Hardy–Weinberg, multiple alleles and estimating allele frequencies

The Hardy–Weinberg principle is the null model of population genetics. It states that, in a large population with random mating and no selection, mutation or migration, allele frequencies do not change from one generation to the next and genotype frequencies settle into the proportions p², 2pq and q² after a single generation of random mating. BIOL-201 used these proportions with a chi-square test as a quality check on genotype data. This lesson goes further in three directions: it extends the principle to loci with many alleles and to sex-linked loci; it shows how to estimate allele frequencies when some genotypes cannot be distinguished, using the expectation–maximization (EM) algorithm on the ABO blood groups; and it replaces the chi-square approximation with an exact test that remains valid for rare alleles and small samples.

Why the proportions arise

The principle follows from a simple observation. Under random mating, the two gene copies that unite to form a zygote are drawn independently from the population’s pool of gametes. If allele A makes up a fraction p of that pool, the chance that both copies are A is p × p, the chance that both are a is q × q, and the chance of one of each is 2pq, the factor of 2 counting the two ways the heterozygote can form. Because each genotype passes on its alleles in proportion to its genotype, the allele frequency in the next gamete pool is p² + pq = p(p + q) = p. Nothing changes.

The assumptions are worth listing because each is the subject of a later lesson. Selection (lesson 3) changes allele frequencies by differential survival and reproduction. Finite population size (lesson 4) causes drift. Mutation (lesson 5) introduces new alleles. Non-random mating and population structure (lesson 7) change genotype frequencies even when allele frequencies stay fixed. Migration (lesson 8) mixes populations with different frequencies. The power of the Hardy–Weinberg principle is that it isolates these forces: any departure from the expected proportions is a clue that one of them is at work, or that the data contain errors.

Multiple alleles

With k alleles at frequencies p_1, p_2, … , p_k, random union of gametes gives homozygotes A_iA_i at frequency p_i² and heterozygotes A_iA_j at frequency 2p_ip_j. There are k homozygote classes and k(k − 1)/2 heterozygote classes. For a microsatellite with 15 alleles, that is 15 homozygotes and 105 heterozygotes, so many genotype classes have expected counts far below one even in a sample of several hundred individuals. This is the main reason the chi-square approximation fails for highly polymorphic loci.

The total expected heterozygosity is 1 − Σp_i², the H_e of lesson 1. With many alleles at similar frequencies, H_e approaches 1, and nearly every individual is heterozygous. This is exactly why forensic laboratories use highly polymorphic STR loci: two unrelated people are very unlikely to share a genotype at any one locus, and the probabilities multiply across loci, as lesson 12 shows.

Sex-linked loci and the approach to equilibrium

For an X-linked locus in a species where females are XX and males XY, females carry two copies and males one. Males inherit their single X from their mothers, while females inherit one X from each parent. If allele frequencies differ between the sexes, say p_f in females and p_m in males, then in the next generation the male frequency equals the previous female frequency, and the female frequency is the average of the previous male and female frequencies:

The difference between the sexes halves and changes sign every generation, so the two frequencies oscillate as they converge on the overall frequency (2p_f + p_m)/3. Unlike an autosomal locus, an X-linked locus therefore needs several generations to reach equilibrium after a disturbance such as the mixing of two populations. At equilibrium, the frequency of a recessive X-linked trait in males equals q, while in females it is q². This is why red–green colour vision deficiency, with q of about 0.08 in many populations of European ancestry, affects about 8 per cent of men but well under 1 per cent of women.

Estimating allele frequencies when genotypes are hidden: the ABO system

Counting alleles is simple when every genotype can be identified. The ABO blood group is the classic case where it cannot. Alleles A and B are codominant to each other and both are dominant to O, so the six genotypes produce only four phenotypes: type A (AA or AO), type B (BB or BO), type AB and type O (OO). Given counts of the four phenotypes, how should the three allele frequencies be estimated?

The answer is maximum likelihood, and the most convenient way to find the maximum is the EM algorithm. The idea is to alternate two steps. In the expectation step, we use the current guess of allele frequencies to split each ambiguous phenotype into its expected genotype counts: of the type A people, a fraction p_A²/(p_A² + 2p_Ap_O) are expected to be AA and the rest AO. In the maximization step, we count alleles in the reconstructed genotypes, exactly as if they had been observed, to get new frequency estimates. Repeating the two steps increases the likelihood at every iteration and converges to the maximum likelihood estimates.

import numpy as np

# invented phenotype counts for a sample of 1,000 blood donors
nA, nB, nAB, nO = 421, 106, 37, 436
n = nA + nB + nAB + nO

pA, pB, pO = 1/3, 1/3, 1/3                     # starting guess
for it in range(1, 101):
    # E step: split ambiguous phenotypes into expected genotype counts
    nAA = nA * pA**2 / (pA**2 + 2*pA*pO); nAO = nA - nAA
    nBB = nB * pB**2 / (pB**2 + 2*pB*pO); nBO = nB - nBB
    # M step: count alleles
    new = np.array([2*nAA + nAO + nAB, 2*nBB + nBO + nAB, 2*nO + nAO + nBO]) / (2*n)
    if np.max(np.abs(new - [pA, pB, pO])) < 1e-10:
        break
    pA, pB, pO = new
print(f"converged after {it} iterations")
print(f"pA = {pA:.4f}, pB = {pB:.4f}, pO = {pO:.4f}")

# expected phenotype counts under these estimates, and a goodness-of-fit check
exp = n * np.array([pA**2 + 2*pA*pO, pB**2 + 2*pB*pO, 2*pA*pB, pO**2])
obs = np.array([nA, nB, nAB, nO])
chi2 = np.sum((obs - exp)**2 / exp)
print("expected A, B, AB, O:", np.round(exp, 1))
print(f"chi-square = {chi2:.2f} on 1 degree of freedom")

converged after 13 iterations
pA = 0.2640, pB = 0.0743, pO = 0.6617
expected A, B, AB, O: [419.1 103.9  39.2 437.8]
chi-square = 0.19 on 1 degree of freedom

The algorithm converges in a handful of iterations. The goodness-of-fit test has only one degree of freedom: four phenotype classes, minus one because the counts sum to n, minus two for the two independent allele frequencies estimated. A quick alternative estimator, Bernstein’s method, uses the square root of the O frequency as p_O, but it does not force the three estimates to sum to 1 and is less efficient. The EM approach generalizes far beyond blood groups: the same logic is used to infer haplotype frequencies from unphased genotypes (lesson 9), to estimate admixture proportions (lesson 8) and to call genotypes from low-coverage sequencing data.

Testing Hardy–Weinberg proportions exactly

The familiar chi-square test compares observed and expected genotype counts. It relies on a large-sample approximation that breaks down when any expected count is small, which is common for rare alleles or highly polymorphic loci. For a biallelic locus there is a better option: the exact test. Conditional on the observed number of each allele, the number of heterozygotes in a sample of n individuals has a known probability distribution under Hardy–Weinberg proportions, first written down by Howard Levene in 1949:

Here n_A and n_B are the allele counts. The p-value is the sum of the probabilities of all heterozygote counts that are as likely as, or less likely than, the observed count. Janis Wigginton and colleagues published a fast algorithm in 2005 that made this test standard for genome-wide SNP data. The code below implements it by enumerating every possible heterozygote count, which is fast enough for any realistic sample.

import numpy as np
from scipy.special import gammaln
from scipy.stats import chi2

def hwe_exact(n_AA, n_AB, n_BB):
    n = n_AA + n_AB + n_BB
    nA = 2*n_AA + n_AB; nB = 2*n - nA
    rare = min(nA, nB)
    het = np.arange(rare % 2, rare + 1, 2)      # heterozygote counts with the right parity
    hom_r = (rare - het) // 2; hom_c = n - het - hom_r
    logp = (gammaln(n+1) - gammaln(hom_r+1) - gammaln(het+1) - gammaln(hom_c+1)
            + het*np.log(2) + gammaln(nA+1) + gammaln(nB+1) - gammaln(2*n+1))
    probs = np.exp(logp)
    p_obs = probs[het == n_AB][0]
    return probs[probs <= p_obs * (1 + 1e-7)].sum()

def hwe_chi2(n_AA, n_AB, n_BB):
    n = n_AA + n_AB + n_BB; p = (2*n_AA + n_AB) / (2*n)
    exp = n * np.array([p**2, 2*p*(1 - p), (1 - p)**2])
    x = np.sum((np.array([n_AA, n_AB, n_BB]) - exp)**2 / exp)
    return chi2.sf(x, 1)

for g in [(57, 36, 7), (88, 10, 2), (95, 3, 2), (40, 20, 40)]:
    print(f"AA, AB, BB = {g}: exact p = {hwe_exact(*g):.3g}, chi-square p = {hwe_chi2(*g):.3g}")

AA, AB, BB = (57, 36, 7): exact p = 0.789, chi-square p = 0.689
AA, AB, BB = (88, 10, 2): exact p = 0.0654, chi-square p = 0.0204
AA, AB, BB = (95, 3, 2): exact p = 0.00265, chi-square p = 2.72e-08
AA, AB, BB = (40, 20, 40): exact p = 7.9e-10, chi-square p = 1.97e-09

For the first genotype set, a common allele in good agreement with expectation, the two tests broadly agree. In the second and third sets, the rare allele appears in only 14 and 7 gene copies. In the second set the chi-square test would declare a significant departure at the 5 per cent level when the exact test does not, and in the third it exaggerates the evidence against Hardy–Weinberg proportions by a large factor because the expected number of BB homozygotes is tiny. The exact test gives a reliable answer. The last set shows a strong heterozygote deficit, which both tests detect easily.

Interpreting departures

A significant departure from Hardy–Weinberg proportions is a symptom, not a diagnosis. A deficit of heterozygotes can result from inbreeding, from population structure (the Wahlund effect of lesson 7), from assortative mating, from selection against heterozygotes or, very often, from technical problems such as null alleles that fail to amplify or from low-coverage sequencing that misses one allele. An excess of heterozygotes can result from selection favouring heterozygotes, from sampling a recent mixture of two populations, from very small numbers of breeders (when allele frequencies happen to differ between the sexes), or from paralogous sequences that are mistakenly merged so that two loci look like one heterozygous locus.

A useful summary of the direction and size of a departure is the fixation index F = 1 − H_o/H_e. Positive F means fewer heterozygotes than expected; negative F means more. Under inbreeding, genotype frequencies become p² + Fpq, 2pq(1 − F) and q² + Fpq. Because inbreeding affects every locus in the genome while technical artefacts and selection affect particular loci, a consistent positive F across many loci points to the mating system or population structure, whereas an extreme F at a handful of loci points to the loci themselves.

In genome-wide studies, each of hundreds of thousands of SNPs is tested, so a p-value threshold of 0.05 would flag thousands of loci by chance alone. Quality-control pipelines therefore use very stringent thresholds, such as p below 10⁻⁶ in controls, and test only in a group that should be in equilibrium. In case–control studies, a departure confined to cases can itself be a sign of association with disease, so filtering on cases would remove real signals.

For a plant population that reproduces partly by selfing at rate S, the equilibrium inbreeding coefficient is F = S/(2 − S). An F of 0.375 implies S of about 0.55, so roughly half of seeds would come from self-fertilization. Estimates of this kind are routine in plant breeding and in conservation genetics of rare plants.