当有了聚类结果(cluster)或注释结果(celltype)后就可以 找不同cluster/celltype间,不同样本间 或者 不同分组间的差异,为后面的 机制探索 或者 样本间/组间异质性研究 提供一些帮助。
但是在对差异结果进行可视化的时候,你是否经常会遇到以下情况,
1)单细胞数据的火山图 中间经常会缺一块 ?Y轴最高处会堆积很多点?
2)除bulk常见的火山图外,文献中是否还有什么其他的展示形式?
3)可以同时展示单细胞各个cluster /celltype的差异情况吗?
一 差异marker分析
Seurat可以通过FindMarkers函数 和 FindAllMarkers函数寻找不同cluster的差异表达基因。
1 FindMarkers函数
1.1 特定cluster/celltype 与其他的差异
是one-others的差异分析方法,由ident.1来制定cluster/celltype ,然今后计算ident.1与其他所有的差异。本例就是Epi与其余的所有celltype找差异marker。
代码语言:javascript复制library(Seurat)
library(tidyverse)
#载入数据
load("sce.anno.RData")
#设置颜色
colour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080","#1E90FF","#7CFC00","#FFFF00",
"#808000","#FF00FF","#FA8072","#7B68EE","#9400D3","#800080","#A0522D","#D2B48C","#D2691E","#87CEEB","#40E0D0","#5F9EA0",
"#FF1493","#0000CD","#008B8B","#FFE4B5","#8A2BE2","#228B22","#E9967A","#4682B4","#32CD32","#F0E68C","#FFFFE0","#EE82EE",
"#FF6347","#6A5ACD","#9932CC","#8B008B","#8B4513","#DEB887")
# find all markers of cluster Epi
cluster1.markers <- FindMarkers(sce2, ident.1 = "Epi", min.pct = 0.25)
head(cluster1.markers, n = 5)
代码语言:javascript复制# p_val avg_log2FC pct.1 pct.2 p_val_adj
# MXRA8 0 -1.5057649 0.024 0.385 0
# MFAP2 0 -0.8893903 0.020 0.279 0
# C1QA 0 -3.0473042 0.048 0.333 0
# C1QC 0 -2.5871212 0.022 0.273 0
# C1QB 0 -3.0683349 0.043 0.316 0
1.2 指定cluster/celltype 之间的差异
示例为计算Epi 和 T B 细胞之间的差异
代码语言:javascript复制cluster2.markers <- FindMarkers(sce2, ident.1 = "Epi", ident.2 = c("T", "B"), min.pct = 0.25)
head(cluster2.markers, n = 5)
FindMarkers函数常用参数:
min.pct: 基因在两个cluster/celltype中任何一个被检测到的百分比阈值,默认0.1 。 logfc.threshold :默认为0.25 thresh.test:设定基因在两个cluster/celltype中基因差异表达量。 test.use :差异marker的检验方式,默认为“wilcox”,可选“bimod”,“roc”,“t”,“negbinom”,“LR”,“MAST”,“DESeq2”等; only.pos:是否仅返回阳性标记(默认为FALSE) 更多参数可以使用 ??FindMarkers 查看。
2 FindAllMarkers 函数
每个cluster/celltype 分别与其他所有cluster/celltype 进行比较,展示各个cluster/celltype 的前2个基因
代码语言:javascript复制sce2.markers <- FindAllMarkers(sce2, only.pos = TRUE, min.pct = 0.25,
logfc.threshold = 0.25)
sce2.markers %>%
group_by(cluster) %>%
top_n(n = 2, wt = avg_log2FC)
代码语言:javascript复制 p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene
1 1.88e-117 1.08 0.913 0.588 2.57e-113 0 LDHB
2 5.01e- 85 1.34 0.437 0.108 6.88e- 81 0 CCR7
3 0 5.57 0.996 0.215 0 1 S100A9
注意FindAllMarkers比FindMarkers的结果多了cluster列,结果为该cluster与其他所有cluster相比的差异基因。
FindAllMarkers 与 FindMarkers 参数类似,但是需要注意only.pos的默认值,only.pos:仅返回阳性标记(默认为TRUE)
二 ,可视化-火山图
1,常规方式绘制火山图
代码语言:javascript复制#自定义阈值
log2FC = 1
padj = 0.05
cluster1.markers$threshold="ns";
cluster1.markers[which(cluster1.markers$avg_log2FC > log2FC & cluster1.markers$p_val_adj <padj),]$threshold="up";
cluster1.markers[which(cluster1.markers$avg_log2FC < (-log2FC) & cluster1.markers$p_val_adj < padj),]$threshold="down";
cluster1.markers$threshold=factor(cluster1.markers$threshold, levels=c('down','ns','up'))
ggplot(data=cluster1.markers, aes(x=avg_log2FC, y=-log10(p_val_adj), color=threshold))
geom_point(alpha=0.8, size=0.8)
geom_vline(xintercept = c(-log2FC, log2FC), linetype=2, color="grey")
geom_hline(yintercept = -log10(padj), linetype=2, color="grey")
#labs(title= ifelse(""==title, "", paste("DEG:", title)))
xlab(bquote(Log[2]*FoldChange))
ylab(bquote(-Log[10]*italic(P.adj)) )
theme_classic(base_size = 14)
scale_color_manual('',labels=c(paste0("down(",table(cluster1.markers$threshold)[[1]],')'),'ns',
paste0("up(",table(cluster1.markers$threshold)[[3]],')' )),
values=c("blue", "grey","red" ) )
guides(color=guide_legend(override.aes = list(size=3, alpha=1)))
更多颜色,参数的设置以及如何添加感兴趣的基因详见ggplot2-plotly|让你的火山图“活”过来 和 火山图|给你geneList,帮我标到火山图上。
回答开始提出的问题1:
(1)因为单细胞自身区别于bulk数据的特异性,大概率会出现很多P值为0或者无限接近于0的基因,绘制常规火山图就会出现Y轴顶出现很多点。
(2)中间的空白是因为FindMarkers的默认参数含还有很多阈值,很多基因被过滤掉了。可以通过调整min.pct,logfc.threshold 和thresh.test等参数得到更多的基因,自行尝试比较。
2, 单细胞火山图展示方式
在生信技能树 两个不同单细胞亚群差异分析,何必一定要做火山图 中提到NC文章 《CD177 modulates the function and homeostasis of tumor-infiltrating regulatory T cells》有以下呈现方式,横坐标使用cluster/celltype间的pct差值,纵坐标使用log2FC ,P值可以使用颜色来表示。
文献阈值参数:Genes labeled have logfold change > 1, Δ Percentage Difference > 20% and adjusted P-value <0.05.
使用ggplot2 点图绘制方式
代码语言:javascript复制library(ggrepel)
cluster1.markers <- cluster1.markers %>%
mutate(Difference = pct.1 - pct.2) %>%
rownames_to_column("gene")
ggplot(cluster1.markers, aes(x=Difference, y=avg_log2FC, color = threshold))
geom_point(size=0.5)
scale_color_manual(values=c( "blue","grey","red") )
geom_label_repel(data=subset(cluster1.markers, avg_log2FC >= 1 & Difference >= 0.2 & p_val_adj <= 0.05),
aes(label=gene), #添加label
color="black", #设置label中标签的颜色
segment.colour = "black",#设置label框的颜色
label.padding = 0.1,
#max.overlaps = 200,
segment.size = 0.3, #框的大小
size=4)
geom_label_repel(data=subset(cluster1.markers, avg_log2FC <= -1 & Difference <= -0.2 & p_val_adj <= 0.05),
aes(label=gene), label.padding = 0.1,
color="black",
segment.colour = "black",
segment.size = 0.3, size=4)
geom_vline(xintercept = 0.0,linetype=2)
geom_hline(yintercept = 0,linetype=2)
theme_classic()
这样就解决了传统火山图 上限一行点的问题,然后可以根据ggplot2的一些参数进行细节的调整和美化。
此外单细胞还有一个特性就是有很多celltype ,那需要像NC文章那样每种cluster绘制一张图吗?是否有方式可以同时展示所有cluster的差异?
3 所有cluster热图
如文献A Spatiotemporal Organ-Wide Gene Expression and Cell Atlas of the Developing Human HeartsciPainter中的Fig2-D所示 ,提出了同时展示所有cluster差异基因的方式
这里介绍一种简单实现的方式,推荐老俊俊封装的scRNAtoolVis-R包,可以一行代码简单出图 jjVolcano 一行代码绘制单细胞火山图,可以去大佬的github star一下哈
支持和鼓励:
代码语言:javascript复制#install.packages('devtools')
#devtools::install_github('junjunlab/scRNAtoolVis')
library(scRNAtoolVis)
jjVolcano(diffData = sce_marker,
log2FC.cutoff = 0.25,
size = 3.5, #设置点的大小
fontface = 'italic', #设置字体形式
#aesCol = c('purple','orange'), #设置点的颜色
tile.col = colour, #设置cluster的颜色
#col.type = "adjustP", #设置矫正方式
topGeneN = 5 #设置展示topN的基因
)
还可以通过polar函数调整形状
代码语言:javascript复制jjVolcano(diffData = sce_marker,
tile.col = colour,
size = 3.5,
topGeneN = 5,
fontface = 'italic',
polar = T)
ylim(-8,10)
更多参数可以去github或者使用??jjVolcano查看 。
◆ ◆ ◆ ◆ ◆
精心整理(含图PLUS版)|R语言生信分析,可视化(R统计,ggplot2绘图,生信图形可视化汇总)