R语言绘制差异火山图示例

ggplot2绘制差异火山图示例

前几篇分别简介了limma、edgeR、DESeq2、EBSeq的差异基因分析方法，对于两组数据间的差异表达基因的可视化，我们一般就以火山图的样式呈现。

本篇简单展示R包ggplot2绘制差异火山图的几个示例。

本篇示例数据的百度盘链接：

https://pan.baidu.com/s/1dkSiKpEYcuwPqRKAurA9_Q

基因表达差异火山图

提到差异火山图，相信很多同学肯定不陌生。因为形似火山（喷发），所以称为火山图。差异火山图最常见于转录组数据的分析中，在基因表达层面，用于展示两组间表达量上调和下调的基因。常规的火山图中主要包含了两个重要信息，差异表达倍数（Fold Change值，简称FC，作图时会对FC进行log转化，根据logFC值的正负判断这些基因的表达量是上调了还是下调了）以及统计学显著性p值（p-value，通常是FDR校正后的p值，根据校正后p值判断基因表达量上调或下调是否具有显著性）。因此在判断差异基因时，与常规的统计学方法相比，除了p值，通常还会考虑差异倍数，即结合这两个统计结果筛选表达量显著上调或下调的基因（一般而言，差异倍数不能太小）。

如下图示例，癌组织与正常组织的基因表达的差异火山图。红色点代表了癌组织中表达量显著上调的基因，蓝色的点代表了表达量显著下调的基因，灰色的点或者p值未达到显著性水平，或者差异倍数太低不具生物学代表性。

作为一种对差异分析结果的可视化呈现方式，差异火山图实质上就是一种散点图。我们只要准备已经计算好的带有Fold Change值以及显著性p值等信息的做图文件，作图就可以了。如上提到，limma、edgeR、DESeq2、EBSeq等方法可计算这些数值。

网盘附件“GSE40290_Volcano_data.txt”是来自某基因表达芯片数据的差异分析结果，包含了两组间基因表达的log2 Fold Change值以及FDR校正后的p值。

我们使用它绘制火山图，一个常见方法如下。

library(ggplot2)
 
##一个展示基因差异表达水平的火山图示例
gene <- read.delim('GSE40290_Volcano_data.txt', sep = '\t', stringsAsFactors = FALSE)
 
#例如这里自定义根据 |log2FC| >= 1 和 adj.P.Val < 0.01 标记差异类型
gene[which(gene$adj.P.Val < 0.01 & gene$logFC <= -1),'sig'] <- 'Down'
gene[which(gene$adj.P.Val < 0.01 & gene$logFC >= 1),'sig'] <- 'Up'
gene[which(gene$adj.P.Val >= 0.01 | abs(gene$logFC) < 1),'sig'] <- 'None'
 
#横轴 log2FC，纵轴 -log10(adj.P.Val)，颜色表示差异
p <- ggplot(gene, aes(x = logFC, y = -log10(adj.P.Val), color = sig)) +
geom_point(alpha = 0.6, size = 1) +
scale_colour_manual(values  = c('red2', 'blue2', 'gray'), limits = c('Up', 'Down', 'None')) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), plot.title = element_text(hjust = 0.5)) +
theme(legend.key = element_rect(fill = 'transparent'), legend.background = element_rect(fill = 'transparent'), legend.position = c(0.9, 0.93)) +
geom_vline(xintercept = c(-1, 1), color = 'gray', size = 0.3) +
geom_hline(yintercept = -log(0.01, 10), color = 'gray', size = 0.3) +
xlim(-5, 5) + ylim(0, 6) +
labs(x = '\nLog2 Fold Change', y = 'Log10 Mean of Normalized Count\n', color = '', title = 'Cancer vs Normal\n')
 
p
 
#ggsave('gene.pdf', p, width = 5, height = 6)
ggsave('gene.png', p, width = 5, height = 6)
##如果想在图中展示基因的名称，例如这里按 p 值把最显著的上调/下调的前 10 个基因标注出
library(ggrepel)
 
up <- subset(gene, sig == 'Up')
up <- up[order(up$adj.P.Val), ][1:10, ]
down <- subset(gene, sig == 'Down')
down <- down[order(down$adj.P.Val), ][1:10, ]
 
#一种样式，借助 ggrepel 包中的函数
p1 <- p + theme(legend.position = 'right') +
geom_text_repel(data = rbind(up, down), aes(x = logFC, y = -log10(adj.P.Val), label = gene_name),
        size = 3,box.padding = unit(0.5, 'lines'), segment.color = 'black', show.legend = FALSE)
p1
 
#ggsave('gene1.pdf', p1, width = 6, height = 6)
ggsave('gene1.png', p1, width = 6, height = 6)
 
#另一种样式，同样借助 ggrepel 包中的函数，不过这种点太多时容易遮挡
p2 <- p + theme(legend.position = 'right') +
geom_label_repel(data = rbind(up, down), aes(label = gene_name), show.legend = FALSE)
p2
 
#ggsave('gene2.pdf', p2, width = 6, height = 6)
ggsave('gene2.png', p2, width = 6, height = 6)

大致就这样吧，这种散点图画起来还是蛮简单的。

OTUs丰度差异火山图

火山图在基因表达数据之外的情形中也经常能够见到。

比方说在微生物组数据中，也常见到火山图，因为微生物组数据的分析方法和转录组数据的方法在某些方面很相似。这里作者使用三张火山图分别表示了植物根际（Rhizosphere）菌群、跟表（Rhizoplane）菌群、根内（Endosphere）菌群相较于土壤（Soil）菌群，在相对丰度上发生显著富集（enriched）和下降（depleted）的OTUs。不过这篇文献中的火山图样式和上面的有点不同，除了log转化后的Fold Change值，另一个坐标轴使用了OTUs的丰度均值（并作了log转化）而不是显著性p值，p值信息隐含在图中以不同颜色表示的差异OTUs中。当然，对于上述那种基因差异表达类型，两轴也可以像这样展示为log2FC和标准化后的count值（对应于这里就是OTUs丰度了）信息。

（Structure, variation, and assembly of the root-associated microbiomes of rice）

网盘附件“OTUs_diff.txt”是来自某细菌16S高通量测序数据的OTUs丰度差异分析结果，包含了两组间OTUs丰度的log2 Fold Change值以及FDR校正后的p值。

我们使用它绘制相似样式的火山图，如下。

##一个展示细菌群落 OTUs 丰度差异水平的火山图示例
otu <- read.delim('OTUs_diff.txt', sep = '\t', stringsAsFactors = FALSE)

#这里把上述数据框“otu”重复上 3 次，假设当它为 3 个不同分组的比较结果吧，仅用于演示分面图
otu1 <- otu
otu1$group <- 'group1 vs 2'
otu2 <- otu
otu2$group <- 'group2 vs 3'
otu3 <- otu
otu3$group <- 'group1 vs 3'
otu <- rbind(otu1, otu2, otu3)

#例如这里自定义根据 |log2FC| >= 1 和 FDR < 0.05 标记差异类型
otu[which(otu$FDR < 0.05 & otu$log2FC <= -1),'sig'] <- 'depleted'
otu[which(otu$FDR < 0.05 & otu$log2FC >= 1),'sig'] <- 'enriched'
otu[which(otu$FDR >= 0.05 | abs(otu$log2FC) < 1),'sig'] <- 'no diff'

#横轴 log10 转化后的标准化后的 OTUs 平均丰度，纵轴 log2FoldChange，颜色表示差异
p <- ggplot(otu, aes(x = log10(baseMean), y = log2FC, color = sig)) +
geom_point(alpha = 0.6, size = 1) +
scale_colour_manual(values = c('red2', 'blue2', 'gray'), limits = c('enriched', 'depleted', 'no diff')) +
theme(panel.grid.major = element_line(color = 'gray', size = 0.2), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent'), legend.position = 'top') +
labs(x = '\nLog10 Mean of Normalized Abundance', y = 'Log2 Fold Change\n', color = '') +
facet_wrap(~group, ncol = 3, scale = 'free')

p

#ggsave('otu.pdf', p, width = 10, height = 3.5)
ggsave('otu.png', p, width = 10, height = 3.5)

这里直接使用了先前差异分析过程中所得的标准化后的OTUs丰度值，如果觉得不合适，也可以替换为原始的均一化后的OTUs丰度值或者更换为相对丰度。

算是一种拓展样式

网盘附件“gene_two.csv”，来自这么一个试验：

对小鼠分别喂食了A药物和B药物，其中喂食A药物的为A组，喂食B药物的为B组，对照组为C组（control）。并在喂药后探究了相关基因的表达量的变化（A相对于C、B相对于C）。

其中，logCPM，CPM标准化后的基因表达量的平均值，已经过log2转化；logFC_A，A组相对于C组基因表达量的log2 Fold Change值；FDR_A，A组相对于C组基因表达程度的显著性p值，已经过FDR校正；logFC_B，B组相对于C组基因表达量的log2 Fold Change值；FDR_B，B组相对于C组基因表达程度的显著性p值，已经过FDR校正。

然后绘制统计图，尽可能在一张图中展示两组的比较信息，如下所示的样式。横纵坐标均为log2 Fold Change值，分别展示A组相对于C组、B组相对于C组的基因表达量的变化。这里同时关注了A药物和B药物对基因表达量的影响，若某基因的表达量与A、B两种药物均存在关联，则选定该基因（可大致分为“A上调、B上调”，“A下调、B上调”，“A上调、B下调”，“A下调、B下调”4种类型）；若仅与A或B之一的药物有关或均不相关，则忽略该基因；并据此在图中标记为不同颜色。同时，根据CPM标准化后的基因表达量的平均值的log2对数取值，将各基因以不同大小的点表示在图中。

下面是做图代码。

##带双 logFC 信息的二维散点图示例
#读取另一个作图数据
two <- read.csv('gene_two.csv')

#标记显著性（默认 p < 0.05）
two[which(two$FDR_A < 0.05 & two$FDR_B < 0.05),'type1'] <- 'sign'
two[which(two$FDR_A >= 0.05 | two$FDR_B >= 0.05),'type1'] <- 'no'

#标记差异倍数（默认 |log2FC| >= 1）
two[which(two$logFC_A <= -1 & two$logFC_B <= -1),'type2'] <- 'a_down.b_down'
two[which(two$logFC_A >= 1 & two$logFC_B <= -1),'type2'] <- 'a_up.b_down'
two[which(two$logFC_A <= -1 & two$logFC_B >= 1),'type2'] <- 'a_down.b_up'
two[which(two$logFC_A >= 1 & two$logFC_B >= 1),'type2'] <- 'a_up_b_up'
two[is.na(two$type2),'type2'] <- 'no'

#合并显著性和差异倍数，用于标记差异基因
two$type3 <- paste(two$type1, two$type2, sep = '.')

#排序，为了使作图时显著的点绘制在前方（减少被遮盖）
two$type3 <- factor(two$type3, levels = c('sign.a_down.b_down', 'sign.a_up.b_down', 'sign.a_down.b_up', 'sign.a_up_b_up', 'no.a_down.b_down', 'no.a_up.b_down', 'no.a_down.b_up', 'no.a_up_b_up', 'sign.no', 'no.no'))
two <- two[order(two$type3, decreasing = TRUE), ]

#ggplot2 作图，点颜色定义为基因差异类型，点大小表示基因表达量 CPM 值
p <- ggplot(two, aes(logFC_A, logFC_B)) +
geom_point(aes(color = type3, size = logCPM), alpha = 0.6, show.legend = FALSE) +
scale_size(range = c(0, 4)) +
scale_color_manual(limits = c('sign.a_down.b_down', 'sign.a_up.b_down', 'sign.a_down.b_up', 'sign.a_up_b_up', 'no.a_down.b_down', 'no.a_up.b_down', 'no.a_down.b_up', 'no.a_up_b_up', 'sign.no', 'no.no'), values = c('red', 'orange', 'purple', 'blue', 'gray', 'gray', 'gray', 'gray', 'gray', 'gray')) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent')) +
geom_vline(xintercept = c(-1, 1), lty = 2) +
geom_hline(yintercept = c(-1, 1), lty = 2) +
labs(x = 'log2FC (group A vs C)', y = 'log2FC (group B vs C)')

p

#在合适的位置添加文字标记（当然，选择 AI、PS 后续添加也很方便）
p <- p +
annotate('text', label = 'A Down\nB Down', -7, -9) +
annotate('text', label = 'A Down\nB Up', -7, 12) +
annotate('text', label = 'A Up\nB Down', 9, -9) +
annotate('text', label = 'A Up\nB Up', 9, 12)

p

#ggsave('gene_two.pdf', p, width = 7, height = 7)
ggsave('gene_two.png', p, width = 7, height = 7)

友情链接