PaperlayerView exact revision on arXiv

Source-linked AI summary

Genome-wide inference of ancestral recombination graphs

Matthew D. Rasmussen, Melissa J. Hubisz, Ilan Gronau, Adam Siepel

arXiv:1306.5110v3q-bio.PE

TL;DR

Existing ARG inference methods have not been fast and accurate enough for large-scale genomic sequence data. The paper introduces ARGweaver, a method with improved scaling that can support analyses of long sequences and provide insights into demographic processes and natural selection.

  • Problem

    Existing ARG inference methods have not been fast and accurate enough for practical use with large-scale genomic sequence data.

  • Method

    ARGweaver introduces an ARG inference method faithful to the sequentially Markov coalescent with improved scaling properties over previous methods.

  • Results

    ARGweaver can sample full sequences megabases in length and provide new insights into demographic processes and natural selection.

  • Takeaways & Limitations

    The method supports using ARG inference to study demographic processes and natural selection in large-scale genomic data.

  • Takeaways & Limitations

    Different methods have somewhat different strengths and weaknesses depending on the application.

Abstract

from arXiv · show

The complex correlation structure of a collection of orthologous DNA sequences is uniquely captured by the "ancestral recombination graph" (ARG), a complete record of coalescence and recombination events in the history of the sample. However, existing methods for ARG inference are computationally intensive, highly approximate, or limited to small numbers of sequences, and, as a consequence, explicit ARG inference is rarely used in applied population genomics. Here, we introduce a new algorithm for ARG inference that is efficient enough to apply to dozens of complete mammalian genomes. The key idea of our approach is to sample an ARG of n chromosomes conditional on an ARG of n-1 chromosomes, an operation we call "threading." Using techniques based on hidden Markov models, we can perform this threading operation exactly, up to the assumptions of the sequentially Markov coalescent and a discretization of time. An extension allows for threading of subtrees instead of individual sequences. Repeated application of these threading operations results in highly efficient Markov chain Monte Carlo samplers for ARGs. We have implemented these methods in a computer program called ARGweaver. Experiments with simulated data indicate that ARGweaver converges rapidly to the true posterior distribution and is effective in recovering various features of the ARG for dozens of sequences generated under realistic parameters for human populations. In applications of ARGweaver to 54 human genome sequences from Complete Genomics, we find clear signatures of natural selection, including regions of unusually ancient ancestry associated with balancing selection and reductions in allele age in sites under directional selection. Preliminary results also indicate that our methods can be used to gain insight into complex features of human population structure, even with a noninformative prior distribution.

Author Summary

The ARG captures the complex correlations among orthologous sequences, but previous inference methods were not practical for large-scale genomic data. The authors introduce ARGweaver, which scales ARG inference to long sequences and entire mammalian genomes, enabling analyses of demography and natural selection.

  • Motivation: An ARG records coalescence and recombination events, representing the changing genealogies and correlation structure of sampled sequences.This structure underlies analyses of population sizes, divergence, gene flow, selection, recombination, and genotype/phenotype associations.
  • Motivation: Previous ARG-inference methods were not fast and accurate enough for practical use with large-scale genomic sequence data.The resulting computational limitation restricted explicit ARG inference in applied population genomics.
  • Method and scope: ARGweaver has improved scaling properties and can analyze sequences megabases in length.The program implements the new algorithm and supports sampling full ARGs for entire mammalian genome sequences with a large computer cluster.
  • Results: ARGweaver performs well in simulation experiments and provides insights into demographic processes and natural selection from real human genome data.The authors report applications to human sequences that reveal information about reduced diversity near functional elements, ancient polymorphisms, and selection-related allele ages.

Introduction

The introduction frames the ARG as a complete but difficult-to-infer representation of genome-wide genealogical relationships. It motivates a new SMC-faithful algorithm that improves scaling while retaining access to detailed ARG features relevant to demography and selection.

  • ARG representation: An ARG records coalescence and recombination events and specifies a complete genealogy at each genomic position.It is therefore richer than commonly used low-dimensional summaries such as the site frequency spectrum, principal components, haplotype maps, and identity-by-descent spectra.
  • ARG representation: Accurate ARGs could simplify estimation of recombination rates, ancestral population sizes, divergence times, gene flow, and selective sweeps.The ARG is presented as an information-rich data structure for population genomic analysis.
  • Inference challenge: Existing explicit ARG methods are computationally intensive, approximate, or limited to small samples, leaving them rarely used in applied population genomics.SMC-based HMM methods require complete characterization of the state space, restricting them to small numbers of samples and reducing power for recent demographic and selection-related features.
  • New approach: The authors introduce an SMC-faithful ARG-inference algorithm with improved scaling properties based on sampling an ARG of n chromosomes conditional on an ARG of n−1 chromosomes.The approach combines the scalability of large-sample methods with explicit ARG characterization.
  • New approach: Explicit ARG characterization enables examination of local coalescence times, allele ages, and gene-tree topologies.These quantities are used to study demographic processes and natural selection across the genome.
  • Applications: Applied to 54 human sequences, ARGweaver provides insight into reduced diversity near functional elements, balancing selection, and direct versus indirect selection on allele age.The authors also report promise for addressing questions about human population structure and selection.

The Sequentially Markov Coalescent

The SMC approximates the coalescent-with-recombination as a first-order Markov process along the genome, generating local trees and recombination breakpoints. The model can be extended with sequence data to generate complete ARGs and aligned DNA sequences.

  • Model definition: The SMC is the starting model for generating a sequence of local genealogical trees and recombination events along a genomic segment.It approximates the full coalescent-with-recombination distribution by simplifying dependencies between neighboring local trees and breakpoints.
  • Discretization: The DSMC discretizes the continuous model, using exhaustive nonoverlapping genomic intervals and a discretized treatment of coalescence and recombination events.This formulation generalizes the pairwise discretized sequentially Markov coalescent to multiple genomes.
  • Model representation: The model represents each nonrecombining interval with a local tree and each breakpoint between intervals with a recombination event.Trees have node ages measured in generations ago, while breakpoints lie along the continuous genomic segment.
  • Generative process: The SMC generates the initial tree from the standard coalescent, then iteratively samples breakpoints, recombination points, and subsequent tree structures.The process repeats until the genomic segment is covered.
  • Generative process: Retaining sampled recombination points turns the sequence of local trees into a complete ARG.The local trees can also generate aligned DNA sequences corresponding to their leaves, making the augmented SMC a generative model for ARGs and sequence data.

The Discretized Sequentially Markov Coalescent

The DSMC approximates the sequentially Markov coalescent by discretizing genomic positions and time, enabling ARGs and sequence alignments to be modeled across multiple genomes.

  • Model definition: The DSMC is a discrete-in-space-and-time approximation of the sequentially Markov coalescent.It generalizes the model to multiple genomes and approaches the SMC as spatial and temporal resolution increase.
  • Time discretization: Time is partitioned into K intervals with events assigned to K + 1 discrete time points, distributed uniformly on a logarithmic scale.The final time point is user-specified, and the discretization is finest near the ARG leaves.
  • ARG representation: The model represents an ARG using local trees T_n and recombination events R_n at individual nucleotide positions.Each local block has an integral length, and recombinations occur between adjacent nucleotides with at most one recombination per pair of sites.
  • Assumptions: The model assumes at most one recombination between adjacent sites, a simplification expected to have at most a minor effect during inference.The stated effect depends on the sparsity of variant sites in most datasets.
  • Sequence generation: The DSMC can generate aligned DNA sequences by sampling ancestral alleles and mutating them along branches of the corresponding local trees.This work uses a Jukes-Cantor mutation model, although other mutation models can be substituted.
  • Sequence likelihood: Although recombination events define the ARG, the likelihood of sequence data given the ARG depends only on the local trees.This separates the complete ARG representation from the tree-dependent sequence-data likelihood.

The Threading Problem

Threading samples the ancestry of one additional sequence conditional on an ARG for the other sequences. The conditional problem becomes an HMM that can be solved exactly under the DSMC assumptions.

  • Motivation: Direct HMM inference is impractical because the number of possible local trees, and therefore the state space, grows super-exponentially with sequence number.Earlier approaches consequently accommodated only small numbers of sequences.
  • Conditional sampling: Threading samples a new ARG G_n conditional on G_n−1, the combined alignment D_n, and model parameters Θ under the DSMC.The new ARG augments the existing graph with the additional coalescence and recombination events for the nth sequence.
  • Two-step threading: The new sequence is represented by coalescence threading Y and recombination threading Z, sampled in two conditional steps.First Y is sampled while integrating over new recombination locations; then Z is sampled conditional on Y.
  • HMM formulation: The coalescence-threading distribution forms an HMM whose states are possible attachment values y_i along the sequence.The locations of new recombinations are integrated into the transition probabilities.
  • Inference: Dynamic programming with stochastic HMM traceback samples the coalescence threading exactly under the DSMC.The model retains Markovian dependencies despite heterogeneous state spaces and normalization structure.

Markov chain Monte Carlo sampling

Threading provides building blocks for ARG samplers, including sequential initialization, single-sequence Gibbs sampling, and subtree-based Metropolis-Hastings updates. Subtree threading addresses poor mixing as sample size grows.

  • Sampler framework: Repeated threading operations form Markov chain Monte Carlo samplers that explore the posterior distribution of ARGs.The threading operation can be used as a building block for full-ARG sampling.
  • Sequential sampling: Sequential sampling constructs an initial ARG by starting with one sequence and threading additional sequences one at a time.The resulting ARG is not a valid posterior sample because each intermediate graph conditions only on the sequences added so far.
  • Single-sequence Gibbs sampling: Single-sequence Gibbs sampling removes and rethreads sequences, producing a valid Gibbs sampler under the DSMC assumptions.Sequences may be rethreaded randomly or in round-robin fashion, and the Markov chain is ergodic under the stated reachability argument.
  • Mixing limitation: Single-sequence rethreading mixes poorly for large samples because it changes external branches inefficiently and rarely rearranges deep internal structure.The mixing problem worsens as the number of sequences grows.
  • Subtree sampling: Subtree threading resamples internal branches and perturbs the deep structure of the ARG by threading subtrees with arbitrary numbers of leaves.A branch graph and dynamic programming identify subtree paths with good continuity properties for efficient sampling.
  • Subtree sampling: ∼40% acceptance for typical human data accompanies substantially improved mixing from generalized subtree threadings.Because the number of removable internal branches varies with ARG structure, the method uses Metropolis-Hastings rather than a unit Hastings ratio.

ARGweaver Program and Visualization

ARGweaver implements threading-based samplers and summarizes posterior ARG samples with genome-browser statistics and leaf traces. Leaf traces make recombination patterns visible through line-order and spacing changes.

  • Program: ARGweaver implements single-sequence Gibbs sampling, subtree Metropolis-Hastings sampling, and sequential sampling for initial ARG construction.The program supports configurable sampling iterations and recording frequency.
  • Performance: A single 1 Mb threading in an ARG of 20 sequences with 20 time steps takes about 1 second.The program is written in C++ and Python and is described as reasonably well optimized.
  • Posterior summaries: Posterior samples are summarized using TMRCA, total branch length, and derived-allele age at each genomic position, with posterior means and 95% credible intervals.These statistics can be visualized in genome-browser tracks alongside other genomic information.
  • Leaf traces: A leaf trace draws one line per haploid sequence, ordering lines by local genealogy and spacing them proportionally to their TMRCAs.Line order and spacing change where recombinations occur.
  • Leaf traces: Recombination hotspots appear as dense clusters of vertical lines, whereas cold spots appear as long blocks of parallel lines.These visual patterns provide an immediate comparison of recombination intensity along the genome.

Simulation Study

ARGweaver converged rapidly and recovered global and local ARG features across simulated datasets, although performance depended on sequence count and showed model-related bias at low µ/ρ.

  • The single-sequence Gibbs sampler converged rapidly for small samples but became inadequate beyond roughly 6–8 sequences.
  • Subtree threading with Metropolis-Hastings addressed the larger-sample convergence problem, allowing the number of sequences to increase to 20 or more.
  • 500 iterations brought 20-sequence, 1-Mb analyses within about 500?
  • ARGweaver recovered simulated global ARG features with fairly high accuracy, while estimate variance decreased as µ/ρ increased from 1 to 6.
  • At low µ/ρ, recombination counts were slightly underestimated and joint probabilities overestimated, consistent with DSMC model misspecification.
  • ARGweaver accurately recovered TMRCA and local-tree topology, detected simulated recombination hotspots with 100 sequences, and estimated allele ages better than allele frequency alone.
  • Its posterior uncertainty for local trees closely matched actual correctness, whereas Margarita often overestimated branch confidence by 10–20%.

Analysis of Real Data

Analyses of human genomes used ARG-derived ancestry and allele-age patterns to identify selection signals, ancient haplotypes, and population structure, while also revealing uncertainty in some interpretations.

  • π and ARG-based TMRCA both reached minimum values in coding exons and increased with distance from exon boundaries.
  • Genic regions showed reduced coalescence rates but similar relative depths of middle and extreme events, consistent with background selection dominating genes.
  • Four HLA-region blocks ranked among the genome’s top twenty TMRCA blocks, supporting unusually ancient ancestry associated with balancing selection.
  • ARG-based TMRCA estimates for shared human–chimpanzee haplotypes ranged from roughly 150,000 to 250,000 generations and were consistent with pre-divergence origins.
  • Allele ages differed across functional classes, and ARGweaver extracted selection-related information beyond the site-frequency spectrum alone.
  • Population patterns recovered major African–Eurasian and within-region splits, while admixture analyses identified gene flow involving Maasai, African American, Mexican, and Puerto Rican populations.
  • The posterior ARG distribution contained useful population-structure information, but inferred networks could not reconstruct the precise topology of the Latino, European, and African subnetwork.

Discussion

The discussion presents explicit genome-wide ARG inference as a scalable framework for studying genealogies, natural selection, and demographic history, while identifying computational, modeling, and data-related limitations. It emphasizes that posterior ARG samples enable estimation of diverse ARG-derived quantities, but broader demographic and selection applications require further development.

  • Limitations and future directions: Subtree threading improves sampler mixing, but extending the framework to full parametric multipopulation demographic models remains technically challenging and may cause slow convergence and long runtimes.The authors suggest that suitable heuristics and approximations may make joint inference feasible.
  • Scalable ARG inference: Explicit ARG inference reaches complete mammalian genomes for the first time, improving scalability over methods restricted to small sequence sets or requiring severe approximations.The approach is based on an efficient threading operation and is implemented with MCMC.
  • ARG-derived quantities: Posterior sampling of complete ARGs enables estimation of quantities such as TMRCA, allele ages, identity-by-descent regions, recombination events, and coalescence times.This extends inference beyond properties directly described by model parameters.
  • Limitations and future directions: The method requires MCMC rather than direct likelihood or maximum-likelihood estimation, and long runtimes plus unwieldy collections of ARG samples may limit practical usefulness.Additional constraints include possible biases from the SMC approximation, coarse time discretization, phasing errors, and high variance in allele-age estimates at individual sites.
  • Natural selection: Genome-wide genealogical characterization reveals ancient TMRCAs near balancing selection, reduced TMRCAs near selected regions, and reduced allele ages under direct or linked selection.These analyses also suggest potential for distinguishing protein-coding genes from partial selective sweeps using the relative TMRCA halflife statistic.
  • Demographic inference: ARG-based summaries may integrate mutation and linkage-disequilibrium patterns for demographic inference, but the noninformative panmictic prior shows limited value for estimating divergence times, ancestral population sizes, and gene flow.The approach shows promise for reconstructing population phylogenies under this prior.

Methods

The DSMC discretizes coalescence and recombination histories onto predefined time points, while threading samples recombination and re-coalescence events using branch-based distributions. The resulting HMM procedure can be sampled in O(nK2L) and is correct under the stated DSMC/SMC assumptions.

  • Discretization scheme: The DSMC requires all coalescence and recombination events to occur at K + 1 predefined time points.The time points run from s0 = 0 to a user-specified maximum sK.
  • Discretization scheme: Logarithmic spacing places finer time resolution near the ARG leaves, where more events are expected.Increasing δ clusters time points more densely near the leaves; the stated settings are sK = 200,000 generations and δ = 10.
  • Time rounding: Time rounding assigns coalescences to designated points and recombinations to the next most recent point, with special handling near the root.The root heuristic redistributes some probability mass to the root so relevant trees remain reachable by threading.
  • Discretization caveat: Discretization can shrink short branches to zero and distort data generation or inference, so the method samples events first and then chooses active branches heuristically.Simulations indicate that this zero-length-branch solution works fairly well in practice.
  • Recombination sampling: Recombination points are sampled by choosing a time interval in proportion to total branch length and then selecting an active branch.The SMC restricts recombination points to times no older than the root of the existing tree.
  • Re-coalescence sampling: The re-coalescence distribution accounts for whether the re-coalescence time is later than or equal to the recombination time and uses rates from the tree excluding the detached branch.The detached branch is reattached through an SPR operation, producing the new local tree.
  • HMM threading: The threading procedure samples the target conditional distribution correctly for every sequence position and has computational cost O(nK2L).Correctness follows by induction on suffixes of the sequence alignment.

Figure Legends

The figure legends depict ARGs as sequences of local trees connected by recombination events and explain threading as adding one sequence through an existing n−1-sequence ARG. They also describe simulations and genomic analyses comparing inferred genealogical features with truth and selection-associated patterns.

  • ARG representation: An ARG is represented by local trees along sequence positions, with recombination events connecting changes between neighboring trees.Together, local trees and recombinations provide a complete description of the ARG.
  • ARG representation: The SMC approximates the full coalescent-with-recombination by making each local tree conditionally independent of previous trees given its predecessor.The DSMC further places events precisely at predefined time points.
  • Threading: Threading removes one sequence from an ARG, samples its coalescence points across existing local trees, and then samples recombination points between adjacent trees.This two-step process constructs new local trees and reconciles them across breakpoints.
  • Threading: The graphical model integrates recombination variables to obtain an HMM and samples coalescent threadings efficiently through the reduced model.The recombination variables are then sampled conditionally on the coalescence threading.
  • Simulation evaluation: Simulations compare true and inferred log joint probability, recombination counts, and total ARG branch length across 100 data sets.Posterior means and credible intervals summarize sampled ARGs after burn-in.
  • Simulation evaluation: For a simulated 1-Mb segment, posterior mean TMRCA is plotted against true TMRCA with 95% credible intervals.The example uses 20 sequences generated with µ/ρ = 6.
  • Human genomic analysis: Across 17,845 protein-coding genes, nucleotide diversity, TMRCA, and relative TMRCA half-life are summarized for 26 haploid European genomes.Noncoding estimates use 300-bp sliding windows, while coding analyses focus on fourfold-degenerate sites.
  • Human genomic analysis: Sweeps show stronger reductions in nucleotide diversity, whereas genes show stronger reductions in TMRCA and sweeps stronger reductions in relative TMRCA half-life.These signatures are described as consistent with a dominant influence from background selection rather than hitchhiking in protein-coding genes.

Supplementary Figure Legends

The supplementary results describe simulation settings, estimator accuracy, and genomic patterns recovered by ARGweaver, including recombination, ancestry, allele age, and selection signals.

  • Simulation settings: Simulation comparisons use coalescent-with-recombination, SMC, and DSMC models across mutation-to-recombination ratios and time discretizations.DSMC analyses consider K = 40, 20, and 10 time intervals, while several experiments vary µ/ρ from 1 to 6.
  • Simulation accuracy: ARGweaver estimates generally remain accurate across tested settings, although recombination counts show a slight under-estimation, especially at low µ/ρ.The under-estimation is attributed to approximations inherent in the DSMC.
  • Recombination-rate inference: In simulations, inferred recombination rates correlate with true rates and reveal hotspots, but hotspot magnitudes are dampened by the uniform prior.Only recombinations producing changes in tree topology are included in the plotted true rate.
  • Posterior calibration: ARGweaver branch posterior probabilities more accurately reflect branch correctness than branch-inference frequency from Margarita.The comparison concerns local tree branch posterior probabilities and their probability of correctness.
  • Selection signals: Regions predicted to be under balancing selection show larger TMRCAs, while strongly selected annotation classes exhibit sharper allele-age reductions as frequency increases.Linked selection can produce comparable or larger age reductions at low allele frequencies, whereas direct selection tends to dominate at high frequencies.

Supplementary Information: Text S1

The supplementary information lists the authors and their Cornell University affiliation.

  • Matthew D. Rasmussen, Melissa J. Hubisz, Ilan Gronau, and Adam Siepel are listed as authors.
  • The authors are affiliated with Cornell University.
  • The listed affiliation is the Department of Biological Statistics and Computational Biology.

Supplementary Methods

The supplementary methods derive efficient HMM transition calculations for threading branches through local trees, covering recombination and recoalescence configurations. They reduce individual transition evaluation from quadratic to linear time in K and reduce forward-table complexity from O(n^2K^2) to O(nK^2).

  • Transition cases: The HMM transition probabilities are simplified by separating configurations with no old recombinations from those with old recombinations.Each major case is further divided into three subcases based on the thread’s previous and current coalescence points.
  • State representation: The threading state records a branch and time, while recombination states record the recombination branch and time.The new branch being threaded is denoted v; single-sequence threading assumes v is external, with subtree threading handled similarly.
  • No-old-recombination case: With no old recombinations, transitions distinguish recoalescence to different branches, the same branch at different times, and the same branch at the same time.The allowed recombination locations and time ranges depend on the relevant branches and coalescence times.
  • Old-recombination case: With an old recombination, deterministic, recombination-point, and recoalescence-point cases determine the possible next thread states.The deterministic case has one nonzero transition, whereas the other cases enumerate states compatible with recombination and recoalescence geometry.
  • Computational efficiency: O(K) time evaluates individual transition probabilities instead of the naive O(K^2) calculation.The improvement follows from simplifying the relevant transition sums, each containing no more than K + 1 terms.
  • Computational efficiency: The forward table takes O(nK^2) rather than O(n^2K^2) time, reducing the algorithm’s total running time.The methods support Gibbs and generalized Metropolis-Hastings Markov chain Monte Carlo samplers through repeated threading operations.

Supplementary Data Analysis

The simulations compare discretized SMC behavior with SMC and CwR models across recombination rates and population sizes. The models produce highly similar distributions, indicating minimal effects of the discretization scheme at realistic human-population parameters.

  • Recombination-count comparisons: At all recombination rates, DSMC, SMC, and CwR produce very similar distributions of recombination counts.The distributions are essentially indistinguishable at lower recombination rates, with only a slight DSMC excess at higher rates.
  • Recombination-count comparisons: The DSMC is not highly sensitive to the number of time intervals K, although coarse discretization most increases high-rate recombination events.
  • Population-size comparisons: The number of segregating sites increases approximately linearly with effective population size N under all models.The comparison uses a single µ/ρ ratio of 2 and population sizes ranging from 10,000 to 30,000 individuals.
  • Overall interpretation: The discretization scheme has at most a minimal effect on measurable mutation and recombination patterns at realistic human-population parameter values.This result supports the adequacy of the model for inference in the tested setting.
Loading 1306.5110v3…