火山图在线绘制-R绘图|ggplot2火山图的绘制

前两期公众号给小伙伴们推荐了两个做基因差异分析网页交互式小工具,很多小伙伴私信我说,用网页小工具做差异分析,结果文件中的火山图默认只展示差异倍数(即logFC)最大的前十个基因。

如果现在要根据自己的需求想展示具有特殊含义的基因(如和癌症发生增长相关,和免疫相关等等)应该怎样做呢?话不多说,直接上工具。

1.打开网址:

点击“火山图快速绘制工具”

2.输入数据:

火山图的输入数据是limma、DESeq2等差异基因分析软件教程得出的结果。 输入文件为所有基因的差异分析结果且必须包含3列内容 :第一列为Gene_symbol,第二列为差异倍数logFC,第三列为FDR(即adj.P.Val),如下图所示。特别注意;用limma包做差异分析时,输出的结果文件包含7列数据, 需要用Excel表格手动调整为3列 。

3.输入参数

4.结果展示:

上一期 R绘图|ggplot2散点图的绘制 简单介绍了散点图在高通量数据展示上的作用,以及如何绘制?散点图在数据展示上存在局限,只能体现基因的差异幅度,并不能体现统计学意义。因此,在高通量文章中,还有一种较为全面的展示数据特点的工具――火山图,可同时直观体现数据差异幅度和统计学意义。

如下图,火山图体现出一组数据(Test/Con)之间的差异幅度和统计学意义分布。X轴代表log2(FC);Y轴代表-log10(q value),灰色代表无差异基因,红色代表上调基因,绿色代表下调基因。X轴的取值可以是FC,也可以是log2处理后的值。X轴的取值可以是q value,也可以是p value(这两个值都属于统计值,只是使用不同方法得出)。在可视化上,形似火山,火山喷发得越高代表统计值越有意义;火山喷发得越外扩代表差异幅度越大。

今天呢,来简单演示下火山图的绘制方法,有数据的小伙伴可以试试!数据可以是高通量检测下所有基因、蛋白、代谢物、甲基化位点信息及其差异倍数和统计学数值。

扩展:

火山图其实是种可视化形式,它的数据来源不一定非得就是FC和P值,只要数据存在两个特征值(连续性数值),其中一个值有正负属性即可。比如:代谢组学数据的FC和VIP值亦可以用火山图展示;GSEA分析结果中通路的P值和NES值也可以用火山图展示。

更多内容可关注公共号“YJY技能修炼”~~~

往期回顾

R绘图|ggplot2散点图的绘制

R绘图|pheatmap热图绘制――基础篇

R绘图|pheatmap热图绘制――中阶篇

R绘图|pheatmap热图绘制――高阶篇

火山图(Volcano Plot)是做RNA-Seq分析的时候特别常用的一张图,因为它可以非常清晰的展示出哪些基因是我们真正想要的显著性的基因。

标准的火山图常用于展示显著差异表达的基因,这里有两个关键词:显著是指P<0.05,差异表达一般我们按照Fold Change(倍数变化)>=2.0作为标准。

当我们拿到基因表达的P值和Fold Change后,为了用火山图展示结果,

一般需要把倍数进行Log2的转化

,比如某基因在实验组表达水平是对照组的4倍,log2(4)=2,同样的如果是1/4,也就是0.25,转换后的结果就是-2。

同样的道理,

对P值进行-log10的转化

,-log10(0.05)约等于1.30103,由于P值越小表示越显著,所以我们进行-log10(P value)转化后,转化值越大表示差异约显著,比如-log10(0.001)=3 > -log10(0.01)=2 > -log10(0.05)=1.30。

上面看不懂没关系,举例说明:

RNA-Seq测序时会有一个对照组的样本,一个处理组的样本,处理组和对照组相比较,会发现某个基因表达量或变高了或变低了。

用处理组的表达量除以对照组就会得到一个倍数,这个倍数就叫做 Fold Change ,

差异表达一般我们按照Fold Change(倍数变化)>=2.0作为标准。

为了用火山图更直观的显示结果我们要将横轴的 Fold Change 进行Log2处理,对纵轴的P值进行 -log10 的处理。

log2(fold change)=0 也就是图中间的位置,就相当于: 处理组fpkm/对照组fpkm=1 也就是说对照组和处理组之间的基因表达量是没有变化的。

所以说出现在中间位置的点是表达量没有发生变化的基因,越往右边走就说明这些基因的表达量在处理组比对照组大,越往左走说明相对于对照组来说,处理组基因表达量变低。

在该图中,不显著的被设置成了粉色,显著的设置成了绿色。

Y轴是 cuffdiff 算的统计检验的 P-value 取了log10之后又取了一个负数,所以说 -log10(P-value) 越高代表着越显著,也就是说越往左上角和越往右上角的点是越显著的点。

我们用来自于 cuffdiff 的输出文件 gene_ 画图,读入R命令:

注:读表用 函数

读入后,看看这个表张什么样子:

前面我们交代过了 fold change=处理组fpkm/对照组fpkm

所以:

我们可以看到我们这里重新赋值求的 log2_foldchange ,原来的表格里有 _change 这一列,为什么还要重新求呢?因为原来表格里求反了求的是:对照组fpkm/处理组fpkm,我们这里改正下。

我们看下计算好的 log2_foldchange 里面有很多的 Inf 和 NaN ,这是怎么回事?

这是因为有些control组基因的FPKM是0,除数是零,是不符合数学规则的,所以会出现这种情况。

我们需要处理下:

control_FPKM == 0 判断一下control_FPKM的值是否等于0

把 control_FPKM == 0 即除数是0的全都筛选出来

我们可以看到经过 log2_foldchange[control_FPKM == 0 ] 操作,把 log2_foldchange 里的所有 NaN / Inf 选出来了,之后我们把筛选出的 NaN / Inf 强行赋值为0

查看新生成的 log2_foldchange 发现还有一些 -Inf

这是因为当 treat_FPKM=0 即分母是0时,求log2就等于负无穷 -Inf

所以,同理,我们把 treat_FPKM 也处理一下,把 log2_foldchange 里所有的 -Inf 赋值为0,这时我们就把画图过程中的异常点过滤掉了。

此时 log2_foldchange 里所有的数都是有效数字了,原来 NaN / Inf 的地方都变成 0 啦~

前面我们也交代过Y轴是 p-value 取 log10 ,这时得到的是一个负值,我们再把它变成整数所以乘以 -1 即:

此时,X,Y轴都准备好了,我们开始作图

这张图虽然是火山图的样子了,但肯定是不符合paper的要求的,我们还需要对它进行精修,至于怎么精修,明天讲~