Efficient computation of the joint sample frequency spectra for multiple populations

Abstract

A wide range of studies in population genetics have employed the sample frequency spectrum (SFS), a summary statistic which describes the distribution of mutant alleles at a polymorphic site in a sample of DNA sequences and provides a highly efficient dimensional reduction of large-scale population genomic variation data. Recently, there has been much interest in analyzing the joint SFS data from multiple populations to infer parameters of complex demographic histories, including variable population sizes, population split times, migration rates, admixture proportions, and so on. SFS-based inference methods require accurate computation of the expected SFS under a given demographic model. Although much methodological progress has been made, existing methods suffer from numerical instability and high computational complexity when multiple populations are involved and the sample size is large. In this paper, we present new analytic formulas and algorithms that enable accurate, efficient computation of the expected joint SFS for thousands of individuals sampled from hundreds of populations related by a complex demographic model with arbitrary population size histories (including piecewise-exponential growth). Our results are implemented in a new software package called momi (MOran Models for Inference). Through an empirical study we demonstrate our improvements to numerical stability and computational complexity.

Keywords: coalescent; demographic inference; population genetics; sum-product algorithm.

Figure 1

A sample path of the coalescent truncated at time τ. Star symbols denote mutations, while ℳ^τ denotes the set of leaves under those mutations. $T_k^{\tau }$ denotes the waiting time in the interval [0, τ) while there are k lineages.

Figure 2

The coalescent with killing for the genealogy in Figure 1. Note that 𝒦_τ is a marked partition, with the blocks killed by mutations in [0, τ) being specially marked.

Figure 3

A demographic history with a pulse migration event (left), and its corresponding directed graph (right).

Figure 4

Average computation time of the joint SFS. For each combination of the sample size n and the number 𝒟 of populations (with $\frac{n}{𝒟}$ samples per population), we generated 20 random datasets, each under a demographic history that is a random binary tree. The expected joint SFS for the resulting segregating sites were then computed using our method (momi) and that of Chen (2012). In the top row, we plot the average runtime per joint SFS entry, and in the bottom row, the average amount of time needed to precompute the truncated SFS for every subpopulation within each demographic history. (a) Runtime results plotted separately for each method in a linear scale. Note the y-axis is on a different scale for each row. (b) Runtime results with the axes on a log-log scale, so that shorter runtimes are visible.

Figure 5

Numerical stability of the two algorithms. The plot compares the numerical values returned by our method (momi) and Chen’s method, for the simulations described in Figure 4. The dashed red line represents the identity y = x. To adequately illustrate the full numerical range (both positive and negative) of values encountered in the simulations, we applied the transformation z ↦ sign(z) log(1 + |z|) to the values of each method in order to produce the scatter plot. The two methods agree for n ≤ 64, while Chen’s method is extremely unstable for larger n.

Figure 6

Comparison of SFS values computed by ∂a∂i and momi, varying the sample size n per population and the number G of grid points per population in ∂a∂i. We used the 3-population out-of-Africa demography inferred in Gutenkunst et al. (2009), but modified to have no gene flow (migration). The x-axis is f_momi(x) the value computed by momi, the y-axis is the absolute value of the ratio $\left|\frac{f_{\partial a\partial i}(\mathbf{x})}{f_{\mathit{\text{momi}}}(\mathbf{x})}\right|$, and the color gives the sign of f_∂a∂i(x) (all values returned by momi were positive). f_∂a∂i generally converges to f_momi(x) as G increases, but for large values of n and G, ∂a∂i appears to diverge, and returns some negative values.

Figure 7

Runtime of momi and ∂a∂i to compute the results in Figure 6.