Population Structure and Eigenanalysis

Picture the old maps of human variation that came out of the early days of genome-wide data. Principal components analysis gave us those axes where Europe ran left to right, Africa down to up, and populations fell into clouds. Useful and vivid, but a little hand-wavy. Patterson, Price, and Reich took that map and bolted a measuring instrument onto it. They transformed PCA from a pretty picture into a testable framework, one that can tell you not just where the axes point, but whether an axis is real structure or just noise. They uncovered a sharp threshold in detectability. For a fixed amount of data, tiny divergences are practically invisible; nudged just above a line, they pop right out. Before we dive into the statistics, here’s the practical promise. This approach is fast. It scales to hundreds of thousands of markers and thousands of people without needing exotic hardware. In their benchmarks, multiplying a one hundred person by one hundred thousand marker matrix by its transpose took about four seconds on a good workstation. A bigger run, two thousand people and half a million markers, was estimated at about two and a half hours on a single processor. It’s also deliberately agnostic about crisp categories. Individuals, not populations, are the rows of the matrix, and the goal isn’t to force everyone into boxes. It’s to find the major axes along which the data vary and ask, formally, which axes mean something. That formal part rides on a random-matrix view of genotype data. Take a big rectangular matrix M: rows are individuals, columns are markers, and each entry is a normalized genotype count. Now form X by multiplying M by its transpose. X is a sample covariance—think of it as how each person correlates with every other across the genome. If there’s no population structure, this X looks like a classic Wishart matrix. Its eigenvalues, the numbers that tell you how much variance lies along each principal axis, should all sit near a noise band. If there is real structure, the top few eigenvalues get big. Johnstone showed that, after a precise centering and scaling, the largest eigenvalue under the null follows the Tracy–Widom distribution, a universal limit that shows up for many kinds of random matrices. Crucially for genetics, you don’t need the genotype entries to be Gaussian for this to hold. By results of Soshnikov on universality, as long as columns are independent and nicely behaved, you still get Tracy–Widom in the limit. What does that buy you? A simple test. Compute the eigenvalues of X. Take the largest one, apply Johnstone’s centering and scaling—those are recipe-like steps that turn the raw eigenvalue into a standardized statistic—and read off a p-value from the Tracy–Widom curve. If it’s small, you’ve got a significant axis of structure. If it’s not, you stop. If it is, you keep going: test the second eigenvalue, then the third, in sequence. After each significant axis, you adjust the effective dimensions of the Wishart model and repeat. Johnstone’s theory states that this sequential testing is conservative for an honest-to-goodness Wishart, and in genetic data, the tail of the distribution—the part that tells you about significance—fits well. Inside that engine is a biological model that links what we care about to what those eigenvalues are doing. Imagine K populations descending from a common ancestor. Let P be the ancestral allele frequency at a site, and p the vector of current population frequencies. Conditional on P, the mean is P along every coordinate, and the covariance has a simple form: P times one minus P, multiplied by a matrix B that captures how much the populations have diverged. If B has rank r—there are r independent directions along which populations have split—then r of the eigenvalues in the theoretical covariance stand out and grow with sample size. The rest don’t. That’s exactly the world Tracy–Widom is designed for: a handful of big signals and a sea of small ones. Here’s where the phase change comes in. In the spiked-covariance picture, there’s a neat threshold first analyzed by Baik, Ben Arous, and Péché. Let the ratio of markers to individuals be n over m, and define c so that c squared equals n divided by m. There’s a number, one plus one over c, that acts like a wall. If the leading population-structure eigenvalue of the true covariance is below that wall, the top sample eigenvalue behaves just like noise and the Tracy–Widom statistic stays modest. Cross that wall, and the leading eigenvalue peels away from the pack; the Tracy–Widom statistic shoots up, and detection becomes easy. Translated into population genetics, that wall is a detectability threshold in terms of the fixation index, F S T, and data size. Patterson and colleagues show, for example, that with a data budget on the order of two to the twentieth genotype entries, the phase change falls around an FST of about 0.0025. Below that, good luck. Just above, you’re in business. An interesting corollary is that above the threshold, if you hold the total data size fixed, increasing the number of individuals tends to grow the Tracy–Widom statistic, so more sampling pays off. None of this matters if the inputs are messy, so they make the preprocessing concrete. Start with the raw genotype matrix C—individuals by markers—with entries zero, one, or two counting the alternate allele. Center each column by subtracting its mean, and scale it so that, under Hardy–Weinberg equilibrium and no structure, every marker has comparable variance. That gives you M, the matrix you’ll actually analyze. This normalization stabilizes variance across markers and makes drift-driven changes comparable. For microsatellites, which behave differently, a simpler mean subtraction can be used. The test itself needs one more ingredient: the effective number of markers. Dense arrays are loaded with linkage disequilibrium, or LD—nearby markers travel together—so the number of effectively independent columns is smaller than the raw count. Rather than trying to specify that a priori, they estimate it from the data, using the moments of the eigenvalues. In practice, they rescale the eigenvalues so they sum to the effective sample size on the individual side, then solve a simple equations system to back out an effective marker count. That plugs into the Tracy–Widom calibration and keeps the test honest when LD is present. Linkage disequilibrium is more than a nuisance—it can wreck eigenvectors. A long block in high linkage disequilibrium can dominate a principal component and make a single region look like global structure. The fix they propose is refreshingly direct. For each marker, predict it from a small window of previous markers, then analyze the residuals instead of the raw genotypes. In other words, strip out the local, predictable part, and feed what’s left to PCA. In simulations, even when they dialed linkage disequilibrium all the way up—blocks in perfect linkage—the residualization pulled the Tracy–Widom statistic back into the expected range. With no linkage disequilibrium, the correction barely moves the needle. The practical advice is to LD-prune aggressively—drop one of each pair in very tight linkage disequilibrium—and, if you’re working with dense panels, residualize within short windows before testing. Quality control extends beyond linkage disequilibrium. Missing data can masquerade as structure, particularly if certain genotypes are systematically called as missing in some labs or batches. A neat diagnostic here is to build a matrix of missingness—zero where a call exists, one where it doesn’t—and run the same eigenanalysis. If that matrix shows structure, your genotype PCA likely will too, for the wrong reason. Closely related individuals are another trap: family structure produces its own principal components. The fix is familiar—remove close relatives up front. And check the inputs carefully. Markers with high error rates or strange call distributions can inflate apparent structure. How well does all this work on real data? On the Human Genome Diversity Project’s African samples, the top eigenvalue translated to a Tracy–Widom statistic around 46, with a p-value below one in a trillion. The second axis also came in decisively, around 6.7 with a p-value below one in a trillion, and the third failed—its p-value hovered around 0.74—so two axes pass and then you stop. In a Northern Thai, Chinese, and Japanese trio, the separation was even crisper: a top Tracy–Widom statistic near 92 with a p-value below one in a trillion, a second around 31 with the same tiny p-value, and no significant third axis. Shriver’s multi-population dataset told a similar story, with leading Tracy–Widom statistics in the 70 to 100 range and vanishingly small p-values for the first few axes, then a slide into the noise. It’s the pattern the theory predicts: a small number of big, interpretable directions, then a long tail of nothing much. Interpretation matters. Those large axes map onto demographic events—splits and major clines—more than they do onto sampling quirks. A lovely illustration comes from admixture. Suppose populations A, B, and D diverged with modest differentiation—think an FST on the order of 0.005—and a fourth group, C, formed recently as an A–B admixture. In their simulations, the admixture fractions in C followed a skewed Beta distribution, averaging about 70 percent from A. The result? Two significant eigenvalues, not three. C didn’t create a new axis; its individuals lie along the line between A and B in PCA space. In real data, when admixture is recent and the panels are informative, the first principal component can correlate almost perfectly with ancestry proportions—correlations around 0.995 aren’t unusual. The catch is admixture linkage disequilibrium. Long ancestry blocks inflate many eigenvalues at once and can push the Tracy–Widom calibration around. In an African-American dataset with a highly informative ancestry panel and very strong differentiation between source populations—FST around 0.58—the Tracy–Widom theory struggled. That’s a sign to reach for admixture-aware tools like STRUCTURE, ADMIXMAP, or ANCESTRYMAP if your goal is local ancestry or fine-grained blocks. For global ancestry and a quick read on structure, PCA remains a strong default. There’s a quiet but important contrast here with heuristic rules you might have seen, like the broken-stick model that compares each eigenvalue to what you’d expect by slicing a stick at random. Broken-stick doesn’t care how many markers you typed, so on huge datasets it can miss real structure. The Tracy–Widom approach incorporates the effective marker count into the calibration and gives you an actual p-value. The testing is clean: if an eigenvalue isn’t significant, you don’t keep looking. You stop. How does this sit alongside STRUCTURE, the popular model-based clustering tool? They view the methods as complements. STRUCTURE tries to assign individuals probabilistically to K clusters and to estimate K. PCA, or eigenanalysis in their language, outputs positions along axes and tests for the presence of structure. When the structure is obvious, both methods tend to agree. When it’s marginal, neither has an easy job. A pragmatic link between them is this: use the number of significant principal components as a guide to K in STRUCTURE—often, one more than the count of significant axes is a sensible starting point. Because PCA is so fast—SMARTPCA is the implementation they provide—you can iterate quickly: check for structure, adjust for it in an association study, or plan a deeper model-based run only where it’s warranted. There’s one more idea to carry forward, and it’s the organizing principle behind the whole framework. Detectability isn’t smooth. It has a ledge. Below a certain combination of fixation index and data size, you won’t see structure no matter how clever your algorithm. Step just over it, and it becomes easy to detect, with p-values plunging into the floor. That mental model helps in three concrete ways. It tells you how many people and markers you need to see a given level of divergence. It explains why sampling design—who you include and in what proportions—can make clines hard or easy to spot without changing the underlying axes. It keeps you honest about the limits of what PCA can and cannot see, especially in the face of linkage disequilibrium and recent admixture. If you’re heading into a new dataset, the playbook looks like this. Clean the inputs and thin tight linkage disequilibrium. Center and normalize by marker so drift is on the same scale across the genome. Consider a short-window residualization if your panel is dense. Estimate the effective number of markers from the eigenvalue moments. Then let Tracy–Widom call balls and strikes on the axes, one by one. Cross-check with simple supervised tests—analysis of variance on labels, if you have them—and with a missingness eigenanalysis to make sure your batch effects aren’t talking. If you need crisp clusters or local ancestry, bring STRUCTURE or related tools into the mix. If you want speed, a principled test for structure, and coordinates you can regress on, start with PCA. The big win here is intellectual economy. Cavalli-Sforza’s maps meet modern random-matrix theory, and out comes a framework that is both interpretable and fast. It gives you numbers where before you had rules of thumb. It tells you when an axis is real. It quantifies that moment when a weak signal crosses a threshold and turns into a discovery.