Bayesian Nonparametric Inference of Population Size Changes from Sequential Genealogies.

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 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 (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 )—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 ; Wilton ). Because of these features, in this work we rely on the SMC′ to model local genealogies with recombination.

Figure 1

SMC′ hidden Markov model for inferring population size trajectories, drawn according to Rasmussen to highlight notation specific to our study. (A) Observed sequence data in a segment of length L from five individuals. Three loci are shown delimited by recombination breakpoints and . Only the derived mutations at polymorphic sites are shown. (B) Corresponding local genealogies for each locus i. The five sampled individuals are depicted as solid black circles. Local genealogies have a Markovian degree-1 dependency. Each intercoalescent time (the time interval between coalescent events denoted as open circles) provides information about past population size (number of solid gray circles at a given time point). Moving from left to right after recombination breakpoint the pruning location is selected from genealogy and the pruned branch is regrafted back on the genealogy (solid blue circle). The coalescent event of depicted as a solid red circle in is deleted. This creates the next genealogy . The process continues until L. At L, the population size trajectory (depicted as a black curve superimposed on ) can be inferred.

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 from a sample of size 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 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 (Sheehan ; Rasmussen ), 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 ). 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 ; Lan ). 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 ; 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 ).

Methods: SMC′ Calculations

Following notation similar to that in Rasmussen (Table 1), a realization of the embedded SMC′ chain consists of a set of m local genealogies recombination breakpoints at chromosomal locations and pruning locations where indicates the time of the recombination event and the branch where recombination happened in genealogy (Figure 1). Genealogy corresponds to the genealogy of n sequences that contains the set of timed ancestral relationships among the n individuals for the chromosomal segment . Genealogy corresponds to the genealogy of the same n sequences for the chromosomal segment for . Finally, denotes the time when two of j lineages coalesce in genealogy measured in units of generations before present.

Table 1

Notation for the SMC′ model used in this work

Symbol	Description
Parameters
ρ	Recombination rate per site per generation
N(t)	Effective population size trajectory with time measured in units of N_0 generations
τ	Hyperparameter that controls the smoothness of the log-Gaussian process prior on N(t)
Notation specific to SMC′ chain
L	Length of observed sequences
b_i	Chromosomal location of the ith recombination breakpoint
m	No. local genealogies corresponding to m−1 recombination events
s_i+1=b_i+1−b_i	Segment length for local genealogy i
g_i	Local genealogy for the segment (b_i−1,b_i]
Notation specific to local genealogy
n	Sample size or no. sequences
l_i	Total tree length of local genealogy g_i
A^i(t)	Piecewise constant function of the number of ancestral lineages at time t in local genealogy g_i
tji	Coalescent time in genealogy g_i when two of j lineages coalesce. A^i(tji−)=j; A^i(tji+)=j−1.
t^i=(tni,tn−1i,…,t2i)	Vector of coalescent times of genealogy g_i
p_i=(u_i,w_i)	Pruning location along local genealogy g_i
u_i	Time when the recombination event happened along the height of the genealogy g_i
w_i	Lineage on genealogy g_i−1 where the recombination event happened
w^′_i	New lineage added on genealogy g_i where the recombination event happened
tnewi	Coalescent time in genealogy g_i when the lineage w_i coalesces
tdeli	Coalescent time in genealogy g_i−1 that no longer exists in genealogy g_i
c_i	Lineage on genealogy g_i that coalesces with lineage w^′_i
Fj,ki	No. free lineages in local genealogy g_i that do not coalesce in the time interval (tj+1i,tki)
I^i(t)	Piecewise constant function that takes values in {0, 1, 2} indicating no. ancestral lineages at time t in genealogy g_i where the pruning event would produce a visible transition to g_i+1
Discretization
d	No. change points at which N(t) is estimated
x=(x_1,…,x_d)	Times at which N(t) is estimated

Using uppercase letters to denote random variables, the evolution of the SMC′ process along chromosomal segments is governed by a point process that represents the random locations of recombination breakpoints. We use for to denote the segment lengths for each local genealogy, with Let be the chain that records the local genealogies, and let represent the chain that records the pruning locations (time and branch) on G. The sequence has the following conditional independence relation:Thus, given a chain of local genealogies, pruning locations, and recombination breakpoints, the joint transition probability to a new genealogy, pruning location, and locus length can be expressed as the product of the locus-length probability conditioned on the current genealogy (Expression 1, above), the pruning location probability conditioned on the current genealogy (Expression 2, above), and the transition probability of the new genealogy conditioned on the current genealogy and pruning location (Expression 3, above).

Complete data transition densities

Consider the chain of local genealogies with recombination breakpoints at . According to the SMC′ process, the first local genealogy follows the standard coalescent densitywhere and are the set of coalescent times in local genealogy . The piecewise constant function denotes the number of ancestral lineages present at time t in genealogy that is, with .

Given a current local genealogy the distribution of the length of the current locus depends on the current state of the SMC′ chain through the local genealogy’s total tree length (the sum of all branch lengths in ) and the recombination rate per site per generation ρ:At recombination breakpoint a new local genealogy is generated that depends on the previous local genealogy and the population size trajectory (Figure 1). To generate we first randomly choose a pruning location (consisting of a pruning time and a lineage ) uniformly along . At pruning location we add a new lineage and coalesce it further in the past at time with some lineage, (Figure 2). We then delete the lineage’s segment from to (the coalescent time of lineage ). The transition density to a new genealogy at recombination breakpoint is thenwhere denotes the total tree length of .

Figure 2

Schematic representation of SMC′ transitions given a recombination breakpoint at location (indicated as an arrow in each panel). (A) Visible transition. We uniformly sample the pruning location from at time along some branch and we add a new branch at and regraft it (dashed black line). The new branch coalesces with some branch at time . We then delete branch and the coalescent time to generate genealogy . Any pruning time along the branch (shown in green) would have produced the same visible transition from to . (B) Invisible transition. We uniformly sample the pruning location add a new branch at and regraft it. The new branch coalesces with itself (dashed black line), that is, and then the segment of is deleted. If any pruning location along the green branches would have produced the same invisible transition.

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

An invisible transition () occurs when . Given the pruning location an invisible transition occurs when and the random variable indicating the lineage that coalesces with lineage takes the value . The probability of an invisible transition is given byThus, the joint transition probability to an invisible event with pruning location given is

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 we need to average over all possible pruning locations along . By comparing the two genealogies and in Figure 2A, we know that corresponds to the lineage some time along or, equivalently, along . In general, comparison of and may not provide complete information to identify the lineage that was pruned. When the children of the node corresponding to and the children of the node corresponding to 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.

We introduce a function equal to the number of possible lineages at time t where the pruning location along would produce a visible transition to . is a piecewise constant function that takes the values in depending on whether the pruning location can happen in or 2 branches at time t. In the example in Figure 2A,For a general piecewise constant function the marginal visible transition density to a new genealogy is

Invisible transitions:

To compute the marginal transition probabilities for invisible events, we must average over all possible pruning locations . Consider the example in Figure 2B and choosing a pruning time () along . To have an invisible transition, the coalescing branch must be the same pruning branch . In Figure 2B the new coalescent time can happen along five lineages in the interval three lineages in the interval and two lineages in the interval . To generalize this calculation, we introduce the quantity with which denotes the number of lineages in that are free (do not coalesce), in the time segment with The time interval includes the interval of pruning up to the interval of self-coalescence . Thus, if the pruning time happens at time an invisible transition with new coalescent time can happen along free lineages.

In Figure 2B, happened in the time interval . If the new coalescent time happens in the interval along the same (unknown) pruning branch, then this invisible transition has probabilitywith

Now consider the same example of Figure 2B but with an unknown pruning time . The joint event where recombination occurs at pruning time and coalescent time occurs in the interval and this results in an invisible transition has probabilitywhere denotes the double integral expression in Equation 9 for ease of notation.

An invisible transition would also result if and along the same (unknown) pruning branch; in Figure 2B, this can happen along three lineages, so and this event has probabilityIf we continue considering the cases where and or we have and Then, the joint probability of an invisible event and isFor the cases when and the new coalescent time falls in another coalescent interval we need to compute the following: the joint probability of and no coalescence in the interval the probability of no coalescence in any of the intermediate coalescent intervals and the probability of coalescing at Then,represents the probability that the pruning location is at time and the new lineage coalesces at time with lineage . Overall, the marginal transition probability to an invisible event is

The likelihood of the embedded SMC′ chain

Instead of having a complete realization of the embedded SMC′ chain of m local genealogies and pruning locations at recombination breakpoints we assume that our data (unless otherwise noted) consist only of m local genealogies at recombination breakpoints from a chromosomal segment of length L (including visible and invisible events). Note that our observed data are not sequence data. More specifically, our observed data areThen, the observed data likelihood iswhere is the survival function in state . Equation 13 is factored into terms that depend on alone and ones that depend on ρ alone. The terms that depend on ρ, given by Equation 5, depend on the data only through total tree lengths and locus lengths . By the factorization theorem for sufficient statistics, local tree lengths and locus lengths are sufficient for inferring ρ. Moreover, recombination locations do not provide information about .

Methods: Inference

Current coalescent-based methods that infer a population size trajectory from whole-genome data assume is a piecewise constant function with change points (Li and Durbin 2011; Sheehan ; Rasmussen ; Schiffels and Durbin 2014). That is,Equation 14 presents two challenges. The first challenge lies in the specification of the change points: the narrower an interval is, the higher the probability that we do not observe coalescent times in that interval; further, the fewer observed coalescent times in an interval, the greater the uncertainty is of the estimate (if the estimate even exists). The second challenge lies in the specification of the time window if is set too far in the past, we might not have enough data to accurately estimate for .

To solve the first challenge, Li and Durbin (2011) and Rasmussen distribute the d change points evenly on a logarithmic scale,where κ is specified by the user. Schiffels and Durbin (2014) propose discretizing time according to the quantiles of the exponential distribution, where λ is the rate of an exponential distribution. Schiffels and Durbin (2014) model the time to the most recent coalescent event and set However, this equation is not directly applicable here because we use all coalescent events for 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

For our Bayesian methodology, we assume the log-Gaussian process prior on the population size trajectory,where denotes a Gaussian process with mean function 0 and inverse covariance function with precision parameter τ. For computational convenience, we use Brownian motion as our prior for since its inverse covariance matrix is sparse. We place a Gamma prior on the precision parameter τ, . Assuming that recombination rate ρ is known, the posterior distribution of model parameters (Figure 3) is thenThe first two factors on the right side of Equation 17, detailed in Equations 8 and 11, involve integration over an infinite-dimensional random function (Equation 16). We approximate the integralby the Riemann sum over a partition of the integration interval. That is,for and for is a representative value of in the interval in our implementation, we set with . This way, we discretize our time window in d evenly spaced segments with the maximum time to the most common ancestor observed in the sequence of local genealogies, and approximate by a piecewise linear function evaluated at .

Figure 3

Structure of our Bayesian model for inferring population size trajectories from a realization of the SMC′ process at recombination breakpoints. Hyperparameter τ controls the smoothness of the log-Gaussian process prior on . Local genealogies depend on and form a Markov chain of degree 1. Given the current local genealogy we sample the location of the new recombination breakpoint and a pruning location on genealogy . The new genealogy depends on and .

We condition on the set of m local genealogies (assuming pruning locations are not known) to generate posterior samples for the vector and τ and use these posterior samples to infer at where . Updating the vector and τ separately is not recommended because of their strong dependency (Lan ). Therefore, we update ( jointly in an MCMC sampling algorithm, using splitHMC (Shahbaba ; Lan ). splitHMC updates all model parameters jointly and it can be extended to a full inferential framework that is directly applicable to sequence data. The splitHMC method relies on Hamiltonian dynamics to propose a new state of the model parameters jointly with a higher acceptance rate than simple methods such as random-walk Metropolis (Neal 2009). splitHMC relies on our ability to calculate the log-likelihood of the observed data and the gradient vector of the log-likelihood (i.e., the score function). The log-likelihoods of the observed data are approximated via sums of the form in Equation 18. We approximate the score function with respect to by applying Fisher’s identity,where, at each iteration in the MCMC, expectation is calculated using the current value of (see Appendix).

Alternatively, one can update in the MCMC algorithm, using the elliptical slice sampler (Murray ) 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 with measures of uncertainty

We assume that the population size trajectory 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 ) to find the maximum-likelihood estimator of . The complete data for inferring are then the set of local genealogies and the set of pruning locations . For the invisible transitions, we also need to know the new coalescent times where denotes the set of indexes of invisible transitions.

The complete data log-likelihood is thenThe EM algorithm starts by initializing the population size trajectory to a piecewise constant function with change points with arbitrarily chosen vector . At the kth iteration of the algorithm we setThe conditional expectation in Equation 20 is conditional on the observed data defined in Equation 12. Let be the ordered set of time points corresponding to the change points and the coalescent time points of local genealogy i. If the transition from to is visible, we replace the jth time point by where j corresponds to the index such that . For ease of notation, we denote the number of time intervals by Letbe an indicator function that takes the value of 1 when the jth interval contains a coalescent time of the first genealogy . Then, the log density of the first genealogy isLetbe an indicator function that takes the value of 1 when the new coalescent time of genealogy i happens in the corresponding time interval and let the adjusted interval length beThen, the augmented transition density can be expressed aswhere and are the vectors with and elements. For the EM algorithm we need to compute the conditional expected vectors and The details of these calculations are in the Appendix.

We use the Fisher information matrix to compute approximate standard errors of and use these standard errors together with asymptotic normality of maximum-likelihood estimators to produce confidence intervals for log population size piecewise trajectories. We compute the observed Fisher information matrix following Louis (1982),where is the gradient and is the Hessian of the complete-data log-likelihood with respect to . This requires the calculation of conditional cross-product means and conditional second moments described in File S7.

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 ); see File S3 for details of these simulations. We assumed that all individuals were sampled at time

Table 2

Simulated demographic scenarios

Demographic model	N(t)
Constant population size	N(t)=1
Exponential growth followed by constant size	N(t)={1,for t∈(0,0.1),exp[−10(t−0.1)],for t∈(0.1,∞).
Population bottleneck	N(t)={1,for t∈(0,0.3),0.1,for t∈(0.3,0.5),1,for t∈(0.5,∞).

The argument t denotes time measured in units of generations.

We compared our point estimates with the truth for each demographic model, using the sum of relative errors (SRE),where is the estimated population size trajectory at time . We compute SRE at equally spaced time points . Second, we compute the mean relative width (MRW) aswhere corresponds to the upper limit and corresponds to the lower limit of . For EM estimates, corresponds to the 95% confidence interval estimated using the observed Fisher information; for Bayesian GP estimates, corresponds to the 95% Bayesian credible interval (BCI) of . To measure how well these intervals cover the truth, we compute the envelope measure (ENV) in the following way:We compute SRE, MRW, and ENV for at equally spaced time points.

For our Bayesian GP estimates, we estimate at time points, unless stated otherwise.

The parameters of the Gamma prior on the GP precision parameter τ were set to 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 with set to be the maximum observed coalescent time. For the cases where we consider only one genealogy (), 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 ; Palacios and Minin 2013).

Figure 5

Inference of population size trajectories for a pair of individuals (). Simulated data under constant population size (A), exponential and constant trajectory (B), and a bottleneck (C). We show estimates from and local genealogies. Dashed lines show the true trajectories, blue lines and light blue areas represent EM point estimates and 95% confidence areas, and red lines and pink areas represent Bayesian GP posterior medians and 95% BCIs. Boxplots of the TMRCA are shown at the bottom of each plot.

Figure 6

Inference of population size trajectories for Simulated data under constant population size (A), exponential and constant trajectory (B), and a bottleneck (C). We show estimates from genealogy, local genealogies, and local genealogies. Dashed lines show the true trajectories, blue lines and light blue areas represent EM point estimates and 95% confidence areas, and red lines and pink areas represent Bayesian GP posterior medians and 95% BCIs. Boxplots of the TMRCA are shown at the bottom of each plot.

Figure 7

Inference of population size trajectories for Simulated data under constant population size (A), exponential and constant trajectory (B), and bottleneck (C). We show estimates from genealogy, local genealogies, and local genealogies. Dashed lines show the true trajectories, blue lines and light blue areas represent EM point estimates and 95% confidence areas, and red lines and pink areas represent Bayesian GP posterior medians and 95% BCIs. Boxplots of the TMRCA are shown at the bottom of each plot.

Sensitivity of EM estimates of 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 (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.

Figure 4

Sensitivity to parameter discretization. Population size trajectories estimated from one simulated genealogy () of 100 individuals with a constant population size are compared. We show true trajectories as dashed lines. (A) Bayesian GP estimates at and 200 equally spaced time points. (B) EM estimates of a piecewise constant trajectory with change points and and 100 (Equation 15). (C) EM estimates of a piecewise constant trajectory with change points and and 500 (Equation 15). Point estimates are shown as solid black lines. 95% percent credible intervals and 95% confidence intervals are shown by gray areas.

Table 3

Summary statistics for simulation results depicted in Figure 4

Simulation of a single genealogy with n=100	SRE	MRW	ENV (%)
MLE d=5, κ=1	41.80	14.76	100.0
MLE d=5, κ=10	41.05	2.98	100.0
MLE d=5, κ=100	57.12	1.72	100.0
MLE d=10, κ=10	47.93	16.08	100.0
MLE d=10, κ=100	61.77	3.91	100.0
MLE d=10, κ=500	31.52	3.60	100.0
Bayesian GP d=50	6.98	1.88	100.0
Bayesian GP d=100	5.52	2.15	100.0
Bayesian GP d=200	4.96	1.70	100.0

SRE is the sum of relative errors (Equation 23), MRW is the mean relative width of the 95% BCI (Equation 24), and ENV is the envelope measure (Equation 25). Values in boldface type indicate best performance.

Comparing methods for estimating

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 () and using Equation 15 (similar to Li and Durbin 2011 and Rasmussen ). 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 () and 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 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 the Bayesian GP method accurately recovers sudden changes as shown in the bottleneck scenario. Although our Bayesian GP prior on is Brownian motion (which is not differentiable at any point), our Bayesian GP recovers smooth curves (Figure 5B).

Table 4

Summary of simulation results depicted in Figure 5

A. Simulations with n = 2
	SRE			MRW			ENV
Simulation and method	m = 100	m = 500	m = 1000	m = 100	m = 500	m = 1000	m = 100 (%)	m = 500 (%)	m = 1000 (%)
Const. EM	39.80	41.78	38.60	0.98	0.26	0.08	31.3	28.0	19.3
Const. GP	30.60	4.25	3.04	0.49	0.33	0.22	100.0	100.0	100.0
Exp. EM	64.68	25.70	33.70	0.91	0.16	0.12	42.0	26.0	6.6
Exp. GP	28.38	32.70	26.76	2.04	0.45	0.33	100.0	56.0	50.6
Bottle. EM	48.48	46.51	127.70	0.43	0.45	1.37	40.6	30.0	34.0
Bottle. GP	33.76	45.14	223.58	3.44	6.84	17.13	98.0	94.6	94.6
B. Simulations with n = 20
	SRE			MRW			ENV
	m = 1	m = 100	m = 1000	m = 1	m = 100	m = 1000	m = 1 (%)	m = 100 (%)	m = 1000 (%)
Const. EM	60.87	121.30	25.60	2.28	2.16	0.23	100.0	37.7	39.3
Const. GP	31.74	3.94	13.22	1.06	0.70	0.36	100.0	100.0	100.0
Exp. EM	40.97	40.66	40.22	3.11	0.37	0.19	100.0	38.6	19.3
Exp. GP	25.35	27.03	65.61	3.53	1.56	0.42	100.0	100.0	39.3
Bottle. EM	147.93	78.40	78.20	6.98	0.81	68.4	66.0	78.6	49.33
Bottle. GP	68.93	78.2	50.92	2.74	2.47	1.47	92.0	79.3	78.6
C. Simulations with n = 100
	SRE			MRW			ENV
	m = 1	m = 100	m = 1000	m = 1	m = 100	m = 1000	m = 1 (%)	m = 100 (%)	m = 1000 (%)
Const. EM	41.05	220.85	43.41	2.98	4.93	0.99	100.0	35.3	48.0
Const. GP	5.52	34.78	12.17	2.15	1.49	0.47	100.0	100.0	89.3
Exp. EM	76.86	40.22	27.63	3.23	0.81	0.13	87.3	42.0	14.0
Exp. GP	114.53	25.82	26.42	3.57	1.55	0.83	100.0	100.0	100.0
Bottle. EM	194.77	59.54	127.68	3.95	1.08	0.85	84.0	51.3	45.3
Bottle. GP	90.27	44.14	42.68	6.98	2.62	1.74	100.0	94.7	96.0

SRE is the sum of relative errors calculated as in (23), MRW is the mean relative width of the 95% BCI as Const: Constant population simulation scenario. Exp: Exponential growth followed by constant population size simulation scenario and Bottle: Population bottleneck. defined in (24), and ENV is the envelope measure calculated as in (25). Values in boldface type indicate best performance for each demographic model and sample size.

Table 4A shows the performance statistics for the estimates of 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 and (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 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 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 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 to improves accuracy: point estimates are closer to the truth (SRE in Table 4, A–C) and credible intervals cover the truth completely (ENV of ).

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

Under the SMC′, marginally at each locus along the chromosome, a local genealogy is a realization of Kingman’s n-coalescent (Kingman 1982), a continuous-time Markov chain taking its values in the set of sequences of partitions of the label set . A local genealogy g of n individuals includes labeled topology and coalescent times . The state space of a local genealogy is then and the cardinality of the set is . However, only the set of ordered coalescent times carries information about . For a single locus, the set of coalescent times provides sufficient statistics for inferring (see Proof in the Appendix). A natural question that follows is whether the coalescent times corresponding to the set of local genealogies are sufficient statistics for inferring under the SMC′ model. We find that the sufficient statistics for inferring under the SMC′ model are the coalescent times, when taken together with local ranked tree shapes (tree with no labels but ranked coalescent events). For a single locus, the set of coalescent times together with the ranked tree shape corresponds to a realization of Tajima’s n-coalescent. Tajima’s n-coalescent (Tajima 1983) is a continuous-time Markov chain taking its values in the set of ranked tree shapes [also called histories, evolutionary relationships, or vintaged and sized coalescent (Sainudiin )]. The state space of Tajima’s local genealogy is then and the cardinality of the set corresponds to the sequence of Euler zigzag numbers whose first 10 elements are (Disanto and Wiehe 2013). The probability of getting a particular type of ranked tree shape of n samples (Tajima 1983) is given bywhere c is the number of cherries, defined as branching events that lead to exactly two leaves.

We defined transition densities in terms of coalescent times and quantities (see Methods: SMC′ calculations). The set of all 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 under the SMC′ model (see Appendix). These observations are crucial for inferring 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 sequences, there are possible labeled topologies while only 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 ). We used ARGweaver (Rasmussen ) 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 and a recombination rate of (Rasmussen ) (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 = 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 ; Li and Durbin 2011; Rasmussen ; Sheehan ; Schiffels and Durbin 2014). In File S5, Figure S4A, we show the estimates of instead of and time measured in units of 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.

Figure 8

Inference of human population size trajectories for Green solid line and green areas represent the posterior median and 95% BCI for the European population (CEU) and blue solid line and blue areas represent the posterior median and 95% BCI for the Yoruban population (YRI). Time is measured in years in the past, assuming a generation length of 25 years and a reference diploid population of 10,000 individuals. The x-axis is log transformed.

Assessing the effect of using genealogies inferred with ARGweaver

We simulated sequence data, used ARGweaver for inferring a set of local genealogies, and used our method on those genealogies to obtain estimates of . To this end, we took the sequences of the first 1000 local genealogies of individuals simulated with MaCS as described in section 3 of File S5. We then generated the sequence lengths ( for the locus corresponding to ) as in Equation 5,where is the tree length of in units of generations and is the current population size. In our simulations, we set . To simulate sequence data of length over genealogy we used Seq-Gen (Rambaut and Grassly 1997) implemented in the R package phyclust (Chen 2011) from the Jukes–Cantor mutation model (Jukes and Cantor 1969) with mutation rate . We then used ARGweaver to infer a sample of local genealogies with the same corresponding parameters and with 200 change points for discretization of time. Figure 9 shows three estimations of effective population size trajectories for our three simulation scenarios. Figure 9, A–C, left, shows our GP-based estimates from 1000 simulated genealogies from MaCS; Figure 9, A–C, center, shows our GP-based estimates from a realization of local genealogies obtained from ARGweaver; and Figure 9, A–C, right, shows our GP-based estimates correcting the number of lineages used in our calculations, replacing by in our likelihood calculations. We find that our estimates are not only noisier using this approach but also biased.

Figure 9

Assessing the effect of a two-step inferential approach. (A–C) Left, our GP-based estimates when 1000 local genealogies are known; center, our GP-based estimates when ARGweaver estimated local genealogies are used for inference; right, our corrected GP-based estimates when ARGweaver estimated local genealogies are used for inference.

Discussion

In this article, we propose a Gaussian-process-based Bayesian nonparametric method for estimating effective population size trajectories 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 ) 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 ; 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 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 . 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 ( instead of 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 follows a log Brownian motion process. This allows us to model 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 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 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 ; Lan ). 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 from a single genealogy than elliptical slice sampling or regular Hamiltonian Monte Carlo sampling (Lan ). 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 ; Rasmussen ).

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 .

We used ARGweaver (Rasmussen ) 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 . 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 ; Gronau ; Sheehan , 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 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 .

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 . 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.