1640738509(1).png
转录组差异基因分析完之后就可以进行后续的分析了,但在此之前重要的是对差异结果的可视化,可视化办法有很多,例如热图,我们之前已经详细的说过如何做了,这一节着重说一下火山图的制作来表示差异基因的结果。
我们提供三种火山图的制作方法,每种做法都是CNS文章中出现的,争取一步到位,做的图就是发文章的图,让自己的可视化更好。
绘制火山图需要的文件是差异分析结果,包含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")
其他更多的精彩内容请移步访问我的个人公众号,谢谢支持!