Site Frequency Spectrum (SFS) Overview

• Site Frequency Spectrum (SFS) is a summary measure that records the number of polymorphic sites with a given count of derived alleles, reflecting the underlying genealogical branch lengths.
• Under the standard neutral coalescent, the expected SFS follows a 1/i law, illustrating that rare alleles are more common and linking mutation counts to branch-length decompositions.
• Extensions of the SFS encompass analysis of linked mutations, multiple populations, triallelic sites, and non-Kingman genealogies, thereby capturing diverse evolutionary scenarios.

The site frequency spectrum (SFS) is a summary of DNA sequence variation that records how many polymorphic sites have a derived allele at each possible frequency in a sample. For a sample of $n$ sequences, the unfolded SFS is the vector $(\xi_1,\dots,\xi_{n-1})$, where $\xi_i$ is the number of sites at which the mutant allele appears in exactly $i$ copies. Under the standard neutral coalescent with infinite-sites mutation, the SFS links mutational counts to the branch-length decomposition of the genealogy and yields the classical expectation $E[\xi_i]=\theta/i$, but a large literature has extended this object to linked sites, triallelic sites, multiple populations, multiple-merger coalescents, spatial genealogies, branching populations, and gene-presence/absence models (Rogers et al., 2021, Ferretti et al., 2016).

1. Definitions and genealogical representation

For a sample of size $n$, the unfolded SFS consists of $n-1$ frequency classes, with $\xi_1$ counting singletons, $\xi_2$ doubletons, and so on up to $\xi_{n-1}$ (Rogers et al., 2021). In a large population, allele frequency may instead be treated as continuous: if $(\xi_1,\dots,\xi_{n-1})$0 is the derived-allele frequency, then $(\xi_1,\dots,\xi_{n-1})$1 denotes the density of mutations at frequency $(\xi_1,\dots,\xi_{n-1})$2 along a sequence of length $(\xi_1,\dots,\xi_{n-1})$3 (Ferretti et al., 2016). When ancestral and derived states are not known, one works with the folded SFS, where counts $(\xi_1,\dots,\xi_{n-1})$4 and $(\xi_1,\dots,\xi_{n-1})$5 are merged as $(\xi_1,\dots,\xi_{n-1})$6 for $(\xi_1,\dots,\xi_{n-1})$7 (Shuai, 2023).

Under the infinite-sites model, the SFS is a Poisson randomization of branch lengths. If $(\xi_1,\dots,\xi_{n-1})$8 denotes the total length of all branches subtending exactly $(\xi_1,\dots,\xi_{n-1})$9 leaves in the genealogy, then mutations occurring on such branches are observed in exactly $\xi_i$0 sampled sequences, and $\xi_i$1 under mutation rate $\xi_i$2 per generation per lineage (Rogers et al., 2021). The same representation appears in more general genealogical settings: in coalescing Brownian motion, for example, $\xi_i$3 denotes the total length of all branches supporting exactly $\xi_i$4 leaves, and conditional on the tree the number of mutations inherited by exactly $\xi_i$5 individuals is Poisson with mean $\xi_i$6 (Shuai, 2023).

This branch-length view makes the SFS a genealogical rather than merely tabular object. It encodes how much tree length sits near the tips, at intermediate depths, or near the root, and consequently it is sensitive to the coalescent law, demographic history, linkage structure, mutation model, and even whether the genomic feature under study is always present in all sampled genomes (Baumdicker, 2014).

2. Standard neutral theory and the invariance property

Under selective neutrality, random mating, constant population size, the Kingman coalescent, and the infinite-sites mutation model, the expected unfolded SFS has the classical form

$\xi_i$7

with $\xi_i$8 in the diploid scaling used in the coalescent note of Rogers and Wooding (Rogers et al., 2021). Equivalently,

$\xi_i$9

so rare derived alleles are expected to be more numerous than common derived alleles (Rogers et al., 2021). Summing over classes gives the expected total number of segregating sites,

$i$0

A striking feature of this result is its invariance under sample-size increase. If one moves from sample size $i$1 to $i$2, then the expectations of the existing classes $i$3 do not change; one only appends the new class $i$4 (Rogers et al., 2021). Rogers and Wooding explain this through coalescent recursions for the expected number of sites at frequency $i$5 at the end of the interval with $i$6 lineages, $i$7, and show that for fixed $i$8, $i$9 does not depend on $E[\xi_i]=\theta/i$0 once $E[\xi_i]=\theta/i$1 (Rogers et al., 2021). In genealogical terms, the expected total branch length subtending exactly $E[\xi_i]=\theta/i$2 leaves is independent of the overall sample size beyond $E[\xi_i]=\theta/i$3, so adding more samples creates new higher-frequency classes without changing the expected amount of branch length supporting the old ones (Rogers et al., 2021).

A common misconception is that the $E[\xi_i]=\theta/i$4 law is a universal neutral baseline. The literature represented here is more restrictive: the formula is specific to the infinite-sites, neutral, constant-size, panmictic model, and the later sections show multiple distinct ways in which the spectrum departs from $E[\xi_i]=\theta/i$5 once those assumptions are relaxed [(Rogers et al., 2021); (Baumdicker, 2014)].

3. Extensions beyond the one-site spectrum

The one-site SFS admits several higher-dimensional generalizations. For two linked neutral sites without recombination, the expected neutral two-site spectrum (2-SFS) records either the population density $E[\xi_i]=\theta/i$6 or the sample counts $E[\xi_i]=\theta/i$7 of pairs of sites with derived frequencies $E[\xi_i]=\theta/i$8 or counts $E[\xi_i]=\theta/i$9 (Ferretti et al., 2016). A central decomposition separates nested mutation pairs, in which one mutation lies inside the background of the other, from disjoint pairs, in which the derived alleles never co-occur on the same chromosome: $n$0 This decomposition captures haplotype structure that is not identifiable from $n$1 alone (Ferretti et al., 2016).

Conditioning on a focal mutation of known frequency yields a conditional SFS of linked sites. In the no-recombination setting, linked mutations relative to the focal can be strictly nested, co-occurring, enclosing, complementary, or strictly disjoint, and exact formulas are available both in finite samples and in the population limit (Ferretti et al., 2016). These expressions provide a neutral reference for linked variation around structural variants, introgressed haplotypes, or chromosomal inversions. In the inversion application, the inversion frequency itself acts as the focal-site frequency, and the expected spectrum of neutral mutations inside the inverted region is exactly the conditional linked-site SFS (Ferretti et al., 2016).

A different extension allows up to two mutation events at a site and thereby produces a triallelic SFS. In that setting, the sample configuration $n$2 of ancestral allele $n$3 and two derived alleles $n$4 depends on whether the two mutation events are nested or nonnested on the underlying genealogy (Jenkins et al., 2013). Under general variable population size, the triallelic SFS can be written in closed form in terms of the second moments $n$5 of inter-coalescent times, in contrast to the usual biallelic SFS, which depends on first moments $n$6 (Jenkins et al., 2013). This suggests a sharper sensitivity to demographic history, and the numerical comparisons in that work indeed show larger Kullback–Leibler divergences between demographic models for the triallelic spectrum than for the diallelic one (Jenkins et al., 2013).

4. Alternative genealogies and nonclassical spectra

When the genealogy is not Kingman, the SFS changes qualitatively. In $n$7-coalescents, multiple ancestral lineages can merge at once, with merger rates determined by a finite measure $n$8 on $n$9 (Birkner et al., 2013). The expected SFS no longer has the simple $n-1$0 form; instead it is computed from recursions involving the block-counting process, the expected sojourn times $n-1$1, and the probability $n-1$2 that a branch at level $n-1$3 subtends $n-1$4 leaves (Birkner et al., 2013). For regularly varying $n-1$5-coalescents with exponent $n-1$6, including the Bolthausen–Sznitman coalescent, branch lengths of fixed order and hence fixed-frequency SFS entries have asymptotics very different from Kingman: $n-1$7 so

$n-1$8

for the Bolthausen–Sznitman case (Diehl et al., 2018). These asymptotics imply an extreme excess of low-frequency variants relative to Kingman.

Spatial expansion can also alter the spectrum. For genealogies modeled by coalescing Brownian motion on the line or circle, the total branch length $n-1$9 supporting exactly $\xi_1$0 leaves satisfies

$\xi_1$1

for each fixed $\xi_1$2 (Shuai, 2023). Thus, for fixed $\xi_1$3, the $\xi_1$4-class branch length is typically of order $\xi_1$5, in strong contrast to the order-one Kingman behavior for fixed frequency classes (Shuai, 2023).

Branching-population models replace coalescent sampling formulas by forward-time growth formulas. In a supercritical birth–death process with mutations only at birth, the large-sample SFS for a uniform sample of size $\xi_1$6 taken at time $\xi_1$7 has

$\xi_1$8

under $\xi_1$9, so the large-$\xi_2$0 shared-mutation spectrum has the $\xi_2$1 form familiar from exponential-growth branching models despite the underlying time-inhomogeneous ancestral reproduction rate (Shuai, 2024). For whole-population SFS in exponentially growing branching processes, almost sure limits are available both at fixed times and at hitting times of large population size, and in the pure-birth case the limiting per-individual SFS is proportional to $\xi_2$2 (Gunnarsson et al., 2023).

Tumor models sharpen this picture. In a birth–death branching process with cell division rate $\xi_2$3, death rate $\xi_2$4, net growth $\xi_2$5, and extinction probability $\xi_2$6, the expected SFS of the total population at the fixed-time scale $\xi_2$7 defined by $\xi_2$8 is

$\xi_2$9

(Gunnarsson et al., 2021). For large $\xi_{n-1}$0 this behaves like $\xi_{n-1}$1, whereas as $\xi_{n-1}$2 with fixed $\xi_{n-1}$3 it approaches $\xi_{n-1}$4, producing a transition between a $\xi_{n-1}$5 spectrum and a $\xi_{n-1}$6 spectrum that reflects cell viability (Gunnarsson et al., 2021). In a different rescued-population model with rare resistance mutations, fixed-$\xi_{n-1}$7 classes among resistant cells satisfy

$\xi_{n-1}$8

showing how rescue dynamics alter the normalization while preserving the branching-process kernel $\xi_{n-1}$9 for small families (Bonnet et al., 2023).

5. Coupled genomic features, linked markers, and dispensable genes

The classical SFS presumes that all sampled sequences carry homologous material at the locus under study. Dispensable genes violate that assumption. In the infinitely many genes plus infinitely many sites framework, each gene can be gained, lost, and, while present, accumulate site mutations (Baumdicker, 2014). The joint gene-and-site spectrum $(\xi_1,\dots,\xi_{n-1})$00 counts mutated sites whose gene is present in exactly $(\xi_1,\dots,\xi_{n-1})$01 sampled genomes and whose mutation is present in exactly $(\xi_1,\dots,\xi_{n-1})$02 of those $(\xi_1,\dots,\xi_{n-1})$03 carriers (Baumdicker, 2014). From this, the conditional within-gene SFS for a gene present in exactly $(\xi_1,\dots,\xi_{n-1})$04 out of $(\xi_1,\dots,\xi_{n-1})$05 individuals is derived explicitly; for $(\xi_1,\dots,\xi_{n-1})$06,

$(\xi_1,\dots,\xi_{n-1})$07

and even when $(\xi_1,\dots,\xi_{n-1})$08, the expectation differs from the classical $(\xi_1,\dots,\xi_{n-1})$09 formula if the gene-loss rate $(\xi_1,\dots,\xi_{n-1})$10 (Baumdicker, 2014).

This has immediate inferential consequences. Watterson’s estimator and Tajima’s $(\xi_1,\dots,\xi_{n-1})$11, when naively applied to dispensable genes, are both biased downward; the paper gives the explicit expectations and shows

$(\xi_1,\dots,\xi_{n-1})$12

(Baumdicker, 2014). The same work reports that the normalized SFS of dispensable genes shows an excess of low-frequency variants relative to the classical neutral SFS and that Tajima’s $(\xi_1,\dots,\xi_{n-1})$13 is shifted toward negative values even under strict neutrality, especially for rare genes and large $(\xi_1,\dots,\xi_{n-1})$14 (Baumdicker, 2014). A common interpretation of negative Tajima’s $(\xi_1,\dots,\xi_{n-1})$15 as evidence for purifying selection or expansion is therefore not portable to dispensable genes without adjusting the neutral baseline.

Linked markers generate a different kind of coupling. In the no-recombination 2-SFS framework, the SFS conditional on a focal neutral mutation gives a natural null model for linked variation around structural variants, introgressed segments, and inversions (Ferretti et al., 2016). The same paper notes an immediate reinterpretation in ancestral-recombination-graph language: because recombination events follow a Poisson process on branches like mutations, the conditional SFS of linked sites can describe the distribution of recombination events along a sequence (Ferretti et al., 2016).

6. Statistical uses, computation, and fundamental limits

Linear functionals of the SFS underpin mutation-rate estimation and neutrality testing. Under the standard coalescent, any statistic of the form

$(\xi_1,\dots,\xi_{n-1})$16

with $(\xi_1,\dots,\xi_{n-1})$17 and $(\xi_1,\dots,\xi_{n-1})$18 is an unbiased estimator of $(\xi_1,\dots,\xi_{n-1})$19 (Hobolth et al., 2021). Watterson’s estimator, the average pairwise-difference estimator, Fay and Wu’s $(\xi_1,\dots,\xi_{n-1})$20, Zeng et al.’s $(\xi_1,\dots,\xi_{n-1})$21, and unnormalized versions of Tajima’s $(\xi_1,\dots,\xi_{n-1})$22 all fit this framework (Hobolth et al., 2021). Multivariate phase-type theory provides an explicit probability generating function for the full SFS under the standard coalescent, and many 0–1 or non-negative integer weighted linear statistics are discrete phase-type distributed (Hobolth et al., 2021). Statistics with mixed-sign weights, such as Tajima’s $(\xi_1,\dots,\xi_{n-1})$23, are generally not discrete phase-type distributed, but their characteristic functions remain explicit and can be numerically inverted (Hobolth et al., 2021). The same framework is implemented in the R package phasty (Hobolth et al., 2021).

For structured demographies, the expected joint SFS across multiple populations can be computed efficiently from truncated single-population spectra and dynamic programming on a demographic DAG. New recursions for the truncated SFS reduce the core computation to $(\xi_1,\dots,\xi_{n-1})$24 per subpopulation, and these ideas are implemented in momi for joint SFS inference under arbitrary population-size histories, including piecewise exponential growth (Kamm et al., 2015). In tree-shaped demographies, forward-time Moran-model representations and FFT-based convolutions further reduce the cost of computing the expected joint spectrum (Kamm et al., 2015).

For $(\xi_1,\dots,\xi_{n-1})$25-coalescents, exact finite-$(\xi_1,\dots,\xi_{n-1})$26 recursions for expectations, variances, and covariances of SFS components extend the classical formulas of Fu from Kingman to multiple-merger genealogies (Birkner et al., 2013). The normalized expected SFS

$(\xi_1,\dots,\xi_{n-1})$27

is interpreted there as the probability that a randomly chosen mutation is seen in $(\xi_1,\dots,\xi_{n-1})$28 copies, and it becomes the basis of pseudo-likelihood inference for parameters such as the Beta-coalescent $(\xi_1,\dots,\xi_{n-1})$29 or the Eldon–Wakeley point-mass parameter $(\xi_1,\dots,\xi_{n-1})$30 (Birkner et al., 2013).

The SFS is also the basis of a large class of neutrality tests, but its informativeness has hard limits. A decomposition of neutrality tests into waiting-time and topology components shows that statistics such as Tajima’s $(\xi_1,\dots,\xi_{n-1})$31 and Fay and Wu’s $(\xi_1,\dots,\xi_{n-1})$32 depend directly on a measure of tree balance that is largely determined by the root balance of the genealogy (Ferretti et al., 2015). At the same time, there are information-theoretic lower bounds on what SFS-based demographic inference can achieve. For population histories with bottlenecks, the minimax error for estimating the size history from $(\xi_1,\dots,\xi_{n-1})$33 independent segregating sites is at least of order $(\xi_1,\dots,\xi_{n-1})$34, and the bound does not depend on the SFS dimension $(\xi_1,\dots,\xi_{n-1})$35 (Terhorst et al., 2015). This implies that, for fixed $(\xi_1,\dots,\xi_{n-1})$36, increasing the number of sampled individuals does not improve the minimax lower bound. A plausible implication is that linkage-aware or genealogy-aware summaries are necessary when one wishes to resolve demographic features that the SFS compresses too aggressively.

The contemporary view of the SFS is therefore twofold. It remains a central and highly tractable summary of mutational variation, with exact formulas, recursion schemes, and matrix methods available for many models (Rogers et al., 2021, Hobolth et al., 2021, Kamm et al., 2015). But its classical $(\xi_1,\dots,\xi_{n-1})$37 form, its standard inferential interpretations, and even its attainable accuracy are all model-dependent, and the modern literature treats the SFS less as a universal law than as a family of genealogy-dependent spectra whose shape changes in systematic and diagnostically useful ways across coalescent, branching, spatial, and pangenomic settings [(Ferretti et al., 2016); (Baumdicker, 2014); (Diehl et al., 2018); (Terhorst et al., 2015)].