De novo assembly of viral quasispecies using overlap graphs

Abstract

A viral quasispecies, the ensemble of viral strains populating an infected person, can be highly diverse. For optimal assessment of virulence, pathogenesis, and therapy selection, determining the haplotypes of the individual strains can play a key role. As many viruses are subject to high mutation and recombination rates, high-quality reference genomes are often not available at the time of a new disease outbreak. We present SAVAGE, a computational tool for reconstructing individual haplotypes of intra-host virus strains without the need for a high-quality reference genome. SAVAGE makes use of either FM-index-based data structures or ad hoc consensus reference sequence for constructing overlap graphs from patient sample data. In this overlap graph, nodes represent reads and/or contigs, while edges reflect that two reads/contigs, based on sound statistical considerations, represent identical haplotypic sequence. Following an iterative scheme, a new overlap assembly algorithm that is based on the enumeration of statistically well-calibrated groups of reads/contigs then efficiently reconstructs the individual haplotypes from this overlap graph. In benchmark experiments on simulated and on real deep-coverage data, SAVAGE drastically outperforms generic de novo assemblers as well as the only specialized de novo viral quasispecies assembler available so far. When run on ad hoc consensus reference sequence, SAVAGE performs very favorably in comparison with state-of-the-art reference genome-guided tools. We also apply SAVAGE on two deep-coverage samples of patients infected by the Zika and the hepatitis C virus, respectively, which sheds light on the genetic structures of the respective viral quasispecies.

Figure 1.

An overview of the workflow and algorithms of SAVAGE. (A) The three stages of SAVAGE. Each assembles sequences into longer sequences. For clarity, we assign different names to the sequences output by each stage: contigs, maximally extended contigs, and master contigs, respectively. (B) Principle of overlap graph construction and distinction among the reads between errors and shared mutations. (C) Each stage has two steps: first, the overlap graph construction; second, assembly. This panel summarizes the differences in each step between the three stages. During overlap graph-based assembly, steps 4 to 6 are repeated iteratively until there are no edges left in the overlap graph.

Figure 2.

Target genome fraction recovered per strain for all 20,000× benchmarks, stratified by strain frequency.

Figure 3.

Performance of SAVAGE-de-novo and SAVAGE-b-ref, depending on pairwise distance and mixture ratio. (A) Target genome fraction recovered (%) considering all maximally extended contigs ≥500 bp. (B) Overall mismatch rate (%) considering all maximally extended contigs ≥500 bp. (C) Relative error of estimated frequency for the minor strain (%). Frequency estimates were computed using Kallisto, and only assemblies containing exactly two maximally extended contigs longer than 4000 bp were evaluated.

Figure 4.

Edge criteria. For an overlap to become an edge in the overlap graph, it must satisfy three criteria. First, the overlap length l must be at least the minimal overlap length L. Second, the overlap quality score QS(R₁, R₂) must be at least the minimal score δ. For overlaps involving paired-end reads, we require both l₁ ≥ L and l₂ ≥ L, and, analogously, $QS(R_{1_a},R_{2_a})\ge \delta$ and $QS(R_{1_b},R_{2_b})\ge \delta$. Finally, we only accept overlaps where the sequence orientations of a paired-end read agree: either both sequences in forward orientation, or both sequences in reverse orientation.

Figure 5.

Algorithmic details. (A) Read orientations: Given an edge u→v with orientations (−,+). Then, if u is labeled +, the induced label for v is −, while if u is labeled −, the induced label for v is +. This procedure leads to a vertex labeling in O(V) time. (B) Transitive edges: An edge u→w is called nontransitive, shown in black, if there is no vertex v such that there are edges u→v, v→w. It is called single transitive, shown in green, if for all vertices v such that there are edges u→v, v→w, one of the edges is nontransitive. It is called double transitive, shown in red, if there is a vertex v with edges u→v, v→w which are both transitive. (C) Read clustering by cliques (top) or by pairs (bottom). (D) Error correction: When a consensus sequence is constructed from a cluster of reads, the extremities are removed.