|
單細(xì)胞常見的可視化方式有DimPlot，F(xiàn)eaturePlot ，DotPlot ，VlnPlot 和 DoHeatmap幾種 ，Seurat均可以實(shí)現(xiàn)，但文獻(xiàn)中的圖大多會精美很多。比如 驚艷umap圖: scRNA復(fù)現(xiàn)|所見即所得，和Cell學(xué)umap，plot1cell完成驚艷的細(xì)胞注釋umap圖； DimPlot美化 scRNA分析 | 定制 美化FeaturePlot 圖，你需要的都在這, DotPlot美化scRNA分析| 和SCI學(xué) 定制化聚類點(diǎn)圖（Dotplot ），含二行代碼出圖方式， DoHeatmap 熱圖：scRNA分析| DoHeatmap 美化，dittoSeq ，scillus 一行代碼出圖，你PICK誰？ 本次介紹Seurat 以及 ggplot2繪制，優(yōu)化堆疊小提琴圖的方法。 仍然使用之前注釋過的sce.anno.RData數(shù)據(jù) ，后臺回復(fù) anno 即可獲取。 library(Seurat)library(tidyverse)
load("sce.anno.RData")head(sce2,2)
1，基礎(chǔ)VlnPlot圖
首先計(jì)算marker基因，然后使用seurat的DoHeatmap 函數(shù)繪制初始熱圖
all_markers <- FindAllMarkers(object = sce2)top5 <- all_markers %>% group_by(cluster) %>% top_n(5, avg_log2FC)###少量基因VlnPlot(sce2, features = c("CD3D","SPP1"))### 所有marker 基因VlnPlot(sce2, features = top5$gene)
當(dāng)展示少量基因時(shí)候，很清晰 。但是更常見的時(shí)候需要同時(shí)展示各個(gè)cluster/celltype的marker gene ，這時(shí)候就會看不清晰。
2，Seurat-堆疊VlnPlot圖Seurat的VlnPlot函數(shù)中stack 參數(shù)可以實(shí)現(xiàn)堆疊小提琴圖，flip 是否翻轉(zhuǎn) #Seurat 的stack 函數(shù)a <- VlnPlot(sce2, features = top5$gene, stack = TRUE, sort = TRUE) + theme(legend.position = "none") + ggtitle("Identity on y-axis")# flip 翻轉(zhuǎn)b <- VlnPlot(sce2, features = top5$gene, stack = TRUE, sort = TRUE, flip = TRUE) + theme(legend.position = "none") + ggtitle("Identity on x-axis")
a + b
3，Seurat-優(yōu)化顏色，大小，方向自定義顏色，是否排序，主題等信息更是和前面的一樣，直接添加theme信息即可。 注意如果想要每種cluster/celltype是一種顏色的話使用split.by參數(shù)。
my36colors <-c('#E5D2DD', '#53A85F', '#F1BB72', '#F3B1A0', '#D6E7A3', '#57C3F3', '#476D87', '#E95C59', '#E59CC4', '#AB3282', '#23452F', '#BD956A', '#8C549C', '#585658', '#9FA3A8', '#E0D4CA', '#5F3D69', '#C5DEBA', '#58A4C3', '#E4C755', '#F7F398', '#AA9A59', '#E63863', '#E39A35', '#C1E6F3', '#6778AE', '#91D0BE', '#B53E2B', '#712820', '#DCC1DD', '#CCE0F5', '#CCC9E6', '#625D9E', '#68A180', '#3A6963', '#968175')
VlnPlot(sce2, features = top_marker$gene, stack = TRUE, sort = TRUE, cols = my36colors, split.by = "celltype" , #每種cluster 一個(gè)顏色 flip = TRUE) + theme(legend.position = "none") + ggtitle("Identity on x-axis")
Seurat的堆疊小提琴圖其實(shí)已經(jīng)可以了，當(dāng)然也可以使用ggplot2進(jìn)行更多的自定義。
1，提取，轉(zhuǎn)化數(shù)據(jù)首先使用FetchData提取出marker gene的表達(dá)量，celltype /seurat_clusters（寬數(shù)據(jù)），然后轉(zhuǎn)為ggplot2讀取的長數(shù)據(jù)類型 。 此外對照上述的圖，可以看到celltype /seurat_clusters一個(gè)表達(dá)量值，而FetchData得到的是每個(gè)cell 的表達(dá)量，因此還需要計(jì)算每種cluster的基因均值。 vln.dat=FetchData(sce2,c(top_marker$gene,"celltype","seurat_clusters"))
vln.dat$Cell <- rownames(vln.dat)#寬轉(zhuǎn)長vln.dat.melt <- reshape2::melt(vln.dat, id.vars = c("Cell","seurat_clusters"), measure.vars = top_marker$gene, variable.name = "gene", value.name = "Expr") %>% group_by(seurat_clusters,gene) %>% #分組 mutate(fillcolor=mean(Expr)) #計(jì)算均值
2，ggplot2 繪制-核心 ggplot(vln.dat.melt, aes(factor(seurat_clusters), Expr, fill = gene)) + geom_violin(scale = "width", adjust = 1, trim = TRUE) + facet_grid(rows = vars(gene), scales = "free", switch = "y")
3，ggplot2 繪制-優(yōu)化上述是ggplot2繪制堆疊小提琴圖的核心代碼，可以做很多調(diào)整 （1）主題（大小，顏色），legend 等 （2）“翻轉(zhuǎn)”（使用aes調(diào)整橫縱坐標(biāo)） p1 <- ggplot(vln.dat.melt, aes(gene, Expr, fill = gene)) + geom_violin(scale = "width", adjust = 1, trim = TRUE) + scale_y_continuous(expand = c(0, 0), position="right", labels = function(x) c(rep(x = "", times = length(x)-2), x[length(x) - 1], "")) + facet_grid(rows = vars(seurat_clusters), scales = "free", switch = "y") + scale_fill_manual(values = my36colors) + theme_cowplot(font_size = 12) + theme(legend.position = "none", panel.spacing = unit(0, "lines"), plot.title = element_text(hjust = 0.5), panel.background = element_rect(fill = NA, color = "black"), plot.margin = margin(7, 7, 0, 7, "pt"), strip.background = element_blank(), strip.text = element_text(face = "bold"), strip.text.y.left = element_text(angle = 0), axis.title.x = element_blank(), axis.ticks.x = element_blank(), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, color = "black") ) + ggtitle("Feature on x-axis with annotation") + ylab("Expression Level")p1
（3）添加基因的分組/注釋 A：添加分組，注釋
假設(shè)知道m(xù)arker gene的通路，也可以添加上（為了美觀先隱藏p1中的橫坐標(biāo)基因標(biāo)簽）
#隱藏axis.text.xp2 <- ggplot(vln.dat.melt, aes(gene, Expr, fill = gene)) + geom_violin(scale = "width", adjust = 1, trim = TRUE) + scale_y_continuous(expand = c(0, 0), position="right", labels = function(x) c(rep(x = "", times = length(x)-2), x[length(x) - 1], "")) + facet_grid(rows = vars(seurat_clusters), scales = "free", switch = "y") + scale_fill_manual(values = my36colors) + theme_cowplot(font_size = 12) + theme(legend.position = "none", panel.spacing = unit(0, "lines"), plot.title = element_text(hjust = 0.5), panel.background = element_rect(fill = NA, color = "black"), plot.margin = margin(7, 7, 0, 7, "pt"), strip.background = element_blank(), strip.text = element_text(face = "bold"), strip.text.y.left = element_text(angle = 0), axis.title.x = element_blank(), axis.ticks.x = element_blank(), axis.text.x = element_blank() #隱藏 ) + ggtitle("Feature on x-axis with annotation") + ylab("Expression Level")p2
B：構(gòu)建注釋信息-基因分組信息
這里通路是隨便寫的，僅為示例，并不是該marker gene 在的通路。 # Create grouping infodf <- data.frame(x = levels(vln.dat.melt$gene), group = c("A","A","B","B","B","B","B","C","C","C","D","D","D", "D","D","D","D","D"), stringsAsFactors = FALSE)df$x <- factor(df$x, levels = levels(vln.dat.melt$gene))df$group <- factor(df$group)#可以修改 注釋 展示的名字levels(df$group) = c("ECM-receptor interaction", "PI3K-Akt signaling pathway", "MAPK signaling pathway", "Cell adhesion molecules")#設(shè)置顏色color <- c("cyan", "pink", "green", "darkorange")
# guides() is used to specify some aesthetic parameters of legend keyp3 <- ggplot(df, aes(x = x, y = 1, fill = group)) + geom_tile() + theme_bw(base_size = 12) + scale_fill_manual(values = my36colors) + scale_y_continuous(expand = c(0, 0)) + guides(fill = guide_legend(direction = "vertical", label.position = "right", title.theme = element_blank(), keyheight = 0.5, nrow = 2)) + theme(legend.position = "bottom", legend.justification = "left", legend.margin = margin(0,0,0,0), legend.box.margin = margin(-10,5,0,0), panel.spacing = unit(0, "lines"), panel.background = element_blank(), panel.border = element_blank(), plot.background = element_blank(), plot.margin = margin(0, 7, 7, 7, "pt"), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, color = "black"), axis.title.y = element_blank(), axis.ticks.y = element_blank(), axis.text.y = element_blank()) + xlab("Feature")p3
C：拼圖收工 # Use plot_grid to join plotsplot_grid(p2, p3, ncol = 1, rel_heights = c(0.78, 0.22), align = "v", axis = "lr")
參考資料： https://github.com/ycl6/StackedVlnPlot 精心整理（含圖PLUS版）|R語言生信分析，可視化（R統(tǒng)計(jì)，ggplot2繪圖，生信圖形可視化匯總） RNAseq純生信挖掘思路分享？不，主要是送你代碼?。ńㄗh收藏）
|
