R语言绘制火山图

火山图在SCI文章中经常出现，其实火山图就是个散点图，通过不同颜色的散点来表示基因的差异。通常横坐标用log2（fold change）表示，差异越大的基因分布在两端，纵坐标用-log10（pvalue）表示，T检验显著性P值的负对数。通常差异倍数越大的基因T检验越显著，所以往往关注左上角和右上角的值。

今天我们通过R来绘制火山图，先导入我们的数据和R包

library(ggplot2)
library(ggrepel)
bc<-read.csv("E:/r/test/huoshantu.csv",sep=',',header=TRUE)

我们来看一下数据（公众号回复：火山图数据，可以获得数据），Gene为基因的名字，log2FoldChange为已经取了对数的FoldChange，pvalue为P值，padj是校正过后的P值。我们先要把P值转换一下，取它的对数，这样画图起来比较对称好看。
bc$log10pvalue=-log10(bc$pvalue)
转换号数据我们就可以做图啦，我们先做个基础的火山图，其实就是散点图。
ggplot(bc, aes(log2FoldChange,log10pvalue))+
geom_point()

非常简单就画好了，一点都不难把，还可以调整X轴的范围和散点的大小或者颜色

ggplot(bc, aes(log2FoldChange,log10pvalue))+
  geom_point(size = 2)+
  scale_x_continuous(limits = c(-2.5, 2.5))

生成一个新指标对基因进行分类，分为3类，log2FoldChange>=1和P小于0.05分为上调，log2FoldChange<=-1和P小于0.05分为下调，剩下的分为未调整。我们可以使用我们的if语句来轻松做到

bc$expression<-ifelse(bc$log2FoldChange >= 1 & bc$pvalue < 0.05,"Up-regulated",
                      ifelse(bc$log2FoldChange <=-1 & bc$pvalue < 0.05,"Down-regulated","Unchanged"))

生成指标expression后就可以进行不同颜色的基因标注

ggplot(bc, aes(log2FoldChange,log10pvalue))+
  geom_point(size = 2,aes(color = expression))+
  scale_x_continuous(limits = c(-2.5, 2.5))

还可以加上标线

ggplot(bc, aes(log2FoldChange,log10pvalue))+
  geom_point(size = 2,aes(color = expression))+
  scale_x_continuous(limits = c(-2.5, 2.5))+
  geom_hline(yintercept=-log10(0.05),linetype=4)+
  geom_vline(xintercept=c(-1,1),linetype=4)

左侧线外的点表示下调超过2倍的基因，右侧线外的点表示上调超过2倍的基因，大于横轴的点表示有显著先的基因。这样我们就把基因区别开来了，看起来一目了然。
还可以加上X轴和Y轴的标签

ggplot(bc, aes(log2FoldChange,log10pvalue))+
  geom_point(size = 2,aes(color = expression))+
  scale_x_continuous(limits = c(-2.5, 2.5))+
  geom_hline(yintercept=-log10(0.05),linetype=4)+
  geom_vline(xintercept=c(-1,1),linetype=4)+
  xlab(expression("log"[2]*" fold change"))+
  ylab(expression("-log"[10]*" p-value"))

接下来我们来对最显著的基因进行标注，上调和下调各标5个，使用dplyr包的管道函数可以轻松帮我们做到

library(dplyr)
data<- bind_rows(
  bc %>%
    filter(expression == 'Up-regulated') %>%
    arrange(pvalue,desc(abs(log2FoldChange))) %>%
    head(5),
  bc %>%
    filter(expression == 'Down-regulated') %>%
    arrange(pvalue,desc(abs(log2FoldChange))) %>%
    head(5)
)

如果不想用管道函数，普通方法也可以实现，就是麻烦点

a<-subset(bc,expression=="Up-regulated")#####a为bc中expression=="Up-regulated"的子集
a1<-head(a[order(abs(a[,2]),decreasing=TRUE),],5)
b<-subset(bc,expression=="Down-regulated")#####b为bc中expression=="Down-regulated"的子集
b1<-head(b[order(abs(b[,2]),decreasing=TRUE),],5)
data1<-rbind(a1,b1)

生成了数据data后就可以进一步绘图了，先把原来绘图函数等于g1

g1<-ggplot(bc, aes(log2FoldChange,log10pvalue))+
  geom_point(size = 2,aes(color = expression))+
  scale_x_continuous(limits = c(-2.5, 2.5))+
  geom_hline(yintercept=-log10(0.05),linetype=4)+
  geom_vline(xintercept=c(-1,1),linetype=4)+
  xlab(expression("log"[2]*" fold change"))+
  ylab(expression("-log"[10]*" p-value")) 

其实就是在原来函数的基础上加上一个geom_label_repel函数进行标签标识

g1+geom_label_repel(data = data,
                       aes(log2FoldChange, -log10(pvalue), label = Gene),
                       size = 2)  

还可以调整标签的字体、位置、颜色大小等，我这里就
参考文献

1. https://www.it610.com/article/1276835499339694080.htm
2. https://mp.weixin.qq.com/s/ApZ6CtnixN9nD4rOTMzwtA
3. http://www.360doc.com/content/17/0725/00/19913717_673898707.shtml

更多精彩文章请关注公众号：零基础说科研