Bio differential expression de visualization

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.

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.

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.

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 LFC; always set max.overlaps = Inf when labeling >10 genes -- the ggrepel default (10) silently drops labels. EnhancedVolcano accepts max.overlaps directly in 1.12+; version 1.10-1.11 has the older maxoverlapsConnectors argument (default 15); for either, falling back to options(ggrepel.max.overlaps = Inf) at the top of the script also works. For full ggplot2 customization (color schemes, faceting, label-set engineering), see data-visualization/volcano-and-ma-plots.

library(EnhancedVolcano)

res_apeglm <- lfcShrink(dds, coef = 'condition_treated_vs_control', type = 'apeglm')

EnhancedVolcano(res_apeglm,
    lab = rownames(res_apeglm),
    x = 'log2FoldChange', y = 'pvalue',
    pCutoff = 0.05, FCcutoff = 1,
    title = 'Treatment vs Control',
    subtitle = 'Shrunken LFC (apeglm); unshrunken Wald p',
    max.overlaps = Inf)

PCA on VST/rlog (Never on Raw Counts)

Goal: Show sample clustering by condition; detect batch effects, swaps, outliers.

Approach: Variance-stabilize first (VST or rlog), THEN PCA. Raw counts make PC1 = library size; log(counts+1) makes PC1 = mean expression. Neither carries biological signal until variance is stabilized.

vsd <- vst(dds, blind = FALSE)
plotPCA(vsd, intgroup = c('condition', 'batch'))

pca_df <- plotPCA(vsd, intgroup = c('condition', 'batch'), returnData = TRUE)
percentVar <- round(100 * attr(pca_df, 'percentVar'))

library(ggplot2)
ggplot(pca_df, aes(PC1, PC2, color = condition, shape = batch)) +
    geom_point(size = 4) +
    xlab(paste0('PC1: ', percentVar[1], '% variance')) +
    ylab(paste0('PC2: ', percentVar[2], '% variance')) +
    theme_bw()

library(limma)
plotMDS(cpm(y, log = TRUE), col = as.numeric(group), pch = 16)

blind=TRUE (default for vst()) re-estimates dispersions ignoring the design -- appropriate for unbiased QC ("are samples consistent independent of design?"). blind=FALSE uses the fitted dispersions -- appropriate for downstream visualization where the design is settled. Modern DESeq2 vignette recommends blind=FALSE for any plot after the model is fit.

PCA pattern	Interpretation	Action
Clear separation by condition on PC1 or PC2	Strong biological signal	Proceed
Separation by batch, not condition	Batch effect dominates	Include batch in design; DO NOT subtract before DE (see batch-correction Nygaard 2016)
One sample far from its group	Outlier or swap	Check library QC; sex check; somalier
Condition signal on PC3+, not PC1-PC2	Subtle effect	May still find DE; review dispersion plot
Two distinct sample clusters not explained by metadata	Hidden covariate	Investigate processing date, lane, machine

Sample Distance Heatmap (for QC)

library(pheatmap)
vsd <- vst(dds, blind = TRUE)
sd <- dist(t(assay(vsd)))
mat <- as.matrix(sd)
ann <- data.frame(condition = colData(dds)$condition,
                  row.names = colnames(dds))
pheatmap(mat, annotation_col = ann, annotation_row = ann,
         clustering_distance_rows = sd, clustering_distance_cols = sd,
         color = colorRampPalette(c('white', 'steelblue'))(100),
         main = 'Sample distance (vst blind)')

The diagonal should be dark; within-group samples should cluster. A within-group sample distant from its peers is a candidate for sample swap.

Top-DE Heatmap and the Row-Scaling Trap

Goal: Show expression patterns of significant genes across samples for results figure.

Approach: Use vst(blind=FALSE), select top genes (by adjusted p-value or MAD-robust variance), choose scaling deliberately.

library(pheatmap)

sig <- rownames(subset(res, padj < 0.01))[1:50]
vsd <- vst(dds, blind = FALSE)
mat <- assay(vsd)[sig, ]

mat_scaled <- t(scale(t(mat)))

ann_col <- data.frame(condition = colData(dds)$condition,
                      batch     = colData(dds)$batch,
                      row.names = colnames(mat))

pheatmap(mat_scaled, annotation_col = ann_col,
         show_rownames = FALSE,
         clustering_distance_rows = 'correlation',
         clustering_distance_cols = 'correlation',
         color = colorRampPalette(c('blue', 'white', 'red'))(100),
         main = 'Top 50 DE genes (z-scored per gene)')

scale='row' (z-score per gene) is the conventional choice for "show me patterns". It DESTROYS absolute expression level information -- a gene at 5-7 with mean 6 looks identical to a gene at 10-1000. For pattern detection: correct. For QC heatmaps showing batch shifts: WRONG -- use scale='none' on assay(vsd).

Top-variable-gene selection robustness:

library(matrixStats)
vars_mad <- rowMads(assay(vsd))
top500 <- order(vars_mad, decreasing = TRUE)[1:500]

rowMads (median absolute deviation) is outlier-robust; rowVars is dominated by single-outlier-sample genes. For exploratory PCA of "top variable genes", MAD selection avoids artifacts.

Per-gene Plot

plotCounts(dds, gene = 'GENE_NAME', intgroup = 'condition')

d <- plotCounts(dds, gene = 'GENE_NAME', intgroup = c('condition','batch'),
                returnData = TRUE)
library(ggplot2)
ggplot(d, aes(x = condition, y = count, color = batch)) +
    geom_jitter(width = 0.1, size = 3) +
    scale_y_log10() +
    ggtitle('GENE_NAME') +
    theme_bw()

With n=3, the boxplot is misleading (3 points per box). Prefer geom_jitter over geom_boxplot at small n.

UpSet for Multi-set Comparisons

For >3 DE gene sets (e.g., contrasts treated_drugA, treated_drugB, treated_drugC each vs control), Venn diagrams become unreadable. UpSet (Lex et al. 2014 IEEE Trans Vis Comput Graph 20:1983) scales:

library(UpSetR)
upset(fromList(list(drugA = sig_drugA, drugB = sig_drugB, drugC = sig_drugC)))

Per-Method Failure Modes

Volcano with unshrunken LFC -- inflated story

Trigger: ggplot(res_df, aes(x=log2FoldChange, ...)) without lfcShrink(); extreme dots at the corners are low-count genes.

Mechanism: Unshrunken MLE LFCs are dominated by very-low-count genes whose log ratios are noisy. The visual top-left and top-right corners look impressive but are artifacts.

Symptom: Top genes by abs(LFC) are obscure low-count genes; reviewer asks "why are these the top hits?"

Fix: res_apeglm <- lfcShrink(dds, coef=..., type='apeglm'); plot from res_apeglm. Label axis "shrunken log2 fold change (apeglm)".

ggrepel max.overlaps silently drops labels

Trigger: geom_text_repel(data = top30, aes(label = gene)); only 10 labels render.

Mechanism: Default max.overlaps = 10; warning printed but easily missed in a knitr/Quarto render.

Symptom: Reviewer asks "where is gene X?"; it was in top30 but did not render.

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

PCA shows batch, not condition

Trigger: plotPCA(vsd, intgroup='batch') cleanly separates batches; intgroup='condition' does not separate.

Mechanism: Batch variance exceeds condition variance.

Symptom: Treatment effect looks weak; DE p-values inflated if batch not in design.

Fix: Include batch in design (design = ~ batch + condition). DO NOT use removeBatchEffect then re-do DE on corrected counts (Nygaard 2016 cardinal sin -- see batch-correction). For VISUALIZATION only, removeBatchEffect is OK.

Heatmap row-scaling hid a sample-level shift

Trigger: QC heatmap with scale='row' looks consistent within group; downstream PCA shows clear sample outlier.

Mechanism: z-score per gene removes per-sample additive shifts. A sample that's globally inflated 1.5x looks identical to peers after row scaling.

Symptom: "The heatmap looked fine but PCA shows a problem."

Fix: For QC heatmaps, use scale = 'none' on assay(vsd) directly. For result heatmaps after QC is clean, scale = 'row' is the appropriate choice for pattern emphasis.

Top-N-by-rowVars dominated by single-outlier-sample genes

Trigger: "Top 500 variable genes" PCA shows a striped pattern, one or two samples driving the spread.

Mechanism: rowVars is squared-deviation; one outlier sample of one gene inflates that gene's "variance" massively.

Symptom: Top variable gene list includes many genes where N-1 samples are flat and one sample is extreme.

Fix: matrixStats::rowMads() for MAD-based selection; or genefilter::rowQ().

Common errors

Error / symptom	Cause	Fix
plotPCA reports only 2 PCs	DESeq2 plotPCA is hard-coded to PC1/PC2	Use prcomp(t(assay(vsd))) and plot any pair
PCA cloud collapses to one point	Forgot to log-transform; raw counts plotted	vst(dds) first
All MA-plot points red	alpha set too high or sig-flag bug	Verify alpha; check padj vs pvalue in flag
pheatmap complains "infinite values"	NA / Inf in scaled matrix; gene with zero variance	Remove zero-variance rows before scaling
Volcano axis labels obscured	Default ggplot theme too compact	theme_bw(base_size = 14)
plotCounts says gene not found	Wrong ID type (symbol vs Ensembl)	Match rownames(dds) exactly
vst() errors with very low gene count post-filter	Default nsub=1000 exceeds available genes	Lower nsub (e.g., vst(dds, nsub=500))

References

• Anders S, Huber W. 2010. Differential expression analysis for sequence count data. Genome Biol 11(10):R106. doi:10.1186/gb-2010-11-10-r106
• Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15(12):550. doi:10.1186/s13059-014-0550-8
• 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
• Robinson MD, McCarthy DJ, Smyth GK. 2010. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26(1):139-140. doi:10.1093/bioinformatics/btp616
• Lex A, Gehlenborg N, Strobelt H, Vuillemot R, Pfister H. 2014. UpSet: Visualization of Intersecting Sets. IEEE Trans Vis Comput Graph 20(12):1983-1992. doi:10.1109/TVCG.2014.2346248
• Schurch NJ et al. 2016. How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA 22(6):839-851. doi:10.1261/rna.053959.115
• Nygaard V, Rødland EA, Hovig E. 2016. Methods that remove batch effects while retaining group differences may lead to exaggerated confidence in downstream analyses. Biostatistics 17(1):29-39. doi:10.1093/biostatistics/kxv027

Related Skills

• deseq2-basics - Generates the dds / res objects plotted here; vst/rlog choice
• edger-basics - Generates y / qlf for plotMD, plotBCV, plotMDS
• de-results - p-value histogram, padj=NA diagnosis informs what to plot
• batch-correction - removeBatchEffect for visualization only (never as DE input)
• expression-matrix/normalization - VST vs rlog vs log-CPM mechanics
• data-visualization/volcano-and-ma-plots - Full custom volcano/MA with apeglm + ggrepel
• data-visualization/dimensionality-reduction-plots - PCA, UMAP, t-SNE customization
• data-visualization/heatmaps-clustering - pheatmap and ComplexHeatmap recipes
• data-visualization/upset-plots - UpSet plot customization