In the previous post, we saw how pooling populations can make well-genotyped variants fail Hardy–Weinberg equilibrium (HWE) tests. How can we account for that structure while still detecting departures that warrant investigation?
Kwong et al. address this with RUTH1, a method that incorporates ancestry and genotype uncertainty into HWE testing. RUTH uses principal components to predict individual allele frequencies, then tests genotype proportions against those expectations.
Principal component analysis
Imagine describing each person by their genotypes at thousands of variants. Each variant is a coordinate: 0, 1, or 2 copies of a chosen allele. A person becomes a point in a space with thousands of dimensions. Principal component analysis (PCA) finds a few directions through that cloud that capture its largest patterns of variation.
With just two measurements, you could draw the cloud on paper. Rotate a line through its center and project every point onto it. The direction with the greatest spread of projected points is PC1. PC2 captures the greatest remaining spread along a direction perpendicular to PC1. PCA applies the same idea in many dimensions.
The PCA algorithm
- Arrange the genotypes in a matrix with individuals as rows and variants as columns. Each entry gᵢⱼ counts copies of the allele chosen for variant j.
- Center each variant by subtracting its mean genotype, 2pⱼ, where pⱼ is its sample allele frequency. A common genetic PCA convention also divides by √[2pⱼ(1 − pⱼ)], its HWE standard deviation. This puts variants with different allele frequencies on a comparable scale. Using the HWE variance for scaling does not establish that a variant is in HWE.
- Find the directions of greatest variation using a singular value decomposition (SVD) of the standardized matrix, or an equivalent eigenvector calculation. The directions are ordered by how much variation they capture and are mutually perpendicular.
- Project each person onto those directions to obtain their PC scores: their positions along the new axes. Plotting PC1 against PC2 gives the familiar genetic PCA scatterplot.
Call the standardized matrix Z and the weight of variant j in component k, vⱼₖ. Each PC score is the sum of the standardized genotypes multiplied by these weights:
\[ \begin{aligned}Z_{ij} &= \frac{g_{ij}-2p_j}{\sqrt{2p_j(1-p_j)}} \\ \mathrm{PC}_{ik} &= \sum_j Z_{ij}v_{jk}\end{aligned} \]
Monomorphic variants have a zero denominator and are excluded. This describes a common PCA setup; score normalization and handling of missing genotypes vary across implementations.
Projection variance
Consider a direction through the center of the cloud. Represent it by a vector v: one weight per variant. For person i, the projection onto that direction is sᵢ = Σⱼ Zᵢⱼvⱼ. Here N is the number of individuals and Σ means to sum over the indicated index. Because each column of Z is centered, these scores also have mean zero. Their variance measures how widely the people spread out along the proposed axis.
\[ \begin{aligned}V(\mathbf v) &= \frac{\sum_i s_i^2}{N-1} \\ \mathbf v_1 &= \underset{\lVert\mathbf v\rVert=1}{\operatorname{argmax}}\; V(\mathbf v)\end{aligned} \]
Read the second line as “choose the unit-length direction v that gives the largest variance.” The length constraint matters: doubling every weight would double the scores without changing the direction, artificially quadrupling their variance.
What argmax means
Max gives the largest value; argmax gives the input that produces it. For the numbers [2, 5, 1], max is 5 and argmax selects the second position. In the PCA equation, the input is an entire direction vector, so argmax returns that vector rather than a position in a list.
Power iteration and eigenvectors
Argmax specifies what we want, but does not tell the computer how to find it. One way is power iteration. Start with a trial unit vector, project the people onto it, and use those projections to update the direction:
- Calculate the scores s = Zv.
- Calculate u = Zᵀs. The transpose Zᵀ turns the scores back into one number per variant: uⱼ = Σᵢ Zᵢⱼsᵢ. Each person’s vector is weighted by their signed projection.
- Normalize: replace v with u / ‖u‖, then repeat.
The update repeatedly multiplies the trial vector by ZᵀZ. An eigenvector is a direction that this multiplication stretches without turning; its eigenvalue is the stretch factor. Components along eigenvectors with larger eigenvalues grow faster, so repeated updates favor the direction with the largest eigenvalue.
The covariance matrix C = ZᵀZ / (N − 1) has the same eigenvectors. Each eigenvalue of C equals the variance of scores along its unit eigenvector. PC1 therefore has the largest eigenvalue. Each later component maximizes variance while remaining perpendicular to all earlier components.
Singular value decomposition computes these directions directly from Z without needing to form ZᵀZ. Power iteration is a way to understand the calculation; it is not a claim about which PCA implementation supplied RUTH’s inputs.
| Term | What it describes |
|---|---|
| Variant loading | A variant’s weight in a PC direction. |
| PC score | One person’s coordinate along that direction. |
| Explained variance | How much of the standardized matrix’s variation that component captures; its fraction is its eigenvalue divided by the sum of eigenvalues. |
PC1 does not inherently mean a particular ancestry, and its sign can be reversed without changing the result. In genetic data, ancestry-related allele-frequency differences often create strong patterns across many variants, so leading PCs can describe ancestry structure. Relatives, batch effects, and unusual genomic regions can also influence the axes.
For ancestry PCA, variant QC and LD pruning help prevent errors or many correlated variants in one region from dominating the result. See the PLINK PCA documentation.2
Principal components as predictors of individual allele frequencies
PCA summarizes patterns across many variants into a few coordinates per person. RUTH takes those coordinates as fixed inputs and fits an allele-frequency relationship for each variant it tests. RUTH’s fitted coefficients β are a separate part of its model from the PCA loadings v.
The PCA plot’s “percent variance explained” refers to variation across the genotype matrix. The variance in the Wahlund identity below refers to variation in allele probabilities at one variant. RUTH connects them by using PC scores to predict those probabilities.
Hardy–Weinberg expectations under population structure
Let p be the frequency of allele A and q = 1 − p. Under HWE, the genotype probabilities are p² for AA, 2pq for Aa, and q² for aa. Heterozygosity follows the concave parabola:
\[ H(p)=2p(1-p) \]
HT is the heterozygosity expected from the total pooled allele frequency. HS averages the separate HWE expectations, traditionally within subpopulations; here we average across individuals with their own allele probabilities pᵢ.
For samples with different allele frequencies, evaluating this curve at their average frequency gives a value at least as large as averaging their separate heights on the curve:
\[ H_T-H_S=2\operatorname{Var}(p_i) \]
In this identity and the HWE illustration below, variance uses denominator N and individuals have equal weight.
Individual-specific allele frequencies
At each variant, RUTH estimates an individual-specific allele frequency (ISAF), written p̂ᵢ, from person i’s PC scores. The hat indicates an estimate. β₀ is the intercept, and each coefficient βₖ describes how the predicted frequency changes with PC k. K is the number of PCs included. These coefficients are fitted separately for each variant.
\[ \hat p_i=\operatorname{clip}\!\left(\beta_0+\sum_{k=1}^{K}\beta_k\mathrm{PC}_{ik},\;\varepsilon,\;1-\varepsilon\right) \]
ε (epsilon) is a small positive minimum allele probability. The linear prediction can fall outside the range of valid probabilities, so clip replaces values below ε with ε and values above 1 − ε with 1 − ε. This also avoids probabilities of exactly zero or one. The paper used ε = 1/(4N), where N is the number of individuals, in its evaluation; we use the same choice below.
The estimated p̂i is the allele probability for someone at person i’s ancestry coordinates. Many samples inform this relationship. With one PC, β₀ = 0.5 and β₁ = 0.2 give predicted frequencies of 0.3, 0.5, and 0.7 at PC scores −1, 0, and 1, respectively.
Testing for HWE departures after ancestry adjustment
RUTH represents the remaining departure from HWE with θ (theta), which the paper calls an inbreeding coefficient. It measures the relative reduction in heterozygosity compared with the ancestry-adjusted HWE expectation. RUTH fits one θ per variant, shared across individuals.
θ = 0 means no departure from that expectation. Positive θ means a heterozygote deficit; negative θ means an excess. Despite its name, this coefficient does not identify the cause: inbreeding, genotyping errors, or population structure missing from the model can all produce a positive estimate. It is not a measure of an individual person’s inbreeding.
Genotype probabilities under the model
For positive θ, probability shifts from heterozygotes equally to the two homozygote types. A heterozygote carries one A and one a, and an equal mixture of AA and aa has the same allele balance. This increases homozygosity while preserving each individual’s modeled allele frequencies pᵢ and qᵢ = 1 − pᵢ.
The heterozygote probability decreases by 2θpᵢqᵢ, and each homozygote probability increases by θpᵢqᵢ:
| Genotype | HWE: θ = 0 | Allow a departure |
|---|---|---|
| aa | qᵢ² | qᵢ² + θpᵢqᵢ |
| Aa | 2pᵢqᵢ | 2pᵢqᵢ − 2θpᵢqᵢ |
| AA | pᵢ² | pᵢ² + θpᵢqᵢ |
For example, with p = q = 0.5 and θ = 0.2, heterozygote probability falls from 50% to 40%, while each homozygote probability rises from 25% to 30%. θ = 0.2 therefore means a 20% reduction relative to expected heterozygosity, a drop of 10 percentage points.
For negative θ, the transfer reverses. Its permitted range depends on the allele frequencies so that all genotype probabilities remain nonnegative. The illustration below focuses on positive θ. RUTH does not explicitly model family relationships or person-specific inbreeding.
The model’s average heterozygosity is Hmodel = (1 − θ)HS. Its gap from the pooled HWE expectation separates into two terms:
\[H_T-H_{\mathrm{model}}=2\operatorname{Var}(p_i)+\theta H_S\]The first term comes from variation in allele probabilities; the second is the additional departure from the ancestry-adjusted expectation. This identity describes expected heterozygosity; the tests below use genotype evidence to assess whether θ differs from zero.
One-PC illustration
21 equally weighted illustrative individuals sit at fixed coordinates from −1 to +1. Here β₀ = 0.5 and ε = 1/(4N). You set the slope and departure; real RUTH estimates parameters from genotype evidence. No genotypes are sampled or p-values calculated here.
Allele frequency predicted from PC1
Ancestry-adjusted heterozygosity
\[\begin{aligned}\bar p &= \frac1N\sum_i p_i \\ H_S &= \frac1N\sum_i 2p_i(1-p_i)\end{aligned}\]The gold square averages the green points’ horizontal and vertical coordinates. It falls below the curve when predicted frequencies differ. The orange point evaluates the curve at the mean frequency. The blue point shows (1 − θ)Hₛ; its distance below the gold square is the additional deficit θHₛ.
Try slope = 0: the ancestry-related gap vanishes. Increase θ: an additional deficit appears even after ancestry is accounted for. The curve shows the HWE baseline; θ lowers the blue point and the third bar while leaving the allele frequencies unchanged.
Genotype likelihoods and uncertainty
A hard call selects aa, Aa, or AA. Sequencing evidence may support several possibilities, particularly at low depth. A genotype likelihood Lᵢ(G) = P(dataᵢ | G) measures how well genotype G explains the reads. The posterior probability P(G | dataᵢ) instead describes support for a genotype after combining those reads with the model. Here G ranges over aa, Aa, and AA.
For one individual, RUTH combines the evidence with its modeled genotype probabilities:
\[ P(\mathrm{data}_i\mid p_i,\theta)=\sum_G L_i(G)P(G\mid p_i,\theta) \]
Combining model probabilities with genotype evidence
This illustration fixes θ at 0. The evidence sliders set relative likelihoods for each genotype. Their values need not sum to one; only their ratios matter for the posterior probabilities below.
| Genotype | Model probability | × Evidence | Contribution | Posterior |
|---|
The marginal likelihood is the sum of contributions. Dividing each contribution by that sum gives the posterior genotype probabilities shown here. Scaling all three likelihoods by the same constant leaves those posteriors unchanged.
Parameter estimation and significance testing
The null hypothesis is θ = 0: HWE after adjustment for the modeled ancestry. The alternative allows θ to differ from zero. The tests assess whether allowing that departure improves the explanation of the genotype evidence enough to exceed what sampling variation would plausibly produce.
Across individuals, RUTH adds the log of each marginal likelihood. It estimates the frequency coefficients using an expectation–maximization (EM) algorithm: alternate between updating support for the possible genotypes and updating parameters to improve the likelihood.
| Test | What it measures |
|---|---|
| Score | Fit the null model (θ = 0). Measure the likelihood’s slope in the θ direction, scaled by its information. |
| Likelihood ratio | Fit both θ = 0 and θ free. Re-estimate the frequency coefficients under each; compare the best log-likelihoods. |
Write the total log-likelihood as ℓ(β, θ). The likelihood-ratio statistic TLRT is twice the improvement from allowing θ to vary; max means the best log-likelihood found by fitting the indicated parameters:
\[ T_{\mathrm{LRT}}=2\left[\max_{\boldsymbol\beta,\theta}\ell(\boldsymbol\beta,\theta)-\max_{\boldsymbol\beta}\ell(\boldsymbol\beta,0)\right] \]
Both tests use an asymptotic χ² distribution with one degree of freedom for p-values. The slope control above holds the frequency relationship fixed for explanation; RUTH’s LRT fits it again under the alternative.
Results and limitations
Kwong et al. evaluated RUTH using simulations and chromosome 20 data from 1000 Genomes and TOPMed. Their comparisons showed benefits from accounting for ancestry and genotype uncertainty, especially at low sequencing depth.
The main caution is reduced power for rare variants with heterozygote deficits. Sparse genotype counts can make the chi-square approximation used by both tests unreliable. The score test can lose further power when θ is far from zero; the likelihood-ratio test avoids that particular disadvantage but retains the approximation limitation.
The authors also did not explicitly model family relationships or individual-level inbreeding, and did not evaluate admixture effects or imputed genotypes. A PC-based frequency model can leave some population structure unexplained.3
For QC, a significant result flags a departure but does not establish its cause. A nonsignificant result, especially at a rare variant, may reflect limited power. Interpret either alongside genotype counts, the direction and magnitude of departure, and other sample and variant QC evidence. Simulations matched to the dataset’s ancestry structure, sample size, allele frequencies, and genotype uncertainty can help assess false-positive rates under HWE and power to detect departures.
References
Kwong et al. (2021), Robust, flexible, and scalable tests for Hardy–Weinberg equilibrium across diverse ancestries. Model and tests: page 3; method overview: Figure 3; evaluation: Figures 1–4. RUTH software.↩︎
PLINK 2.0 documentation: population stratification and principal component analysis.↩︎
Kwong et al. (2021), Discussion, pages 11–12, for the authors’ discussion of power, asymptotic approximations, and modeling limitations.↩︎