R绘图|ggplot2火山图的绘制

上一期R绘图|ggplot2散点图的绘制简单介绍了散点图在高通量数据展示上的作用，以及如何绘制？散点图在数据展示上存在局限，只能体现基因的差异幅度，并不能体现统计学意义。因此，在高通量文章中，还有一种较为全面的展示数据特点的工具——火山图，可同时直观体现数据差异幅度和统计学意义。

如下图，火山图体现出一组数据（Test/Con）之间的差异幅度和统计学意义分布。X轴代表log2(FC)；Y轴代表-log10(q value)，灰色代表无差异基因，红色代表上调基因，绿色代表下调基因。X轴的取值可以是FC，也可以是log2处理后的值。X轴的取值可以是q value，也可以是p value（这两个值都属于统计值，只是使用不同方法得出）。在可视化上，形似火山，火山喷发得越高代表统计值越有意义；火山喷发得越外扩代表差异幅度越大。

image

今天呢，来简单演示下火山图的绘制方法，有数据的小伙伴可以试试！数据可以是高通量检测下所有基因、蛋白、代谢物、甲基化位点信息及其差异倍数和统计学数值。

常规火山图的绘制

首先清除环境，安装并加载所需要的R包

rm(list = ls()) #清除环境内存
#install.packages("ggplot2") #已经安装的，无需重复安装
library(ggplot2)#加载ggplot2包
library(ggrepel)#用来添加名字

载入数据，设置阈值和标题

data = read.table("volcano.txt", header=TRUE)
head(data,10)#查看data数
##           Gene.Symbol       FC     LogFC     P.Value        P
## 1          uc008bia.2 3.825804 1.9357630 0.002615940 2.582372
## 2          uc029xmk.1 2.338380 1.2255091 0.037124123 1.430344
## 3          uc029roz.1 2.338380 1.2255091 0.037124123 1.430344
## 4                 Tnf 1.982197 0.9871004 0.002971479 2.527027
## 5               jokla 1.961320 0.9718248 0.001067592 2.971595
## 6          uc029wtq.1 1.932778 0.9506760 0.194224499 0.711696
## 7                Ccl3 1.821079 0.8647932 0.003162985 2.499903
## 8           Mir466f-2 1.789048 0.8391922 0.043636948 1.360146
## 9  ENSMUSG00000088978 1.782004 0.8335004 0.021773938 1.662063
## 10      NONMMUG022964 1.709356 0.7734525 0.002087741 2.680323
logFC_cutoff <- with(data,mean(abs( LogFC)) + 2*sd(abs( LogFC)))#设置阈值,常规咱们可直接用abs(FC>2),即logFC=1去设置
logFC_cutoff #查看数据
## [1] 0.1812603
data$sig = as.factor(ifelse(data$P.Value < 0.05 & abs(data$LogFC) > logFC_cutoff,
                              ifelse(data$LogFC > logFC_cutoff ,'UP','DOWN'),'NOT'))#添加sig这栏up,down,not的行数据
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(data[data$sig =='UP',]) ,
                    '\nThe number of down gene is ',nrow(data[data$sig =='DOWN',])) #设置标题，体现阈值、上下调基因数量
this_tile
## [1] "Cutoff for logFC is 0.181\nThe number of up gene is 980\nThe number of down gene is 395"

ggplot2绘制火山图

g = ggplot(data=data, aes(x=LogFC, y=-log10(P.Value), color=sig)) +
  geom_point(alpha=0.4, size=2.5) +
  theme_bw(base_size=15)+
  theme(panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.line = element_line(colour = "black"))+
  xlab("LogFC") + ylab("-Lgp") +
  ggtitle( this_tile ) +
  theme(plot.title = element_text(size=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red'))
print(g)

image.png

添加特定基因名称

g1=g+geom_text_repel(data=subset(data,abs(LogFC)>1.5), aes(label=Gene.Symbol),col="green",alpha = 0.5) #logFC大于1.5的基因标记上名字
print(g1)

image.png

添加分区线

g2<- g+ geom_hline(yintercept=-log10(0.05),colour="black", linetype="dashed") +  #加一条y=-1*log10(0.05)的直线，它的颜色为black，线条类型为dashed即虚线
     geom_vline(xintercept=c(-0.18,0.18),colour="black", linetype="dashed")      #加一条x=0.18和-0.18的直线，它的颜色为black，线条类型为dashed即虚线
print(g2)

image.png

美化版火山图的绘制


m=ggplot(data=data, aes(x=LogFC, y =-log10(P.Value))) +
#将数据分成四象限，对每个象限的数据进行颜色和大小的定义
  geom_point(data=subset(data,abs(data$LogFC) <= 0.18),aes(size=abs(LogFC)),color="gray",alpha=0.5) +
  geom_point(data=subset(data,data$P.Value >= 0.05&abs(data$LogFC) > 0.18),aes(size=abs(LogFC)),color="gray",alpha=0.5) +
  geom_point(data=subset(data,data$P.Value<0.05 & data$LogFC > 0.1826),aes(size=abs(LogFC)),color="red",alpha=0.5) +
  geom_point(data=subset(data,data$P.Value<0.05 & data$LogFC < -0.18),aes(size=abs(LogFC)),color="darkgreen",alpha=0.5) +
  ## 添加分割线（横向+竖向）
  geom_hline(yintercept = -log10(0.05),lty=4,lwd=0.6,alpha=0.8)+
  geom_vline(xintercept = c(0.18,-0.18),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"))+ 
  ## 修改坐标轴
  labs(x="Log2 (fold change)",y="-Log10 (p-value)")+
  ## 去掉图注
  theme(legend.position='none')+
  ## 添加文本信息
  geom_text_repel(data=subset(data, abs(LogFC) > 1), aes(label=Gene.Symbol),color="black",alpha = 0.8)
print(m)

image.png

ggsave(m,file="meihua_volcano.pdf",width = 8,height = 6)#保存为pdf格式

扩展：火山图其实是种可视化形式，它的数据来源不一定非得就是FC和P值，只要数据存在两个特征值（连续性数值），其中一个值有正负属性即可。比如：代谢组学数据的FC和VIP值亦可以用火山图展示；GSEA分析结果中通路的P值和NES值也可以用火山图展示。
更多内容可关注公共号“YJY技能修炼”~~~

往期回顾
R绘图|ggplot2散点图的绘制
 R绘图|pheatmap热图绘制——基础篇
 R绘图|pheatmap热图绘制——中阶篇
 R绘图|pheatmap热图绘制——高阶篇