Clustering and Heatmaps

Clustering is the assignment of objects — here, samples (e.g., patients, cell types, treatment conditions) or features (e.g., genes, proteins, metabolites) — into groups such that objects within a group are more similar to each other than to objects in other groups. It is unsupervised: no label ("tumor" vs. "normal", "responder" vs. "non-responder") is given to the algorithm in advance. The groups it finds are candidate structure in the data, to be interpreted, or discarded, afterward. Two ingredients make clustering possible. First, a distance or dissimilarity metric: a rule that converts a pair of expression profiles (each a vector of thousands of numbers) into a single number describing how far apart they are. Second, a clustering algorithm that uses those pairwise distances to build groups. The two dominant algorithm families are hierarchical clustering, which builds a nested tree of clusters at every scale simultaneously, and partitioning methods such as k-means, which split the data into a fixed number k of non-overlapping groups directly. Hierarchical clustering is the method almost always paired with heatmaps, because its output — a tree, or dendrogram — gives a natural ordering for rows and columns that groups similar things adjacently, exactly what a heatmap needs to be readable. A heatmap is then simply the data matrix rendered as a grid of colored cells, with rows and columns reordered according to their dendrograms and the two trees drawn along the matrix edges. The color encodes a value — usually not raw expression, but a standardized value (a z-score: how many standard deviations above or below that gene's own mean, across samples, this particular reading falls). Standardizing row-wise is what allows genes with wildly different absolute expression levels (a housekeeping gene expressed at 10,000 units, a transcription factor expressed at 20) to sit on the same red-blue color scale and be compared meaningfully by shape of variation rather than magnitude.

Step 1 — Choose and compute a distance metric. Every clustering decision begins here, and the choice changes the answer. For two expression profiles $x$ and $y$ across $n$ genes:

Euclidean distance treats each gene as a coordinate axis and measures straight-line distance:

This is sensitive to absolute scale and magnitude — two samples with the same pattern of ups and downs but different overall amplitude will be scored as far apart.

1 − Pearson correlation measures whether two profiles rise and fall together, ignoring absolute magnitude:

Correlation distance is the standard choice for gene-expression heatmaps precisely because biologists usually care about co-regulation (shape of the pattern), not raw signal intensity, which is confounded by technical factors like RNA amount and probe efficiency.

Manhattan distance ($\sum_i |x_i - y_i|$) is more robust to outlier genes than Euclidean because it does not square the differences.

Step 2 — Build the tree with agglomerative hierarchical clustering. The dominant algorithm is agglomerative (bottom-up): 1. Start with every sample as its own cluster of size one. 2. Compute all pairwise distances. 3. Merge the two closest clusters into one. 4. Recompute the distance from this new cluster to every remaining cluster, using a linkage rule. 5. Repeat steps 3–4 until only one cluster (containing everything) remains. The linkage rule is the crucial, under-appreciated choice governing how a merged cluster's distance to others is defined:

The same distance matrix, run through single versus complete linkage, can produce visibly different trees — this is not a minor implementation detail, it is a modeling choice with real consequences for what "clusters" appear.

Step 3 — Cut or read the tree. The dendrogram encodes a nested hierarchy: at low height, only near-identical profiles are merged; at high height, everything eventually merges into one root. The height (or branch length) at which two branches join represents the dissimilarity at which those two groups were judged "the same." A discrete cluster assignment is obtained by cutting the tree at some chosen height — but the tree itself does not tell you where to cut. Choosing 2 versus 4 versus 8 clusters is an analyst decision, often guided by heuristics like the Statistics of High-Dimensional Data concept of within-cluster variance explained, or by simply looking for visually long branches ("the biggest jumps").

Step 4 — Reorder the matrix and color it. The heatmap draws the expression matrix with rows and columns permuted to match the leaf order of their respective dendrograms — genes that clustered together sit in adjacent rows, samples that clustered together sit in adjacent columns. Before coloring, each gene's row is typically standardized (subtract its mean, divide by its standard deviation, across samples) so that color encodes relative pattern (how far above or below this gene's own typical level) rather than absolute abundance. A diverging color scale — traditionally red (high) to black or white to blue (low) — is centered at zero so that "no change" is visually neutral.

Common alternative: k-means and other partitioning. k-means clustering, by contrast, requires the analyst to fix the number of clusters k in advance, then iteratively (1) assigns each point to its nearest of k cluster centers, and (2) recomputes each center as the mean of its assigned points, repeating until assignments stop changing. It scales better to very large numbers of samples (single-cell datasets with hundreds of thousands of cells routinely use k-means or graph-based variants like Louvain clustering rather than full hierarchical clustering, which requires computing and storing all pairwise distances — an $O(n^2)$ cost that becomes prohibitive past tens of thousands of objects). But k-means does not naturally produce the nested tree a heatmap's dendrogram needs, which is why hierarchical clustering remains the default for heatmap figures even when other methods are used for the primary analysis.

The result that made hierarchical clustering and heatmaps central to biology — rather than a niche statistical technique — was the discovery of intrinsic molecular subtypes of breast cancer by Charles Perou, Therese Sørlie, and colleagues, published in Nature (2000) and PNAS (2001).

What was done. The team used cDNA microarrays to measure the expression of roughly 8,000 genes across 65 surgical breast tumor samples (plus paired samples before and after chemotherapy from some patients, and normal tissue controls). Critically, they did not begin with a hypothesis that "tumors fall into four subtypes." They computed pairwise correlation distances between all samples across the full gene set, and applied average-linkage hierarchical clustering to both genes and samples, then displayed the reordered matrix as a heatmap — precisely the pipeline described above. They also identified an "intrinsic gene set" — genes whose expression varied more between different patients' tumors than between repeated samples from the same tumor — to focus the clustering on genes likely to reflect stable tumor biology rather than sampling noise.

What was observed. The sample dendrogram did not produce a random-looking tree. Instead, it split cleanly into a small number of major branches, each corresponding to a block of co-varying genes visible as a colored rectangle in the heatmap. These branches lined up with genes known from independent biology: one cluster of tumors expressed estrogen-receptor pathway genes strongly (later named luminal), one cluster overexpressed ERBB2 (HER2-enriched), and one cluster expressed genes otherwise associated with basal/myoepithelial breast cells and lacked estrogen-receptor and ERBB2 expression (basal-like). Repeated samples from the same tumor, taken at different times or from different biopsy sites, reliably clustered next to each other rather than scattering randomly across the tree — direct evidence that the clustering was picking up a stable, biological signal rather than technical noise.

Why this counted as proof, not just a pretty picture. Two follow-up pieces of evidence elevated this from "an interesting pattern in a heatmap" to an accepted biological classification. First, Sørlie et al. (2001, 2003) showed the same subtype structure re-emerged in independent patient cohorts profiled on different arrays — the clusters were reproducible, not an artifact of one dataset. Second, and decisively, the subtypes had different clinical outcomes: patients with basal-like tumors had significantly worse relapse-free and overall survival than luminal-A patients, a finding replicated across cohorts and later confirmed with immunohistochemistry (ER, PR, HER2, and basal cytokeratin staining) as a cheaper clinical surrogate for full expression profiling. A heatmap pattern that also predicts who lives and who dies is no longer merely descriptive — it is a validated biological and clinical entity, now embedded in the WHO classification of breast cancer and in treatment guidelines (e.g., basal-like/triple-negative tumors do not benefit from anti-HER2 or anti-estrogen therapy and are treated differently as a direct consequence of this clustering result).

The negative control that matters equally: Luria–Delbrück-style caution isn't the model here, but a parallel warning is. Clustering algorithms will always output clusters, even applied to pure random noise with no true group structure — this is a mathematical certainty of the algorithm, not evidence of biology. The reason breast-cancer subtyping is trusted is not that clustering "worked" once, but that the specific clusters found were (a) reproducible across independent cohorts, (b) aligned with independently known biology (ER, HER2 status), and (c) prognostically meaningful. Any heatmap lacking these three checks should be read as a hypothesis, not a finding.

Example 1 — computing distances by hand and seeing why the metric matters. Four samples, expression of one gene measured across four conditions (values in arbitrary fluorescence units):

Compare S1 and S4 by Euclidean distance:

Compare S1 and S2:

By Euclidean distance, S1 and S2 look far closer than S1 and S4 — dominated by the fact that Gene B has much larger absolute values than Gene A, so differences in Gene B dominate the sum of squares. This is exactly the scale-sensitivity problem: Gene B, simply by having bigger numbers, drowns out Gene A's contribution entirely, even though Gene A might be the biologically interesting one. Now standardize each gene (z-score across the four samples) before computing distance. Gene A has mean 6.75, SD ≈ 8.6; Gene B has mean 21, SD ≈ 13.6. After standardizing, both genes contribute comparably to distance, and the resulting hierarchical clustering — say, using average linkage — will group samples by relative pattern across both genes rather than being dictated by whichever gene happens to have larger raw units. This single worked calculation is why the standard pipeline standardizes rows before clustering and coloring: skipping this step is one of the most common ways a heatmap misleads.

Example 2 — reading a dendrogram cut and clinical stakes. Suppose a hierarchical clustering of 12 tumor samples (average linkage, 1 − Pearson correlation distance) yields a dendrogram whose top-level split (highest merge height, dissimilarity = 0.95) separates 5 samples from 7. Within the 7-sample branch, a second, much lower split (dissimilarity = 0.35) further divides it into groups of 4 and 3. An analyst must decide: report 2 clusters (5 vs. 7), or 3 clusters (5 vs. 4 vs. 3)? The height gap between the top split (0.95) and the next split (0.35) is large — this "big jump, then small jumps" pattern is one common heuristic for choosing where to cut: it suggests that the 5-vs-7 split reflects a real, strong difference (worth calling two clusters), while the split within the 7-sample group is comparatively weak and might simply reflect finer, noisier substructure rather than a distinct biological subtype. Say gene-level annotation of the two major branches shows the 5-sample branch is uniformly high for a basal marker (KRT5) and low for the estrogen receptor gene (ESR1), while the 7-sample branch shows the reverse. This external, independent biological annotation is what turns "two branches in a tree" into a defensible claim of "basal-like vs. luminal tumors" — exactly the reasoning that validated the Perou/Sørlie subtypes. Absent such an independent anchor, a clinician should treat the 5-vs-7 split as a candidate grouping requiring confirmation (e.g., immunohistochemistry, an independent cohort), not as an established diagnosis — the dendrogram alone never proves clinical relevance.

Charles Perou trained as a molecular biologist and, working with the laboratory of Patrick Brown and David Botstein at Stanford in the late 1990s — the same group that helped develop the cDNA microarray itself — turned the newly invented technology on human breast tumors rather than model organisms. This was a deliberately clinical application of a tool initially built for basic yeast and cell-biology research. Perou's key methodological contribution was the "intrinsic gene set": recognizing that to find biologically stable tumor subtypes, you should first filter for genes whose expression varies more between patients than within repeated samples of the same patient's tumor — a subtle but crucial statistical move that reduced noise before clustering. Working with Therese Sørlie (then a visiting researcher from Norway) and the wider Stanford/Norway collaboration, Perou co-authored the 2000 Nature paper "Molecular portraits of human breast tumours" and the 2001 PNAS paper establishing that the resulting subtypes predicted survival. He has continued to lead refinement of the intrinsic subtype classification (the PAM50 gene panel used clinically today derives directly from this work) and remains active in breast-cancer genomics at UNC Chapel Hill.

Clinical oncology. The PAM50 test, a 50-gene expression assay directly descended from the intrinsic-subtype clustering work, is used in hospitals to classify breast tumors and guide chemotherapy versus hormone-therapy decisions — a heatmap-discovered structure is now a billed clinical diagnostic.

Single-cell genomics. Modern single-cell RNA-sequencing atlases cluster tens of thousands to millions of cells to discover novel cell types and states (e.g., the Human Cell Atlas project); heatmaps of marker genes across identified clusters remain the standard way researchers communicate and validate that a cluster corresponds to a known or novel cell type.

Drug discovery and repurposing. Connectivity Map–style approaches cluster the gene-expression signatures of thousands of drug treatments, grouping compounds whose transcriptional effects resemble each other — enabling repurposing hypotheses (a compound clustering near a known anti-inflammatory may share its mechanism).

Microbiome studies. Clustering of microbial community composition across samples (e.g., gut samples from healthy versus diseased cohorts) using Bray-Curtis or UniFrac distances — close relatives of the distance metrics described here — routinely appears as clustered heatmaps in microbiome and public-health research on conditions from inflammatory bowel disease to malnutrition.

Agriculture and plant breeding. Expression or genotype clustering of crop germplasm collections helps breeders identify genetically distinct subpopulations for trait mapping and to avoid redundant crosses between near-identical lines.

Conservation genetics. Hierarchical clustering of genetic distance matrices between populations informs conservation units and management boundaries — the same tree-building math, applied to allele frequencies instead of gene expression.

Ask the tutor: If I re-run hierarchical clustering on the same expression data but change only the linkage method, how different can the resulting clusters legitimately be, and how would I tell whether my biological conclusion is robust to that choice?