---
title: "Publication-ready volcano plots with enhanced colouring and labeling"
author: "Kevin Blighe"
date: "`r Sys.Date()`"
package: "`r packageVersion('EnhancedVolcano')`"
output:
  html_document:
    toc: true
    toc_depth: 3
    number_sections: true
    theme: united
    highlight: tango
fig_width: 7
bibliography: library.bib
vignette: >
    %\VignetteIndexEntry{Publication-ready volcano plots with enhanced colouring and labeling}
    %\VignetteEngine{knitr::rmarkdown}
    %\VignetteEncoding{UTF-8}
    %\usepackage[utf8]{inputenc}
---

# Introduction

Volcano plots represent a useful way to visualise the results of differential
expression analyses. Here, we present a highly-configurable function that
produces publication-ready volcano plots [@EnhancedVolcano]. EnhancedVolcano
will attempt to fit as many transcript names in the plot window as possible,
thus avoiding 'clogging' up the plot with labels that could not otherwise
have been read.


```{r, echo=FALSE}

    library(knitr)

    opts_chunk$set(tidy = TRUE, message = FALSE, warning = FALSE)

```

# Installation

## 1. Download the package from Bioconductor

```{r getPackage, eval=FALSE}

    if (!requireNamespace("BiocManager", quietly = TRUE))

        install.packages("BiocManager")

        BiocManager::install("EnhancedVolcano")

```

Note: to install development version:

```{r getPackageDevel, eval=FALSE}

    devtools::install_github("kevinblighe/EnhancedVolcano")

```

## 2. Load the package into R session

```{r Load, message=FALSE}

    library(EnhancedVolcano)

```


# Quick start

For this example, we will follow the tutorial (from Section 3.1) of [RNA-seq workflow: gene-level
exploratory analysis and differential expression](http://master.bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html). Specifically, we will load the 'airway' data, where different airway smooth muscle cells were treated with dexamethasone.

```{r}

    library(airway)

    library(magrittr)

    data("airway")

    airway$dex %<>% relevel("untrt")

```

Conduct differential expression using DESeq2 in order to create 2 sets of results:

```{r}

    library("DESeq2")

    dds <- DESeqDataSet(airway, design = ~ cell + dex)
    dds <- DESeq(dds, betaPrior=FALSE)
    res1 <- results(dds,
        contrast = c("dex","trt","untrt"))
    res1 <- lfcShrink(dds,
        contrast = c("dex","trt","untrt"), res=res1)
    res2 <- results(dds,
        contrast = c("cell", "N061011", "N61311"))
    res2 <- lfcShrink(dds,
        contrast = c("cell", "N061011", "N61311"), res=res2)

```

## Plot the most basic volcano plot

For the most basic volcano plot, only a single data-frame or -matrix of
test results is required, containing transcript names, log2FC, and adjusted
or unadjusted P values. The default cut-off for log2FC is >|2|; the default
cut-off for P value is 0.05.

```{r ex1, fig.height = 8, fig.width = 6, fig.cap = "Plot the most basic volcano plot."}

    EnhancedVolcano(res1,

        lab = rownames(res1),

        x = "log2FoldChange",

        y = "pvalue")

```

# Advanced features

Virtually all aspects of an EnhancedVolcano plot can be configured for the
purposes of accommodating all types of statistical distributions and labelling
preferences. EnhancedVolcano will only attempt to label genes that pass the
thresholds that you set for statistical significance, i.e., 'pCutoff' and
'FCcutoff'. In addition, it will only label as many of these that can
reasonably fit in the plot space. The user can optionally supply a vector of
transcript names (as 'selectLab') that s/he wishes to label in the plot.

## Modify cut-offs for log2FC and P value; add title; adjust point and label size

The default P value cut-off of 0.05 may be too relaxed for most studies, which
may therefore necessitate increasing this threshold by a few orders of magnitude.
Equally, the log2FC cut-offs may be too stringent, given that moderated 'shrunk'
estimates of log2FC differences in differential expression analysis can now be
calculated.

In this example, we also modify the point and label size, which can help to
improve clarity where many transcripts went into the differential expression analysis.

```{r ex2, fig.height = 8, fig.width = 6, fig.cap = "Modify cut-offs for log2FC and P value; add title; adjust point and label size."}

    EnhancedVolcano(res2,

        lab = rownames(res2),

        x = "log2FoldChange",

        y = "pvalue",

        pCutoff = 10e-12,

        FCcutoff = 1.5,

        transcriptPointSize = 1.5,

        transcriptLabSize = 3.0,

        title = "N061011 versus N61311")

```

## Adjust colour and alpha for point shading

The default colour scheme may not be to everyone's taste. Here we make it such that
only the transcripts passing both the log2FC and P value thresholds are coloured red,
with everything else black. We also adjust the value for 'alpha', which controls the
transparency of the plotted points: 1 = 100% opaque; 0 = 100% transparent.

```{r ex3, fig.height = 8, fig.width = 6, fig.cap = "Adjust colour and alpha for point shading."}

    EnhancedVolcano(res2,

        lab = rownames(res2),

        x = "log2FoldChange",

        y = "pvalue",

        pCutoff = 10e-12,

        FCcutoff = 1.5,

        transcriptPointSize = 1.5,

        transcriptLabSize = 3.0,

        title = "N061011 versus N61311",

        col=c("black", "black", "black", "red3"),

        colAlpha = 1)

```

## Adjust axis limits

The x-axis limits for log2FC defaults to the max and min of the log2FC values passed
to EnhancedVolcano. This can often render the plot asymmetrical; so, the user may
wish to set these axis limits to the same absolute values, e.g., c(-8, 8). One can
also modify the y-axis limits, but this should be a less common occurrence.


```{r ex4, fig.height = 8, fig.width = 9, fig.cap = "Adjust axis limits."}

    EnhancedVolcano(res2,

        lab = rownames(res2),

        x = "log2FoldChange",

        y = "pvalue",

        pCutoff = 10e-12,

        FCcutoff = 1.5,

        transcriptPointSize = 1.5,

        transcriptLabSize = 3.0,

        title = "N061011 versus N61311",

        colAlpha = 1,

        xlim = c(-8, 8),

        ylim = c(0, -log10(10e-32)))

```

## Adjust cut-off lines

The lines that are drawn to indicate cut-off points are also modifiable. The
parameter 'cutoffLineType' accepts the following values: "blank", "solid",
"dashed", "dotted", "dotdash", "longdash", and "twodash". The colour and
thickness of these can also be modified with 'cutoffLineCol' and
'cutoffLineWidth'. To disable the lines, set either cutoffLineType="blank"
or cutoffLineWidth=0.

```{r ex5, fig.height = 8, fig.width = 9, fig.cap = "Adjust cut-off lines."}

    EnhancedVolcano(res2,

        lab = rownames(res2),

        x = "log2FoldChange",

        y = "pvalue",

        pCutoff = 10e-12,

        FCcutoff = 1.5,

        transcriptPointSize = 1.5,

        transcriptLabSize = 3.0,

        title = "N061011 versus N61311",

        colAlpha = 1,

        xlim = c(-8, 8),

        ylim = c(0, -log10(10e-32)),

        cutoffLineType = "twodash",

        cutoffLineCol = "red3",

        cutoffLineWidth = 1.5)

```

## Adjust legend position, size, and text

The position of the legend can also be changed to "left" or "right"
(and stacked vertically), or "top" or "bottom" (stacked horizontally). The
legend text, label size, and icon size can also be modified.

```{r ex6, fig.height = 8, fig.width = 12, fig.cap = "Adjust legend position, size, and text."}

    EnhancedVolcano(res2,

        lab = rownames(res2),

        x = "log2FoldChange",

        y = "pvalue",

        pCutoff = 10e-12,

        FCcutoff = 1.5,

        transcriptPointSize = 1.5,

        transcriptLabSize = 3.0,

        colAlpha = 1,

        cutoffLineType = "twodash",

        cutoffLineCol = "red4",

        cutoffLineWidth = 1.0,

        legend=c("NS","Log (base 2) fold-change","P value",
            "P value & Log (base 2) fold-change"),

        legendPosition = "right",

        legendLabSize = 14,

        legendIconSize = 5.0)

```

Note: to make the legend completely invisible, specify:

```{r eval=FALSE}

legend=c("","","",""), legendLabSize=-1, legendIconSize=-1

```

## Plot adjusted p-values

Volcano plots do not have to be plotted with nominal (unadjusted P values).
Simply provide a column name relating to adjusted P values and you can also
generate a volcano with these. In this case, the cutoff for the P value 
then relates to the adjusted P value. Here, we also modify the axis titles
by supplying an expression via the bquote function.

```{r ex7, fig.height = 8, fig.width = 6, fig.cap = "Plot adjusted p-values."}

    EnhancedVolcano(res2,

        lab = rownames(res2),

        x = "log2FoldChange",

        y = "padj",

        xlab = bquote(~Log[2]~ "fold change"),

        ylab = bquote(~-Log[10]~adjusted~italic(P)),

        pCutoff = 0.0001,

        FCcutoff = 1.0,

        xlim=c(-6,6),

        transcriptLabSize = 3.0,

        colAlpha = 1,

        legend=c("NS","Log2 FC","Adjusted p-value",
            "Adjusted p-value & Log2 FC"),

        legendPosition = "bottom",

        legendLabSize = 10,

        legendIconSize = 3.0)

```

## Fit more labels by adding connectors

In order to maximise free space in the plot window, one can fit more
transcript labels by adding connectors from labels to points, where
appropriate. The width and colour of these connectors can also be
modified with widthConnectors and colConnectors, respectively.

The result may not always be desirable as it can make the plot
look overcrowded.

```{r ex8, fig.height = 8, fig.width = 6.5, fig.cap = "Fit more labels by adding connectors."}

    EnhancedVolcano(res2,

        lab = rownames(res2),

        x = "log2FoldChange",

        y = "padj",

        xlab = bquote(~Log[2]~ "fold change"),

        ylab = bquote(~-Log[10]~adjusted~italic(P)),

        pCutoff = 0.0001,

        FCcutoff = 2.0,

        xlim = c(-6,6),

        transcriptLabSize = 3.0,

        colAlpha = 1,

        legend=c("NS","Log2 FC","Adjusted p-value",
            "Adjusted p-value & Log2 FC"),

        legendPosition = "bottom",

        legendLabSize = 10,

        legendIconSize = 3.0,

        DrawConnectors = TRUE,

        widthConnectors = 0.2,

        colConnectors = "grey30")

```

## Only label key transcripts

In many situations, people may only wish to label their key transcripts
/ transcripts of interest. One can therefore supply a vector of these
transcripts via the 'selectLab' parameter, the contents of which have
to also be present in the vector passed to 'lab'. In addition, only
those transcripts that pass both the cutoff for log2FC and P value will
be labelled.

```{r ex9, fig.height = 8, fig.width = 6, fig.cap = "Only label key transcripts."}

    EnhancedVolcano(res2,

        lab = rownames(res2),

        x = "log2FoldChange",

        y = "padj",

        selectLab = c("ENSG00000106565","ENSG00000187758"),

        xlab = bquote(~Log[2]~ "fold change"),

        ylab = bquote(~-Log[10]~adjusted~italic(P)),

        pCutoff = 0.0001,

        FCcutoff = 2.0,

        xlim = c(-6,6),

        transcriptPointSize = 1.8,

        transcriptLabSize = 5.0,

        colAlpha = 1,

        legend=c("NS","Log2 FC","Adjusted p-value",
            "Adjusted p-value & Log2 FC"),

        legendPosition = "bottom",

        legendLabSize = 10,

        legendIconSize = 3.0)

```

## Modify border and remove gridlines

One can also modify the border by adding a complete border around the
entire plot, and also its width and colour. The gridlines can also be
removed.

```{r ex10, fig.height = 8, fig.width = 6, fig.cap = "Modify border and remove gridlines"}

    EnhancedVolcano(res2,

        lab = rownames(res2),

        x = "log2FoldChange",

        y = "padj",

        selectLab = c("ENSG00000106565","ENSG00000187758"),

        xlab = bquote(~Log[2]~ "fold change"),

        ylab = bquote(~-Log[10]~adjusted~italic(P)),

        pCutoff = 0.00000000000000001,

        FCcutoff = 2.0,

        xlim = c(-6,6),

        transcriptLabSize = 3.0,

        colAlpha = 1,

        legend=c("NS","Log2 FC","Adjusted p-value",
            "Adjusted p-value & Log2 FC"),

        legendPosition = "bottom",

        legendLabSize = 10,

        legendIconSize = 3.0,

        DrawConnectors = FALSE,

        border = "full",

        borderWidth = 1.5,

        borderColour = "black",

        gridlines.major = FALSE,

        gridlines.minor = FALSE)

```


## Plot multiple volcanos on the same page

One can also plot multiple volcanos on the same plot via the use
of the grid and gridExtra packages.

```{r ex11, fig.height = 8, fig.width = 11, fig.cap = "Plot multiple volcanos on the same page."}

    p1 <- EnhancedVolcano(res1,

        lab = rownames(res1),

        x = "log2FoldChange",

        y = "pvalue",

        pCutoff = 10e-24,

        FCcutoff = 2.0,

        transcriptLabSize = 2.5,

        colAlpha = 1,

        legendPosition = "bottom",

        legendLabSize = 10,

        legendIconSize = 3.0)


    p2 <- EnhancedVolcano(res2,

        lab = rownames(res2),

        x = "log2FoldChange",

        y = "padj",

        selectLab = c("ENSG00000106565","ENSG00000187758"),

        xlab = bquote(~Log[2]~ "fold change"),

        ylab = bquote(~-Log[10]~adjusted~italic(P)),

        pCutoff = 0.0001,

        FCcutoff = 2.0,

        xlim = c(-6,6),

        transcriptLabSize = 5.0,

        colAlpha = 1,

        legend=c("NS","Log2 FC","Adjusted p-value",
            "Adjusted p-value & Log2 FC"),

        legendPosition = "bottom",

        legendLabSize = 10,

        legendIconSize = 3.0)


    library(gridExtra)

    library(grid)

    grid.arrange(p1, p2, ncol=2, top="EnhancedVolcano")

    grid.rect(gp=gpar(fill=NA))

```

## Override colouring scheme with custom key-value pairs

In certain situations, one may wish to over-ride the default colour scheme
with their own colour-scheme, such as colouring transcripts by pathway or
group. This can be achieved by supplying a named vector as colOverride.

In this example, we just wish to colour all transcripts with log2FC > 2 as
'high' and those with log2FC < -2 as 'low'.

```{r ex12, fig.height = 8, fig.width = 11, fig.cap = "Override colour"}

    # create custom key-value pairs for high, low, mid expression
    keyvals <- rep("black", nrow(res2))
    names(keyvals) <- rep("Mid / NA", nrow(res2))

    keyvals[which(res2$log2FoldChange > 2)] <- "gold"
    names(keyvals)[which(res2$log2FoldChange > 2)] <- "High"

    keyvals[which(res2$log2FoldChange < -2)] <- "royalblue"
    names(keyvals)[which(res2$log2FoldChange < -2)] <- "Low"

    p1 <- EnhancedVolcano(res2,

        title = "Over-ride",

        lab = rownames(res2),

        x = "log2FoldChange",

        y = "padj",

        selectLab = rownames(res2)[which(names(keyvals) %in% c("High", "Low"))],

        xlab = bquote(~Log[2]~ "fold change"),

        ylab = bquote(~-Log[10]~adjusted~italic(P)),

        pCutoff = 0.00000000000000001,

        FCcutoff = 1.0,

        xlim = c(-6,6),

        transcriptLabSize = 2.0,

        colAlpha = 1,

        legendPosition = "top",

        transcriptPointSize = 1.2,

        legendLabSize = 15,

        legendIconSize = 5.0,

        DrawConnectors = TRUE,

        widthConnectors = 0.3,

        colConnectors = "grey50",

        border = "partial",

        borderWidth = 1.5,

        borderColour = "black",

        gridlines.major = TRUE,

        gridlines.minor = FALSE,

        colOverride = keyvals)


    p2 <- EnhancedVolcano(res2,

        title = "No over-ride",

        lab = rownames(res2),

        x = "log2FoldChange",

        y = "padj",

        selectLab = rownames(res2)[which(names(keyvals) %in% c("High", "Low"))],

        xlab = bquote(~Log[2]~ "fold change"),

        ylab = bquote(~-Log[10]~adjusted~italic(P)),

        pCutoff = 0.00000000000000001,

        FCcutoff = 1.0,

        xlim = c(-6,6),

        transcriptLabSize = 2.0,

        colAlpha = 1,

        legendPosition = "top",

        transcriptPointSize = 1.2,

        legendLabSize = 15,

        legendIconSize = 5.0,

        DrawConnectors = FALSE,

        widthConnectors = 0.3,

        colConnectors = "grey50",

        border = "full",

        borderWidth = 1.0,

        borderColour = "black",

        gridlines.major = TRUE,

        gridlines.minor = FALSE,

        colOverride = NULL)


    library(gridExtra)

    library(grid)

    grid.arrange(p1, p2, ncol=2, top="EnhancedVolcano")

    grid.rect(gp=gpar(fill=NA))

```

# Acknowledgments

The development of *EnhancedVolcano* has benefited from contributions
and suggestions from:

Sharmila Rana,

[Myles Lewis](https://www.qmul.ac.uk/whri/people/academic-staff/items/lewismyles.html),

Luke Dow - Assistant Professor at Weill Cornell Medicine

# Session info

```{r}

sessionInfo()

```

## References

@EnhancedVolcano
