Browse code

First revision

liuwd15 authored on 06/10/2020 08:06:07
Showing 1 changed files
... ...
@@ -13,27 +13,26 @@
13 13
 #' @seealso \code{geom_smooth} for plotting smooth lines, \code{facet_wrap} for faceting genes.
14 14
 #'
15 15
 #' @importFrom ggplot2 ggplot aes_string theme theme_bw geom_smooth geom_line labs facet_wrap
16
-#' @importFrom SummarizedExperiment rowData colData
16
+#' @importFrom SummarizedExperiment assay rowData colData
17 17
 #' @export
18 18
 #'
19 19
 #' @examples
20 20
 #' data(zh.data)
21
-#' zh <- CreateTomo(zh.data)
22
-#' zh <- Normalize(zh)
23
-#' LinePlot(zh,
21
+#' zh <- createTomo(zh.data)
22
+#' linePlot(zh,
24 23
 #'  c("ENSDARG00000002131", "ENSDARG00000003061", "ENSDARG00000076075", "ENSDARG00000076850"))
25 24
 #'
26 25
 #' # Do not smooth lines.
27
-#' LinePlot(zh,
26
+#' linePlot(zh,
28 27
 #'  c("ENSDARG00000002131", "ENSDARG00000003061", "ENSDARG00000076075", "ENSDARG00000076850"), span=0)
29 28
 #'
30 29
 #' # Plot genes in different facets.
31
-#' LinePlot(zh,
30
+#' linePlot(zh,
32 31
 #'  c("ENSDARG00000002131", "ENSDARG00000003061", "ENSDARG00000076075", "ENSDARG00000076850"),
33 32
 #'  facet=TRUE)
34
-LinePlot <- function(object, genes, matrix='normalized', facet=FALSE, span=0.3)
33
+linePlot <- function(object, genes, matrix='normalized', facet=FALSE, span=0.3)
35 34
 {
36
-    exp_matrix <- GetMatrix(object, matrix)
35
+    exp_matrix <- assay(object, matrix)
37 36
     exp_genes <- subset(exp_matrix, rowData(object)$gene %in% genes)
38 37
     exp_genes_df <- NULL
39 38
     sections <- factor(colData(object)$section, levels=colData(object)$section)
... ...
@@ -45,7 +44,7 @@ LinePlot <- function(object, genes, matrix='normalized', facet=FALSE, span=0.3)
45 44
     ylab <- c('Count', 'Normalized count', 'Scaled expression')[c('count', 'normalized', 'scaled') == matrix]
46 45
 
47 46
     g <- ggplot(exp_genes_df, aes_string(x='section', y='value', group='gene', color='gene')) +
48
-        labs(x='', y=ylab) +
47
+        labs(y=ylab) +
49 48
         theme_bw() +
50 49
         theme(legend.position='top', axis.text.x=element_text(angle=90, hjust=1, vjust=0.5))
51 50
 
... ...
@@ -78,25 +77,24 @@ LinePlot <- function(object, genes, matrix='normalized', facet=FALSE, span=0.3)
78 77
 #'
79 78
 #' @examples
80 79
 #' data(zh.data)
81
-#' zh <- CreateTomo(zh.data)
82
-#' zh <- Normalize(zh)
83
-#' zh <- TSNE(zh)
80
+#' zh <- createTomo(zh.data)
81
+#' zh <- runTSNE(zh)
84 82
 #' # Plot TSNE embeddings.
85
-#' EmbedPlot(zh)
83
+#' embedPlot(zh)
86 84
 #'
87 85
 #' # Plot UMAP embeddings.
88
-#' zh <- UMAP(zh)
89
-#' EmbedPlot(zh, method="UMAP")
86
+#' zh <- runUMAP(zh)
87
+#' embedPlot(zh, method="UMAP")
90 88
 #'
91 89
 #' # Color sections by kmeans cluster labels.
92
-#' zh <- KMeans(zh, 3)
93
-#' EmbedPlot(zh, group="kmeans_cluster")
94
-EmbedPlot <- function(object, group='section', method='TSNE')
90
+#' zh <- kmeansClust(zh, 3)
91
+#' embedPlot(zh, group="kmeans_cluster")
92
+embedPlot <- function(object, group='section', method='TSNE')
95 93
 {
96 94
     if(!group %in% names(colData(object)))
97 95
     {
98 96
         group <- 'section'
99
-        cat("Unknown parameter 'group' for plotting sections!")
97
+        message("Unknown parameter 'group' for plotting sections!")
100 98
     }
101 99
     colData(object)[[group]] <- as.character(colData(object)[[group]])
102 100
 
... ...
@@ -113,7 +111,7 @@ EmbedPlot <- function(object, group='section', method='TSNE')
113 111
         }
114 112
         else
115 113
         {
116
-            cat("Function 'TSNE' must be run before plotting TSNE embeddings.\n")
114
+            message("Function 'runTSNE' must be run before plotting TSNE embeddings.\n")
117 115
         }
118 116
     }
119 117
     else if(method == 'UMAP')
... ...
@@ -129,7 +127,7 @@ EmbedPlot <- function(object, group='section', method='TSNE')
129 127
         }
130 128
         else
131 129
         {
132
-            cat("Function 'UMAP' must be run before plotting UMAP embeddings.\n")
130
+            message("Function 'runUMAP' must be run before plotting UMAP embeddings.\n")
133 131
         }
134 132
     }
135 133
     else if(method == 'PCA')
... ...
@@ -145,12 +143,12 @@ EmbedPlot <- function(object, group='section', method='TSNE')
145 143
         }
146 144
         else
147 145
         {
148
-            cat("Function 'PCA' must be run before plotting PC embeddings.\n")
146
+            message("Function 'runPCA' must be run before plotting PC embeddings.\n")
149 147
         }
150 148
     }
151 149
     else
152 150
     {
153
-        cat("Unknown embeddings!\n")
151
+        message("Unknown embeddings!\n")
154 152
     }
155 153
 }
156 154
 
... ...
@@ -165,26 +163,25 @@ EmbedPlot <- function(object, group='section', method='TSNE')
165 163
 #'
166 164
 #' @return A \code{ggplot} object.
167 165
 #'
168
-#' @importFrom SummarizedExperiment assays rowData
166
+#' @importFrom SummarizedExperiment rowData
169 167
 #' @importFrom ggplot2 ggplot aes_string geom_point theme theme_bw
170 168
 #'
171 169
 #' @export
172 170
 #'
173 171
 #' @examples
174 172
 #' data(zh.data)
175
-#' zh <- CreateTomo(zh.data)
176
-#' zh <- Normalize(zh)
177
-#' peak_genes <- FindPeakGene(zh)
178
-#' zh <- TSNE(zh, peak_genes$gene)
173
+#' zh <- createTomo(zh.data)
174
+#' peak_genes <- findPeakGene(zh)
175
+#' zh <- runTSNE(zh, peak_genes$gene)
179 176
 #' # Color genes by peak centers.
180
-#' GeneEmbedPlot(zh, peak_genes)
177
+#' geneEmbedPlot(zh, peak_genes)
181 178
 #'
182 179
 #' # Color genes by peak starts.
183
-#' GeneEmbedPlot(zh, peak_genes, group="start")
180
+#' geneEmbedPlot(zh, peak_genes, group="start")
184 181
 #'
185 182
 #' # Do not color genes.
186
-#' GeneEmbedPlot(zh, peak_genes["gene"])
187
-GeneEmbedPlot <- function(object, gene.df, group='center', method='TSNE')
183
+#' geneEmbedPlot(zh, peak_genes["gene"])
184
+geneEmbedPlot <- function(object, gene.df, group='center', method='TSNE')
188 185
 {
189 186
     if(group %in% names(gene.df))
190 187
         gene.df[[group]] <- as.factor(gene.df[[group]])
... ...
@@ -206,7 +203,7 @@ GeneEmbedPlot <- function(object, gene.df, group='center', method='TSNE')
206 203
         }
207 204
         else
208 205
         {
209
-            cat("Function 'TSNE' must be run for input genes before plotting TSNE embeddings.\n")
206
+            message("Function 'TSNE' must be run for input genes before plotting TSNE embeddings.\n")
210 207
         }
211 208
     }
212 209
     else if(method == 'UMAP')
... ...
@@ -227,7 +224,7 @@ GeneEmbedPlot <- function(object, gene.df, group='center', method='TSNE')
227 224
         }
228 225
         else
229 226
         {
230
-            cat("Function 'UMAP' must be run for input genes before plotting UMAP embeddings.\n")
227
+            message("Function 'UMAP' must be run for input genes before plotting UMAP embeddings.\n")
231 228
         }
232 229
     }
233 230
     else if(method == 'PCA')
... ...
@@ -248,12 +245,12 @@ GeneEmbedPlot <- function(object, gene.df, group='center', method='TSNE')
248 245
         }
249 246
         else
250 247
         {
251
-            cat("Function 'PCA' must be run for input genes before plotting PC embeddings.\n")
248
+            message("Function 'runPCA' must be run for input genes before plotting PC embeddings.\n")
252 249
         }
253 250
     }
254 251
     else
255 252
     {
256
-        cat("Unknown embeddings!\n")
253
+        message("Unknown embeddings!\n")
257 254
     }
258 255
 }
259 256
 
... ...
@@ -268,39 +265,39 @@ GeneEmbedPlot <- function(object, gene.df, group='center', method='TSNE')
268 265
 #'
269 266
 #' @return A \code{ggplot} object.
270 267
 #'
268
+#' @importFrom SummarizedExperiment assay
271 269
 #' @importFrom reshape2 melt
272
-#' @importFrom ggplot2 ggplot aes_string geom_raster scale_fill_gradient2 theme element_blank element_text
270
+#' @importFrom RColorBrewer brewer.pal
271
+#' @importFrom ggplot2 ggplot aes_string geom_raster scale_fill_distiller theme element_blank element_text
273 272
 #'
274 273
 #' @export
275 274
 #'
276 275
 #' @examples
277 276
 #' data(zh.data)
278
-#' zh <- CreateTomo(zh.data)
279
-#' zh <- Normalize(zh)
280
-#' zh <- Scale(zh)
277
+#' zh <- createTomo(zh.data)
281 278
 #'
282 279
 #' # Plot some genes.
283
-#' ExpHeatmap(zh,
280
+#' expHeatmap(zh,
284 281
 #'  c("ENSDARG00000002131", "ENSDARG00000003061", "ENSDARG00000076075", "ENSDARG00000076850"))
285 282
 #'
286 283
 #' # Plot peak genes.
287
-#' peak_genes <- FindPeakGene(zh)
288
-#' ExpHeatmap(zh, peak_genes$gene)
284
+#' peak_genes <- findPeakGene(zh)
285
+#' expHeatmap(zh, peak_genes$gene)
289 286
 #'
290 287
 #' # Remove gene names if too many genes are in the heatmap.
291
-#' ExpHeatmap(zh, peak_genes$gene, size=0)
292
-ExpHeatmap <- function(object, genes, matrix='scaled', size=5)
288
+#' expHeatmap(zh, peak_genes$gene, size=0)
289
+expHeatmap <- function(object, genes, matrix='scaled', size=5)
293 290
 {
294 291
     genes <- rev(as.character(genes))
295
-    exp_matrix <- GetMatrix(object, matrix)[genes, ]
296
-    if(matrix == 'normalized')
292
+    exp_matrix <- assay(object, matrix)[genes, ]
293
+    if(matrix %in% c('count', 'normalized'))
297 294
         exp_matrix <- log10(exp_matrix + 1)
298 295
     genes_df <- melt(exp_matrix, varnames=c('gene','section'))
299 296
     exp_name <-  c('Log10(Count+1)', 'Log10(Normalized count+1)', 'Scaled expression\n(Z score)')[c('count', 'normalized', 'scaled') == matrix]
300 297
 
301 298
     g <- ggplot(genes_df, aes_string(x='section', y='gene', fill='value')) +
302 299
         geom_raster() +
303
-        scale_fill_gradient2(low='magenta', mid='black', high='yellow', name=exp_name) +
300
+        scale_fill_distiller(palette='RdYlBu', name=exp_name) +
304 301
         theme(axis.line=element_blank(),
305 302
               axis.text.x=element_text(angle=90, hjust=1, vjust=0.5),
306 303
               axis.ticks=element_blank(),
... ...
@@ -324,28 +321,27 @@ ExpHeatmap <- function(object, genes, matrix='scaled', size=5)
324 321
 #'
325 322
 #' @return A \code{ggplot} object.
326 323
 #'
324
+#' @importFrom SummarizedExperiment assay
327 325
 #' @importFrom stats cor
328 326
 #' @importFrom reshape2 melt
329
-#' @importFrom ggplot2 ggplot aes_string geom_raster scale_fill_gradientn labs theme element_blank element_text
330
-#' @importFrom colorRamps rgb.tables
327
+#' @importFrom RColorBrewer brewer.pal
328
+#' @importFrom ggplot2 ggplot aes_string geom_raster scale_fill_distiller labs theme element_blank element_text
331 329
 #'
332 330
 #' @export
333 331
 #'
334 332
 #' @examples
335 333
 #' data(zh.data)
336
-#' zh <- CreateTomo(zh.data)
337
-#' zh <- Normalize(zh)
338
-#' zh <- Scale(zh)
339
-#' CorHeatmap(zh)
334
+#' zh <- createTomo(zh.data)
335
+#' corHeatmap(zh)
340 336
 #'
341 337
 #' # Use Spearman correlation coefficients.
342
-#' CorHeatmap(zh, cor.method='spearman')
338
+#' corHeatmap(zh, cor.method='spearman')
343 339
 #'
344 340
 #' # Set max correlation coefficients to 0.3.
345
-#' CorHeatmap(zh, max.cor=0.3)
346
-CorHeatmap <- function(object, matrix='scaled', max.cor=0.5, cor.method='pearson')
341
+#' corHeatmap(zh, max.cor=0.3)
342
+corHeatmap <- function(object, matrix='scaled', max.cor=0.5, cor.method='pearson')
347 343
 {
348
-    exp_matrix <- GetMatrix(object, matrix=matrix)
344
+    exp_matrix <- assay(object, matrix)
349 345
     cor_matrix <- cor(exp_matrix, method=cor.method)
350 346
     cor_matrix[cor_matrix > max.cor] <- max.cor
351 347
     cor_df <- melt(cor_matrix, varnames=c('section1','section2'))
... ...
@@ -354,7 +350,7 @@ CorHeatmap <- function(object, matrix='scaled', max.cor=0.5, cor.method='pearson
354 350
 
355 351
     g <- ggplot(cor_df, aes_string(x='section1', y='section2', fill='value')) +
356 352
         geom_raster() +
357
-        scale_fill_gradientn(colors=rgb.tables(10), name=cor_name) +
353
+        scale_fill_distiller(palette='RdYlBu', name=cor_name) +
358 354
         labs(x='', y='') +
359 355
         theme(axis.line=element_blank(),
360 356
               axis.text.x=element_text(angle=90, hjust=1, vjust=0.5),
... ...
@@ -377,39 +373,39 @@ CorHeatmap <- function(object, matrix='scaled', max.cor=0.5, cor.method='pearson
377 373
 #'
378 374
 #' @details This method can create a pure heatmap or a heatmap with side bar. If you prefer a pure heatmap, input a \code{gene.df} with a single column of gene names.
379 375
 #' However, you may want to show additional information of genes with a side bar, and the grouping information should be saved as additional column(s) of \code{gene.df}, and declared as \code{group}.
380
-#' By default, you can use the output by \code{FindPeakGene} as input \code{gene.df}. Peak genes will be grouped by their centers on the side bar.
376
+#' By default, you can use the output by \code{findPeakGene} as input \code{gene.df}. Peak genes will be grouped by their centers on the side bar.
381 377
 #'
382 378
 #' @return A \code{ggplot} object.
383 379
 #'
380
+#' @importFrom SummarizedExperiment assay
384 381
 #' @importFrom stats cor
385 382
 #' @importFrom grDevices rainbow
386 383
 #' @importFrom reshape2 melt
387
-#' @importFrom ggplot2 ggplot aes_string geom_raster scale_fill_gradientn labs theme element_blank element_text ggplot_build scale_color_manual guides annotation_raster guide_legend coord_cartesian
384
+#' @importFrom RColorBrewer brewer.pal
385
+#' @importFrom ggplot2 ggplot aes_string geom_raster scale_fill_distiller labs theme element_blank element_text ggplot_build scale_color_manual guides annotation_raster guide_legend coord_cartesian
388 386
 #'
389 387
 #' @export
390 388
 #'
391 389
 #' @examples
392 390
 #' data(zh.data)
393
-#' zh <- CreateTomo(zh.data)
394
-#' zh <- Normalize(zh)
395
-#' zh <- Scale(zh)
391
+#' zh <- createTomo(zh.data)
396 392
 #'
397 393
 #' # Correlation heatmap for all peak genes.
398
-#' peak_genes <- FindPeakGene(zh)
399
-#' GeneCorHeatmap(zh, peak_genes)
394
+#' peak_genes <- findPeakGene(zh)
395
+#' geneCorHeatmap(zh, peak_genes)
400 396
 #'
401 397
 #' # Use Spearman correlation coefficients.
402
-#' GeneCorHeatmap(zh, peak_genes, cor.method="spearman")
398
+#' geneCorHeatmap(zh, peak_genes, cor.method="spearman")
403 399
 #'
404 400
 #' # Group genes by peak start.
405
-#' GeneCorHeatmap(zh, peak_genes, group="start")
401
+#' geneCorHeatmap(zh, peak_genes, group="start")
406 402
 #'
407 403
 #' # Plot without side bar.
408
-#' GeneCorHeatmap(zh, data.frame(
404
+#' geneCorHeatmap(zh, data.frame(
409 405
 #'  gene=c("ENSDARG00000002131", "ENSDARG00000003061", "ENSDARG00000076075", "ENSDARG00000076850")))
410
-GeneCorHeatmap <- function(object, gene.df, group='center', matrix='scaled', size=5, cor.method='pearson')
406
+geneCorHeatmap <- function(object, gene.df, group='center', matrix='scaled', size=5, cor.method='pearson')
411 407
 {
412
-    exp_matrix <- GetMatrix(object, matrix=matrix)
408
+    exp_matrix <- assay(object, matrix)
413 409
     genes <- as.character(gene.df$gene)
414 410
     cor_matrix <- cor(t(exp_matrix[genes, ]), method=cor.method)
415 411
     cor_df <- melt(cor_matrix, varnames=c('gene1', 'gene2'))
... ...
@@ -430,7 +426,7 @@ GeneCorHeatmap <- function(object, gene.df, group='center', matrix='scaled', siz
430 426
 
431 427
     g <- ggplot(cor_df, aes_string(x='gene1', y='gene2', fill='value')) +
432 428
         geom_raster() +
433
-        scale_fill_gradientn(colors=rgb.tables(10), name=cor_name)  +
429
+        scale_fill_distiller(palette='RdYlBu', name=cor_name)  +
434 430
         labs(x='', y='') +
435 431
         theme(axis.line=element_blank(),
436 432
               axis.text.y=element_blank(),
... ...
@@ -451,9 +447,10 @@ GeneCorHeatmap <- function(object, gene.df, group='center', matrix='scaled', siz
451 447
         y_max <- y_min + 0.03 * y_range
452 448
 
453 449
         g <- g + geom_point(aes_string(color='group'), alpha=0, shape=15, size=5) +
454
-            scale_color_manual(name=group, values=color_legend)  +
450
+            scale_color_manual(name=group, values=color_legend) +
455 451
             guides(color=guide_legend(override.aes=list(alpha=1))) +
456
-            annotation_raster(t(color_legend[as.character(group_genes)]), -Inf, Inf, ymin=y_min, ymax=y_max) +
452
+            annotation_raster(t(color_legend[as.character(group_genes)]),
453
+                              -Inf, Inf, ymin=y_min, ymax=y_max) +
457 454
             coord_cartesian(ylim = c(0, y_max), clip='off')
458 455
     }
459 456
 
Browse code

Edit files to pass BiocCheck

liuwd15 authored on 15/09/2020 02:53:28
Showing 1 changed files
... ...
@@ -33,31 +33,31 @@
33 33
 #'  facet=TRUE)
34 34
 LinePlot <- function(object, genes, matrix='normalized', facet=FALSE, span=0.3)
35 35
 {
36
-  exp_matrix <- GetMatrix(object, matrix)
37
-  exp_genes <- subset(exp_matrix, rowData(object)$gene %in% genes)
38
-  exp_genes_df <- NULL
39
-  sections <- factor(colData(object)$section, levels=colData(object)$section)
40
-  for(i in 1:nrow(exp_genes))
41
-  {
42
-    exp_gene_df <- data.frame(gene=genes[i], section=sections, value=exp_genes[i,])
43
-    exp_genes_df <- rbind(exp_genes_df, exp_gene_df)
44
-  }
45
-  ylab <- c('Count', 'Normalized count', 'Scaled expression')[c('count', 'normalized', 'scaled') == matrix]
36
+    exp_matrix <- GetMatrix(object, matrix)
37
+    exp_genes <- subset(exp_matrix, rowData(object)$gene %in% genes)
38
+    exp_genes_df <- NULL
39
+    sections <- factor(colData(object)$section, levels=colData(object)$section)
40
+    for(i in seq_len(nrow(exp_genes)))
41
+    {
42
+        exp_gene_df <- data.frame(gene=genes[i], section=sections, value=exp_genes[i,])
43
+        exp_genes_df <- rbind(exp_genes_df, exp_gene_df)
44
+    }
45
+    ylab <- c('Count', 'Normalized count', 'Scaled expression')[c('count', 'normalized', 'scaled') == matrix]
46 46
 
47
-  g <- ggplot(exp_genes_df, aes_string(x='section', y='value', group='gene', color='gene')) +
48
-    labs(x='', y=ylab) +
49
-    theme_bw() +
50
-    theme(legend.position='top', axis.text.x=element_text(angle=90, hjust=1, vjust=0.5))
47
+    g <- ggplot(exp_genes_df, aes_string(x='section', y='value', group='gene', color='gene')) +
48
+        labs(x='', y=ylab) +
49
+        theme_bw() +
50
+        theme(legend.position='top', axis.text.x=element_text(angle=90, hjust=1, vjust=0.5))
51 51
 
52
-  if(span > 0)
53
-    g <- g + geom_smooth(method='loess', span = span, se=F)
54
-  else
55
-    g <- g + geom_line(size=1)
52
+    if(span > 0)
53
+        g <- g + geom_smooth(method='loess', span=span, se=FALSE)
54
+    else
55
+        g <- g + geom_line(size=1)
56 56
 
57
-  if(facet)
58
-    g <- g + facet_wrap(~gene, scales = 'free') + theme(legend.position='none')
57
+    if(facet)
58
+        g <- g + facet_wrap(~gene, scales='free') + theme(legend.position='none')
59 59
 
60
-  return(g)
60
+    return(g)
61 61
 }
62 62
 
63 63
 #' Embedding plot for sections
... ...
@@ -93,65 +93,65 @@ LinePlot <- function(object, genes, matrix='normalized', facet=FALSE, span=0.3)
93 93
 #' EmbedPlot(zh, group="kmeans_cluster")
94 94
 EmbedPlot <- function(object, group='section', method='TSNE')
95 95
 {
96
-  if(!group %in% names(colData(object)))
97
-  {
98
-    group <- 'section'
99
-    cat("Unknown parameter 'group' for plotting sections!")
100
-  }
101
-  colData(object)[[group]] <- as.character(colData(object)[[group]])
102
-
103
-  if(method == 'TSNE')
104
-  {
105
-    if(all(c('TSNE1','TSNE2') %in% names(colData(object))))
96
+    if(!group %in% names(colData(object)))
106 97
     {
107
-      g <- ggplot(data.frame(colData(object)), aes_string(x='TSNE1', y='TSNE2', color=group)) +
108
-        geom_point() +
109
-        geom_text_repel(aes_string(label='section')) +
110
-        theme_bw() +
111
-        theme(legend.position='none')
112
-      return(g)
98
+        group <- 'section'
99
+        cat("Unknown parameter 'group' for plotting sections!")
113 100
     }
114
-    else
115
-    {
116
-      cat("Function 'TSNE' must be run before plotting TSNE embeddings.\n")
117
-    }
118
-  }
119
-  else if(method == 'UMAP')
120
-  {
121
-    if(all(c('UMAP1','UMAP2') %in% names(colData(object))))
101
+    colData(object)[[group]] <- as.character(colData(object)[[group]])
102
+
103
+    if(method == 'TSNE')
122 104
     {
123
-      g <- ggplot(data.frame(colData(object)), aes_string(x='UMAP1', y='UMAP2', color=group)) +
124
-        geom_point(aes_string()) +
125
-        geom_text_repel(aes_string(label='section')) +
126
-        theme_bw() +
127
-        theme(legend.position='none')
128
-      return(g)
105
+        if(all(c('TSNE1','TSNE2') %in% names(colData(object))))
106
+        {
107
+            g <- ggplot(data.frame(colData(object)), aes_string(x='TSNE1', y='TSNE2', color=group)) +
108
+                geom_point() +
109
+                geom_text_repel(aes_string(label='section')) +
110
+                theme_bw() +
111
+                theme(legend.position='none')
112
+            return(g)
113
+        }
114
+        else
115
+        {
116
+            cat("Function 'TSNE' must be run before plotting TSNE embeddings.\n")
117
+        }
129 118
     }
130
-    else
119
+    else if(method == 'UMAP')
131 120
     {
132
-      cat("Function 'UMAP' must be run before plotting UMAP embeddings.\n")
121
+        if(all(c('UMAP1','UMAP2') %in% names(colData(object))))
122
+        {
123
+            g <- ggplot(data.frame(colData(object)), aes_string(x='UMAP1', y='UMAP2', color=group)) +
124
+                geom_point(aes_string()) +
125
+                geom_text_repel(aes_string(label='section')) +
126
+                theme_bw() +
127
+                theme(legend.position='none')
128
+            return(g)
129
+        }
130
+        else
131
+        {
132
+            cat("Function 'UMAP' must be run before plotting UMAP embeddings.\n")
133
+        }
133 134
     }
134
-  }
135
-  else if(method == 'PCA')
136
-  {
137
-    if(all(c('PC1','PC2') %in% names(colData(object))))
135
+    else if(method == 'PCA')
138 136
     {
139
-      g <- ggplot(data.frame(colData(object)), aes_string(x='PC1', y='PC2', color=group)) +
140
-        geom_point(aes_string()) +
141
-        geom_text_repel(aes_string(label='section')) +
142
-        theme_bw() +
143
-        theme(legend.position='none')
144
-      return(g)
137
+        if(all(c('PC1','PC2') %in% names(colData(object))))
138
+        {
139
+            g <- ggplot(data.frame(colData(object)), aes_string(x='PC1', y='PC2', color=group)) +
140
+                geom_point(aes_string()) +
141
+                geom_text_repel(aes_string(label='section')) +
142
+                theme_bw() +
143
+                theme(legend.position='none')
144
+            return(g)
145
+        }
146
+        else
147
+        {
148
+            cat("Function 'PCA' must be run before plotting PC embeddings.\n")
149
+        }
145 150
     }
146 151
     else
147 152
     {
148
-      cat("Function 'PCA' must be run before plotting PC embeddings.\n")
153
+        cat("Unknown embeddings!\n")
149 154
     }
150
-  }
151
-  else
152
-  {
153
-    cat("Unknown embeddings!\n")
154
-  }
155 155
 }
156 156
 
157 157
 #' Embedding plot for genes
... ...
@@ -186,75 +186,75 @@ EmbedPlot <- function(object, group='section', method='TSNE')
186 186
 #' GeneEmbedPlot(zh, peak_genes["gene"])
187 187
 GeneEmbedPlot <- function(object, gene.df, group='center', method='TSNE')
188 188
 {
189
-  if(group %in% names(gene.df))
190
-    gene.df[[group]] <- as.factor(gene.df[[group]])
191
-  if(method == 'TSNE')
192
-  {
193
-    if(all(c('TSNE1','TSNE2') %in% names(rowData(object))))
189
+    if(group %in% names(gene.df))
190
+        gene.df[[group]] <- as.factor(gene.df[[group]])
191
+    if(method == 'TSNE')
194 192
     {
195
-      gene.df$TSNE1 <- rowData(object)[as.character(gene.df$gene), 'TSNE1']
196
-      gene.df$TSNE2 <- rowData(object)[as.character(gene.df$gene), 'TSNE2']
197
-      g <- ggplot(gene.df, aes_string(x='TSNE1', y='TSNE2')) +
198
-        theme_bw()
193
+        if(all(c('TSNE1','TSNE2') %in% names(rowData(object))))
194
+        {
195
+            gene.df$TSNE1 <- rowData(object)[as.character(gene.df$gene), 'TSNE1']
196
+            gene.df$TSNE2 <- rowData(object)[as.character(gene.df$gene), 'TSNE2']
197
+            g <- ggplot(gene.df, aes_string(x='TSNE1', y='TSNE2')) +
198
+                theme_bw()
199 199
 
200
-      if(group %in% names(gene.df))
201
-        g <- g + geom_point(aes_string(color=group))
202
-      else
203
-        g <- g + geom_point()
200
+            if(group %in% names(gene.df))
201
+                g <- g + geom_point(aes_string(color=group))
202
+            else
203
+                g <- g + geom_point()
204 204
 
205
-      return(g)
205
+            return(g)
206
+        }
207
+        else
208
+        {
209
+            cat("Function 'TSNE' must be run for input genes before plotting TSNE embeddings.\n")
210
+        }
206 211
     }
207
-    else
208
-    {
209
-      cat("Function 'TSNE' must be run for input genes before plotting TSNE embeddings.\n")
210
-    }
211
-  }
212
-  else if(method == 'UMAP')
213
-  {
214
-    if(all(c('UMAP1','UMAP2') %in% names(rowData(object))))
212
+    else if(method == 'UMAP')
215 213
     {
216
-      gene.df$UMAP1 <- rowData(object)[as.character(gene.df$gene), 'UMAP1']
217
-      gene.df$UMAP2 <- rowData(object)[as.character(gene.df$gene), 'UMAP2']
218
-      g <- ggplot(gene.df, aes_string(x='UMAP1', y='UMAP2')) +
219
-        theme_bw()
214
+        if(all(c('UMAP1','UMAP2') %in% names(rowData(object))))
215
+        {
216
+            gene.df$UMAP1 <- rowData(object)[as.character(gene.df$gene), 'UMAP1']
217
+            gene.df$UMAP2 <- rowData(object)[as.character(gene.df$gene), 'UMAP2']
218
+            g <- ggplot(gene.df, aes_string(x='UMAP1', y='UMAP2')) +
219
+                theme_bw()
220 220
 
221
-      if(group %in% names(gene.df))
222
-        g <- g + geom_point(aes_string(color=group))
223
-      else
224
-        g <- g + geom_point()
221
+            if(group %in% names(gene.df))
222
+                g <- g + geom_point(aes_string(color=group))
223
+            else
224
+                g <- g + geom_point()
225 225
 
226
-      return(g)
227
-    }
228
-    else
229
-    {
230
-      cat("Function 'UMAP' must be run for input genes before plotting UMAP embeddings.\n")
226
+            return(g)
227
+        }
228
+        else
229
+        {
230
+            cat("Function 'UMAP' must be run for input genes before plotting UMAP embeddings.\n")
231
+        }
231 232
     }
232
-  }
233
-  else if(method == 'PCA')
234
-  {
235
-    if(all(c('PC1','PC2') %in% names(rowData(object))))
233
+    else if(method == 'PCA')
236 234
     {
237
-      gene.df$PC1 <- rowData(object)[as.character(gene.df$gene), 'PC1']
238
-      gene.df$PC2 <- rowData(object)[as.character(gene.df$gene), 'PC2']
239
-      g <- ggplot(gene.df, aes_string(x='PC1', y='PC2')) +
240
-        theme_bw()
235
+        if(all(c('PC1','PC2') %in% names(rowData(object))))
236
+        {
237
+            gene.df$PC1 <- rowData(object)[as.character(gene.df$gene), 'PC1']
238
+            gene.df$PC2 <- rowData(object)[as.character(gene.df$gene), 'PC2']
239
+            g <- ggplot(gene.df, aes_string(x='PC1', y='PC2')) +
240
+                theme_bw()
241 241
 
242
-      if(group %in% names(gene.df))
243
-        g <- g + geom_point(aes_string(color=group))
244
-      else
245
-        g <- g + geom_point()
242
+            if(group %in% names(gene.df))
243
+                g <- g + geom_point(aes_string(color=group))
244
+            else
245
+                g <- g + geom_point()
246 246
 
247
-      return(g)
247
+            return(g)
248
+        }
249
+        else
250
+        {
251
+            cat("Function 'PCA' must be run for input genes before plotting PC embeddings.\n")
252
+        }
248 253
     }
249 254
     else
250 255
     {
251
-      cat("Function 'PCA' must be run for input genes before plotting PC embeddings.\n")
256
+        cat("Unknown embeddings!\n")
252 257
     }
253
-  }
254
-  else
255
-  {
256
-    cat("Unknown embeddings!\n")
257
-  }
258 258
 }
259 259
 
260 260
 #' Expression heatmap
... ...
@@ -291,26 +291,26 @@ GeneEmbedPlot <- function(object, gene.df, group='center', method='TSNE')
291 291
 #' ExpHeatmap(zh, peak_genes$gene, size=0)
292 292
 ExpHeatmap <- function(object, genes, matrix='scaled', size=5)
293 293
 {
294
-  genes <- rev(as.character(genes))
295
-  exp_matrix <- GetMatrix(object, matrix)[genes, ]
296
-  if(matrix == 'normalized')
297
-    exp_matrix <- log10(exp_matrix + 1)
298
-  genes_df <- melt(exp_matrix, varnames=c('gene','section'))
299
-  exp_name <-  c('Log10(Count+1)', 'Log10(Normalized count+1)', 'Scaled expression\n(Z score)')[c('count', 'normalized', 'scaled') == matrix]
294
+    genes <- rev(as.character(genes))
295
+    exp_matrix <- GetMatrix(object, matrix)[genes, ]
296
+    if(matrix == 'normalized')
297
+        exp_matrix <- log10(exp_matrix + 1)
298
+    genes_df <- melt(exp_matrix, varnames=c('gene','section'))
299
+    exp_name <-  c('Log10(Count+1)', 'Log10(Normalized count+1)', 'Scaled expression\n(Z score)')[c('count', 'normalized', 'scaled') == matrix]
300 300
 
301
-  g <- ggplot(genes_df, aes_string(x='section', y='gene', fill='value')) +
302
-    geom_raster() +
303
-    scale_fill_gradient2(low='magenta', mid='black', high='yellow', name=exp_name) +
304
-    theme(axis.line=element_blank(),
305
-          axis.text.x=element_text(angle=90, hjust=1, vjust=0.5),
306
-          axis.ticks=element_blank(),
307
-          panel.background=element_blank())
308
-  if(size == 0)
309
-    g <- g + theme(axis.text.y=element_blank())
310
-  else
311
-    g <- g + theme(axis.text.y=element_text(size=size))
301
+    g <- ggplot(genes_df, aes_string(x='section', y='gene', fill='value')) +
302
+        geom_raster() +
303
+        scale_fill_gradient2(low='magenta', mid='black', high='yellow', name=exp_name) +
304
+        theme(axis.line=element_blank(),
305
+              axis.text.x=element_text(angle=90, hjust=1, vjust=0.5),
306
+              axis.ticks=element_blank(),
307
+              panel.background=element_blank())
308
+    if(size == 0)
309
+        g <- g + theme(axis.text.y=element_blank())
310
+    else
311
+        g <- g + theme(axis.text.y=element_text(size=size))
312 312
 
313
-  return(g)
313
+    return(g)
314 314
 }
315 315
 
316 316
 #' Correlation heatmap of sections
... ...
@@ -345,23 +345,23 @@ ExpHeatmap <- function(object, genes, matrix='scaled', size=5)
345 345
 #' CorHeatmap(zh, max.cor=0.3)
346 346
 CorHeatmap <- function(object, matrix='scaled', max.cor=0.5, cor.method='pearson')
347 347
 {
348
-  exp_matrix <- GetMatrix(object, matrix=matrix)
349
-  cor_matrix <- cor(exp_matrix, method=cor.method)
350
-  cor_matrix[cor_matrix > max.cor] <- max.cor
351
-  cor_df <- melt(cor_matrix, varnames=c('section1','section2'))
352
-  # c('Pearson r', 'Kendall τ', 'Spearman ρ')
353
-  cor_name <- expression("Pearson"~r,"Kendall"~tau,"Spearman"~rho)[c('pearson', 'kendall', 'spearman') == cor.method]
348
+    exp_matrix <- GetMatrix(object, matrix=matrix)
349
+    cor_matrix <- cor(exp_matrix, method=cor.method)
350
+    cor_matrix[cor_matrix > max.cor] <- max.cor
351
+    cor_df <- melt(cor_matrix, varnames=c('section1','section2'))
352
+    # c('Pearson r', 'Kendall τ', 'Spearman ρ')
353
+    cor_name <- expression("Pearson"~r,"Kendall"~tau,"Spearman"~rho)[c('pearson', 'kendall', 'spearman') == cor.method]
354 354
 
355
-  g <- ggplot(cor_df, aes_string(x='section1', y='section2', fill='value')) +
356
-    geom_raster() +
357
-    scale_fill_gradientn(colors=rgb.tables(10), name=cor_name) +
358
-    labs(x='', y='') +
359
-    theme(axis.line=element_blank(),
360
-          axis.text.x=element_text(angle=90, hjust=1, vjust=0.5),
361
-          axis.ticks=element_blank(),
362
-          panel.background=element_blank())
355
+    g <- ggplot(cor_df, aes_string(x='section1', y='section2', fill='value')) +
356
+        geom_raster() +
357
+        scale_fill_gradientn(colors=rgb.tables(10), name=cor_name) +
358
+        labs(x='', y='') +
359
+        theme(axis.line=element_blank(),
360
+              axis.text.x=element_text(angle=90, hjust=1, vjust=0.5),
361
+              axis.ticks=element_blank(),
362
+              panel.background=element_blank())
363 363
 
364
-  return(g)
364
+    return(g)
365 365
 }
366 366
 
367 367
 #' Correlation heatmap of genes
... ...
@@ -409,53 +409,53 @@ CorHeatmap <- function(object, matrix='scaled', max.cor=0.5, cor.method='pearson
409 409
 #'  gene=c("ENSDARG00000002131", "ENSDARG00000003061", "ENSDARG00000076075", "ENSDARG00000076850")))
410 410
 GeneCorHeatmap <- function(object, gene.df, group='center', matrix='scaled', size=5, cor.method='pearson')
411 411
 {
412
-  exp_matrix <- GetMatrix(object, matrix=matrix)
413
-  genes <- as.character(gene.df$gene)
414
-  cor_matrix <- cor(t(exp_matrix[genes, ]), method=cor.method)
415
-  cor_df <- melt(cor_matrix, varnames=c('gene1', 'gene2'))
416
-  cor_name <- expression("Pearson"~r,"Kendall"~tau,"Spearman"~rho)[c('pearson', 'kendall', 'spearman') == cor.method]
417
-  # plot sidebar if genes are grouped
418
-  plot_sidebar <- group %in% names(gene.df)
412
+    exp_matrix <- GetMatrix(object, matrix=matrix)
413
+    genes <- as.character(gene.df$gene)
414
+    cor_matrix <- cor(t(exp_matrix[genes, ]), method=cor.method)
415
+    cor_df <- melt(cor_matrix, varnames=c('gene1', 'gene2'))
416
+    cor_name <- expression("Pearson"~r,"Kendall"~tau,"Spearman"~rho)[c('pearson', 'kendall', 'spearman') == cor.method]
417
+    # plot sidebar if genes are grouped
418
+    plot_sidebar <- group %in% names(gene.df)
419 419
 
420
-  if(plot_sidebar)
421
-  {
422
-    group_genes <- gene.df[[group]]
423
-    names(group_genes) <- gene.df$gene
424
-    cor_df$group <- as.factor(group_genes[cor_df$gene1])
425
-    # Colors for column side bar
426
-    n_groups <- length(unique(group_genes))
427
-    color_legend <- rainbow(n_groups)
428
-    names(color_legend) <- unique(group_genes)
429
-  }
420
+    if(plot_sidebar)
421
+    {
422
+        group_genes <- gene.df[[group]]
423
+        names(group_genes) <- gene.df$gene
424
+        cor_df$group <- as.factor(group_genes[cor_df$gene1])
425
+        # Colors for column side bar
426
+        n_groups <- length(unique(group_genes))
427
+        color_legend <- rainbow(n_groups)
428
+        names(color_legend) <- unique(group_genes)
429
+    }
430 430
 
431
-  g <- ggplot(cor_df, aes_string(x='gene1', y='gene2', fill='value')) +
432
-    geom_raster() +
433
-    scale_fill_gradientn(colors=rgb.tables(10), name=cor_name)  +
434
-    labs(x='', y='') +
435
-    theme(axis.line=element_blank(),
436
-          axis.text.y=element_blank(),
437
-          axis.ticks=element_blank(),
438
-          panel.background=element_blank(),
439
-          legend.key=element_blank(),
440
-          legend.box='horizontal')
441
-  if(size == 0)
442
-    g <- g + theme(axis.text.x=element_blank())
443
-  else
444
-    g <- g + theme(axis.text.x=element_text(angle=90, hjust=1, size=size))
431
+    g <- ggplot(cor_df, aes_string(x='gene1', y='gene2', fill='value')) +
432
+        geom_raster() +
433
+        scale_fill_gradientn(colors=rgb.tables(10), name=cor_name)  +
434
+        labs(x='', y='') +
435
+        theme(axis.line=element_blank(),
436
+              axis.text.y=element_blank(),
437
+              axis.ticks=element_blank(),
438
+              panel.background=element_blank(),
439
+              legend.key=element_blank(),
440
+              legend.box='horizontal')
441
+    if(size == 0)
442
+        g <- g + theme(axis.text.x=element_blank())
443
+    else
444
+        g <- g + theme(axis.text.x=element_text(angle=90, hjust=1, size=size))
445 445
 
446
-  if(plot_sidebar)
447
-  {
448
-    gbuild <- ggplot_build(plot = g)
449
-    y_range <- diff(x = gbuild$layout$panel_params[[1]]$y.range)
450
-    y_min <- max(gbuild$layout$panel_params[[1]]$y.range) + 0.01 * y_range
451
-    y_max <- y_min + 0.03 * y_range
446
+    if(plot_sidebar)
447
+    {
448
+        gbuild <- ggplot_build(plot = g)
449
+        y_range <- diff(x = gbuild$layout$panel_params[[1]]$y.range)
450
+        y_min <- max(gbuild$layout$panel_params[[1]]$y.range) + 0.01 * y_range
451
+        y_max <- y_min + 0.03 * y_range
452 452
 
453
-    g <- g + geom_point(aes_string(color='group'), alpha=0, shape=15, size=5) +
454
-      scale_color_manual(name=group, values=color_legend)  +
455
-      guides(color=guide_legend(override.aes=list(alpha=1))) +
456
-      annotation_raster(t(color_legend[as.character(group_genes)]), -Inf, Inf, ymin=y_min, ymax=y_max) +
457
-      coord_cartesian(ylim = c(0, y_max), clip='off')
458
-  }
453
+        g <- g + geom_point(aes_string(color='group'), alpha=0, shape=15, size=5) +
454
+            scale_color_manual(name=group, values=color_legend)  +
455
+            guides(color=guide_legend(override.aes=list(alpha=1))) +
456
+            annotation_raster(t(color_legend[as.character(group_genes)]), -Inf, Inf, ymin=y_min, ymax=y_max) +
457
+            coord_cartesian(ylim = c(0, y_max), clip='off')
458
+    }
459 459
 
460
-  return(g)
460
+    return(g)
461 461
 }
Browse code

Initial commit

liuwd15 authored on 14/09/2020 10:09:23
Showing 1 changed files
1 1
new file mode 100644
... ...
@@ -0,0 +1,461 @@
1
+#' Line plot for expression traces
2
+#'
3
+#' Plot expression traces for genes across sections in a \code{SummarizedExperiment} object.
4
+#'
5
+#' @param object A \code{SummarizedExperiment} object.
6
+#' @param genes A character vector of gene names for plotting expression traces.
7
+#' @param matrix Character, must be one of \code{"count"}, \code{"normalized"}, or \code{"scaled"}.
8
+#' @param facet Logical. Plot the expression trace of each gene in a facet if it is \code{TRUE}.
9
+#' @param span Numeric, the amount of smoothing for the default loess smoother. Smaller numbers produce wigglier lines, larger numbers produce smoother lines. Set it to 0 for non-smoothing lines.
10
+#'
11
+#' @return A \code{ggplot} object.
12
+#'
13
+#' @seealso \code{geom_smooth} for plotting smooth lines, \code{facet_wrap} for faceting genes.
14
+#'
15
+#' @importFrom ggplot2 ggplot aes_string theme theme_bw geom_smooth geom_line labs facet_wrap
16
+#' @importFrom SummarizedExperiment rowData colData
17
+#' @export
18
+#'
19
+#' @examples
20
+#' data(zh.data)
21
+#' zh <- CreateTomo(zh.data)
22
+#' zh <- Normalize(zh)
23
+#' LinePlot(zh,
24
+#'  c("ENSDARG00000002131", "ENSDARG00000003061", "ENSDARG00000076075", "ENSDARG00000076850"))
25
+#'
26
+#' # Do not smooth lines.
27
+#' LinePlot(zh,
28
+#'  c("ENSDARG00000002131", "ENSDARG00000003061", "ENSDARG00000076075", "ENSDARG00000076850"), span=0)
29
+#'
30
+#' # Plot genes in different facets.
31
+#' LinePlot(zh,
32
+#'  c("ENSDARG00000002131", "ENSDARG00000003061", "ENSDARG00000076075", "ENSDARG00000076850"),
33
+#'  facet=TRUE)
34
+LinePlot <- function(object, genes, matrix='normalized', facet=FALSE, span=0.3)
35
+{
36
+  exp_matrix <- GetMatrix(object, matrix)
37
+  exp_genes <- subset(exp_matrix, rowData(object)$gene %in% genes)
38
+  exp_genes_df <- NULL
39
+  sections <- factor(colData(object)$section, levels=colData(object)$section)
40
+  for(i in 1:nrow(exp_genes))
41
+  {
42
+    exp_gene_df <- data.frame(gene=genes[i], section=sections, value=exp_genes[i,])
43
+    exp_genes_df <- rbind(exp_genes_df, exp_gene_df)
44
+  }
45
+  ylab <- c('Count', 'Normalized count', 'Scaled expression')[c('count', 'normalized', 'scaled') == matrix]
46
+
47
+  g <- ggplot(exp_genes_df, aes_string(x='section', y='value', group='gene', color='gene')) +
48
+    labs(x='', y=ylab) +
49
+    theme_bw() +
50
+    theme(legend.position='top', axis.text.x=element_text(angle=90, hjust=1, vjust=0.5))
51
+
52
+  if(span > 0)
53
+    g <- g + geom_smooth(method='loess', span = span, se=F)
54
+  else
55
+    g <- g + geom_line(size=1)
56
+
57
+  if(facet)
58
+    g <- g + facet_wrap(~gene, scales = 'free') + theme(legend.position='none')
59
+
60
+  return(g)
61
+}
62
+
63
+#' Embedding plot for sections
64
+#'
65
+#' Scatter plot for sections with two-dimenstional embeddings in a \code{SummarizedExperiment} object. Each point stands for a section.
66
+#'
67
+#' @param object A \code{SummarizedExperiment} object.
68
+#' @param group Character, a variable in slot \code{meta} defining the groups of sections. Sections in the same group have same colors.
69
+#' @param method Character, the embeddings for scatter plot. Must be one of \code{"TSNE"}, \code{"UMAP"}, or \code{"PCA"}.
70
+#'
71
+#' @return A \code{ggplot} object.
72
+#'
73
+#' @importFrom SummarizedExperiment colData
74
+#' @importFrom ggplot2 ggplot aes_string geom_point theme theme_bw
75
+#' @importFrom ggrepel geom_text_repel
76
+#'
77
+#' @export
78
+#'
79
+#' @examples
80
+#' data(zh.data)
81
+#' zh <- CreateTomo(zh.data)
82
+#' zh <- Normalize(zh)
83
+#' zh <- TSNE(zh)
84
+#' # Plot TSNE embeddings.
85
+#' EmbedPlot(zh)
86
+#'
87
+#' # Plot UMAP embeddings.
88
+#' zh <- UMAP(zh)
89
+#' EmbedPlot(zh, method="UMAP")
90
+#'
91
+#' # Color sections by kmeans cluster labels.
92
+#' zh <- KMeans(zh, 3)
93
+#' EmbedPlot(zh, group="kmeans_cluster")
94
+EmbedPlot <- function(object, group='section', method='TSNE')
95
+{
96
+  if(!group %in% names(colData(object)))
97
+  {
98
+    group <- 'section'
99
+    cat("Unknown parameter 'group' for plotting sections!")
100
+  }
101
+  colData(object)[[group]] <- as.character(colData(object)[[group]])
102
+
103
+  if(method == 'TSNE')
104
+  {
105
+    if(all(c('TSNE1','TSNE2') %in% names(colData(object))))
106
+    {
107
+      g <- ggplot(data.frame(colData(object)), aes_string(x='TSNE1', y='TSNE2', color=group)) +
108
+        geom_point() +
109
+        geom_text_repel(aes_string(label='section')) +
110
+        theme_bw() +
111
+        theme(legend.position='none')
112
+      return(g)
113
+    }
114
+    else
115
+    {
116
+      cat("Function 'TSNE' must be run before plotting TSNE embeddings.\n")
117
+    }
118
+  }
119
+  else if(method == 'UMAP')
120
+  {
121
+    if(all(c('UMAP1','UMAP2') %in% names(colData(object))))
122
+    {
123
+      g <- ggplot(data.frame(colData(object)), aes_string(x='UMAP1', y='UMAP2', color=group)) +
124
+        geom_point(aes_string()) +
125
+        geom_text_repel(aes_string(label='section')) +
126
+        theme_bw() +
127
+        theme(legend.position='none')
128
+      return(g)
129
+    }
130
+    else
131
+    {
132
+      cat("Function 'UMAP' must be run before plotting UMAP embeddings.\n")
133
+    }
134
+  }
135
+  else if(method == 'PCA')
136
+  {
137
+    if(all(c('PC1','PC2') %in% names(colData(object))))
138
+    {
139
+      g <- ggplot(data.frame(colData(object)), aes_string(x='PC1', y='PC2', color=group)) +
140
+        geom_point(aes_string()) +
141
+        geom_text_repel(aes_string(label='section')) +
142
+        theme_bw() +
143
+        theme(legend.position='none')
144
+      return(g)
145
+    }
146
+    else
147
+    {
148
+      cat("Function 'PCA' must be run before plotting PC embeddings.\n")
149
+    }
150
+  }
151
+  else
152
+  {
153
+    cat("Unknown embeddings!\n")
154
+  }
155
+}
156
+
157
+#' Embedding plot for genes
158
+#'
159
+#' Scatter plot for genes with two-dimenstional embeddings in a \code{SummarizedExperiment} object. Each point stands for a gene.
160
+#'
161
+#' @param object A \code{SummarizedExperiment} object.
162
+#' @param gene.df Data.frame. The first column must be a vector of gene names, and has the name \code{"gene"}. Additional columns in \code{gene.df} can be used to set the colors of genes.
163
+#' @param group Character, a column name in \code{gene.df} defining the groups of genes. Genes in the same group have same colors.
164
+#' @param method Character, the embeddings for scatter plot. Must be one of \code{"TSNE"}, \code{"UMAP"}, or \code{"PCA"}.
165
+#'
166
+#' @return A \code{ggplot} object.
167
+#'
168
+#' @importFrom SummarizedExperiment assays rowData
169
+#' @importFrom ggplot2 ggplot aes_string geom_point theme theme_bw
170
+#'
171
+#' @export
172
+#'
173
+#' @examples
174
+#' data(zh.data)
175
+#' zh <- CreateTomo(zh.data)
176
+#' zh <- Normalize(zh)
177
+#' peak_genes <- FindPeakGene(zh)
178
+#' zh <- TSNE(zh, peak_genes$gene)
179
+#' # Color genes by peak centers.
180
+#' GeneEmbedPlot(zh, peak_genes)
181
+#'
182
+#' # Color genes by peak starts.
183
+#' GeneEmbedPlot(zh, peak_genes, group="start")
184
+#'
185
+#' # Do not color genes.
186
+#' GeneEmbedPlot(zh, peak_genes["gene"])
187
+GeneEmbedPlot <- function(object, gene.df, group='center', method='TSNE')
188
+{
189
+  if(group %in% names(gene.df))
190
+    gene.df[[group]] <- as.factor(gene.df[[group]])
191
+  if(method == 'TSNE')
192
+  {
193
+    if(all(c('TSNE1','TSNE2') %in% names(rowData(object))))
194
+    {
195
+      gene.df$TSNE1 <- rowData(object)[as.character(gene.df$gene), 'TSNE1']
196
+      gene.df$TSNE2 <- rowData(object)[as.character(gene.df$gene), 'TSNE2']
197
+      g <- ggplot(gene.df, aes_string(x='TSNE1', y='TSNE2')) +
198
+        theme_bw()
199
+
200
+      if(group %in% names(gene.df))
201
+        g <- g + geom_point(aes_string(color=group))
202
+      else
203
+        g <- g + geom_point()
204
+
205
+      return(g)
206
+    }
207
+    else
208
+    {
209
+      cat("Function 'TSNE' must be run for input genes before plotting TSNE embeddings.\n")
210
+    }
211
+  }
212
+  else if(method == 'UMAP')
213
+  {
214
+    if(all(c('UMAP1','UMAP2') %in% names(rowData(object))))
215
+    {
216
+      gene.df$UMAP1 <- rowData(object)[as.character(gene.df$gene), 'UMAP1']
217
+      gene.df$UMAP2 <- rowData(object)[as.character(gene.df$gene), 'UMAP2']
218
+      g <- ggplot(gene.df, aes_string(x='UMAP1', y='UMAP2')) +
219
+        theme_bw()
220
+
221
+      if(group %in% names(gene.df))
222
+        g <- g + geom_point(aes_string(color=group))
223
+      else
224
+        g <- g + geom_point()
225
+
226
+      return(g)
227
+    }
228
+    else
229
+    {
230
+      cat("Function 'UMAP' must be run for input genes before plotting UMAP embeddings.\n")
231
+    }
232
+  }
233
+  else if(method == 'PCA')
234
+  {
235
+    if(all(c('PC1','PC2') %in% names(rowData(object))))
236
+    {
237
+      gene.df$PC1 <- rowData(object)[as.character(gene.df$gene), 'PC1']
238
+      gene.df$PC2 <- rowData(object)[as.character(gene.df$gene), 'PC2']
239
+      g <- ggplot(gene.df, aes_string(x='PC1', y='PC2')) +
240
+        theme_bw()
241
+
242
+      if(group %in% names(gene.df))
243
+        g <- g + geom_point(aes_string(color=group))
244
+      else
245
+        g <- g + geom_point()
246
+
247
+      return(g)
248
+    }
249
+    else
250
+    {
251
+      cat("Function 'PCA' must be run for input genes before plotting PC embeddings.\n")
252
+    }
253
+  }
254
+  else
255
+  {
256
+    cat("Unknown embeddings!\n")
257
+  }
258
+}
259
+
260
+#' Expression heatmap
261
+#'
262
+#' Heatmap for gene expression across sections in a \code{SummarizedExperiment} object.
263
+#'
264
+#' @param object A \code{SummarizedExperiment} object.
265
+#' @param genes A vector of character, the name of genes to plot heatmap.
266
+#' @param matrix Character, must be one of \code{"count"}, \code{"normalized"}, or \code{"scaled"}.
267
+#' @param size Character, the size of gene names. Set it to 0 if you do not want to show gene names.
268
+#'
269
+#' @return A \code{ggplot} object.
270
+#'
271
+#' @importFrom reshape2 melt
272
+#' @importFrom ggplot2 ggplot aes_string geom_raster scale_fill_gradient2 theme element_blank element_text
273
+#'
274
+#' @export
275
+#'
276
+#' @examples
277
+#' data(zh.data)
278
+#' zh <- CreateTomo(zh.data)
279
+#' zh <- Normalize(zh)
280
+#' zh <- Scale(zh)
281
+#'
282
+#' # Plot some genes.
283
+#' ExpHeatmap(zh,
284
+#'  c("ENSDARG00000002131", "ENSDARG00000003061", "ENSDARG00000076075", "ENSDARG00000076850"))
285
+#'
286
+#' # Plot peak genes.
287
+#' peak_genes <- FindPeakGene(zh)
288
+#' ExpHeatmap(zh, peak_genes$gene)
289
+#'
290
+#' # Remove gene names if too many genes are in the heatmap.
291
+#' ExpHeatmap(zh, peak_genes$gene, size=0)
292
+ExpHeatmap <- function(object, genes, matrix='scaled', size=5)
293
+{
294
+  genes <- rev(as.character(genes))
295
+  exp_matrix <- GetMatrix(object, matrix)[genes, ]
296
+  if(matrix == 'normalized')
297
+    exp_matrix <- log10(exp_matrix + 1)
298
+  genes_df <- melt(exp_matrix, varnames=c('gene','section'))
299
+  exp_name <-  c('Log10(Count+1)', 'Log10(Normalized count+1)', 'Scaled expression\n(Z score)')[c('count', 'normalized', 'scaled') == matrix]
300
+
301
+  g <- ggplot(genes_df, aes_string(x='section', y='gene', fill='value')) +
302
+    geom_raster() +
303
+    scale_fill_gradient2(low='magenta', mid='black', high='yellow', name=exp_name) +
304
+    theme(axis.line=element_blank(),
305
+          axis.text.x=element_text(angle=90, hjust=1, vjust=0.5),
306
+          axis.ticks=element_blank(),
307
+          panel.background=element_blank())
308
+  if(size == 0)
309
+    g <- g + theme(axis.text.y=element_blank())
310
+  else
311
+    g <- g + theme(axis.text.y=element_text(size=size))
312
+
313
+  return(g)
314
+}
315
+
316
+#' Correlation heatmap of sections
317
+#'
318
+#' Heatmap pf correlation coefficients between any two sections in a \code{SummarizedExperiment} object.
319
+#'
320
+#' @param object A \code{SummarizedExperiment} object.
321
+#' @param matrix Character, must be one of \code{"count"}, \code{"normalized"}, or \code{"scaled"}.
322
+#' @param max.cor Numeric, correlation coefficients bigger than \code{max.cor} are set to \code{max.cor}. It is used to clearly show small correlation coefficients.
323
+#' @param cor.method Character, the method to calculate correlation coefficients. must be one of \code{"pearson"}, \code{"kendall"}, or \code{"spearman"}.
324
+#'
325
+#' @return A \code{ggplot} object.
326
+#'
327
+#' @importFrom stats cor
328
+#' @importFrom reshape2 melt
329
+#' @importFrom ggplot2 ggplot aes_string geom_raster scale_fill_gradientn labs theme element_blank element_text
330
+#' @importFrom colorRamps rgb.tables
331
+#'
332
+#' @export
333
+#'
334
+#' @examples
335
+#' data(zh.data)
336
+#' zh <- CreateTomo(zh.data)
337
+#' zh <- Normalize(zh)
338
+#' zh <- Scale(zh)
339
+#' CorHeatmap(zh)
340
+#'
341
+#' # Use Spearman correlation coefficients.
342
+#' CorHeatmap(zh, cor.method='spearman')
343
+#'
344
+#' # Set max correlation coefficients to 0.3.
345
+#' CorHeatmap(zh, max.cor=0.3)
346
+CorHeatmap <- function(object, matrix='scaled', max.cor=0.5, cor.method='pearson')
347
+{
348
+  exp_matrix <- GetMatrix(object, matrix=matrix)
349
+  cor_matrix <- cor(exp_matrix, method=cor.method)
350
+  cor_matrix[cor_matrix > max.cor] <- max.cor
351
+  cor_df <- melt(cor_matrix, varnames=c('section1','section2'))
352
+  # c('Pearson r', 'Kendall τ', 'Spearman ρ')
353
+  cor_name <- expression("Pearson"~r,"Kendall"~tau,"Spearman"~rho)[c('pearson', 'kendall', 'spearman') == cor.method]
354
+
355
+  g <- ggplot(cor_df, aes_string(x='section1', y='section2', fill='value')) +
356
+    geom_raster() +
357
+    scale_fill_gradientn(colors=rgb.tables(10), name=cor_name) +
358
+    labs(x='', y='') +
359
+    theme(axis.line=element_blank(),
360
+          axis.text.x=element_text(angle=90, hjust=1, vjust=0.5),
361
+          axis.ticks=element_blank(),
362
+          panel.background=element_blank())
363
+
364
+  return(g)
365
+}
366
+
367
+#' Correlation heatmap of genes
368
+#'
369
+#' Heatmap of correlation coefficients between any two queried genes in a \code{SummarizedExperiment} object.
370
+#'
371
+#' @param object A \code{SummarizedExperiment} object.
372
+#' @param gene.df Data.frame. The first column must be a vector of gene names, and has the name \code{"gene"}. Additional columns in \code{gene.df} can be used to set the colors of genes.
373
+#' @param group Character, a column name in \code{gene.df} defining the groups of genes. Genes in the same group have same colors on the side bar.
374
+#' @param matrix Character, must be one of \code{"count"}, \code{"normalized"}, or \code{"scaled"}.
375
+#' @param size Numeric, the size of gene names. Set it to 0 if you do not want to show gene names.
376
+#' @param cor.method Character, the method to calculate correlation coefficients. must be one of \code{"pearson"}, \code{"kendall"}, or \code{"spearman"}.
377
+#'
378
+#' @details This method can create a pure heatmap or a heatmap with side bar. If you prefer a pure heatmap, input a \code{gene.df} with a single column of gene names.
379
+#' However, you may want to show additional information of genes with a side bar, and the grouping information should be saved as additional column(s) of \code{gene.df}, and declared as \code{group}.
380
+#' By default, you can use the output by \code{FindPeakGene} as input \code{gene.df}. Peak genes will be grouped by their centers on the side bar.
381
+#'
382
+#' @return A \code{ggplot} object.
383
+#'
384
+#' @importFrom stats cor
385
+#' @importFrom grDevices rainbow
386
+#' @importFrom reshape2 melt
387
+#' @importFrom ggplot2 ggplot aes_string geom_raster scale_fill_gradientn labs theme element_blank element_text ggplot_build scale_color_manual guides annotation_raster guide_legend coord_cartesian
388
+#'
389
+#' @export
390
+#'
391
+#' @examples
392
+#' data(zh.data)
393
+#' zh <- CreateTomo(zh.data)
394
+#' zh <- Normalize(zh)
395
+#' zh <- Scale(zh)
396
+#'
397
+#' # Correlation heatmap for all peak genes.
398
+#' peak_genes <- FindPeakGene(zh)
399
+#' GeneCorHeatmap(zh, peak_genes)
400
+#'
401
+#' # Use Spearman correlation coefficients.
402
+#' GeneCorHeatmap(zh, peak_genes, cor.method="spearman")
403
+#'
404
+#' # Group genes by peak start.
405
+#' GeneCorHeatmap(zh, peak_genes, group="start")
406
+#'
407
+#' # Plot without side bar.
408
+#' GeneCorHeatmap(zh, data.frame(
409
+#'  gene=c("ENSDARG00000002131", "ENSDARG00000003061", "ENSDARG00000076075", "ENSDARG00000076850")))
410
+GeneCorHeatmap <- function(object, gene.df, group='center', matrix='scaled', size=5, cor.method='pearson')
411
+{
412
+  exp_matrix <- GetMatrix(object, matrix=matrix)
413
+  genes <- as.character(gene.df$gene)
414
+  cor_matrix <- cor(t(exp_matrix[genes, ]), method=cor.method)
415
+  cor_df <- melt(cor_matrix, varnames=c('gene1', 'gene2'))
416
+  cor_name <- expression("Pearson"~r,"Kendall"~tau,"Spearman"~rho)[c('pearson', 'kendall', 'spearman') == cor.method]
417
+  # plot sidebar if genes are grouped
418
+  plot_sidebar <- group %in% names(gene.df)
419
+
420
+  if(plot_sidebar)
421
+  {
422
+    group_genes <- gene.df[[group]]
423
+    names(group_genes) <- gene.df$gene
424
+    cor_df$group <- as.factor(group_genes[cor_df$gene1])
425
+    # Colors for column side bar
426
+    n_groups <- length(unique(group_genes))
427
+    color_legend <- rainbow(n_groups)
428
+    names(color_legend) <- unique(group_genes)
429
+  }
430
+
431
+  g <- ggplot(cor_df, aes_string(x='gene1', y='gene2', fill='value')) +
432
+    geom_raster() +
433
+    scale_fill_gradientn(colors=rgb.tables(10), name=cor_name)  +
434
+    labs(x='', y='') +
435
+    theme(axis.line=element_blank(),
436
+          axis.text.y=element_blank(),
437
+          axis.ticks=element_blank(),
438
+          panel.background=element_blank(),
439
+          legend.key=element_blank(),
440
+          legend.box='horizontal')
441
+  if(size == 0)
442
+    g <- g + theme(axis.text.x=element_blank())
443
+  else
444
+    g <- g + theme(axis.text.x=element_text(angle=90, hjust=1, size=size))
445
+
446
+  if(plot_sidebar)
447
+  {
448
+    gbuild <- ggplot_build(plot = g)
449
+    y_range <- diff(x = gbuild$layout$panel_params[[1]]$y.range)
450
+    y_min <- max(gbuild$layout$panel_params[[1]]$y.range) + 0.01 * y_range
451
+    y_max <- y_min + 0.03 * y_range
452
+
453
+    g <- g + geom_point(aes_string(color='group'), alpha=0, shape=15, size=5) +
454
+      scale_color_manual(name=group, values=color_legend)  +
455
+      guides(color=guide_legend(override.aes=list(alpha=1))) +
456
+      annotation_raster(t(color_legend[as.character(group_genes)]), -Inf, Inf, ymin=y_min, ymax=y_max) +
457
+      coord_cartesian(ylim = c(0, y_max), clip='off')
458
+  }
459
+
460
+  return(g)
461
+}