Chapter 4 Principal components analysis
4.1 Motivation
We use principal components analysis (PCA) to clean up and compact the log-normalized expression matrix. Consider each gene as a dimension of our dataset where the cells are the observations, i.e., each cell’s expression profile defines its location in the high-dimensional expression space. PCA discovers axes in this high-dimensional space that capture the largest amount of variation (Pearson 1901). Each principal component (PC) corresponds to an axis in this space, where the earliest PCs capture the dominant factors of heterogeneity in our data. The idea is to use the first few PCs to approximate our original dataset9. Similarly, the Euclidean distances between cells in the PC space are approximations of the corresponding distances in the original dataset. PCA effectively compresses multiple genes into a much smaller number of dimensions (a.k.a., “eigengenes”, from Langfelder and Horvath (2007)), which allows us to use a much smaller matrix in downstream steps like clustering.
4.2 Getting the top PCs
Our assumption is that biological processes affect multiple genes in a coordinated manner. This means that the earlier PCs are likely to represent biological structure as more variation can be captured by considering the correlated behavior of many genes. In contrast, random technical or biological noise is expected to affect each gene independently. There is unlikely to be an axis that can capture random variation across many genes, meaning that noise should mostly be concentrated in the later PCs. By retaining the earlier PCs, we can focus on the biological signal while removing random noise. To demonstrate, we’ll pull out our favorite mouse brain dataset from Zeisel et al. (2015):
library(scRNAseq)
sce.zeisel <- ZeiselBrainData()
is.mito.zeisel <- rowData(sce.zeisel)$featureType=="mito"
# Performing some QC to set up the dataset prior to normalization.
library(scrapper)
sce.qc.zeisel <- quickRnaQc.se(sce.zeisel, subsets=list(MT=is.mito.zeisel), altexp.proportions="ERCC")
sce.qc.zeisel <- sce.qc.zeisel[,sce.qc.zeisel$keep]
# Computing log-normalized expression values.
sce.norm.zeisel <- normalizeRnaCounts.se(sce.qc.zeisel, size.factors=sce.qc.zeisel$sum)
# Computing the variances.
sce.var.zeisel <- chooseRnaHvgs.se(sce.norm.zeisel, more.var.args=list(use.min.width=TRUE))We run the PCA on our HVG-filtered log-normalized expression matrix to compute the top 25 PCs. This returns a matrix of “PC scores”, i.e., the coordinates for each cell in the new low-dimensional space. As discussed in Chapter 3, we restrict this step to the top HVGs to reduce the impact of random noise. While PCA is robust to noise, too much of it may cause the earlier PCs to ignore meaningful structure (Johnstone and Lu 2009).
## Mode FALSE TRUE
## logical 16006 4000sce.pca.zeisel <- runPca.se(sce.var.zeisel, features=hvgs.zeisel, number=25)
dim(reducedDim(sce.pca.zeisel, "PCA"))## [1] 2866 25We then use the PC score matrix in place of the original log-expression matrix in downstream steps like clustering (Chapter 6) and visualization (Chapter 5). The former is much smaller than the latter as the relevant heterogeneity is compacted into a smaller number of dimensions. This enables more efficient calculation of distances between cells and reduces the impact of random noise.
# This performs graph-based clustering on the PC score matrix. See the linked
# chapters for more in-depth discussion of this step.
table(clusterGraph.se(sce.pca.zeisel, reddim.type="PCA")$clusters)##
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 288 256 204 396 189 320 184 176 198 118 220 87 155 754.3 How many PCs?
The million dollar question is, how many of the top PCs should we retain for downstream analyses? Using more PCs will retain more biological signal at the cost of including more unintersting noise that might mask that signal. As with the number of HVGs, it is hard to determine whether an “optimal” choice exists here. Sure, technical variation is almost always uninteresting, but there is no straightforward way to automatically determine which aspects of biological variation are relevant to a particular scientific question. For example, heterogeneity within a population might be interesting when studying continuous processes like metabolic flux or differentiation potential, but could be considered noise in applications that only aim to distinguish between distinct cell types.
Most practitioners will simply set \(d\) to a “reasonable” but arbitrary value, typically ranging from 10 to 50. This is often satisfactory as the later PCs explain so little variance that their inclusion or omission has no major effect. For example, in the Zeisel dataset, very few PCs explain more than 1% of the variance in the entire dataset (Figure 4.1). Choosing between, say, 20 and 40 PCs would not even amount to 5 percentage points’ worth of difference in variance. In fact, the main consequence of using more PCs is simply an increase in compute time as more dimensions need to be processed, but most PC-related calculations are fast enough that this is not a practical concern.
sce.more.zeisel <- runPca.se(sce.var.zeisel, features=hvgs.zeisel, number=50)
pca.meta <- metadata(sce.more.zeisel)$PCA
percent.var <- pca.meta$variance.explained / pca.meta$total.variance * 100
plot(percent.var, log="y", xlab="PC", ylab="Variance explained (%)", type="b")
Figure 4.1: Percentage of variance explained by successive PCs in the Zeisel dataset, shown on a log-scale.
If we really must try to guess the “best” number of PCs10, here are a few approaches:
- We can choose the elbow point in the scree plot (Figure 4.1), e.g., using the
findElbowPoint()function from the PCAtools package. The assumption is that there should be a sharp drop in the percentage of variance explained when we move past the last PC corresponding to biological structure. However, the ideal cut-off can be difficult to gauge when there are sources of weaker biological variation. - We can keep the number of PCs that cumulatively explain variance equal to the sum of the biological components among the HVGs. This relies on the decomposition of each gene’s variance into biological and technical components (see Chapter 3). In practice, the separation between biological and technical variation is usually not so clear-cut; some technical noise will still contribute to the earlier PCs while some biological variation will still be present in the later PCs.
- We can use random matrix theory to select an appropriate number of PCs. This might involve the Marchenko-Pastur limit (Shekhar et al. 2016), Horn’s parallel analysis (Horn 1965), or the Gavish-Donoho threshold for optimal reconstruction (Gavish and Donoho 2014) (see relevant functions in PCAtools). Each of these methods has its own limitations, e.g., requirement for i.i.d. noise.
But if we’re really concerned about the number of PCs, it’s probably just better to repeat the analysis with different number of PCs. This allows us explore other perspectives of the data at different trade-offs between biological signal and technical noise.
4.4 Blocking on uninteresting factors
Larger datasets typically consist of multiple blocks of cells with uninteresting differences between them, e.g.. batch effects, variability between donors. We don’t want to waste our top PCs on capturing these differences - instead, we want our PCA to focus on the biological structure within each block. To demonstrate, let’s look at a dataset consisting of two plates of wild-type and oncogene-induced 416B cells (Lun et al. 2017). Differences in expression due to the plate of origin are obviously technical and should be ignored. To make life more exciting, we will also consider the oncogene induction status to be an uninteresting experimental factor11 that should not be allowed to dominate the PCA.
library(scRNAseq)
sce.416b <- LunSpikeInData("416b")
# Combining the plate of origin and oncogene induction status into a single
# blocking factor of 'uninteresting' variation.
plate.416b <- sce.416b$block
pheno.416b <- ifelse(sce.416b$phenotype == "wild type phenotype", "WT", "induced")
sce.416b$block <- factor(paste0(pheno.416b, "-", plate.416b))
# Computing the QC metrics.
library(scrapper)
is.mito.416b <- which(any(seqnames(rowRanges(sce.416b)) == "MT"))
sce.qc.416b <- quickRnaQc.se(
sce.416b,
subsets=list(MT=is.mito.416b),
altexp.proportions="ERCC",
block=sce.416b$block
)
sce.qc.416b <- sce.qc.416b[,sce.qc.416b$keep]
# Computing log-normalized expression values.
sce.norm.416b <- normalizeRnaCounts.se(
sce.qc.416b,
size.factors=sce.qc.416b$sum,
block=sce.qc.416b$block
)
# Choosing the top HVGs after blocking on the uninteresting factors.
sce.var.416b <- chooseRnaHvgs.se(
sce.norm.416b,
more.choose.args=list(top=1000), # just picking a cool-looking number of top genes here.
block=sce.norm.416b$block
)We specify block= to instruct runPca.se() to focus on the variation within each block.
Briefly, this involves centering each block at the origin, finding the axes of largest variation among the residuals,
and then projecting the expression values for each cell onto these axes to obtain that its PC scores.
By examining the variance of the residuals, we reduce the contribution of uninteresting differences between blocks;
conversely, we provide more opportunities for the top PCs to capture biological variation within each block.
For example, blocking removes some of the shift between the induced and wild-type subpopulations on the first two PCs (Figure 4.2).
is.hvg.416b <- rowData(sce.var.416b)$hvg
sce.pca.416b <- runPca.se(sce.var.416b, features=is.hvg.416b, number=20)
sce.block.416b <- runPca.se(
sce.var.416b,
features=is.hvg.416b,
number=20,
block=sce.var.416b$block
)
library(scater)
gridExtra::grid.arrange(
plotReducedDim(sce.pca.416b, dimred="PCA", colour_by="block") +
ggtitle("Without blocking"),
plotReducedDim(sce.block.416b, dimred="PCA", colour_by="block") +
ggtitle("Blocked"),
ncol=2
)
Figure 4.2: First two PCs for the 416B dataset, before and after blocking on uninteresting experimental factors. Each point represents a cell, colored by its combination of experimental factors.
It’s worth mentioning that blocking in runPca.se() is not guaranteed to remove the differences between blocks.
If we’re lucky, blocking alone is sufficient if all of the inter-block differences are orthogonal to the major biological variation,
such that taking the first few PCs will focus on the latter and remove the former.
In practice, inter-block differences tend to have some biological component that will be preserved within the first few PCs.
For example, oncogene induction should have some effect on proliferation decisions in 416B cells,
so if there is any variation in proliferative pathway activity within each block,
the differences between blocks would still be captured in the top PCs.
If such differences are uninteresting, their removal requires some additional effort (see Chapter 8) prior to downstream steps like clustering.
Nonetheless, blocking is still helpful as it eliminates at least some of these differences and preserves more biological signal in the top PCs.
For completeness’ sake: technically speaking, we could attempt to remove the differences between blocks explicitly by setting center.scores.by.block=TRUE.
This will center the PC scores within each block at the origin, and is effectively equivalent to deriving the PC scores from the residuals.
(Careful readers will note that the default behavior of block= is to compute only the axes of variation from the residuals, not the scores themselves.)
Centering each block removes any systematic difference in location between blocks but is only correct in very limiting circumstances,
e.g., assuming all blocks have the same subpopulation composition and the difference between blocks is consistent for all cell subpopulations.
Such assumptions may be appropriate in some situations - say, when technical replicates are involved - but are not generally applicable.
In our 416B example, the two subpopulations are now forced together (Figure 4.3) for better or worse.
sce.resid.416b <- runPca.se(
sce.var.416b,
features=is.hvg.416b,
number=20,
block=sce.var.416b$block,
more.pca.args=list(center.scores.by.block=TRUE)
)
library(scater)
plotReducedDim(sce.resid.416b, dimred="PCA", colour_by="block")
Figure 4.3: First two PCs for the 416B dataset with blocking, where PC scores are computed from residuals. Each point represents a cell, colored by its combination of experimental factors.
As with HVGs, we should only use block= for experimental factors that are not interesting.
If we were interested in the effects of oncogene induction, we should not block on it to ensure that the PCA can capture the associated changes in expression.
Sometimes, though, it is not obvious whether something is “interesting” or not,
as we may wish to ignore some biological differences to obtain a consistent set of clusters across treatment conditions, tissues, etc.
Check out Chapter 8 for a more detailed discussion.
4.5 Visualizing the PCs
We might as well touch on another common use of PCA, which is visualization of high-dimensional data. This is used in a variety of fields and applications including bulk RNA-seq, but is not as effective for scRNA-seq data. If we’re lucky, our population structure is simple enough that the first two PCs capture most of the relevant biology (Figures 4.2 and 4.3). However, in most cases, relevant biological heterogeneity is spread throughout 10-50 PCs that are much harder to visualize. For example, examination of the top 4 PCs is still insufficient to resolve all subpopulations identified by Zeisel et al. (2015) (Figure 4.4).
library(scater)
plotReducedDim(
sce.pca.zeisel,
dimred="PCA",
ncomponents=4,
colour_by="level1class"
)
Figure 4.4: PCA plot of the first 4 PCs in the Zeisel brain data. Each point is a cell, coloured according to the annotation provided by the original authors.
The problem here is that PCA is a linear technique, i.e., only variation along a line in high-dimensional space is captured by each PC. As such, it cannot efficiently represent high-dimensional differences in the first 2 PCs. If the first PC is devoted to resolving the biggest difference between subpopulations, and the second PC is devoted to resolving the next biggest difference, then the remaining differences will not be visible in the plot. That said, PCA is still useful as the top PCs are often used as input to more sophisticated algorithms for dimensionality reduction (Chapter 5).
Session information
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.4 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] scater_1.41.2 ggplot2_4.0.3
## [3] scuttle_1.23.1 ensembldb_2.37.3
## [5] AnnotationFilter_1.37.0 GenomicFeatures_1.65.0
## [7] AnnotationDbi_1.75.2 scrapper_1.7.3
## [9] scRNAseq_2.27.0 SingleCellExperiment_1.35.2
## [11] SummarizedExperiment_1.43.0 Biobase_2.73.2
## [13] GenomicRanges_1.65.1 Seqinfo_1.3.0
## [15] IRanges_2.47.2 S4Vectors_0.51.6
## [17] BiocGenerics_0.59.10 generics_0.1.4
## [19] MatrixGenerics_1.25.0 matrixStats_1.5.0
## [21] BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 jsonlite_2.0.0 magrittr_2.0.5
## [4] ggbeeswarm_0.7.3 gypsum_1.9.0 farver_2.1.2
## [7] rmarkdown_2.31 BiocIO_1.23.3 vctrs_0.7.3
## [10] memoise_2.0.1 Rsamtools_2.29.0 RCurl_1.98-1.19
## [13] htmltools_0.5.9 S4Arrays_1.13.0 BiocBaseUtils_1.15.1
## [16] AnnotationHub_4.3.2 curl_7.1.0 BiocNeighbors_2.7.2
## [19] Rhdf5lib_2.1.0 SparseArray_1.13.2 rhdf5_2.57.3
## [22] sass_0.4.10 alabaster.base_1.13.1 bslib_0.11.0
## [25] alabaster.sce_1.13.0 httr2_1.3.0 cachem_1.1.0
## [28] GenomicAlignments_1.49.1 lifecycle_1.0.5 pkgconfig_2.0.3
## [31] rsvd_1.0.5 Matrix_1.7-6 R6_2.6.1
## [34] fastmap_1.2.0 digest_0.6.39 irlba_2.3.7
## [37] ExperimentHub_3.3.1 RSQLite_3.53.3 beachmat_2.29.0
## [40] labeling_0.4.3 filelock_1.0.3 httr_1.4.8
## [43] abind_1.4-8 compiler_4.6.1 bit64_4.8.2
## [46] withr_3.0.3 S7_0.2.2 BiocParallel_1.47.0
## [49] viridis_0.6.5 DBI_1.3.0 HDF5Array_1.41.0
## [52] alabaster.ranges_1.13.0 alabaster.schemas_1.13.0 rappdirs_0.3.4
## [55] DelayedArray_0.39.3 rjson_0.2.23 tools_4.6.1
## [58] vipor_0.4.7 otel_0.2.0 beeswarm_0.4.0
## [61] glue_1.8.1 h5mread_1.5.0 restfulr_0.0.17
## [64] rhdf5filters_1.25.2 grid_4.6.1 gtable_0.3.6
## [67] BiocSingular_1.29.0 ScaledMatrix_1.21.0 XVector_0.53.0
## [70] ggrepel_0.9.8 BiocVersion_3.24.0 pillar_1.11.1
## [73] dplyr_1.2.1 BiocFileCache_3.3.0 lattice_0.22-9
## [76] rtracklayer_1.73.0 bit_4.6.0 tidyselect_1.2.1
## [79] Biostrings_2.81.6 knitr_1.51 gridExtra_2.3.1
## [82] bookdown_0.47 ProtGenerics_1.45.0 xfun_0.60
## [85] UCSC.utils_1.9.0 lazyeval_0.2.3 yaml_2.3.12
## [88] evaluate_1.0.5 codetools_0.2-20 cigarillo_1.3.1
## [91] tibble_3.3.1 alabaster.matrix_1.13.0 BiocManager_1.30.27
## [94] cli_3.6.6 jquerylib_0.1.4 dichromat_2.0-1
## [97] Rcpp_1.1.2 GenomeInfoDb_1.49.1 dbplyr_2.6.0
## [100] png_0.1-9 XML_3.99-0.23 parallel_4.6.1
## [103] blob_1.3.0 bitops_1.1-0 viridisLite_0.4.3
## [106] alabaster.se_1.13.0 scales_1.4.0 purrr_1.2.2
## [109] crayon_1.5.3 rlang_1.3.0 cowplot_1.2.0
## [112] KEGGREST_1.53.6References
Gavish, M., and D. L. Donoho. 2014. “The optimal hard threshold for singular values is \(4/\sqrt {3}\).” IEEE Trans. Inf. Theory 60 (8): 5040–53.
Horn, J. L. 1965. “A Rationale and Test for the Number of Factors in Factor Analysis.” Psychometrika 30 (2): 179–85.
Johnstone, I. M., and A. Y. Lu. 2009. “On consistency and sparsity for principal components analysis in high dimensions.” J. Am. Stat. Assoc. 104 (486): 682–93.
Langfelder, P., and S. Horvath. 2007. “Eigengene networks for studying the relationships between co-expression modules.” BMC Syst. Biol. 1 (November): 54.
Lun, A. T. L., F. J. Calero-Nieto, L. Haim-Vilmovsky, B. Gottgens, and J. C. Marioni. 2017. “Assessing the reliability of spike-in normalization for analyses of single-cell RNA sequencing data.” Genome Res. 27 (11): 1795–1806.
Pearson, K. 1901. “On lines and planes of closest fit to systems of points in space.” Lond. Edinb. Dubl. Phil. Mag. 2 (11): 559–72.
Shekhar, K., S. W. Lapan, I. E. Whitney, N. M. Tran, E. Z. Macosko, M. Kowalczyk, X. Adiconis, et al. 2016. “Comprehensive classification of retinal bipolar neurons by single-cell transcriptomics.” Cell 166 (5): 1308–23.
Zeisel, A., A. B. Munoz-Manchado, S. Codeluppi, P. Lonnerberg, G. La Manno, A. Jureus, S. Marques, et al. 2015. “Brain structure. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq.” Science 347 (6226): 1138–42.
And in fact, PCA does provide the optimal low-rank approximation, based on its relation to the SVD and the Eckart-Young-Mirsky theorem. This is arguably one of the rare cases of mathematical rigor in this entire book.↩︎
Probably because a reviewer asked us to.↩︎
Which is not entirely unreasonable. Say that we want to identify matching cell states across the wild-type and induced populations. In such cases, we would want to ignore the induction effect so that the matching states will cluster together.↩︎