How to Make and Read a PCA Plot for RNA-seq Data

Oak
rna-seq pca deseq2 quality-control r
A studio-lit iridescent feathered creature on a black plinth, surrounded by two drifting clusters of small glowing beads.

By the end of this you will have a scree plot, a two-dimensional scatter of your samples with the percent variance printed on each axis, a table of which metadata columns explain which principal components, and a ranked list of the genes driving each one. You need gene-level counts (a salmon/kallisto quantification directory or a featureCounts matrix), a sample metadata table with one row per library, and R with DESeq2, tximport, ggplot2, and limma installed. Everything below assumes bulk RNA-seq with somewhere between 6 and a few hundred libraries, which covers most personal longitudinal designs (one person, repeated blood draws) as well as ordinary case-control studies. A PCA plot answers one question well: how far apart are my samples in overall expression space, and does that distance line up with anything I recorded? That is a quality control question before it is a biology question, and it should be the first figure you make after quantification 1.

1. Build the count matrix and the metadata table together

PCA is only as informative as the metadata you can color it by, so assemble both objects in the same script and keep the column order locked. If you quantified with Salmon, import to gene level with tximport and a transcript-to-gene map:

library(tximport); library(DESeq2); library(readr)

meta <- read_csv("samples.csv")           # sample_id, draw_date, batch, rin, condition, ...
meta$batch <- factor(meta$batch)
files <- file.path("quants", meta$sample_id, "quant.sf")
names(files) <- meta$sample_id
stopifnot(all(file.exists(files)))

txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ 1)

Use countsFromAbundance = "lengthScaledTPM" so that changes in transcript-length composition between samples do not masquerade as expression changes. The design = ~ 1 is deliberate: for an unsupervised look you do not want the design influencing the dispersion estimates used by the transform. Include in samples.csv every nuisance variable you have, even the ones you assume are irrelevant: library prep batch, flow cell, sequencing date, RIN (RNA integrity number, a 1–10 degradation score from the Bioanalyzer), input mass, draw time of day, and for whole blood, the globin depletion method. You will regret each one you left out at step 5.

2. Filter low counts, then transform

Raw counts are unusable for PCA. Counts are heteroscedastic: a gene averaging 20,000 reads has an enormous raw variance compared to a gene averaging 20, so principal components computed on counts are dominated by a handful of highly expressed genes (in whole blood, hemoglobin and a few ribosomal genes) regardless of biology. Log transformation overcorrects in the other direction, inflating the apparent variance of low-count genes where a jump from 1 to 4 reads is noise. The variance-stabilizing transform in DESeq2 fits the mean-variance relationship across all genes and applies a transform that flattens it:

keep <- rowSums(counts(dds) >= 10) >= max(3, floor(0.2 * ncol(dds)))
dds  <- dds[keep, ]

vsd <- vst(dds, blind = TRUE, nsub = 1000)   # or rlog(dds, blind = TRUE) if n < 30
mat <- assay(vsd)

We use vst by default. It is fast on hundreds of samples and, in our experience, gives PCA geometry nearly identical to rlog for libraries of reasonable depth. Use rlog when you have fewer than about 30 samples and their library sizes vary by more than roughly fourfold, which is where rlog handles shrinkage more gracefully. blind = TRUE means the transform ignores the design, which is what you want for an honest-to-the-data QC plot. If your counts came from featureCounts rather than tximport, build the DESeqDataSet with DESeqDataSetFromMatrix and the same filter.

There is a deeper point here that is easy to skip. Sequencing counts are compositional: a library reports proportions of a fixed read budget, not absolute molecule counts, so one gene rising forces every other gene’s share down. Compositional data analysis handles this with log-ratio transforms such as the centered log-ratio, and a CLR-transformed matrix is a defensible alternative input to PCA 2. In practice VST and CLR produce similar sample geometry unless a few genes dominate the library, which does happen in whole blood without globin depletion. If your top 10 genes hold more than about 40% of reads, run the PCA both ways and compare.

3. Compute the principal components

The one-liner in DESeq2 gets you a plot immediately:

plotPCA(vsd, intgroup = c("batch", "condition"), ntop = 500)

Under the hood plotPCA selects the ntop genes with highest row variance, transposes the matrix so samples are rows, and calls prcomp with centering on and scaling off. Do it yourself once so you control each choice:

rv     <- matrixStats::rowVars(mat)
select <- order(rv, decreasing = TRUE)[seq_len(min(2000, length(rv)))]
pca    <- prcomp(t(mat[select, ]), center = TRUE, scale. = FALSE)

pv <- pca$sdev^2 / sum(pca$sdev^2)
round(100 * pv[1:8], 1)

Three choices deserve an opinion. First, ntop: the default of 500 is a filter, and a tight one. We run the PCA at 500, 2000, and all genes, and check that the sample ordering is stable. If the picture changes qualitatively between 500 and all genes, your structure is carried by a small gene set and you should look at it directly rather than trusting the scatter. Second, center = TRUE is mandatory; without it PC1 is the mean expression vector and tells you nothing. Third, scale. = FALSE is right for a VST matrix, because the transform already put genes on a comparable variance footing, and scaling would promote near-constant genes to equal standing with genuinely variable ones. Scale only if you are feeding in something with heterogeneous units, which is the situation in multi-omic integration where proteins, metabolites, and transcripts live on different scales 3.

4. Plot it, with the variance on the axes and the aspect ratio honest

library(ggplot2)
df <- data.frame(PC1 = pca$x[,1], PC2 = pca$x[,2],
                 pca$x[,3:4], meta[match(rownames(pca$x), meta$sample_id), ])

ggplot(df, aes(PC1, PC2, color = batch, shape = condition)) +
  geom_point(size = 3) +
  ggrepel::geom_text_repel(aes(label = sample_id), size = 2.5) +
  coord_fixed(ratio = sqrt(pv[2] / pv[1])) +
  labs(x = sprintf("PC1 (%.1f%%)", 100 * pv[1]),
       y = sprintf("PC2 (%.1f%%)", 100 * pv[2])) +
  theme_bw()

The coord_fixed call matters more than it looks. If PC1 holds 60% of the variance and PC2 holds 8%, drawing them on equal-length axes visually inflates PC2 by a factor of about 2.7 and invites you to interpret scatter that is essentially noise. Setting the ratio to sqrt(pv2/pv1) renders distances in the plot proportional to distances in the data.

Read the plot in this order. Look at the scree values first: if PC1 is 70% and PC2 is 5%, you have one dominant axis and the second dimension is decoration. Look for samples sitting far from every other point, which are usually technical failures rather than interesting biology. Then ask whether any metadata variable separates the clusters, and only then ask whether your variable of interest does. A PCA plot is a map of between-sample distances, and its main use is showing you which experimental factor dominates those distances 1. Absolute coordinates and the sign of each component are arbitrary: prcomp can flip PC1 when you add a sample, so never compare signs across runs.

5. Find out what each component means

This is the step most tutorials omit, and it is where the plot stops being decorative. Regress each principal component on each metadata column and report the R²:

vars <- c("batch", "rin", "lib_size", "draw_date", "condition", "neutrophil_pct")
pc_meta_r2 <- sapply(vars, function(v) {
  sapply(1:6, function(k) {
    fit <- lm(pca$x[, k] ~ df[[v]])
    summary(fit)$r.squared
  })
})
round(pc_meta_r2, 2)   # rows PC1..PC6, columns = metadata

Any cell above roughly 0.5 is a component you can name. In whole blood, PC1 very often tracks cell-type composition (neutrophil fraction) rather than anything you manipulated, and in degraded samples PC1 tracks RIN. Follow up by inspecting the loadings, the per-gene weights that define the component:

top1 <- sort(pca$rotation[, 1])
head(names(top1), 25); tail(names(top1), 25)

If the extremes of PC1 are hemoglobin subunits, you are looking at globin carryover. If they are mitochondrial or histone genes, suspect degradation or a cell-stress artifact. Running a gene set enrichment on the loading vector turns an anonymous axis into an interpretable one, and pcaExplorer provides exactly this workflow interactively, pairing the sample scatter with gene loadings and functional annotation of the components 4. PCA is one member of a broader family of matrix factorizations whose latent factors can be read as biological programs, and treating the loadings as the object of interest rather than a byproduct is the productive stance 5.

6. Decide what to do about batch

If a nuisance variable explains a leading component, you have two options and they are not interchangeable. For the model, put the variable in the design (design = ~ batch + condition) and let DESeq2 account for it during testing. For the picture only, subtract it:

library(limma)
mat_adj <- removeBatchEffect(assay(vsd), batch = vsd$batch,
                             design = model.matrix(~ vsd$condition))

Passing design preserves the effect you care about while removing the batch shift. Plot the adjusted matrix to confirm the batch separation collapses, and label the figure as batch-adjusted. Never feed a removeBatchEffect output into DESeq() or any count-based test; the output is no longer counts, and the standard errors will be wrong because the adjustment is treated as known rather than estimated. If batch is perfectly confounded with your condition, no adjustment rescues it and the PCA is telling you the experiment cannot answer the question.

7. Know when PCA is the wrong tool

PCA finds orthogonal directions of maximum variance under a linear model, so it fails when the structure you care about is nonlinear, when it lives in a low-variance direction, or when a few outliers dominate the covariance matrix. Three neighboring questions come up constantly.

MDS versus PCA: limma::plotMDS and edgeR compute, by default, a distance between each pair of samples using the top 500 genes that distinguish that specific pair, then embed those distances. Because the gene set differs per pair, an MDS plot is not a projection and its axes have no loadings, though it is often more sensitive to structure carried by different gene sets in different sample groups. Classical MDS on Euclidean distances is mathematically equivalent to PCA; the limma default is not classical MDS. Run both. Agreement is reassuring, disagreement tells you the variable structure is heterogeneous.

Clustering: PCA does not cluster. It rotates the coordinate system, and you perceive clusters by eye. If you want cluster assignments, run k-means or hierarchical clustering on the first few PCs, which both denoises and speeds things up, but choose the number of components from the scree plot rather than defaulting to two. k-means partitions samples into groups; PCA produces continuous coordinates. They answer different questions and are frequently used together.

PCA versus PCR: principal component regression uses the leading PCs as predictors in a regression. It is a modeling method, unrelated to the polymerase chain reaction, and unrelated to what you are doing here.

For nonlinear structure, UMAP or t-SNE on the top 20 PCs will show more, at the cost of meaningless global distances and axes you cannot interpret. For multi-omic tables (transcripts plus proteins plus metabolites), single-table PCA on a concatenated matrix is a poor default because the block with the most features dominates, and joint factorization methods such as MOFA, DIABLO, or multi-block PLS are designed for the problem 67.

8. For longitudinal data from one person, decenter the individual

When every sample comes from the same person across months, PC1 is no longer between-person biology, so the leading component usually captures whichever technical or physiological axis moves most: sequencing batch, draw time, acute inflammation, or cell composition. Two practices help. First, plot the trajectory rather than a cloud: connect points in date order with geom_path and color by draw date, which turns drift into a visible line and makes a single aberrant draw obvious. Second, if you also have samples from other people, center within person (subtract each person’s gene means) before the PCA so the components describe within-person change rather than between-person differences. Quality control pipelines built for repeated-measures omics apply exactly this logic, running PCA-style projections at the sample level as a screen for outlying runs before any statistics 8.

Common problems

One sample sits far from all others on PC1. Check its library size, duplication rate, and RIN before touching biology. Recompute the PCA with that sample dropped: if the remaining structure reorganizes completely, the outlier was absorbing the covariance and everything you read from the first plot was about that one library.

PC1 is 90% of the variance and separates two obvious groups. Usually a batch, a protocol change, or a tissue difference. Confirm with the metadata regression in step 5 before celebrating.

The plot changes when you change ntop. Expected, and informative. Report the value you used, and check stability at 500, 2000, and all genes.

Replicates that should be identical are spread widely. Look at the axis percentages first. If PC2 holds 4%, the spread may be nothing at all, which is why the fixed aspect ratio matters.

PCA on TPM or FPKM instead of VST. This works less well than it appears. Both are ratio measures sensitive to a few dominant genes, and neither stabilizes variance across the expression range. If TPM is all you have, use log2(TPM + 1) and filter to genes with reasonable expression, and be aware you are looking at a compositional quantity 2.

Interpreting a component as a biological finding without checking loadings. A component is a direction, not a result. Name it from its loadings and its metadata associations, or leave it unnamed.

Finally, none of this is a clinical readout. A PCA plot of your own transcriptomes shows how your samples relate to each other, not whether anything is wrong with you. Any question about health status, symptoms, or what to do next belongs with a physician who can see the full picture.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Keunhong Son, Sungryul Yu, Wonseok Shin, et al. A Simple Guideline to Assess the Characteristics of RNA-Seq Data. BioMed Research International, 2018. https://doi.org/10.1155/2018/2906292 ↩ ↩2

  2. Thomas P Quinn, Ionas Erb, Greg Gloor, et al. A field guide for the compositional analysis of any-omics data. GigaScience, 2019. https://doi.org/10.1093/gigascience/giz107 ↩ ↩2

  3. Dhivyaa Rajasundaram, Joachim Selbig. More effort — more results: recent advances in integrative ‘omics’ data analysis. Current Opinion in Plant Biology, 2016. https://doi.org/10.1016/j.pbi.2015.12.010 ↩

  4. Federico Marini, Harald Binder. pcaExplorer: an R/Bioconductor package for interacting with RNA-seq principal components. BMC Bioinformatics, 2019. https://doi.org/10.1186/s12859-019-2879-1 ↩

  5. Genevieve L. Stein-O’Brien, Raman Arora, Aedin C. Culhane, et al. Enter the Matrix: Factorization Uncovers Knowledge from Omics. Trends in Genetics, 2018. https://doi.org/10.1016/j.tig.2018.07.003 ↩

  6. Irene Sui Lan Zeng, Thomas Lumley. Review of Statistical Learning Methods in Integrated Omics Studies (An Integrated Information Science). Bioinformatics and Biology Insights, 2018. https://doi.org/10.1177/1177932218759292 ↩

  7. Aurelia Morabito, Giulia De Simone, Roberta Pastorelli, et al. Algorithms and tools for data-driven omics integration to achieve multilayer biological insights: a narrative review. Journal of Translational Medicine, 2025. https://doi.org/10.1186/s12967-025-06446-x ↩

  8. David J. Degnan, Kelly G. Stratton, Rachel Richardson, et al. pmartR 2.0: A Quality Control, Visualization, and Statistics Pipeline for Multiple Omics Datatypes. Journal of Proteome Research, 2023. https://doi.org/10.1021/acs.jproteome.2c00610 ↩