Heatmaps and Clustering

Before you start

You've done From Raw Counts to Differential Expression, so you have a dds object and a results table.
A little R, and a free Google account to run the Colab notebook.

Learning objectives

By the end of this lesson you will be able to: build a heatmap of expression in R, understand how clustering groups genes and samples, and choose scaling and gene selection so the figure shows real structure rather than noise.

First rule: never heatmap raw counts

It's tempting to throw your count matrix straight into a heatmap, but it will mislead you. Raw counts span an enormous range (some genes have tens of reads, others tens of thousands), and they depend on how deeply each sample was sequenced. A heatmap of raw counts just shows you a handful of very high-count genes and hides everything else.

The fix is a variance-stabilizing transformation. It puts every gene on a comparable, roughly even scale so that real patterns, not raw magnitude, drive the colours. DESeq2 gives you one in a single line.

Decode the jargonVariance-stabilizing transformation (vst)

A transformation that flattens the relationship between a gene's average count and its variability, so highly expressed genes stop dominating. The result is values you can safely cluster, plot, and compare across samples. DESeq2's vst() is fast and works well for most datasets; for very small ones, rlog() is a gentler alternative.

▶ Code along in Google Colab

Set the runtime to R (Runtime, then Change runtime type, then R), and run each block as you read. The DESeq2 install takes about 3 to 5 minutes the first time, which is completely normal and only happens once per session.

Open the R notebook in Colab →

1Set up and transform

Rebuild the dataset from the previous lessons, then transform it. We'll also install pheatmap, a small, friendly package for drawing heatmaps.

library(DESeq2)
install.packages("pheatmap")   # small CRAN package, installs quickly
library(pheatmap)

base <- "https://raw.githubusercontent.com/owkin/PyDESeq2/main/datasets/synthetic/"
counts <- as.matrix(read.csv(paste0(base, "test_counts.csv"), row.names = 1))
metadata <- read.csv(paste0(base, "test_metadata.csv"), row.names = 1)
metadata$condition <- factor(metadata$condition)
metadata <- metadata[colnames(counts), , drop = FALSE]

dds <- DESeqDataSetFromMatrix(counts, metadata, design = ~ condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "B", "A"))

# the transformed values we'll use for every plot below
vsd <- vst(dds, blind = FALSE)
# tip: on a very small dataset, use rlog(dds) instead if vst complains

2Do my samples cluster by condition? (the QC heatmap)

Before trusting any result, ask a basic question: do replicates of the same condition look alike? We answer it by measuring the distance between every pair of samples and drawing those distances as a heatmap.

sampleDists <- dist(t(assay(vsd)))   # distance between samples
mat <- as.matrix(sampleDists)

pheatmap(mat,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         main = "Sample-to-sample distances")

How to read it: darker (smaller distance) means two samples are more similar. You want samples of the same condition to sit together in tight blocks, with a clear split between groups. If one sample clusters with the wrong group, that's a red flag worth investigating: a possible sample swap, a batch effect, or an outlier. This single plot catches problems that would otherwise quietly ruin your conclusions.

3The same question, even faster: PCA

A PCA plot squeezes all the genes down to two axes that capture the biggest sources of variation, so you can see the grouping in one glance. DESeq2 has it built in.

plotPCA(vsd, intgroup = "condition")

Samples of the same condition should cluster together and the two conditions should separate, ideally along the first axis (PC1). PCA and the distance heatmap answer the same question two ways; seeing both agree is reassuring.

4Which genes separate the groups? (the gene heatmap)

Now the satisfying one. We take the most differentially expressed genes, scale each gene so we see its pattern rather than its absolute level, and let clustering arrange them.

# the 30 most significant genes
top <- head(order(res$padj), 30)
mat <- assay(vsd)[top, ]

# z-score each gene across samples (centre and scale per row)
mat <- t(scale(t(mat)))

# label the columns by condition
ann <- as.data.frame(colData(vsd)[, "condition", drop = FALSE])

pheatmap(mat,
         annotation_col = ann,
         show_rownames = FALSE,
         main = "Top 30 genes (z-scored)")

Read it as blocks of colour. A clean result shows two opposing blocks: a set of genes high (one colour) in condition A and low in B, and another set doing the reverse. The dendrogram across the top should group the samples by condition on its own, which is a lovely independent confirmation that your differential-expression hits are real and coherent.

Decode the jargonZ-score scaling

For each gene, we subtract its mean and divide by its standard deviation across samples. This rescales every gene to the same footing so the heatmap shows relative change (up or down versus that gene's own average) instead of which genes are simply loud. Without it, a few high-expression genes would wash out the pattern you care about.

🔶 Level up: prettier and bigger heatmaps

pheatmap has many options worth exploring: custom colour palettes, splitting rows or columns into groups with cutree_rows, and adding multiple annotation tracks. For complex, publication-grade figures, the ComplexHeatmap Bioconductor package is the community standard once you outgrow the basics.

Check your understanding

Why transform the counts with vst before making a heatmap?

Raw counts span a huge range and depend on sequencing depth, so a few high-count genes would dominate the colours. The variance-stabilizing transform puts genes on a comparable scale, so real patterns drive the heatmap instead of raw magnitude.

What does the sample-distance heatmap tell you?

Whether replicates of the same condition are actually similar to each other. Tight same-condition blocks with a clean split between groups is healthy. A sample sitting with the wrong group hints at a swap, a batch effect, or an outlier.

Why z-score each gene before the gene heatmap?

So the heatmap shows each gene's relative change (up or down versus its own average) rather than its absolute expression level. Otherwise a handful of very highly expressed genes would dominate and hide the pattern.

Sources & further reading

Love MI, Huber W, Anders S. DESeq2 (variance-stabilizing transformation used before clustering). Genome Biology 15:550, 2014. doi:10.1186/s13059-014-0550-8
Conesa A, et al. A survey of best practices for RNA-seq data analysis. Genome Biology 17:13, 2016. doi:10.1186/s13059-016-0881-8

Last reviewed: June 2026.

After z-scoring a gene, what does a value of 0 mean for one sample?

A z-score of 0 means that sample sits exactly at the gene's mean. Positive values are above average (higher than usual) and negative values are below, which is why heatmap colours read as relative up or down.

What does hierarchical clustering do to the rows and columns of a gene heatmap?

Clustering groups genes with similar expression patterns together and samples that look alike together, so blocks of co-regulated genes and related samples become visible as clear patches of colour.

Next in Track 3

Heatmaps and Clustering That Mean Something