Applied - Heatmaps and Clustermaps

Core Principle: Clustered heatmaps combine hierarchical distance-based clustering of samples and genes with standardized color-mapped intensity matrices to reveal co-expression patterns.

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

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

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

Exercise

seaborn and scipy are both absent from this sandbox, so build the clustermap's clustering step yourself: z-score the rows, take the correlation distance, and run average linkage with NumPy. Doing it by hand is what makes the three decisions in the lesson visible.

Write cluster(text), returning a list of output lines. The first input line holds the number of clusters k. The second is the header, the literal gene followed by the sample names. Every later line is a gene name and its values. Z-score each row using the population standard deviation; a row with zero standard deviation cannot be scaled, so set it to all zeros and count it. With Z the z-scored matrix over m samples, the correlation matrix is Z @ Z.T / m, which gives a constant row a correlation of zero with everything; the distance is 1 - r.

Cluster with average linkage: repeatedly merge the two clusters with the smallest mean pairwise distance between their members until k clusters remain, breaking ties in favour of the pair earliest in the current cluster list. Treat a k at or above the number of genes as "no merging at all", and a k below one as one. Print genes=<n> constant=<n>, then one line per cluster as cluster<i>=<names>, where the member names are sorted ascending, joined with commas, and the clusters themselves are numbered from one in ascending order of their first member name.