Linkage disequilibrium and effective size

Drift builds association between nearby sites, recombination breaks it, and the balance tells you the population size.

Drift makes LD, recombination removes it

Two biallelic sites on one chromosome are in linkage equilibrium when their alleles combine at random. The module starts a Wright-Fisher population of N diploids in exact equilibrium. Each generation, 2N gametes are drawn from the parents. That draw is a sample, so by chance some allele pairs end up together more often than the allele frequencies predict. The difference is linkage disequilibrium (LD), measured here as r² = D²/(pA qA pB qB).

Recombination works against this. Sites c apart (the recombination fraction) are separated in a fraction c of gametes. Tightly linked sites keep their LD for many generations, and loosely linked ones lose it almost at once. The population settles at a balance set by ρ = 4Nc. Sved (1971) gave the result: r² ≈ 1/(1 + 4Nc).

Take N = 200 and two sites 1 cM apart, so c = 0.01 and 4Nc = 8. Then r² ≈ 1/9 = 0.11. For unlinked sites Hill (1981) gives r² ≈ 1/(3N), which is 0.0017. A sample of n chromosomes adds roughly 1/n on top of both numbers, 0.005 for the 200 chromosomes sampled at the defaults. The model has no mutation, so alleles are also drifting to fixation, and the sites left polymorphic slowly thin out.

r² falls with distance

Figure 1. Haplotype r² for 16,110 pairs of the 180 sites that pass the 5% frequency filter, after 200 generations at N = 200, seed 2. The ringed dots are bin means: 0.262 at 0.3 cM, 0.121 at 1.0 cM, 0.061 at 2.2 cM. The solid curve is the exact expectation at generation 200, the dashed one is Sved's equilibrium, and the floor at 0.0050 is sampling noise.

Tight pairs take longer to settle

The balance is not instant. The exact recursion for the gamete pool is x(t) = (1 − c)² [(1 − 1/2N) x(t−1) + 1/2N] + c²/2N, starting from x = 0. Its fixed point is Sved's formula for small c and Hill's for c = ½. The approach takes about 1/(2c + 1/2N) generations, because the (1 − c)² factor shrinks old LD by 2c per generation.

For the tightest class in the module, 0.33 cM apart, that is about 110 generations. For a pair 50 cM apart (c near 0.32) it is about 2. This is why the run length matters. A short run has reached equilibrium for loose pairs and not for tight ones, and it is why estimates from tight pairs reflect the founding population, which had no LD, more than the present one.

LD builds up at three distances

Figure 2. Mean r² of the whole population through 200 generations for a tight class (0.33 cM, ρ = 2.7), a moderate one (5 cM, ρ = 38) and a loose one (50 cM, ρ = 253). Dashed lines are the recursion. At generation 200 the observed means are 0.261, 0.017 and 0.006, against 0.231, 0.026 and 0.005 expected.

Exercise

Stop early

Open the module, type the seed 2, set Gen to 20 and run. Compare the binned r² near 0.3 cM and the estimate of N with the default run.

What to expect. Tight pairs have not built their LD: the 0.3 cM bin is at 0.052, against 0.262 after 200 generations. Loose pairs are already at equilibrium, so the LD-Ne estimate is 213 against the true 200. The tight bins hold no information about N yet.

Open the walkthrough step →

Reading N from LD

At equilibrium E[r²] ≈ f(c)/N, where f(c) = [(1 − c)² + c²] / [2c(2 − c)] and f = 1/3 for unlinked sites. Invert that and N ≈ f(c)/r². LD-Ne does this over all pairs with c ≥ 0.2, at least about 26 cM apart, which equilibrate in a few generations. It works on unphased genotypes of S sampled individuals, as real data are, so it first removes the sampling floor: r²′ = [r² − 1/(S − 1)] (S − 1)/(S − 2). Then N̂ = Σ f(c) / Σ r²′ over the pairs.

At seed 2 the estimate is N̂ = 201, from 9,243 loose pairs, with an approximate 95% interval of 134 to 396. The true value is 200. The interval comes from a delete-one-block jackknife over 20 blocks of sites, computed on 1/N̂, which is why its upper end can be infinite. At seed 1 the same parameters give 285, with an interval of 195 to 525. A single run scatters widely.

The weight f(c) matters. Pairs 26 cM to 100 cM apart on one chromosome are loosely linked, not unlinked. Treating them as unlinked, with f = 1/3, gives 115 here, a bias of nearly half.

N from each distance bin

Figure 3. Each distance bin inverted for N with the recursion (filled dots) and with Sved's formula (hollow dots), against the true N = 200. The band beyond about 26 cM is the LD-Ne estimate of 201 with its interval of 134 to 396. Tight bins overestimate N, because they still remember a founding population with no LD.

Exercise

Shorten the map, widen the data

With the seed 2, set Map (M) to 0.2 and run. Then set it to 0.3. Then set it back to 1 and set Sites to 1000 and Sample to 200.

What to expect. At 0.2 M no pair has c ≥ 0.2, so there is no LD-Ne estimate. At 0.3 M there are 467 loose pairs, N̂ = 159, and the interval is 46 to infinity. With 1000 sites and 200 individuals on 1 M there are 100,602 loose pairs and N̂ = 217 with an interval of 188 to 256.

Open the walkthrough step →

Summary

Drift creates LD every generation and recombination removes a fraction c of it, so r² settles near 1/(1 + 4Nc).
For unlinked sites r² is about 1/(3N), which is the signal LD-Ne reads.
Tight pairs need about 1/(2c + 1/2N) generations to reach the balance. Loose pairs need a few.
N̂ = Σ f(c) / Σ r²′ uses pairs with c ≥ 0.2 and removes the sampling floor 1/(S − 1). Reading them as unlinked is biased low.
At the defaults and seed 2, N̂ = 201 (134 to 396) for a true N of 200, and seed 1 gives 285. One run is one draw.
More sites and more individuals narrow the interval. Longer runs do not, because there is no mutation to replace the variation drift removes.