让你的火山图更丰富

下面是去年实习生的分享

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

官方全面的教程：https://github.com/kevinblighe/EnhancedVolcano

示例差异基因数据

代码语言：javascript复制

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值三部分信息的差异结果就可以用于绘制火山图。

安装、加载包

代码语言：javascript复制

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

基本绘制

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

代码语言：javascript复制

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

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

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

1、标题修改

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

代码语言：javascript复制

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

2、阈值修改

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

代码语言：javascript复制

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值

代码语言：javascript复制

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的图标大小

代码语言：javascript复制

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()会显示部分具有显著意义的基因名。我们也可以自定义修改显示哪些基因的标签

代码语言：javascript复制

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

ggplot2 legend point position title

0 人点赞