生物信息学习的正确姿势
NGS系列文章包括NGS基础、高颜值在线绘图和分析、转录组分析 （Nature重磅综述|关于RNA-seq你想知道的全在这）、ChIP-seq分析 （ChIP-seq基本分析流程）、单细胞测序分析 (重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程)、DNA甲基化分析、重测序分析、GEO数据挖掘（典型医学设计实验GEO数据分析 (step-by-step)）、批次效应处理等内容。
火山图是散点图的一种，它将统计测试中的统计显著性量度（如p value）和变化幅度相结合，从而能够帮助快速直观地识别那些变化幅度较大且具有统计学意义的数据点（基因等）。常应用于转录组研究，也能应用于基因组，蛋白质组，代谢组等统计数据。
所以关注火山图（其它类型图也是），先理解每个点是什么（点代表基因、样品、通路或其它的，这个认识可以来自于常识，更准确的是看作者的描述），然后看横轴代表什么、纵轴代表什么，再看图例中展示的其他信息，如颜色、大小和形状分别代表什么。这些都理顺了，图理解就不难了。
如图一：
每个点代表一个检测到的基因。
横轴和纵轴用于固定点在空间的位置。
一般横轴是Log2(fold change)，点越偏离中心，表示差异倍数越大。
纵轴是-Log 10 (adjusted P-value)，点越靠图的顶部表示差异越显著。
点的大小和颜色也可以表示更多的属性，如下图中点的颜色标记其对应的基因是上调, 下调还是无差异。
大小也可用于展示基因表达的平均丰度，一般我们关注表达水平较高且差异较大的基因用于后续的分析和验证。
图一（图源：易生信PPT）
火山图理解常见的几个问题
但没想到，在我们易生信培训过程中，对火山图的问题还是比较多的，我们一个个的说一下。
什么是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<0.05为统计检验显著的常规标准。
什么是adjusted P-value?
这里面就涉及到一个统计学问题了。做差异基因检测时，要对成千上万的基因分别做差异统计检验。统计学家认为做这么多次的检验，本身就会引入假阳性结果，需要做一个多重假设检验校正。
这个校正怎么做呢？最简单粗暴的方法是每一次统计检验获得的P-value都乘以总的统计检验的次数获得adjusted P-value (这就是Bonferroni correction）。
但这样操作太严苛了，很容易降低统计检出力，找不到有差异的基因。后续又有统计学家提出相对不这么严苛的计算方法，如holm, hochberg, hommel, BH, BY, fdr等。BH是我们比较常用的一个校正方法，获得的值是假阳性率 FDR （false discovery rate）。
FDR筛选时就可以不用遵循0.05这个标准了。我们可以设置FDR<0.05表示我们容许数据中存在至多5%假阳性率；FDR<0.1表示我们对假阳性率的容忍度至多是10%。当然如果说我们设置FDR<0.5，即数据中最多可能有一半是假阳性就说不过去了。
同样为什么做-Log 10转换呢？
因为FDR值是0-1之间，数值越小越是统计显著，也越是我们关注的。-Log 10 (adjusted P-value)转换后正好是反了多来，数值越大越显著，而且以10为底很容易换算回去。
理解完这些之后，再来看火山图。
整体来看，基因有上调就有下调，图整体是以X=0的垂线左右对称的。如果数据中大部分点都是上调或下调，成偏态分布时，需考虑标准化步骤没有处理好，或数据存在批次效应，导致数据存在系统偏差。
图的左上角和右上角是差异基因集中的地方，也是我们关注的重点。
图一中左侧的火山图还展示了基因表达的平均丰度，即基因在所有样品中表达的均值。一般变化倍数大、平均表达也比较高的基因会更可信，更适合后期实验检测，否则就算变化倍数再大，表达低的基因也难以被检测到。
番外：
差异倍数fold change还有另外一种处理方式。假如有两个样品A和B。如果某个基因在A中表达比较高，则计算fold change是用A/B; 。如果某个基因在B中表达比较高，则计算fold change是用B/A，然后乘以-1; gtools::foldchange是这么操作的。
adjusted P-value, q value, fdr一般代表相同的含义，都是多重假设检验校正后的P-value,可能的区别就在于校正算法的不同。
几个代表性火山图
火山图虽然用的多，但其实能提供的信息算不上多，一般是在上面标记一些关注的基因的名字，然后在正文中做下描述。标记基因名字的方式也比较多，图二中左图的颜色标示是一个不错的选择。
图二（图源：易生信PPT）
图二右图来自2017年发表在Cell的一篇文章-Epigenetic Therapy Ties MYC Depletion to Reversing Immune Evasion and Treating Lung Cancer。
一排火山图放在一起是不是很有气势，更主要的是展示了5种疫苗诱导的差异基因数目显著不同，在图上红点多少展示出的视觉冲击还是优于图标中的数字表示的，更容易留下直观的印象。个人觉得是一个很有特色的火山图案例。
图三
图三来自文章Edwards, J., et al. (2015). PNAS Fig. 2A。
这是一篇16S分析文章较系统的作品，两年被引用147次，推荐阅读。上面的火山图展示了水稻根不同生态位相对于土壤中显著差异的OTU，横坐标是相对丰度平均值（Log10 转换），纵坐标是Log10(fold change)，整体类似于图一中的左图，只是转换了X和Y轴变量。
图四
火山图就是散点图，点的颜色可展示代表性属性。
图四来源— https://arxiv.org/pdf/1103.3434.pdf ：
第6号染色体上的探针/基因用红色标记，在基因注释中带有“细胞因子”的探针/基因用蓝色标记。
增强火山图之在基本火山图的基础上，标注有变量名。
上图共有64102个变量，绿色的点的|log2FC|>1，蓝色的点是在P value <0.0001的基础上，校正得到的符合要求的点。红色的点是满足了以上两点要求的变量。
如有雷同数据，可大胆参照模仿，更多增强火山图见：
传送门（代码）➷➷
火山图绘制 （新版更强大）
最简单的绘制方法是使用我们的在线网站——imageGP(http://www.ehbio.com/ImageGP/)。
还有我们的往期推文：
https://blog.sciencenet.cn/blog-118204-1256841.html
上一篇：
WGCNA分析，简单全面的最新教程（可以在线做了）下一篇：
GO、GSEA富集分析一网打进
