Bio Differential Expression De Visualization

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

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