画火山图
发布网友
发布时间：2022-12-04 07:45
共1个回答
热心网友
时间：2023-02-06 00:53
火山图（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_exp.diff 画图，读入R命令：
注：读表用 read.table 函数
读入后，看看这个表张什么样子：
前面我们交代过了 fold change=处理组fpkm/对照组fpkm
所以：
我们可以看到我们这里重新赋值求的 log2_foldchange ，原来的表格里有 log2.fold_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的要求的，我们还需要对它进行精修，至于怎么精修，明天讲~
