Inference Under a Wright-Fisher Model Using an Accurate Beta Approximation

Abstract

The large amount and high quality of genomic data available today enable, in principle, accurate inference of evolutionary histories of observed populations. The Wright-Fisher model is one of the most widely used models for this purpose. It describes the stochastic behavior in time of allele frequencies and the influence of evolutionary pressures, such as mutation and selection. Despite its simple mathematical formulation, exact results for the distribution of allele frequency (DAF) as a function of time are not available in closed analytical form. Existing approximations build on the computationally intensive diffusion limit or rely on matching moments of the DAF. One of the moment-based approximations relies on the beta distribution, which can accurately describe the DAF when the allele frequency is not close to the boundaries (0 and 1). Nonetheless, under a Wright-Fisher model, the probability of being on the boundary can be positive, corresponding to the allele being either lost or fixed. Here we introduce the beta with spikes, an extension of the beta approximation that explicitly models the loss and fixation probabilities as two spikes at the boundaries. We show that the addition of spikes greatly improves the quality of the approximation. We additionally illustrate, using both simulated and real data, how the beta with spikes can be used for inference of divergence times between populations with comparable performance to an existing state-of-the-art method.

Keywords: Wright-Fisher; beta; divergence times; linear evolutionary pressures; pure genetic drift.

Figure 1

Fit of the beta and beta with spikes approximations for a population of size $2N=200$. (A–F) The true discrete DAF as given by the Wright-Fisher model under pure genetic drift and the corresponding discretized beta (A–C) and beta with spikes (D–F) approximations. The distributions are conditional on an initial frequency $x_0=0.2$ and for different time points $t/2N=0.035$ (A, D), $t/2N=0.115$ (B, E), and $t/2N=0.195$ (C, F), where t is the number of discrete generations that the population has evolved. (G, H) The Hellinger distance between the true DAF and the beta (G) and beta with spikes (H) as a function of $x_0$ and $t/2N$. Each row and column corresponds to specific values of the scaled parameters $4Na$ and $4Nb$. The distances corresponding to the distributions from A–F are marked with arrows. The discretization procedure and Hellinger distance are detailed in File S1.

Figure 2

History of three populations in the present. Ancestral population 5 splits in populations 3 and 4, which further splits in populations 2 and 1. For each SNP i and present population $j\in \left\{1,2,3\right\}$, the data consist of the sample size $n_{ij}$ and allele count $z_{ij}$. The branch length between populations k and j is given as ${\left(t/2N\right)}_{k\to j}$ and represents the scaled number of generations that population j evolved since the split from the ancestral population k. The unknown allele frequencies of each population are denoted as $X_{ij}$, with $1\le j\le 5$.

Figure 3

Inference of divergence times for simulation scenarios I (A) and II (B). The figure shows box plots summarizing the inferred lengths for the four branches of the tree, indicated at the top of each column in black. The inferred lengths are plotted for beta (B), beta with spikes (BS), and Kim Tree (KT) (Gautier and Vitalis 2013). The (true) simulated length of each branch is plotted as a horizontal line. Each plot is scaled relative to the corresponding simulated branch length τ, with the limits of the y-axis set to $\left[\tau \cdot 0.1,\tau \cdot 1.5\right]$.

Figure 4

Inference of divergence times for the chimpanzee exome data. The figure shows box plots summarizing the inferred lengths using 50 data sets with 10,000 SNPs that were randomly sampled from the full data set. The corresponding tree branches are indicated at the top of each plot in black. The inferred lengths are plotted for beta (B), beta with spikes (BS), and Kim Tree (KT) (Gautier and Vitalis 2013). The nonsolid lines indicate the inferred lengths when running the methods on the full data set of 42,064 SNPs. The populations at the leaves are Eastern (E), Central (C), and Western (W). Each plot is scaled relative to the corresponding branch length τ inferred by beta with spikes on the full data set. The limits of the y-axis are set to $\left[\tau \cdot 0.05,\tau \cdot 1.5\right]$.