RNA-seq学习：No.12 EnhancedVolcano绘

EnhancedVolcano包可根据差异分析结果，基于ggplot2绘图结构，方便地绘制美观的火山图，下面根据自己的理解小结下基本用法。

• 官方全面的教程：https://github.com/kevinblighe/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,]

library('DESeq2')
dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds <- DESeq(dds, betaPrior=FALSE)
res <- results(dds,
               contrast = c('dex','trt','untrt'))
res <- lfcShrink(dds,
                 contrast = c('dex','trt','untrt'), res=res, type = 'normal')
res <- as.data.frame(res)
head(res)
#             baseMean log2FoldChange      lfcSE       stat       pvalue         padj
# TSPAN6   710.0931707    -0.37807189 0.09851236 -3.8404448 0.0001228116 0.0009522932
# TNMD       0.0000000             NA         NA         NA           NA           NA
# DPM1     521.2572396     0.19826365 0.10931684  1.8155169 0.0694445184 0.1910397405
# SCYL3    237.6068046     0.03234467 0.13821470  0.2371917 0.8125081096 0.9118161375
# C1orf112  58.0358739    -0.08835419 0.25056704 -0.3194810 0.7493618190 0.8773885438
# FGR        0.3194343    -0.08459224 0.15186225 -0.3948862 0.6929268648           NA

如上，只要包含包含基因名、差异倍数、P值三部分信息的差异结果就可以用于绘制火山图。

• 安装、加载包

# BiocManager::install('EnhancedVolcano')
library(EnhancedVolcano)

基本绘制

如下代码，需要分别交代基因名；x轴为差异倍数；y轴为P值

EnhancedVolcano(res,
    lab = rownames(res),
    x = 'log2FoldChange',
    y = 'pvalue')

如下图结果，基本绘制了不错的火山图。

EnhancedVolcano()也提供了很多调整的参数，可供优化选择

1、标题修改

• title =主标题
• subtitle = 副标题，默认为 "EnhancedVolcano"
• caption = 图注，默认为基因总数统计

EnhancedVolcano(res,
                lab = rownames(res),
                x = 'log2FoldChange',
                y = 'pvalue',
                title = 'Disease versus Normal',
                subtitle = NULL,
                caption = NULL)

2、阈值修改

• 差异倍数：pCutoff = 10e-32
• p值：FCcutoff = 2

EnhancedVolcano(res,
                lab = rownames(res),
                x = 'log2FoldChange',
                y = 'pvalue',
                pCutoff = 10e-32,
                FCcutoff = 3)

3、点point的修改

• 点的大小：pointSize = 2
• 点的不透明度(0~1)：colAlpha = 1/2
• 点的颜色col = c("grey30", "forestgreen", "royalblue", "red2")，分别对应NS，仅差异倍数，仅P值，差异倍数与P值

EnhancedVolcano(res,
                lab = rownames(res),
                x = 'log2FoldChange',
                y = 'pvalue',
                pointSize = 1,
                colAlpha = 0.3,
                col = c("black","black","black","red"))

4、legend的修改

• legendLabels=修改legend的标签内容
• legendLabSize = 14 修改legend的标签大小
• legendPosition = "top"修改legend的位置
• legendIconSize = 5 修改legend的图标大小

EnhancedVolcano(res,
                lab = rownames(res),
                x = 'log2FoldChange',
                y = 'pvalue',
                legendLabels=c('Not sig.',
                               'Log (base 2) FC',
                               'p-value',
                               'p-value & Log (base 2) FC'),
                legendPosition = 'right',
                legendLabSize = 5,
                legendIconSize = 5.0)

可使用ggplot2的语法 + theme(legend.position="none") 设置取消legend

5、设置point label

• 如上面的图，EnhancedVolcano()会显示部分具有显著意义的基因名。我们也可以自定义修改显示哪些基因的标签

p1 = EnhancedVolcano(res,
                lab = rownames(res),
                x = 'log2FoldChange',
                y = 'pvalue',
                selectLab = c('VCAM1','KCTD12','ADAM12',
                              'CXCL12','CACNB2','SPARCL1',
                              'DUSP1','SAMHD1','MAOA'),
                drawConnectors = TRUE)
p2 = EnhancedVolcano(res,
                     lab = rownames(res),
                     x = 'log2FoldChange',
                     y = 'pvalue',
                     col = c("black","black","black","red"),
                     selectLab = "") +
  theme(legend.position="none")
library(patchwork)
p1 | p2