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}}\),

\[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

Since \(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

\(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

The same geometric form, so the expectation 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 reads:

\[= \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.

Examples from corn domestication

Maize was domesticated from teosinte ~9,000 years ago through a bottleneck (Wright et al. 2005; Hufford et al. 2012):

• Genome-wide: the bottleneck lost diversity and skewed the SFS — a demographic signature present everywhere.
• At domestication loci (e.g. tb1, teosinte branched1): a selective sweep drove a far stronger local loss of diversity and strongly negative \(D\), standing out above the genome-wide background.

Examples from corn domestication

Significant evidence for a bottleneck in the history of corn domestication from Tajima’s \(D\) from Wright et al. (2005)

Examples from corn domestication

Do you expect more rare variants or intermediate variants based on this tree?

Neighbor joining tree of maize and teosinte from Hufford et al. (2012)

Examples from corn domestication

Signatures of selection recognizable from multiple lines of evidence:

Sweep for domestication loci on chromosomes one from Hufford et al. (2012)

Selection drags its neighbors

Selection acts on one site, but linkage spreads the consequences to neutral neighbors — the dominant force shaping diversity along a genome.

• Selective sweep (hitchhiking): a favored mutation sweeps to fixation, dragging its haplotype and erasing nearby diversity — a valley in \(\pi\) (Maynard Smith and Haigh 1974).
• Background selection: purging of linked deleterious mutations also removes neutral variation, depressing diversity in low-recombination regions (Charlesworth, Morgan, and Charlesworth 1993).

How deep and how wide?

\(T_{\text{fix}}\) is the number of generations for a new beneficial allele to sweep from a single copy to fixation. \(c\) is the recombination rate — the probability per generation of a crossover between the selected site and a neutral neighbor.

\[ T_{\text{fix}} \approx \frac{2\ln(2N_e)}{s}, \qquad \text{swept width} \sim \frac{s}{c}. \]

• Stronger selection \(s\) ⟹ faster fixation ⟹ deeper valley (less time for recombination to rescue neutral variants).
• More recombination \(c\) ⟹ narrower valley (neighbors escape the sweep).
• Diversity recovers with time \(T\) after the sweep.

How deep and how wide?

Recent sweep drops \(\pi\) 9-fold at the centre; 4,000 generations later it has refilled to within 45% of the background.

Hard, soft, and partial sweeps

• Hard sweep — a single new mutation fixes; one haplotype dominates. This is the strongest signal but kind of rare for natural plant populations.
• Soft sweep — selection on standing variation or recurrent mutation. Multiple haplotypes can carry the adaptive allele and the variation around the target makes a softer peak.
• Partial sweep — there is an adaptive allele but it remains at an intermediate frequency and does not reach fixation.

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.

Summary — key equations

\[ \theta = 4N_e\mu \qquad \mathbb E[\xi_i] = \frac{\theta}{i} \qquad D = \frac{\pi - \theta_W}{\sqrt{\widehat{\mathrm{Var}}(\pi - \theta_W)}} \qquad T_{\text{fix}} \approx \frac{2\ln(2N_e)}{s} \]

Summary - Review questions

• 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 Tajima’s \(D\) has difficulty distinguishing between them?
• How can you differentiate demographic effects from selection?
• Linked selection can affect upstream and downstream diversity, but what affects the width and depth of the affected region?