Applied - Differential Expression Volcano Plot

Core Principle: Volcano plots display differential gene expression by graphing effect size ($\log_2$ fold change) against statistical significance ($-\log_{10}$ adjusted $p$-value), highlighting biologically significant targets.

Uses: Matplotlib · Seaborn · Colormaps and Color · Matplotlib Artists · Saving Figures · DataFrame

The data

A DESeq2/edgeR-style result table: gene, baseMean, log2FoldChange, pvalue, padj.

import numpy as np, pandas as pd, matplotlib.pyplot as plt, seaborn as sns

de = pd.read_csv("deseq2_results.tsv", sep="\t")

de = de.assign(
    neglog10p=lambda d: -np.log10(d["padj"].clip(lower=1e-300)),
    sig=lambda d: (d["padj"] < 0.05) & (d["log2FoldChange"].abs() > 1),
    direction=lambda d: np.select(
        [(d["padj"] < 0.05) & (d["log2FoldChange"] > 1),
         (d["padj"] < 0.05) & (d["log2FoldChange"] < -1)],
        ["up", "down"], default="ns"),
)

Three details in those four lines:

• .clip(lower=1e-300): -log10(0) is inf, which breaks the axis limits and silently drops points. Very small p-values do underflow to 0 in float64.
• padj, not pvalue: plotting raw p-values on a 20,000-gene test is not a defensible figure.
• DESeq2 sets padj to NA for genes it filtered out. Do not fill those with 1: that makes untested genes look tested-and-not-significant. Drop them and say how many.

n_na = de["padj"].isna().sum()
de = de.dropna(subset=["padj"])
print(f"{n_na:,} genes excluded by independent filtering")

The plot

sns.set_theme(style="ticks", context="paper")
import matplotlib as mpl
mpl.rcParams.update({"pdf.fonttype": 42, "savefig.bbox": "tight"})

PALETTE = {"ns": "#CCCCCC", "up": "#C44E52", "down": "#4C72B0"}

fig, ax = plt.subplots(figsize=(4.5, 4))

for cat in ["ns", "down", "up"]:            # grey FIRST so hits draw on top
    sub = de[de["direction"] == cat]
    ax.scatter(sub["log2FoldChange"], sub["neglog10p"],
               s=6, c=PALETTE[cat], alpha=0.6 if cat == "ns" else 0.9,
               edgecolors="none", rasterized=True, label=f"{cat} ({len(sub):,})")

ax.axhline(-np.log10(0.05), ls="--", lw=0.8, c="grey", zorder=0)
ax.axvline(-1, ls="--", lw=0.8, c="grey", zorder=0)
ax.axvline(1, ls="--", lw=0.8, c="grey", zorder=0)

ax.set(xlabel=r"log$_2$ fold change",
       ylabel=r"$-$log$_{10}$ adjusted $p$")
ax.spines[["top", "right"]].set_visible(False)
ax.legend(frameon=False, fontsize=7, markerscale=2, loc="upper left")

rasterized=True is not optional. A vector PDF with 20,000 individually-stored scatter points is enormous and crashes Illustrator; rasterizing that one Artist keeps the text and axes as vectors. See Saving Figures.

Draw order matters. Plotting "ns" first means the significant points land on top rather than being buried under grey.

Labelling the top genes

top = pd.concat([
    de.query("direction == 'up'").nlargest(8, "neglog10p"),
    de.query("direction == 'down'").nlargest(8, "neglog10p"),
])

try:
    from adjustText import adjust_text
    texts = [ax.text(r.log2FoldChange, r.neglog10p, r.gene, fontsize=6)
             for r in top.itertuples()]
    adjust_text(texts, ax=ax,
                arrowprops=dict(arrowstyle="-", lw=0.4, color="grey"))
except ImportError:
    for r in top.itertuples():
        ax.annotate(r.gene, (r.log2FoldChange, r.neglog10p), fontsize=6,
                    xytext=(3, 3), textcoords="offset points")

adjustText repels overlapping labels iteratively. Without it, your 16 labels become an unreadable pile in exactly the region you most want to read.

Label ~10–20 genes, not 100. The point of labelling is to direct attention; labelling everything directs it nowhere.

Capping the axes honestly

XCAP, YCAP = 6, 50
de["x_plot"] = de["log2FoldChange"].clip(-XCAP, XCAP)
de["y_plot"] = de["neglog10p"].clip(upper=YCAP)
de["capped"] = (de["log2FoldChange"].abs() > XCAP) | (de["neglog10p"] > YCAP)

# draw capped points as triangles so it is visible
capped = de[de["capped"]]
ax.scatter(capped["x_plot"], capped["y_plot"], marker="^", s=20,
           c=[PALETTE[d] for d in capped["direction"]], edgecolors="black", lw=0.3)

A few extreme genes otherwise compress everything else into a horizontal smear. Capping is legitimate; capping silently is not. Use a distinct marker and state it in the caption.

Saving

fig.savefig("volcano.pdf", dpi=300, bbox_inches="tight")
fig.savefig("volcano.png", dpi=300, bbox_inches="tight")
plt.close(fig)

Both formats every time: PDF for the manuscript, PNG for the slide.

A reusable function

The ax=None pattern is what makes it composable into a multi-panel figure:

fig, axes = plt.subplots(1, 3, figsize=(13, 4), sharey=True, layout="constrained")
for ax, (name, sub) in zip(axes, contrasts.items()):
    volcano(sub, ax=ax, title=name)

See pyplot vs Object-Oriented API.

The companion plots

A volcano alone is not enough for a methods-complete figure.

<Axes: title={'center': 'Volcano'}, xlabel='log$_2$ fold change', ylabel='$-$log$_{10}$ adj. $p$'>

The p-value histogram is the most diagnostic plot in DE analysis and the most often omitted. Expected shape: uniform, with a spike near 0 for the true positives.

• Uniform with a spike at 0: healthy.
• Uniform with no spike: no signal.
• U-shaped or a spike near 1: model misspecification, usually unmodelled covariates or variance overestimation.
• Bumpy / multimodal: often a batch effect or a discreteness artefact.

Look at this before you look at the volcano. A beautiful volcano built on a pathological p-value distribution is a beautiful wrong figure.

Common mistakes

• -log10(0) = inf. Clip first.
• Raw p-values instead of adjusted.
• Filling DESeq2 NA padj with 1.
• Not rasterizing, producing a 200 MB PDF.
• Labelling too many genes.
• Silent axis capping.
• Significant points buried under grey. Draw order.
• Only fold-change thresholds, no significance (or vice versa).
• Skipping the p-value histogram.
• Red/green colouring: the worst possible choice for colourblind readers. See Colormaps and Color.

See also

Common Plot Types · Matplotlib Artists · Colormaps and Color · Saving Figures · Applied - Expression Matrices with NumPy · Applied - Heatmaps and Clustermaps · Statistical Estimation in seaborn

Exercise

seaborn is not installed in this sandbox, so build the volcano with matplotlib alone, and assert on the numbers the Axes carries rather than on pixels. Use the Agg backend and never call show.

Write volcano(text), returning a list of output lines. The first input line holds three whitespace-separated values: the absolute log2 fold-change threshold, the adjusted p-value threshold, and the number of genes to label. The rest is a tab-separated table with the header gene, baseMean, log2FoldChange, pvalue, padj; a padj of NA means the gene was filtered by independent filtering and must be dropped before anything else. Classify each remaining gene as up when the fold change is at or above the threshold and padj is below the p threshold, down when it is at or below the negative threshold and equally significant, and ns otherwise. Plot each class as its own scatter, set the x limits symmetrically to plus and minus round(max(abs(log2FoldChange)) + 0.5, 2), and label the most significant genes, ranked by padj ascending then absolute fold change descending then gene name ascending.

Print plotted=<n> up=<n> down=<n> ns=<n>, then xlim=<low>,<high> reading the limits back from the Axes and rounding each to two decimals, then ymax=<v> for the largest plotted value of -log10(padj) rounded to two decimals, then labels=<g1,g2,...> in the ranked order, or labels=none when no gene is labelled. Clip padj at 1e-300 before taking the logarithm so a zero cannot produce infinity. If every gene was dropped, print the first line with all counts zero and then no genes to plot, and nothing else.