# Heatmap Clustering Doesn't Match Your PCA? The most common cause of inconsistent heatmap and PCA clustering is that the two figures were never computed from the same matrix. A typical session builds a PCA on variance-stabilized counts from all genes, then builds a heatmap on the top 500 most variable genes, row-scaled, clustered with a different distance metric and a different linkage rule. Both figures are correct. They summarize different objects, so they tell different stories. The first thing to check is whether both figures derive from one transformed matrix, and whether the heatmap's column tree was computed before or after row scaling. In pheatmap, the `scale` argument is applied before clustering, so the sample tree describes z-scored genes even when the colors appear to describe the untransformed values. This single detail explains a large share of the reports of heatmap and PCA clustering looking inconsistent. Some definitions, since the rest of this guide depends on them. [Principal component analysis](/blog/data-analysis/principal-component-analysis-explained) (PCA) treats each sample as a point in gene-expression space and projects those points onto new axes, the principal components, which capture decreasing shares of the total variance. A clustered heatmap groups samples with hierarchical clustering, which computes all pairwise distances between samples and then merges the closest samples and clusters step by step under a linkage rule, drawing the result as a dendrogram, a tree in which branch length encodes distance. The two methods summarize the same data with different mathematics, so they will never match exactly. They should agree on the dominant structure. When they disagree, one of a small number of pipeline differences is responsible. This guide walks the causes of inconsistent heatmap and PCA clustering from most common to rarest, explains the mechanism behind each failure, gives a diagnostic sequence you can follow in minutes, and ends with prevention practices. This mismatch is one of the most frequently reported RNA-seq errors in exploratory analysis, and it is almost never a software bug. ## Why Heatmap and PCA Clustering Look Inconsistent PCA takes a samples-by-genes matrix, centers it, finds the directions of maximal variance, and plots each sample by its scores on the first two components. It is a global summary, dominated by whatever genes carry the most variance. Hierarchical clustering computes a distance between every pair of samples using a distance metric, then merges clusters under a linkage rule until one tree remains. It is sensitive to the metric, the linkage, and the exact values that were clustered. Because the two methods weight the matrix differently, disagreement has a structural basis even before any pipeline mistake enters. The causes of inconsistent heatmap and PCA clustering fall into seven buckets, ordered here by how often they explain the problem: 1. The two figures use different transformations of the count matrix. 2. The two figures use different gene subsets. 3. The heatmap row-scales the matrix before clustering, which changes the sample geometry. 4. The two figures use different distance metrics or linkage rules. 5. Outlier samples or low-count genes inject noise into one view. 6. Batch effects appear in one view and not the other. 7. Display artifacts, including arbitrary branch order, compressed color scales, and reading the wrong tree. One expectation to set before troubleshooting: agreement means the strongest pattern, usually the first two principal components together with the top-level splits of the dendrogram, tells the same story. Perfect identity between a scatter plot and a tree is impossible and not the goal. ## The First Check: Are Both Figures Using the Same Matrix? ![Decision path asking whether the PCA and heatmap matrices match in dimensions, sample names and gene names, with fixes for each mismatch.](/article-images/illustrations/heatmap-pca-clustering-inconsistent-same-matrix-check.webp) Before touching any parameters, verify the inputs. Print the dimensions and sample names of the matrix fed to each figure. ```r dim(mat_pca) # genes x samples fed to prcomp dim(mat_heat) # genes x samples fed to pheatmap identical(colnames(mat_pca), colnames(mat_heat)) identical(rownames(mat_pca), rownames(mat_heat)) ``` If the dimensions differ, you have your answer. Three mismatches account for most cases: 1. Different transformations. The PCA came from `vst()`, the heatmap from log2 counts per million, or one of them from raw counts. 2. Different gene subsets. The PCA used every gene, the heatmap used differentially expressed genes or the top variable genes. 3. Different samples. The heatmap quietly dropped two outliers that the PCA kept, often because the figures were made months apart by different scripts. Also recover the exact calls used to build each figure. If the heatmap call contains `scale = "row"`, the column tree was computed on z-scored genes, which is a different geometry from the one a default PCA uses. In practice, most bioinformatics troubleshooting sessions for this problem end at this step, because the two figures were built from different objects. ## Troubleshooting Table: Symptom to Fix This table is the core of the guide. Find the symptom that matches your figures, confirm the cause, and apply the fix. | Symptom | Likely cause | How to confirm | Fix | |--|--|--|--| | PCA separates groups, heatmap colors look mixed | Different gene subsets in the two figures | Print `dim()` of both matrices | Rebuild both from one matrix and one gene set | | Heatmap splits a group that PCA keeps together | Row scaling changed the sample geometry | Recompute sample distances with and without scaling | Match PCA scaling to the tree, or pass a precomputed tree | | PCA separates groups, heatmap looks flat | Color scale compressed by extreme z-scores | Check `quantile()` of the scaled matrix | Cap values at plus and minus 3, use symmetric breaks | | One sample sits alone on a long branch | Outlier library or single-linkage chaining | Compare library sizes, mapping rates, largest distances | Fix or remove with a documented QC reason, avoid single linkage | | One figure separates by batch, the other by condition | Batch effect plus different gene sets | Color both figures by batch | Model batch in the design, correct only visualization matrices | | Two heatmaps of the same data disagree | Branch rotation, arbitrary order | Compute the cophenetic correlation | Compare cluster membership, not branch order | | Clusters make no sense at all | Raw counts, TPM, or a transposed matrix | Check transformation provenance and `dim()` | Normalize, log-transform, verify orientation | | The tree has thousands of leaves | Reading the gene tree instead of the sample tree | Check which axis the tree is drawn on | The sample tree is the column tree | | A few samples swap between figures | Near-ties in the distance matrix | Print the distances among those samples | Acceptable, report cluster membership rather than order | ## Cause 1: Different Transformations of the Count Matrix Raw counts are the integer numbers of reads mapped to each gene. They carry library size, the total number of reads sequenced per sample, which can differ several-fold between lanes for purely technical reasons. Normalization removes that component. DESeq2 estimates a size factor per sample with the median-of-ratios method, and edgeR applies trimmed mean of M values, usually abbreviated TMM. After normalization, values still span orders of magnitude, so analysis proceeds on a log scale with a small pseudocount added to avoid the logarithm of zero, for example `prior.count = 2` in the edgeR `cpm()` function. Without the log transform, a handful of highly expressed genes dominates every distance, and a gene that doubles from 50 to 100 counts is invisible next to one that drifts from 10,000 to 11,000. Transcripts per million (TPM) rescales within each sample only. It cannot correct composition shifts between samples, where a few very highly expressed genes consume reads and deflate everything else, so count-based normalized values are preferred for sample-level comparisons. That preference is a long-standing consensus in the RNA-seq methods literature. Variance-stabilizing transformations go one step further. The DESeq2 package provides `vst()` and `rlog()`, which transform normalized counts so that variance stops depending on mean expression, exactly the property that makes Euclidean distances and PCA behave well across expression ranges. The DESeq2 documentation recommends `vst()` for large experiments and `rlog()` for small ones, roughly 30 samples or fewer. The edgeR equivalent for visualization is log2 counts per million (logCPM). ```r library(DESeq2) vsd <- vst(dds, blind = FALSE) mat_vst <- assay(vsd) # edgeR alternative library(edgeR) y <- DGEList(counts = counts) y <- calcNormFactors(y) mat_logcpm <- cpm(y, log = TRUE, prior.count = 2) ``` Pick one transformation, compute it once, and feed both figures from it. A PCA on vst values and a heatmap on logCPM will usually still agree on strong structure. A PCA on vst values and a heatmap on raw counts will not, because the raw-count distances measure library size and total expression rather than expression pattern. ## Cause 2: Different Gene Subsets A PCA on all genes and a heatmap on a selected subset answer different questions. Differentially expressed (DE) genes are chosen for mean differences between conditions, not for variance structure, and a 50-gene heatmap of a clean signature can show perfect separation while a genome-wide PCA shows batch as PC1. Both figures are correct summaries of their own inputs. Most heatmaps use the most variable genes, because those carry the structure worth displaying. DESeq2's `plotPCA()` function uses the 500 most variable genes by default, so a heatmap built on those same 500 genes is the natural partner figure. If your PCA used 18,000 genes and your heatmap used 500, disagreement is expected rather than suspicious. Filter low-count genes before selecting the subset, because genes with near-zero counts contribute distances dominated by zeros: ```r keep <- rowSums(counts(dds, normalized = TRUE) >= 10) >= 3 mat <- assay(vsd)[keep, ] rv <- apply(mat, 1, var) top <- names(sort(rv, decreasing = TRUE))[1:500] ``` The filter shown here, a count of at least 10 in at least 3 samples, is a common heuristic rather than a rule. edgeR offers `filterByExpr()` as a systematic alternative. Whatever you choose, apply it once and use the resulting matrix for every figure. ## Cause 3: Row Scaling Changes Sample Distances ![Worked calculation showing that after row scaling 500 genes, Euclidean distance between two samples equals the square root of 2G(1 minus r).](/article-images/illustrations/heatmap-pca-clustering-inconsistent-row-scaled-distance.webp) This cause is the most technically subtle and the most frequently missed. Row scaling replaces each gene's values with z-scores, meaning the gene's value minus the gene's mean across samples, divided by the gene's standard deviation. pheatmap applies this scaling before clustering, so the column tree is computed on the z-scored matrix. ComplexHeatmap clusters on whatever matrix you pass, so scaling is explicit there, but the same principle applies. The consequence is exact, not approximate. After row scaling G genes, the Euclidean distance between two samples equals the square root of 2G(1 minus r), where r is the Pearson correlation between the two samples across those genes. With 500 genes and a correlation of 0.9, that distance is 10.0. With a correlation of 0, it is 31.6. The column tree of a row-scaled heatmap is therefore a correlation tree, whether you intended one or not. A default PCA does not scale genes. `prcomp()` defaults to `center = TRUE, scale. = FALSE`, and DESeq2's `plotPCA()` also centers without scaling. That PCA is dominated by genes with the largest absolute variance and retains magnitude information. So a row-scaled heatmap and a default PCA live in different geometries, and they can disagree even when both use the identical 500 genes. The practical symptoms follow directly. A sample with a uniform downward shift, for example lower input or degraded RNA, shows perfectly normal colors in a row-scaled heatmap, because every gene is divided by its own spread, yet it separates from its group in an unscaled PCA. Pattern-only differences behave the reverse, visible in the scaled heatmap and diluted in the unscaled PCA. Decide which question you are asking, then make both figures answer it: - Absolute expression levels. No row scaling, Euclidean distance, PCA with `scale. = FALSE`. - Relative expression patterns. Row scaling, and PCA with `scale. = TRUE` or a correlation distance, so both views share the correlation geometry. You can also decouple display from clustering by passing a precomputed tree: ```r z <- t(scale(t(mat[top, ]))) # row-scaled, for colors # tree on the unscaled matrix, magnitude included column_hc <- hclust(dist(t(mat[top, ])), method = "ward.D2") pheatmap(mat[top, ], scale = "row", cluster_cols = column_hc) # tree on the scaled matrix, pattern only column_hc <- hclust(dist(z), method = "ward.D2") pheatmap(mat[top, ], scale = "row", cluster_cols = column_hc) ``` Either choice is defensible. Silence about the choice is what produces figures that appear to contradict each other. ## Cause 4: Distance Metric and Linkage Mismatch A distance metric defines how different two samples are. Euclidean distance on log counts includes magnitude, so a sample that is uniformly 2-fold lower looks far away. Pearson correlation distance, computed as `1 - cor()`, ignores magnitude and mean level, so it measures pattern only. Spearman correlation, the rank-based version, additionally resists a few extreme genes that would otherwise dominate the correlation. Linkage defines how clusters merge. Complete linkage joins on the farthest pair of members and produces compact clusters. Average linkage uses the mean pairwise distance and is a common compromise. Ward linkage merges the pair that minimizes within-cluster variance and tends to produce balanced trees, which can visually exaggerate separation. Single linkage joins on the closest pair and produces chaining, long straggly branches that one outlier can drag across the plot, and it should be avoided for expression data. One R-specific trap deserves its own sentence. In `hclust()`, `method = "ward.D2"` implements Ward's minimum-variance criterion on ordinary Euclidean distances, while `method = "ward.D"` expects squared distances, so `ward.D2` is the correct choice with `dist()` output. Using `ward.D` on unsquared distances silently produces a tree that does not correspond to the Ward criterion anyone expects. ```r d_euc <- dist(t(mat), method = "euclidean") d_cor <- as.dist(1 - cor(mat, method = "pearson")) d_spr <- as.dist(1 - cor(mat, method = "spearman")) hc <- hclust(d_euc, method = "ward.D2") plot(hc, labels = as.character(batch)) ``` Match the metric to the question and use the same metric in both figures. Euclidean pairs naturally with a default PCA. Correlation distance pairs with a PCA computed on scaled values. Mixing a correlation tree with an unscaled PCA, or a Euclidean tree with a scaled PCA, is a reliable way to manufacture disagreement. ## Cause 5: Outliers and Low-Count Noise An outlier is a sample whose technical quality or biological state places it far from every other sample. One bad library, with low library size, low mapping rate, few detected genes, or degraded RNA reflected in a low RNA Integrity Number (RIN), can dominate PC1 or sit alone on a long branch, and a heatmap built with single linkage will chain other samples onto it. Quantify before you remove anything: ```r libsize <- colSums(counts(dds)) detected <- colSums(counts(dds) > 10) sort(libsize) sort(detected) d <- dist(t(mat)) sort(as.matrix(d)[, 1], decreasing = TRUE)[1:5] ``` If one sample shows a library size several-fold below the rest and the fewest detected genes, the tree and the PCA are both telling you about RNA quality, not biology. The fixes, in order of preference, are to return to the raw data and repair the library if possible, to remove the sample with a documented quality-control reason, and to reduce its leverage with Spearman correlation or average linkage. Never delete a sample without a reason you can state in one sentence, because reviewers will ask for exactly that sentence. Low-count genes inject the complementary noise. A gene with counts of 0, 0, 3, 0 across samples contributes distance dominated by ties and rounding, so filter before clustering, as shown in Cause 2. ## Cause 6: Batch Effects and Confounded Structure A batch effect is a systematic technical difference between groups of samples processed together, caused by lane, library preparation date, reagent lot, or operator. Because technical variation often exceeds biological signal, PCA frequently places batch on PC1. If your heatmap uses DE genes from a model that included batch, it can show biology while the PCA shows batch, and the two figures will appear to contradict each other even though both are correct. Diagnose by coloring both figures with the batch variable: ```r pc <- prcomp(t(mat), center = TRUE, scale. = TRUE) plot(pc$x[, 1], pc$x[, 2], col = as.numeric(batch), pch = 19, xlab = "PC1", ylab = "PC2") legend("topright", legend = levels(factor(batch)), col = seq_along(levels(factor(batch))), pch = 19) ``` If samples arrange by batch in the PCA, include batch in the statistical design, for example `design = ~ batch + condition` in DESeq2, so that condition effects are estimated within batches. For visualization matrices, the limma function `removeBatchEffect()` produces batch-adjusted values suitable for heatmaps and PCA, and the ComBat function in the sva package is the other standard option. Batch-adjusted values belong in figures, not in differential expression testing, which the design already handles. One hard limit cannot be fixed by any transform. If batch is fully confounded with condition, meaning all wild-type samples sit in batch 1 and all knockout samples sit in batch 2, no computational method can separate the two, and any figure will show their mixture. The honest fixes are more samples or a merged analysis across batches, and the honest figure caption says exactly what the design can and cannot distinguish. ## Cause 7: Display Artifacts and Reading Errors The rarest causes are not statistical at all. They live in how trees are drawn and read. ### Branch order is arbitrary A dendrogram is defined only up to rotation. At every merge, the two branches could be drawn in either order, and heatmap functions rotate branches so that the tree does not cross the matrix. Two heatmaps of the same distance matrix can therefore look different while encoding identical structure. Compare cluster membership after cutting the tree, for example with `cutree(hc, k = 3)`, or check how faithfully the tree represents the distances: ```r cor(as.vector(cophenetic(hc)), as.vector(d)) ``` The `cophenetic()` function returns the distances as represented by the tree, and a high correlation with the original distances means the tree is faithful. Branch order itself carries no information. ### Reading the wrong tree A clustered heatmap carries two trees. The row tree clusters genes, the column tree clusters samples. A surprising number of reported inconsistencies dissolve once the reader looks at the column tree. If your tree has thousands of leaves, you are looking at the gene tree. ### Color scale compression A few extreme z-scores, for example a gene at plus 12 in one sample, stretch the color scale so that the bulk of the cells, sitting between minus 2 and plus 2, all receive nearly the same color. Real structure becomes invisible even though the tree still shows it. Cap the values and use symmetric breaks: ```r z <- t(scale(t(mat[top, ]))) z[z > 3] <- 3 z[z < -3] <- -3 pheatmap(z, cluster_cols = column_hc, border_color = NA) ``` Capping is a display decision and should be stated in the caption, because it changes what the colors claim to show. ### Magnitude hidden by row scaling Row scaling erases absolute level, so a sample that is uniformly lower than the rest shows no color difference at all. If the biology you care about is a global shift in output, the heatmap is the wrong figure for it, and the unscaled PCA is the right one. This is not an artifact to fix but a property to respect. ## How to Fix Heatmap and PCA Clustering Inconsistency: A Diagnostic Sequence Work through these steps in order. Each one either resolves the problem or eliminates a cause, and changing one parameter at a time is what makes bioinformatics troubleshooting conclusive rather than guesswork. 1. **Inventory the inputs.** Print `dim()` and compare `colnames()` for the matrix behind each figure. Any mismatch ends the investigation. 2. **Fix the transformation.** Compute one normalized, transformed matrix, preferably `vst()` output, and save it as a single object. 3. **Fix the gene set.** Apply one low-count filter and one selection rule, for example the top 500 variable genes, and subset once. 4. **Recompute the PCA from the heatmap matrix exactly**, including whatever scaling the heatmap applies, so both figures describe the same object. 5. **Quantify agreement numerically** rather than by eye, using the code below. 6. **Check the variance explained.** If PC1 and PC2 together explain 30 percent of the variance, the scatter plot discards 70 percent of the structure, and disagreement with any tree is unsurprising. 7. **Inspect the extreme distances** and the samples involved in them, along with library size and QC metrics. 8. **Change one parameter at a time**, record it, and re-run both figures. ```r pc <- prcomp(t(mat), center = TRUE, scale. = TRUE) d <- dist(t(mat)) cor(as.vector(d), as.vector(dist(pc$x[, 1:2]))) # PCA geometry vs distances cor(as.vector(cophenetic(hc)), as.vector(d)) # tree vs distances summary(pc)$importance[2, 1:2] # variance explained by PC1 and PC2 ``` As a rule of thumb, a correlation above roughly 0.8 between the PCA distances and the clustering distances means the figures describe the same geometry. Below that, one of the causes in this article is active, and the sequence above will locate it. ## Worked Example: One Matrix, Two Views That Agree The following script builds both figures from one object so that agreement is guaranteed by construction. It assumes a DESeq2 object `dds` with columns `condition` and `batch` in its column metadata. ```r library(DESeq2) library(pheatmap) ## one transformation, computed once vsd <- vst(dds, blind = FALSE) ## one gene set: drop low counts, keep the 500 most variable genes keep <- rowSums(counts(dds, normalized = TRUE) >= 10) >= 3 mat <- assay(vsd)[keep, ] top <- names(sort(apply(mat, 1, var), decreasing = TRUE))[1:500] ## PCA on exactly those samples and genes, unscaled genes pc <- prcomp(t(mat[top, ]), center = TRUE, scale. = FALSE) var_pc <- 100 * pc$sdev^2 / sum(pc$sdev^2) plot(pc$x[, 1], pc$x[, 2], col = as.numeric(condition), pch = 19, xlab = paste0("PC1 (", round(var_pc[1]), "%)"), ylab = paste0("PC2 (", round(var_pc[2]), "%)")) ## sample tree on the same matrix, same geometry column_hc <- hclust(dist(t(mat[top, ])), method = "ward.D2") ## heatmap: row scaling for colors only, tree passed in explicitly pheatmap(mat[top, ], scale = "row", cluster_cols = column_hc, clustering_method = "ward.D2", annotation_col = as.data.frame(colData(dds)[, c("condition", "batch")]), border_color = NA) ``` These figures agree because the column tree and the PCA both use Euclidean distances among samples computed on the same 500 genes, and the tree is passed to the heatmap explicitly so that row scaling affects colors only. If you want the tree to reflect expression pattern instead of magnitude, compute it on the row-scaled matrix and switch the PCA to `scale. = TRUE`. Either way, the caption should state the transformation, the gene set, the distance, and the linkage. *Safety note: this workflow is fully computational and involves no reagents or biological hazards. If the count data derive from clinical specimens, work from de-identified data and follow your institutional biosafety and data governance policy.* ## Common Mistakes and Limitations - **Comparing two different questions.** A heatmap of DE genes and a PCA of all genes are not two views of one analysis. Calling their difference a contradiction is the single most common mistake in this area. - **Trusting branch order.** Left-to-right position in a dendrogram is arbitrary. Only cluster membership and branch length carry meaning. - **Reading colors as distances.** Colors show per-gene z-scores. The distances live in the tree. A flat-looking heatmap with a clear tree is not a contradiction. - **Ignoring variance explained.** If the first two components capture a small share of the variance, the scatter plot is a weak summary, and demanding that a tree match it exactly is unreasonable. - **Clustering on raw counts or TPM.** Raw counts measure library size as much as biology, and TPM cannot correct composition shifts between samples. - **Row scaling and then claiming magnitude differences.** After row scaling, a uniformly lower sample is invisible by construction. - **Treating either figure as a test.** Neither PCA nor a heatmap is a statistical test of group separation. For formal claims, use a permutation-based test such as PERMANOVA, implemented as `adonis2()` in the vegan package, or a cross-validated classifier. - **Forgetting that parameters are choices.** Transformation, filter, distance, and linkage are analyst decisions that shape both figures. Report them, and interpret the results with someone who knows the biology of the system. ## Prevention Practices 1. **Build one analysis matrix at the start** and derive every figure from it. One object, one transformation, one gene filter. 2. **Script everything.** No manual subsetting, no clicking in the plot pane, no figures made months apart by different scripts without reconciliation. 3. **Record the parameters in every caption**: transformation, gene set and size, distance metric, linkage, and any capping applied to the color scale. 4. **Pin your software versions** with `sessionInfo()` or an environment manager such as renv, so the figures can be regenerated exactly. 5. **Run the numeric agreement check** from the diagnostic sequence before any submission or presentation, not after a reviewer asks. 6. **Predefine the gene selection rule** in the project README, so the top-500 choice is a decision rather than a default nobody examined. ## Frequently Asked Questions ### Why does my heatmap show clusters that my PCA does not? Because the two figures usually summarize different matrices, most often a heatmap built on a small set of row-scaled genes and a PCA built on all genes. Rebuild both from one transformed matrix and the dominant structure will usually match. ### Should PCA and the heatmap use the same gene set? Yes, if you want the figures to agree. A heatmap of the 500 most variable genes answers a different question than a PCA on 15,000 genes. Either restrict both figures to the same set or state each figure's gene set in its caption. ### Is it wrong to cluster samples on row-scaled heatmap values? It is not wrong, but it changes the question. Row scaling makes Euclidean distance between samples a function of Pearson correlation, so the column tree matches correlation-based PCA rather than a default PCA on log counts. ### Which distance metric should I use for RNA-seq sample clustering? Euclidean distance on variance-stabilized log counts is the default and pairs naturally with PCA. Use Pearson correlation distance when expression pattern matters more than magnitude, and Spearman correlation when you want robustness to a few extreme genes. ### Why do my heatmap colors look flat even though PCA separates the groups? A few extreme z-scores compress the color scale, so most cells land in a narrow band. Cap values at plus and minus 3 with symmetric breaks and the structure usually reappears. ### Can batch effects make heatmap and PCA disagree? Yes, when the two figures use different gene sets that respond differently to batch. Color both figures by batch, include batch in the statistical design, and use batch-corrected values only for visualization. ### Does the dendrogram branch order mean anything? No. Branch order within a cluster is arbitrary, and heatmap functions rotate branches to avoid crossings. Compare cluster membership after cutting the tree, not the left-to-right order. ### How many genes should I use for the heatmap? There is no universal number. A common choice is the 500 to 2,000 most variable genes after low-count filtering, large enough to carry the structure you claim and small enough to render legibly.