Volcano Plots for RNA-seq: Fold Change, P-Value Cutoffs, and Labeling Key Genes

You ran DESeq2 on your RNA-seq counts and have a results table with 20,000 rows. Most genes are noise; a few hundred are biology. A volcano plot is how you triage them at a glance — effect size on x, significance on y, with the genes worth labeling sitting in the upper corners.

This guide covers the choices that actually matter: which fold-change column to plot, which p-value to use (raw vs. adjusted, and why it matters), how to set thresholds without p-hacking, and how to label a crowded plot without it turning into a wall of text. Code is in R using ggplot2 and EnhancedVolcano, but the principles apply to any plotting tool.

What a Volcano Plot Shows

The x-axis is log₂ fold change (log2FC) — positive = upregulated in treatment, negative = downregulated. Log₂ is the convention because a value of 1 means “doubled” and −1 means “halved,” symmetric around zero. Linear fold change would compress downregulation into 0–1 and stretch upregulation to infinity, which is unreadable.

The y-axis is −log₁₀(p), so smaller p-values plot higher up. A gene with adjusted p = 10⁻⁶ sits at y = 6; a gene with p = 0.05 sits at y = 1.3. The transformation matters because raw p-values cluster between 0 and 1 and would be uninformative on a linear axis. The visual shape — two “wings” rising from the center — gives the plot its name.

Genes in the upper-left and upper-right corners are the candidates: large effect size AND high significance. Genes at the bottom (low significance) or near x = 0 (no effect) are not interesting regardless of which corner of the table they occupy.

Which P-Value to Plot: Raw or Adjusted

This is the most important decision and the one most tutorials get wrong. RNA-seq tests 15,000–25,000 genes simultaneously. Raw p-values are not interpretable at that scale — with 20,000 tests at α = 0.05, you expect 1,000 false positives by chance alone.

Plot the adjusted p-value (DESeq2's padj, edgeR's FDR, limma's adj.P.Val) on the y-axis. These columns are Benjamini-Hochberg-adjusted by default in all three tools, controlling the false discovery rate rather than the family-wise error rate. The interpretation is “of the genes I call significant, what fraction are expected to be false positives?” — which is the question a biologist asks. For more on why FDR is the right choice for genome-scale testing instead of Bonferroni, see Tukey, Bonferroni, and FDR compared.

Setting Thresholds Without P-Hacking

The two conventional cutoffs are |log2FC| ≥ 1 (i.e., 2-fold change in either direction) and padj < 0.05. They're conventions, not rules — and they should be chosen before looking at your data, then justified in the methods section.

A few sharper alternatives to the defaults:

• Stricter FDR for exploratory work: padj < 0.01 if you plan to validate ~50 genes by qPCR — tighter cutoff buys you a smaller, higher-confidence shortlist.
• Effect-size threshold by biology, not convention: in a tightly regulated pathway, |log2FC| ≥ 0.5 (1.4-fold) may be the meaningful threshold; in a knockout vs. WT comparison, |log2FC| ≥ 2 (4-fold) is more appropriate. Don't import the |log2FC| ≥ 1 default from a paper studying a different system.
• Shrunken log fold changes for low-count genes: DESeq2's lfcShrink() with type = "apeglm" shrinks fold-change estimates for low-count genes toward zero, removing the “wings” of high-FC, low-confidence genes near the bottom of the plot. Use shrunken estimates for visualization; use raw estimates for downstream pathway analysis.

The threshold is a methodological choice that goes in your figure legend, not a knob you twist after seeing results. Twisting it until your favorite gene crosses the line is p-hacking by another name.

Building the Plot in R: ggplot2 Approach

The minimum viable volcano plot, assuming a DESeq2 results object res:

library(ggplot2)
library(dplyr)

df <- as.data.frame(res) %>%
  filter(!is.na(padj)) %>%
  mutate(
    significant = padj < 0.05 & abs(log2FoldChange) >= 1,
    direction = case_when(
      significant & log2FoldChange > 0 ~ "Up",
      significant & log2FoldChange < 0 ~ "Down",
      TRUE ~ "NS"
    )
  )

ggplot(df, aes(x = log2FoldChange, y = -log10(padj), color = direction)) +
  geom_point(alpha = 0.6, size = 1.2) +
  scale_color_manual(values = c("Up" = "#D55E00", "Down" = "#0072B2", "NS" = "#BBBBBB")) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "grey50") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey50") +
  labs(x = "log₂ fold change", y = "−log₁₀ adjusted p-value", color = NULL) +
  theme_classic(base_size = 9) +
  theme(legend.position = "top")

The filter(!is.na(padj)) step matters: DESeq2 sets padj to NA for genes that fail independent filtering (low mean count, extreme outliers). Plotting NAs as −log₁₀(NA) produces silent gaps. The colorblind-safe Okabe-Ito palette (orange/blue/grey) is the recommended default for Nature and PLOS figures.

Labeling Genes Without Overplotting

The default approach — geom_text() — produces unreadable plots the moment you have more than ~10 labels. Use ggrepel to avoid overlap, and label selectively:

library(ggrepel)

top_genes <- df %>%
  filter(significant) %>%
  arrange(padj) %>%
  slice_head(n = 15)

# Add geom_text_repel layer to the plot from the previous block
... +
  geom_text_repel(
    data = top_genes,
    aes(label = gene_symbol),
    size = 2.5,
    box.padding = 0.4,
    max.overlaps = 20
  )

EnhancedVolcano: The One-Liner Alternative

If you want a publication-ready volcano plot in three lines, the EnhancedVolcano Bioconductor package handles thresholds, coloring, and label-repulsion in one call:

library(EnhancedVolcano)

EnhancedVolcano(res,
  lab = rownames(res),
  x = "log2FoldChange",
  y = "padj",
  pCutoff = 0.05,
  FCcutoff = 1,
  title = "Treatment vs. Control",
  subtitle = "DESeq2 differential expression"
)

EnhancedVolcano uses ggrepel internally and applies the standard four-color scheme (NS, log2FC pass only, p pass only, both pass). It's the right choice for a quick figure; ggplot2 from scratch is the right choice when you need precise control over typography, color, or annotation for journal submission.

Volcano vs. MA Plot: When to Use Which

An MA plot puts log₂ mean expression on x and log₂ fold change on y. It's the diagnostic plot — you check for a horizontal smear around y = 0 (good normalization) and you spot trumpet-shaped fold-change inflation at low counts (a sign that shrinkage is needed).

A volcano plot is the selection plot — you read it to pick candidates. They answer different questions and you should usually make both. MA plot first (is the analysis well-normalized?), volcano plot second (which genes are interesting?). Don't substitute one for the other.

What Reviewers Catch

• Raw p-values on the y-axis. Use adjusted p-values for any genome-scale experiment.
• No threshold lines. Dashed lines at the chosen log2FC and padj cutoffs let the reader see where your filter is.
• Missing axis transformation note. The figure legend should say “y-axis: −log₁₀(adjusted p-value).” Don't assume the reader infers it.
• Colorblind-unfriendly palette. Red-green is the worst possible choice for ~8% of male readers. Use Okabe-Ito or viridis.
• Inconsistent labeling. Either label by rank (top N by padj) or by curated list, not a mix. State which in the legend.

Once the plot is right, the methods section sentence reads: “Differentially expressed genes were visualized as a volcano plot of log₂ fold change against −log₁₀ Benjamini-Hochberg-adjusted p-value, with thresholds at |log₂FC| ≥ 1 and padj < 0.05.” Specific, reproducible, defensible — the bar reviewers check.

For the upstream design decision — how many biological replicates you need to detect a given fold change at adequate power — the sample size and power calculator handles two-group comparisons and ANOVA designs. RNA-seq-specific power tools (PROPER, RNASeqPower) are more accurate for count data but the t-test approximation is a reasonable first pass for budget planning.