Population structure with PCA
Principal components turn allele-frequency differences between populations into axes, as long as the differences are large enough to see.
Read first: Genetic drift
Each person is a point
A genotype table has one row per individual and one column per SNP, each cell 0, 1 or 2 copies of the alternate allele. Principal component analysis (PCA) finds the directions in which the rows differ most. The first direction is PC1, the next PC2, and so on. Each individual gets a score on each PC, and plotting PC1 against PC2 puts every person on a map.
The module simulates that table. K populations split from one ancestor and drift apart by the same amount F, which is also their expected pairwise FST. For each SNP the ancestor has frequency p₀ between 0.05 and 0.95, and population k gets a frequency with mean p₀ and variance F p₀(1 − p₀). Genotypes are drawn under Hardy-Weinberg from those frequencies. All SNPs are unlinked.
Before the PCA, each SNP is centred and scaled: x = (g − 2p̂) / √(p̂(1 − p̂)), where p̂ is the sample frequency. This is Patterson scaling. It gives every SNP the same weight regardless of its frequency. PC scores are then the eigenvectors of the individual-by-individual matrix XXᵀ/m, with m the number of SNPs.
Three populations, two structure axes
K populations give K − 1 axes
Three populations form a triangle, and a triangle is flat: two axes describe it. In general K populations need K − 1 PCs. With K = 2 and no admixture the same design gives one structure axis, which explains 6.3% of the variance at 30 per population. Every other PC is noise.
Noise PCs are not zero. With n individuals and m SNPs, the noise eigenvalues spread over a band, the Marchenko-Pastur bulk, whose upper edge sits at σ²(1 + √y)² with y = n/m. Here y = 0.106 and the edge holds 1.57% of the variance. A structure axis shows up only if its eigenvalue clears the edge. The cutoff is the Baik-Ben Arous-Peche (BBP) transition.
For K equal populations the cutoff is F/(1 − F) = (K/2)/√(nm). Patterson's rule for two populations is FST ≈ 1/√(nm), which is 0.0031 for this sample. The module computes the cutoff for the actual design, admixed individuals included: 0.0048. Below that value PCA cannot see the populations, however real they are.
Two spikes above the bulk
Exercise
Go below the threshold
Open the module with the defaults (176 individuals, 2000 SNPs), type the seed 1, set FST to 0.001 and Admix % to 0, and run. Read PC1 in the summary strip and the scree annotations.
What to expect. No spike is detectable. PC1 explains 1.06% of the variance, under the bulk edge of 1.08%, and the three clusters dissolve. The BBP cutoff for this design is 0.0027, above the true 0.001. The ancestry RMSE is 0.505, against 0.069 at the defaults.
Open the walkthrough step →Distance reads FST, and mixtures fill the triangle
Above the threshold, distances between population centroids carry information. McVean (2009) showed that for small F the squared distance between two population centroids is about 8 FST + 4/s, where s is the number of individuals per population. With F = 0.05 and s = 30 that is 0.40 + 0.13 = 0.53, a distance of 0.73. The run gives 0.72. The 4/s term is sampling noise in the centroids, and it shrinks as you add people.
An admixed individual inherits a share q from one population and 1 − q from another. Its allele frequency is the same mix, so its PC score lands on the line between the two population centroids, at the same fraction. In the map, the 16 admixed individuals fall between Pop 1 and Pop 2. Admixture proportions come from a second method, an ADMIXTURE-style maximum-likelihood fit started from the PCA, and the module compares it with the true proportions.
Structure that is smooth in space looks different. In the isolation-by-distance model, demes sit on a line and exchange migrants with their neighbours at rate m. Nearby demes share ancestry and distant ones do not, so the genotype matrix varies continuously along the line. The leading PCs of such a matrix are waves along the line: PC1 is a single gradient and PC2 a cosine that bends it into an arch.
A line of demes makes an arch
Exercise
Tighten and loosen the line
Choose IBD · line in the module (the preset) and type the seed 1. Run once, then set m to 0.01 and Gen to 300 and run, then set m to 0.5 with Gen back at 100 and run.
What to expect. The preset gives PC1 6.8%, PC2 3.4%, FST 0.068 and |r(PC1, x)| = 0.98. With m = 0.01 and 300 generations the demes diverge: PC1 14.3%, PC2 9.9%, FST 0.32. With m = 0.5 migration keeps them mixed: PC1 3.3%, PC2 1.8%, FST 0.020, and two spikes remain.
Open the walkthrough step →Summary
- PCA scores each individual on axes of genotype variation after Patterson scaling; K populations need K − 1 axes.
- Noise fills a Marchenko-Pastur bulk with upper edge σ²(1 + √y)², y = n/m. An axis counts only if it clears that edge.
- The detection cutoff for K equal populations is F/(1 − F) = (K/2)/√(nm). Below it the populations are invisible to PCA.
- Squared centroid distance is about 8 FST + 4/s, so distance on the map reads FST.
- Admixed individuals fall between the centroids of their source populations, in proportion to their ancestry.
- Isolation by distance gives gradients and an arch in PC2. Neither one is evidence of admixture.