Principal Component Analysis: An Introduction

Suppose you measure $p$ variables (SNP genotypes, gene expression levels, skull dimensions) on $n$ individuals. Each individual is a point in $p$-dimensional space. If the variables are correlated — which biological variables almost always are — the cloud of points is not spherical; it is elongated, like a cigar or a pancake, because variation along some combinations of variables is much larger than along others. PCA finds the orientation of that cigar. Formally, a principal component is a linear combination of the original (usually centered, and often standardized) variables,

chosen so that (1) the variance of $PC_1$ across individuals is as large as possible, and (2) the weight vector $\mathbf{a} = (a_1,\dots,a_p)$ has unit length, $\sum a_i^2 = 1$, so the solution is not made arbitrarily large by simply inflating the coefficients. The second component, $PC_2$, is defined the same way subject to the additional constraint that it be uncorrelated with (geometrically, orthogonal to) $PC_1$. This continues until there are as many PCs as original variables, at which point the total variance in the data has been fully re-partitioned — no information is lost, only reorganized, if you keep all $p$ components. The crucial theorem that makes PCA computable is that these variance-maximizing directions are exactly the eigenvectors of the covariance matrix $\mathbf{\Sigma}$ (or the sample covariance matrix $\mathbf{S}$), and the variance captured by each PC is exactly the corresponding eigenvalue $\lambda_i$. Because $\mathbf{S}$ is symmetric, its eigenvectors are automatically orthogonal, which is precisely the uncorrelated-axes property PCA requires — this is not a coincidence but the reason PCA is defined via eigen-decomposition rather than by some more complicated search. The eigenvalues sum to the total variance in the data (the trace of $\mathbf{S}$), so $\lambda_i / \sum_j \lambda_j$ is the proportion of total variance explained by $PC_i$ — the number quoted on every PCA plot axis label.

Step 1: Center (and possibly scale) the data. For each variable, subtract its mean across individuals so that every variable has mean zero. If variables are on very different scales (e.g., height in cm mixed with a hormone concentration in pg/mL), each is also divided by its standard deviation, converting the covariance matrix into a correlation matrix; otherwise PCA would be dominated by whichever variable happens to have large numeric units, an artifact of measurement scale rather than biology.

Step 2: Compute the covariance matrix. For $p$ centered variables measured on $n$ individuals, the sample covariance matrix is

a $p \times p$ symmetric matrix whose diagonal entries are variances of each variable and whose off-diagonal entries are covariances between pairs of variables. This matrix is the entire input to PCA — PCA never looks at the raw data again once $\mathbf{S}$ (or the correlation matrix) is formed.

Step 3: Eigen-decompose the covariance matrix. Solve $\mathbf{S}\mathbf{a} = \lambda\mathbf{a}$ for eigenvectors $\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_p$ and eigenvalues $\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_p \geq 0$. Each eigenvector, called a loading vector, gives the recipe (the weights $a_i$) for building one principal component out of the original variables. In genomic applications this decomposition is usually done via Singular Value Decomposition of the (centered) data matrix directly, which is numerically more stable than forming $\mathbf{S}$ explicitly and is mathematically equivalent.

Step 4: Project the data. Each individual's PC score is obtained by taking the dot product of their (centered) trait vector with each loading vector:

These scores are the new coordinates plotted on a PCA scatterplot; an individual far to the right on the PC1 axis has a trait combination strongly aligned with the PC1 loading vector.

Step 5: Decide how many PCs to keep. Because $\lambda_i$ measures variance captured, a scree plot of eigenvalue against PC number typically shows a steep drop followed by a long, shallow tail; the "elbow" is often used as an informal cutoff, keeping only the PCs before the tail flattens, on the reasoning that the tail represents noise rather than structured signal.

The key intuition from this figure generalizes without change to 500,000 dimensions: the covariance matrix defines an ellipsoid (in $p$ dimensions) whose principal axes are exactly the PCs, ordered from longest (most variance) to shortest.

In genomic datasets the "signal" PCs (here, PC1–PC3, together explaining 77% of variance) typically correspond to interpretable axes such as ancestry, sex, or a strong batch effect, while the long flat tail is diffuse noise from individual SNP genotyping error or expression measurement variability — this is exactly the noise floor that motivates dimensionality reduction before downstream statistical testing in Genome-Wide Association Statistics.

The canonical demonstration that PCA recovers real, interpretable biological structure — not just an arbitrary rotation of noise — is the population-genetic work initiated by Cavalli-Sforza's group in the 1970s–90s and brought to genome-wide data by Alkes Price et al. (2006, Nature Genetics, "Principal components analysis corrects for stratification in genome-wide association studies") and, most vividly, John Novembre et al. (2008, Nature, "Genes mirror geography within Europe").

What was done. Novembre and colleagues genotyped roughly 3,000 Europeans (the POPRES dataset) at about 500,000 SNPs across the genome. They built the individual-by-SNP genotype matrix (each entry coded as 0, 1, or 2 copies of a reference allele), computed the sample covariance matrix across individuals, and performed PCA — technically an SVD-based variant adapted for genotype data (the EIGENSOFT/smartpca method developed by Price, Patterson, and Reich). No geographic coordinates or country labels were given to the algorithm; only genotypes were used as input.

What was observed. The scatterplot of PC1 versus PC2 scores reproduced, with startling fidelity, a rotated and slightly distorted map of Europe: individuals from Portugal and Spain clustered in the southwest of the plot, Finns and other Scandinavians clustered in the northeast, and individuals from intermediate countries (France, Germany, Italy, the Balkans) filled in the intermediate space in roughly the correct relative geographic positions. The correlation between an individual's PC1–PC2 position and their true country-of-origin coordinates was strong enough that the authors could predict an individual's country of origin to within a few hundred kilometers, purely from genotype-derived PC scores.

What it proved. This result demonstrated that most human genetic variation within a continental group is structured by isolation by distance — genetic similarity decays smoothly with geographic distance because migration and mating are geographically local — and that this smooth, continuous gradient (rather than discrete "races" or clusters) is precisely the kind of variance-maximizing linear structure that PCA is built to detect. It also validated the earlier, more theoretical proposal of Price et al. (2006) that the top few PCs of genome-wide genotype data are an efficient empirical proxy for ancestry, directly justifying their use as covariates to correct for population stratification — a confounding structure that, left uncorrected, produces false-positive associations in genome-wide association studies (a problem foreshadowing the methodological core of Genome-Wide Association Statistics). Related work by Cavalli-Sforza, Menozzi and Piazza in The History and Geography of Human Genes (1994) had shown similar continent-scale gradients using classical PCA on allele-frequency data decades earlier, and Reich and colleagues later extended the same logic to ancient DNA, using PCA to project ancient samples onto axes built from present-day populations and thereby infer ancestry proportions and migration history.

Example 1 — population structure from a toy SNP dataset. Suppose you genotype 6 individuals at 4 SNPs (coded as allele dosage 0, 1, or 2), and — to keep the arithmetic tractable — the individuals come from two weakly diverged populations, A (individuals 1–3) and B (individuals 4–6):

Column means are (1.0, 1.0, 1.0, 0.833). Centering subtracts these from each column; e.g., individual 1's centered SNP1 value is $0 - 1.0 = -1.0$. Without grinding through the full $4\times4$ eigen-decomposition by hand, the qualitative structure is instructive: SNP1 and SNP2 both take low values in Pop A and high values in Pop B, so they are strongly positively correlated across individuals; SNP3 takes high values in Pop A and low values in Pop B, so it is strongly negatively correlated with SNP1 and SNP2. SNP4 shows no consistent population difference — it varies almost independently of population. A PCA on this covariance structure will produce a PC1 loading vector with large positive weight on SNP1 and SNP2, large negative weight on SNP3, and near-zero weight on SNP4 — because that combination is exactly the linear combination whose value differs most between the two clusters, hence has maximal variance across all six individuals. Projecting each individual onto this PC1 gives strongly negative scores for individuals 1–3 and strongly positive scores for individuals 4–6, cleanly separating the two populations on a single axis — this is the mechanism, in miniature, behind the Novembre et al. European map: real, replicated allele-frequency differences between geographically or ancestrally structured groups are precisely the kind of correlated variation that maximizes variance along a shared axis, so PCA discovers population structure as a side effect of maximizing variance, without ever being told population labels. If instead you standardized each SNP (dividing by its standard deviation) before running PCA, SNP4 — despite carrying almost no population signal — would receive a small but nonzero PC1 loading (it is not perfectly independent in a sample of only 6), and doing PCA on unstandardized allele counts is in fact the more common and defensible choice in population genetics, because SNP allele frequency and variance are mathematically linked (variance $=2p(1-p)$ under Hardy–Weinberg), and this relationship is precisely what many ancestry-informative PCA methods intentionally rely on.

Example 2 — a morphometric PCA with an asymmetric size/shape signal. A field biologist measures 5 skull traits (in mm) on a sample of turtle skulls spanning both sexes and a wide size range: skull length (SL), skull width (SW), jaw length (JL), orbit diameter (OD), and interorbital distance (IOD). The raw trait correlation matrix (correlation, not covariance, is used here because traits are on the same unit but on different natural scales of variability) looks like this:

Notice the asymmetry: SL, SW, JL, and IOD are all strongly positively correlated with one another (0.74–0.91), but OD (orbit diameter) is only weakly correlated with the others (0.22–0.34). This is a classic biological pattern — most skull dimensions scale together with overall body/head size, but eye size is comparatively size-invariant across individuals (a phenomenon well documented in vertebrate cranial allometry, since neural and sensory tissue often grows more slowly than the rest of the skull, a pattern called negative allometry). PCA on this correlation matrix typically yields:

PC1's loadings would be approximately (0.49, 0.47, 0.48, 0.20, 0.46) — nearly equal, large positive weights on SL, SW, JL, IOD, and a much smaller weight on OD — confirming PC1 is a size axis: an individual with a high PC1 score is simply a large-skulled turtle across almost every dimension measured. PC2's loadings would be dominated by OD with an opposite sign from the rest, roughly (−0.25, −0.22, −0.20, 0.90, −0.18) — a shape axis contrasting relatively large eyes against a comparatively small skull, capturing individuals whose orbit size deviates from what overall skull size alone would predict. A biologist plotting PC1 (x-axis) against PC2 (y-axis) would see the population spread mainly left-to-right by size, with residual up-down scatter reflecting the size-independent, orbit-driven variation — precisely the kind of size/shape decomposition that motivated Sewall Wright's early use of principal-axis methods on turtle shell measurements in the 1930s and that remains standard practice in modern geometric morphometrics.

Principal Component Analysis was invented independently twice: by Karl Pearson in 1901, in a paper titled "On Lines and Planes of Closest Fit to Systems of Points in Space," which posed PCA geometrically as finding the best-fitting line or plane through a cloud of points by minimizing the perpendicular (orthogonal) distances of points to that line — a formulation mathematically equivalent to maximizing variance along the line. The method was independently rediscovered and given its modern statistical formulation, name, and computational algorithm by Harold Hotelling in 1933, in "Analysis of a Complex of Statistical Variables into Principal Components," where he framed the problem explicitly in terms of extracting uncorrelated components ranked by variance from a set of correlated psychological test scores. Its biological uptake was slower: it entered population genetics substantially through Luigi Luca Cavalli-Sforza, whose decades-long project mapping human allele-frequency variation onto geography (culminating in the 1994 monograph The History and Geography of Human Genes, with Paolo Menozzi and Alberto Piazza) used principal-components maps of Europe and the world as a central visual and analytical tool decades before genome-wide SNP data existed, and it was modernized for genome-scale data by Nick Patterson, Alkes Price, and David Reich in the mid-2000s, whose EIGENSOFT software and stratification-correction methodology made PCA a mandatory preprocessing step in essentially every genome-wide association study conducted since.

Genome-wide association studies (GWAS). The top 5–20 PCs of genome-wide genotype data are routinely included as covariates in association-testing regressions to absorb population stratification, preventing ancestry-correlated allele frequency differences from masquerading as disease associations — a direct methodological bridge to Genome-Wide Association Statistics.

Single-cell and bulk transcriptomics. PCA is the near-universal first step in single-cell RNA-seq pipelines (e.g., Seurat, Scanpy), used to reduce tens of thousands of genes to 10–50 PCs before clustering cells into putative cell types or feeding the reduced space into nonlinear embeddings such as t-SNE or UMAP.

Ancient DNA and human population history. Ancient genomes are commonly "projected" onto PCs built from present-day reference populations, letting researchers infer the ancestry composition and migration history of archaeological samples relative to modern genetic variation (used extensively in Reich-lab ancient-DNA studies of Eurasian population turnover).

Conservation genetics. PCA on SNP or microsatellite data from endangered populations detects cryptic population structure, inbreeding-driven substructure, and admixture between subspecies, informing management units and translocation decisions.

Agricultural and plant breeding genomics. PCA of genotyped breeding-line panels corrects for population structure in genomic-selection and marker-association models in crops such as maize, rice, and wheat, analogous to its GWAS role in humans.

Quantitative morphology and paleontology. PCA of landmark or linear-measurement matrices remains the standard tool for summarizing shape variation across fossil specimens, living species comparisons, and developmental growth series, separating size-related from shape-related variation as in the turtle-skull example above.

Public health surveillance. PCA of pathogen genomic surveillance data (e.g., influenza or SARS-CoV-2 sequence variation) helps visualize antigenic drift and lineage relationships, informing vaccine strain selection.

Ask the tutor: How do I decide whether a cluster I see on a PCA plot of my RNA-seq or SNP data reflects real biology versus a technical batch effect?