Volcano Plots in RNA-seq: How to Build and Correctly Read Differential Expression Results

A volcano plot is the standard first-look figure for an RNA-seq differential expression result: one point per gene, plotted so that genes with the largest and most confidently supported changes sit toward the upper corners and everything else collapses toward the middle. It is quick to read, which is exactly what makes it easy to misread. The two mistakes that do the most damage are using the wrong p-value on the y-axis and picking thresholds that look reasonable but were never actually derived from the data. Both are avoidable once you understand what is actually being plotted.

What a volcano plot actually plots

Each point is one tested gene (or transcript). The axes are:

• x-axis: log₂ fold change (log2FC) — the estimated effect size, i.e. how much a gene’s expression differs between conditions on a log₂ scale. A value of 1 means roughly a 2-fold increase; -1 means roughly a 2-fold decrease. This comes straight from the differential expression model’s coefficient for your comparison of interest.
• y-axis: -log₁₀(adjusted p-value) — the statistical confidence that the fold change is real and not noise, transformed so that stronger evidence plots higher. A gene with an adjusted p-value of 0.01 sits at y = 2; one at 0.0001 sits at y = 4.

The genes worth following up sit in the upper-left and upper-right corners: large effect size and strong statistical support. Everything is measured against the model’s own estimate of noise, which is exactly why the choice of which p-value goes on that y-axis is not a cosmetic detail — see the next section.

The single most common and serious error: raw p-values instead of adjusted p-values

An RNA-seq differential expression call is not one test — it is one test per gene, typically 15,000-25,000 of them in a single comparison. At a raw significance threshold of 0.05 and no correction, roughly 5% of genes with truly zero effect will still cross that line by chance alone. Test 20,000 genes and that is on the order of 1,000 false positives, before a single real signal is counted. This is exactly why differential expression tools apply a multiple-testing correction — almost universally the Benjamini-Hochberg (BH) procedure, which controls the false discovery rate (FDR) rather than the per-gene error rate.

DESeq2’s results() function reports this BH-adjusted value in a column literally named padj, alongside the raw, uncorrected pvalue column. edgeR’s topTags() reports the same BH-adjusted quantity in a column named FDR. Both packages compute the raw p-value and the adjusted one side by side specifically so you cannot accidentally reach for the wrong one — but a volcano plot built by hand, or copied from a template, will silently transform whichever column you hand it. Wire up pvalue instead of padj by habit or by column-order mistake, and the plot looks identical in shape but is now showing hundreds of genes as “significant” that would not survive correction. It is one of the easiest RNA-seq mistakes to make and one of the hardest to catch by eye, because a volcano plot built on raw p-values is not obviously wrong — it is just wrong.

The fix is procedural, not statistical: before plotting, confirm which column feeds the y-axis, by name, every time. If you are scripting this, hard-code the column name (padj, FDR, or whatever your pipeline’s equivalent is) rather than a positional index, so a reordered results table cannot swap it silently.

Choosing defensible thresholds, not arbitrary round numbers

The conventional volcano plot cutoff — adjusted p-value below 0.05 and |log2FC| above 1 (i.e. 2-fold) — is a convention, not a derivation. Both halves deserve a second look before you adopt them unexamined:

• The significance threshold should reflect what you can afford downstream, not just habit. DESeq2’s own results() function defaults its internal independent-filtering optimization to an adjusted-p-value cutoff of 0.1 via its alpha argument — not the 0.05 many analysts assume by default — specifically because independent filtering is tuned to maximize the number of genes passing at whatever FDR you actually intend to use. If you plan to validate hits individually and each false positive costs real bench time, a tighter FDR (1% or lower) is defensible; if the plot is a screen feeding a much larger downstream filter, a looser one may be fine. State the number you used and why, rather than defaulting to 0.05 because it is familiar.
• The fold-change threshold is more often misapplied than the significance one. Filtering a gene list post hoc by “adjusted p < 0.05 AND |log2FC| > 1” treats the two as independent screens, but they are not measuring independent things — a gene with a huge fold change and very few reads can still fail on significance, and a gene with a tiny but extremely consistent fold change can pass significance easily at high sequencing depth. The more defensible approach is to test against the fold-change threshold directly, inside the model, rather than filtering afterward: DESeq2’s results() accepts an lfcThreshold argument that runs a modified Wald test against that boundary rather than against zero, and edgeR provides the equivalent TREAT approach via glmTreat() (McCarthy & Smyth, 2009). Both report a p-value that already accounts for the threshold you care about, instead of asking you to intersect two separately-computed columns.

Either way, pick the numbers before looking at how many genes they keep — adjusting a threshold until the gene list “looks right” is p-hacking by another name, just applied to fold change instead of p-value.

Building a volcano plot from DESeq2 or edgeR output

For visualization and ranking specifically (not for the significance test itself), shrink the fold-change estimates first. Raw log2FC estimates for low-count genes are noisy and inflate outward on their own, producing a characteristic flared base of low-count genes with enormous fold changes and no significance — that flare is a display artifact, not biology. DESeq2’s vignette addresses this directly, under “Log fold change shrinkage for visualization and ranking,” recommending lfcShrink() with the apeglm estimator for exactly this purpose — use the shrunken estimate for the x-axis position, but the unshrunken model’s p-values for the significance test itself.

A minimal, correct R example using DESeq2 and ggplot2:

library(DESeq2)
library(ggplot2)

# res <- results(dds, alpha = 0.05)  # unshrunken, for the significance test
res_shrunk <- lfcShrink(dds, coef = 'condition_treated_vs_control', type = 'apeglm')

df <- as.data.frame(res_shrunk)
df$sig <- with(df, !is.na(padj) & padj < 0.05 & abs(log2FoldChange) > 1)

ggplot(df, aes(x = log2FoldChange, y = -log10(padj), color = sig)) +
  geom_point(alpha = 0.6, size = 1, na.rm = TRUE) +
  scale_color_manual(values = c('grey70', 'firebrick'), guide = 'none') +
  geom_vline(xintercept = c(-1, 1), linetype = 'dashed', color = 'grey40') +
  geom_hline(yintercept = -log10(0.05), linetype = 'dashed', color = 'grey40') +
  labs(x = 'log2 fold change (shrunken)', y = expression(-log[10]~'adjusted p-value')) +
  theme_minimal()

Two details worth keeping even outside this exact example: filter or handle NA values in padj before plotting (DESeq2 sets padj to NA for genes removed by independent filtering or flagged as outliers — plotting them as zero or dropping them silently are both wrong), and use -log10() only after confirming the column is the adjusted p-value, per the section above. If you would rather not hand-roll the plot, the Bioconductor package EnhancedVolcano wraps this same logic — correct axis choice, threshold lines, labeling — with sensible defaults, and takes a standard DESeq2/edgeR results data frame directly.

How to correctly read a volcano plot

Three misreadings account for most of the wrong conclusions drawn from this figure:

Far right or left, but not high up, is not significant

A point sitting at the extreme edge of the x-axis draws the eye immediately — it looks like the biggest hit on the plot. If that same point sits low on the y-axis, it is not statistically significant, regardless of how large the fold change looks. This is the single most common misreading of a volcano plot: treating horizontal position alone as evidence. A large estimated fold change with weak statistical support usually means the gene had low counts, high variance across replicates, or too few replicates for the model to distinguish a real effect from noise — and unshrunk estimates make this worse, which is exactly the flare artifact described above. The gene may still be biologically interesting, but the plot itself has not told you it changed; more replicates or a targeted follow-up assay would be needed to say that.

Near the center but very high up is real, but may be small

The mirror-image error: a gene near x = 0 but far up the y-axis is highly significant but has a small effect size. At sufficient sequencing depth and replicate count, even a modest, consistent fold change can produce an extremely small adjusted p-value. Statistically real does not automatically mean biologically important — a 10% change detected with overwhelming confidence is still a 10% change. This is precisely the case the fold-change-threshold discussion above exists to handle: deciding, before you look, how large an effect has to be to matter for your question.

Shape of the point cloud is informative, not decorative

A genuinely bimodal or heavily right-skewed cloud, or a hard “wall” of points at one exact -log10(p) value, is usually a data or pipeline artifact rather than biology — a batch effect not accounted for in the model’s design formula, an outlier sample inflating dispersion estimates, or (the wall pattern specifically) p-values that were floored or ties introduced by an upstream filtering step. Treat an unusual overall shape as a QC prompt to go back to the model, not as a plotting quirk to ignore.

Volcano plot vs. MA plot

The volcano plot’s closest relative is the MA plot, which puts log2 fold change on the y-axis against mean normalized count (not significance) on the x-axis. DESeq2’s own plotMA() function colors points by adjusted p-value rather than plotting it as an axis. The two are complementary, not interchangeable: a volcano plot ranks genes by strength of statistical evidence, while an MA plot shows whether fold-change estimates are shrinking appropriately as counts get smaller — genuinely low-count, low-confidence genes should visibly compress toward zero fold change in a well-shrunk MA plot. If a volcano plot shows large fold changes surviving at low significance, the MA plot is often the faster way to see whether that is a count-depth problem across the whole dataset or isolated to a handful of genes.

Frequently asked questions

Should I use padj or FDR on the y-axis if my pipeline reports both a raw and an adjusted value under different names?

Use whichever column your tool documents as the multiple-testing-corrected value — padj in DESeq2, FDR in edgeR’s topTags(), or the equivalent in another package. Confirm this from the tool’s own documentation for the exact function version you are using rather than assuming a column position or name carries over between tools.

Why do some genes have no point on the plot at all?

Genes with an NA adjusted p-value — because they were removed by independent filtering (very low mean count, unlikely to reach significance at any fold change) or flagged as count outliers by Cook’s distance — have no defined y-coordinate and should be excluded from the plot rather than plotted at y = 0, which would misrepresent them as tested-and-not-significant when they were never meaningfully testable.

Is 2-fold (log2FC of 1) always the right cutoff?

No. It is a common convention, not a rule. The right cutoff depends on what fold change is biologically meaningful for your system and what you can afford to validate; testing directly against a chosen threshold with lfcThreshold (DESeq2) or glmTreat() (edgeR) is more defensible than filtering an arbitrary round number after the fact.

Can I make a volcano plot from limma-voom or another tool’s output?

Yes — the plot only needs a fold-change column and an adjusted-p-value column, and limma’s topTable() reports both (logFC and adj.P.Val, again BH by default). The axis logic and the misreadings above apply identically regardless of which tool produced the numbers.

My plot has a symmetric “flare” of extreme fold changes near the bottom. Is that a problem?

Usually yes, if you plotted unshrunken log2 fold changes. Low-count genes produce noisy, inflated fold-change estimates in either direction, which is exactly what shrinkage (lfcShrink()/apeglm in DESeq2) is designed to pull back toward zero for genes the data cannot support a large estimate for. See the building section above.

References

Love, M.I., Huber, W., Anders, S. (2014). “Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.” Genome Biology 15:550. The DESeq2 Bioconductor vignette documents the padj/pvalue distinction, the default alpha argument, and the lfcShrink()/apeglm shrinkage workflow cited above. McCarthy, D.J., Smyth, G.K. (2009). “Testing significance relative to a fold-change threshold is a TREAT.” Bioinformatics 25(6):765-771 — the method implemented in edgeR’s glmTreat(). For the underlying statistics of differential expression modeling, thresholding, and replicate requirements referenced throughout this guide, see CASRAI’s differential gene expression analysis guide and the broader RNA-seq experimental design and analysis guide; for the multiple-testing statistics behind the y-axis specifically, see Bonferroni correction and multiple comparisons, q-values and the positive false discovery rate, and the p-value dictionary entry. For how a p-value differs from a confidence interval more generally, see p-value vs confidence interval.