干货丨如何绘制新颖好看的火山图?

火山图（Volcano plot）是一种比较“远古”的一种散点图，广泛应用于转录组、蛋白组等组间差异分析结果的展示。图表的横轴一般展示差异倍数的变化，而纵轴表示差异分析结果的可靠性。通过找到显著差异的基因，研究者可以进一步做后续的实验。

在生物信息分析中，火山图是非常常见的一种数据展示形式。由于火山图可以非常清晰的展示出哪些基因在不同样本中是具有差异表达显著性的基因，因此在生物医学中，常应用于病例和对照组的转录组研究，也能应用于基因组，蛋白质组，代谢组等统计数据。

常见火山图类型：

那如何做出并读懂这类图，下面进入实战部分吧！

01数据的准备

绘制火山图的方法非常多，常用的有：使用R语言中的ggplot2包、EnhancedVolcano包，Origin等绘制高颜值火山图，下面我们介绍‘使用R语言中的ggplot2包’方法绘制火山图。

首先，安装并且载入相关的R包

install.packages(c("ggrepel", "dplyr", "gt"))

# 载入R包
library(dplyr) # 用于数据处理
library(ggrepel) # 给图片贴上标签
library(gt) # 制作表格

接下来准备作图的数据。

今天用到的RNA-seq数据是经{DESeq2}包处理后的，数据的格式是文本text，可直接导入到R中：

mydata <- read.table("https://gist.githubusercontent.com/stephenturner/806e31fce55a8b7175af/raw/9479acb809fae09aa50ea7df54a1199b3f1ffa11/results.txt", header = TRUE, sep = "")

等数据下载完毕，再查看数据的前六行， 同样是使用代码:

# 查看数据
head(mydata) %>% # 查看数据的前六行
gt() # 制作成表格

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

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

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

◆ 纵轴是-Log 10 (adjusted P-value)，点越靠图的顶部表示差异越显著。

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

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

02制作基础版本的火山图

现在万事俱备，开始作图：

volc_plot <- ggplot(mydata, aes(log2FoldChange, -log10(pvalue))) +# 将p值进行-log10转化
geom_point(size = 0.4) + # 调整点的大小
xlab(expression("log"[2]*" fold change")) + # 修饰横坐标的题目
ylab(expression("-log"[10]*" p-value")) + # 修饰纵坐标的题目
scale_x_continuous(limits = c(-2.5, 2.5)) + # 调整横坐标的范围
theme_minimal()

volc_plot

基本的火山图已经搞定！接下来根据基因的表达水平进行着色。

03根据基因的表达水平进行着色

首先，我们设定一下基因表达水平的规则：

将“fold change”≥2[即log2(2) = 1]并且p＜0.05的基因标记为：“上调”； 将“fold change” ≤0.5[即log2(0.5) = -1]并且p＜0.05的基因标记为：“下调”； 剩下其他的基因则标记为：“不变”。

数据部分可以这么操作：

# 数据处理
mydata1 <- mydata %>%
mutate(expression = case_when(log2FoldChange >= 1 & pvalue < 0.05 ~ "Up-regulated", # 上调
log2FoldChange <= -1 & pvalue < 0.05 ~ "Down-regulated", # 下调
TRUE ~ "Unchanged")) # 不变

# 查看数据
head(mydata1) %>% # 查看前六行
gt() # 制作表格

从上可知，新加了一个名为expression的变量，用于标注基因的表达情况。

接下来，根据expression的不同，将基因进行着色：

volc_plot1 <- ggplot(mydata1, aes(log2FoldChange, -log10(pvalue))) +
geom_point(size = 1.4, aes(color = expression)) + # 根据expression水平进行着色
xlab(expression("log"[2]*" fold change")) + # 修饰x轴题目
ylab(expression("-log"[10]*" p-value")) + # 修饰y轴题目
scale_x_continuous(limits = c(-2.5, 2.5)) +
scale_color_manual(values = c("steelblue", "grey", "red")) + # 添加三种颜色
theme_minimal()

volc_plot1

左侧蓝色的点为下调的基因（并且有统计学意义），右侧红色的点为上调的基因（并且有统计学意义），而中间灰色的指的是表达水平差异不显著的基因。

现在终于明白，为啥将p值进行-log10转化，将“fold change”进行log2转化了。因为转化后，只要看图的左上角（蓝色）和右上角（红色），就马上可以识别出下调和上调的基因！

04分别标注出上调或下调最显著的5个基因

因为将所有基因的名字都标注起来是不可能的，整个画面会像一团浆糊！所以小编只标注上调或下调最显著的前5个基因，共10个，可以这么做：

top_10 <-bind_rows(
mydata1 %>%
filter(expression == 'Up-regulated') %>%
arrange(pvalue, desc(abs(log2FoldChange))) %>%
head(5),
mydata1 %>%
filter(expression == 'Down-regulated') %>%
arrange(pvalue, desc(abs(log2FoldChange))) %>%
head(5)
)

top_10 %>% gt()

接下来，使用top_10里面的基因名称标注火山图：

volc_plot2 <- volc_plot1 +
geom_label_repel(data = top_10,
aes(log2FoldChange, -log10(pvalue), label = Gene),
size = 2)
volc_plot2

火山图的本质就是一张散点图，我们只要理清了绘图逻辑，无论用什么软件，都可以绘制出漂亮的火山图！