Bayesian Nonparametric Inference of Population Size Changes from Sequential Genealogies

Abstract

Sophisticated inferential tools coupled with the coalescent model have recently emerged for estimating past population sizes from genomic data. Recent methods that model recombination require small sample sizes, make constraining assumptions about population size changes, and do not report measures of uncertainty for estimates. Here, we develop a Gaussian process-based Bayesian nonparametric method coupled with a sequentially Markov coalescent model that allows accurate inference of population sizes over time from a set of genealogies. In contrast to current methods, our approach considers a broad class of recombination events, including those that do not change local genealogies. We show that our method outperforms recent likelihood-based methods that rely on discretization of the parameter space. We illustrate the application of our method to multiple demographic histories, including population bottlenecks and exponential growth. In simulation, our Bayesian approach produces point estimates four times more accurate than maximum-likelihood estimation (based on the sum of absolute differences between the truth and the estimated values). Further, our method’s credible intervals for population size as a function of time cover 90% of true values across multiple demographic scenarios, enabling formal hypothesis testing about population size differences over time. Using genealogies estimated with ARGweaver, we apply our method to European and Yoruban samples from the 1000 Genomes Project and confirm key known aspects of population size history over the past 150,000 years.

FOR a single nonrecombining locus, neutral coalescent theory predicts the set of timed ancestral relationships among sampled individuals, known as a gene genealogy (Kingman 1982; Hudson 1983, 1990; Tajima 1983). In the coalescent model with variable population size, the rate at which two lineages have a common ancestor (or coalesce) is a function of the population size in the past. Here we denote the population size trajectory by $N(t),$ where t is time in the past, and use the term local genealogy to describe ancestral relationships at one nonrecombining locus. When analyzing multilocus sequences, a single local genealogy will not represent the full history of the sample. Instead, the set of ancestral relationships and recombination events among a sample of multilocus sequences can be represented by a graph, known as the ancestral recombination graph (ARG), which depicts the complex structure of neighboring local genealogies and results in a computationally expensive model for inferring $N(t)$ (Griffiths and Marjoram 1997; Wiuf and Hein 1999).

Recent studies have leveraged approximations for the coalescent with recombination—the sequentially Markov coalescent (SMC) (McVean and Cardin 2005) and its variant SMC′ (Marjoram and Wall 2006; Chen et al. 2009)—both of which model local genealogies as a continuous-time Markov process along sequences (Figure 1). The difference between the SMC and the SMC′ is that the SMC models only the class of recombination events that alter local genealogies of the sample; in general, the SMC′ is a better approximation to the ARG than the SMC (Chen et al. 2009; Wilton et al. 2015). Because of these features, in this work we rely on the SMC′ to model local genealogies with recombination.

Under the coalescent and SMC′ models, population size trajectories and sequence data are separated by two stochastic processes: (i) a state process that describes the relationship between the population size trajectory and the set of local genealogies and (ii) an observation process that describes how the hidden local genealogies are observed through patterns of nucleotide diversity in the sequence data. The observation process includes mutation and genotyping error while the state process models coalescence. Population size trajectories are then inferred from sequence data, using these coalescent-based hidden Markov models. In this study, we restrict attention to the state process and present a novel Bayesian approach for inferring population size trajectories from local genealogies. We solve a number of key modeling and inference problems and thus provide a basis for developing efficient algorithms to infer population parameters from sequence data directly.

Whole-genome inference of population size trajectories has been hampered by the enormous state space of local genealogies for large sample sizes. The pioneering pairwise sequentially Markov coalescent (PSMC) method of Li and Durbin (2011) employed the SMC to infer $N(t)$ from a sample of size 2 (⁠$n=2$⁠). In this method, time is discretized and the population size trajectory is piecewise constant. Subsequent methods for samples larger than 2 similarly rely on the discretization of time. The natural extension of the PSMC to $n>2$ is the multiple sequentially Markovian coalescent (MSMC) (Schiffels and Durbin 2014). However, the MSMC models only the most recent coalescent event of the sample; thus MSMC’s estimation of population sizes is limited to very recent times. Other recent methods propose efficient ways of exploring the state space of hidden genealogies for $n>2$ (Sheehan et al. 2013; Rasmussen et al. 2014), yet also rely on discretizing the state space of local genealogies and assume a piecewise constant trajectory of population sizes.

Gaussian process-based Bayesian inference of population size trajectories has proved to be a powerful and flexible nonparametric approach when applied to a single local genealogy (Palacios and Minin 2013; Lan et al. 2015). The two main advantages of the Gaussian process (GP)-based approach are (i) it does not require a specific functional form of the population size trajectory (such as constant or exponential growth) and (ii) it does not require an arbitrary specification of change points in a piecewise constant or linear framework.

In this article, we overcome the limitations of existing methods—discretizing time, assuming a piecewise constant trajectory, and reporting only point estimates for past population sizes—by introducing a Bayesian nonparametric approach with a GP to model the population size trajectory as a continuous function of time. More specifically, we model the logarithm of the population size trajectory a priori as a Gaussian process (the log ensures our estimates are positive). As mentioned above, we assume that local gene genealogies are known. For our Bayesian approach, we develop a Markov chain Monte Carlo (MCMC) method to sample from the posterior distribution of population sizes over time. Our MCMC algorithm uses the recently developed algorithm, split Hamiltonian Monte Carlo (splitHMC) (Shahbaba et al. 2014; Lan et al. 2015). To compare our Bayesian GP-based estimation of population size trajectories with a piecewise constant maximum-likelihood-based estimation (e.g., Li and Durbin 2011; Sheehan et al. 2013; Schiffels and Durbin 2014), we implement the expectation-maximization (EM) algorithm within our framework and compute the observed Fisher information to obtain confidence intervals of the maximum-likelihood estimates.

Finally, we address a key problem for inference of population size trajectories under sequentially Markov coalescent models: the efficient computation of transition densities needed in the calculation of likelihoods. Here, we express the transition densities of local genealogies in terms of local ranked tree shapes (Tajima 1983) and coalescent times and show that these quantities are statistically sufficient for inferring population size trajectories either from sequence data directly or from the set of local genealogies. The use of ranked tree shapes allows us to exploit the state process of local genealogies efficiently since the space of ranked tree shapes has a smaller cardinality than the space of labeled topologies (Sainudiin et al. 2014).

Methods: SMC′ Calculations

Following notation similar to that in Rasmussen et al. (2014) (Table 1), a realization of the embedded SMC′ chain consists of a set of m local genealogies $(g0,g1,…,gm−1),$  $m−1$ recombination breakpoints at chromosomal locations $(b1,b2,…,bm−1),$ and $m−1$ pruning locations $(p1,p2,…,pm−1),$ where $pi=(ui,wi)$ indicates the time of the recombination event $ui$ and the branch $wi$ where recombination happened in genealogy $gi−1$ (Figure 1). Genealogy $g0$ corresponds to the genealogy of n sequences that contains the set of timed ancestral relationships among the n individuals for the chromosomal segment $(0,b1]$⁠. Genealogy $gi$ corresponds to the genealogy of the same n sequences for the chromosomal segment $(bi,bi+1]$ for $i=1,2,…,m−2$⁠. Finally, $tji$ denotes the time when two of j lineages coalesce in genealogy $gi,$ measured in units of generations before present.

Notation for the SMC′ model used in this work

Complete data transition densities

This generative process for local genealogies can result in a visible transition, where a genealogy $gi$ is different from $gi−1$ (Figure 2A), or an invisible transition, where $gi$ is identical to $gi−1$ (Figure 2B).

Transition densities averaged over unknown pruning locations

Even though we assume that local genealogies are known, to build inferential frameworks for sequence data in the future, we do not wish to make the same assumption about pruning locations. Thus, we average over pruning locations to obtain marginal transition densities between genealogies for both visible and invisible transitions.

Visible transitions:

To compute the marginal visible transition density to a new genealogy $gi={gi−1∖{tdeli}∪{tnewi,(w′i,ci)}},$ we need to average over all possible pruning locations $pi=(ui,wi)$ along $gi−1$⁠. By comparing the two genealogies $gi−1$ and $gi$ in Figure 2A, we know that $pi$ corresponds to the lineage $wi$ some time along $(0,t4i−1)$ or, equivalently, along $(0,tdeli)$⁠. In general, comparison of $gi−1$ and $gi$ may not provide complete information to identify the lineage that was pruned. When the children of the node corresponding to $tdel$ and the children of the node corresponding to $tnew$ are the same, pruning different branches can lead to the same transition. We enumerate all cases of incomplete information for visible transitions in Supporting Information, File S1, and File S2.

Invisible transitions:

To compute the marginal transition probabilities for invisible events, we must average over all possible pruning locations $pi$⁠. Consider the example in Figure 2B and choosing a pruning time (⁠$ui$⁠) along $gi−1$⁠. To have an invisible transition, the coalescing branch $Ci$ must be the same pruning branch $Wi$⁠. In Figure 2B the new coalescent time $Tnewi$ can happen along five lineages in the interval $(0,t5i−1),$ three lineages in the interval $(t5i−1,t4i−1),$ and two lineages in the interval $(t4i−1,t3i−1)$⁠. To generalize this calculation, we introduce the quantity $Fj,ki$ with $n≥j≥k≥2,$ which denotes the number of lineages in $gi$ that are free (do not coalesce), in the time segment $(tj+1i,tki),$ with $tn+1i=0.$ The time interval $(tj+1i,tki)$ includes the interval of pruning $(tj+1i,tji)$ up to the interval of self-coalescence $(tk+1i,tki)$⁠. Thus, if the pruning time happens at time $Ui∈(tj+1i,tji),$ an invisible transition with new coalescent time $Tnewi∈(tk+1i,tki)$ can happen along $Fj,ki$ free lineages.

The likelihood of the embedded SMC′ chain

Methods: Inference

In the following sections, we first present our Bayesian nonparametric method and then develop a maximum-likelihood method under a piecewise constant trajectory so we can directly compare an EM-based method to our Bayesian nonparametric method.

Gaussian-process-based Bayesian nonparametric estimation of $N(t)$

Alternatively, one can update $N(t)$ in the MCMC algorithm, using the elliptical slice sampler (Murray et al. 2010) with a fixed value of τ (perhaps estimated from previous studies or from a preliminary run from the split Hamiltonian Monte Carlo algorithm). The advantage of using the elliptical slice sampler over the split Hamiltonian Monte Carlo is purely computational (the elliptical slice sampler does not require calculation of the score function).

Maximum-likelihood estimation of $N(t)$ with measures of uncertainty

We assume that the population size trajectory $N(t)$ is defined as in Equation 14. The standard coalescent density (Equation 4) and the transition densities defined in Equations 8 and 11 are tractable, so calculation of the likelihood (Equation 13) is tractable. However, maximization of the likelihood function cannot be performed analytically because pruning locations are missing. We implement an EM algorithm (Dempster et al. 1977) to find the maximum-likelihood estimator of $N=(N1,…,Nd)$⁠. The complete data $Yc$ for inferring $N(t)$ are then the set of local genealogies $g0,…,gm−1$ and the set of pruning locations $p1…,pm−1$⁠. For the invisible transitions, we also need to know the new coalescent times ${tnewi}i∈ℐ,$ where $ℐ⊂{1,2,…,m−1}$ denotes the set of indexes of invisible transitions.

Data availability

The R code for all simulation studies and analysis of sequence data conducted in this article are publicly available at http://ramachandran-data.brown.edu/.

Results

We simulated 1000 local genealogies of 2, 20, and 100 individuals from each of the three different demographic models described in Table 2, using MaCS (Chen et al. 2009); see File S3 for details of these simulations. We assumed that all individuals were sampled at time $t=0.$

Simulated demographic scenarios

For our Bayesian GP estimates, we estimate $N(xi)$ at $d=100$ time points, unless stated otherwise.

The parameters of the Gamma prior on the GP precision parameter τ were set to $α=β=0.001,$ reflecting our lack of prior information about the smoothness of the population size trajectory.

For our EM estimates, we used different discretizations based on Equation 15 and varying the number of change points d and κ over the fixed interval $(0,xd)$ with $xd$ set to be the maximum observed coalescent time. For the cases where we consider only one genealogy (⁠$m=1$⁠), the EM approach becomes standard maximum-likelihood estimation.

We summarize our posterior inference and compare our Bayesian GP method to the EM method in Figure 5, Figure 6, and Figure 7. The population size trajectory is log-transformed for ease of visualization and for direct comparison with other methods (Minin et al. 2008; Palacios and Minin 2013).

Sensitivity of EM estimates of $N(t)$ to discretization

In Figure 4, we show our Bayesian GP and EM estimates of a constant population size trajectory from a single genealogy of 100 individuals with different discretizations. We find that our Bayesian GP point estimates depicted in Figure 4A recover the truth (dashed line) almost perfectly with less uncertainty than the EM (Figure 4, B and C). Comparing our Bayesian GP estimates with different discretizations [50, 100, and 200 equally spaced time points (Figure 4A)], we find that increasing the number of time points improves inference (Table 3) but that the differences between estimates among the three discretizations are marginal (Figure 4A). In contrast, we show that different grid definitions alter the EM estimates (Figure 4B). It is not clear how to define a good strategy for the definition of the grid for the EM method, even for the simple model of constant population size. For example, increasing κ from 100 to 500 with 5 change points (Figure 4B) does not improve estimation. Increasing the number of change points does not necessarily improve the estimates either, for example, increasing the number of change points from 5 to 10 for $κ=10$ (Figure 4, B and C). EM grid sensitivity is persistent even when the number of genealogies increases; Figure S2 in File S4 shows that the best definition of change points when our data consist of 1000 local genealogies of 100 individuals is 10 evenly distributed change points.

Summary statistics for simulation results depicted in Figure 4

Comparing methods for estimating $N(t)$

Figure 5 shows the estimated population size trajectories when the number of samples is two for the three different demographic scenarios and varying the number of local genealogies (100, 500, and 1000 local genealogies). For constant and exponential growth, our EM method assumes a piecewise constant trajectory of 10 change points (⁠$d=10$⁠) and $κ=1,$ using Equation 15 (similar to Li and Durbin 2011 and Rasmussen et al. 2014). For the bottleneck scenario, some of the intervals did not have coalescent events; hence, for this case we assumed a piecewise constant trajectory of 5 change points (⁠$d=5$⁠) and $κ=1$ for constructing our EM estimates. We show the boxplots of the time to the most recent common ancestor (TMRCA) at the bottom of each plot in Figure 5, which indicate the uncertainty expected in our estimates.

Both approaches, EM and Bayesian GP, show narrower confidence and credible intervals at the center of the distribution of the TMRCA, particularly during the bottleneck in Figure 5C.

For the constant population size model in Figure 5A, our Bayesian GP considerably outperforms our EM estimates. This is not surprising since a priori  $logN(t)$ has mean 0 in our Bayesian approach (Equation 16). Moreover, EM confidence intervals cover the truth only ∼30% of the time, while the GP method covers 100% of the truth (Table 4A). Despite placing a mean-0 prior on $logN(t),$ the Bayesian GP method accurately recovers sudden changes as shown in the bottleneck scenario. Although our Bayesian GP prior on $logN(t)$ is Brownian motion (which is not differentiable at any point), our Bayesian GP recovers smooth curves (Figure 5B).

Summary of simulation results depicted in Figure 5

Table 4A shows the performance statistics for the estimates of $N(t)$ in Figure 5. In general, our Bayesian GP has wider credible intervals than the EM confidence intervals but these credible intervals cover the true trajectory better than the EM confidence intervals in all cases (MRW and ENV in Table 4). Our Bayesian GP estimates also generally have smaller sums of relative errors (SRE in Table 4). Under the bottleneck scenario, our Bayesian GP produces greater sums of relative errors than does the EM, but our Bayesian GP estimates recover the truth more accurately than the EM during the bottleneck.

Figure 6 and Figure 7 show our estimates when $n=20$ and $n=100$ (Table 4, B and C, gives performance statistics). In general, our GP-based estimates have smaller SRE and larger ENV than the EM-based estimates and hence, the MRW is usually wider in the GP-based estimates, accurately reflecting the uncertainty of the estimates. As expected, increasing the number of loci (m) generally decreases the width of the confidence and credible intervals of our estimates (MRW). Although this is generally true for EM estimates as well, EM estimates have very low coverage of the truth (MRE in Table 4) when the number of loci increases.

Sampling more individuals vs. sequencing more loci

Figure 5, Figure 6, and Figure 7 show our estimates for $n=2,20,$ and 100 sampled individuals across varying numbers of loci. Performance of EM estimates depends strongly on the definition of the grid, so we focus here on the Bayesian GP estimates. We find that increasing the number of loci decreases uncertainty of our estimates and allows us to infer $N(t)$ farther back in time. Increasing the number of samples does not necessarily increase the performance of our GP estimates (File S6). For example, under the bottleneck scenario, we are able to detect the bottleneck fairly accurately even for two samples with $m=1000$ local genealogies. This is because most TMRCAs observed under the bottleneck scenario occur during the bottleneck (Figure 5, Figure 6, and Figure 7), regardless of the sample size. In contrast, in our exponential growth scenario, increasing the number of samples from $n=2$ to $n=100$ improves accuracy: point estimates are closer to the truth (SRE in Table 4, A–C) and credible intervals cover the truth completely (ENV of $100%$⁠).

Sequential Tajima’s genealogies are sufficient statistics under the SMC′

We defined transition densities in terms of coalescent times and $Fi,j$ quantities (see Methods: SMC′ calculations). The set of all $Fi,j$ quantities from a local genealogy forms a triangular matrix: an F matrix. We show that (i) F matrices are in bijection with ranked tree shapes and (ii) the set of local Tajima’s genealogies has sufficient statistics for inferring $N(t)$ under the SMC′ model (see Appendix). These observations are crucial for inferring $N(t)$ from sequence data directly. Coalescent-based inference from sequence data relies on marginalization over the hidden state space of genealogies. In the Appendix, we show that the state space needed is the space of local Tajima’s genealogies, as opposed to the space of local Kingman’s genealogies. For $n=10$ sequences, there are $2,571,912,000$ possible labeled topologies while only $7936$ possible ranked tree shapes.

Application to human data

We applied our method to a 2-Mb region on chromosome 1 (187,500,000–189,500,000) with no genes from five Yorubans from Ibadan, Nigeria (YRI) and five Utah residents of central European descent (CEU) from the 1000 Genomes pilot project (1000 Genomes Project Consortium 2012) and previously analyzed for the same purpose (Sheehan et al. 2013). We used ARGweaver (Rasmussen et al. 2014) to obtain a sample path of local genealogies for the two populations (YRI and CEU). The parameters used were 200 change points, a mutation rate of $μ=1.26×10−8,$ and a recombination rate of $ρ=1.6×10−8$ (Rasmussen et al. 2014) (File S5). We note that ARGweaver assumes the SMC process and our method assumes the SMC′ process. Moreover, our inference is based on a single sample of the SMC process with known pruning times. Our ARGweaver set of local genealogies is discretized at 200 time points and our GP-based inference is influenced by this discretization. In Figure 8 we show our estimates of past Yoruban (in blue) and European population sizes (in green). The two population size trajectories experience a series of bottlenecks and overlap until ∼100 KYA, assuming a diploid reference population size of $N0$ = 10,000 and a generation time of 25 years. In Figure 8 we recover an out-of-Africa bottleneck that starts ∼100 KYA and ends ∼30 KYA in the European population. These results are consistent with previously published results (Gronau et al. 2011; Li and Durbin 2011; Rasmussen et al. 2011; Sheehan et al. 2013; Schiffels and Durbin 2014). In File S5, Figure S4A, we show the estimates of $logN(t)$ instead of $N(t)$ and time measured in units of $N0$ generations (as in Figure 5, Figure 6, and Figure 7). We note that this two-step procedure of inferring local genealogies with ARGweaver and then using our method introduces biases and ignores genealogical uncertainty. In File S5, we correct for some of the bias caused by using this two-step procedure and show that our inferred population size trajectory remains valid for the recent past.

Assessing the effect of using genealogies inferred with ARGweaver

Discussion

In this article, we propose a Gaussian-process-based Bayesian nonparametric method for estimating effective population size trajectories $N(t)$ from a sequence of local genealogies, accounting for recombination. Under a variety of simulated demographic scenarios and sampling designs, our method recovers the truth with better precision and accuracy than a maximum-likelihood approach (Figure 5, Figure 6, and Figure 7). We apply our method to genealogies estimated using ARGweaver (Rasmussen et al. 2014) for European and African samples in the 1000 Genomes; this application to real data recovers the known features of the out-of-Africa bottleneck (Figure 8).

Several recent approaches have emerged for inferring population size trajectories from multiple whole-genome sequences using the SMC (Li and Durbin 2011; Sheehan et al. 2013; Schiffels and Durbin 2014). However, current SMC-based methods rely on maximum-likelihood inference (EM) of both a discretized parameter space and a discretized state space to gain computational tractability, and incur the costs of reduced accuracy and biased estimates. Although in principle the EM approach and the Bayesian nonparametric approach approximate $N(t)$ similarly—by either a piecewise constant or a piecewise linear function—the Bayesian nonparametric approach is not affected by increasing the number of parameters (or change points) in the estimation of $N(t)$⁠. For comparison with existing methods, we implemented an EM approach to infer population size trajectories from a sequence of local genealogies and we note that increasing the number of loci may actually increase the bias of the EM estimates (Figure 5, Figure 6, and Figure 7). For example, in simulation, our EM approach incorrectly detects the initial period of the simulated bottleneck (⁠$∼0.8N0$ instead of $0.5N0$ generations ago) with narrow confidence intervals (Figure 7C).

Using Bayesian GP for inferring population size trajectories offers many advantages over the EM approach. Similar to Palacios and Minin’s (2013) approach to inference from a single genealogy, we a priori assume that $N(t)$ follows a log Brownian motion process. This allows us to model $N(t)$ as a continuous positive function. The main advantage of using a Brownian motion process is that its inverse covariance function is a sparse matrix that allows for fast computations. Since the likelihood function involves integration over $N(t),$ this integral is approximated by the Riemann sum over a regular grid of points. The finer the grid is, the better the approximation. We find that our method performs well for inferring $N(t)$ at 100 change points in all our examples and, more importantly, results are not sensitive to the number of change points used in the analysis (Figure 4). Our Bayesian approach relies on MCMC for inference from the posterior distribution of model parameters. Because population sizes at different grid points are correlated, we adapt the recently developed MCMC technique splitHMC for jointly sampling all model parameters (Shahbaba et al. 2014; Lan et al. 2015). splitHMC is a Metropolis sampling algorithm that efficiently proposes states that are distant from current states with high acceptance rates. It has been shown to be more efficient in inferring $N(t)$ from a single genealogy than elliptical slice sampling or regular Hamiltonian Monte Carlo sampling (Lan et al. 2015). However, splitHMC relies on calculating the score function at every single iteration. Because the pruning time in each local genealogy is unknown, we calculate the score function via Fisher’s formula.

In simulations, we find that our algorithm scales well with hundreds of individuals; our computational bottleneck is in the number of local genealogies. We envision that extending the current methodology to inference from sequence data directly will require a strategy for sampling shorter genomic segments. This would be a probabilistic alternative to arbitrarily choosing segment lengths (Sheehan et al. 2013; Rasmussen et al. 2014).

Under the SMC model, every recombination event along the genome translates to a new coalescent event for the sample under study, so increasing the number of loci results in more realizations of the coalescent process. The longer the segments are and the larger the number of samples taken, the greater the chance of observing variation due to recombination. This fact makes it hard to define a sampling strategy: Longer genomes or larger sample sizes? We show that increasing the number of local genealogies improves precision of our Bayesian GP estimates (Figure 5, Figure 6, and Figure 7). However, resolution into the past from contemporaneous sequences highly depends on the actual population size trajectory $N(t)$⁠.

We used ARGweaver (Rasmussen et al. 2014) to generate two samples of contiguous local genealogies corresponding to a 2-Mb region of chromosome 1 for five Europeans (CEU) and five Africans (YRI) from the 1000 Genomes Project; this genomic region is free of genes and was also analyzed in Sheehan et al. (2013). Taking these two samples of local genealogies as our data (4186 local genealogies for CEU and 6247 local genealogies for YRI), we were able to use our Bayesian GP method to infer Yoruban and European effective population size trajectories (Figure 8). We find an out-of-Africa bottleneck that began ∼100 KYA and ended ∼30 KYA in the European population, consistent with Li and Durbin (2011); Rasmussen et al. (2011); Gronau et al. (2011); Sheehan et al. (2013), and Schiffels and Durbin (2014). We note that our estimates are based on a single sample of local genealogies and thus ignore genealogical uncertainty. Moreover, we generated our data from the posterior distribution of local genealogies, using ARGweaver at 200 time intervals, so our GP-based approach cannot fully detect sudden changes that may occur between the discretized times. In addition, ARGweaver assumes an SMC prior model on local genealogies and our GP-based method assumes the SMC′ process; the lack of invisible recombination events in ARGweaver’s genealogies will bias inference (as shown in simulated data in Figure 9).

The natural next extension for our method presented in this study is to infer $N(t)$ from sequence data directly and not from the set of local genealogies. Our MCMC approach allows us to extend the current methodology in a Bayesian hierarchical framework where the SMC′ process would be used as a prior distribution over local genealogies. The work we present here suggests a combination of ARGweaver accommodating SMC′ and GP priors would result in an efficient method for inferring population size trajectories from sequence data directly. In addition, our model can be easily modified to model a variable recombination rate along chromosomal segments and to jointly infer variable recombination rates and $N(t)$⁠.

Finally, we show that, under the SMC′ model, local ranked tree shapes and coalescent times correspond to a set of local Tajima’s genealogies; these Tajima’s genealogies are sufficient statistics for inferring $N(t)$⁠. Under the SMC′ model, the state space needed for inferring population size trajectories from sequence data is that of a sequence of local Tajima’s genealogies. This lumping, or reduction of the original SMC′ process, will allow more efficient inference from sequence data directly.

Current methods for inferring population size trajectories make trade-offs to analyze whole genomes that limit both biological understanding of sudden population size changes and the ability to test hypotheses regarding population size changes. This work represents a critical set of theoretical results that lay the groundwork for efficient estimation of detailed histories from sequence data with measures of uncertainty.

Appendix

Discretization

For both our Bayesian method and our EM method, we assume that $N(t)$ is a piecewise linear (or piecewise constant) function with d change points. Let $xi={x1i,x2i,…,xd+n−1i}$ be the ordered set of time points corresponding to the change points $x1,…,xd$ and the coalescent time points $ti$ of local genealogy i. Then, we calculate all the factors needed for the observed data likelihood (Equation 13) and the complete data likelihood (Equation 19).

Expectation-Maximization Algorithm

E step

M step

Observed Score Function for Split Hamiltonian Monte Carlo

Our Bayesian approach relies on splitHMC sampling from the posterior distribution of model parameters. This method requires the calculation of the observed score function. We use Fisher’s identity and calculate the observed score function as the conditional expected complete score function. The lth element of $∇ℒobs$ is

Sufficient Statistics Under SMC′

Here, we first formally show that under the standard coalescent and for a single locus, the set of coalescent times has sufficient statistics for inferring $N(t);$ that is, information about the topology is irrelevant for inference of $N(t)$⁠. We then investigate the properties of our F quantities (see Methods: SMC′ Calculations) needed for the calculation of the transition densities (Equation 11) and show through a series of propositions that the F quantities and coalescent times are the sufficient statistics for inferring $N(t)$ under the SMC′ process.

Proposition 1. For a single locus, the set of coalescent times are sufficient statistics for inferring $N(t)$⁠.

Proof. This can be proved using the factorization theorem. The marginal density of a local genealogy (Equation 3) has a unique factor that depends on $N(t)$ and g only through $tn,…,t2$⁠. The values of $A(t)$ are induced by the natural order of the coalescent times.

▪

Let F denote a lower triangular matrix of size $n×n$ with the $Fi,j$ entry the number of lineages that do not coalesce in the time interval $(ti+1,tj),$ as defined in Methods: SMC′ Calculations and with the following properties:

1. $Fi,1=0$ for all $i=1,…,n$ (The first column contains 0’s for completion).

2. $Fi,j=0$ for all $j>i$ (lower triangular matrix).

3. $Fi,i=i$ for all $i≥2$ (the diagonal corresponds to the number of lineages at each intercoalescent interval).

4. $Fi,i−1=i−2$ for all $i≥2$ (at each intercoalescent interval, we lose two free lineages, so the second diagonal correspond to the number of lineages minus two).

5. For $j<n−1,$ the last row of F is defined according to

   $Fn,j−1={Fn,j−2,with probability p=(Fn,j2)(j2),Fn,j−1,with probability p=Fn,j(j−Fn,j)(j2),Fn,j,with probability p=(j−Fn,j2)(j2).$
6. Let c denote the number of cherries; then

   $c=∑j=2n1{Fn,j−Fn,j−1=2}.$

7. For $i<n$ and $j<i−1,$ if $Fn,j−1=Fn,j−2,$ then $Fi,j−1=Fi,j−2.$

8. Let $vi$ denote the set of lineages in the intercoalescent interval $(ti,ti−1)$ with direct descendant internal nodes. The lineage labels correspond to the label of the coalescent time, when the direct descendant internal node was created. That is, the lineage created at $tn$ has label n: $vn={n};$ the lineage created at $ti$ has label i. Let $|vi|$ denote the size of the set $vi$⁠. Note that $1≤|vi|≤c$ and

   $|vi|=∑j=in1{Fn,j−Fn,j−1=2}−∑j=in1{Fn,j−Fn,j−1=0}.$
9. For $i<n$ and $j<i−1,$ if $Fn,j−1=Fn,j−1,$ then at time $tj,$ there is a coalescence between a singleton and a lineage in the set $vj$⁠. Let $aj$ be the lineage selected uniformly at random from $vj;$ then

   $Fi,j−1={Fi,j−1if i>ajFi,j−2if j<i≤aj.$
10. For $i<n$ and $j<i−1,$ if $Fn,j−1=Fn,j,$ then at time $tj,$ there is a coalescence between two lineages $aj1$ and $aj2$ from the set $vj$⁠. Let $aj1$ denote the minimum and $aj2$ the maximum of the two lineages selected; then

    $Fi,j−1={Fi,jif i>aj2Fi,j−1if aj1<i≤aj2Fi,j−2if j<i≤aj1.$

We show the correspondence between a ranked tree shape and the F matrix in the example of Figure A1. The first row and the first column are set to 0, and the first two diagonals are known with probability 1: $Fi,i=i$ and $Fi,i−1=Fi,i−2$ for $i>1.$ In our example, $n=5$ and so the first diagonal corresponds to $(0,2,3,4,5)$ and the second diagonal corresponds to $(0,1,2,3)$⁠. The last row, $F5,$ contains 0, followed by the number of branches that do not coalesce in the time intervals $(t6,t2),$  $(t6,t3),$  $(t6,t4),$ and $(t6,t5)$ corresponding to $(0,0,2,3,5)$⁠.

Proposition 2.  There is a bijection between the set of ranked tree shapes $Gn$  and  $ℱ,$  the set of  F  matrices.

▪

Note that the entries of the F matrix correspond to the same quantities needed to express the transition density of an invisible event (Equation 11). We claim that the sequence of coalescent times sets $t0,t1,…,tm−1$ and $F0,F1,…,Fm−1$ matrices corresponding to the ranked tree shapes of local genealogies $g0,g1,…,gm−1$ are sufficient statistics to infer $N(t)$ under the SMC′ process. We prove this through Propositions  3–6.

Proposition 3.  The probability density of Tajima’s genealogy is proportional, up to a combinatorial factor, to the probability density of Kingman’s genealogy.

Proposition 4.  The marginal visible transition density from a local Kingman’s genealogy $gi−1$  to  $Gi$  is proportional to the marginal visible transition density from  the corresponding local Tajima’s genealogy  $gi−1T$  to  $GiT$⁠.

▪

Proposition 5.  The marginal invisible transition density from a local Kingman’s genealogy $gi−1$  to  $Gi$  is equal to the marginal invisible transition density from the  corresponding local Tajima’s genealogy  $gi−1T$  to  $GiT$⁠.

▪

Proposition 6.  The likelihood of partially observed embedded SMC ′ chain of  local Kingman’s genealogies is proportional, up to a combinatorial factor, to the  likelihood of partially observed embedded SMC ′ chain of the corresponding local  Tajima’s  genealogies.

Proof. The proof follows from Propositions 3–5 needed to express the likelihood of partially observed embedded SMC′ chain (Equation 13).

▪

Acknowledgments

We thank Amandine Veber for her valuable suggestions and comments. We thank Shiwei Lan for his suggestions for speeding up the MCMC sampling scheme and Sara Sheehan and Melissa J. Hubisz for helpful discussions. We also thank two referees for suggestions that improved this article. J.A.P. acknowledges scholarship from Consejo Nacional de Ciencia y Tecnología Mexico to pursue her research work. This research is supported in part by National Science Foundation CAREER award DBI-1452622 (to S.R.). S.R. is a Pew Scholar in the Biomedical Sciences, supported by The Pew Charitable Trusts, and an Alfred P. Sloan Research Fellow.

Footnotes

Literature Cited