Heatmaps for Gene Expression: Clustering, Color Scales, and Annotation in R

The default heatmap from heatmap() in base R is the figure that gets you a reviewer comment about “publication-quality figures.” The colors are red-green, the dendrogram is on top of the rows, the labels overflow the margins, and the scaling is whatever R guessed. None of those defaults are right, and each of them is a separate choice you have to make deliberately.

This guide covers the four choices that determine whether a gene expression heatmap reads as signal or as static: scaling (rows, columns, or none), clustering method, color palette, and row/column ordering. Code uses pheatmap and ComplexHeatmap in R; the principles transfer to Python's seaborn or Morpheus.

When a Heatmap Is the Right Plot

Heatmaps work when you have:

• A matrix — rows are genes (or samples), columns are samples (or conditions), and the cell value is an expression measurement.
• 20–200 features to show. Fewer than ~20 and a bar plot or line plot is more readable. More than ~500 and individual rows become a 1-pixel smear — switch to a different visualization (PCA, UMAP, or pathway-level summary).
• A clustering or grouping you want to surface. If the rows are pre-ordered (e.g., by fold change) and you're just showing values, a heatmap is showy but adds little over a sorted table.

If you have an arbitrary RNA-seq results table with 20,000 genes, a heatmap is the wrong starting point — build a volcano plot first to pick the candidates, then heatmap the 50–100 you want to show.

Choice 1: Scaling (the choice that matters most)

Raw expression values span orders of magnitude. A highly expressed housekeeping gene at TPM = 5,000 and a low-expressed transcription factor at TPM = 5 cannot be color-encoded on the same scale — the TF row will be uniformly black and you'll see nothing. Scaling fixes this, but the scaling choice changes what the figure shows:

• Row scaling (z-score per gene): subtract the row mean, divide by the row SD. Each gene's color now shows how each sample compares to that gene's mean. This is the right default for “which samples upregulate which genes.” The trade-off: absolute expression level is hidden — a uniformly low-expressed gene and a uniformly high one can look identical.
• Column scaling: same z-score, but per sample. Rarely useful for gene expression — usually applied when columns are heterogeneous metabolites or proteins measured at different scales.
• No scaling: use only when values are already comparable — log₂ fold changes (already a relative measure), already-normalized expression, or correlation matrices. Plotting raw counts with no scaling produces an uninterpretable figure.
• Log-transform before scaling: for raw counts or TPM/FPKM, apply log2(x + 1) first, then row-scale. Gene expression is approximately log-normal because expression results from multiplicative processes; the log puts values on a scale where SD is biologically meaningful.

Choice 2: Clustering Method

Hierarchical clustering on rows and/or columns groups similar profiles together. Three parameters — distance metric, linkage method, and what to cluster — each have defaults you should change deliberately:

Distance metric measures how different two profiles are:

• Euclidean distance — default. Works when values are scaled (z-scores). Bad for raw expression because it's dominated by the highest-expressed gene.
• 1 − Pearson correlation — better for gene expression. Treats two genes that go up and down together as similar regardless of overall expression level. Use when scaling won't be applied.
• 1 − Spearman correlation — robust to outlier samples. Use when one or two samples have unusually high/low expression that you don't want to drive the clustering.

Linkage method determines how clusters merge:

• complete (the pheatmap default) — merges based on the farthest pair. Produces compact, tight clusters. Sensitive to outliers.
• average (UPGMA) — merges based on the mean pairwise distance. The standard in microarray analysis. A reasonable default for gene expression.
• Ward.D2 — minimizes within-cluster variance. Produces balanced, even-sized clusters. Often the best choice when you want to read distinct cluster blocks from the heatmap.
• single — chains, produces long stringy clusters. Almost never the right choice for gene expression.

For most RNA-seq heatmaps, the combination that produces readable figures is row-scaling + Pearson distance + Ward.D2 linkage. The pheatmap defaults (no scaling, Euclidean distance, complete linkage) are not those.

Choice 3: Color Palette

The red-green default fails for ~8% of male readers (the most common form of color blindness). It also visually overweights the extremes, making moderate effects look stronger than they are. Three palettes that work for gene expression:

• Blue-white-red — for row-scaled data centered on zero. The diverging palette signals “negative below zero, positive above.” ColorBrewer's RdBu reversed is the standard implementation.
• Viridis — perceptually uniform, colorblind-safe, monotonic. The right choice for non-scaled data (e.g., log2 expression values) where a single low-to-high direction is the message.
• Magma / Plasma — variants of viridis with different hue ranges. Magma works well for dark-background presentations; plasma for print.

The ColorBrewer reference is the canonical source for diverging and sequential palettes. Avoid rainbow palettes (they're not perceptually uniform) and any palette where the brightest color is yellow (yellow recedes against white backgrounds — the very feature you want to highlight disappears).

Choice 4: Row and Column Ordering

Default ordering is alphabetical (gene symbols A–Z), which is essentially random. Three principled alternatives:

1. Hierarchical clustering (the heatmap default). Genes and samples are reordered so similar profiles sit next to each other. The dendrogram shows the relationship structure.
2. Manual / biological ordering. When the rows are pathway-curated (e.g., “the 15 cell-cycle genes” in a known biological order), preserve that order. Don't cluster — the order carries information.
3. By value. Sort rows by mean expression, by fold change, or by a sample-specific value. Useful for “top genes by criterion X” figures where the reader wants to read down the list.

For samples (columns), almost always cluster — the dendrogram on top reveals replicate structure and batch effects. If two biological replicates of the same condition cluster apart with replicates of a different condition between them, you have a batch effect to address before publishing the figure.

Building the Heatmap in R

pheatmap is the most readable R package for ~80% of gene expression heatmaps. A working example with the choices above:

library(pheatmap)

# Assuming `counts` is a gene-by-sample matrix and `sample_info` is a data.frame
# with sample annotations (treatment, timepoint, etc.)

log_counts <- log2(counts + 1)
top_var_genes <- order(apply(log_counts, 1, var), decreasing = TRUE)[1:100]
mat <- log_counts[top_var_genes, ]

pheatmap(
  mat,
  scale = "row",
  clustering_distance_rows = "correlation",
  clustering_distance_cols = "correlation",
  clustering_method = "ward.D2",
  color = colorRampPalette(rev(RColorBrewer::brewer.pal(11, "RdBu")))(50),
  annotation_col = sample_info[, c("treatment", "timepoint"), drop = FALSE],
  show_rownames = TRUE,
  show_colnames = TRUE,
  fontsize_row = 6,
  fontsize_col = 8,
  border_color = NA,
  cellwidth = 14,
  cellheight = 7
)

The annotation_col argument adds colored bars above the columns identifying each sample's treatment and timepoint — essential for reading clustering against the experimental design. border_color = NA removes gridlines (the default white grid is visually noisy at high cell counts). cellwidth and cellheight force a printable aspect ratio rather than letting the device size dictate it.

When to Reach for ComplexHeatmap

pheatmap stops being enough when you need:

• Multiple side annotations (e.g., gene category bars on the left AND fold change bars on the right)
• Split heatmaps (rows or columns grouped into separate blocks with gaps)
• Heatmaps combined with bar plots, point plots, or density plots in the same figure panel
• Programmatic generation across many comparisons with consistent formatting

The ComplexHeatmap Bioconductor package (Gu, 2022) is the standard for these. It's harder to use than pheatmap for a basic plot but unmatched once the figure has multiple annotation tracks. The package is now the most-cited heatmap tool in genomics papers.

Spot-Check Your Heatmap Before Submission

• Color bar legend has numeric ticks, not just “Low/High.”
• Scale is labeled — “Row z-score,” “log₂ TPM,” “log₂ fold change.” Don't make the reader infer.
• Replicates cluster with replicates. If they don't, the figure exposes a batch effect — address it before publishing.
• Row labels are legible at print size. 6-point font for 100+ genes is usable; 3-point isn't. If you can't read it at print size, reduce row count.
• Color is colorblind-safe. Test with a simulator (Coblis, Sim Daltonism) before submission.
• Clustering method is stated in the legend. “Hierarchical clustering with Pearson distance and Ward.D2 linkage” tells the reader what they're seeing.

Methods sentence: “The 100 most variable genes after log₂(x + 1) transformation were row-scaled and clustered using Pearson correlation distance with Ward.D2 linkage, visualized with pheatmap (v1.0.12) in R 4.3.1.” Specific, reproducible, defensible — the same bar reviewers apply to the rest of the methods section.