PCA and Sample Clustering for RNA-Seq in R (DESeq2)

Q: VST or rlog — which transformation should I use?

Both stabilize the variance for visualization. vst() is fast and the default for typical studies; rlog() is more robust on very small datasets (a handful of samples) but much slower on large ones. For QC plots either is fine — use VST unless you have very few samples.

Q: Should blind be TRUE or FALSE in vst()?

For quality control use blind = TRUE: the transformation ignores the design, giving an unbiased look at whether samples group as expected. For transformed values you will feed into downstream modelling, use blind = FALSE so the design-related variance is preserved.

Check your samples before differential expression — variance-stabilize, run PCA, and cluster samples to spot grouping, outliers, and batch effects

Bioinformatics

A practical guide to exploratory quality control for bulk RNA-seq in R with DESeq2. Before testing for differential expression, check that your samples make sense: apply a variance-stabilizing transformation (VST), draw a PCA plot to see whether samples separate by condition (and to catch outliers or batch effects), and a sample-to-sample distance heatmap to see how they cluster. Read the percent-variance axes in plain language. On the airway dataset.

Published

June 29, 2026

Modified

July 7, 2026

Key takeaways

Look at your samples before you test them. A quick PCA + clustering check catches mislabelled samples, outliers, and batch effects before they wreck your differential expression.
Variance-stabilize first. Raw counts have variance that grows with the mean, so a few loud genes would dominate. DESeq2’s VST (vst()) puts all genes on a comparable scale for the visual steps.
PCA projects the samples onto the directions of greatest variation. If your biological groups separate along PC1 or PC2, the signal is strong; if a batch separates them instead, you have a problem to model.
A sample-to-sample distance heatmap is the complementary view: it clusters samples so you can see which group together and spot any that sit alone.
These are QC steps, not the test — they run on transformed counts; the differential expression test runs on the raw counts (DESeq2 lesson).

Introduction

Before asking which genes changed, ask a simpler question: do my samples make sense? Sequencing experiments are full of avoidable disasters — a swapped label, a failed library, a batch processed on a different day. A two-minute exploratory check catches them before they corrupt every downstream result.

The check has two standard views, both built on a variance-stabilized version of the counts: a PCA plot (do samples separate by condition?) and a sample-to-sample distance heatmap (how do they cluster?). This guide runs both in R with DESeq2 on the airway dataset — airway smooth-muscle cells from four donors, each treated and untreated with dexamethasone.

Start from the count data

We build the DESeq2 object from the raw counts — the same starting point as the DESeq2 lesson.

library(airway)
library(DESeq2)

data("airway")
airway$dex <- relevel(airway$dex, ref = "untrt")
dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds <- dds[rowSums(counts(dds)) >= 10, ]
dds

class: DESeqDataSet 
dim: 22369 8 
metadata(2): '' version
assays(1): counts
rownames(22369): ENSG00000000003 ENSG00000000419 ... ENSG00000273487
  ENSG00000273488
rowData names(10): gene_id gene_name ... seq_coord_system symbol
colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
colData names(9): SampleName cell ... Sample BioSample

Stabilize the variance (VST)

Raw counts are a poor input for PCA and clustering: their variance grows with the mean, so a handful of highly expressed genes would dominate purely because they are big. The variance-stabilizing transformation (vst) fixes this — it puts low- and high-count genes on a comparable scale, like a smarter log2.

vsd <- vst(dds, blind = TRUE)     # blind = TRUE: ignore the design, a pure QC view
head(assay(vsd))[, 1:4]

                SRR1039508 SRR1039509 SRR1039512 SRR1039513
ENSG00000000003   9.464955   9.085088   9.613950   9.349298
ENSG00000000419   8.963799   9.271415   9.154078   9.199673
ENSG00000000457   8.211542   8.117807   8.051019   8.165650
ENSG00000000460   6.610550   6.652804   6.153025   6.507349
ENSG00000000971  11.652764  12.018678  12.366560  12.641598
ENSG00000001036  10.495488  10.263000  10.563194  10.407888

We set blind = TRUE so the transformation does not use the experimental design — exactly what you want for an unbiased quality-control look. (For downstream modelling you would use blind = FALSE.) For very small studies, rlog() can be more robust; for typical and larger experiments, VST is faster and the standard choice.

PCA — do samples separate by condition?

Principal component analysis projects each sample onto the directions that capture the most variation. PC1 is the biggest source of variance, PC2 the next. If your biological grouping lines up with PC1 or PC2, the effect is strong relative to the noise.

library(ggplot2)

pcaData <- plotPCA(vsd, intgroup = c("dex", "cell"), returnData = TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))

ggplot(pcaData, aes(PC1, PC2, color = dex, shape = cell)) +
  geom_point(size = 4) +
  scale_color_viridis_d(end = 0.85) +
  labs(x = paste0("PC1: ", percentVar[1], "% variance"),
       y = paste0("PC2: ", percentVar[2], "% variance"),
       color = "Treatment", shape = "Donor") +
  coord_fixed() +
  theme_minimal()

PCA plot of the variance-stabilized counts. PC1 (41% of variance) on the x-axis separates the untreated samples (left) from the treated samples (right); PC2 (26%) on the y-axis spreads the four donor cell lines. Points are coloured by treatment and shaped by donor, showing a clean treatment split along PC1.

Read the plot. PC1 captures 41% of the variance and cleanly separates treated from untreated samples (left vs right) — the dexamethasone effect is the single biggest source of variation, a reassuring sign before testing. PC2 (26%) captures donor-to-donor variation (most visibly one donor high, another low), the person-to-person differences we control for with cell in the design. No sample sits off on its own, so there are no obvious outliers.

If a batch separates your samples on PC1

The danger sign is when PC1 lines up with a technical variable — sequencing day, lane, kit — instead of your biology. That means a batch effect is the dominant signal. The fix is to add the batch to your design (~ batch + condition) so DESeq2 accounts for it, not to ignore it. PCA is how you catch it.

Sample-to-sample distance heatmap

PCA summarizes the samples in two dimensions; a distance heatmap shows the full pairwise picture. Compute the distances between samples on the VST data and cluster them — samples from the same group should sit together.

library(pheatmap)

sampleDists <- dist(t(assay(vsd)))
mat <- as.matrix(sampleDists)
rownames(mat) <- paste(vsd$dex, vsd$cell, sep = " · ")
colnames(mat) <- NULL

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

Heatmap of sample-to-sample distances on the variance-stabilized counts, with samples hierarchically clustered on both axes. Darker cells are more similar; one donor's two samples (N080611) cluster together and split off from the rest first, while treatment drives the split among the remaining samples.

The dendrogram tells a subtler story than PCA did. One donor (N080611) dominates: its treated and untreated samples cluster together and split off from the rest first; only among the remaining six does treatment then drive the split. Why the difference from PCA? The PCA used the 500 most variable genes — where the treatment effect stands out — while these distances use all genes, so a single strongly-different donor can dominate the overall picture. That the two views emphasize different structure is itself informative: it confirms the donor effect is real and large, which is exactly why the design controls for it with ~ cell + dex. Always read both — they answer slightly different questions.

Common issues

Samples don’t separate by condition. Either the effect is genuinely small (GSEA on the next step may still find coordinated shifts), or a stronger nuisance variable dominates. Colour the PCA by other variables (batch, RIN, donor) to find what PC1 is really tracking.

One sample sits far from the rest. A likely outlier — a failed library or a swap. Check its library size and QC metrics; consider dropping it, but never silently. Document the decision.

You ran PCA on raw or plain-log counts. Don’t. Without VST (or rlog), high-count genes dominate the variance and the plot reflects sequencing depth more than biology. Always transform first.

Frequently asked questions

VST or rlog — which transformation should I use?

Both stabilize the variance for visualization. vst() is fast and the default for typical studies; rlog() is more robust on very small datasets (a handful of samples) but much slower on large ones. For QC plots either is fine — use VST unless you have very few samples.

Should blind be TRUE or FALSE in vst()?

For quality control use blind = TRUE: the transformation ignores the design, giving an unbiased look at whether samples group as expected. For transformed values you will feed into downstream modelling, use blind = FALSE so the design-related variance is preserved.

Why not just run PCA on the raw counts?

Because count variance grows with the mean: a few highly expressed genes would dominate the principal components, and the plot would reflect sequencing depth more than biology. The VST puts all genes on a comparable scale first, so PCA sees real structure.

What if PC1 separates by batch instead of condition?

That is a batch effect dominating your data. Don’t ignore it — add the batch variable to the DESeq2 design (~ batch + condition) so the test accounts for it. PCA is exactly the tool that reveals the problem before you test.

Test your understanding

Your turn: which variable drives PC2?

The PCA above colours by treatment. Re-plot it colouring by donor (cell) instead, and decide which principal component the donor effect lines up with. Is it PC1 or PC2?

Hint

Reuse pcaData; just change the color = aesthetic to cell. Look at which axis spreads the donors apart.

Solution

ggplot(pcaData, aes(PC1, PC2, color = cell, shape = dex)) +
  geom_point(size = 4) +
  coord_fixed() + theme_minimal()

The donors spread along PC2, not PC1 — treatment is the dominant axis (PC1, 41%), donor the secondary one (PC2, 26%). That is exactly why the design is ~ cell + dex: it removes the donor variance so the test can see the treatment effect cleanly.

Quick check: why transform counts before PCA?

A. To make the test more significant. B. Because raw-count variance grows with the mean, so big genes would dominate the components. C. PCA cannot handle integers.

Show answer

B. Count variance scales with expression level, so without a variance-stabilizing transform a few highly expressed genes drive the principal components and the plot tracks sequencing depth, not biology.

Conclusion

You ran the standard exploratory check for bulk RNA-seq: variance-stabilize the counts (vst), then look at the samples two ways — a PCA plot (here PC1 captured 41% of the variance and cleanly split treated from untreated, with donor variation on PC2) and a sample-distance heatmap (which, using all genes, was dominated by one donor — a reminder that the two views can emphasize different structure). The practical rule: always do this before differential expression — it is the cheapest way to catch a swapped label, an outlier, or a batch effect, and it tells you whether your effect is strong enough to trust. With the samples checked, you are ready to test for differential expression.

References

Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12), 550.
Love, M. I., Anders, S., & Huber, W. Analyzing RNA-seq data with DESeq2 — the Bioconductor vignette (data transformations & quality assessment).
Himes, B. E., et al. (2014). RNA-Seq transcriptome profiling identifies CRISPLD2 as a glucocorticoid responsive gene in airway smooth muscle cells (the airway dataset). PLoS ONE, 9(6), e99625.

Reuse

CC BY-NC-SA 4.0

Citation

BibTeX citation:

@online{2026,
  author = {},
  title = {PCA and {Sample} {Clustering} for {RNA-Seq} in {R}
    {(DESeq2)}},
  date = {2026-06-29},
  url = {https://www.datanovia.com/learn/bioinformatics/bulk-rna-seq/pca-and-sample-clustering},
  langid = {en}
}

For attribution, please cite this work as:

“PCA and Sample Clustering for RNA-Seq in R (DESeq2).” 2026. June 29. https://www.datanovia.com/learn/bioinformatics/bulk-rna-seq/pca-and-sample-clustering.