-
Notifications
You must be signed in to change notification settings - Fork 18
/
PlotGeneSet.R
148 lines (121 loc) · 7.53 KB
/
PlotGeneSet.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#' @title Plot the Expression Profiles of a Gene Set
#' @description
#' This function simply visualizes the gene expression profiles of
#' a defined subset of genes stored in the input \code{ExpressionSet}.
#' @param ExpressionSet a standard PhyloExpressionSet or DivergenceExpressionSet object.
#' @param gene.set a character vector storing the gene ids for which gene expression profiles shall be visualized.
#' @param get.subset a logical value indicating whether or not an \code{ExpressionSet} subset of the selected \code{gene.set} should be retuned.
#' @param use.only.map a logical value indicating whether instead of a standard \code{ExpressionSet} only a \code{Phylostratigraphic Map} or \code{Divergene Map} is passed to the function.
#' @param colors colors for gene expression profiles. Default: \code{colors = NULL}, hence default colours are used.
#' @param plot.legend a logical value indicating whether gene ids should be printed as legend next to the plot.
#' @param y.ticks a numeric value specifying the number of ticks to be drawn on the y-axis.
#' @param digits.ylab a numeric value specifying the number of digits shown for the expression levels on the y-axis.
#' @param ... additional parameters passed to \code{\link{matplot}}.
#' @author Hajk-Georg Drost
#' @details
#'
#' This function simply visualizes or subsets the gene expression levels of a set of genes
#' that are stored in the input \code{ExpressionSet}.
#'
#' @seealso \code{\link{SelectGeneSet}}, \code{\link{PlotEnrichment}}, \code{\link{DiffGenes}}
#' @examples
#' data(PhyloExpressionSetExample)
#'
#' # the best parameter setting to visualize this plot:
#' # png("test_png.png",700,400)
#' PlotGeneSet(ExpressionSet = PhyloExpressionSetExample,
#' gene.set = PhyloExpressionSetExample[1:5, 2],
#' type = "l",
#' lty = 1,
#' lwd = 4,
#' xlab = "Ontogeny",
#' ylab = "Expression Level")
#'
#' # dev.off()
#'
#' # In case you would like to work with the expression levels
#' # of selected genes you can specify the 'get.subset' argument:
#'
#' PlotGeneSet(ExpressionSet = PhyloExpressionSetExample,
#' gene.set = PhyloExpressionSetExample[1:5, 2],
#' get.subset = TRUE)
#'
#'
#' # get a gene subset using only a phylostratihraphic map
#' ExamplePSMap <- PhyloExpressionSetExample[ , 1:2]
#'
#' PlotGeneSet(ExpressionSet = ExamplePSMap,
#' gene.set = PhyloExpressionSetExample[1:5, 2],
#' get.subset = TRUE,
#' use.only.map = TRUE)
#'
#' @export
PlotGeneSet <- function(ExpressionSet,
gene.set,
get.subset = FALSE,
use.only.map = FALSE,
colors = NULL,
plot.legend = TRUE,
y.ticks = 6,
digits.ylab = 4, ... ){
ExpressionSet <- as.data.frame(ExpressionSet)
if (!use.only.map)
is.ExpressionSet(ExpressionSet)
GeneSubSet.indixes <- stats::na.omit(match(tolower(gene.set), tolower(ExpressionSet[ , 2])))
if (length(GeneSubSet.indixes) == 0)
stop ("None of your input gene ids could be found in the ExpressionSet.", call. = FALSE)
if (length(GeneSubSet.indixes) != length(gene.set))
warning ("Only ",length(GeneSubSet.indixes), " out of your ", length(gene.set), " gene ids could be found in the ExpressionSet.", call. = FALSE)
GeneSubSet <- ExpressionSet[GeneSubSet.indixes , ]
ncols <- ncol(GeneSubSet)
if(!is.null(colors)){
if (length(colors) != length(gene.set))
stop ("The number of colors and the number of genes do not match.", call. = FALSE)
}
# http://www.compbiome.com/2010/12/r-using-rcolorbrewer-to-colour-your.html
if (is.null(colors))
colors <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(8,"Dark2"))(length(GeneSubSet.indixes))
if(!get.subset){
# define arguments for different graphics functions
plot.args <- c("type","lwd","col","cex.lab","main","xlab","ylab")
axis.args <- c("las", "cex.axis")
legend.args <- c("border","angle","density","box.lwd","cex")
dots <- list(...)
ellipsis.names <- names(dots)
ylim.range <- range(min(GeneSubSet[ , 3:ncols]),max(GeneSubSet[ , 3:ncols]))
if((length(ellipsis.names[grep("ylab",ellipsis.names)]) > 0) | (length(ellipsis.names[grep("xlab",ellipsis.names)]) > 0)){
if(plot.legend)
graphics::par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = TRUE)
do.call(graphics::matplot,c(list(x = t(GeneSubSet[ , 3:ncols]),
col = colors[1:length(GeneSubSet.indixes)],
axes = FALSE),
dots[!is.element(names(dots),c(axis.args,legend.args))]))
} else {
if(plot.legend)
graphics::par(mar = c(5.1, 4.1, 4.1, 8.1), xpd = TRUE)
do.call(graphics::matplot,c(list(x = t(GeneSubSet[ , 3:ncols]),
col = colors[1:length(GeneSubSet.indixes)],
xaxt = "n",
xlab = "Ontogeny",
ylab = "Expression Level",
axes = FALSE ),
dots[!is.element(names(dots),c(axis.args,legend.args))]))
}
do.call(graphics::axis,c(list(side = 1,at = 1:(ncols-2),
labels = names(ExpressionSet)[3:ncols]),dots[!is.element(names(dots),c(plot.args,legend.args))]))
do.call(graphics::axis,c(list(side = 2,at = format(seq(ylim.range[1],ylim.range[2],length.out = y.ticks),digits = digits.ylab),
labels = format(seq(ylim.range[1],ylim.range[2],length.out = y.ticks),digits = digits.ylab)),
dots[!is.element(names(dots),c(plot.args,legend.args))]))
if(plot.legend){
do.call(graphics::legend, c(list("topright",
inset = c(-0.2,0),
legend = GeneSubSet[ , 2],
title = "Genes",
col = colors[1:length(GeneSubSet.indixes)],
lwd = 4,
bty = "n"), dots[!is.element(names(dots),c(plot.args,axis.args))]))
}
}
if (get.subset)
return(GeneSubSet)
}