The sequential coalescent with recombination model

The coalescent-with-recombination (CwR) process, and the graph structures that are the realisations of the process, was first described by Hudson [42], and was given an elegant mathematical description by Griffiths [43], who named the resulting structure the Ancestral Recombination Graph (ARG). Like the coalescent process, these models run backwards in evolutionary time and consider the entire sequence at once, making it difficult to use them for inference on whole genomes. The first model of the CwR process that evolves sequentially rather than in the evolutionary time direction was introduced by Wiuf and Hein [44], opening up the possibility of inference over very long sequences. Like Griffiths’ process, the Wiuf-Hein algorithm operates on an ARG-like graph, but it is more efficient as it does not include many of the non-observable recombination events included in Griffiths’ process. The Sequential Coalescent with Recombination Model (SCRM) [45] further improved efficiency by modifying Wiuf and Hein’s algorithm to operate on a local genealogy rather than an ARG-like structure. Besides the “local” tree over the observed samples, this genealogy includes branches to non-contemporaneous tips that correspond to recombination events encountered earlier in the sequence. Recombinations on these “non-local” branches can be postponed until they affect observed sequences, and can sometimes be ignored altogether, leading to further efficiency gains while the resulting sample still follows the exact CwR process. An even more efficient but approximate algorithm is obtained by culling some non-local branches. In the extreme case of culling all non-local branches the SCRM approximation is equivalent to the SMC’ model [46, 47]. With a suitable definition of “current state” (i.e., the local tree including all non-local branches) these are all Markov processes, and can all be used in the Markov jump particle filter; here we use the SCRM model with tunable accuracy as implemented in [45].

The state space $$ \mathbb{T}$$ of the Markov process is the set of all possible genealogies at a given locus. The probability measure of a complete realisation x can be written as

Mutations follow a Poisson process whose rate at s depends on the state x_s via μ(s)B(x_s) where μ(s) is the mutation rate at s per nucleotide and per generation. Mutations are not observed directly, but their descendants are; a complete observation is represented by a set $$ y={\left\{(s_j,A_j)\right\}}_{j=1,\dots ,\left|y\right|}\in \mathcal{Y}$$ where s_j ∈ [1, L) is the locus of mutation j, and A_j ∈ {0, 1}^S are the wildtype (0) and alternative (1) alleles observed in the S samples. The conditional probability measure of the observations y given a realisation x is

Particle filters

Particle filters methods, also known as Sequential Monte Carlo (SMC) [22], generate samples from complex probability distributions with high-dimensional latent variables. An SMC method uses importance sampling (IS) to approximate a target distribution using weighted random samples (particles) drawn from a tractable distribution. We briefly review the discrete-time case. Suppose that particles {(x⁽ⁱ⁾, w⁽ⁱ⁾)}_i=1,…,N, approximate a distribution with density p(x), such that

Note that the algorithm can be seen as a recipe to transform a sample from P(X) to a sample from P(X)P(Y|X)/P(Y) = P(X|Y), that is, an application of Bayes’ theorem. Following this interpretation we will refer to P(X) as the prior distribution, and P(X|Y) as the posterior.

The algorithm generates an approximation to $$ P\left(x_{1··s}\right|y_{1··s})$$ rather than $$ P\left(x_s\right|y_{1··s})$$, but we follow [22] in calling it a particle filter algorithm instead of a smoothing algorithm (although our use of fixed-lag distributions for parameter estimation is a partial smoothing operation).

The marginal likelihood can be estimated (although with high variance, see [51]) by setting the weights to $$ N^{-1}{\sum }_i{w}_s^{\left(i\right)}$$ rather than N⁻¹ when particles are resampled. This makes the weights asymptotically normalized, so that (3) becomes E_P(X,Y=y)[f] ≈ ∑_i w⁽ⁱ⁾ f(x⁽ⁱ⁾), and $$ P(Y=y)=\int P(x,y)\text{d}x=E_{P(X,Y=y)}\left[1\right]\approx {\sum }_i{w}_L^{\left(i\right)}$$.

Algorithm 1 Particle filter

Input: $$ y_{1··L}$$

Output: Particles $$ \left\{(x_{1··L}^{\left(i\right)},w_L^{\left(i\right)})\right\}$$ approximating $$ P\left(x_{1··L}\right|y_{1··L})$$

 $$ w_0^{\left(i\right)}\leftarrow 1/N$$, $$ x_0^{\left(i\right)}\leftarrow \varnothing (i=1,\dots ,N)$$

 For s from 0 to L − 1

  Loop invariant: $$ \left\{(x_{1··s}^{\left(i\right)},w_s^{\left(i\right)})\right\}\sim p\left(x_{1··s}\right|y_{1··s})$$

  If ESS < N/2:

   Resample, with replacement, $$ \left\{x_{1··s}^{\left(i\right)}\right\}$$ proportional to $$ \left\{w_s^{\left(i\right)}\right\}$$

   $$ w_s^{\left(i\right)}\leftarrow N^{-1}(i=1,\dots ,N)$$

  For i from 1 to N:

   Sample $$ x_{s+1}^{\left(i\right)}\sim q\left(x_{s+1}\right|x_s^{\left(i\right)})$$

   $$ w_{s+1}^{\left(i\right)}\leftarrow w_s^{\left(i\right)}\frac{p\left(x_{s+1}^{\left(i\right)}\right|x_s^{\left(i\right)})}{q\left(x_{s+1}^{\left(i\right)}\right|x_s^{\left(i\right)})}g\left(y_{s+1}\right|x_{s+1}^{\left(i\right)})$$.

Continuous-time and -space Markov jump processes

For the hidden process we now consider Markov jump processes, which have as realisations piecewise constant functions $$ x:[1,L)\mapsto \mathbb{T}$$ where $$ \mathbb{T}$$ is the state space of the Markov process. Recall that in the model we consider, $$ \mathbb{T}$$ is the space of rooted genealogical trees with branch lengths. Let $$ (\mathcal{X},{\mathcal{F}}_x,{\pi }_x)$$ be a probability space, where $$ \mathcal{X}={\mathbb{T}}^{[1,L)}$$ is the space of possible realisations of the hidden stochastic process X = {X_s}_{s ∈ [1, L)}, $$ {\mathcal{F}}_x\subset \mathcal{P}\left(\mathcal{X}\right)$$ is the σ-algebra of events, and π_x(X) is the probability measure on $$ \mathcal{X}$$ induced by the stochastic process X. See the Appendix (“Conditional distributions and the Markov property”) for some remarks on how to define a Markov model when the phase space $$ \mathbb{T}$$ is uncountable.

The complete model is defined by specifying the observation process. We consider models where observations Y are generated by a Poisson process whose intensity at time (i.e. locus) s depends on X_s [a Cox process, see e.g. 52]. The space of observations $$ \mathcal{Y}$$ consists of finite subsets of [1, L) × M, where M is a discrete set of potential events, each of which may occur at some s ∈ [1, L). For a full observation $$ y=(({\tilde{s}}_1,m_1),\dots ,({\tilde{s}}_k,m_k))\in \mathcal{Y}$$ we write |y| ≔ k for the number of events in y. Writing λ(y) for the Lebesgue measure (ds)^|y|, the emission distribution π(Y|X = x) has a density r(y|x) relative to λ(y). For Cox processes this density has the form

The absence of events in an interval s ∈ [a, b) is also informative about the latent variable through the exponential factor in (8). In practice however, not all intervals may have been observed, so that events may or may not have occurred in these intervals. Assuming that the “observation process” is independent of the Markov jump process X, such unobserved intervals can simply be left out of integral (8).

Some more notation is needed to describe the Markov jump process version of algorithm 1. As above π_x denotes the prior distribution of the latent variable X, and ξ_x denotes the proposal distribution, both Markov processes on $$ \mathcal{X}$$, playing the role of p(x) and q(x) in the discrete case. We write $$ a:b$$ for the interval $$ [a,b)\subset \mathbb{R}$$, and α^{a: b} for the restriction of a measure or function α to $$ a:b$$; similarly y_a:b ≔ y∩([a, b) × M) and X_a:b ≔ {X_s}_s∈[a,b). The particle filter algorithm uses the notation (dα/dβ)(x) for distributions α and β to denote their Radon-Nikodym derivative: the ratio of their density functions with respect to a common reference measure, evaluated at x. To simplify notation we write the Radon-Nikodym derivative of two conditional distributions $$ \alpha \left(X\right|\mathcal{G})$$ and $$ \beta \left(X\right|\mathcal{G})$$ at x as $$ (\text{d}\alpha /\text{d}\beta )\left(x\right|\mathcal{G})$$, and we also do not explicitly restrict distributions to their appropriate intervals when this is clear from the context, so that we write for example $$ (\text{d}\pi /\text{d}\lambda )\left(y_{s_j:s_{j+1}}\right|X_{s_j:s_{j+1}}=x)$$ instead of $$ (\text{d}{\pi }^{s_j:s_{j+1}}\left(Y\right|X_{s_j:s_{j+1}}=x)/\text{d}{\lambda }^{s_j:s_{j+1}})\left(y_{s_j:s_{j+1}}\right)$$. With this notation we can formulate Algorithm 2.

Algorithm 2 Particle filter for Markov jump processes

Input: $$ y_{1:L}\in \mathcal{Y}$$; waypoints 1 = s₀ < s₁ < … < s_K = L.

Output: Particles $$ \left\{(x_{1:L}^{\left(i\right)},w_L^{\left(i\right)})\right\}$$ approximating the posterior distribution π(X|Y = y_1:L)

  $$ w_1^{\left(i\right)}\leftarrow N^{-1}$$, $$ x_{1:1}^{\left(i\right)}\leftarrow \varnothing $$ (i = 1, …, N)

 For j from 0 to K − 1

  Loop invariant: $$ \left\{(x_{1:s_j}^{\left(i\right)},w_{s_j}^{\left(i\right)})\right\}\approx \pi \left(X_{1:s_j}\right|Y_{1:s_j}=y_{1:s_j})$$

  If $$ ESS\left(\left\{w_{s_j}^{\left(i\right)}\right\}\right)<N/2$$:

   Resample $$ \left\{x_{1:s_j}^{\left(i\right)}\right\}$$ with probabilities proportional to $$ \left\{w_{s_j}^{\left(i\right)}\right\}$$

   $$ w_{s_j}^{\left(i\right)}\leftarrow N^{-1}$$ (i = 1, …, N)

  For i from 1 to N:

   Sample $$ x_{s_j:s_{j+1}}^{\left(i\right)}\sim {\xi }_x\left(X_{s_j:s_{j+1}}\right|X_{s_j}=x_{s_j}^{\left(i\right)})$$

   $$ {\displaystyle w_{s_{j+1}}^{\left(i\right)}\leftarrow w_{s_j}^{\left(i\right)}\frac{\text{d}{\pi }_x}{\text{d}{\xi }_x}\left(x_{s_j:s_{j+1}}^{\left(i\right)}\right|X_{s_j}=x_{s_j}^{\left(i\right)})\frac{\text{d}\pi }{\text{d}\lambda }\left(y_{s_j:s_{j+1}}\right|X_{s_j:s_{j+1}}=x_{s_j:s_{j+1}}^{\left(i\right)})}$$.

The choice of waypoints s₁, …, s_K is discussed below; in particular they need not be the same as the event loci $$ {\tilde{s}}_1,\dots ,{\tilde{s}}_{\left|y\right|}$$ of the observation y. Note that there is no initialization step; instead, initially $$ x_{1:1}^{\left(i\right)}=\varnothing $$, and the first sample will be drawn from ξ conditioned on an empty set, i.e. the unconditional distribution. The loop invariant holds when j = 0 since $$ 1:s_0=\varnothing $$. As with Algorithm 1 it is possible to estimate the likelihood density π_θ(y_1:L) by replacing the factors N⁻¹ with $$ N^{-1}{\sum }_i{w}_{s_j}^{\left(i\right)}$$; then the likelihood density w.r.t. λ(dy) = (ds)^|y| is approximated by $$ {\sum }_i{w}_L^{\left(i\right)}$$.

Note that by the nature of Markov jump processes, particles that start with identical latent variables have a positive probability of remaining identical after a finite time. Combined with resampling, this causes a considerable number of particles to have one or more identical siblings. For computational efficiency we represent such particles once, and keep track of their multiplicity k. When evolving a particle with multiplicity k > 1, we increase the exit rate k-fold, and when an event occurs one particle is spawned off while the remaining k − 1 continue unchanged.

Using lookahead to improve the particle filter

At the jth iteration, Algorithm 2 uses data up to waypoint s_j to build particles approximating $$ \pi \left(X_{1:s_j}\right|Y_{1:s_j}=y_{1:s_j})$$. This is reasonable as $$ \pi \left(X_{1:s_j}\right|y_{1:s_j})$$ is independent of data beyond s_j. However, not all particles are equally important for approximating subsequent posteriors, which suggests to emphasise particles that will be relevant in future at the expense of those relevant only to $$ \pi \left(X_{1:s_j}\right|y_{1:s_j})$$. This echoes the justification of resampling: although resampling increases the variance of the approximation to the current partial posterior, the variance at subsequent iterations by increasing the number of particles that are likely to contribute to future distributions. For discrete-time models p(X_1:n|y_1:n), the Auxiliary Particle Filter (APF) [40] implements this intuition by targeting a resampling distribution [53], which includes a “lookahead” factor $$ \tilde{p}\left(y_{i+1}\right|x_i)$$ approximating the probability of observing data y_i+1 given the current state x_i. Importance sampling is used to keep track of the desired distribution p(X_1:i|y_1:i).

In the continuous-time context it is natural to look an arbitrary distance ahead. Similar to APF, the lookahead distribution can be conditioned on the current state only, and must be an approximation of the true distribution. It should be heavy-tailed with respect to the true distribution to ensure that the variance of the estimator remains finite [22], which implies that the distribution should not depend on data too far beyond s; what is “too far” depends on how well the lookahead distribution approximates the true distribution.

The lookahead distribution is only evaluated on a fixed observation y, and is used to quantify the plausibility of a current state $$ x_s^{\left(i\right)}$$, rather than to define a distribution over y. For this reason we call it a lookahead likelihood. In fact, for correctness of the algorithm it is not necessary that this likelihood derives from a probability distribution. We define the lookahead likelihood as a family of functions $$ h^s\left(y_{s:L}\right|x_s):{\mathcal{Y}}_{s:L}\times \mathbb{T}\to \mathbb{R}$$, and an associated family of unnormalized distributions $$ {\tilde{\pi }}^s(x_{1:s},y_{1:L})={\pi }^{1:s}(x_{1:s},y_{1:s})h^s\left(y_{s:L}\right|x_s){\lambda }^{s:L}\left(y_{s:L}\right)$$ on $$ {\mathcal{X}}_{1:s}\times \mathcal{Y}$$. The functions h^s can be chosen arbitrarily, except that h^s(⋅, x_s)λ^s:L must be absolutely continuous w.r.t. π^s:L(⋅|X_s = x_s) to ensure that importance sampling is justified. The lookahead Algorithm 3 keeps track of two sets of weights, which together with a single set of samples form two sets of particles that approximate the resampling and target distributions.

Algorithm 3 Markov-jump particle filter with lookahead

Input: $$ y_{1:L}\in \mathcal{Y}$$; waypoints 1 = s₀ < s₁ < … < s_K = L.

Output: Particles $$ \left\{(x_{1:L}^{\left(i\right)},w_L^{\left(i\right)})\right\}$$ approximating π(X|Y_1:L = y_1:L)

 $$ w_1^{\left(i\right)}\leftarrow 1/N$$, $$ v_1^{\left(i\right)}\leftarrow 1/N$$, $$ x_{1:1}^{\left(i\right)}\leftarrow \varnothing $$ (i = 1, …, N)

 For j from 0 to K − 1

  Loop invariant: $$ \left\{(x_{1:s_j}^{\left(i\right)},w_{s_j}^{\left(i\right)})\right\}\approx {\pi }^{1:s_j}\left(X_{1:s_j}\right|Y_{1:s_j}=y_{1:s_j})$$

  Loop invariant: $$ \left\{(x_{1:s_j}^{\left(i\right)},v_{s_j}^{\left(i\right)})\right\}\approx {\tilde{\pi }}^{s_j}\left(X_{1:s_j}\right|Y=y_{1:L})$$

  If $$ ESS\left(\left\{v_{s_j}^{\left(i\right)}\right\}\right)<N/2$$:

   Resample $$ \left\{x_{1:s_j}^{\left(i\right)}\right\}$$ with probabilities proportional to $$ \left\{v_{s_j}^{\left(i\right)}\right\}$$

   $$ w_{s_j}^{\left(i\right)}\leftarrow N^{-1}{w}_{s_j}^{\left(i\right)}/v_{s_j}^{\left(i\right)}$$ (i = 1, …, N)

   $$ v_{s_j}^{\left(i\right)}\leftarrow N^{-1}$$ (i = 1, …, N)

  For i from 1 to N:

   Sample $$ x_{s_j:s_{j+1}}^{\left(i\right)}\sim {\xi }_x\left(X_{s_j:s_{j+1}}\right|X_{s_j}=x_{s_j}^{\left(i\right)})$$

   $$ {\displaystyle w_{s_{j+1}}^{\left(i\right)}\leftarrow w_{s_j}^{\left(i\right)}\frac{\text{d}{\pi }_x}{\text{d}{\xi }_x}\left(x_{s_j:s_{j+1}}^{\left(i\right)}\right|X_{s_j}=x_{s_j}^{\left(i\right)})\frac{\text{d}\pi }{\text{d}\lambda }\left(y_{s_j:s_{j+1}}\right|X_{s_j:s_{j+1}}=x_{s_j:s_{j+1}}^{\left(i\right)})}$$

   $$ {\displaystyle v_{s_{j+1}}^{\left(i\right)}\leftarrow v_{s_j}^{\left(i\right)}\frac{\text{d}{\pi }_x}{\text{d}{\xi }_x}\left(x_{s_j:s_{j+1}}^{\left(i\right)}\right|X_{s_j}=x_{s_j}^{\left(i\right)})\frac{\text{d}{\tilde{\pi }}^{s_{j+1}}}{\text{d}{\tilde{\pi }}^{s_j}}\left(y_{s_j:L}\right|X_{s_j:s_{j+1}}=x_{s_j:s_{j+1}}^{\left(i\right)})}$$

(see Appendix, “Proof of Algorithm 3”.) To implement the lookahead particle filter we need a tractable approximate likelihood of future data given a current genealogy. To do this we simplify the full likelihood, and ignore all data except for a digest of singletons and doubletons that are informative of the topology and branch lengths near the tips of the genealogy—in particular, singletons are informative of terminal branch lengths, and doubletons identify the existence of nodes with precisely two descendants (“cherries”). This digest consists of the distance s_i to the nearest future singleton for each haploid sequence, and ≤ n/2 mutually consistent cherries c_k = (a_k, b_k) identified by their two descendants a_k, b_k, together with loci $$ s_k^{\prime }\le s_k^{″}$$ where their first and last supporting doubleton were observed (Fig 1a). Under some simplifying assumptions we derive an approximation of the likelihood $$ h^s(\left\{t_i\right\},\{c_k,s_k^{\prime },s_k^{″}\}|x_s)$$ of the current genealogy given these data; see Appendix (“A lookahead likelihood”) for details.

Fig 1

a. Example of data digest. Lines represent genomes of 6 lineages, circles observed genetic variants. Of the data shown, one singleton (yellow) and five doubletons (red) contribute to the digest. Cherry c₃ is supported by a single doubleton; r does not contribute because the mutation patterns p and q are incompatible with c₃. Similarly, p does not contribute because it is incompatible with c₂ and c₃. b. Partial genealogy (unbroken lines) over 6 lineages. Open circles and arrowheads represent potential recombination and coalescence events that would change the terminal branch length for lineage 1 (t,u), and remove cherry c₃ (x,y).

Choosing waypoints

The choice of waypoints s_j can significantly impact the performance of the algorithm: choosing too few increases the variance of the approximation, and choosing too many slows down the algorithm without increasing its accuracy. Waypoints determine where the algorithms might perform a resampling step. A high density of waypoints is therefore always acceptable, but a low density may result in particle degeneracy. Choosing a waypoint at every event ensures that any weight variance induced at these sites is mitigated, but there is still the opportunity for weight variance to build up between events.

If ξ_x = π_x, particle weights diverge only because different particles $$ (x_{1:s}^{\left(i\right)},w_s^{\left(i\right)})$$ experience a different total intensity $$ r\left(x_s^{\left(i\right)}\right)$$ of observed events. If ESS₀ is the current estimated sample size, then under some assumptions, along an interval of length L where no events occur we have

To apply this to our situation, assume a panmictic population with constant diploid effective population size N_e. The variance of the total coalescent branch length in a sample of n individuals is $$ {\left(4N_e\right)}^2{\sum }_{i=1}^{n-1}{i}^{-2}$$ [54]. The variance of total mutation intensity σ² is obtained by multiplying this by μ², since the rate of mutations on the coalescent tree is μ times the total branch length. Rewriting this in terms of the heterozygosity θ = 4N_e μ, and approximating the sum with $$ {\sum }_{i=1}^{\infty }{i}^{-2}={\pi }^2/6$$ gives σ² = θ² π²/6, and a minimum waypoint distance of $$ 1/\sqrt{2{\sigma }^2}=\sqrt{3}/\pi \theta \approx 1/2\theta $$.

Because the assumptions mentioned above are in practice only met approximately, this minimum waypoint density should be taken as a guide; breakdown of the assumptions can be monitored by tracking the ESS, increasing the density of waypoints if necessary.

Parameter inference

Parameters can be inferred by stochastic expectation maximization (SEM), which involves maximizing the expected log likelihood over the posterior distribution of the latent variable. The probability density for a Poisson process is $$ \frac{1}{c!}{\theta }^c{e}^{-q\theta }$$, where c is the event count, and θ is the rate of events per unit of “opportunity” q, measured in units of time or space or some combination of them. The expected log likelihood c log θ − qθ (ignoring constants) is maximized for θ = c/q, where c and q are the expected event count and opportunity. We consider Markov jump processes X_s with parameters θ and distribution

To evaluate the expectations above we do not use the complete set of events in the full realisations x, since resampling causes early parts of x to become degenerate due to “coalescences” of the particle’s history of events along the sequence, which would lead to high variance of the estimates. Using only the most recent events is also problematic as these have not been shaped by many observations and mostly follow the prior π_x(x|θ), resulting in highly biased estimates. Smoothing techniques such as two-filter smoothing [55] cannot be used here since finite-time transition probabilities are intractable. For discrete-time models fixed-lag smoothing is often effective [39]. For our model the optimal lag depends on the epoch, as the age of tree nodes strongly influence their correlation distance. For each epoch we determine the correlation distance empirically, and for the lag we use this distance multiplied by a factor α; we obtain good results with α = 1.

Particularly in cases where some event types are rare, Variational Bayes can improve on EM by iteratively estimating posterior distributions rather than point estimates of θ. A tractable algorithm is obtained if the joint posterior π(x, θ|y)dxdθ is approximated as a product of two independent distributions over x and θ, and an appropriate prior over θ is chosen. For the Poisson example above, combining a Γ(θ|α₀, β₀) prior with the likelihood θ^c e^−qθ results in a Γ(θ|α₀ + c, β₀ + q) posterior. Similarly, with this choice the Variational Bayes approximation results in an inferred posterior distribution of the form

Explicitly, for model (1) the parameters $$ {\theta }^{\prime }=({\rho }_{\mathsf{\text{EM}}},C_{\mathsf{\text{EM}}})$$ maximising $$ \mathbb{E}\left[\text{log}{\pi }_x\left(x\right|{\theta }^{\prime })\right]$$, where the expectation is taken over the posterior x ∼ π(x|y, θ)dx as approximated by Algorithm 3, is

Simulation study

We implemented the model and algorithm above in a Python/C++ program SMCSMC (Sequential Monte Carlo for the Sequentially Markovian Coalescent) and assessed it on simulated and real data.

To investigate the effect of the lookahead particle filter, we simulated four 50 megabase (Mb) diploid genomes under a constant population-size model (N_e = 10, 000, μ = 2.5 × 10⁻⁸ and ρ = 10⁻⁸, both per generation and per site, generation time g = 30 years). We inferred population sizes N_e through evolutionary time, defined as the inverse of twice the instantaneous coalescent rate, as a piecewise constant function across 9 epochs (with boundaries at 400, 800, 1200, 2k, 4k, 8k, 20k, 40k and 60k generations) using particle filters Algorithms 2 and 3, as well as a recombination rate, which was taken to be constant through evolutionary time (and along the genome). Although recombination rate can be inferred, we here focus on the accuracy of the inferred N_e through evolutionary time. Observations are often available as unphased genotypes, and we assessed both algorithms using phased and unphased data, using the same simulations for both. Experiments were run for 15 EM iterations and repeated 15 times (Fig 2a).

Fig 2

Accuracy of population size inferences in simulated data.

Shown are true population sizes (black) and median inferred population sizes across 15 independent runs (blue); shaded areas denote quartiles and full extent. a Impact of lookahead, phasing and number of particles on the bias in population size estimates for recent epochs, for data simulated under a constant population size model. b Inference in the “zigzag” model on phased data using lookahead and 30, 000 particles, comparing inference using stochastic Expectation Maximization (SEM) and Variational Bayes (VB).

On phased data (Fig 2a, top rows), N_e values inferred without lookahead show a strong positive bias in recent epochs, corresponding to a negative bias in the inferred coalescence rate. Increasing the number of particles reduces this bias somewhat. By contrast, the lookahead filter shows no discernable bias on these data, even for as little as 1, 000 particles. On unphased data (Fig 2a, bottom rows), the default particle filter continues to work reasonably well; in fact the bias appears somewhat reduced compared to phased data analyses, presumably because integrating over the phase makes the likelihood surface smoother, reducing particle degeneracy. By contrast, the lookahead particle filter shows an increased bias on these data compared to the default implementation. This is presumably because of the reliance of the lookahead likelihood on the distance to the next singleton; this statistic is much less informative for unphased data, making the lookahead procedure less effective, and even counterproductive for early epochs.

We next investigated the impact of using Variational Bayes instead of stochastic EM, using the lookahead filter on phased data. We simulated four 2 gigabase (Gb) diploid genomes using human-like evolutionary parameters (μ = 1.25 × 10⁻⁸, ρ = 3.5 × 10⁻⁸, g = 29, N_e(0) = 14312) under a “zigzag” model similar to that used in [16] and [18], and inferred N_e across 37 approximately exponentially spaced epochs; see Appendix (“Implementation Details”). Both approaches give accurate N_e inferences from 2, 000 years up to 1 million years ago (Mya); other experiments indicate that population sizes can be inferred up to 10 Mya (but see Fig 3b). The upwards bias in the most recent epochs is reduced considerably by the Variational Bayes approach compared to SEM (Fig 2b), although some bias remains.

Fig 3

Population size inferences by SMCSMC on four diploid samples.

Left, three human populations (CEU, CHB, YRI), together with inferences from msmc using 1, 2 and 4 diploid samples. Right, simulated populations resembling CEU and YRI population histories. All inferences (SMCSMC, msmc) were run for 20 iterations.

Inference on human subpopulations

We applied SMCSMC to three sets of four phased diploid samples, of Northern European (CEU), Han Chinese (CHB) and Yoruban (YRI) ancestry respectively, from the 1000 Genomes project. For comparison we also ran msmc [16] inferring on the same data, and on subsets of 2 and 1 diploid samples. Inferences show good agreement where msmchas power (Fig 3). Since the inferences show some variation particularly in more recent epochs, we simulated data under a demographic model closely resembling CEU and YRI ancestry as inferred by multiple methods (see Appendix, “Implementation Details”), and we inferred population sizes using SMCSMC and msmc as before. This confirmed the accuracy of SMCSMC inferences from about 5,000 to 5 million years ago, while inferences in more recent epochs show more variability. A representative comparison of run times is provided in Table 1.

Table 1

Table lists means ± one standard deviation across 10 independent runs in a high performance compute environment. Note that due to parallel execution of SMCSMC (146 genomic chunks) and msmc (8 cores), wall clock time was considerably less than the total CPU time.

Runtimes (total CPU time, hours) for analyzing one or two diploid human genomes using msmc (40 EM iterations), and SMCSMC (15 Variational Bayes iterations).

Algorithm	2 haploids	4 haploids
msmc	5.2±0.5	107.3±18.7
SMCSMC 5,000 particles	109.2±5.7	277±15
SMCSMC 10,000 particles	219±11	673±32

Conditional distributions and the Markov property

Here we outline how to define a conditional distribution $$ \pi (·|\mathcal{G})$$ given a distribution π on $$ \mathcal{X}$$ and a conditioning subset $$ \mathcal{G}\subset \mathcal{X}$$ of measure 0. Suppose $$ {\mathcal{G}}_{\tau }$$ is a family of subsets of $$ \mathcal{X}$$ so that $$ {\cup }_{\tau }{\mathcal{G}}_{\tau }=\mathcal{X}$$. A particular subset $$ {\mathcal{G}}_{\tau }$$ for a fixed τ plays the role of the conditioning event B in the standard definition P(A|B) = P(A ∩ B)/P(B). It can be shown that, under some conditions, there exists an essentially unique family of measures $$ {\pi }_{{\mathcal{G}}_{\tau }}$$ and a measure μ so that $$ {\pi }_{{\mathcal{G}}_{\tau }}$$ is concentrated on $$ {\mathcal{G}}_{\tau }$$, $$ {\pi }_{{\mathcal{G}}_{\tau }}\left(\mathcal{X}\right)=1$$ for all τ, and $$ E_{\pi }\left[f\right]=\int \int f\left(x\right){\pi }_{{\mathcal{G}}_{\tau }}\left(\text{d}x\right)\mu \left(\text{d}\tau \right)$$ for well-behaved functions f [60], making it possible to define the conditional expectation as $$ E_{\pi }[f|\mathcal{G}]=E_{{\pi }_{\mathcal{G}}}[f]={\displaystyle \int f}(x){\pi }_{\mathcal{G}}(\text{d}x).$$ Using this, the Markov property of π can be expressed in terms of conditional expectations:

Proof of Algorithm 3

The algorithm is proved by induction on j. For j = 0 the loop invariant holds, while for j = K it implies the output condition. Suppose the loop invariant is true for some j. If ESS < N/2, assume w.l.o.g. that $$ v^{\left(i\right)}=v_{s_j}^{\left(i\right)}$$ are normalized, let i_k be the index of the kth new particle, $$ {\widehat{v}}^{\left(k\right)}=N^{-1}$$ and $$ {\widehat{w}}^{\left(k\right)}=N^{-1}{w}^{\left(i_k\right)}/v^{\left(i_k\right)}$$ be its weights, and write $$ {\pi }^j=\pi \left(X_{1:s_j}\right|Y_{1:s_j}=y_{1:s_j})$$, $$ {\tilde{\pi }}^j={\tilde{\pi }}^{s_j}\left(X_{1:s_j}\right|Y_{1:s_j}=y_{1:s_j})$$, then

After sampling $$ x_{s_j:s_{j+1}}^{\left(i\right)}\sim {\xi }_x\left(X_{s_j:s_{j+1}}\right|X_{s_j}=x_{s_j}^{\left(i\right)})$$, the particles $$ \left\{(x_{1:s_j}^{\left(i\right)},w_{s_j}^{\left(i\right)})\right\}$$ approximate $$ \pi \left(X_{1:s_j}\right|Y=y_{1:s_j}){\xi }_x\left(X_{s_j:s_{j+1}}\right|X_{s_j})$$. To make this distribution absolutely continuous w.r.t. $$ \pi (X_{1:s_{j+1}},Y_{s_j:s_{j+1}}|Y_{1:s_j})$$, multiply it with the constant measure $$ {\lambda }_{s_j:s_{j+1}}\left(y_{s_j:s_{j+1}}\right)$$; any measure will do as long as it has a density w.r.t. $$ {\lambda }_{s_j:s_{j+1}}$$ and is independent of X. Taking the Radon-Nikodym derivative of these two distributions gives

Particle weight variance

To derive a criterion on the waypoints that limits the effect of weight variance build-up, let R(s) = f(X_s) be the stochastic variable that measures the instantaneous rate of occurrence of emission events for a particular (random) particle X, and let $$ W\left(s\right)=W_0\text{exp}(-{\int }_0^{s}R\left(u\right)\text{d}u)$$ be that particle’s time-dependent weight; the dependence on W on X is not written explicitly. Note that the expression for W(s) is valid as long as no events have occurred in the interval [0, L). We assume that R(s) is time-homogeneous, that it can be approximated by a Gaussian process, that particles are drawn from the equilibrium distribution, and that W₀ and R(s) are independent. Write 〈V(X)〉 ≔ ∫V(X)dπ(X) for the expectation of V over π(X). Writing $$ R\left(s\right)=\mu +\tilde{R}\left(s\right)$$ where μ = 〈R(s)〉 is the mean event rate (which is independent of s by assumption), then

In practice particles will not be drawn from the equilibrium distribution π_x(X), but from the joint distribution on X and Y conditioned on observations y. However, for most likelihoods conditioning will reduce the variance of R as observations tend to constrain the distribution of likely particles, making this a conservative assumption. The other assumption that is likely not met is that R(t) is a Gaussian process; it is less clear whether making this approximation will in practice be conservative.

The sequential coalescent with recombination process

In formula (1), if s is a recombination point, x_s is the genealogy just left of the recombination point and includes the infinite branch from the root, so that b_u(x_s) = 1 for u above the root.

The measure (1) describes the CwR process exactly as long as x encodes both the local genealogy and the non-local branches used by the SCRM algorithm. In practice the SCRM algorithm prunes some of these branches, and we use (1) on the pruned x.

Note that we take the view that the realisation x encodes not only the sequence of genealogies x_s but also the number of recombinations |x| (some of which may not change the tree), their loci $$ s_j=s_j^x$$, and the recombination and coalescence times $$ {\nu }_j^x$$ and $$ {\tau }_j^x$$. This information is also kept in the implementation of the algorithm, and is used to calculate the sufficient statistics required for inference of the coalescence and recombination rates.

Variational Bayes for Markov jump processes

We consider hidden Markov models where the latent variable follows a Markov jump process over $$ x\in \mathcal{X}$$, that with respect to a suitable measure dxdy admits a probability density of the form

A Variational Bayes approach approximates the true joint posterior density π(x, θ|y) ∝ π_xy(x, y|θ)π_θ(θ), where π_θ is a prior on the parameters, with a probability density ϕ(x, θ) that is easier to work with (here the constant of proportionality implied by “∝” hides a constant density λ(y)). Following Hinton and van Camp [61] and MacKay [62], we choose to constrain ϕ by requiring it to factorize as ϕ(x, θ) = ϕ_x(x)ϕ_θ(θ), and we choose to optimize it by minimizing the Kullback-Leibler divergence KL(ϕ||π), also referred to as the variational free energy [63],

A lookahead likelihood

Let s_i denote the distance along the genome to the nearest future singleton in each sequence, and let c_k = (a_k, b_k) be ≤ n/2 mutually consistent cherries with loci $$ s_k^{\prime }\le s_k^{″}$$ of their first and last supporting doubleton. To simplify notation we assume that the current locus is 0 (Fig 1a).

Note that recombinations result in a change of a terminal branch length (TBL) if either the recombination occurred in the branch itself and the new lineage does not coalesce back into it, or the recombination occurred outside the branch and the new lineage coalesces into it (Fig 1b). To compute the likelihood that the first singleton in lineage i occurs at locus s_i, we assume that all TBLs are equal to l_i, and that coalescences occur before l_i. Then, the total rate of events that change the TBL i is

To approximate the likelihood of the doubleton data, note that a node c with precisely two descendants (a, b) (a “cherry”) at height l changes if a recombination occurs in either branch a or b and the new lineage coalesces out, or a recombination occurs outside of a and b and coalesces into either (Fig 1b). Again assuming that all TBLs are l and coalescences occur before l, the total rate of change is $$ 2l\rho \frac{n-2}{n}+(n-2)l\rho \frac{2}{n}=4l\rho \frac{n-2}{n}:={\rho }_C$$. When a cherry changes, we assume that the new cherry is drawn from the equilibrium distribution. To calculate the probablity of observing c = (a, b) at equilibrium, assume that a tree supports 1 ≤ k ≤ n/2 cherries. The branches of c are among the 2k branches subtended by the tree’s k cherries with probability $$ \frac{2k}{n}\frac{2k-1}{n-1}$$, and a is paired with b with probability $$ \frac{1}{2k-1}$$. Since k has mean n/3 if n ≥ 3 [64], the probability of observing (a, b) at equilibrium is $$ \frac{2}{3(n-1)}$$. We approximate the likelihood of a doubleton by 0 if the c is not in the tree, and by 1 if it is. Then, the likelihood of observing c_k = (a_k, b_k) at the last known locus $$ s_k^{″}$$ conditional on the tree currently containing c_k is

These likelihoods show good performance, but result in some negative bias in inferred population size for recent epochs. We traced this to the lack of correlation between l_i and $$ l_i^{\prime }$$, requiring a single very recent coalescence to explain a long segment devoid of singletons, rather than allowing for the possibility of several correlated coalescences each in slightly earlier epochs. To model correlations, we averaged the likelihood above over ρ′ = ρ and ρ′ = ρ/2 each weighted with probability 1/2. This effectively removed the negative bias.

To deal with missing data, we reduce μ proportionally to the missing segment length and the number of lineages missing. For unphased mutation data, singletons and doubletons can still be extracted, and are greedily assigned to compatible lineages. The likelihoods are also similarly calculated, by greedily assigning cherries to observed doubletons. Unphased singletons can result from mutations on either of the individual’s alleles; the likelihood term uses the sum of the two branch lengths for that individual to calculate the expected rate of unphased singletons.

Implementation details

While x_1:s refers to the entire sequence of genealogies along the sequence segment 1: s, storing this sequence would require too much memory. Instead we only store the most recent genealogy x_s (including non-local branches where appropriate), which is sufficient to simulate subsequent genealogies using the SCRM algorithm. To implement epoch-dependent lags when harvesting sufficient statistics, we do store records of the events (recombinations, coalescences and migrations) that changed x along the sequence, as well as the associated opportunities, for each particle and each epoch; this implicitly stores the full ARG. To avoid making copies of potentially many event records when particles are duplicated at resampling, these are stored in a linked list, and are shared by duplicated particles where appropriate, forming a tree structure. Records are removed dynamically after contributing to the summary statistics, and when particles fail to be resampled, ensuring that memory usage is bounded.

The likelihood calculations involve many evaluations of the exponential function, often for small exponents. We use the continued-fraction approximation $$ e^x\approx 1+2x/(2-x+\frac{1}{6}{x}^2)$$ for |x| < 0.03, with relative error bounded by 10⁻¹⁰ [65].

Table 2 shows the commands to generate the data for the three simulation experiments. Epoch boundaries for N_e inference in generations for the zigzag experiment were defined by taking interval boundaries −14312log(1−i/256)/2, i = 0, …, 255, merging intervals according to the pattern 4 * 1 + 7 * 2 + 8 * 5 + 7 * 13 + 1 * 15 + 8 * 11 + 1 * 3 (37 epochs; see [14]). For the real data experiments, epochs boundaries for the 32 epochs were logarithmically spaced from 133 to 133016 generations ago, using generation time g = 29 years, without merging intervals (command line option -P 133 133016 31*1).

Table 2

Commands to generate simulation data.

Experiment	Command
zigzag	scrm 8 1 -N0 14312 -t 1431200 -r 400736 2000000000 -eN 0 1 -eG 0.000582262 1318.18 -eG 0.00232905 -329.546 -eG 0.00931619 82.3865 -eG 0.0372648 -20.5966 -eG 0.149059 5.14916 -eN 0.596236 0.1 -seed 1 -T -L -p 10 -l 300000
CEU	scrm 8 1 -N0 14312 -t 1789000 -r 500920 2500000000 -eN 0 10.4807 -eG 0.00120468 214.8965 -eG 0.0180702 -14.15827 -eG 0.180702 1.33255 -eG 1.084212 -0.563414 -eN 2.71053 2.096143 -seed 1 -T -L -p 10 -l 300000
YRI	scrm 8 1 -N0 14312 -t 1789000 -r 500920 2500000000 -eN 0 10.4807 -eG 0.00120468 502.8635 -eG 0.00542106 0 -eG 0.0451755 -5.89189 -eG 0.180702 1.33255 -eG 1.084212 -0.563414 -eN 2.71053 2.096143 -seed 1 -T -L -p 10 -l 300000