bio-differential-expression-de-visualization

What this skill does

## Version Compatibility

Reference examples tested with: DESeq2 1.42+, edgeR 4.0+, limma 3.58+, ggplot2 3.5+, pheatmap 1.0+, RColorBrewer 1.1+, ggrepel 0.9+, EnhancedVolcano 1.20+, matrixStats 1.2+

Before using code patterns, verify installed versions match. If versions differ:
- 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.

# DE Visualization

**"Make the standard DE figure panel"** -> Use built-in functions or thin wrappers to produce diagnostic plots (dispersion, p-value histogram, PCA, sample distance) and result plots (MA, volcano, heatmap of top DE genes, per-gene counts), interpreted as diagnostics of the underlying model.

## Scope

This skill covers DE-specific built-in plots and immediate wrappers. For richer customization:
- Custom volcano/MA with apeglm-shrunken LFC and ggrepel labelling -> `data-visualization/volcano-and-ma-plots`
- PCA / UMAP / t-SNE customization -> `data-visualization/dimensionality-reduction-plots`
- Heatmap customization and ComplexHeatmap recipes -> `data-visualization/heatmaps-clustering`

## The Single Most Important Modern Insight -- A volcano with shrunken LFC compresses the cloud, but the p-values are unchanged

`lfcShrink()` pulls noisy estimates toward zero. On the volcano, that pulls genes horizontally toward the center. But the y-axis (`-log10(pvalue)`) is the **unshrunken Wald p-value** -- shrinkage does NOT recompute p-values (Zhu, Ibrahim, Love 2019 *Bioinformatics* 35:2084). A naive reader sees fewer extreme dots and concludes "fewer genes are significant". Wrong: the same genes are significant; the effect sizes are smaller and more honest.

Always label the volcano x-axis "shrunken log2 fold change (apeglm)" and note the y-axis comes from the unshrunken Wald test. The whole point of the apeglm volcano is the honest effect-size axis; if a publication shows an unshrunken volcano, it is showing inflated effects from low-count noise.

The MA plot has its own version of this: shrinkage flattens the left side (low-mean, formerly extreme LFC) and barely touches the right (high-mean, well-estimated LFC). That asymmetry is the visual signature of working shrinkage.

## Plot Taxonomy

| Plot | Diagnostic OR result | Built-in function | What it tests |
|------|---------------------|-------------------|---------------|
| Dispersion plot | Diagnostic | `plotDispEsts(dds)` (DESeq2), `plotBCV(y)` (edgeR) | Mean-dispersion trend fit quality |
| p-value histogram | Diagnostic | None; use ggplot2 | Null calibration, hidden batch, over-correction |
| PCA on VST/rlog | Diagnostic + result | `plotPCA(vsd, intgroup=...)` (DESeq2), `plotMDS()` (edgeR via limma) | Sample clustering, batch effects, outliers |
| Sample distance heatmap | Diagnostic | `pheatmap` on `dist(t(assay(vsd)))` | Within-group consistency, sample swaps |
| MA plot | Diagnostic + result | `plotMA(res)` (DESeq2), `plotMD(qlf)` (edgeR) | Normalization sanity, LFC vs mean |
| Volcano | Result | ggplot2 wrapper; `EnhancedVolcano` | Top-effect, top-significance gene story |
| Top-DE heatmap | Result | `pheatmap` on `assay(vsd)[sig_genes,]` | Per-gene pattern across conditions |
| `plotCounts` per gene | Result | `plotCounts(dds, gene, intgroup)` | Per-gene biology |

## Decision Tree by Scenario

| Scenario | Recommended approach |
|----------|---------------------|
| PCA for unbiased QC | `vst(dds, blind = TRUE)`; ask "do samples group as expected without design influence?" |
| PCA for results figure | `vst(dds, blind = FALSE)`; design is settled, accept its influence on dispersion |
| n < 30, library sizes vary >4x | `rlog(dds, blind = FALSE)` instead of vst |
| n > 30 | `vst()`; rlog impractical |
| Volcano | Plot shrunken LFC on x, unshrunken p-value on y; label both axes |
| Sample distance heatmap | `vst(blind = TRUE)`; tells if a sample is the wrong group regardless of design |
| Top-DE heatmap, want to see PATTERN | `scale = 'row'` (z-score per gene) |
| Top-DE heatmap, want to see ABSOLUTE LEVEL | `scale = 'none'` on `assay(vsd)`; otherwise weak signal looks strong |
| Top-variable-gene selection | `matrixStats::rowMads(assay(vsd))` instead of `rowVars` -- MAD is outlier-robust |
| n = 3, top genes in volcano | Note Schurch 2016 finding: 20-40% of true positives missed; treat as exploratory |
| Many groups, comparing DE sets | UpSet plot (Lex 2014); Venn drowns above 3 sets |

## Dispersion Diagnostic (Run This First)

**Goal:** Verify the dispersion-mean trend was fit acceptably before trusting any results.

**Approach:** `plotDispEsts(dds)` (DESeq2) or `plotBCV(y)` (edgeR) shows gene-wise (black/blue), fitted trend (red), and final shrunken (blue) dispersions vs mean.

```r
plotDispEsts(dds)

plotBCV(y)
```

| Pattern | Meaning | Action |
|---------|---------|--------|
| Cloud follows trend; final shrunken estimates pulled toward red curve | Healthy fit | Proceed |
| Red trend nowhere near the gene-wise cloud | Parametric trend failed | `DESeq(dds, fitType = 'local')` or `fitType = 'mean'` |
| Many gene-wise dispersions FAR ABOVE the trend | Outlier or unmodeled batch genes | Inspect rather than trust QL F-test alone |
| Final estimates much lower than gene-wise everywhere | Excessive shrinkage; sample too small or trend too flat | Check `useEM`, robust hyperparameter setting |
| BCV decreases monotonically with mean | Correct in edgeR | Default trend |

A plot inspected before trusting results is worth a hundred lines of statistical safeguards.

## P-value Histogram (Run This Second)

**Goal:** Detect model misspecification or hidden batch before reporting any gene list.

**Approach:** Histogram of raw p-values; under a correctly specified null, uniform with a spike near zero.

```r
library(ggplot2)
ggplot(res_df, aes(x = pvalue)) +
    geom_histogram(bins = 50, fill = 'steelblue', color = 'white') +
    labs(x = 'P-value', y = 'Frequency', title = 'P-value distribution') +
    theme_bw()
```

| Shape | Meaning | Action |
|-------|---------|--------|
| Uniform + spike at 0 | Correctly specified | Proceed |
| U-shape (spikes at 0 AND 1) | Anti-conservative; hidden batch or unmodeled covariate | Add the missing covariate; re-fit |
| Depleted near 0, spike near 1 | Conservative; over-modeled or wrong dispersion | Simplify model; check dispersion plot |
| Spike only at p = 1 | Discrete artifact from very-low-count genes | Pre-filter more aggressively |

## MA Plot (LFC vs Mean)

**Goal:** Inspect the relationship between LFC and mean expression for normalization correctness and shrinkage effect.

**Approach:** `plotMA` (DESeq2) or `plotMD` (edgeR). Always pick `ylim` deliberately; default can flatten the signal.

```r
plotMA(res, ylim = c(-5, 5), main = 'MA plot (unshrunken)')

res_apeglm <- lfcShrink(dds, coef = 'condition_treated_vs_control', type = 'apeglm')
plotMA(res_apeglm, ylim = c(-5, 5), main = 'MA plot (apeglm-shrunken)')

plotMD(qlf, main = 'edgeR MD plot')
abline(h = c(-1, 1), col = 'blue', lty = 2)
```

| Pattern | Meaning |
|---------|---------|
| Symmetric cloud centered at LFC = 0 | Correct normalization |
| Cloud median clearly above or below 0 | Normalization failed (TMM/RLE assumption violated) -- see normalization skill |
| Funnel widening at low mean | Expected (low counts noisier) |
| Dramatic up/down asymmetry | Possibly real (large biological perturbation), possibly normalization failure -- cross-check |
| Discrete horizontal bands at low mean | Low-count artifacts; pre-filter more aggressively |

The apeglm-shrunken MA visually flattens the left side; the post-shrinkage cloud should be tighter at low means.

## Volcano with Shrunken LFC

**Goal:** Show effect size vs significance with honest fold changes.

**Approach:** Use a built-in renderer (EnhancedVolcano for quick publication-quality output) on shrunken LFCs. Always plot shrunken LF