跟着 Hindawi 学画图：漂亮的火山图

原创 JunJunLab 老俊俊的生信笔记 2022-08-15

收录于合集

#100天生信/数据科学自我挑战 100 个

#R绘图 104 个

学习自己喜欢的是一种态度

1、来源

在一篇文献里看到了几个火山图，还挺好看的，但这个期刊没怎么听过：

下面是火山图：

火山图是来自 3 个不同数据集的差异结果，分别对应 GSE775、GSE19322、GSE97494，作者也提供了数据。下面我们来重复这些火山图。

2、分析和观察

我们挑第二个火山图来测试，近距离观察一下细节：

注意：
点的大小不是按值的属性映射的！而是有 3 种大小的点，所以得靠图层叠加分别改变点的大小。
点具有透明度属性。
基因类型分为 5 类：超级显著上调和下调的，显著上调和下调的，不显著变化的基因。
以上类型基因的阈值选定标准。
X、Y 轴标签文本有下标和斜体。
带有特定基因的标签。

3、绘图

分析完以后画图就变得比较清晰明了了，首先我们得对数据的基因进行分类，分为 5 类。

3.1、基因分类

正常的火山图分为 3 类即可，一般 ifelse 函数就可以快速分类，但是对于比较多的分类，这个函数写多了容易看不懂，阅读性不高，我们可以使用 case_when 函数分类，代码清晰明了：

# 加载R包
library(ggplot2)
library(ggrepel)
library(ggnewscale)
library(tidyverse)

# 设置工作目录
setwd('c:/Users/admin/Desktop/')

# 读取数据
vo <- read.table('GSE19322.txt',header = T)

# 查看数据
head(vo,3)
    Gene     logFC  AveExpr         t     P.Value   adj.P.Val         B
1 Cacnb2 -2.712469 5.420691 -8.048005 0.000017900 0.000317613  3.328462
2 Synpo2 -2.427160 7.858508 -4.183735 0.002221717 0.006269240 -1.719538
3 Gm3579 -2.348745 7.581799 -3.879394 0.003541202 0.008918279 -2.203795

由图可以看出 超级（上调或下调）基因 的 -log10 adj.pvalue >= 3 且 log2FC > 2 或 < -2，我们按照 5 类基因的阈值来筛选：

# 计算-log10 adj.P.Val
vo$Log10.adj.Pvalue <- -log10(vo$adj.P.Val)

# 分类
vo$type <- case_when(vo$logFC > 2 & vo$Log10.adj.Pvalue >= 3 ~ 'Super-Up',
                     vo$logFC < -2 & vo$Log10.adj.Pvalue >= 3 ~ 'Super-Down',
                     vo$logFC >= 1 & vo$adj.P.Val < 0.05 ~ 'Up',
                     vo$logFC <= -1 & vo$adj.P.Val < 0.05 ~ 'Down',
                     vo$adj.P.Val >= 0.05 ~ 'NoSig',
                     abs(vo$logFC) < 1 ~ 'NoSig'
)

# 查看分类数量
table(vo$type)
      Down      NoSig Super-Down   Super-Up         Up
       437      12451          9         17        134

3.2、绘制第一层

我们先绘制第一层 不显著基因 的图层：

# 绘制不显著基因图层
nosig <- ggplot(vo,aes(x = logFC,y = Log10.adj.Pvalue)) +
  # 不显著基因图层
  geom_point(data = vo %>% filter(type == 'NoSig' ),
             # 点大小
             size = 1,
             color = 'grey',alpha = .5) +
  # 修改颜色
  new_scale('color') +
  scale_color_manual(values = c('Up' = '#FF7878','Down' = '#B3E283'),
                     name = 'Normal DEGs') +
  # 横坐标刻度及范围调整
  scale_x_continuous(limits = c(-5,5),breaks = seq(-5,5,1)) +
  # 水平虚线
  geom_hline(yintercept = c(-log10(0.05),3),
             size = 0.7,
             color = 'grey',
             lty = 'dashed') +
  # 竖直虚线
  geom_vline(xintercept = c(-2,-1,1,2),
             size = 0.7,
             color = 'grey',
             lty = 'dashed') +
  theme_bw(base_size = 18) +
  theme(panel.grid = element_blank()) +
  # 添加x轴标签
  xlab(bquote(Log[2] ~ '(fold change)')) +
  # 添加y轴标签
  ylab(bquote(-Log[10] ~ italic('adj.Pvalue')))

nosig

3.3、绘制第二层

这一层我们添加 显著性变化 的基因上去：

# 绘制显著基因图层
sig <- nosig +
  # 显著基因图层
  geom_point(data = vo %>% filter(type == 'Up' | type == 'Down'),
             # 点大小
             size = 3,
             aes(color = type),alpha = .5) +
  # 修改颜色
  new_scale('color') +
  scale_color_manual(values = c('Super-Up' = '#D32626','Super-Down' = '#438A5E'),
                     name = 'Super Anno')

sig

3.4、绘制第三层

接下来添加超级变化的基因图层上去，点要更大一些：

# 绘制超级显著基因图层
supersig <- sig +
  # 显著超级显著基因图层
  geom_point(data = vo %>% filter(type == 'Super-Down'| type == 'Super-Up'),
             # 点大小
             size = 5,
             aes(color = type),alpha = .6)

supersig

3.5、添加基因名

最好添加上感兴趣的基因名即可：

# 添加基因名
anno_gene <- supersig +
  geom_text_repel(data = vo[which(vo$Gene %in% c('Cxcl2','Ctla2a','Il6','Hspa1b','Hspa1a')),],
                  aes(label = Gene),
                  # 标签颜色大小
                  color = '#BE6283',size = 5,
                  # 标签微调
                  nudge_y = -0.1,nudge_x = 0.2,
                  # 标签箭头
                  arrow = arrow(angle = 0,ends = 'last',length = unit(0.5,'cm'))) +
  labs(title = 'GSE19322')

anno_gene

完成，这里我们用同样的方法把数据替换一下就可以画出其它两个火山图了。

3.6、基因名

基因名，添加标签的时候运行一下下面，然后替换成相应的变量即可：

# 基因名
gene_set1 <- c('Hspa1a','Hspa1b','Inhbb','Cxcl2','Gp49a',
               'Ptgs2','Cxcl5','Ccl2','S100a9','S100a8')
gene_set2 <- c('Cxcl2','Ctla2a','Il6','Hspa1b','Hspa1a')
gene_set3 <- c('Il1r2','Timp1','Arg1','Spp1','Chil3',
               'Cxcl2','Slfn4','Retnlg')