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

What this skill does

## 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 |

```r
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.

```r
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',
                TRUE ~ 'NS'
            ),
            neg_log10_p = -log10(pvalue)
        )

    if (is.null(label_genes)) {
        label_genes <- res %>%
            filter(significance != 'NS') %>%
            mutate(rank_score = -log10(pvalue) * abs(log2FoldChange)) %>%
            arrange(desc(rank_score)) %>%
            head(top_n) %>%
            pull(gene)
    }
    res$label <- ifelse(res$gene %in% label_genes, res$gene, '')

    okabe_ito <- c(Up = '#D55E00', Down = '#0072B2', NS = '#999999')

    ggplot(res, aes(log2FoldChange, neg_log10_p, color = significance)) +
        geom_point(alpha = 0.6, size = 1.3) +
        scale_color_manual(values = okabe_ito, name = NULL) +
        geom_vline(xintercept = c(-lfc_threshold, lfc_threshold),
                   linetype = 'dashed', color = 'grey40', linewidth = 0.3) +
        geom_hline(yintercept = -log10(fdr), linetype = 'dashed',
                   color = 'grey40', linewidth = 0.3) +
        geom_text_repel(aes(label = label), color = 'black', size = 3,
                        max.overlaps = Inf, box.padding = 0.4, segment.size = 0.2,
                        min.segment.length = 0) +
        labs(x = expression(log[2]~'fold change (shrunken)'),
             y = expression(-log[10]~italic(p))) +
        theme_classic(base_size = 10) +
        theme(panel.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

```r
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.