A sample in the present

A sample of 10 individuals in the present, \(t_0\).

Time runs left to right, into the past.

One generation back

Each individual draws a parent uniformly at random from the previous generation, independently of every other individual.

This generational sampling assumption is the Wright-Fisher model. It would be appropriate for a monoecious annual plant population.

Two alleles that share a parent

\[P(\text{share an ancestor 1 generation ago}) = \frac{1}{2N}\]

The first allele picks some parent; the second picks the same one with probability \(1/2N\).

Two alleles that do not

\[P(\text{do not share an ancestor}) = 1 - \frac{1}{2N}\]

Following two lineages back

Waiting for the first coalescence

\(P(\text{two alleles coalesce in generation } i \mid N)\)

\[P(t=i) = \left( 1 - \frac{1}{2N} \right)^{i-1}\left( \frac{1}{2N} \right)\]

That is a geometric distribution

Geometric with \(p = 1/2N\):

\[P(t \mid N) = \left( 1 - \frac{1}{2N} \right)^{i-1}\left( \frac{1}{2N} \right)\]

\[\boxed{\;\mathbb{E}[t] = \frac{1}{p} = 2N\;}\]

Not coalescing for a while

\(P(\text{two alleles do not coalesce in } i \text{ generations})\)

\[P(t > i) = \left( 1 - \frac{1}{2N} \right)^{i}\]

An exponential approximation

\[P(t > i) = \left( 1 - \frac{1}{2N} \right)^{i}\]

Since \(\left( 1 - \frac{1}{2N} \right) \approx e^{-\frac{1}{2N}}\), \(\qquad P(t > i) \approx e^{-\frac{i}{2N}}\)

Rescaling time

\[P(t > i) \approx e^{-\frac{i}{2N}} \approx e^{-T}\]

Measuring time in units of \(2N\) generations removes \(N\) from the expression entirely.

Coalescent units

\[P(t > i) \approx e^{-T}\]

\[T = \frac{\text{number of generations}}{2N}\]

One coalescent unit \(=2N\) generations.

More than two alleles

Everything so far followed a pair. A real sample has \(n\) lineages, and any pair of them can coalesce.

One lineage

Two lineages

Three lineages

Three alleles, one generation

\(P(\text{3 alleles do not coalesce in 1 generation})\)

\[= \left( 1 - \frac{1}{2N} \right) \times \left( 1 - \frac{2}{2N} \right)\]

The second lineage must miss the first; the third must miss both.

Two alleles

\(P(\text{2 alleles do not coalesce in 1 generation})\)

\[= 1 - \frac{1}{2N}\]

Three alleles

\(P(\text{3 alleles do not coalesce in 1 generation})\)

\[= \left( 1 - \frac{1}{2N} \right) \times \left( 1 - \frac{2}{2N} \right)\]

\(n\) alleles

\(P(n \text{ alleles do not coalesce in 1 generation})\)

\[= \left( 1 - \frac{1}{2N} \right) \times \left( 1 - \frac{2}{2N} \right) \times \cdots \times \left( 1 - \frac{n-1}{2N} \right)\]

Dropping the small terms

• Expanding the product leaves terms that look like \(1/N\) terms that have \(1/N^2\).
• Because we expect \(n \ll N\), we can drop the terms that have \(1/N^2\) from the expansion.

\[1 - \sum_{i=1}^{n-1}\frac{i}{2N} + \sum_{i<j}\frac{i}{2N}\cdot\frac{j}{2N} - \cdots\]

\[\approx 1 - \frac{1 + 2 + \cdots + (n-1)}{2N}\]

The complement

\[P(n \text{ alleles do}) = 1 - P(n \text{ alleles do not})\]

\[= \frac{1 + 2 + \cdots + (n-1)}{2N}\]

Summing the series

Given: \[ 1 + 2 + \cdots + (n-1) = \frac{n(n-1)}{2}, \] \[P(n \text{ alleles do}) = \frac{n(n-1)}{4N}\]

Waiting time for \(n\) alleles

We can now develop an expression for the probability that \(n\) alleles coalesce in an arbitrary generation:

\[ P(n \text{ alleles coalesce in the } i\text{th generation}) \]

\[= \left[ 1 - \frac{n(n-1)}{4N} \right]^{i-1} \frac{n(n-1)}{4N}\]

Geometric again

This fits nicely into the expectation of a geometric distribution again. So, the expected waiting time is the reciprocal of the success probability:

\[\boxed{\;\mathbb{E}[t] = \frac{4N}{n(n-1)}\;}\]

The combinatoric form

\(\frac{n(n-1)}{2}\) is the number of pairs, so the same result can be written as:

\[= \left[ 1 - \binom{n}{2}\left( \frac{1}{2N} \right) \right]^{i-1} \binom{n}{2}\frac{1}{2N}\]

Each of the \(\binom{n}{2}\) pairs coalesces at rate \(1/2N\).

Expectation, combinatoric form

\[\mathbb{E}[t] = \frac{2N}{\binom{n}{2}}\]

Identical to \(4N/n(n-1)\) — the pair-counting view just makes the mechanism explicit.

To continuous time

As before, the geometric waiting time is well approximated by an exponential once \(N\) is large.

The per-pair rate

A single pair coalesces at rate

\[T = \frac{1}{2N}\]

per generation.

Counting the ways

There are \(n-1\) coalescent events to get from \(n\) lineages to one, and \(\binom{n}{2}\) ways for the next one to happen.

The \(j\)th coalescence

With \(\binom{j}{2} = \frac{j(j-1)}{2}\) pairs available when \(j\) lineages remain:

\[f(T_j) = \frac{j(j-1)}{2}\exp\left\{ -\frac{j(j-1)}{2}T_j \right\}\]

All coalescences on a genealogy \(G\)

The intervals are independent, so the density of the whole set of waiting times is a product:

\[f(T \mid G) = \prod_{j=2}^{n} \frac{j(j-1)}{2}\exp\left\{ -\frac{j(j-1)}{2}T_j \right\}\]

Expectations for the TMRCA

Each interval is exponential, so its expectation is the reciprocal of its rate:

\[\mathbb{E}[T_j] = \frac{2}{j(j-1)}\]

and the time to the most recent common ancestor is the sum of the intervals:

\[\mathbb{E}[T_{MRCA}] = \mathbb{E}(T_n + T_{n-1} + \cdots + T_2)\]

Collapsing the sum

\[\mathbb{E}[T_{MRCA}] = \sum_{j=2}^{n} \frac{2}{j(j-1)} = 2\sum_{j=2}^{n}\left( \frac{1}{j-1} - \frac{1}{j} \right)\]

The sum telescopes:

\[\mathbb{E}[T_{MRCA}] = 2\left( 1 - \frac{1}{n} \right) \approx 2 \;\text{ coalescent units } = 4N \text{ generations}\]

Important conclusion

\(\mathbb{E}[T_{MRCA}] converges to 2\) coalescent units as \(n\) grows

This means that adding more samples does not add much information about deeper coalescent events. Ten lineages should already capture \(2(1-1/10)=1.8\) of the maximum 2.

From waiting times to branch length

So far we have the depth of the tree. That can be useful for estimating split times for speciation later on, but for now we only care about counting segregating sites. They occur on branches, so what we need is total branch length.

While \(j\) lineages remain, there are \(j\) branches accumulating length, each for a duration \(T_j\). So

\[\mathbb{E}[L] = \sum_{j=2}^{n} j \; \mathbb{E}[T_j]\]

This takes the form of the nth coalescent before.

The sum collapses again

Substituting \(\mathbb{E}[T_j] = \dfrac{4N}{j(j-1)}\) generations:

\[\mathbb{E}[L] = \sum_{j=2}^{n} j \times \frac{4N}{j(j-1)} = 4N\sum_{j=2}^{n} \frac{1}{j-1} = 4N\sum_{i=1}^{n-1} \frac{1}{i}\]

\[\boxed{\;\mathbb{E}[L] = 4N a_n, \qquad a_n = \sum_{i=1}^{n-1}\frac{1}{i}\;}\]

The \(j\) in the numerator cancels the \(j\) in \(j(j-1)\) — more lineages means more branches, but proportionally shorter intervals.

\(a_n\) grows like \(\log n\)

• Here \(n\) is logarithmic on the x axis, \(a_n\) is linear on the y axis.
• Going from \(n=2\) to \(n=10\) nearly triples the tree length.
• Going from \(n=10\) to \(n=200\) — twenty times the sequencing — only doubles it.

Thus, we do not emphasize sampling many individuals within a population.

Mutations land on the tree

Under the infinite-sites model, mutations arrive along branches at rate \(\mu\) per generation, independently of the genealogy.

So the number of segregating sites is just a thinned version of the total branch length:

\[\mathbb{E}[S] = \mu\,\mathbb{E}[L] = 4N\mu\,a_n\]

The definition of \(\theta\)

Writing \(\theta = 4N\mu\), the population mutation parameter:

\[\boxed{\;\mathbb{E}[S] = \theta\,a_n\;}\]

\(\theta\) compounds mutation and population size. A large population mutating slowly and a small population mutating fast produce the same spectrum of variation. Your organismal expertise is needed to separate the two.

Watterson’s estimator

Rearranging \(\mathbb{E}[S] = \theta a_n\) gives an estimator you can compute from a sequence alignment (Watterson 1975):

\[\boxed{\;\hat\theta_W = \frac{S}{a_n}\;}\]

• Unbiased by construction: \(\mathbb{E}[\hat\theta_W] = \theta\).
• Needs only \(S\), the count of variable sites, and \(n\).
• Per site, divide by sequence length.

Definition and neutral expectation

The SFS counts segregating sites by the number of copies \(i\) of the derived allele in a sample of \(n\). Under the neutral, constant-size coalescent:

\[ \boxed{\;\mathbb E[\xi_i] = \frac{\theta}{i}, \qquad i = 1, \dots, n-1\;} \]

• A mutation on a branch subtending \(i\) tips appears \(i\) times
• rare variants (singletons, \(i=1\)) should dominate the SFS
• \(\pi\) and \(\theta_W\) are different summary statistics from the same SFS: \(\theta_W \propto \sum_i \xi_i\) (counts sites equally), \(\pi \propto \sum_i i(n-i)\,\xi_i\) (weights intermediate frequencies).

Definition and neutral expectation

Simulated: \(n=50\), 2 Mb, constant \(N_e=10^4\)

Folded vs. unfolded

• Unfolded SFS needs an outgroup to polarize ancestral vs. derived (often difficult in practice).
• Folded SFS uses the minor-allele count \(\eta_i = \xi_i + \xi_{n-i}\)

A complication of the unfolded SFS is that errors in polarization can create a heavy tail, which looks like a historical bottleneck or a completed selective sweep.

Polarization errors in action

• Introduce 25% polarization error and the unfolded spectrum grows a fake high-frequency-derived tail.
• The folded spectrum is unchanged in both panels.

Construction

Both \(\pi\) and \(\theta_W\) estimate \(\theta\), so their difference reflects when evolution deviates from neutrality. Tajima (Tajima 1989) derived the standardization constant that allows the difference to be interpreted as a Z-score, where anything less than about \(-2\) or greater than about \(+2\) is a significant departure from neutrality:

\[ D = \frac{\pi - \theta_W}{\sqrt{\widehat{\mathrm{Var}}(\pi - \theta_W)}} . \]

• \(\pi\) is affected most by intermediate-frequency variants; \(\theta_W\) counts all segregating sites.
• So Tajima’s \(D\) asks: is the SFS skewed towards rare variants (\(\pi < \theta_W\)) or intermediate variants (\(\pi > \theta_W\))?

Interpreting the sign

“Bottleneck” is not one signature

20 replicates each; \(n=50\), 2 Mb.

Demography vs. Selection

• Demography (bottlenecks, expansions, structure) acts on every locus — it shifts the genome-wide SFS.
• Selection acts on specific loci and their linked neighborhoods — a local departure against the genome-wide background.
• ⟹ Build the genome-wide distribution of \(D\) (or \(F_{ST}\)), then flag outliers.

30 replicates, 50 kb windows, \(s=0.05\). Both lines are the same demography — only the sweep differs.

Computing Tajima’s \(D\) from an SFS

The sweep signature is an SFS skewed to rare variants. Compute \(D\) directly from a folded/unfolded spectrum:

TajimaD <- function(sfs) {              # sfs = c(#singletons, #doubletons, ...)
  n  <- length(sfs) + 1                 # number of sampled chromosomes
  ss <- sum(sfs)                        # segregating sites
  a1 <- sum(1 / seq_len(n - 1)); a2 <- sum(1 / seq_len(n - 1)^2)
  b1 <- (n + 1) / (3 * (n - 1)); b2 <- 2 * (n^2 + 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)
  theta_pi <- sum(sapply(seq_along(sfs), function(i) i*(n-i)*sfs[i])) / choose(n, 2)
  theta_w  <- ss / a1
  (theta_pi - theta_w) / sqrt(e1*ss + e2*ss*(ss-1))
}

Sweep vs. neutral — try it

sweep   <- c(20, 2, 1, 1, 1)   # rare-variant excess  → D < 0
neutral <- c( 6, 1, 2, 2, 4)   # flatter spectrum     → D ≈ 0
TajimaD(sweep)     # strongly negative
TajimaD(neutral)   # near zero

• The sweep SFS piles up singletons ⟹ \(\pi \ll \theta_W\) ⟹ negative \(D\).
• Change the counts and watch \(D\) move — this is the engine behind the genome scans in deck 05.

From one population to two

The SFS counts sites by the number of copies of an allele in one sample. The joint SFS counts them in two samples at once:

\[\xi_{i,j} = \text{number of sites with } i \text{ copies in population A and } j \text{ copies in B}\]

• Add up over B and you get back A’s ordinary SFS: the 1D spectra are the margins.
• Axes (\(i = 0\) or \(j = 0\)): private variants, segregating in one population only.
• Interior: shared variants, inherited from the common ancestor or carried across by gene flow.
• Opposite corners \((0, n_B)\) and \((n_A, 0)\): fixed differences between the populations.

What shapes a joint SFS

Simulated with msprime: 10 diploids per population, \(N_e = 10^4\), 20 Mb. The recent split is 0.1 coalescent units ago; the old split is 2.

Folding, again

Three populations, three pairwise spectra

• \(K\) populations give \(\binom{K}{2}\) pairwise spectra, and fastsimcoal2 fits them all together (Excoffier et al. 2013, 2021).
• Below: the paper’s split times and pre-crash sizes, with no recent crash; only gene flow differs.

Simulated with msprime (Baumdicker et al. 2022) on the mouse lemur north-vs-south tree. Top: no gene flow. Bottom: gene flow between Ambavala and the ancestor of the two southern populations.

The likelihood of an SFS

\[\ln L(\Theta) = \sum_{k} m_k \ln p_k(\Theta)\]

• \(m_k\): observed number of sites in cell \(k\) of every pairwise spectrum, including monomorphic sites.
• \(p_k(\Theta)\): probability that a site lands in cell \(k\) under parameters \(\Theta\). fastsimcoal2 estimates it by coalescent simulation, then climbs towards the best \(\Theta\).
• Every site is counted once in each pairwise spectrum, and nearby sites are linked, so this is a composite likelihood: good for finding \(\hat\Theta\), overconfident about uncertainty.
• Monomorphic sites plus \(\mu = 1.52\times10^{-8}\) (from a mouse lemur pedigree) turn \(\theta\) into individuals and generations; 3.5 years per generation gives years. As before: \(\theta\) compounds \(N_e\) and \(\mu\), and organismal knowledge separates them.

Case study: Goodman’s mouse lemur

• A small, forest-dependent primate of Madagascar’s eastern rainforests and isolated Central Highland forest patches (Tiley et al. 2022).
• Was its forest fragmented by paleoclimate (the LGM, 19–26.5 ka) or by people (~2 ka)?
• North vs. south: Ambavala (north) vs. Ambatovy and Tsinjoarivo (south). RADseq, with spectra estimated from genotype likelihoods.

Observed folded joint spectra and Hudson’s \(F_{ST}\) for each pair.

Eleven hypotheses

Blue arrows: gene flow. Dashed line: gene flow stops (models 8 and 10) or resumes (model 9). Red line: a recent size change in all three populations. Grey: the best model.

• Built hierarchically: which branches exchanged genes, then when, then a recent size change.
• Each model is a .tpl file (the history, written backward in time) plus an .est file (free parameters and search ranges).

Choosing among models

The Akaike information criterion (Akaike 1974) penalizes each free parameter; Akaike weights turn AIC differences into model probabilities:

\[\text{AIC} = 2k - 2\ln\hat L \qquad\qquad w_i = \frac{e^{-\Delta_i/2}}{\sum_j e^{-\Delta_j/2}}\]

What the best model says

Model 10 estimates and 95% CIs from Tiley et al. (2022) (Table S15), at 3.5 years per generation.

• North and south split ~277 ka; at the estimates, north–south gene flow ends about when Tsinjoarivo and Ambatovy split, ~7.8 ka (after the LGM).
• A 20- to 45-fold decline in Tsinjoarivo and Ambavala ~1.9 ka, near human arrival, but the CI spans 0.1–12 ka. Ambatovy did not decline.
• Tajima’s \(D\) tracks how small each population is now: Tsinjoarivo (smallest) \(+0.11\), Ambavala \(-0.20\), Ambatovy (no decline) \(-0.43\).

Winning is not fitting

Simulate the best model at its estimates and compare it with the data:

• Without the crash, the same tree gives \(F_{ST} \approx 0.035\) for the two southern populations: drift in a small population drives fast differentiation.
• AIC only ranks the models you thought of. Simulation asks whether the winner can actually produce the data.

From fastsimcoal2’s units to individuals and years

fastsimcoal2 reports what the data pin down, given the \(\mu\) fixed in the .tpl file:

• Drift is measured in coalescent units: \(T/N = T/2N_e\), as earlier in this lecture.
• We refitted model 10 with fastsimcoal 2.8 (10 replicates; practical §9–10). Its crash is more recent than Table S15’s: 332 generations, which add 0.097 units of drift to Tsinjoarivo, more than the 4,988 generations before it (0.074).
• Every branch is under 1 unit long, even the ~300,000-year-old one, so much ancestral variation is still shared.

Rescaling is not refitting

The SFS fixes \(N\mu\), \(T\mu\) and \(Nm\), not \(N\), \(T\) and \(m\) separately:

\[\mu \to c\,\mu \quad\Longrightarrow\quad N \to N/c, \qquad T \to T/c, \qquad m \to c\,m \qquad \text{(same likelihood)}\]

• Unchanged: fold declines, the order of events, \(Nm\), drift in coalescent units. Changed: individuals and years, which also scale with the generation time \(g\).

Event dates from our model-10 refit across illustrative mutation rates (lines) and generation times.

Rescaling in R

# our best model-10 replicate, exactly as fastsimcoal2 wrote it
est <- read.table("data/model10_best_fit/10.bestlhoods", header = TRUE)

mu <- 1.52e-8; g <- 3.5
est$Tau_Crash * g                     # years since the size change               ~1,160
est$Theta_NodeS / 2                   # diploid Ne of the southern ancestor       ~73,500
2 * est$Theta_NodeS * mu              # theta per site = 4 Ne mu                  ~0.0045
est$Tau_Crash / est$Theta_TsinCurr    # drift since the crash (coalescent units)  ~0.10

• Double \(\mu\): the first two numbers halve, while \(\theta\) and the drift do not change. They are what the data measured.

Caveats

• Composite likelihood: ~30 M linked sites make \(\Delta\)AIC and naive confidence intervals overconfident. The paper used a block bootstrap (blocks of 10,000 sites).
• Linked selection shifts the SFS genome-wide too: background selection (Charlesworth, Morgan, and Charlesworth 1993) can mimic a change in \(N_e\). Demographic models assume neutral sites.
• Identifiability: with biallelic SNPs, population size and split time can trade off against each other.
• Fold consistently: the model’s spectrum must be folded exactly as the data were (e.g. within each pair, as ANGSD does; fastsimcoal2’s --foldedSFS).
• Absolute \(N_e\) differs between methods (stairway plots vs. fastsimcoal2 in this study). Trust the pattern of change more than the exact numbers.

Conceptual review

• Under what conditions does \(\pi = \theta_W = \theta = 4N_e\mu\)?
• What does \(\theta_W\) depend on?
• What does Tajima’s \(D\) tell us about the skew of the SFS?
• What is a potential risk of the unfolded SFS?
• What are demographic scenarios where summary statistics like Tajima’s \(D\) have difficulty?
• How would a recent split and an old split with ongoing gene flow look different in a joint SFS?
• What additional information do you need to calibrate a demographic model to time in years?
• If a new study revises the mutation rate, which parameter estimates change, and which stay the same?