Heatmap maker with hierarchical clustering

Row z-scores. Each value becomes $z=(x-\overline{x})/s$, where $\overline{x}$ and $s$ are the mean and sample standard deviation ($n-1$) of its row, the same as R’s scale() and pheatmap’s scale = "row". Every row then has mean 0 and SD 1, so the colours show whether a gene is high or low in a sample relative to its own average, not its absolute level. Rows with no variation get $z=0$.

Distances. Euclidean distance $d=\sqrt{\sum (x_i-y_i{)}^2}$ groups rows with similar values; $1-\text{Pearson }r$ (0 for identical patterns, 1 for unrelated, 2 for opposite) groups rows with the same shape of profile regardless of level, the measure used by Eisen et al. (1998). On z-scored rows the two are closely related: the squared Euclidean distance between z-scored rows is $2(n-1)(1-r)$, so they rank pairs of rows in the same order and usually give very similar trees. With missing cells, distances use the positions both rows have, and Euclidean distance is scaled by $\sqrt{p/m}$ as in R’s dist().

Average linkage (UPGMA). Start with every row as its own cluster; merge the two closest; set the distance from the merged cluster to every other as the size-weighted mean of the two old distances, $d(A\cup B,C)$ $=(|A|d(A,C)+|B|d(B,C))$ $/(|A|+|B|)$; repeat until one cluster is left. The dendrogram joins clusters at the height of their merge, and the heatmap rows follow its leaf order. Merge heights and leaf order match SciPy’s linkage(method="average") and dendrogram on the example.

Colours. Sequential data (expression, intensity) use viridis or cividis, which are perceptually uniform and readable with colour-vision deficiency and in greyscale; data around a centre (z-scores, log fold changes) use the ColorBrewer RdBu or PuOr diverging scales, centred so 0 is white. Rainbow and red–green scales distort the data and exclude colour-blind readers (Crameri et al. 2020). Sequential scales run from the smallest to the largest value; diverging scales are symmetric around the centre.

Sources: Eisen MB, Spellman PT, Brown PO, Botstein D. Cluster analysis and display of genome-wide expression patterns. Proc Natl Acad Sci USA 1998;95:14863–14868, doi:10.1073/pnas.95.25.14863. Crameri F, Shephard GE, Heron PJ. The misuse of colour in science communication. Nat Commun 2020;11:5444, doi:10.1038/s41467-020-19160-7.

The built-in matrix is illustrative log₂ expression of 16 genes in 4 control and 4 hypoxia samples: hypoxia-induced genes (VEGFA, CA9, SLC2A1, PGK1, BNIP3, LDHA, NDRG1), proliferation genes that fall (MKI67, CCNB1, TOP2A, CDK1, PLK1) and four housekeeping genes.

1. Z-score a row. VEGFA: 6.1, 6.3, 5.9, 6.2, 8.4, 8.9, 8.6, 8.2. Mean 7.325, SD 1.302, so $z=-0.94,-0.79,-1.09,-0.86,$ $0.83,1.21,0.98,0.67$. PGK1 is about 3 log₂ units higher but gets similar z-scores, so the two genes look alike on the heatmap.
2. Cluster the rows. With z-scores, Euclidean distance and average linkage, the seven hypoxia genes form one tight cluster (all joined by height 0.45) and the five proliferation genes another (by 0.14). The housekeeping genes vary by only 0.1–0.4 log₂ units, but z-scoring stretches that noise to an SD of 1, so each attaches to the group its small trend follows: GAPDH (11.5 → 11.9) joins the hypoxia genes at 0.67, ACTB and RPLP0 join the proliferation genes at 3.06, and TBP joins the hypoxia side at 3.47; the last merge is at 4.94. This is the usual caveat of row scaling: remove flat genes (low variance) before plotting.
3. Cluster the columns. The samples split into the four controls and the four hypoxia samples, joined at height 7.09; within groups the first merge is Ctrl_1 with Ctrl_4 at 0.47.
4. Without z-scores, Euclidean clustering groups genes partly by level instead: the highly expressed ACTB and GAPDH join at 1.25 and then take in RPLP0, PGK1 and LDHA, so the hypoxia genes are split by how strongly they are expressed. That is why expression heatmaps are usually row-scaled.