bio-data-visualization-volcano-and-ma-plots

name	bio-data-visualization-volcano-and-ma-plots
description	Build volcano and MA plots from differential-expression / association results with LFC shrinkage, FDR-adjusted thresholds, sensible label placement, and axis-truncation conventions. Covers EnhancedVolcano, ggplot2, matplotlib, and the apeglm/ashr/normal shrinkage decision. Use when visualizing differential-expression results (RNA-seq, ChIP-seq, ATAC-seq, proteomics) or any per-feature effect-size + p-value table.
origin	openai4s
category	bioskills/data-visualization
metadata	{"tool_type":"mixed","primary_tool":"ggplot2","third_party":{"name":"GPTomics/bioSkills","repository":"https://github.com/GPTomics/bioSkills","commit":"d91ed3d563019e649dc854c56ccd62551359488a","license":"MIT"}}

Version Compatibility

Reference examples tested with: DESeq2 1.42+, EnhancedVolcano 1.20+, ggplot2 3.5+, ggrepel 0.9.5+, matplotlib 3.8+, numpy 1.26+, adjustText 1.1+, apeglm 1.28+, ashr 2.2+.

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Volcano and MA Plots

"Plot differential-expression results" -> Place per-feature shrunken effect estimate on the x-axis and a significance measure (-log10 padj or -log10 p) on the y-axis. The decision space spans which effect estimate (raw vs shrunken), which significance measure (raw p vs adjusted vs s-value), how to encode categories (color by direction, not by gradient), how to label (top-N is rarely informative), and how to handle the tail (extreme p compresses the plot).

• R: EnhancedVolcano::EnhancedVolcano(), ggplot2 + ggrepel, DESeq2::plotMA()
• Python: matplotlib.scatter with adjustText, sanbomics.tools.volcano, custom seaborn

The Single Most Important Modern Insight -- Plot Shrunken LFC

Raw log2 fold change from DESeq2 / edgeR is the maximum-likelihood estimate and inflates wildly at low counts. A gene with 2 vs 0 reads gets log2FC = Inf; one with 4 vs 1 gets log2FC = 2 with a huge standard error. A naive volcano labels these as "top hits" purely because they have extreme estimates, not because they have a real signal.

DESeq2::lfcShrink() applies an empirical-Bayes prior to pull noisy low-count LFCs toward zero while leaving well-estimated genes essentially untouched. Since DESeq2 v1.28, the default shrinkage is type='apeglm' (Zhu, Ibrahim, Love 2019 Bioinformatics 35:2084) which uses a Cauchy prior — heavy enough to preserve large real effects, sharp enough at zero to deflate noise. Plot the shrunken LFC. The unshrunken LFC is a misleading effect estimate for ranking, labeling, or thresholding.

A complementary modern alternative is the s-value (Stephens 2017 Biostatistics 18:275): the local false sign rate — probability that the sign of the effect is wrong. s-values rank genes by how confident the sign is, which is what a volcano plot is implicitly trying to communicate. Where padj answers "is the effect non-zero," s answers "do we know which direction."

Shrinkage Method Selection

Method	Prior	Best for	Fails when
apeglm (Zhu 2019)	Cauchy	Default; preserves large LFCs, deflates noise	Requires coef=; no support for contrast=
ashr (Stephens 2017)	Mixture of normals	Comparisons requiring contrast=; supports s-values	Slightly more aggressive shrinkage of medium effects
normal (DESeq2 original)	Zero-centered normal	Legacy reproducibility only	Over-shrinks large effects; deprecated in current DESeq2 vignette
Unshrunken MLE	None	NEVER for volcano/MA plots	Low-count genes dominate the tails with no real signal

library(DESeq2)
dds <- DESeq(dds)
res_apeglm <- lfcShrink(dds, coef = 'condition_treated_vs_control', type = 'apeglm')
res_ashr <- lfcShrink(dds, contrast = c('condition', 'treated', 'control'), type = 'ashr')
# ashr also returns svalue column (Stephens 2017 local false sign rate)

For edgeR users: topTags() already provides moderated p-values but does not shrink LFC. Use glmTreat() for a moderated test against a non-zero LFC threshold; this is the edgeR equivalent of the shrunken-LFC philosophy.

Decision Tree by Scenario

Scenario	Recommended	Why
Bulk RNA-seq with DESeq2	lfcShrink(type='apeglm') + plot on padj threshold	apeglm is the default since DESeq2 v1.28
Non-default contrast required	lfcShrink(type='ashr')	apeglm requires coef=, ashr accepts contrast=
Single-cell pseudobulk DE	lfcShrink(type='apeglm') + raster the plot	sc datasets often >20k genes; raster prevents PDF lag
Proteomics (limma/MSstats)	Already moderated; plot logFC vs adj.P.Val directly	limma's empirical-Bayes already shrinks
Microarray (limma)	topTable() adj.P.Val + logFC	Same — limma is the original shrunken-LFC method
ATAC/ChIP differential peaks	lfcShrink(type='apeglm') on DESeq2/DiffBind	Treat peaks as features identically to genes
Want to rank by sign-confidence	ashr's svalue, NOT padj	Stephens 2017 — sign-aware ranking
Many comparisons in one figure	Faceted MA plot with shared y-axis	MA scales better than volcano for >6 panels

Volcano with ggplot2 + ggrepel

Goal: Plot shrunken LFC vs -log10 p-value, color by significance class, label genes of interest with non-overlapping repulsion.

Approach: Compute a categorical significance variable from padj AND |LFC| thresholds; pre-select labels (genes of interest OR top-N by combined rank) before plotting; use ggrepel::geom_text_repel with max.overlaps = Inf to guarantee every selected label appears.

library(ggplot2)
library(ggrepel)
library(dplyr)

volcano_plot <- function(res, fdr = 0.05, lfc_threshold = 1, label_genes = NULL, top_n = 10) {
    res <- as.data.frame(res) %>%
        tibble::rownames_to_column('gene') %>%
        mutate(
            significance = case_when(
                is.na(padj) ~ 'NS',
                padj < fdr & log2FoldChange > lfc_threshold ~ 'Up',
                padj < fdr & log2FoldChange < -lfc_threshold ~ 'Down',
                  
            
            neg_log10_p  log10pvalue
        

     label_genes 
        label_genes  res 
            filtersignificance   
            mutaterank_score  log10pvalue  log2FoldChange 
            arrangedescrank_score 
            headtop_n 
            pullgene
    
    reslabel  ifelseresgene  label_genes resgene 

    okabe_ito  Up   Down   NS  

    ggplotres aeslog2FoldChange neg_log10_p color  significance 
        geom_pointalpha   size   
        scale_color_manualvalues  okabe_ito name   
        geom_vlinexintercept  lfc_threshold lfc_threshold
                   linetype   color   linewidth   
        geom_hlineyintercept  log10fdr linetype  
                   color   linewidth   
        geom_text_repelaeslabel  label color   size  
                        max.overlaps   box.padding   segment.size  
                        min.segment.length   
        labsx  
             y  italicp 
        theme_classicbase_size   
        themepanel.grid  element_blank

Key design choices encoded above:

• Colors from Okabe-Ito (Wong 2011 Nat Methods 8:441) — CVD-safe categorical palette
• max.overlaps = Inf because the ggrepel default of 10 silently drops labels with no error (see [[api_gotchas]])
• Threshold line on -log10(fdr) matches the adjusted p threshold; drawing it on raw p creates a meaningless line
• rank_score = -log10(pvalue) * abs(log2FoldChange) selects labels that are both significant AND have non-trivial effect; pure top-N-by-p selects high-count genes with tiny effects

EnhancedVolcano -- Production Use and Its Gotchas

library(EnhancedVolcano)
EnhancedVolcano(res,
    lab = rownames(res),
    x = 'log2FoldChange',
    y = 'padj',                    # use padj NOT pvalue for the threshold line
    pCutoff = 0.05,
    FCcutoff = 1,
    selectLab = c('TP53', 'MYC', 'BRCA1'),
    drawConnectors = TRUE,
    widthConnectors = 0.3,
    maxoverlapsConnectors = Inf,
    colAlpha = 0.6,
    pointSize = 1.5,
    labSize = 3,
    col = c('grey60', '#0072B2',  
    legendPosition  

Gotcha 1: selectLab filters through pCutoff AND FCcutoff. Genes explicitly listed but failing thresholds appear unlabeled, with no warning. This is the most common "why is the gene missing" failure. To force-label regardless of thresholds, pre-shrink the input so the genes pass the cutoff, or switch to a manual ggplot layer.

Gotcha 2: y = 'pvalue' vs y = 'padj'. Many tutorials use raw pvalue for the y-axis, then draw the threshold line at -log10(0.05) — that line corresponds to a raw p < 0.05, not an FDR. Use y = 'padj' so the threshold is meaningful.

Gotcha 3: Asymmetric x-limits hide the diverging null distribution. Use xlim = c(-max(abs(LFC)), max(abs(LFC))) so the plot is symmetric around zero.

MA Plot -- The Underused Diagnostic

The MA plot (log2-mean vs log2-fold-change) is the original RNA-seq diagnostic (Dudoit 2002 JASA). It exposes the abundance-dependent variance structure that the volcano hides:

• A fan-shaped MA plot with extreme LFCs concentrated at low baseMean indicates inadequate shrinkage
• A horizontal "stripe" of significant genes at one LFC value indicates batch confound with treatment
• An asymmetric distribution (more Up than Down) at the low-count end indicates library-size normalization failure

library(DESeq2)
plotMA(res_apeglm, alpha = 0.05, ylim = c(-5, 5))
# alpha colors significant points; ylim clips for readability without losing the gene

import matplotlib.pyplot as plt
import numpy as np

def ma_plot(res, fdr=0.05, ax=None):
    if ax is None:
        fig, ax = plt.subplots(figsize=(6, 5))
    sig = (res['padj'] < fdr) & res['padj'].notna()
    ax.scatter(np.log10(res.loc[~sig, 'baseMean']), res.loc[~sig, 'log2FoldChange'],
               c='#999999', s=4, alpha=0.4, rasterized=True)
    ax.scatter(np.log10(res.loc[sig, 'baseMean']), res.loc[sig, 'log2FoldChange'],
               c='#D55E00', s=6, alpha=0.7, rasterized=True)
    ax.axhline(0, color='black', linewidth=0.5)
    ax.set_xlabel(r'$\log_{10}$ mean normalized count')
    ax.set_ylabel(r'$\log_2$ fold change (shrunken)')
    return ax

rasterized=True is critical for >5000 points — vector scatter creates 5MB+ PDFs that crash Illustrator.

Per-Method Failure Modes

Unshrunken LFC plotted as volcano

Trigger: Calling results(dds) and plotting log2FoldChange directly, without lfcShrink().

Mechanism: ML estimate has infinite-variance tails at low counts; one read difference produces log2FC = Inf.

Symptom: "Top hits" by |LFC| are all genes with baseMean < 5; biologically interesting genes with moderate LFC are hidden in the noise cloud.

Fix: lfcShrink(dds, coef=..., type='apeglm') for the default case; type='ashr' if contrast= is needed.

Raw p threshold line drawn on adjusted axis

Trigger: Drawing geom_hline(yintercept = -log10(0.05)) on a plot whose y-axis is -log10(padj).

Mechanism: padj < 0.05 corresponds to FDR < 5% control, NOT to raw p < 0.05. The drawn line is at the wrong y-value relative to the data.

Symptom: Visible "significant" points sit below the FDR line; the legend says "FDR < 0.05" but the line doesn't separate them correctly.

Fix: Be explicit: if y-axis is padj, use -log10(fdr) as the threshold line value AND label it "FDR threshold." If y-axis is raw p, an FDR threshold cannot be drawn as a horizontal line — the FDR threshold moves per gene.

Top-N-by-p selects low-effect-size hits

Trigger: head(arrange(res, pvalue), 20) to choose labels.

Mechanism: With large N, the smallest p-values belong to high-count, low-variance, biologically-boring genes (housekeeping). Effect size and statistical confidence are not the same thing.

Symptom: Labels are GAPDH, ACTB, B2M — never the gene that drives the biology.

Fix: Rank by -log10(p) * abs(log2FoldChange) (geometric average of the two axes), OR pre-specify labels of interest from prior knowledge.

ggrepel max.overlaps silently drops labels

Trigger: Default max.overlaps = 10; 30 genes labeled in code; only 10 render.

Mechanism: ggrepel emits a warning ("18 unlabeled data points (too many overlaps)") but no error. In a Quarto/Rmd render the warning is buried in the log.

Symptom: Reviewer asks "where is gene X?"; the label was specified in code but did not render.

Fix: geom_text_repel(..., max.overlaps = Inf) or options(ggrepel.max.overlaps = Inf) at the top of the script.

Extreme p-values compress the upper axis

Trigger: Genes with p = 1e-200 or smaller (common in cancer datasets) push the y-axis maximum to 200; all biologically meaningful genes pile up at the bottom.

Mechanism: -log10 expands the tail; one ultra-significant gene visually dominates.

Symptom: Volcano looks like an Eiffel Tower with most genes squished near y = 0-20.

Fix: Cap the y-axis (ylim = c(0, 50) and use coord_cartesian so capped points stay in the data but render at the edge) OR transform with sqrt(-log10(p)) to compress the tail OR split the y-axis with ggbreak.

lfcShrink(type='normal') on a modern DESeq2

Trigger: Following old tutorials that pre-date DESeq2 v1.28 when apeglm became the default.

Mechanism: The 'normal' prior over-shrinks large real effects toward zero.

Symptom: Volcano looks "too clean" — genuine 8-fold changes appear as 2-3 fold.

Fix: Use type='apeglm' (default since v1.28) or type='ashr'. Vignette removed 'normal' from recommendations.

EnhancedVolcano's selectLab filters by thresholds

Trigger: Listing genes in selectLab that have padj > pCutoff.

Mechanism: Source code applies pCutoff AND FCcutoff filter to selectLab membership; silently drops any that don't pass.

Symptom: Specific genes requested in selectLab do not appear; no warning.

Fix: Pre-shrink (so genes pass) or build the labeled subset manually with ggrepel and add as an annotation layer.

Reconciliation: When Methods Disagree

Pattern	Likely cause	Action
apeglm and ashr give different "top hits"	Different prior shapes; medium-effect, medium-count genes are most sensitive	Both are valid; pick one and document. apeglm is the DESeq2 default and the published recommendation
EnhancedVolcano shows fewer points than ggplot	EnhancedVolcano drops padj = NA (DESeq2 independent filtering)	Confirm by counting is.na(res$padj); to include them, set NA padj to 1 before plotting
Volcano has many "significant" genes but MA plot shows them all at low baseMean	Unshrunken LFC; the volcano is showing fold-change noise	Re-plot with shrunken LFC; the MA-plot fan is the diagnostic
Forest of horizontal stripes in MA at integer LFC	Pseudocount-induced quantization in low-count genes	Increase normalization-method aggressiveness OR filter low-count genes upstream
Half the genes have padj = NA	DESeq2 independent filtering (Bourgon-Gentleman-Huber 2010 PNAS) excluded them as low-mean	This is correct behavior; do NOT set independentFiltering = FALSE to hide it. Report the NA count

Operational rule: the volcano plot is read in this priority order — (1) is LFC shrunken? (2) is the y-axis padj or raw p? (3) does the threshold line match the axis? (4) are labels selected by combined rank? If any of these is wrong, the plot is misleading regardless of how attractive it looks.

Quantitative Thresholds

Threshold	Value	Source
Default LFC cutoff for "biologically relevant"	|log2FC| > 1 (2-fold)	Convention; sensitive analyses use 0.58 (1.5-fold) for subtle effects
Default FDR cutoff	padj < 0.05	Benjamini-Hochberg 1995 JRSS-B 57:289
Stricter cutoff for unbiased screens	padj < 0.01	Reduces false positives in unbiased genome-wide analyses
Relaxed cutoff for exploratory / hypothesis-generating	padj < 0.10 or 0.20	Acceptable for follow-up enrichment, NOT for "hits"
s-value cutoff (Stephens 2017)	s < 0.005 corresponds approximately to padj < 0.05	Stephens 2017 Biostatistics 18:275
Raster threshold	>5000 points	Vector PDF crashes Illustrator; raster scatter, keep axes vector
ggrepel max.overlaps	Set to Inf for publication	Default 10 silently drops labels
Volcano y-axis cap	-log10(p) > 50 typically warrants capping	Visual compression of biologically meaningful genes

Common Errors

Error / symptom	Cause	Solution
All "top hits" have baseMean < 5	Unshrunken LFC	lfcShrink(type='apeglm')
Threshold line doesn't separate colored from grey points	y = pvalue but threshold drawn at FDR	Switch y to padj OR redraw line at FDR-equivalent p
Labeled gene does not appear in EnhancedVolcano	selectLab filtered by pCutoff/FCcutoff	Pre-shrink or build labels manually
Volcano renders as a flat horizontal cloud	Extreme p (e.g., 1e-200) dominates y-axis	Cap with coord_cartesian(ylim = c(0, 50)) or use sqrt transform
PDF crashes Illustrator	Vector scatter of >10000 points	Set rasterized = TRUE in geom_point or matplotlib scatter
ggrepel labels 10 of 30 selected genes	Default max.overlaps = 10	geom_text_repel(max.overlaps = Inf)
Volcano "significant" gene count differs from DESeq2 summary	EnhancedVolcano drops padj = NA	Set NA padj to 1 or document the discrepancy
Up and Down counts asymmetric for a balanced experiment	Library-size normalization failure	Re-run DESeq2 with estimateSizeFactors(type='poscounts') for sparse data

Anticipated Reviewer Pushback

Pushback	Standard response
"Why is this LFC shrunken? Show me the unshrunken."	Shrunken LFC is the recommended estimate for ranking and visualization (Zhu 2019). Unshrunken LFC inflates at low counts and gives misleading rank. Unshrunken is available in the supplementary table
"Why padj < 0.05 not p < 0.05?"	padj controls FDR via Benjamini-Hochberg. Raw p < 0.05 across 20000 genes yields ~1000 false positives by chance; padj < 0.05 caps the expected false-positive rate at 5% of called hits
"Why are X gene and Y gene not labeled?"	Labels selected by combined rank (-log10(p) * |LFC|) or pre-specified gene list. List of all significant genes is in supplementary table T1
"The volcano looks too clean / too sparse."	Color encodes 3 categories (Up/Down/NS), not a gradient. Gradient encoding implies a continuous interpretation of significance which is invalid — significance is a threshold decision
"Why is the x-axis asymmetric?"	Asymmetric x-axis reflects the asymmetry of the data. If symmetry is preferred for visual interpretation, use xlim = c(-X, X) with X = max(|LFC|)

References

• Anders S, Huber W. 2010. Differential expression analysis for sequence count data. Genome Biol 11:R106.
• Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. JRSS-B 57:289-300.
• Bourgon R, Gentleman R, Huber W. 2010. Independent filtering increases detection power for high-throughput experiments. PNAS 107:9546-9551.
• Dudoit S, Yang YH, Callow MJ, Speed TP. 2002. Statistical methods for identifying differentially expressed genes in replicated cDNA microarray experiments. Stat Sin 12:111-139.
• Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15:550.
• Stephens M. 2017. False discovery rates: a new deal. Biostatistics 18(2):275-294. doi:10.1093/biostatistics/kxw041
• Wong B. 2011. Points of view: Color blindness. Nat Methods 8(6):441. doi:10.1038/nmeth.1618
• Zhu A, Ibrahim JG, Love MI. 2019. Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. Bioinformatics 35(12):2084-2092. doi:10.1093/bioinformatics/bty895

Related Skills

• differential-expression/de-results - Filter and rank DE result tables before plotting
• differential-expression/deseq2-basics - Run DESeq2 to produce the input results object
• differential-expression/de-visualization - DESeq2 / edgeR built-in plot helpers
• data-visualization/distribution-plots - Boxplot / raincloud follow-up for specific gene panels
• data-visualization/color-palettes - Okabe-Ito and CVD-safe palette selection
• data-visualization/ggplot2-fundamentals - Underlying grammar of graphics
• pathway-analysis/go-enrichment - Functional enrichment from the gene lists produced