画火山图

Posted 2023-04-02

参考技术A

火山图（Volcano Plot）是做RNA-Seq分析的时候特别常用的一张图，因为它可以非常清晰的展示出哪些基因是我们真正想要的显著性的基因。

标准的火山图常用于展示显著差异表达的基因，这里有两个关键词：显著是指P<0.05，差异表达一般我们按照Fold Change(倍数变化)>=2.0作为标准。
当我们拿到基因表达的P值和Fold Change后，为了用火山图展示结果， 一般需要把倍数进行Log2的转化 ，比如某基因在实验组表达水平是对照组的4倍，log2（4）=2，同样的如果是1/4，也就是0.25，转换后的结果就是-2。

同样的道理， 对P值进行-log10的转化 ，-log10（0.05）约等于1.30103，由于P值越小表示越显著，所以我们进行-log10（P value）转化后，转化值越大表示差异约显著，比如-log10（0.001）=3 > -log10(0.01)=2 > -log10(0.05)=1.30。

上面看不懂没关系，举例说明：

RNA-Seq测序时会有一个对照组的样本，一个处理组的样本，处理组和对照组相比较，会发现某个基因表达量或变高了或变低了。 用处理组的表达量除以对照组就会得到一个倍数，这个倍数就叫做 Fold Change ， 差异表达一般我们按照Fold Change(倍数变化)>=2.0作为标准。 为了用火山图更直观的显示结果我们要将横轴的 Fold Change 进行Log2处理，对纵轴的P值进行 -log10 的处理。
log2(fold change)=0 也就是图中间的位置，就相当于： 处理组fpkm/对照组fpkm=1 也就是说对照组和处理组之间的基因表达量是没有变化的。 所以说出现在中间位置的点是表达量没有发生变化的基因，越往右边走就说明这些基因的表达量在处理组比对照组大，越往左走说明相对于对照组来说，处理组基因表达量变低。
在该图中，不显著的被设置成了粉色，显著的设置成了绿色。

Y轴是 cuffdiff 算的统计检验的 P-value 取了log10之后又取了一个负数，所以说 -log10(P-value) 越高代表着越显著，也就是说越往左上角和越往右上角的点是越显著的点。

我们用来自于 cuffdiff 的输出文件 gene_exp.diff 画图，读入R命令：

注：读表用 read.table 函数
读入后，看看这个表张什么样子：

前面我们交代过了 fold change=处理组fpkm/对照组fpkm
所以：

我们可以看到我们这里重新赋值求的 log2_foldchange ，原来的表格里有 log2.fold_change 这一列，为什么还要重新求呢？因为原来表格里求反了求的是：对照组fpkm/处理组fpkm，我们这里改正下。
我们看下计算好的 log2_foldchange 里面有很多的 Inf 和 NaN ，这是怎么回事？

这是因为有些control组基因的FPKM是0，除数是零，是不符合数学规则的，所以会出现这种情况。
我们需要处理下：
control_FPKM == 0 判断一下control_FPKM的值是否等于0
把 control_FPKM == 0 即除数是0的全都筛选出来

我们可以看到经过 log2_foldchange[control_FPKM == 0 ] 操作，把 log2_foldchange 里的所有 NaN / Inf 选出来了，之后我们把筛选出的 NaN / Inf 强行赋值为0

查看新生成的 log2_foldchange 发现还有一些 -Inf

这是因为当 treat_FPKM=0 即分母是0时，求log2就等于负无穷 -Inf
所以，同理，我们把 treat_FPKM 也处理一下，把 log2_foldchange 里所有的 -Inf 赋值为0，这时我们就把画图过程中的异常点过滤掉了。

此时 log2_foldchange 里所有的数都是有效数字了，原来 NaN / Inf 的地方都变成 0 啦~

前面我们也交代过Y轴是 p-value 取 log10 ，这时得到的是一个负值，我们再把它变成整数所以乘以 -1 即：

此时，X,Y轴都准备好了，我们开始作图

这张图虽然是火山图的样子了，但肯定是不符合paper的要求的，我们还需要对它进行精修，至于怎么精修，明天讲~

快来看看如何使用R语言绘制一张漂亮的火山图

“

火山图（Volcano Plot）是差异基因可视化非常好的工具，本文将向大家展示如何用R语言来实现火山图的绘制

”

在对基因芯片数据、测序数据进行差异表达分析之后，我们通常会选择绘制火山图（Volcano Plot）来对各基因的表达情况进行可视化，正如下图所示：

绘制火山图的工具有很多，本质上就是画一张散点图，然后可以根据自己的需要进行标注。上面这张图是使用R语言的ggplot2包绘制的，除此之外我们还可以用GraphPad Prism软件绘制，或者使用Excel也是可以，当然还有一些在线工具可以使用。

本篇推文，我就来教大家如何使用R语言绘制如上图所示的火山图，同时为了发扬偷懒的精神，方便日后使用，我最后还把绘图过程封装成了一个函数，这样再次使用的时候只用输入数据、简单调整参数，就可以很方便地得到一张好看的火山图啦！

好，废话不多说，我们这就开始吧！

首先准备绘图所需数据，我们需要一个数据框，里面包含了每个基因的表达变化情况，如下图所示：

（此为数据集GSE118370经limma包分析的结果：Tumor vs. Normal，命名为DEG）

我们还需要定义一下上下调基因：以logFC>1且adj.P.Val<0.05为显著上调基因，logFC<-1且adj.P.Val<0.05为显著下调基因，其余为表达变化不显著的基因，并在DEG中新增一列，命名为Change（logFC代表对变化倍数进行log2转换，adj.P.Val为校正后的P值，筛选标准更为严格）。

##给DEG添加一列，标注各基因上下调情况
DEG$Change <- as.factor(ifelse(DEG$adj.P.Val < 0.05 & abs(DEG$logFC) > 1,
                               ifelse(DEG$logFC > 1, 'Up', 'Down'),'NOT'))

加上Change这一列后，就可以开始绘图了。

library(ggplot2)
ggplot(data = DEG, aes(x = logFC, y = -log10(adj.P.Val),colour = Change))+
  geom_point(alpha = 0.5)

怎么样，最简单版的火山图就绘制好了，接下来就是对图片元素进行修饰，例如加上基因的名字、添加水平竖直分割线、调整背景等，具体的修饰过程参见下述代码：

##ggrepel包用来给散点添加标签
library(ggrepel)

##为了给散点添加标签，我们还需要给DEG再加上一列，为基因的名称（其实就是每行的名字）
##把DEG按logFC绝对值从大到小排序
DEG <- DEG[order(abs(DEG$logFC),decreasing = T),]

##新增加一列，为基因的名字
DEG$anno_name <- rownames(DEG)

##只取表达变化前10的基因进行标注，所以其余的都设置为NA
DEG$anno_name[11:nrow(DEG)] <- NA

##设置一个标题
title <- paste0('cutoff for logFC is ',1,
                '
cutoff for adj.P.Val is ',0.05,
                '
the number of up gene is ',nrow(DEG[DEG$Change=='Up',]),
                '
the number of down gene is ',nrow(DEG[DEG$Change=='Down',]))

##现在开始绘图
ggplot(data = DEG, aes(x = logFC, y = -log10(adj.P.Val),colour = Change))+
  geom_point(alpha = 0.5)+
  scale_colour_manual(values = c('blue','grey10','red'))+ #自定义定义散点颜色
  geom_text_repel(aes(label=anno_name),show.legend = F,segment.colour = 'black')+ #给散点添加标签
  xlab('log2(FC)')+ylab('-log10(adj.P.Val)')+ggtitle(title)+ #设置坐标轴及主标题
  geom_hline(yintercept = -log10(0.05),lty=4,lwd=0.6,alpha=0.8)+ #设置水平分割线
  geom_vline(xintercept = c(1,-1),lty=4,lwd=0.6,alpha=0.8)+ #设置竖直分割线
  theme_bw()+theme(panel.border = element_blank(), #设置背景
                   panel.grid.major = element_blank(),
                   panel.grid.minor = element_blank(),
                   axis.line = element_line(colour = "black"))

这样，最开始那张好看的火山图就绘制完成了，大家可以根据自己的喜好对图中元素进行修改。最后放一个大招：把整个绘图过程定义成一个函数，这样我们运行函数之后就只用短短一句话就可以生成火山图啦！

plot_vocalno <- function(DEG,logfccutoff=1,pcutoff=0.05,show.number=10,p=c('adj.p','p')){
##加载绘图所需包
library(ggrepel);library(ggplot2)

##根据输入选择P.Value或者adj.P.Val
if(p[1]=='p'){p <- 'P.Value'} else {
if(p[1]=='adj.p'){p <- 'adj.P.Val'} else {
  print('请输入正确的p！');break}}

#Change列要根据设置的cutoff阈值而定
DEG$Change = as.factor(ifelse(DEG[, colnames(DEG) %in% p[1]] < pcutoff & abs(DEG$logFC) > logfccutoff,
                              ifelse(DEG$logFC > logfccutoff, 'Up', 'Down'),'NOT'))

##还要规定一下show.number值
show.number <- floor(show.number)

##设置后续画图的标签列，以下代码是为了避免那些logFC符合但p值不符合的
if(show.number >= 1){
  DEG$anno_name <- rownames(DEG) #把行名设为打标签的一列
  DEG <- DEG[order(abs(DEG$logFC),decreasing = T),] #按logFC绝对值从大到小进行排序
  k <- 0
  for(i in 1:nrow(DEG)){
    if(DEG$Change[i]=='NOT') {DEG$anno_name[i] <- NA} else {k <- k+1}
    if(k==show.number) break
  }
  DEG$anno_name[{i+1} : nrow(DEG)] <- NA #只要前面的的show.number个，剩下的都为NA
} else {
  DEG$anno_name <- NA #当show.number为0时，整列都设为NA
}

##设置标题
title <- paste0('cutoff for logFC is ',logfccutoff,
                '
cutoff for ',p,' is ',pcutoff,
                '
the number of up gene is ',nrow(DEG[DEG$Change=='Up',]),
                '
the number of down gene is ',nrow(DEG[DEG$Change=='Down',]))

##自定义颜色（因为有时候没有Up或者Down，所以要选择一下）
manualcolor <- c('blue','grey10','red')
if(sum(DEG$Change %in% 'Up')==0) manualcolor <- c('blue','grey10')
if(sum(DEG$Change %in% 'Down')==0) manualcolor <- c('grey10','red')
if(sum(DEG$Change %in% 'Up')==0 & sum(DEG$Change %in% 'Down')==0) manualcolor <- c('grey10')

##根据p的选择，要选择一下纵坐标轴的映射列
if(p=='P.Value'){
  plot_tmp <- ggplot(data = DEG,aes(x=logFC,y=-log10(P.Value),colour=Change))
} else {
  plot_tmp <- ggplot(data = DEG,aes(x=logFC,y=-log10(adj.P.Val),colour=Change))
}

##开始绘图
plot_tmp+
  geom_point(alpha=.5,size=1)+ #画散点
  scale_colour_manual(values = manualcolor)+ #自定义定义散点颜色
  geom_text_repel(aes(label=anno_name),show.legend = F,segment.colour = 'black')+ #给散点添加标签
  xlab('log2(FC)')+ylab(paste0('-log10(',p,')'))+ggtitle(title)+ #设置坐标轴及主标题
  geom_hline(yintercept = -log10(pcutoff),lty=4,lwd=0.6,alpha=0.8)+ #设置水平分割线
  geom_vline(xintercept = c(logfccutoff,-logfccutoff),lty=4,lwd=0.6,alpha=0.8)+ #设置竖直分割线
  theme_bw()+theme(panel.border = element_blank(), #设置背景
                   panel.grid.major = element_blank(),
                   panel.grid.minor = element_blank(),
                   axis.line = element_line(colour = "black"))
}

参数解释：

1. DEG：差异表达分析结果，格式如前图；

2. logfccutoff：即为logFC的筛选阈值，默认为1；

3. pcutoff：即为P值的筛选阈值，默认为0.05；

4. show.number：设置散点标签的个数；

5. p：设置选择adj.P.Val还是P.Value，默认选择adj.P.Val。

（如发现代码有错误之处或不足，欢迎指正！）

—END—