Accelerated Bayesian inference of population size history from recombining sequence data

Abstract

I present PHLASH, a new Bayesian method for inferring population history from whole genome sequence data. PHLASH is population history learning by averaging sampled histories: it works by drawing random, low-dimensional projections of the coalescent intensity function from the posterior distribution of a PSMC-like model, and averaging them together to form an accurate and adaptive size history estimator. On simulated data, PHLASH tends to be faster and have lower error than several competing methods including SMC++, MSMC2, and FITCOAL. Moreover, it provides a full posterior distribution over population size history, leading to automatic uncertainty quantification of the point estimates, as well to new Bayesian testing procedures for detecting population structure and ancient bottlenecks. On the technical side, the key advance is a novel algorithm for computing the score function (gradient of the log-likelihood) of a coalescent hidden Markov model: when there are $M$ hidden states, the algorithm requires $𝒪\left(M^2\right)$ time and $𝒪\left(1\right)$ memory per decoded position, the same cost as evaluating the log-likelihood itself using the naïve forward algorithm. This algorithm is combined with a hand-tuned implementation that fully leverages the power of modern GPU hardware, and the entire method has been released as an easy-to-use Python software package.

Figure D.1:

Simulation results for models 1–3. For $n=1000$, shaded regions indicate 95% credible intervals. For all other $n$, the solid (dotted) lines represent the median (min/max) across all simulation replicates.

Figure D.2:

Simulation results for models 4–6. See Figure D.1 for additional information.

Figure D.3:

Simulation results for models 7–9. See Figure D.1 for additional information.

Figure D.4:

Simulation results for models 10–12. See Figure D.1 for additional information.

Figure D.5:

$N_e$ estimates for archaic humans. Solid lines are posterior medians, and blue shaded area shows the 95% credible intervals for the filtered data. “Unfiltered” refers to the raw genotypes extracted from the Wohns et al. (2022) dataset, while “filtered” are the estimates obtained after running the filtering procedure described in the main text.

Figure 1:

Parallel algorithm for evaluating the likelihood of an HMM. The exact algorithm evaluates the likelihood sequentially across each observation (grey blocks in figure). The parallel algorithm breaks the observation sequence up into a set of overlapping blocks which can be evaluated in parallel. Blue blocks are used to “warm up” the filtering distribution, and red blocks contribute to the likelihood calculation.

Figure 2:

Example output of the algorithm for the “sawtooth” demography (Schiffels and Durbin, 2014; Adrion et al., 2020). $n=1000$ diploid samples were simulated according to the size history shown in black. The posterior distribution returned by PHLASH is summarized by its median (dark blue line), central 95% point-wise credible interval (shaded blue region), and 95% simultaneous credible band (dashed green line).

Figure 3:

$L_2$ versus total variation error for the African3Epoch_S16 demography when n = 1. phlash has large mean-squared error relative to msmc2 and smc++ (a), but it mainly occurs at very recent times, when there is low density of coalescence. Total variation error (b–c) accounts for this by measuring the $L_1$ distance between the true and inferred density functions.

Figure 4:

Benchmarking results. (a) Time needed to evaluate score function as a function of batch size. The sequence length was $L={10}^5$. Custom algorithm is compared forward- and reverse-mode automatic differentiation of the accelerated/linear-time forward algorithm using Jax (Bradbury et al., 2018). For $n=1000$, reverse-mode automatic differentiation exhausted the available GPU memory and was unable to be evaluated. (b) Mean CPU time and (c) peak memory usage for the various methods, with standard error bars.

Figure 5:

Accuracy of likelihood approximations. Panel (a) is a box-whisker plot of the relative error, defined as $\left|{\text{loglik}}_{\text{exact}}-{\text{loglik}}_{\text{parallel}}\right|/{\text{loglik}}_{\text{exact}}$, for various settings of the overlap parameter $f$, across all simulations. Panels (b)–(e) compare the results of model fitting using the composite likelihood approximate, versus independent likelihoods (see text for description).

Figure 6:

Posterior distribution of the recombination rate for each of the models shown in Table D.1.

Figure 7:

Results of running phlash on 3,609 genomes. Panel (a): estimates for all populations. Vertical shaded region is 813–930kya, see Section 4.4. (b) Estimates after grouping by super-population. (c-d): Estimates and 95% credible intervals for three (c) modern and (d) archaic subpopulations.

Figure 8:

Posterior inference of population structure. (a) Posterior distribution (median + 95% credible interval) of cross-coalescent rate between the Yoruba (YRI) and Han (CHB) populations. (b) Posterior distribution of non-parametric ratio estimator. (c) Three-population out-of-Africa model used to generate simulated data.

Figure 9:

Analyzing ancient DNA using phlash. a) Histogram of the different mutation types in the ancient samples. Enrichment for $\text{C}\to \text{T}$ mutations are likely the result of deamination, as are $\text{T}\to \text{C}$ mutations after accounting for uncertainty in the ancestral allele. $\text{G}\leftrightarrow \text{A}$ mutations are likely the result of deaminated bases on the complementary strand. b) Size history estimates for the Vindija Neanderthal before and after filtering out probable deaminated bases.

Figure 10:

Searching for signs of an ancient bottleneck in real and simulated human data. (a) Histogram of posterior median of ${\inf }_{t>t_{\text{anc}}}{N}_e\left(t\right)$ for each of the 159 populations considered. (b-c) The estimated size history of the Yoruba/Han population (black line, estimated from real data) was modified to incorporate a bottleneck. Orange, green, and blue correspond to bottleneck strength $\alpha \in \left\{{10}^0,{10}^{-1},{10}^{-2}\right\}$ respectively. Grey shaded region is $\left[813\times {10}^3,930\times {10}^3\right]$ years before present, assuming a generation time of 24 years. Dashed lines represent $\alpha N_e\left(t\right)$, the population size during the bottleneck, and solid lines represent the corresponding phlash estimates. (d) Mirror plot of posterior density of test statistic for Han and Yoruba. Shaded density curves are its distribution in simulated data; bar plots are a histogram of its distribution in real data.