EnhancedVolcano包可根据差异分析结果，基于ggplot2绘图结构，方便地绘制美观的火山图，下面根据自己的理解小结下基本用法。
示例差异基因数据
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
