Question

How to design a Volcano plot

0

Entering edit mode

3.7 years ago

suxxa • 0

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 • 4.7k views

ADD COMMENT • link updated 3.7 years ago by Kevin Blighe 88k • written 3.7 years ago by suxxa • 0

1

Entering edit mode

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).

ADD REPLY • link 3.7 years ago by ATpoint 85k

0

Entering edit mode

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?

ADD REPLY • link 3.7 years ago by suxxa • 0

Ram · Accepted Answer · 2021-03-24

2

Entering edit mode

3.7 years ago

Kevin Blighe 88k

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.

ADD COMMENT • link 3.7 years ago by Kevin Blighe 88k

0

Entering edit mode

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?

ADD REPLY • link 3.7 years ago by suxxa • 0

1

Entering edit mode

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

ADD REPLY • link 3.7 years ago by Kevin Blighe 88k

0

Entering edit mode

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" ...

ADD REPLY • link updated 3.7 years ago by Ram 44k • written 3.7 years ago by suxxa • 0

0

Entering edit mode

Please deactivate selectLab and then try again

ADD REPLY • link 3.7 years ago by Kevin Blighe 88k

0

Entering edit mode

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

ADD REPLY • link 3.7 years ago by suxxa • 0

0

Entering edit mode

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

ADD REPLY • link 3.7 years ago by Kevin Blighe 88k

0

Entering edit mode

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

ADD REPLY • link updated 3.7 years ago by rpolicastro 13k • written 3.7 years ago by suxxa • 0

1

Entering edit mode

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

ADD REPLY • link 3.7 years ago by Kevin Blighe 88k

0

Entering edit mode

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

ADD REPLY • link 3.7 years ago by suxxa • 0

1

Entering edit mode

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

ADD REPLY • link 3.7 years ago by Kevin Blighe 88k

0

Entering edit mode

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

ADD REPLY • link 3.7 years ago by suxxa • 0

1

Entering edit mode

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

ADD REPLY • link 3.7 years ago by Kevin Blighe 88k

0

Entering edit mode

Many thanks Kevin now

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

ADD REPLY • link 3.7 years ago by suxxa • 0

0

Entering edit mode

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

ADD REPLY • link 3.7 years ago by Kevin Blighe 88k

0

Entering edit mode

Sure, many thanks again!

ADD REPLY • link 3.7 years ago by suxxa • 0

1

Entering edit mode

No problem. No more space

ADD REPLY • link 3.7 years ago by Kevin Blighe 88k