下面是去年实习生的分享
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()
也提供了很多调整的参数,可供优化选择
title =
主标题subtitle =
副标题,默认为 "EnhancedVolcano"caption =
图注,默认为基因总数统计EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'pvalue',
title = 'Disease versus Normal',
subtitle = NULL,
caption = NULL)
pCutoff = 10e-32
FCcutoff = 2
EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'pvalue',
pCutoff = 10e-32,
FCcutoff = 3)
pointSize = 2
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"))
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
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