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

1. 前言

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

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

2. 实例绘图

1. 数据读取

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

DEG<-read.table("TCGA-COAD-edgeR-DEG.xls",sep = "\t",check.names = F,header = T)
head(DEG)

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

2. 自行写代码实现

########自己写代码
# 绘制火山图====================================
# 设置pvalue和logFC的阈值
cut_off_pvalue = 0.0001
cut_off_logFC = 2
# 根据阈值分别为上调基因设置‘up’，下调基因设置‘Down’，无差异设置‘None’，保存到change列
# 这里的change列用来设置火山图点的颜色
DEG$Sig = ifelse(DEG$PValue < cut_off_pvalue & 
                abs(DEG$logFC) >= 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<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为底很容易换算回去。

火山图就这么实现完成，您学会了没，读懂火山图的意思了没？记得关注公众号 桓峰基因，未来将推出视频教程，敬请关注！

桓峰基因 生物信息分析，SCI文章撰写及生物信息基础知识学习：R语言学习，perl基础编程，linux系统命令，Python遇见更好的你 42篇原创内容 

Reference:

Qiu C C, Su Q S, Zhu S Y, et al. Identification of Potential Biomarkers and Biological Pathways in Juvenile Dermatomyositis Based on miRNA-mRNA Network. BioMed Research International, 2019, 2019.

Lin X D, Wu Y P, Chen S H, et al. Identification of a five‐mRNA signature as a novel potential prognostic biomarker in pediatric Wilms tumo. Molecular genetics & genomic medicine, 2019: e1032

本文使用 文章同步助手 同步