PaperlayerView exact revision on arXiv

Source-linked AI summary

Bayesian inference of sampled ancestor trees for epidemiology and fossil calibration

Alexandra Gavryushkina, David Welch, Tanja Stadler, Alexei Drummond

arXiv:1406.4573v3q-bio.PE

TL;DR

Standard phylogenetic models assume samples are leaves, limiting analyses when sampled individuals can be ancestors. The paper develops a Bayesian MCMC framework for sampled ancestor trees and applies it to epidemiological and fossil-calibrated data. It identifies sampled ancestors, supports divergence-time inference, and establishes parameter non-identifiability in some transmission models.

  • Problem

    Standard phylogenetic models assume all samples are terminal nodes, preventing direct-ancestor relationships in serially sampled infections and fossil-inclusive phylogenies.

  • Method

    The paper implements reversible-jump Bayesian MCMC in BEAST2 and extends birth-death skyline models to sampled ancestor trees.

  • Results

    The method identifies sampled ancestors in HIV data and infers divergence times from fossil and recent bear taxa while accounting for sampled-ancestor relationships.

  • Takeaways & Limitations

    Sampled ancestor models support infectious-transmission and fossil-calibrated macroevolutionary analyses within a unified Bayesian framework.

  • Takeaways & Limitations

    Transmission birth-death models with sampled ancestors require one parameter to be fixed or strongly constrained because their parameters are non-identifiable.

Abstract

from arXiv · show

Phylogenetic analyses which include fossils or molecular sequences that are sampled through time require models that allow one sample to be a direct ancestor of another sample. As previously available phylogenetic inference tools assume that all samples are tips, they do not allow for this possibility. We have developed and implemented a Bayesian Markov Chain Monte Carlo (MCMC) algorithm to infer what we call sampled ancestor trees, that is, trees in which sampled individuals can be direct ancestors of other sampled individuals. We use a family of birth-death models where individuals may remain in the tree process after the sampling, in particular we extend the birth-death skyline model [Stadler et al, 2013] to sampled ancestor trees. This method allows the detection of sampled ancestors as well as estimation of the probability that an individual will be removed from the process when it is sampled. We show that sampled ancestor birth-death models where all samples come from different time points are non-identifiable and thus require one parameter to be known in order to infer other parameters. We apply this method to epidemiological data, where the possibility of sampled ancestors enables us to identify individuals that infected other individuals after being sampled and to infer fundamental epidemiological parameters. We also apply the method to infer divergence times and diversification rates when fossils are included among the species samples, so that fossilisation events are modelled as a part of the tree branching process. Such modelling has many advantages as argued in literature. The sampler is available as an open-source BEAST2 package (https://github.com/gavryushkina/sampled-ancestors).

Introduction

The paper addresses phylogenetic settings where sampled individuals may be direct ancestors, including serially sampled infections and fossil taxa. It develops a Bayesian MCMC framework for sampled ancestor trees and applies it to HIV and bear data.

  • Motivation: Standard phylogenetic models treat all samples as terminal nodes, which is inappropriate when sampled individuals can remain infectious or ancestral.This issue arises in serially sampled epidemiological and fossil data.
  • Motivation: Birth-death-sampling models allow incomplete sampling but lacked software support for infection after sampling, so applications often ignored sampled ancestors.
  • Motivation: Fossilisation events can be modelled within the tree-process prior to jointly analyse fossil and recent taxa while allowing ancestor-descendant relationships.
  • Computational challenge: Sampled ancestor models require extensions to standard MCMC because sampled nodes on branches create non-binary trees with variable dimensionality.
  • Contribution: The paper implements a reversible-jump MCMC kernel in BEAST2, extends the birth-death skyline model, and analyses serially sampled HIV data and fossil-bearing bear sequences.

Methods

The method models incomplete, serially sampled birth-death processes in which sampled individuals may remain in the process and become sampled ancestors. A piecewise-constant skyline extension covers changing birth, death, sampling, and removal parameters.

  • Transmission model: The transmission birth-death process uses birth rate λ, extinction rate µ, sampling rate ψ, removal probability r, and origin time t_or.Lineages bifurcate, go extinct, or are sampled through time.
  • Sampled ancestor trees: Sampling without removal produces degree-two sampled ancestors, while the reconstructed tree has no origin node because origin time is a model parameter.
  • Tree representation: The genealogy combines a ranked labeled topology with a time vector containing bifurcation, tip, and sampled two-degree-node times.
  • Likelihood: The model likelihood is conditioned on observing at least one sampled individual and is formulated for non-oriented labeled trees.
  • Identifiability: The transmission model is unidentifiable because its likelihood depends on λ−µ−ψ, λψ, and ψ(1−r), rather than λ, µ, ψ, and r independently.
  • Fossilized birth-death process: The fossilized birth-death process sets r to zero and permits present-day sampling with probability ρ, allowing all four parameters λ, µ, ψ, and ρ to be identified.
  • Skyline extension: The skyline extension uses piecewise-constant parameters across time intervals and contains skyline transmission and fossilized birth-death models as special cases.

Markov chain Monte Carlo Operators

The MCMC implementation uses operators that move among sampled ancestor trees while preserving valid topology and timing constraints. Reversible-jump moves additionally change tree dimension by converting sampled ancestors into leaves or vice versa.

  • Operator framework: The operator suite explores sampled ancestor trees with a fixed number of sampled nodes using node heights as continuous state variables.
  • Wilson Balding: The extended Wilson Balding operator prunes and regrafts subtrees, can change the root, and uses reversible-jump updates when tree dimension changes.
  • Wilson Balding: Wilson Balding selects a non-root edge, chooses an eligible older edge or leaf, and proposes a new parent height when attaching to a branch.
  • Hastings ratios: The proposal density combines edge-selection and height-selection probabilities, while shared selection terms cancel in the forward-backward Hastings ratio.
  • Dimension-changing move: A dimension-changing move inserts a node to turn a sampled ancestor into a leaf or replaces a leaf with a sampled ancestor when topology permits.
  • Exchange and height operators: Narrow and wide exchange operators are extended to sampled ancestor trees, with scaling and uniform height moves contributing Hastings ratios α^(k−2) and 1, respectively.

Simulations and empirical data analysis

The authors evaluate parameter and genealogy recovery through simulations and apply the framework to bear fossils and UK HIV-1 data. Analyses use birth-death parameterisations, sequence evolution models, and MCMC estimation of trees and evolutionary parameters.

  • Fossilized birth-death simulations: The fossilized birth-death simulation fixes model parameters, generates trees and GTR sequences, and estimates tree parameters, clock rate, and sampled-ancestor features.
  • Transmission simulations: The transmission simulations draw parameters from broad uniform distributions, discard trees outside the sampled-node range, and report 100 retained trees averaging 53 sampled nodes.
  • Transmission inference: Only three of four birth-death parameters are inferred in transmission-process MCMC, so the fossilisation proportion s is fixed to its true value.
  • Skyline simulations: Skyline transmission simulations vary parameters across two or three intervals while fixing either r, ψ, or the full vector of removal probabilities.
  • Bear data: The bear reanalysis compares BEAST2 and DPPDiv under the same fossilized birth-death model, using a fixed extant topology and strict molecular clocks.
  • HIV data: The UK HIV analysis uses a skyline model without ρ-sampling, one rate shift in 1999, and priors on reproductive number, removal rate, sampling proportion, and removal probability.

Results

The sampled ancestor framework recovered trees and parameters in simulations, supported accurate fossil and HIV analyses, and identified sampled ancestors in epidemiological data.

  • The framework was implemented in BEAST2 and applied to simulations, fossil-calibrated divergence dating, and an HIV dataset.The fossil-bear analysis was also compared with an alternative implementation.
  • Simulation of sampled ancestor models: Simulation studies recovered trees and model parameters from sequence data and sampling times across sampled ancestor birth-death scenarios.The studies covered both sampled ancestor birth-death and skyline processes.
  • Simulation of sampled ancestor models: Fixing one tree-model parameter enabled accurate recovery of the remaining parameters in otherwise non-identifiable model variants.For example, fixing ψ allowed recovery of λ, µ, and r; other fixed-parameter scenarios also produced accurate estimates.
  • Simulation of sampled ancestor models: Fossilized birth-death simulations had worst-case median relative parameter errors of 0.22, while tree-property errors were at most 0.09.True parameters and tree properties fell within 95% HPD intervals at least 95% of the time.
  • Simulation of sampled ancestor models: Transmission birth-death simulations had maximum median relative errors of 0.28 for parameters and 0.06 for tree properties.In the worst case, a parameter or tree property was inside the 95% HPD interval 92% of the time.
  • Application of sampled ancestor Skyline model to HIV dataset: In the HIV-1 dataset, three sampled nodes had posterior sampled-ancestor probabilities of 61%, 59%, and 49%, with Bayes factors of 5.9, 8.7, and 4.2.Other sampled nodes had probabilities below 4%.

Discussion

The sampler supports sampled-ancestor analyses for epidemiological and fossil data, while exposing both practical benefits and identifiable limitations. Its applications include uncertainty-aware fossil phylogenies, direct-ancestor detection, and inference of removal-at-sampling parameters, but transmission models may require parameter constraints.

  • Fossil data: The sampler can analyze fossil and recent taxa jointly to infer divergence times while integrating over fossil attachment times and topological attachment points.Unlike an approach requiring a known extant topology, this implementation can account for uncertainty in the extant phylogeny.
  • Fossil data: When extant phylogeny uncertainty is present, the sampler can account for it, although analytical calculation produces faster mixing when the topology is well resolved.
  • Epidemiological data: Simulation studies show that sampled ancestors can be detected, including sampled individuals that later infected other individuals in epidemiological analyses.
  • Epidemiological data: The removal-at-sampling parameter r can be estimated and often reflects the probability that diagnosed patients remain able to cause further infections.
  • Model limitations: Transmission birth-death models including r are non-identifiable unless one parameter is fixed or strongly constrained by prior information.With fossils, all four tree-process parameters λ, µ, ρ, and ψ can be estimated when time-stamped comparative data are available.
  • Broader significance: The implementation is presented as a first full MCMC sampler for sampled ancestor trees and as a computational basis for future fossil dating and phylodynamics.

Figures

The figures illustrate sampled ancestor trees, reversible-jump tree proposals, simulation estimates, ancestor-identification performance, and applications to bear and HIV data.

  • Figure 1: Figure 1 contrasts a full sampled ancestor birth-death tree with its reconstructed tree, including sampled ancestors and skyline parameter intervals.Nodes A, B, and D are sampled ancestors; skyline parameters change between t0–t1 and t1–t2, with additional sampling attempts at t1 and t2.
  • Figure 2: The extended Wilson Balding operator prunes and reattaches subtrees while allowing sampled ancestor trees to gain or lose nodes.Pruning from a branch followed by attachment to a leaf removes a node, whereas pruning from a node followed by attachment to an edge introduces one.
  • Figure 3: Figure 3 compares median estimates and 95% HPD intervals with true values for tree height and the number of sampled ancestors.The fossilized birth-death simulations evaluate both overall tree height and sampled-ancestor counts.
  • Figure 4: Figure 4 plots relative 95% HPD interval widths for turnover rate ν against tree size in fossilized birth-death simulations.The plot focuses on uncertainty in turnover-rate estimates as simulated trees become larger.
  • Figure 5: Figure 5 compares median estimates and 95% HPD intervals with true values for turnover rate ν and removal probability r in transmission simulations.The two panels show the turnover-rate and removal-probability parameters separately.
  • Figure 6: Figure 6 evaluates sampled-ancestor classification across posterior-probability thresholds using ROC sensitivity-specificity trade-offs, with an optimal threshold of 0.45.The dashed diagonal represents random guessing; curves nearer the upper-left indicate more accurate tests.
  • Figure 7: Figure 7 shows bear divergence-time estimates from DPPDiv and BEAST2 fossilized birth-death analyses, which give the same results.Bars represent 95% HPD intervals and dots represent means for nodes spanning bear outgroups, the living-bear clade, and within-clade divergences.

Tables

The table summarizes Hastings ratios for the extended Wilson Balding operator, whose proposals can alter sampled ancestor tree topology and dimension.

  • Table 1: Table 1 summarizes the Hastings ratio q(g*|g) and q(g|g*) for cases of the extended Wilson Balding operator.The operator handles pruning from branches or nodes and attachment to branches or leaves.

1 Sampled ancestor Skyline model

The sampled ancestor skyline model derives tree densities for birth-death processes with sampled nodes that may remain in the process. Its re-parameterisation reveals that, under specified sampling conditions, one original parameter cannot be identified from the sampled tree.

  • Sampled ancestor skyline process: The SABD skyline process assigns probability densities to reconstructed trees with branching, sampling, and sampled-ancestor events.The derivation uses edge-specific probabilities and traverses the tree from tips to the root.
  • Tree density derivation: The tree density can be expressed through oriented trees and, after accounting for labelings and orientations, yields the density for labeled genealogies.With fixed parameters, the oriented-tree density depends on event times and node counts, not on lineage connectivity.
  • Re-parameterisation: When intermediate sampling proportions are zero, the model can be re-parameterised with 4l parameters instead of the original 4l + 1.The re-parameterisation uses combinations involving d_i, f_i, g_i, h, and k_i.
  • Identifiability: The sampled-ancestor tree likelihood then depends only on the re-parameterised quantities, so one parameter is not identifiable from the sampled tree.This non-identifiability is established for the tree likelihood conditioned on at least one ψ-sampled individual.
  • Fossilised birth-death special case: For the skyline fossilised birth-death case with r = 0 and ρ_l ≠ 0, the re-parameterisation does not reduce the number of parameters.The model has 3l + 1 initial parameters but 4l re-parameterisation parameters.

2 Testing Operators

The authors test sampled-ancestor MCMC operators by comparing sampled tree distributions with analytically calculated probabilities. They evaluate topology sampling using repeated runs and effective sample size.

  • Implementation: The BEAST2 sampled-ancestor add-on implements operators for random walks in sampled ancestor tree space.The implementation is tested against calculations performed in Mathematica.
  • Topology validation: The topology test considers eight non-ranked tree topologies and compares estimated marginal probabilities with exact Mathematica probabilities.Agreement is assessed using standard errors and whether estimates fall within two standard errors of the true values.
  • MCMC diagnostics: Effective sample size is calculated by assigning integers to the eight topologies and applying ESS estimation to the resulting integer sample.The study uses 100 MCMC runs to test operators that preserve sampled-node times.

3 Simulation studies

Simulation studies evaluate sampled-ancestor inference under fixed and prior-drawn parameters, with either trees alone or sequences and dates as data. Performance is summarised using posterior errors, biases, interval widths, and coverage.

  • Simulation design: The simulations include tree-only and sequence-based scenarios for estimating tree parameters, trees, and molecular parameters.Scenario 1 fixes trees during MCMC, whereas Scenario 2 jointly estimates trees and parameters from sequences and sampled-node dates.
  • Simulation design: The study simulates 100 trees in most scenarios, varying whether parameters are fixed or drawn from priors and whether simulations stop by sample size or time of origin.Different scenarios use fixed sampled-node counts, fixed time intervals, or parameter draws from uniform priors.
  • Identifiability assumptions: For models without ρ-sampling, one parameter is fixed to its true value in MCMC because not all parameters can be inferred.The time of origin receives a uniform prior on [0, 1000] in all scenarios.
  • Sequence simulation: Sequences are simulated at 2000 bp under GTR with fixed rates and frequencies, using a strict molecular clock with a fixed substitution rate.The specified substitution-rate parameters are listed explicitly in the simulation setup.
  • Evaluation criteria: Parameter performance is evaluated using posterior medians, errors, relative biases, relative 95% HPD widths, and 95% HPD accuracy.Across 100 runs, the reported summaries are medians of the corresponding statistics, with HPD accuracy counting inclusion of the true value.

4 HIV-1 data analysis

The section notes that some taxon names in Figure 8M are accompanied by accession numbers in a table.

  • Figure 8M: Accession numbers are provided in a table for some taxon names shown in Figure 8M.
  • Figure 8M: The table supplements Figure 8M by linking selected taxon names to accession numbers.
  • Figure 8M: Figure 8M should be read together with the accession-number table for the taxa listed there.
Loading 1406.4573v3…