RNA 6. 差异基因表达之

您所在的位置:网站首页 火山260 RNA 6. 差异基因表达之

RNA 6. 差异基因表达之

2023-03-10 20:02| 来源: 网络整理| 查看: 265

    火山图也是 RNA 数据分析中 SCI 文章中出现频率最高的图形,那么我们该怎么绘制这样一幅高质量可发表的图形呢?这期就为大家细细道来!

1. 前言

火山图是散点图的一种,它将统计测试中的统计显著性量度(如p value)和变化幅度相结合,从而能够帮助快速直观地识别那些变化幅度较大且具有统计学意义的数据点(基因等)。常应用于转录组研究,也能应用于基因组,蛋白质组,代谢组等统计数据。比如文章中出现的火山图,如下:

所以关注火山图(其它类型图也是),先理解每个点是什么(点代表基因、样品、通路或其它的,这个认识可以来自于常识,更准确的是看作者的描述),然后看横轴代表什么、纵轴代表什么,再看图例中展示的其他信息,如颜色、大小和形状分别代表什么。这些都理顺了,图理解就不难了。

2. 实例绘图1. 数据读取

我们这次使用的是基于 TCGA-COAD 表达数据,利用软件包 edgeR 计算出来的差异表达结果,整个过程就不多说了,利用其结果做例子吧,如下:

DEG= cut_off_logFC, ifelse(DEG$logFC> cut_off_logFC ,'Up','Down'),'None')

查看下差异基因个数,上调的基因为 2832 ,下调的基因个数为 1296 ,还算可以,不是很多,说明我们的阈值选择的算合理,如果基因个数过多,cutoff 可以设严格一些,如果过少就设置低一些,根据情况自己筛选即可,如下:

table(DEG$Sig)## ## Down None Up ## 1296 27649 2832

整理好差异基因表格,如下:

head(DEG)## logFC logCPM LR PValue FDR Sig## 1 -5.924830 3.20085238 1457.1137 8.159939e-319 2.592984e-314 Down## 2 -4.223787 2.59792097 1145.7948 3.679815e-251 5.846674e-247 Down## 3 -5.131288 1.71477711 1122.0643 5.289173e-246 5.602468e-242 Down## 4 -4.101967 1.51480635 1085.9423 3.752089e-238 2.980753e-234 Down## 5 -4.295806 3.39558390 1080.7407 5.067873e-237 3.220836e-233 Down## 6 -6.284258 -0.02616284 916.3497 2.739233e-201 1.450744e-197 Down

整理好表格就可以绘图了,自己写段代码,随心改还是挺方便的,如下:

library(ggplot2)ggplot(DEG, aes(x = logFC, y = -log10(PValue), colour=Sig)) + geom_point(alpha=0.4, size=3.5) + scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757"))+ # 辅助线 geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) + geom_hline(yintercept = -log10(cut_off_pvalue), lty=4,col="black",lwd=0.8) + # 坐标轴 labs(x="log2(Fold Change)", y="-log10 (P-value)")+ theme_bw()+ ggtitle("Volcano Plot")+ # 图例 theme(plot.title = element_text(hjust = 0.5), legend.position="right", legend.title = element_blank() ) 3. 软件包 EnhancedVolcano

火山图的绘制仍然使用 EnhancedVolcano 包,因为非常好用,只需要修改几个参数,比如输入矩阵的变量,x, y 变量所使用的列名称即可,输出火山图。

软件包 EnhancedVolcano 安装并加载,如下:

require(EnhancedVolcano)

绘制火山图,如下:

EnhancedVolcano(DEG, lab = rownames(DEG), labSize = 2, x = "logFC", y = "PValue", #selectLab = rownames(lrt)[1:4], xlab = bquote(~Log[2]~ "Fold Change"), ylab = bquote(~-Log[10]~italic(P)), pCutoff = cut_off_pvalue,## pvalue闃堝€? FCcutoff = cut_off_logFC,## FC cutoff xlim = c(-5,5), ylim = c(0,260), colAlpha = 0.6, legendLabels =c("NS","Log2 FC"," P-value", " P-value & Log2 FC"), legendPosition = "top", legendLabSize = 10, legendIconSize = 3.0, pointSize = 1.5, title ="Volcano Plot", subtitle = 'EnhancedVolcano') 3. 火山图解读

关于火山图的结果,内行人一看就懂,而刚接触生信的人就是有一堆的疑惑,迟迟不能自拔的感觉,我特别能理解,记得当时我接触测序学生信的时候,怎么也理解不了基因组测序整个流程,以及 Reads 会有 3' 端和 5' 端之分,哈哈,当我终于明白时,10 已过去,呜呜...

我们对着上面实例图形进行说明,如下:

每个点代表一个检测到的基因;

横轴和纵轴用于固定点在空间的位置;

一般横轴是Log2(fold change),点越偏离中心,表示差异倍数越大;

纵轴是-Log 10 (P-value),点越靠图的顶部表示差异越显著;

点的大小和颜色也可以表示更多的属性,如下图中点的颜色标记其对应的基因是上调, 下调还是无差异;

大小也可用于展示基因表达的平均丰度,一般我们关注表达水平较高且差异较大的基因用于后续的分析和验证。

火山图理解常见的几个问题,公众号生信宝典,如下:

什么是 fold change?翻译成中文是差异倍数,简单来说就是基因在一组样品中的表达值的均值除以其在另一组样品中的表达值的均值。所以火山图只适合展示两组样品之间的比较。

为什么要做 Log 2转换?两个数相除获得的结果 (fold change)要么大于1,要么小于1,要么等于1。这是一句正确的废话吧?那么对应于基因差异呢?简单说,大于1表示上调(可以描述为上调多少倍),小于1表示下调(可以描述为下调为原来的多少分之多少)。大于1可以到多大呢?多大都有可能。小于1可以到多小呢?最小到0。用原始的fold change描述上调方便,描述下调不方便。绘制到图中时,上调占的空间多,下调占的空间少,展示起来不方便。所以一般会做Log 2转换。默认我们都会用两倍差异 (fold change == 2 | 0.5)作为一个筛选标准。Log2转换的优势就体现出来了,上调的基因转换后Log2 (fold change)都大于等于1,下调的基因转换后Log2 (fold change)都小于等于-1。无论是展示还是描述是不是都更方便了。

P-value 都比较熟悉,统计检验获得的是否统计差异显著的一个衡量值,约定成俗的P-value



【本文地址】


今日新闻


推荐新闻


CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3