Courtesy Friederike jared.andrews07
NB - please post questions on Bioconductor Support: https://support.bioconductor.org/
[ direct link to create a post: https://support.bioconductor.org/p/new/post/?tag_val=EnhancedVolcano ]
Yep, there is a way to use EnhancedVolcano to also produce an MA plot.
MA plots plot fold-change versus mean expression (as opposed to fold-change versus -log10 p-values in a Volcano). MA plots have been used since the days of the cDNA array and, while array data was typically already log2-transformed as part of normalisation, RNA-seq data is count-based; thus, we typically use natural log (normalised count + 1)
for RNA-seq data:
1, load airway data
library(airway)
library(magrittr)
data('airway')
airway$dex %<>% relevel('untrt')
2, convert Ensembl genes to gene symbols
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, differential expression with DESeq2
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')
4, Prepare results object
EnhancedVolcano will internally transform the y-axis values by -log10()
; here, we counteract this by creating a new column in the data for base mean.
res$baseMeanNew <- 1 / (10^log(res$baseMean + 1))
Test it:
log(5 + 1)
-log10(1 / (10 ^ log(5 + 1)))
5, generate an EnhancedMA
# define different cell-types that will be shaded
celltype1 <- c('VCAM1','WNT2')
library(EnhancedVolcano)
EnhancedVolcano(res,
lab = rownames(res),
title = 'MA plot',
subtitle = 'EnhancedMA',
x = 'log2FoldChange',
y = 'baseMeanNew',
xlab = bquote(~Log[2]~ 'fold change'),
ylab = bquote(~Log[e]~ 'base mean + 1'),
ylim = c(0, 12),
pCutoff = 1 / (10 ^ 7),
FCcutoff = 2.0,
pointSize = 3.5,
labSize = 4.0,
boxedLabels = TRUE,
colAlpha = 1,
legendLabels = c('NS', expression(Log[2]~FC),
'Mean expression', expression(Mean-expression~and~log[2]~FC)),
legendPosition = 'bottom',
legendLabSize = 16,
legendIconSize = 4.0,
drawConnectors = TRUE,
widthConnectors = 0.75,
# encircle
encircle = celltype1,
encircleCol = 'black',
encircleSize = 2.5,
encircleFill = 'pink',
encircleAlpha = 1/2) + coord_flip()
Excellent work! - thanks!