Applied - Heatmaps and Clustermaps

In one line: the expression heatmap — where the colormap choice, the scaling choice and the clustering choice each independently determine what the reader concludes.

Uses: Seaborn · Colormaps and Color · Figure-level vs Axes-level · Linear Algebra with NumPy · Reshaping with pivot and melt

Basic heatmap

import seaborn as sns, matplotlib.pyplot as plt, numpy as np, pandas as pd
sns.set_theme(style="white", context="paper")

fig, ax = plt.subplots(figsize=(6, 8))
sns.heatmap(mat, cmap="rocket", ax=ax,
            cbar_kws=dict(label=r"log$_2$ CPM"),
            xticklabels=True, yticklabels=True)

sns.heatmap does not cluster. It draws the matrix in the order you give it. sns.clustermap clusters.

Clustermap — the workhorse

cg = sns.clustermap(
    mat,                             # genes × samples
    z_score=0,                       # 0 = z-score ROWS (genes); 1 = columns
    cmap="vlag", center=0,
    figsize=(8, 10),
    method="average", metric="correlation",
    col_colors=col_colors,           # annotation bars
    row_cluster=True, col_cluster=True,
    dendrogram_ratio=(0.1, 0.15),
    cbar_pos=(0.02, 0.85, 0.03, 0.1),
    xticklabels=True, yticklabels=False,
)
cg.ax_heatmap.set(xlabel="", ylabel="")
cg.savefig("clustermap.pdf", dpi=300)

Handles to reach through:

cg.ax_heatmap; cg.ax_row_dendrogram; cg.ax_col_dendrogram; cg.ax_cbar
cg.dendrogram_row.reordered_ind        # ← the clustered gene order
cg.dendrogram_col.reordered_ind
cg.data2d                               # the transformed matrix as plotted

reordered_ind is genuinely useful — export it to reuse the same gene order in another figure, or to write out the clustered table.

The three decisions

1. Scaling

z_score=0        # each ROW centred and scaled → compare patterns ACROSS genes
z_score=1        # each COLUMN scaled → compare across samples
standard_scale=0 # scale rows to [0, 1]
None             # raw values

Without scaling, a heatmap of raw expression shows you which genes are highly expressed — and nothing else. Ribosomal genes are bright, everything else is dark, and every sample looks the same. That is almost never the question.

With z_score=0, you see relative patterns per gene, which is usually the question ("which samples are high for this gene?"). But it also amplifies noise in genes with tiny absolute variation — a gene ranging 0.1–0.15 CPM z-scores into a dramatic-looking pattern.

The defensible combination: filter to genes with meaningful variation first, then z-score.

top_var = mat.var(axis=1).nlargest(500).index
cg = sns.clustermap(mat.loc[top_var], z_score=0, cmap="vlag", center=0)

State the scaling in the caption. "Row z-scored log₂ CPM, top 500 variable genes" is a complete description; "expression heatmap" is not.

2. Colormap

cmap="vlag", center=0        # z-scores, log fold changes, correlations
cmap="rocket"                 # raw expression, counts, coverage
cmap="vlag", vmin=-1, vmax=1  # correlation, full range

Diverging with center=0 for anything signed. Without center=, if your z-scores run −2 to +6 the colormap's white point lands at +2, and every unchanged gene renders as up-regulated. See Colormaps and Color.

sns.clustermap(mat, z_score=0, cmap="vlag", center=0, vmin=-3, vmax=3)
sns.heatmap(mat, robust=True)      # clip to the 2nd–98th percentile

Clipping stops one extreme value from consuming the whole colour range. robust=True does it automatically — useful, and worth mentioning in the caption.

3. Clustering

method="average"    # UPGMA — the usual default
method="ward"       # compact clusters; requires Euclidean
method="complete"   # maximum linkage; tight clusters
metric="correlation"  # pattern similarity, magnitude-insensitive
metric="euclidean"    # magnitude matters
metric="cosine"

correlation for expression patterns — two genes with the same shape at different magnitudes cluster together, which is usually what "co-expressed" means. euclidean for already-z-scored data (where it is nearly equivalent).

Note method="ward" is only valid with metric="euclidean"; seaborn will let you combine it with others and the result is not meaningful.

Hierarchical clustering always produces a dendrogram, whether or not there is structure. Random data yields a confident-looking tree. Do not read cluster boundaries out of a heatmap without independent support.

Annotation bars

condition_lut = dict(zip(meta["condition"].unique(),
                         sns.color_palette("colorblind", meta["condition"].nunique())))
batch_lut = dict(zip(meta["batch"].unique(),
                     sns.color_palette("Set2", meta["batch"].nunique())))

col_colors = pd.DataFrame({
    "condition": meta["condition"].map(condition_lut),
    "batch": meta["batch"].map(batch_lut),
}, index=meta.index)

cg = sns.clustermap(mat, z_score=0, cmap="vlag", center=0, col_colors=col_colors)

from matplotlib.patches import Patch
handles = [Patch(facecolor=c, label=l) for l, c in condition_lut.items()]
cg.ax_col_dendrogram.legend(handles=handles, ncol=3, loc="center",
                            bbox_to_anchor=(0.5, 1.1), frameon=False)

Annotation bars are the most informative part of the figure. If samples cluster by batch rather than by condition, the coloured bar shows it immediately — and that is a result about your experiment, not your biology. Always include batch alongside condition.

Correlation heatmap for QC

corr = np.corrcoef(log_cpm, rowvar=False)          # sample × sample
corr = pd.DataFrame(corr, index=samples, columns=samples)

cg = sns.clustermap(corr, cmap="vlag", center=0, vmin=0, vmax=1,
                    figsize=(8, 8), col_colors=col_colors, row_colors=col_colors,
                    cbar_kws=dict(label="Pearson r"))

Do this before any analysis. Replicates should correlate ~0.95+; a sample that correlates poorly with its own group is a swap, a failure, or a mislabelling. Catching it here costs ten minutes; catching it after the paper is submitted costs rather more.

Note rowvar=False — with genes as rows, the default rowvar=True would compute a 20,000 × 20,000 gene correlation matrix (3.2 GB). See Linear Algebra with NumPy.

Big matrices

# 20,000 genes is unreadable — subset
top = mat.loc[mat.var(axis=1).nlargest(100).index]

sns.clustermap(top, z_score=0, yticklabels=True,
               figsize=(8, 14),
               cbar_kws=dict(label="z-score"))

Above ~100 rows, gene labels are illegible. Options: subset to a panel or the top-variable genes; drop the labels (yticklabels=False) and treat it as a texture plot; or split into several figures by module.

A heatmap of 20,000 unlabelled rows communicates "there is structure" and nothing more. That is sometimes the point — but say so, rather than pretending it is a gene-level figure.

Composing into a panel

clustermap is figure-level and creates its own Figure — it cannot take ax=. To put a heatmap in a panel, use sns.heatmap with a pre-computed order:

from scipy.cluster.hierarchy import linkage, dendrogram
from scipy.spatial.distance import pdist

Z = linkage(pdist(mat, metric="correlation"), method="average")
order = dendrogram(Z, no_plot=True)["leaves"]

fig, axes = plt.subplots(1, 2, figsize=(12, 6))
sns.heatmap(mat.iloc[order], cmap="vlag", center=0, ax=axes[0], yticklabels=False)
volcano(de, ax=axes[1])

See Figure-level vs Axes-level.

Common mistakes

• No scaling, so the plot shows only expression magnitude.
• Diverging colormap without center=0.
• Reading clusters as real without independent support.
• ward with a non-Euclidean metric.
• np.corrcoef orientation producing a 3.2 GB matrix.
• Too many rows to label, but labelled anyway.
• No annotation bars, hiding batch effects.
• Not stating the scaling and metric in the caption.
• ax= on clustermap → TypeError.
• Clustering on unfiltered data, where low-expression noise dominates.

See also

Colormaps and Color · Figure-level vs Axes-level · Seaborn · Linear Algebra with NumPy · Applied - Expression Matrices with NumPy · Applied - Quality Control Plots