The coalescent and neutral diversity statistics

Classical population genetics looks forward in time: given allele frequencies now, what will they be in the next generation? With DNA sequence data, it is often more natural to look backward. A sample of sequences has a genealogy, a tree describing how their lineages merge into common ancestors as we trace them into the past. Coalescent theory, introduced by John Kingman in 1982 and developed by Richard Hudson, Fumio Tajima and others, describes the shape of that genealogy under neutral evolution. It is the mathematical foundation of modern population genomics: it explains why the site frequency spectrum of lesson 1 is proportional to 1/i, provides two estimators of the population mutation rate, and gives a test of neutrality, Tajima’s D, that compares them. It also makes simulation of genomic data fast enough to fit complex demographic models.

Two lineages

Take two gene copies sampled from a Wright–Fisher population of 2N copies. In the previous generation, each copy had a parent chosen at random from the 2N copies, so the chance that both came from the same parent, that their lineages coalesce, is 1/(2N). If they do not coalesce, the same reasoning applies one generation further back. The number of generations until coalescence therefore follows a geometric distribution with mean 2N. When N is large it is well approximated by an exponential distribution with mean 2N generations. It is conventional to measure coalescent time in units of 2N generations (or 2N_e when the population is not ideal), so that the expected time for two lineages to coalesce is 1.

This simple result already explains the relationship between diversity and the mutation rate. Two sequences differ at every site where a mutation occurred on either branch since their common ancestor. The total branch length connecting them is twice the coalescence time, on average 4N generations, so with μ neutral mutations per site per generation, the expected number of differences per site is 4Nμ = θ. Nucleotide diversity, π, is therefore an estimator of θ, as lesson 5 assumed.

Many lineages

With k lineages, any of the k(k − 1)/2 pairs could coalesce. In coalescent time units, each pair coalesces at rate 1, so the waiting time until the first coalescence among k lineages is exponential with rate k(k − 1)/2 and mean 2/[k(k − 1)]. After that event, k − 1 lineages remain, and the process repeats until a single lineage, the most recent common ancestor (MRCA) of the sample, is reached. Two features follow immediately. The last stage, when only two lineages remain, is the longest, with expected duration 1. And adding more sequences to a sample mostly adds short branches near the tips:

in units of 2N generations. The expected time to the MRCA of a sample of 10 is already 1.8 units, 90 per cent of the value for an infinitely large sample. Sampling more individuals barely changes how far back the genealogy reaches. This is why even a small sample of genomes contains information about the deep history of a population, and why large samples are mainly useful for studying rare, recent variation.

Mutations on the genealogy and the site frequency spectrum

Under neutrality, mutations do not affect the genealogy; they simply fall on its branches at a constant rate, like raindrops on a tree. A mutation on a branch is inherited by every sampled sequence below that branch. Under the infinite-sites model, each mutation produces one segregating site, and its derived allele count is the number of descendants of the branch it fell on. This view makes the site frequency spectrum easy to understand. Mutations on the external branches leading to single tips produce singletons; mutations on branches subtending two tips produce doubletons; and so on. Fu showed in 1995 that the expected total length of branches with i descendants is proportional to 1/i, which is the 1/i spectrum of lesson 1.

Demographic history changes the shape of the genealogy and therefore the spectrum. In a population that has grown rapidly, coalescence rates were low in the recent past, when the population was large, and high further back, when it was small. The genealogy is star-like, with long external branches and short internal ones, producing an excess of singletons. A bottleneck or a population subdivided into isolated groups does the opposite: lineages coalesce quickly within groups or during the bottleneck, but the last few lineages take a long time to merge, giving long internal branches and an excess of intermediate-frequency variants.

Two estimators of θ

The coalescent gives two estimators of θ that use different aspects of the data. The first is π itself: the average number of pairwise differences, which has expectation θ (per site or per locus). The second, proposed by G. A. Watterson in 1975, uses the number of segregating sites S. Because the expected total branch length is 2Σ1/i in units of 2N generations, and mutations accumulate at rate θ/2 per unit of branch length, the expected number of segregating sites is θ times a_n = Σ1/i, summed from i = 1 to n − 1:

Both estimators have expectation θ under the standard neutral model with constant size, but they weight variants differently. Watterson’s estimator counts every segregating site equally, so it is sensitive to rare variants, which make up most segregating sites. π weights each site by the number of pairwise differences it creates, c(n − c) for a derived allele in c of n sequences, so intermediate-frequency variants contribute most.

Tajima’s D

In 1989 Fumio Tajima proposed comparing the two estimators. His statistic D is the difference π − θ_W divided by an estimate of its standard deviation under the neutral model:

Here e_1 and e_2 are constants that depend only on the sample size. Under the standard neutral model, D has a mean near zero and lies roughly between −2 and 2 most of the time. Negative D means an excess of rare variants, as produced by population growth, by a recent selective sweep (which wipes out variation and leaves new, rare mutations accumulating afterwards) or by purifying selection removing deleterious variants before they reach intermediate frequency. Positive D means an excess of intermediate-frequency variants, as produced by balancing selection, by a population bottleneck or by population structure.

The table shows the central difficulty in using D. Demography and selection can produce the same signal. The way out is to compare loci: demographic history affects the whole genome, while selection affects particular regions. A genome-wide average D of −0.5 in a population that has expanded, with a handful of genes showing D below −2, is more informative than any single value. Human populations outside Africa generally show more negative genome-wide D than many African populations because of the combined effects of the out-of-Africa bottleneck and subsequent rapid growth.

import numpy as np
rng = np.random.default_rng(322)

def tajima_D(counts, n):
    counts = np.asarray(counts); S = len(counts)
    i = np.arange(1, n)
    a1, a2 = np.sum(1/i), np.sum(1/i**2)
    b1 = (n + 1) / (3*(n - 1)); b2 = 2*(n*n + n + 3) / (9*n*(n - 1))
    c1 = b1 - 1/a1; c2 = b2 - (n + 2)/(a1*n) + a2/a1**2
    e1 = c1/a1; e2 = c2/(a1**2 + a2)
    pi = np.sum(counts*(n - counts)) / (n*(n - 1)/2)
    return pi, S/a1, (pi - S/a1) / np.sqrt(e1*S + e2*S*(S - 1))

def coalescent(n, theta):
    # returns derived-allele counts for each mutation on one neutral genealogy
    lineages = [frozenset([k]) for k in range(n)]
    counts = []
    while len(lineages) > 1:
        k = len(lineages)
        t = rng.exponential(2 / (k*(k - 1)))      # time in units of 2N generations
        for lin in lineages:                       # mutations at rate theta/2 per branch
            counts += [len(lin)] * rng.poisson(theta/2 * t)
        a, b = sorted(rng.choice(k, 2, replace=False), reverse=True)
        merged = lineages.pop(a) | lineages.pop(b)
        lineages.append(merged)
    return counts

n, theta = 10, 5.0
print("worked example:  pi = %.2f, theta_W = %.2f, D = %.2f" %
      tajima_D([1]*6 + [2]*3 + [3]*2 + [5, 7, 9], n))

D_vals, S_vals, pi_vals = [], [], []
for _ in range(4000):
    c = coalescent(n, theta)
    if len(c) > 1:
        p, w, d = tajima_D(c, n); D_vals.append(d); S_vals.append(len(c)); pi_vals.append(p)
print(f"neutral simulations: mean S = {np.mean(S_vals):.2f} (expected {theta*np.sum(1/np.arange(1,n)):.2f}), "
      f"mean pi = {np.mean(pi_vals):.2f} (expected {theta:.2f})")
print(f"mean D = {np.mean(D_vals):.2f}; 2.5% and 97.5% quantiles = {np.percentile(D_vals, 2.5):.2f}, {np.percentile(D_vals, 97.5):.2f}")

worked example:  pi = 4.42, theta_W = 4.95, D = -0.49
neutral simulations: mean S = 14.31 (expected 14.14), mean pi = 5.06 (expected 5.00)
mean D = -0.06; 2.5% and 97.5% quantiles = -1.74, 1.70

The worked example gives a D well inside the range of neutral simulations, so the slight excess of rare variants is not significant. The simulations also confirm the coalescent expectations for S and π. Simulating a null distribution like this, ideally under a demographic model fitted to genome-wide data, is how modern analyses judge whether a particular locus is unusual. Specialized software such as msprime simulates whole chromosomes with recombination for thousands of genomes in seconds, using the same principles as this short function.

Recombination and the ancestral recombination graph

So far each locus has had a single genealogy. Recombination means that different segments of a chromosome can have different genealogies, because, going back in time, a recombination event splits a lineage into two parents. The full history of a recombining chromosome is called the ancestral recombination graph. Nearby sites share almost the same tree; distant sites have nearly independent trees. This has two practical consequences. First, each independent stretch of the genome provides a new sample of the coalescent process, so a single diploid genome contains information about many genealogies. Second, because coalescence times vary along the genome, the density of heterozygous sites in one genome can be used to reconstruct how N_e changed over time.

The pairwise sequentially Markovian coalescent (PSMC) method of Heng Li and Richard Durbin, published in 2011, does exactly this. In regions where the two copies of a chromosome in one person coalesced long ago, heterozygous sites are dense; where they coalesced recently, sites are sparse. Converting the distribution of coalescence times into a history of N_e shows, for example, a reduction in the effective size of non-African human populations roughly 50,000 to 100,000 years ago, consistent with the out-of-Africa bottleneck. PSMC has been applied to hundreds of species, from mammoths to birds, and newer methods use many genomes at once to resolve more recent history. Converting coalescent time units into years requires assumed values for the mutation rate and generation time, so absolute dates carry considerable uncertainty.