转录组不求人系列(十): NCS级别的火山图，总有一款适合你！

1640738509(1).png

转录组差异基因分析完之后就可以进行后续的分析了，但在此之前重要的是对差异结果的可视化，可视化办法有很多，例如热图，我们之前已经详细的说过如何做了，这一节着重说一下火山图的制作来表示差异基因的结果。

我们提供三种火山图的制作方法，每种做法都是CNS文章中出现的，争取一步到位，做的图就是发文章的图，让自己的可视化更好。

image

绘制火山图需要的文件是差异分析结果，包含FC和P值：

image

第一种

第一种方法比较简单，读入文件，横轴为log2FC，纵轴为-Log10Pvalue，直接用plot函数画图即可；

setwd("E:/生物信息学/转录组火山图")
df <- read.csv("DEGs_trans.csv",row.names = 1) #导入数据，第一列作为行名
plot(df$log2FoldChange, -log10(df$padj), 
     col="#00000033", pch=19,cex=2.5,
     xlab=paste("log2 (fold change)"),
     ylab="-log10 (P_value)")+title("DEGs")

image

初步画好了散点图，接着加入阈值条件，显示差异基因：

#筛选上下调
up <- subset(df, df$padj < 0.01 & df$log2FoldChange > 5)
down <- subset(df, df$padj < 0.01 & df$log2FoldChange < as.numeric(-5))
#绘制上下调
points(up$log2FoldChange, -log10(up$padj), col=1, bg = brewer.pal(9, "YlOrRd")[6], pch=21, cex=2.5)
points(down$log2FoldChange, -log10(down$padj), col = 1, bg = brewer.pal(11,"RdBu")[9], pch = 21,cex=2.5)
#加上阈值线
abline(h=-log10(0.01),v=c(-5,5),lty=2,lwd=1)

image

这样一副火山图就基本完成了，还有一个最后的步骤就是显示感兴趣的基因名称，plot画图可以用text函数添加，只需要找到基因的坐标即可。

text(-5.301223,-log10(3.260000e-64),"OGFOD2")

image

其他更多的精彩内容请移步访问我的个人公众号，谢谢支持！