Displaying points which fall outside of the plot window in EnhancedVolcano
3
0
Entering edit mode
8 weeks ago
os306 ▴ 10

Hello,

I am using DESeq2 for a differential gene expression analysis of bulk RNAseq data. I am attempting to generate a volcano plot using EnchancedVolcano and was wondering whether it is possible to plot points which fall outside of the figure window. For example, the function plotMA displays points which fall out of the window as open triangles pointing either up or down. Is it possible to do this using EnhancedVolcano (or some other method?)

Thanks very much

EnhancedVolcano plot volcano • 866 views
0
Entering edit mode
0
Entering edit mode

Thanks very much for taking the time to reply. Sorry if this is a silly question, but I had a look at that part of the vignette and I still wasn’t sure how to achieve what I’m looking for (i.e. open triangles to represent outliers that fall outside of the plot window). Would be very grateful for any guidance!

1
Entering edit mode

There is no equivalent function in EnhancedVolcano for this. I would suggest that you simply modify the limits of the x axis via xlim, or, if you want to genuinely reproduce this functionality in EnhancedVolcano, you could achieve it by modifying the underlying log [base 2] fold change values in your input object for these points, and also change their shape to be open triangles.

It's obviously quite difficult to account for every eventuality, and the EnhancedVolcano function / package is already 'parameter overloaded'. I am tentative about doing any further modifications.

0
Entering edit mode

Hi, thanks for your reply. I totally get it. Regardless, it's a fantastic package and very intuitive, particularly for a novice like myself. Thank you for making it!

1
Entering edit mode

No problem. You could try the solution by ATpoint (below) and adapt it to make a volcano, but, then, one can begin to see how there is a learning curve with ggplot2's 'unique' format.

4
Entering edit mode
7 weeks ago

ATpoint made me try this! Unfortunately, it is not as easy with EnhancedVolcano, and I may suggest a customised way via ggplot2, as per ATpoint's code.

Credit to Carlo Yague for the left and right arrow codes.

From the vignette (https://github.com/kevinblighe/EnhancedVolcano):

library(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') p1 <- EnhancedVolcano(res, lab = rownames(res), x = 'log2FoldChange', y = 'pvalue', drawConnectors = TRUE, legendPosition = 'none')  Now, for the new volcano, where absolute log2 FC 2.5 is our cut-off for display purposes: res.new <- res cutoff <- 2.5 # identify genes that pass the threshold of 2.5 up <- rownames(subset(res, log2FoldChange >= cutoff)) up down <- rownames(subset(res, log2FoldChange <= cutoff * -1)) down # impute log2 FCs in the object to be max 2.5 or min -2.5 res.new$log2FoldChange <- ifelse(res.new$log2FoldChange >= cutoff, cutoff, ifelse(res.new$log2FoldChange <= cutoff * -1, cutoff * -1,
res.new$log2FoldChange)) max(res.new$log2FoldChange, na.rm = TRUE)
min(res.new$log2FoldChange, na.rm = TRUE) # custom shapes for the points customshape <- rep(19, nrow(res.new)) names(customshape) <- rep('group1', nrow(res.new)) customshape[which(rownames(res.new) %in% up)] <- -9658 names(customshape)[which(rownames(res.new) %in% up)] <- 'group2' customshape[which(rownames(res.new) %in% down)] <- -9668 names(customshape)[which(rownames(res.new) %in% down)] <- 'group3' # custom sizes for the points customsize <- rep(2.0, nrow(res.new)) customsize [which(rownames(res.new) %in% up)] <- 8 customsize [which(rownames(res.new) %in% down)] <- 8 p2 <- EnhancedVolcano(res.new, lab = rownames(res.new), x = 'log2FoldChange', y = 'pvalue', drawConnectors = TRUE, shapeCustom = customshape, legendPosition = 'none', pointSize = customsize, xlim = c(-2.5,2.5))  Now draw the 2 plots: cowplot::plot_grid(p1, p2, ncol = 2)  ADD COMMENT 1 Entering edit mode Thank you, once again, for taking the time to do this. This, and the package from ATpoint below, are exactly what I was looking for. ADD REPLY 2 Entering edit mode 7 weeks ago ATpoint 53k I also wrote a little package a while ago that can do this type of plot, both for Volcanos and MA-plots. It allows the user to set the x- and ylimits and will plot points beyond these limits as trianges, and also allows to color points when you provide cutoffs for x/y-axis and p-values. You can install as described here: https://github.com/ATpoint/vizzy#installation Then type this below to get the help section for the plots. library(vizzy) ?ggMAplot  For Volcanos, say you want to plot a certain data range at the x-axis (here -2 to 3), with color highlighting significant genes: library(vizzy) library(DESeq2) #/ example data: dds <- DESeq2::DESeq(DESeq2::makeExampleDESeqDataSet(5000,10)) res <- DESeq2::results(dds) %>% data.frame %>% na.omit res$baseMean<-log2(res$baseMean+1) theme_set(theme_classic(base_size = 12.5)) #/ plot: ggMAplot(xval = res$log2FoldChange, yval = -log10(res$pvalue), pval = res$pvalue,
title = "DE results", subtitle = "Volcano plot", preset = "volcano", xlim=c(-2, 3),
x.ablines = c(-log2(1.5), log2(1.5)))


0
Entering edit mode

This is exactly what I was looking for. I'm very grateful to you and Kevin for taking the time to help me.

0
Entering edit mode

Hi there, just to follow up on this, I used your package to generate my volcano plots and got exactly what I was looking for. I now wish to add labels to the genes on the plot (similar to what EnhancedVolcano does). Can I do this using ggrepel()? I wasn't sure how to implement that since vizzy is a ggplot wrapper. Once again, thanks for your patience and guidance!

1
Entering edit mode

I just added a basic option for this. If you install the most recent version and then look at the last two code examples in the Example section and the last paragraph in the Details this should basically be what you need. Something like this:

As you will see in the examples you will need to customly add these points via e.g. ggrepel. I just added an option to directly build the gene names into the ggplot object so that you can easily use + to add further options as you like. I am not an advocate for adding too many custom options to a plotting wrapper but rather add a general option which users then can customize with whatever ggplot extension they like. See examples section.

0
Entering edit mode

Thank you, this is brilliant.

3
Entering edit mode
7 weeks ago
ATpoint 53k

Here is a basic function that you can adapt to make a custom MA-plot because you mentioned plotMA. It might look a bit intimidating at first but learning ggplot is a super handy and important skill in bioinformatics/R. I tried to comment each line of code, please comment if you need help. You could wrap this code into a function that allows to easily customize it.

library(ggplot2)
library(dplyr)
library(DESeq2)

#/ simulate data:
set.seed(2021)
dds <- DESeq2::DESeq(DESeq2::makeExampleDESeqDataSet(5000,10))
res <- DESeq2::results(dds) %>% data.frame %>% na.omit
res$baseMean<-log2(res$baseMean+1)
res$padj[abs(res$log2FoldChange) > log2(1.5)] <- 0.001

#/ see the range of fold changes:
summary(res\$log2FoldChange)

#/ and say you want to plot from -1 to 1:
my_ylimits <- c(-1, 1)

#/ set a nice theme:
theme_set(theme_minimal(base_size = 15))

#/ transform data and prepare the plot:
res %>%

# define what is significant and what not (p<0.05)
levels=c("signif", "nonsignif")),

# if beyond the ylimits mark as "outlier"
outlier=factor(case_when(log2FoldChange > max(my_ylimits) ~ "yes",
log2FoldChange < min(my_ylimits) ~ "yes",
TRUE ~ "no"),
levels=c("yes", "no")),

# set fold changes beyond ylimits to exactly the ylimit
log2FoldChange=case_when(log2FoldChange > max(my_ylimits) ~ max(my_ylimits),
log2FoldChange < min(my_ylimits) ~ min(my_ylimits),
TRUE ~ log2FoldChange)) %>%

# plot colored by signif, and shaped by outlier
ggplot(aes(x=baseMean, y=log2FoldChange, color=signif, shape=outlier)) +
geom_point() +

# make outliers triangles and all others dots:
scale_shape_manual(values=c(17, 20)) +

#/ make signif red and others grey
scale_colour_manual(values=c("firebrick", "grey50")) +

#/ and remove the part of the legend saying outlier or not as we do not need it:
guides(shape=FALSE)


0
Entering edit mode

Thank you for this, much appreciated. I will attempt to adapt this to a volcano plot created in ggplot2 and will comment here if I get stuck.