做完单细胞差异基因分析(FindMarkers/FindAllmarkers)之后,按照常规流程绘制出来的火山图看上去会很奇怪。
1、为什么火山图顶部聚集了很多基因?
这是由于单细胞有别于bulk的特性,会出现两组细胞之间的p值过于显著出现或接近0的情况,那么在取-log10之后的值会很大。
2、为什么火山图中间出现了空白
这是因为在进行FindMarkers时默认设置了一定的阈值,可以通过修改参数的阈值来修改自己的火山图(比如下边的示例代码中的min.pct和logfc.threshold参数)
markers <- FindMarkers(scRNA_tumor,ident.1 = 'acitve', only.pos = FALSE,
min.pct = 0.25, logfc.threshold = 0.25)
常规火山图的代码:
#增加一列显著基因,待会标名字
markers$sign <- ifelse(markers$p_val_adj < 0.005 & abs(markers$avg_log2FC) > 2,
rownames(markers),NA)
#当然也可以自定义的,随机的
k <- c("TP53","CD34","CD68")
markers <- markers %>%
mutate(Sign = ifelse(rownames(markers) %in% k, rownames(markers), NA))
#自定义阈值
log2FC = 0.585
padj = 0.05
colnames(markers)
library(ggplot2)
library(ggrepel) # 这个R包用于添加基因名称
ggplot(data = markers, aes(x = avg_log2FC, y = -log10(p_val_adj))) +
geom_point(alpha=0.4, size=3.5, aes(color=change)) +
#scale_color_manual(values=c("blue", "grey","red"))+
scale_color_manual('change',labels=c(paste0("down(",table(markers$change)[[1]],')'),
'ns',
paste0("up(",table(markers$change)[[3]],')' )),
values=c("blue", "grey","red" ))+
geom_vline(xintercept=c(-log2FC,log2FC),lty=4,col="black",linewidth=0.8) +
geom_hline(yintercept = -log10(padj),lty=4,col="black",linewidth=0.8) +
geom_label_repel(aes(label = sign ), fontface="bold",
color="black", box.padding=unit(0.35, "lines"),
point.padding=unit(0.5, "lines"), segment.colour = "grey50") +
theme_bw()
那么除了不画常规的火山图,还可以采用另一种方式进行绘制,在生信技能树推文中就提供了新的绘制方式。
横坐标 使用不同cluster细胞的的pct差值,纵坐标 使用log2FC ,P值 用颜色来表示。
针对自己的数据也进行了尝试
新的火山图绘制代码
# 单细胞火山图展示方式-new
# Findmarkers参数设置为0!
library(ggrepel)
markers <- markers %>%
mutate(Difference = pct.1 - pct.2) %>%
rownames_to_column("gene")
#画图
ggplot(markers, aes(x=Difference, y=avg_log2FC, color = change)) +
geom_point(size=1.2) +
scale_color_manual('change',labels=c(paste0("down(",table(markers$change)[[1]],')'),
'ns',
paste0("up(",table(markers$change)[[3]],')' )),
values=c("#0d1b46", "grey","tomato" ))+
geom_label_repel(data=subset(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.2,
segment.size = 0.3, #框的大小
size=4,
max.overlaps = 20) + # 增加 max.overlaps 参数
geom_label_repel(data=subset(markers,
avg_log2FC <= -1 & Difference <= -0.1 & p_val_adj <= 0.05),
aes(label=gene),
label.padding = 0.2,
color="black",
segment.colour = "black",
segment.size = 0.3, size=4,
max.overlaps = 20) +
geom_vline(xintercept = 0,linetype = 2) +
geom_hline(yintercept = 0,linetype = 2) +
labs(x="△ Percentage difference",y="Log-Fold change") +
theme_bw()
ggsave("volcanoplot_new.pdf",width = 9,height = 7)
致谢:感谢曾老师,小洁老师以及生信技能树团队全体成员。
注:若对内容有疑惑或者有发现明确错误的朋友,请联系后台(希望多多交流)。更多内容可关注公众号:生信方舟
- END -