A basic volcano plot
Passing a single data frame colours every gene by significance. By
default a gene is called significant when its adjusted
p-value (FDR) is below 0.05 and its
|log2FoldChange| exceeds 1.
ggvolc turns the output of a differential-expression
(DE) analysis into a polished, publication-ready volcano plot. It
accepts the output of DESeq2, edgeR,
and limma directly, highlights and labels genes of
interest, combines the plot with a gt table, and can render
an interactive version with ggiraph.
The package ships with two data frames. all_genes is a
complete DESeq2-style result table; attention_genes is a
small subset of genes you might want to highlight.
data(all_genes)
data(attention_genes)
head(all_genes)
#> genes baseMean log2FoldChange lfcSE stat pvalue
#> 1 GCR1 7201.5782 2.244064 0.2004959 11.192564 4.434241e-29
#> 2 OPI10 1009.4171 -2.257454 0.2096469 -10.767889 4.880607e-27
#> 3 AGA2 249.1173 3.829474 0.3623263 10.569132 4.143136e-26
#> 4 FIM1_1376 5237.5035 2.550409 0.2560379 9.961059 2.256459e-23
#> 5 HMG1 10838.1037 2.214300 0.2229065 9.933763 2.968371e-23
#> 6 FIM1_3918 2456.8070 2.288243 0.2356228 9.711467 2.694309e-22
#> padj
#> 1 2.153711e-25
#> 2 1.185255e-23
#> 3 6.707736e-23
#> 4 2.739905e-20
#> 5 2.883475e-20
#> 6 2.181043e-19Passing a single data frame colours every gene by significance. By
default a gene is called significant when its adjusted
p-value (FDR) is below 0.05 and its
|log2FoldChange| exceeds 1.
Supply a second data frame and those genes are drawn with a black outline and labelled.
The add_seg = TRUE argument adds dashed guides at the
fold-change and significance thresholds.
You usually do not want to build a separate data frame just to label
a handful of hits. label_top picks the most significant
genes for you, and label_dir controls the direction they
are drawn from.
# the 10 most significant genes overall
ggvolc(all_genes, label_top = 10, add_seg = TRUE, title = "Top 10 hits")
# the top 8 up- and top 8 down-regulated genes
ggvolc(all_genes, label_top = 8, label_dir = "each",
add_seg = TRUE, title = "Top 8 up + 8 down")
label_dir accepts "both" (the default),
"up", "down", and "each".
ggvolc() calls significance on the adjusted p-value
(padj) by default, which is the right cutoff for most DE
workflows. The y-axis and the significance line follow the choice, so
the plot stays internally consistent. Switch to the raw p-value with
sig_col = "pvalue".
ggvolc() is also robust to p-values of exactly
0 (which DESeq2 and edgeR can emit for the strongest
genes): rather than dropping those points, their -log10
value is capped so they stay near the top of the plot.
Point size can encode either the fold change or the significance.
genes_table() composes the volcano plot with a
gt table of the highlighted genes using
patchwork. The result is a single object you can save with
ggplot2::ggsave().
To tabulate the most significant genes without curating them
yourself, pass the full DE table and a top_n. Genes are
ranked by sig_col (padj by default, falling
back to pvalue); dir = "each" takes the top N
up- and down-regulated genes.
p_top <- ggvolc(all_genes, label_top = 10, add_seg = TRUE)
genes_table(p_top, all_genes, top_n = 10)
Column names from all three pipelines are detected and mapped
automatically, so you can pass their output straight in. Gene
identifiers held in the row names (as edgeR and limma often do) are
promoted to a genes column for you.
The package ships edger_genes, an example
topTags()-style table with the gene identifiers in the row
names, so you can try this right away:
data(edger_genes)
head(edger_genes, 3)
#> logFC logCPM PValue FDR
#> GCR1 2.244064 12.8143 4.434241e-29 2.153711e-25
#> OPI10 -2.257454 9.9807 4.880607e-27 1.185255e-23
#> AGA2 3.829474 7.9665 4.143136e-26 6.707736e-23
ggvolc(edger_genes, label_top = 8, add_seg = TRUE, title = "edgeR input")
| Pipeline | Fold change | p-value | adjusted p | expression |
|---|---|---|---|---|
DESeq2 (results()) |
log2FoldChange |
pvalue |
padj |
baseMean |
edgeR (topTags()) |
logFC |
PValue |
FDR |
logCPM |
limma (topTable()) |
logFC |
P.Value |
adj.P.Val |
AveExpr |
You can also skip the conversion and pass the result objects
directly: DESeqResults or DESeqDataSet
(DESeq2), TopTags, DGELRT or
DGEExact (edgeR), and MArrayLM (limma, choose
the coefficient with coef).
ggvolc(DESeq2::results(dds))
ggvolc(edgeR::glmQLFTest(fit, coef = 2))
ggvolc(limma::eBayes(fit), coef = "treatment")de_table() performs the same conversion if you want the
data frame itself.
palette selects a color preset.
ggvolc_palettes() lists them, including colorblind-safe
options and presets from the ltc palettes.
ggvolc_palettes()
#> palette down not_sig up colorblind_safe source
#> 1 default #00798c #808080 #d1495b FALSE ltc: minou
#> 2 okabe_ito #0072B2 #999999 #D55E00 TRUE Okabe & Ito (2008)
#> 3 rdbu #0571b0 #999999 #ca0020 TRUE ColorBrewer RdBu
#> 4 pantone23 #1F2C43 #B8B8B8 #FFB000 TRUE ltc: pantone23
#> 5 seafarer #013D5A #BDBDBD #E4A25B TRUE ltc: seafarer
#> 6 midnight #075056 #B0B0B0 #FF5B04 FALSE ltc: midnight
#> 7 maya #3d5a80 #A9A9A9 #ee6c4d FALSE ltc: maya
#> 8 trio1 #0E7175 #B0B0B0 #FD7901 FALSE ltc: trio1
ggvolc(all_genes, label_top = 8, label_dir = "each", palette = "okabe_ito")
genes_table() reuses the plot’s colors, darkening light
ones so the table text stays readable.
# color named gene sets
ribosomal <- grep("^RP[LS][0-9]", all_genes$genes, value = TRUE)
heat_shock <- grep("^HSP[0-9]|^SSA[0-9]", all_genes$genes, value = TRUE)
ggvolc(all_genes, highlight = list(Ribosomal = ribosomal, `Heat shock` = heat_shock))
To compare contrasts, pass a named list of results.
stress_timecourse holds real DESeq2 results from a fission
yeast stress time course (Leong et al. 2014):
flip = TRUE puts the fold change on the y-axis instead,
which lines it up across the panels:
Each panel compares one time point with 0 minutes. The x-axis (fold
change) is shared; each panel uses its own y-axis height, and
facet_scales = "fixed" shares that too.
A named vector sets the colors directly, and
up_reg_color, down_reg_color and
not_sig_color still override single colors:
With the optional ggiraph
package installed, interactive = TRUE returns a widget
where hovering a point reveals the gene name and its statistics.
See the live, hover-to-inspect version on the package website.