GPTomics/bio-data-visualization-heatmaps-clustering

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.

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.

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.

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$condition,
    Batch     = metadata$batch,
    Age       = anno_barplot(metadata$age),
    col = list(
        Condition = c(Control = '#56B4E9', Treatment = '#D55E00'),
        Batch     = c(A = '#009E73', B = '#0072B2', C = '#CC79A7')
    ),
    annotation_name_gp = gpar(fontsize = 8)
)

# Row metadata
ha_row <- rowAnnotation(
    Pathway = gene_info$pathway,
    LogFC   = anno_barplot(gene_info$log2FC, baseline = 0,
                            gp = gpar(fill = ifelse(gene_info$log2FC > 0,
                                                     '#D55E00', '#0072B2'))),
    col = list(Pathway = c(Metabolism = '#8491B4', Signaling = '#91D1C2'))
)

ht <- Heatmap(mat,
              name = 'Z-score',
              col  = col_fun,
              top_annotation  = ha_col,
              left_annotation = ha_row,
              row_split    = gene_info$pathway,
              column_split = metadata$condition,
              clustering_method_rows    = 'ward.D2',
              clustering_method_columns = 'ward.D2',
              clustering_distance_rows    = 'euclidean',
              clustering_distance_columns = 'euclidean',
              show_row_names = FALSE,
              use_raster = TRUE)          # rasterize cell layer for >2000 rows

draw(ht, merge_legends = TRUE)            # draw() not bare Heatmap()

The draw() requirement (silent failure)

A bare Heatmap(mat) works at the R console because auto-print invokes draw(). Inside for, lapply, function, Quarto/Rmd chunks, or Rscript, a bare Heatmap() produces no output and no error. Always wrap in draw() non-interactively. Only draw() exposes merge_legends, heatmap_legend_side, ht_gap, and padding.

seaborn.clustermap (Python)

import seaborn as sns
import numpy as np
import pandas as pd

# Robust symmetric bounds (1-99% quantile)
vmax = np.quantile(np.abs(df.values[~np.isnan(df.values)]), 0.99)

# col_colors / row_colors for categorical annotations
condition_colors = metadata['condition'].map({'Control': '#56B4E9', 'Treatment': '#D55E00'})
batch_colors = metadata['batch'].map({'A': '#009E73', 'B': '#0072B2', 'C': '#CC79A7'})
col_colors = pd.DataFrame({'Condition': condition_colors, 'Batch': batch_colors})

g = sns.clustermap(df,
                   cmap='RdBu_r', center=0, vmin=-vmax, vmax=vmax,
                   row_cluster=True, col_cluster=True,
                   method='ward',                      # seaborn uses scipy ward, equivalent to R ward.D2
                   metric='euclidean',
                   z_score=0,                          # 0 = rows, 1 = columns
                   col_colors=col_colors,
                   dendrogram_ratio=0.15,
                   cbar_pos=(0.02, 0.8, 0.03, 0.15),
                   figsize=(10, 12),
                   rasterized=True)                    # rasterize the cell layer

standard_scale vs z_score confusion:

• z_score=0 standardizes ROWS to mean 0, SD 1 (most common request)
• z_score=1 standardizes COLUMNS to mean 0, SD 1
• standard_scale=0 rescales ROWS to [0, 1] via (x − min) / (max − min) — NOT z-scoring, compresses outliers nonlinearly
• The two are mutually exclusive — passing both errors

A heatmap published with standard_scale looks like a z-scored heatmap but the color encoding is not interpretable as standard deviations.

OncoPrint -- The Specialized Mutation-Matrix Heatmap

OncoPrint (Cerami 2012 Cancer Discov 2:401; canonical at cBioPortal) is a stylized heatmap for mutation matrices where each cell encodes multiple alteration types via overlapping rectangles. Different from generic heatmaps — see data-visualization/oncoprint-mutation-matrices for the dedicated skill. Mentioned here only to note that ComplexHeatmap::oncoPrint() is the R implementation and inherits all the cluster/annotation machinery of Heatmap().

Per-Method Failure Modes

ward.D used when ward.D2 was intended

Trigger: clustering_method = 'ward' or 'ward.D' (with or without the .D).

Mechanism: R hclust ward.D is a legacy implementation that does NOT use squared distances — it does not implement Ward's minimum-variance criterion (Murtagh-Legendre 2014).

Symptom: Different dendrogram than published papers that say "Ward"; clusters look subtly different; reproducibility issues across R versions.

Fix: Always specify ward.D2 explicitly. For pheatmap and ComplexHeatmap, pass clustering_method_rows = 'ward.D2' and clustering_method_columns = 'ward.D2'.

One outlier compresses the color scale

Trigger: Plotting matrix without quantile clipping; one extreme value dominates the color range.

Mechanism: Default colorRamp2(c(min(mat), 0, max(mat)), ...) is dominated by the outlier — the rest of the matrix renders within a narrow band of pale colors.

Symptom: Heatmap looks "washed out" except for one cell or one column; biological pattern invisible.

Fix: bounds <- quantile(abs(mat), 0.99); col_fun <- colorRamp2(c(-bounds, 0, bounds), c('#0072B2', 'white', '#D55E00')). ComplexHeatmap's colorRamp2 does NOT clip values exceeding the range — they render at the extreme color, which is the intended behavior.

ComplexHeatmap silently produces no output in a script

Trigger: Bare Heatmap(mat) inside a for loop, lapply, function(), Quarto/Rmd code chunk, or Rscript invocation.

Mechanism: Auto-print only happens at the top-level R prompt; in non-interactive contexts the Heatmap object is created but never rendered.

Symptom: No error, no warning, no PDF output. Looks like the script ran successfully.

Fix: Always wrap in draw(): pdf('out.pdf', ...); draw(Heatmap(mat, ...)); dev.off(). Use the draw() call to set legend layout: draw(ht, merge_legends = TRUE, heatmap_legend_side = 'right').

Clustering applied to ordered conditions

Trigger: cluster_columns = TRUE when columns are an ordered sequence (time points, dose levels, treatment stages).

Mechanism: Hierarchical clustering re-orders columns to maximize within-cluster similarity, destroying the time/dose axis.

Symptom: Time-course heatmap with time points scrambled; reader cannot follow temporal pattern.

Fix: cluster_columns = FALSE for ordered conditions. To group while preserving order, use column_split or column_order explicitly.

Z-score on a sparse matrix

Trigger: Row z-scoring a matrix with many zero values (e.g., single-cell expression, sparse peak counts).

Mechanism: (x − mean) / sd is ill-defined when a row is mostly zeros — sd is dominated by the few non-zero values; z-scores explode for the non-zero entries.

Symptom: A few cells render as extreme colors; most cells are washed-out near-zero.

Fix: For single-cell, work with cluster-summarized pseudobulk matrices, not raw single-cell expression. For sparse peak data, filter rows by minimum non-zero count before scaling.

Correlation distance applied to data with batch effect

Trigger: clustering_distance_rows = 'correlation' on a matrix where samples have a strong batch effect.

Mechanism: Correlation preserves co-regulation pattern but is sensitive to global structure. If batch shifts ALL genes up in one batch, correlation distance reads this as "co-regulation."

Symptom: Modules cluster by batch, not by biology.

Fix: Batch-correct (limma::removeBatchEffect or ComBat) before clustering. Confirm with PCA that batch is no longer the dominant axis.

pheatmap's gaps_col interacts with cluster_cols

Trigger: Setting gaps_col = c(5, 10) with cluster_cols = TRUE.

Mechanism: When cluster_cols = TRUE, gaps_col is silently ignored; the dendrogram-determined order has no concept of position.

Symptom: No gaps appear; no warning.

Fix: To use gaps_col, set cluster_cols = FALSE and pre-arrange the columns explicitly. Or in ComplexHeatmap use column_split to combine clustering AND visual gaps.

Raster rendering at low DPI looks pixelated

Trigger: use_raster = TRUE (default for large heatmaps in ComplexHeatmap) with default raster_quality = 1.

Mechanism: Default raster quality is set for screen rendering; the bitmap is upscaled for PDF output, producing blocky cells.

Symptom: Heatmap cells look pixelated at print zoom; diagonal "stair-stepping" on cell boundaries.

Fix: Heatmap(..., use_raster = TRUE, raster_quality = 5) increases the raster resolution. Set raster_device = 'CairoPNG' for transparency support.

Reconciliation: When Methods Disagree

| Pattern | Likely cause | Action | |---------|--------------|--------| | ComplexHeatmap and pheatmap give different dendrograms | pheatmap uses dist() with default method = 'euclidean'; ComplexHeatmap defaults to the same but different clustering distance defaults | Verify both with clustering_distance_rows = 'euclidean', clustering_method_rows = 'ward.D2' explicitly | | Same code, different dendrogram across R versions | R 3.1 renamed 'ward' to 'ward.D' and added 'ward.D2' | Always specify ward.D2 explicitly; never 'ward' | | Z-scored heatmap with extreme colors only in a few cells | Sparse matrix with zero-inflation | Filter low-expression rows; OR shift to robust scaling ((x - median) / mad) | | Modules cluster by batch | Correlation distance picked up batch effect | Batch-correct upstream; verify via PCA | | seaborn clustermap produces different clusters than R | seaborn method='ward' calls scipy.cluster.hierarchy.linkage which IS ward.D2-equivalent; difference is usually metric default | Set metric='euclidean' explicitly in both | | OncoPrint mutual-exclusivity panel appears empty | column_order is being computed by clustering instead of preserved | Pass column_order = ... explicitly to oncoPrint |

Operational rule: a clustered heatmap is reproducible only when scaling, distance, linkage, color bounds, and (for OLO) the seriation method are all explicitly stated. Defaults differ across packages and across versions of the same package.

Quantitative Thresholds

| Threshold | Value | Source | |-----------|-------|--------| | Robust color bound | 1st-99th percentile of |matrix| | Standard publication practice; suppresses single-outlier dominance | | Raster trigger | >2000 rows or >2000 columns | ComplexHeatmap default use_raster = TRUE above 2000 | | OLO practical limit | ~5000 rows | O(n^4) worst case; modern Bar-Joseph implementations faster | | z-score symmetry | bounds around 0 | Z-scores are symmetric by construction | | Minimum non-zero count to z-score | >=3 non-zero values per row | Below this sd is unreliable | | Single-cell pseudobulk threshold | Downsample to <500 cells/group | Otherwise PDF rendering hangs |

Common Errors

| Error / symptom | Cause | Solution | |-----------------|-------|----------| | No PDF output from a script | Bare Heatmap() without draw() | Wrap in draw(); use pdf() / dev.off() explicitly | | Color scale washed out | One outlier dominates | Clip to 1st-99th percentile; symmetric bounds for diverging | | Time-course columns scrambled | cluster_columns = TRUE on ordered data | cluster_columns = FALSE; use column_split | | pheatmap gaps_col ignored | Conflict with cluster_cols = TRUE | Disable clustering OR switch to ComplexHeatmap split | | Dendrogram differs from a paper | Default clustering_method mismatch | Always specify ward.D2; never ward | | seaborn standard_scale interpreted as z-score | Different rescaling functions | Use z_score=0 for row z-scoring; standard_scale is min-max not z | | Heatmap renders pixelated | Default raster_quality = 1 | Set raster_quality = 5 for publication | | Z-score blows up for some rows | Sparse rows; near-zero sd | Filter low-expression rows OR robust scale | | Modules cluster by batch not biology | Batch effect not removed | limma::removeBatchEffect or ComBat upstream |

Anticipated Reviewer Pushback

| Pushback | Standard response | |----------|-------------------| | "Which clustering method?" | Explicit ward.D2 (Murtagh-Legendre 2014) on row-z-scored Euclidean distance. Alternatives evaluated in supplementary | | "How were leaves ordered?" | Optimal Leaf Ordering via seriation::seriate (Bar-Joseph 2001); reduces visually-adjacent dissimilarity | | "Why this color scale?" | Diverging palette symmetric around zero, bounds = 1st-99th percentile of |z|. Robust to outliers per standard practice | | "Why z-score?" | Removes absolute level so co-regulated genes cluster regardless of magnitude. Raw values clustered separately (supplementary) | | "Why row split by pathway?" | Pre-specified gene-set annotation (KEGG/Reactome) to verify clustering recovers known biology, not to bias it | | "Reproducibility across R versions?" | ward.D2 is stable across R 3.1+; clustering_method = 'ward' is not — never used |

References

• Bar-Joseph Z, Gifford DK, Jaakkola TS. 2001. Fast optimal leaf ordering for hierarchical clustering. Bioinformatics 17(suppl 1):S22-S29. doi:10.1093/bioinformatics/17.suppl_1.S22
• Cerami E, Gao J, Dogrusoz U, et al. 2012. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov 2(5):401-404.
• Crameri F, Shephard GE, Heron PJ. 2020. The misuse of colour in science communication. Nat Commun 11:5444.
• Gehlenborg N, Wong B. 2012. Points of view: Heat maps. Nat Methods 9(3):213.
• Gehlenborg N, Wong B. 2012. Points of view: Mapping quantitative data to color. Nat Methods 9(8):769.
• Gu Z, Eils R, Schlesner M. 2016. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 32(18):2847-2849.
• Murtagh F, Legendre P. 2014. Ward's hierarchical agglomerative clustering method: which algorithms implement Ward's criterion? J Classif 31(3):274-295.
• Nuñez JR, Anderton CR, Renslow RS. 2018. Optimizing colormaps with consideration for color vision deficiency to enable accurate interpretation of scientific data. PLOS ONE 13(7):e0199239.

Related Skills

• data-visualization/color-palettes - Sequential and diverging colormap selection
• data-visualization/oncoprint-mutation-matrices - Mutation-matrix heatmap (ComplexHeatmap oncoPrint)
• data-visualization/multipanel-figures - Combine heatmaps into journal layouts
• data-visualization/dimensionality-reduction-plots - PCA / UMAP as alternative views of the same matrix
• differential-expression/de-visualization - Heatmap of top DE genes
• single-cell/markers-annotation - Single-cell dotplot / matrixplot as scRNA alternatives