Source-linked AI summary
Viral population estimation using pyrosequencing
Nicholas Eriksson, Lior Pachter, Yumi Mitsuya, Soo-Yon Rhee, Chunlin Wang, Baback Gharizadeh, Mostafa Ronaghi, Robert W. Shafer, Niko Beerenwinkel
TL;DR
Within-host viral diversity is difficult to quantify from short, error-prone reads, yet matters for disease progression, vaccines, and therapy. The paper develops a reference-guided pipeline combining error correction, minimal graph-based haplotype reconstruction, and EM frequency estimation. Simulations and clonal-sequencing comparisons show effective quantitative reconstruction, while performance is constrained by population diversity and reference alignment.
Problem
Within-host viral diversity matters for disease progression, vaccine design, and antiretroviral therapy, but sequencing errors and short reads complicate population-structure inference.
Method
The method aligns reads, corrects errors locally, constructs a minimal read-consistent haplotype set, and estimates haplotype frequencies by maximum likelihood using EM.
Results
Pyrosequencing effectively reconstructs virus-population diversity, supported by extensive simulations and comparisons with direct clonal sequencing.
Takeaways & Limitations
Pyrosequencing can quantitatively assess the diversity of RNA-virus populations such as HIV.
Takeaways & Limitations
The method works best for suitably diverse populations and becomes less accurate with highly diverse populations or true indels because reference alignment degrades.
Abstract
from arXiv · showhide
The diversity of virus populations within single infected hosts presents a major difficulty for the natural immune response as well as for vaccine design and antiviral drug therapy. Recently developed pyrophosphate based sequencing technologies (pyrosequencing) can be used for quantifying this diversity by ultra-deep sequencing of virus samples. We present computational methods for the analysis of such sequence data and apply these techniques to pyrosequencing data obtained from HIV populations within patients harboring drug resistant virus strains. Our main result is the estimation of the population structure of the sample from the pyrosequencing reads. This inference is based on a statistical approach to error correction, followed by a combinatorial algorithm for constructing a minimal set of haplotypes that explain the data. Using this set of explaining haplotypes, we apply a statistical model to infer the frequencies of the haplotypes in the population via an EM algorithm. We demonstrate that pyrosequencing reads allow for effective population reconstruction by extensive simulations and by comparison to 165 sequences obtained directly from clonal sequencing of four independent, diverse HIV populations. Thus, pyrosequencing can be used for cost-effective estimation of the structure of virus populations, promising new insights into viral evolutionary dynamics and disease control strategies.
Synopsis
Viral population estimation proceeds from aligned pyrosequencing reads through error correction and haplotype assembly to maximum-likelihood frequency estimates.
- Reads are aligned to a reference strain, corrected for errors, and assembled into haplotype candidates before frequency estimation.
- Maximum-likelihood estimates of reconstructed haplotype frequencies constitute the inferred virus population.
1 Introduction
The paper addresses the challenge of reconstructing within-host viral population structure from short, error-prone pyrosequencing reads. It develops a reference-guided computational procedure and applies it to HIV populations.
- Pyrosequencing is faster, less laborious, and cheaper than Sanger sequencing, but produces shorter, more error-prone reads.
- Viral genetic diversity affects disease progression, vaccine design, and antiretroviral drug therapy, motivating quantitative population-structure inference.
- Each read originates from one unknown haplotype, so the task is to reconstruct consistent haplotypes and infer their relative frequencies.
- The procedure aligns reads, corrects errors locally, assembles consistent haplotypes with graph paths, and estimates frequencies by maximum likelihood.
- The method is validated through simulations and comparisons with direct clonal sequencing of four independent HIV populations.
- The methods apply when a reference genome is available, read coverage is sufficient, and haplotypes are sufficiently separated genetically.
2 Results
The paper combines local error correction with statistical testing and clustering to prepare pyrosequencing reads for haplotype reconstruction and frequency estimation.
- The procedure begins with reference alignment and error correction, followed by haplotype reconstruction and frequency estimation.
- Because pyrosequencing errors include indels, substitutions, and ambiguous bases, error correction is necessary before inferring the virus population.
- Error correction: Fixed-width alignment windows are clustered using tests for over-represented mutations and co-occurring mutation pairs.
- Error correction: Reads are clustered into estimated haplotype groups and corrected to quality-score-weighted consensus sequences.
Algorithm 1. (Local error correction)
Local correction and graph-based reconstruction seek a minimal, read-consistent haplotype explanation, while sequencing coverage and error behavior constrain detectable population diversity.
- Algorithm 1. (Local error correction): The local correction algorithm tests candidate mutations and mutation pairs, clusters reads into k groups, corrects them to cluster centers, and outputs corrected reads.
- Algorithm 1. (Local error correction): The correction procedure can leave errors or introduce misc corrections through false positives, false negatives, or misclustering.
- Algorithm 1. (Local error correction): Misclustering becomes more likely as window size grows and reads overlap the window only partially.
- Algorithm 1. (Local error correction): At 2.5 errors per kb, variants below 1% frequency would be erased in a dataset of 10,000 reads.
- Haplotype reconstruction: A completely consistent haplotype is constructible from overlapping reads that agree on their overlaps, and the read graph represents these possibilities.
- Haplotype reconstruction: A source-to-sink path in the read graph corresponds to a completely consistent haplotype, enabling efficient path enumeration.
- Coverage and resolution: With 99% probability, 2,302 reads cover haplotypes at least 5% frequent and 11,508 reads cover haplotypes at least 1% frequent under the stated model.
- Haplotype reconstruction: The reconstruction target is a small haplotype set explaining the reads, found through graph coverage by paths.
Theorem 4. (Minimal cover of the read graph)
The read graph has a well-defined minimal-cover size, even though the corresponding path cover and chain decomposition may not be unique. This size is the smallest number of haplotypes that can explain the reads.
- Every minimal read-graph cover has cardinality equal to the size of the largest vertex set with no paths between its elements.
- A minimal cover can be computed through a maximum-matching problem in an associated bipartite graph.
- The matching problem requires at worst cubic time in the number of irredundant reads.
- The minimal path cover and minimal chain decomposition are generally non-unique, but their cardinality is a well-defined global invariant of the reads.
- The minimal cover size is the smallest number of haplotypes capable of explaining the observed data.
Algorithm 5. (Construction of a minimal set of explaining haplotypes)
Algorithm 5 constructs a minimal set of explaining haplotypes by converting error-corrected reads into graph paths. It can also generate additional candidate paths, while relaxing consistency requirements when coverage is too low.
- Algorithm 5 builds the read graph, computes a minimal chain decomposition, extends chains to source-to-sink paths, and outputs the corresponding haplotypes.
- The algorithm samples candidate haplotypes from minimal explaining sets rather than enumerating all graph paths.The HIV datasets produced up to 10^9 paths, while simulations often produced more than 10^12.
- The generative model represents five genomes with two haplotypes distributed as p = (3/5, 2/5), then samples reads uniformly from those haplotypes.
- When coverage is too low for a completely consistent explaining set, the method adds edges between non-overlapping reads before recomputing a minimal set.
2.3 Haplotype frequency estimation
Haplotype frequency estimation uses a generative model in which reads arise from unknown haplotype frequencies, then fits those frequencies by maximum likelihood with EM. Candidate haplotypes are restricted to a feasible explaining set.
- The candidate set H contains explaining haplotypes rather than all possible haplotypes, making frequency estimation feasible.
- The model draws a haplotype according to p and then samples uniformly from reads consistent with that haplotype.
- A read has conditional probability 1/K when it is consistent with a haplotype and probability 0 otherwise.
- Because haplotypes have equal length, K is independent of both the read and haplotype.
- The method estimates p by maximizing the log-likelihood of the observed read counts.
- The EM algorithm runs in O(|R||H|) per iteration and typically converges within minutes for 5000 reads and 200 candidate haplotypes.
2.4 Simulation results
Simulations evaluated error correction, haplotype reconstruction, frequency estimation, and their combined performance across varying diversity, population size, read depth, and error rates. The procedure generally improved with greater diversity but degraded with sequencing errors and larger populations.
- Error correction: 30-fold error-rate reduction left approximately 0.1 errors per kb after correction, although reduction decreased as populations became larger and more diverse.The post-alignment error rate was about 1–3 errors per kb.
- Haplotype reconstruction: At 3% diversity, one of five haplotypes was reconstructed well; at 7%, chain decomposition exactly reconstructed all five.At 5% diversity, reconstruction was almost correct except for a few small crossovers.
- Frequency estimation: EM frequency estimation significantly outperformed the read-count heuristic at every tested read-set size, with KL divergence decreasing roughly exponentially as reads increased.The improvement in prediction accuracy increased with the number of reads.
- Combined procedure: Combined reconstruction performance increased with diversity and dropped slightly as the number of equal-frequency haplotypes increased, across populations of 5–100 haplotypes.These tests used 10,000 error-free reads and diversities between 3% and 8%.
- Combined procedure: At 0.1 errors per kb, approximately 3,500 reads recovered about 55% of a population with 5% diversity.Performance decreased as errors were introduced, but the combined procedure still recovered much of the original population.
- Combined procedure: At a 0.2 error rate, performance showed a small loss as read count increased because additional reads created more graph paths and erroneous haplotypes received positive probabilities.The minimal path-cover size increased approximately linearly with the number of reads.
2.5 Analysis of HIV samples
The method was evaluated on drug-resistant HIV populations from four patients using paired pyrosequencing and clonal Sanger sequencing. Pyrosequencing closely mirrored clonal sequencing, while also producing additional haplotypes whose status could not be resolved from these data.
- Validation against clonal sequencing: Across four HIV samples, at least 51.8% of the inferred populations matched a clonal haplotype within one amino acid difference.The authors used ϕ1 for this more sensitive comparison instead of ϕ10.
- V11909 analysis: In the V11909 sample, 61 of 65 clonal haplotypes were matched in the inferred population, and inferred frequencies reasonably matched clonal mutation-pattern frequencies.The unmatched count was four haplotypes, or 6.1%.
- Resolution: The sequencing depth was sufficient, under the Lander–Waterman model, to reconstruct with 99% probability haplotypes occurring at frequencies of at least 2.2%.The comparison used pyrosequencing reads from the HIV samples.
- Validation against clonal sequencing: Sanger sequencing yielded 65 clonal sequences, 37 of which were mixtures of two or more clones.This complicates direct interpretation of the clonal reference set.
3 Discussion
The paper frames pyrosequencing population reconstruction as a sequence of computational problems and shows that the resulting methods can estimate viral diversity at useful scale. It also identifies resolution, diversity, read length, and error structure as important boundaries for interpretation and future development.
- Method: The procedure separates error correction, haplotype assembly, and frequency estimation into three computationally manageable tasks.The method combines local clustering, combinatorial haplotype reconstruction, and an EM algorithm for maximum-likelihood frequency estimation.
- Method: A minimal path cover of the read graph provides a computable measure related to population diversity.Consistency requires haplotypes to be constructible from overlapping reads, while parsimony minimizes the explaining set.
- Validation: Extensive simulations and comparisons with direct clonal sequencing show effective quantitative assessment of RNA-virus population diversity.The validation included four independent HIV populations and 165 clonal sequences.
- Resolution and scale: Approximately 10,000 reads can detect variants present at 1% frequency, while 200,000 reads could support sequencing about twenty populations at comparable resolution.The latter process is described as less labor intensive than limiting-dilution clonal sequencing for one population.
- Limitations: The method works best for suitably diverse populations, whereas repeats, true indels, and greater diversity can reduce linkage or alignment accuracy.Low diversity weakens read linkage; sufficiently long repeats can destroy it, and more diverse samples or indels challenge single-reference alignment.
- Limitations: Error-correction resolution is limited by the overall error rate, and haplotype reconstruction is the least immediately practical step for datasets containing 200,000 reads.The reconstruction procedure is at worst cubic in the number of irredundant reads, although redundant-read elimination improves feasibility.
4 Methods
The methods separate error correction, haplotype reconstruction, and frequency estimation into manageable computational tasks. Statistical tests identify distinct haplotypes, graph algorithms construct explaining haplotypes, and an EM algorithm estimates their frequencies.
- Statistical tests for error correction: Column-wise binomial tests identify mutation counts that exceed expectations under a one-haplotype sequencing-error model.The test uses the overlap depth d, error rate ϵ, and significance threshold α; with d≈1000, ϵ=0.0025, and α=0.001, nine mutations are significant.
- Statistical tests for error correction: Fisher’s exact tests evaluate whether mutations in two alignment columns co-occur more often than expected under the one-haplotype null model.The co-occurrence count C follows a hypergeometric distribution, with mutation counts n_u and n_v defining the marginal observations.
- Statistical tests for error correction: The correction rule subtracts two haplotypes when significant mutations occur together on exactly the same reads, preventing over-counting.This handles cases such as two observed mutations that may represent only the haplotypes AA and CC rather than four inferred haplotypes.
- Haplotype reconstruction: The read graph is decomposed into a minimum number of chains using a bipartite maximum-matching construction derived from Dilworth’s theorem.The associated bipartite graph connects reads when one is reachable from another in the acyclic read graph.
- Haplotype frequency estimation: An EM algorithm estimates maximum-likelihood haplotype frequencies by iteratively estimating read-to-haplotype assignments and optimizing the hidden-model likelihood.Simulations used 200 candidate haplotypes and 10 random EM starting points, setting frequencies below 10^-6 to zero.
- Evaluation: The evaluation measures haplotype recovery, Hamming distance, KL divergence between frequency distributions, and the inferred population fraction near an original haplotype.The latter statistic is denoted ϕ_n, where n is the permitted number of amino-acid differences.
Supporting Information
Supporting analyses examine reconstruction accuracy, error correction, and lower bounds on population size across simulated viral populations. They show that multiple chain decompositions can recover haplotypes, error correction substantially reduces errors, and reconstruction resolution depends on reads and errors.
- Haplotype reconstruction: Multiple chain decompositions can recover different haplotypes from the same read graph, allowing all haplotypes to be reconstructed.In the example, one decomposition accurately recovers one haplotype, while another accurately recovers two different haplotypes.
- Haplotype reconstruction: The reconstruction-efficiency analysis measures the average amino-acid Hamming distance between each original haplotype and its closest reconstructed match.It evaluates populations of 10, 20, and 50 haplotypes using up to 1000 candidate haplotypes generated from 10,000 error-free reads.
- Error correction: Error correction reduced simulated pyrosequencing error rates by approximately 30-fold for populations with 4% diversity.The simulations used populations of up to 50 equally frequent haplotypes and initial post-alignment error rates of 3–6 errors per kb.
- Population-size lower bound: The minimal explaining-haplotype cover provides an accurate lower bound for error-free reads, whereas introduced errors make the bound increase approximately linearly with read count.Simulations varied population diversity from 3–7% and summarized variability using interquartile-range error bars over 100 trials.