R绘制火山图 EnhancedVolcano+ggplot

什么是火山图？

火山图其实是一种很形象的叫法，它可以通过关注对象的落点从而直观地展示该对象的所属区域。其通常用于展示差异结果，比如RNA-seq差异基因展示。读懂了“火山”火星喷射的落点横纵坐标意义，就读懂了火山图：

图片1.png

X轴：Log2 Fold Change, Fold Change指样本间表达量的差异倍数；取Log2是为了让差异较大的和差异比较小的数值在视觉上缩小距离 (e.g. 原来2倍的差异等同于1Log2FC)。一般默认取log2FC绝对值大于1为差异基因的筛选标准。

Y轴：-Log(adjust P-value), 对矫正后的P值取负对数(-log)；矫正P值为多重假设检验矫正过的差异显著性P值。由于转录组测序的差异表达分析是对大量的基因表达值进行独立的统计假设检验会存在假阳性问题，因此在进行差异表达分析过程中，采用了公认的Benjamini-Hochberg校正方法对原有假设检验得到的显著性p值（p-value）进行校正，剔除假阳性。

红点: p<0.05且log2FC>1的基因; 蓝点：p<0.05且log2FC< -1的基因; 灰色点：FC<2,即|log2FC|<1。

R绘制火山图

下面就进入主题，用R绘制火山图，我们的教程小白也能看懂哦😊！
R包：Enhanced Volcano Plot的实例操作解说 EnhancedVolcano by Kevin Blighe

1. 安装

  if (!requireNamespace('BiocManager', quietly = TRUE))
    install.packages('BiocManager')
  BiocManager::install('EnhancedVolcano')
  BiocManager::install('DESeq2')
  BiocManager::install("airway")
  BiocManager::install("magrittr")

2. 加载数据集，先安装R包airway，magrittr。

为了快速入门及有效展示，我们参考EnhancedVolcano的官方解说，之后会对部分参数解说，调整参数画出长在自己审美点的火山图。（因为R包EnhancedVolcano是更新的，如果参数有变，请参照报错提示修改代码）

  library(EnhancedVolcano)
  library(airway)
  library(magrittr)
  data('airway')
  airway$dex %<>% relevel('untrt')

  ens <- rownames(airway)
  library(org.Hs.eg.db)
  symbols <- mapIds(org.Hs.eg.db, keys = ens, column = c('SYMBOL'), 
  keytype    = 'ENSEMBL')
  symbols <- symbols[!is.na(symbols)]
  symbols <- symbols[match(rownames(airway), names(symbols))]
  rownames(airway) <- symbols
  keep <- !is.na(rownames(airway))
  airway <- airway[keep,]

3. DESeq2筛选差异基因

  library('DESeq2')
  dds <- DESeqDataSet(airway, design = ~ cell + dex)
  dds <- DESeq(dds, betaPrior=FALSE)
  res <- results(dds, contrast = c('dex','trt','untrt'))

4. EnhancedVolcano绘制

  volcanoplot<-EnhancedVolcano(res,
    lab = rownames(res),
    x = 'log2FoldChange',
    y = 'pvalue')
ggsave("volcanoplot.png",volcanolot)

volcanoplot.png

5. 调整绘图参数

指南中的火山图过于凌乱，调整参数使火山图更加简洁

##收敛结果，会改变log2FC，不改变P指。目的是为了使火山图
res <- lfcShrink(dds, contrast = c('dex','trt','untrt'), res=res, type = 'normal')
#去除P值为NA的基因
res<-na.omit(res)
#筛选差异显著基因，由于该实验中显著且|log2FC|>1的点太多，我们制定更严格的筛选标准，令|log2FC|>2
celltype<-rownames(res[res$padj<0.05&abs(res$log2FoldChange)>2,])
head(celltype)
p1<-EnhancedVolcano(res,
    lab = rownames(res),
    selectLab=celltype,
    axisLabSize=13,
    title="Treat vs Untreat",
    subtitle=NULL,
    titleLabSize=15,
    #caption = paste0('Total = ', nrow(toptable), ' variables'),   
    captionLabSize=11,
    #drawConnectors = T,
    #widthConnectors = 0.2,
    #colConnectors = 'grey50',
    #boxedLabels=T,
    labSize=2.5,  
    x = 'log2FoldChange',
    y = 'padj',
    xlim = c(-5,5),
    pointSize=2,
   #shape = c(6, 6, 19, 16),
    colAlpha=0.7,
    gridlines.major=F,
    gridlines.minor=F,
    border="full",
    borderWidth=1,  
   #boxedLabels=T,
    cutoffLineType="longdash",
    cutoffLineCol="red",
    legendVisible=F, 
    pCutoff=0.05,
    FCcutoff=2,
    #vline=c(-2,2),
    xlab = bquote(~log[2]~ 'fold change'),
    ylab=bquote(~-log[10]~'adjusted p-value')
)

p2<-p1+theme(axis.text.x = element_text(color="black", size=12),
axis.text.y = element_text(color="black", size=12),
plot.title = element_text(hjust = 0.5))
ggsave(p2,file="volcanoplot2.png")

volcanoplot2.png

这样的火山图是不是错落有致，简洁明了，是你想象中的模样😍。

6. ggplot2绘制火山图（使上下调显著差异基因显示不同的颜色）

res$threshold<-as.factor(ifelse(res$log2FoldChange >= 2,'Up',ifelse(res$padj<0.05 & res$log2FoldChange <= -2,'Down','Not')))
res<-data.frame(res)
p3<-ggplot(data=res, aes(x=log2FoldChange, y=-log10(padj), colour=threshold, fill=threshold)) + 
 scale_color_manual(values=c("blue", "grey","red"))+
 geom_point(alpha=0.6,size=2) +
 #xlim(c(-6, 6)) +
 #ylim(c(0, 300)) +
 theme_bw(base_size = 12, base_family = "Times") +
 geom_vline(xintercept=c(-2,2),lty=1,col=c("blue","red"),lwd=0.6)+
 geom_hline(yintercept = -log10(0.05),lty=2,col="black",lwd=0.6)+
 theme(legend.position="right",
 panel.grid=element_blank(),
 legend.title = element_blank(),
 legend.text= element_text(face="bold", color="black",family = "Times", size=8),
 plot.title = element_text(hjust = 0.5),
 axis.text.x = element_text(color="black", size=12),
 axis.text.y = element_text(color="black", size=12),
 axis.title.x = element_text(face="bold", color="black", size=12),
 axis.title.y = element_text(face="bold",color="black", size=12))+
 labs(x="log2 (Fold Change)",y="-log10 (p-value)",title="Treat vs Untreat")

volcanoplot3.png

这版的风格跟EnhancedVolcano的保持一致了，我们就是这么专一😜！
如果你有什么问题，在留言区域留言哦，如果你喜欢我们的文章，请点个赞吧，还可以收藏我们的简书账号，从9月开始，我们会定时更新哦!