This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to design a Volcano plot

Hi,

I'm trying to design something similar to this with EnhancedVolcano to see only the upregulated and downregulated genes.

Is it possibile to do it?

enhancedvolcano
• 7,634 views • link
updated by Kevin Blighe Kevin Blighe • written by suxxa • suxxa

Sure, simply filter the data you give to the plot to only include genes that are below the pvalue cutoff you want.

Say your data are df so something like df[df$FDR<0.05,] or tidyverse-like df %>% filter(FDR<0.05).

• 0 views • link
written by ATpoint ATpoint

Thanks, yes I understood and I'll follow your advice, but how can I colour the upregulated genes in red and the downregulated genes in blue? I'm using this code as a template but I'm not making any improvements...

load("results_Liver.V8.rda") 

EnhancedVolcano(resLFC,
            lab = rownames(resLFC),
            x = 'log2FoldChange',
            y = 'padj',
            title = 'Liver - SEX',
            selectLab = c('VCAM1'),
            xlab = bquote(~Log[2]~ 'fold change'),
            pCutoff = 10e-14,
            pointSize = 4.0,
            labSize = 6.0,
            labCol = 'black',
            labFace = 'bold',
            boxedLabels = TRUE,
            colAlpha = 4/5,
            legendLabels=c('Not sig.','Log (base 2) FC','p-adj',
                           'p-adj & Log (base 2) FC'),
            legendPosition = 'right',
            legendLabSize = 14,
            legendIconSize = 4.0,
            drawConnectors = TRUE,
            widthConnectors = 1.0,
            colConnectors = 'black')

This is the result liver

But I would like to colour differently the upregulated and downregulated genes and delete the pvalue -300 limit because there are too many genes squeezed and not easily seen, so maybe i could put -600 as a limit?

• 0 views • link
written by suxxa • suxxa

1 answer

Hi, the vignette, which I wrote, is quite detailed, if you can take a few moments to go through it so that you can see how to, for example, modify colours: https://bioconductor.org/packages/release/bioc/vignettes/EnhancedVolcano/inst/doc/EnhancedVolcano.html

The value of ~300 on the y-axis relates to p-values of 0. There is not much that EnhancedVolcano can do with p-values of 0. The relevant lines of code are: https://github.com/kevinblighe/EnhancedVolcano/blob/master/R/EnhancedVolcano.R#L274-L294

It has been suggested elsewhere to use geom_jitter() to make these values more aesthetically pleasing. I personally see absolutely no issue with them presented as they currently are [presented] in your figure, provided that you explain why they appear this way.

• 0 views • link
written by Kevin Blighe Kevin Blighe

Hi Kevin, thanks for your work, yes I read your vignette but other errors appeared like:

Error in grid.Call(C_convert, x, as.integer(whatfrom), as.integer(whatto),  : Viewport has zero dimension(s)

I'm using this code:

Keyvals etc...

EnhancedVolcano(resLFC,
            lab = rownames(resLFC),
            x = 'log2FoldChange',
            y = 'pvalue',
            selectLab = rownames(resLFC)[which(names(keyvals) %in% c('high', 'low'))],
            xlab = bquote(~Log[2]~ 'fold change'),
            title = 'Liver - SEX',
            pCutoff = 10e-14,
            FCcutoff = 1.0,
            pointSize = 3.5,
            labSize = 4.5,
            colCustom = keyvals,
            colAlpha = 1,
            legendPosition = 'left',
            legendLabSize = 15,
            legendIconSize = 5.0,
            drawConnectors = TRUE,
            widthConnectors = 1.0,
            colConnectors = 'black',
            arrowheads = FALSE,
            gridlines.major = TRUE,
            gridlines.minor = FALSE,
            border = 'partial',
            borderWidth = 1.5,
            borderColour = 'black')

About the limit value of 300 maybe I explained it wrong, is it possible to increase this limit or it's fixed?

• 0 views • link
written by suxxa • suxxa

Hi, can you confirm that all devices are off? - please run dev.off() until it returns an error. Then, what is the output of:

dim(resLFC)
str(resLFC)
str(keyvals)

Again, regarding the p-values, what can we realistically do with p-values of 0? -log10(0.00) == Infinity; so, we cannot plot these. The issue about these is sourced at that program that you use for differential expression - it calculates p-values of 0.

In order to plot these, my design choice in EnhancedVolcano is to automatically convert these p-values of 0 to 10^-1 the lowest non-zero p-value. For example, if our lowest non-zero p-value is 0.0001, then any p-value of 0 will be converted to 0.00001

You can technically manually modify these values in your results table prior to running EnhancedVolcano, if you wish

However, there is always a hard limit, which may be machine specific. This can be given by:

.Machine$double.xmin
[1] 2.225074e-308
-log10(.Machine$double.xmin)
[1] 307.6527
• 0 views • link
written by Kevin Blighe Kevin Blighe

Hi, I run dev.off() and after 1 null device result I obtained an error and I runned another time the script but obtained the same error Error in grid.Call(C_convert, x, as.integer(whatfrom), as.integer(whatto), : Viewport has zero dimension(s)

Output of dim(resLFC)

  [1] 54496     5

Output of str(resLFC)

Formal class 'DESeqResults' [package "DESeq2"] with 7 slots
  ..@ priorInfo      :List of 4
  .. ..$ type         : chr "apeglm"
  .. ..$ package      : chr "apeglm"
  .. ..$ version      :Classes 'package_version', 'numeric_version'  hidden list of 1
  .. .. ..$ : int [1:3] 1 12 0
  .. ..$ prior.control:List of 7
  .. .. ..$ no.shrink            : int [1:3] 1 2 4
  .. .. ..$ prior.mean           : num 0
  .. .. ..$ prior.scale          : num 0.106
  .. .. ..$ prior.df             : num 1
  .. .. ..$ prior.no.shrink.mean : num 0
  .. .. ..$ prior.no.shrink.scale: num 15
  .. .. ..$ prior.var            : num 0.0112
  ..@ rownames       : chr [1:54496] "ENST00000336079.7" "ENST00000430575.1" "ENST00000253320.8" "ENST00000477725.1" ...
  ..@ nrows          : int 54496
  ..@ listData       :List of 5
  .. ..$ baseMean      : num [1:54496] 1185 930 613 1009 311 ...
  .. ..$ log2FoldChange: Named num [1:54496] -8.3 -7.61 -8.63 -9.42 -7.94 ...
  .. .. ..- attr(*, "names")= chr [1:54496] "ENST00000336079.7" "ENST00000430575.1" "ENST00000253320.8" "ENST00000477725.1" ...
  .. ..$ lfcSE         : Named num [1:54496] 0.126 0.119 0.136 0.152 0.156 ...
  .. .. ..- attr(*, "names")= chr [1:54496] "ENST00000336079.7" "ENST00000430575.1" "ENST00000253320.8" "ENST00000477725.1" ...
  .. ..$ pvalue        : num [1:54496] 0 0 0 0 0 0 0 0 0 0 ...
  .. ..$ padj          : num [1:54496] 0 0 0 0 0 0 0 0 0 0 ...
  ..@ elementType    : chr "ANY"
  ..@ elementMetadata:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  .. .. ..@ rownames       : NULL
  .. .. ..@ nrows          : int 5
  .. .. ..@ listData       :List of 2
  .. .. .. ..$ type       : chr [1:5] "intermediate" "results" "results" "results" ...
  .. .. .. ..$ description: chr [1:5] "mean of normalized counts for all samples" "log2 fold change (MAP): SEX 2 vs 1" "posterior SD: SEX 2 vs 1" "Wald test p-value: SEX 2 vs 1" ...
  .. .. ..@ elementType    : chr "ANY"
  .. .. ..@ elementMetadata: NULL
  .. .. ..@ metadata       : list()
  ..@ metadata       :List of 6
  .. ..$ filterThreshold: Named num 7.85
  .. .. ..- attr(*, "names")= chr "0%"
  .. ..$ filterTheta    : num 0
  .. ..$ filterNumRej   :'data.frame':  50 obs. of  2 variables:
  .. .. ..$ theta : num [1:50] 0 0.0194 0.0388 0.0582 0.0776 ...
  .. .. ..$ numRej: num [1:50] 610 603 595 584 577 572 560 554 538 532 ...
  .. ..$ lo.fit         :List of 2
  .. .. ..$ x: num [1:50] 0 0.0194 0.0388 0.0582 0.0776 ...
  .. .. ..$ y: num [1:50] 611 602 594 586 578 ...
  .. ..$ alpha          : num 0.1
  .. ..$ lfcThreshold   : num 0

Output of > str(keyvals)

    Named chr [1:54496] "royalblue" "royalblue" "royalblue" "royalblue" "royalblue" "royalblue" "royalblue" "royalblue" "gold" ...
 - attr(*, "names")= chr [1:54496] "low" "low" "low" "low" ...
• 0 views • link
updated by Ram Ram • written by suxxa • suxxa

Please deactivate selectLab and then try again

• 0 views • link
written by Kevin Blighe Kevin Blighe

even with #selectLab = rownames(resLFC)[which(names(keyvals) %in% c('high', 'low'))], it doesn't work...

• 0 views • link
written by suxxa • suxxa

Please show exactly how you created keyvals. Also the output of head(resLFC), please

• 0 views • link
written by Kevin Blighe Kevin Blighe 
keyvals <- ifelse(

  resLFC$log2FoldChange < -1.4, 'royalblue',

  ifelse(resLFC$log2FoldChange > 1.4, 'gold',

  'black'))

keyvals[is.na(keyvals)] <- 'black'

names(keyvals)[keyvals == 'gold'] <- 'high'

names(keyvals)[keyvals == 'black'] <- 'mid'

names(keyvals)[keyvals == 'royalblue'] <- 'low'

and for head head

• 0 views • link
updated by rpolicastro rpolicastro • written by suxxa • suxxa

Strange, but I cannot see what you are doing locally. The error sometimes appears [if I recall correctly] when the plot window is too small for the output that is being produced. This could be an issue if you are using RStudio.

..and what is the output of:

rownames(resLFC)[which(names(keyvals) %in% c('high', 'low'))]

That could be the issue because, with this, you are asking it to label too many transcripts.

Try with drawConnectors = FALSE and I think that it will generate

• 0 views • link
written by Kevin Blighe Kevin Blighe

Hi Kevin, many thanks I resolved disabling connectors, even using RStudio. Just another question, how could I clean this plot deleting all the ENST references?

result

output of rawnames

here

• 0 views • link
written by suxxa • suxxa

Hey again, in which way do you wish to tidy these?

• 0 views • link
written by Kevin Blighe Kevin Blighe

Thanks again Kevin, I would like to see only the various coloured dots, if possible...

• 0 views • link
written by suxxa • suxxa

Hi, I think that, in that case, you need to use lab = NULL or lab = NA

• 0 views • link
written by Kevin Blighe Kevin Blighe

Many thanks Kevin now

What if I would like to colour black under the pvalue threshold?

• 0 views • link
written by suxxa • suxxa

In order to do that, you will have to go via the keyvals method

• 0 views • link
written by Kevin Blighe Kevin Blighe

Sure, many thanks again!

• 0 views • link
written by suxxa • suxxa

Log in to answer this question.

More posts like this