The site frequency spectrum and demography

Counting sites by allele count

Sample n chromosomes at many loci. At each variable site, count how many chromosomes carry the derived (new) allele. The site frequency spectrum (SFS) is the number of sites ξᵢ for each count i from 1 to n − 1. The sites with i = 1 are singletons.

In a population of constant size with no selection, the expected spectrum is E[ξᵢ] = θ/i, where θ = 4N₀μ is the mutation rate for the whole locus. With n = 30 and θ = 10 that is 10 singletons, 5 doubletons and 3.3 tripletons per locus. Adding up all classes gives E[S] = θ × aₙ = 10 × 3.96 = 39.6 variable sites. The 1/i shape appears because most of the length of a genealogy lies in short branches near the tips, and a mutation on a tip branch is carried by one chromosome.

The module simulates 2000 independent loci from the coalescent, with one genealogy per locus and no recombination inside it. It draws mutations on the branches with a Poisson count, averages the spectrum over loci, and overlays the exact expectation for the chosen history.

How demography moves the tree

Size change alters how fast lineages merge, which alters the branch lengths. In the growth model the population was smaller in the past, so going back in time lineages merge quickly and the genealogy becomes star-like: long tips and a short trunk. Mutations fall mostly on tips, so singletons increase relative to everything else.

The summary statistics compare two estimates of θ. θ_W = S/aₙ counts variable sites, θ_π is the mean number of pairwise differences, and Tajima's D is their difference scaled by its standard deviation. Rare variants raise S and barely raise θ_π, so an excess of singletons makes D negative. A deficit of rare variants and an excess of intermediate ones makes it positive.

Fay and Wu's H reads the same spectrum differently. It compares θ_π with θ_H, an estimate that weights each class by i², so it responds to derived alleles at high frequency. A selective sweep produces many of them and makes H negative. A bottleneck gives a mean normalised H of −0.46, with 10% of loci below the null interval, so it mimics a sweep in H as it mimics balancing selection in D.

The module gives each statistic a null: a second set of 2000 constant-size loci with the same n and θ. A locus is flagged when its D falls outside the 2.5% and 97.5% quantiles of the null, which are −1.70 and 1.75 here.