转录组数据分析笔记（8）——ggplot2和ggrepel绘制火山图

主要介绍下在用DESeq2得到我们想要的差异表达基因后，如何在R中用ggplot和ggrepel绘制火山图。

1. 前言

前面介绍了怎么用DESeq2做两组样本的差异基因表达分析，以及怎么用dplyr包给DESeq2运行结果增加一列分组信息，我们先看下两个R包运行结束后生成的gene_0_1.csv文件是怎么样的：

A-G列结果是DESeq2跑的，我们用到的只有基因名，log2FoldChange和padj这三列，通过log2FoldChange绝对值大于1，调整后的pvalue也就是padj（即FDR值）小于0.05筛选除差异表达基因，最后加上group列方便查看。

到这里已经可以通过排序找到我们感兴趣的基因了，但是这样的数据不够直观，我们还可以用最著名的绘图R包ggplot2做个火山图。这里需要准备的绘图R包是ggplot2，还有添加标签的R包ggrepel。

1.2 两个注意点

什么是火山图就不多bb了，重要的是知道我们可以从火山图获得两个信息：差异表达倍数（FoldChange值）和统计学显著性的标志p值。为了更方便比较和作图，我们一般用log2FC代替Fold Change值并作为X轴数据，表示两样品（组）间表达量的比值，对其取以2为底的对数即为log2FC，一般默认取log2FC绝对值大于1为差异基因的筛选标准；用FDR（也就是padj值）代替pvalue，并取-log10（FDR）值作为Y轴，FDR是错误发现率，是pvalue值进行校正得到的。

log2FC有正有负很好理解，可能有同学发现，有的基因明明有pvalue值，但是校正之后的FDR值却是NA，如下：

查阅了一些资料，当基因的在所有样本中表达量为0，则两个值都为NA；当read count数较低时，DESeq2进行Independent Filtering过滤了一部分可能造成假阳性的结果，此时padj值为NA。因此，这部分数据在做火山图的时候因为没有对应的Y值也会被过滤掉。

2. 流程代码

这部分需要一点R语言基础，需要知道怎么改自己的参数。假设前面没有用dplyr包做差异筛选（做了也不影响，只是多一列数据而已）只用DESeq2跑了个结果，我们同样可以用ggplot2包做筛选，用ggrepel包美化做标签。继续用前面DESeq2生成的csv文件，流程代码如下：

library("ggplot2")
gene <- read.csv("gene_0_1.csv", stringsAsFactors = FALSE)
gene[which(gene$log2FoldChange <= -1 & gene$padj < 0.05), "GROUP"] <- "DOWN"
gene[which(gene$log2FoldChange >= 1 & gene$padj < 0.05), "GROUP"] <- "UP"
gene[which(abs(gene$log2FoldChange) < 1 | gene$padj >= 0.05), "GROUP"] <- "NOT CHANGE" # |代表或，和linux里的管道是完全不一样的。以上三步新建了一列GROUP，筛选并赋予了三个值。
res <- ggplot(gene,     # ggplot数据来源，这里省略了data = 和mapping = 
       aes(x = log2FoldChange,      # 表示映射关系，就是定义xy
           y = -log10(padj), 
           col = GROUP)) +      # 注意这里定义颜色用col，以GROUP值区分
         geom_point(alpha = 0.5,    # ggplot做散点图，设置点透明度和大小
                    size = 1) + 
         scale_color_manual(values = c("red","blue","grey"),
                            limits = c("UP","DOWN","NOT CHANGE")) +     # 自定义颜色 
         theme(panel.grid = element_blank(),    # 去网格线
               panel.background = element_rect(color = "black",
                                               fill = "transparent"),   # 去背景色，透明
               plot.title = element_text(hjust = 0.5),      # 调整图标标题位置为中间
               legend.key = element_rect(fill = "transparent"),
               legend.background = element_rect(fill = "transparent"),
               legend.position = "right") +     # 设置legend图标
         geom_vline(xintercept = c(-1, 1),
                    color = "gray",
                    size = 0.3) +       # 设置x轴辅助线
         geom_hline(yintercept = -log10(0.05),
                    color = "gray",
                    size = 0.3) +       # 设置y轴辅助线
         labs(x = "log2 Fold Change",
              y = "-log10(FDR) ",
              title = "LD 1 day vs LD 0 day")       # 设置坐标轴标题和火山图标题
res     # 查看结果，plots中可以查看

library("ggrepel")
up <- subset(gene, GROUP == "UP")       # subset从数据框中筛选符合条件的数据
up <- up[order(up$padj), ][1:5, ]       # order升序排序，取前5个
down <- subset(gene, GROUP == "DOWN")
down <- down[order(down$padj), ][1:5, ]
resdata <- res + 
  geom_text_repel(data = rbind(up, down),       # ggrepel特有的函数
                  aes(x = log2FoldChange,
                      y = -log10(padj),
                      label = X ),      # label值指定哪列做标签
                  size = 3,
                  box.padding = unit(0.5, "lines"),
                  segment.color = "#cccccc",
                  show.legend = FALSE)      # 以上都是特有参数
resdata     # 查看结果
ggsave("gene_0_1.png", resdata, width = 10, height = 10)        # 输出结果文件

3. 代码详解

3.1 ggplot2

1
2
3

gene[which(gene$log2FoldChange <= -1 & gene$padj < 0.05), "GROUP"] <- "DOWN"
gene[which(gene$log2FoldChange >= 1 & gene$padj < 0.05), "GROUP"] <- "UP"
gene[which(abs(gene$log2FoldChange) < 1 | gene$padj >= 0.05), "GROUP"] <- "NOT CHANGE"

之前这里稍微解释一下，即使前面没有用dplyr包，用别的方法同样可以筛选差异基因并且新增一列分组数据，万变不离其宗，核心的判断方式是一样的。如果前面学了linux，注意最后 | 这个符号在R里表示或，不要和管道混淆。

ggplot(data = 输入的数据,
      mapping = aes(x = 定义值,
                   y = 定义值,
                   其他参数)) +
      genom_point(参数) +       # 选择作图方法和参数 
    其他设置函数和参数              # 可有可无，美观相关的东西

我总结了一下ggplot最基本的结构，data和mapping是缺省值，可以不写。

输入的数据可以是表格，可以是数据框等等；aes自定义点的映射范围，大小，颜色等等；作图方法有很多，比如点状图是genom_point。自由度很高，能设置的东西也非常之多，只有两点需要注意，同一个函数不同参数用 , 隔开；不同函数用 + 隔开。

中间的设置函数也稍微解释一下：

scale_color_manual() 该函数是R中的一种自定义配色方法，手动把颜色赋值给参数value。我这里将UP的点赋予了红色，DOWN的点赋予蓝色，其他点赋予灰色。

theme() 该函数与主题配置有关，参数非常多，可以选取需要的比如背景色、网格线等等进行设置。这里举个例子，legend是图标，在ggplot中legend有四部分： legend.tittle, legend.text, legend.key和legend.backgroud，而每一个部分都有四种函数可以嵌套（也是是对应4种处理方式）：element_text()绘制标题相关；element_rect()绘制背景相关；element_blank()空主题，对元素不分配绘图空间；element_get()得到当前主题的设置。每个函数还有相应的参数，说起来就没完没了了。。。常用的设置知道就行。

geom_vline() 和 geom_hline() 这两个函数分别设置x轴y轴辅助线，目的是使我的火山图更直观，从图上可以直接看到我的分类依据。

labs() 该函数自定义x轴y轴标签和图标标题。这里提一嘴，火山图标题也是一个注意点，一般是 处理组vs对照组 ，因此前面也说到DESeq2处理数据要注意顺序问题，不然会得出完全相反的结论，在火山图上的表现为与实际火山图呈镜像对称，这也很好理解。

这里看一下res结果，我们可以在plots中看到缩略的预览图：

3.3 ggrepel

在过滤掉4000多个FDR值为NA的点，我们获得了一个不怎么像火山的火山图 ~~（简直丑爆了）~~ ，但是数据还是挺美丽的，上调区域和下调区域都有前后对比差异非常大的基因：log2FC绝对值越大，差异越明显；-log10（FDR）值越大，可信度越高。

但是这个图还有个缺点，我不知道这些代表性的差异点对应什么基因名，我还要回到excel里去筛选排序。因此，我推荐用ggrepel包对火山图进一步美化，加上基因标签，能一眼看到我感兴趣的基因。

这个包的原理和发展咱就不说了，已经是半夜4点了。。。简单介绍下中间用到的函数的结构。

subset() 从数据框中筛选符合条件的数据，我将UP的点和DOWN的点都提取出来。

order() order是升序排序，因为上调和下调的基因都比较多，全部打上基因名标签是不现实也没有意义的。我按照padj列也就是FDR值进行升序排序，取前5个可信度最高的基因打上基因名标签。当基因较少的时候是可以全部打上标签的。

在前面ggplot2作图的基础上，我们加上geom_text_repel()这个特殊的ggreple包函数，这个函数是基于函数geom_label()做的改良，它将标签置于一个方框中，并且每个标签有算法优化不会重叠。该函数的结构与前面的ggplot前半部分类似：

geom_text_repel(data = 输入数据,
                aes(x = 定义值,
                   y = 定义值,
                   label = X ),
                其他参数

这里所有参数设置都是平行的，所以只需要 , 隔开。

我之前说到提取了up和down的数据，这里我把它们rbind一下合在一起，就形成了新的数据框数据，也就是我只对前面排序筛选的上调下调各5个基因打标签。这里注意下aes()这里的 label 是指定标签的，也就是我们这里的基因名，应该用的行名才能和数据一一对应，这里我用X是偶然发现的一个很有趣的事：

前面导入gene数据框的时候，自动把行名加到了第一列成了单独的一列，且该列列名系统定义为X，我们可以进入environment找到gene点开看看这个数据框结构，如下所示：