Applied - Differential Expression Volcano Plot
In one line: the standard DE figure — effect size against significance — built properly, with the details that separate a publishable figure from a default one.
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)isinf, which breaks the axis limits and silently drops points. Very small p-values do underflow to 0 in float64.padj, notpvalue— plotting raw p-values on a 20,000-gene test is not a defensible figure.- DESeq2 sets
padjto 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
def volcano(de, ax=None, fc_thresh=1.0, p_thresh=0.05, n_label=10,
palette=None, title=None):
"""Volcano plot. Returns the Axes so it composes into a panel."""
palette = palette or {"ns": "#CCCCCC", "up": "#C44E52", "down": "#4C72B0"}
if ax is None:
_, ax = plt.subplots(figsize=(4.5, 4))
d = de.dropna(subset=["padj"]).assign(
y=lambda x: -np.log10(x["padj"].clip(lower=1e-300)),
direction=lambda x: np.select(
[(x["padj"] < p_thresh) & (x["log2FoldChange"] > fc_thresh),
(x["padj"] < p_thresh) & (x["log2FoldChange"] < -fc_thresh)],
["up", "down"], default="ns"),
)
for cat in ["ns", "down", "up"]:
s = d[d["direction"] == cat]
ax.scatter(s["log2FoldChange"], s["y"], s=6, c=palette[cat],
alpha=0.6 if cat == "ns" else 0.9, edgecolors="none",
rasterized=True, label=f"{cat} ({len(s):,})")
ax.axhline(-np.log10(p_thresh), ls="--", lw=0.8, c="grey", zorder=0)
for v in (-fc_thresh, fc_thresh):
ax.axvline(v, ls="--", lw=0.8, c="grey", zorder=0)
if n_label:
top = pd.concat([d.query("direction=='up'").nlargest(n_label, "y"),
d.query("direction=='down'").nlargest(n_label, "y")])
for r in top.itertuples():
ax.annotate(r.gene, (r.log2FoldChange, r.y), fontsize=6,
xytext=(3, 3), textcoords="offset points")
ax.set(xlabel=r"log$_2$ fold change", ylabel=r"$-$log$_{10}$ adj. $p$",
title=title)
ax.spines[["top", "right"]].set_visible(False)
ax.legend(frameon=False, fontsize=7, markerscale=2)
return ax
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.
fig, axes = plt.subplots(1, 3, figsize=(13, 4), layout="constrained")
# MA plot — effect vs expression level; reveals normalisation problems
axes[0].scatter(np.log10(de["baseMean"] + 1), de["log2FoldChange"],
s=3, alpha=0.3, c=np.where(de["sig"], "crimson", "lightgrey"),
rasterized=True)
axes[0].axhline(0, c="grey", lw=0.8)
axes[0].set(xlabel=r"log$_{10}$ mean expression", ylabel=r"log$_2$ FC", title="MA")
# p-value histogram — the DE sanity check
axes[1].hist(de["pvalue"].dropna(), bins=50, edgecolor="white")
axes[1].set(xlabel=r"raw $p$", ylabel="genes", title="p-value distribution")
# volcano
volcano(de, ax=axes[2], title="Volcano")
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
padjwith 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