WgLink: reconstructing whole-genome viral haplotypes using L0 + L1-regularization

Many tools can reconstruct viral sequences based on next generation sequencing reads. Although existing tools effectively recover local regions, their accuracy suffers when reconstructing the whole viral genomes (strains). Moreover, they consume significant memory when the sequencing coverage is high or when the genome size is large. We present WgLink to meet this challenge. WgLink takes local reconstructions produced by other tools as input and patches the resulting segments together into coherent whole-genome strains. We accomplish this using an L0 + L1-regularized regression synthesizing variant allele frequency data with physical linkage between multiple variants spanning multiple regions simultaneously. WgLink achieves higher accuracy than existing tools both on simulated and real data sets while using significantly less memory (RAM) and fewer CPU hours. Source code and binaries are freely available at https://github.com/theLongLab/wglink.


INTRODUCTION 25 26
In the viral sequencing projects based on next-generation sequencing (NGS) instruments, 27 researchers typically sample a collection of viruses in tissue rather than targeting individual 28 viruses. As a result, sequencing reads represent a mixture of viral particles reflecting the within-29 host viral population, the constituents of which are called viral quasispecies. Reconstructing 30 haplotypes from this mixture is an important but difficult computational problem. 31 Haplotype reconstruction from pooled NGS data is a prerequisite for many downstream 32 biological analyses focusing on individual strains. This has led bioinformaticians to develop tools 33 facilitating such reconstruction. These tools employ an range of computational and statistical 34 techniques, e.g. graph algorithms (1, 2) regularization (3), Tenser factorization (4), Dirichlet 35 processes (5), and Bayesian methods (6). 36 37 Many tools can effectively reconstruct haplotypes locally at the gene level (4,6). However, the 38 accuracy of global reconstruction, i.e. at the level of the whole viral genome (strain) remains low. 39 Moreover, the techniques for such global reconstruction consume significant memory resources. 40 There is much room for improvement. 41

42
In this work, we present WgLink, a novel tool using regularized regression to stitch together 43 regional haplotypes into global ones. As input, WgLink takes regional haplotypes produced by 44 existing tools, the reference genome of the focal species, and the sequencing reads of a sample. 45 It outputs whole viral genome haplotypes present in the sample together with their within-host 46 frequencies. WgLink's novelty lies both in its design, using a regression incorporating both allele 47 frequencies and physical linkage disequilibrium (PLD) observed in sequencing reads, and its 48 implementation, solving this regression problem using a state-of-the-art tool for extra-fast 0 -49 regularization (7) (Methods). As a result, WgLink generates more accurate viral haplotypes 50 using less memory and fewer CPU hours. 51 52 METHOD 53 54 55 The WgLink algorithm. Based on regional haplotypes generated by some other tool, WgLink 56 first generates all possible whole-genome haplotypes using a breadth-first search 57 (Supplementary Notes). To collect allele frequency and LD information in a consistent file 58 format, WgLink uses BWA (8) to map reads back to the reference genome. Subsequently, 59 WgLink performs a regularized regression to estimate the haplotype frequencies of all 60 candidates. 61 62 As indicated above, the design of this regularized regression contains two key innovations. 63 First is the synthesis of allele frequencies at individual sites and the LD between pairs of sites in 64 the context of a single regression problem. We explain why these innovations lead to the claimed 65 improvements. Many existing tools aiming to reconstruct haplotypes targeting different species 66 use allele frequencies only (3,9,10). More specifically: Letting Y denote the vector of within-67 host frequencies of all genetic variants and = ( 1 , , …, ) denote the potential haplotypes, 68 represented by vectors of 0s and 1s (representing reference and alternate alleles), then the 69 linear regression of ~∑ gives an estimate of the haplotype frequencies 1 , 2 , … , . 70 ( Supplementary Fig. S1). Note that in a regularized form, most of the are 0. The sample size 71 for solving this regression is the number of genetic variants. This model is useful only if there 72 are sufficient polymorphic sites to distinguish the various haplotypes. However, this linear 73 regression approach using allele frequency alone is insufficient for distinguishing haplotypes that 74 differ subtly, e.g. by a few recombination events. For instance, the actual haplotypes and 75 frequencies in Fig. 1, (x1:1/6, x2:2/6, x3:3/6), will have the same allele frequencies as the 76 configuration (x4:1/6, x5:2/6, x6:3/6) and, thus, won't be distinguishable by regression on allele 77 frequency only. As remedy for this insensitivity, WgLink incorporates a quantity derived from 78 pairs of alleles on the same sequencing read (or paired-end reads) that we call Physical Linkage 79 Disequilibrium, or PLD (in contrast to statistical LD in population genetics). More precisely, for 80 any two sites, one sample will be added as follows: the proportion of reads carrying allele 1 on 81 both sites (out of all reads covering both sites) will be added as a Y-value. The corresponding xi 82 value will be set to 1 if the corresponding i-th potential haplotype contains allele 1 on both sites 83 and to 0, otherwise. As such, samples capturing the PLD, i.e., the blue diamonds in

95
The second innovation is adding of a mixture of 0 and 1 regularization to the regression 96 problem. Existing tools use 1 regularization (i.e., LASSO) only (9), which does not discriminate 97 feasible solutions for which the haplotype frequencies sum (correctly) to 1. For example, an 1 -98 regularized regression does not favor a correct solution of 2 haplotypes, each with frequency 99 0.5, over a competing, incorrect solution of 10 haplotypes, all with frequency being 0.1, because 100 of the penalty term being 1.0 in both cases. RegressHaplo attempted to address this using a 101 revised 1 -regularization term. To fundamentally solve this problem, we add an 0 -regularization 102 term, imposing a much stiffer penalty on high-complexity models (Supplementary Notes). 0 -103 regularized regression is a non-convex optimization problem whose solution was, until recently, 104 extremely slow. This changed with the release of L0Learn package (7), an exciting breakthrough 105 in the field, that can solve 0 -regularized regression problems efficiently. Although the L0Learn 106 paper does not present a precise estimate of time-complexity, our benchmarking shows that our 107 L0Learn-based implementation of WgLink is significantly faster than other tools 108 Table S1). 109

Detailed algorithm for the breadth-first search
The local haplotypes produced by the other reconstruction tool are linked into potential global haplotypes using a breadth-first search (BFS) algorithm. As haplotypes are inferred for each tiling window, haplotypes from adjacent windows will overlap over about half of their sequence length. If two overlapping haplotypes match on this overlap, they are stitched together to form a longer haplotype. Since the overlapping regions may not be identical due to errors in the previous steps, two haplotypes may be stitched together as long as the number of mismatches in their overlapping region is less than a user-specified cutoff parameter.
Next, global haplotypes are determined by traversing a tree structure whose nodes are the locally stitched haplotypes. The i-th level of the tree represents the i-th genomic region. We draw an edge between two nodes in neighboring levels if and only if their corresponding haplotypes can be stitched together.
Paths that traverse the tree enumerated using BFS, represent candidate global haplotypes, as stitched together local ones. These are vetted further in WgLink's regularized regression phase.

Real Data, Simulations, and HPC specs
Real data. For benchmarking on real data, we used a public dataset containing a five-strain(1). Specifically, this dataset describes an in vitro generated quasispecies population consisting of Simulations. We used the above five-strain as templates for simulated data. In real world, haplotype abundance is unknown and random. To reflect this, we randomly generated goldstandard frequencies for each haplotype. Based on the gold-standard haplotypes (i.e., the five strains in real data) and their frequencies, DWGSIM (2) (https://github.com/nh13/DWGSIM) was used to generate paired-end reads for a mixture of the haplotypes. We varied three factors in our construction: sequencing error (low = 0.1%, high = 0.4%), read length (low = 150bp, high = 250bp) and average sequencing coverage per base pair (low = 2500X, high = 10000X), giving a total eight combinations (see Supplementary Figures 2 -9). For each combination of simulation parameters, we performed fifty (50) simulation rounds and plotted the averages. To avoid unusually low values, we discarded any haplotype with abundance less than 1%, replacing it with a new simulation.
Data processing. For both real and simulated data, we used BWA (version 0.7.15) (3) and SAMTools (4) (version 1.4) mapping and sorting reads, respectively. We used HXB2, the standard reference genome for HIV studies, as the reference genome.

Parameters and commands for different tools
We used the default parameters of TenSQR, aBayesQR is WgLink. According to the TenSQR(5) and aBayesQR manual(6), specified min_mapping_qual, min_read_length and MEC improvement threshold were used for the real data analysis.

Comparison criteria Estimated number of haplotypes
The estimated numbers of haplotypes are the numbers of predicted haplotypes by different tools (including WgLink and alternative tools). The error in an estimate is simply the difference between the prediction and the actual number.

Reconstruction Rate (two versions)
Reconstruction rate is, by definition, the number of mismatches between the gold-standard halophytes and the reconstructed ones. Since these may differ in number, we need an explicit mapping between gold-standard and reconstructed haplotypes to facilitate the comparison. There are two natural ways constructing such: mapping the gold-standard haplotypes to reconstructed haplotypes or mapping the reconstructed haplotypes to gold-standard haplotypes, corresponding to Reconstruction Rates 1 and 2 (see below) respectively.
Reconstruction Rate 1 is defined as where the , = 1, . . . , , are the gold-standard haplotypes, is the predicted haplotype closest to in terms of minimum Hamming distance, and ℎ is the number of mismatch sites between .

Reconstruction Rate 2
Reconstruction Rate 2 is defined as where the , = 1, . . . , , are the predicted haplotypes, is the gold-standard haplotype closest to in terms of minimum Hamming distance, and ℎ is the number of mismatch sites between .

Jensen-Shannon divergence (JSD)
The Jensen-Shannon divergence (JSD) evaluates the divergence between the frequencies of the simulated and reconstructed haplotype distributions and represents our tool's ability to predict haplotype abundance. The probability distributions A and P represent the frequency distributions gold-standard and predicted haplotypes, respectively. JSD is defined as The quantity D(A||P) is the Kullback-Leibler divergence from A to P, defined by D( ∥ P) = � ( ) log 2 ( ) ( ) .
JSD is symmetric in A and P and takes values in the interval [0,1], the value 0 representing a perfect prediction.

Supplementary Figures
Supplementary Fig. 1. A naïve regression based on allele frequencies only. Colored bars on the left represents reads belongs to three actual haplotypes (x1, x2, and x3), corresponding to the vectors on the right. x4, x5, and x6 are potential haplotypes with actual frequencies 0. (As stated in the main text, without incorporating physical linkage disequilibrium, it can't distinguish many incorrect configurations.)