bio-data-visualization-heatmaps-clustering

What this skill does

## Version Compatibility

Reference examples tested with: ComplexHeatmap 2.18+, pheatmap 1.0.13 (still maintained as of 2025-06), circlize 0.4.16+, seaborn 0.13+, scipy 1.12+, scanpy 1.10+, ggplot2 3.5+.

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.

# Heatmaps and Hierarchical Clustering

**"Make a clustered heatmap"** -> Render an expression / feature matrix as a colored grid with hierarchical-clustering dendrograms, after committing to (a) how to scale the data (row z-score vs raw vs robust), (b) which distance metric (Euclidean vs correlation vs Manhattan), (c) which linkage criterion (ward.D2 vs complete vs average), (d) how to order the leaves (default vs optimal leaf ordering), (e) how to map values to color (sequential vs diverging, robust quantile bounds), and (f) which package can handle the matrix size and annotation complexity.

- R: `ComplexHeatmap::Heatmap()` (modern default), `pheatmap::pheatmap()` (still maintained, simpler API)
- Python: `seaborn.clustermap()`, `scanpy.pl.heatmap()` (single-cell-aware)

## The Single Most Important Modern Insight -- Distance, Linkage, and Scaling Are Three Independent Decisions

A heatmap's dendrograms are produced by three orthogonal choices, each with material biological consequences:

1. **Scaling** decides what "similar" means. Row z-scoring asks "do these genes covary across samples?" — it strips absolute level. Raw values ask "do these genes have similar magnitude AND pattern?" Robust scaling (quantile clip) asks "do these covary after suppressing outliers?"
2. **Distance metric** decides how dissimilar two profiles are. Euclidean on z-scored data ≈ `1 − Pearson` correlation; Manhattan tolerates outliers; correlation distance preserves co-regulation patterns regardless of amplitude.
3. **Linkage** decides how to merge clusters. Ward minimizes within-cluster sum of squares (compact, spherical clusters); complete uses max distance (compact, outlier-sensitive); average is balanced; single chains (almost never what genomics wants).

The biological story changes depending on these choices. A "module" identified with `complete` linkage on Euclidean distance of raw counts is *not the same module* identified with `ward.D2` on correlation distance of z-scored data. Both can be defensible; neither is automatic. **Pick deliberately, document, and verify the clustering against orthogonal evidence before claiming the modules are biological.**

## The ward.D vs ward.D2 Trap (Murtagh-Legendre 2014)

R `stats::hclust` exposes two methods both labeled "Ward": `ward.D` and `ward.D2`. They produce *different* dendrograms on the same data. Only `ward.D2` (squared distances input) implements Ward's actual minimum-variance criterion (Murtagh & Legendre 2014 *J Classif* 31:274). `ward.D` is a historical implementation that does not.

```r
hclust(dist(x), method = 'ward.D')   # NOT Ward's criterion -- legacy
hclust(dist(x), method = 'ward.D2')  # Ward's actual minimum-variance criterion
```

`pheatmap::pheatmap(clustering_method='ward.D2')` and `ComplexHeatmap::Heatmap(clustering_method_rows='ward.D2')` both pass through to `hclust`. Always specify `ward.D2` unless reproducing a paper that used the unlabeled `ward` (which actually called `ward.D` pre-R 3.1).

## Decision Tree by Scenario

| Scenario | Scaling | Distance | Linkage | Why |
|----------|---------|----------|---------|-----|
| Bulk RNA-seq, expression patterns across samples | row z-score | euclidean | ward.D2 | Standard; z-score removes absolute level so co-regulated genes cluster regardless of magnitude |
| Methylation beta values (already bounded [0,1]) | raw (no scale) | euclidean or manhattan | ward.D2 | Beta values are interpretable on absolute scale; scaling would distort |
| Co-expression module discovery | row z-score | correlation (`1 - cor`) | average | WGCNA convention; preserves co-regulation pattern |
| ChIP/ATAC peak intensity across samples | raw log-counts | euclidean | ward.D2 | Peaks are interpretable on absolute scale after log |
| Sample QC (correlation of samples) | column-wise raw | correlation | ward.D2 | The correlation IS the data; don't scale before computing it |
| Methylation array with outliers | raw + clip 1-99% | euclidean | ward.D2 | Outliers dominate Euclidean; robust clip preserves signal |
| Single-cell pseudobulk by cell type | row z-score | euclidean | ward.D2 | Same as bulk; downsample to <500 cells per type for rendering |
| Mutation matrix (binary present/absent) | raw | binary (jaccard) | average or complete | Standard distance for binary data; ward inappropriate |
| Drug response across cell lines | row z-score | spearman correlation | ward.D2 | Drug-rank patterns matter more than absolute IC50 |

## Color Mapping -- The Quietly Most-Important Choice

A heatmap is a color encoding of a matrix. The default linear mapping from data to color is rarely correct:

1. **Diverging data needs symmetric bounds.** For z-scores or log-fold changes, the color bar must be symmetric around zero. `colorRamp2(c(-2, 0, 2), c('#0072B2', 'white', '#D55E00'))` ALWAYS, not `colorRamp2(c(min, mean, max), ...)`.

2. **Robust quantile bounds.** Single outliers compress the entire color scale. Clip at 1st/99th percentile before mapping: `bounds <- quantile(mat, c(0.01, 0.99))`. ComplexHeatmap's `colorRamp2(c(bounds[1], 0, bounds[2]), ...)` is standard. Without this, one outlier sample turns the entire heatmap pale.

3. **Sequential data uses a perceptually-uniform colormap.** viridis, magma, cividis (Nuñez 2018), or batlow (Crameri 2020). NOT jet, NOT rainbow, NOT `colorRampPalette(c('blue','red'))(100)` which has a non-monotonic luminance.

4. **Diverging palettes from Crameri** (`vik`, `roma`) or ColorBrewer (`RdBu`, `BrBG`) are perceptually uniform. Reverse the default direction for log-fold-change (negative = blue, positive = red, by biological convention).

## Optimal Leaf Ordering (Bar-Joseph 2001)

A dendrogram for n leaves has 2^(n-1) consistent linear orderings — only one is the leaf order shown. Default `hclust` gives a deterministic but visually arbitrary ordering. **Optimal Leaf Ordering (OLO)** chooses the consistent ordering that minimizes the sum of distances between adjacent leaves — making visually adjacent rows actually similar, and revealing block structure in the heatmap that the default ordering hides.

```r
library(ComplexHeatmap)
library(seriation)

# OLO via seriation
dist_rows <- dist(mat)
hc_rows <- hclust(dist_rows, method = 'ward.D2')
olo_rows <- seriate(dist_rows, method = 'OLO', control = list(hclust = hc_rows))

Heatmap(mat,
        cluster_rows = as.dendrogram(olo_rows[[1]]),
        cluster_columns = TRUE,
        clustering_method_columns = 'ward.D2')
```

For matrices >2000 rows OLO becomes slow (O(n^4) in worst case; modern implementations are much faster). The trade-off is worth it for publication figures.

## Annotation Tracks -- ComplexHeatmap as the Reference

**Goal:** Render an annotated heatmap with column metadata (condition, batch, age), row metadata (pathway, gene class), and split panels for grouped display.

**Approach:** Define `HeatmapAnnotation` (column) and `rowAnnotation` objects with explicit color lists; render with `Heatmap()` specifying `row_split`/`column_split` for grouped layout; use `draw()` to commit, not bare `Heatmap()`, when running non-interactively.

```r
library(ComplexHeatmap)
library(circlize)

# Robust symmetric color mapping
bounds <- quantile(abs(mat[!is.na(mat)]), 0.99)
col_fun <- colorRamp2(c(-bounds, 0, bounds), c('#0072B2', 'white', '#D55E00'))

# Column metadata
ha_col <- HeatmapAnnotation(
    Condition = metadata$c