Full-Length Envelope Analyzer (FLEA): A tool for longitudinal analysis of viral amplicons

Next generation sequencing of viral populations has advanced our understanding of viral population dynamics, the development of drug resistance, and escape from host immune responses. Many applications require complete gene sequences, which can be impossible to reconstruct from short reads. HIV env, the protein of interest for HIV vaccine studies, is exceptionally challenging for long-read sequencing and analysis due to its length, high substitution rate, and extensive indel variation. While long-read sequencing is attractive in this setting, the analysis of such data is not well handled by existing methods. To address this, we introduce FLEA (Full-Length Envelope Analyzer), which performs end-to-end analysis and visualization of long-read sequencing data. FLEA consists of both a pipeline (optionally run on a high-performance cluster), and a client-side web application that provides interactive results. The pipeline transforms FASTQ reads into high-quality consensus sequences (HQCSs) and uses them to build a codon-aware multiple sequence alignment. The resulting alignment is then used to infer phylogenies, selection pressure, and evolutionary dynamics. The web application provides publication-quality plots and interactive visualizations, including an annotated viral alignment browser, time series plots of evolutionary dynamics, visualizations of gene-wide selective pressures (such as dN/dS) across time and across protein structure, and a phylogenetic tree browser. We demonstrate how FLEA may be used to process Pacific Biosciences HIV env data and describe recent examples of its use. Simulations show how FLEA dramatically reduces the error rate of this sequencing platform, providing an accurate portrait of complex and variable HIV env populations. A public instance of FLEA is hosted at http://flea.datamonkey.org. The Python source code for the FLEA pipeline can be found at https://github.com/veg/flea-pipeline. The client-side application is available at https://github.com/veg/flea-web-app. A live demo of the P018 results can be found at http://flea.murrell.group/view/P018.

Full-length sequences can resolve features that are difficult to assemble from short sequences [8,31]. For instance, Pacific Biosciences SMRT sequences were able to resolve 1.5 kb msg isoforms from Pneumocystis jirovecii, but reads from a 454 instrument could not be assembled correctly [31]. For tracking evolutionary patterns in viral populations, accurately resolving these features provides a more accurate history of the population, which becomes especially relevant when epistatic interactions and linkage between mutations effect phenotypic changes in the pathogen [32][33][34]. For example, studies of HIV env frequently use functional assays to measure the potency with which a given antibody or donor serum neutralizes a specific env strain [35], which requires knowing the full env sequence.
We have developed a pipeline for handling long read HIV env sequencing data from within-host viral populations: the Full-Length Envelope Analyzer (FLEA). FLEA addresses the specific challenges posed by large volumes of such data, e.g., using the sequencing protocols we previously described in Laird Smith et al [36], which also contains an overview of a prototype of FLEA. Here we describe the full pipeline and experimentally demonstrate its ability to resolve populations of closely related variants. FLEA uses state-of-the-art tools and methods at every step and can be accessed through a web browser or on a high-performance cluster. FLEA is readily extensible to other genes and systems. FLEA has recently been used by the authors in two high-profile studies. In [37], we describe how FLEA was used to process PacBio HIV env data from a clinical trial of monoclonal antibody 10-1074. For sequences sampled before and after therapy, FLEA reveals that prior to antibody therapy low-frequency env variants were present with mutations that typically confer resistance to 10-1074. Additionally, when resistance emerges, it emerges multiple times, exploiting many different resistance pathways. FLEA was also used to characterize the longitudinal env population that drove development of a broadly neutralizing antibodies against the apex of the env trimer, sampled from donor PC64 from the Protocol C primary infection cohort [38].
There exist dozens of standalone pipelines developed for analyzing HIV and related sequence data, including longitudinal samples [4,9,12,39]. However, it was necessary to develop a new tool due to HIV env's extensive natural indel variation and the high rate of indels in long PacBio reads, which are especially problematic when any spurious indel in the 2.6kb env amplicon corrupts the reading frame, rendering the sequence uninterpretable. Previous analysis [36] determined that, for PacBio amplicon sequencing of an Env clone (for the set of sequencing and filtering conditions employed therein), 4 out of 5 errors are indels, and these occur more commonly in long homopolymer runs, with the per-base error rate ranging from 1 in 300 for a homopolymer of length 2, but up to nearly 1 in 50 for a homopolymer of length 6. On average, the per-base error rate was around 1 in 200, yielding an average of roughly 13 errors per 2.6kb sequence.
With HIV env, the common strategy of mapping reads to a reference fails because the diversity in variable regions of env, predominantly driven by extensive long insertions and deletions, means that these regions in sampled reads lack homology to those in any heterologous reference sequence, causing alignment-to-reference strategies to fail. Instead, FLEA relies on a finegrained cluster-and-consensus strategy to remove spurious indels from reads. The task is related to Liang et al. (2016), but, rather than distinguishing a small number of variants at 81-91% identity, we must distinguish potentially hundreds of variants that differ by only a handful of bases.
The main contribution of the FLEA pipeline, therefore, is the reconstruction of a population, including accurate inference of relative frequencies of closely-related minority variants, from SMRT sequences alone. In addition, it performs many useful analyses on this population, such as alignment, phylogenetic reconstruction, and selection inference, and provides interactive visualizations for the results. The full pipeline is available as an online resource, or for local installation.

Pipeline
The input to FLEA is a set of FASTQ files from the PacBio RS-II or Sequel. Each set corresponds to one time point, containing circular consensus sequence (CCS) reads, which can be obtained using the "Reads of Insert" protocol on PacBio's SMRTportal or SMRTanalysis tools. Upon completion, the FLEA pipeline produces results as JSON (Javascript Object Notation) files, a standard format for machine (and human-) readable structured data. The logic of FLEA is implemented in Nextflow [40], a workflow framework for deploying parallel pipelines to clusters and clouds.
FLEA consists of multiple sub-pipelines, as shown in Fig 1. Details of the quality and consensus pipelines are depicted in Fig 2. Together, these two pipelines take error-prone CCS reads and convert them into unique high-quality consensus sequences. The alignment pipeline generates a multiple sequence alignment, which is used by multiple methods in the analysis pipeline.

Quality assurance sub-pipeline
The first steps remove low quality reads and filter out common sequencing artifacts. Parameters given in these steps were chosen for full-length HIV envelope sequences from the RS-II or  Sequel platforms. Other reads with different properties (error rates, error models, lengths, homopolymer distributions, etc.) likely require different parameters. All steps are run independently per time point.
1. Filter by error rate. The input FASTQ files contain Phred scores for each base, encoding the probabilities of incorrect base calls. USEARCH [41] is used to remove reads with an expected error rate greater than 1%, computed as the mean of the per-base error probabilities.
2. Trim heads/tails. A fraction of reads from the Laird Smith et al. sequencing protocol contain poly-A or poly-T heads or tails (cause unknown), which can be hundreds of bases long and sometimes contain a small number of other bases. These heads and tails are trimmed with a hidden Markov model (Fig 3) implemented in Pomegranate [42]. The emission probabilities of the model were fixed, and the transitions trained using Baum-Welch. The Viterbi path for each sequence is computed, and bases emitted by head and tail nodes are removed.
3. Filter long runs. Reads with homonucleotide runs longer than 16 bases are discarded. This length was chosen to be twice the length of the longest such run in the LANL HIV database [43].

Filter contaminants and trim reads.
Sample contamination can introduce non-native sequences that interfere with subsequent analyses, so these contaminants must identified and discarded. USEARCH is used to compare reads to a contaminant database and a reference database using usearch_global. Alignments returned from querying the database are then used to trim reads to the gene boundaries. Trimming terminal insertions is vital for the accuracy of downstream tasks, such as length filtering and clustering. The contaminant database contains HXB2 and NL4-3 env, each ubiquitous in labs working with env sequences and a common source of sample contamination. Reads that match with � 98% identity are discarded. Since a 1% error rate cutoff was earlier used, this parameter conservatively ensures that these contaminants are almost certainly identified. The reference database contains thirty-eight sequences representing the major HIV Group M subtypes from the LANL sequence database [43]. Reads with � 70% identity to every sequence in the reference database are discarded. This cutoff is chosen to retain reads remotely similar to HIV Group M while excluding contaminants such as human or bacterial genome reads. If a sample is from SIV, or from a non group-M HIV+ donor, then more appropriate reference sequences should be added to the database. 5. Filter by length. By default, sequences shorter than 90% or longer than 110% of the length of the reference sequence are discarded. However, sequences with large deletions are frequently observed in HIV. These likely represent replication incompetent envelopes, and their reduced length can cause them to be dramatically oversampled due to PCR length bias. Users who want to include these species in their analyses should modify these parameters.
Reads that pass this quality assurance phase have low expected error rates and no homonucleotide runs, are within 70% identity of at least one reference sequence, are (after trimming) no more than 10% different in length than a reference sequence, and do not match the contaminant database. We refer to these sequences as quality-controlled sequences (QCS).
Consensus sub-pipeline for variant identification. Even for highly diverse populations, unique reads in a sequencing run outnumber the true unique variants, predominantly due to sequencing errors. The problem is far more significant in long reads than in short reads, precluding the use of amplicon denoising strategies used to reduce error rates in short read sequencing [44]. To accommodate this effect, the next phase of the FLEA pipeline clusters and combines QCS reads, attempting to infer the true variants in each time point. It also attempts to detect and correct frameshift errors.
All of the following tasks are run separately for each time point, yielding sets of unique inframe consensus sequences. We refer to these sequences as high-quality consensus sequences (HQCS).
1. Cluster. USEARCH is used with the cluster_fast command to generate clusters with 99% nucleotide identity. This parameter approximates the 1% error cutoff used in the error rate filtering step, so that pairwise distances of sequences in the same cluster are consistent with the sequencing error. cluster_fast runs in a single pass, so it is sensitive to input order. Sequences are sorted from lowest to highest quality according to expected error rate; experiment suggests that this order yields better results (Tables A, B, and C in S1 Text).
2. Select and subsample clusters. Clusters with fewer than three members are discarded, because they are too small to de-noise by majority consensus. Clusters with more than 50 members are subsampled to the top 50 with the lowest expected error rate to speed up the multiple sequence alignment step.
3. Align and consensus. MAFFT [45] is used to align each cluster. The consensus sequence of each alignment is computed.

Frame correction
In-frame consensus sequences from all time points are collected into a USEARCH database for frame correction. usearch_global is then used to align each out-of-frame sequence to its top hit. The nucleotide alignment is used to correct incomplete codons: short insertions (1 or 2 base pairs) are discarded, and single deletions are replaced with the aligned base. Sequences with longer insertions or deletions are discarded. All changes are logged, so that the user can identify the sequences that have been corrected.
6. Copy numbers The number of sequences per cluster provides an estimate of the relative abundance of that HQCS in the population. Those numbers are further augmented by adding sequences orphaned by cluster filtering and HQCS dereplication. usearch_global is used to assign each QCS to its nearest HQCS. The number of sequences accrued by each HQCS is interpreted as its copy number.
Alignment sub-pipeline. The HQCSs from all time points are combined into a single file, translated to protein sequences, and aligned using MAFFT. A Python script then transfers the gaps from each aligned protein sequence to the corresponding nucleotide sequences to produce a codon-level nucleotide multiple sequence alignment of all unique variants from all time points.
Analysis sub-pipeline. The analyses used in FLEA take as input the two outputs of the alignment phase: a codon multiple sequence alignment of all unique HQCS sequences from all time points, and their associated copy numbers. These data are used for the following analyses.
1. Time point metrics. HyPhy [46] scripts are used to compute evolutionary metrics (total, dN, and dS divergence and diversity) and phenotypic metrics (protein length, potential Nlinked glycosylation sites, isoelectric point) for each annotated region (e.g., V1, MPER) in the amplicon for each time point.
2. Early consensus. The early consensus is inferred by taking the copy-number-weighted codon consensus of the codon-aligned HQCSs from the earliest time point. By including gaps, this consensus sequence is already aligned with the rest of the multiple sequence alignment. This strategy is acceptable for primary infection studies from single founders with very low early diversity, in which case the consensus sequence should closely match the most recent common ancestor.
3. Reference coordinates. MAFFT is used to assign HXB2 [47] coordinates to the gapped early consensus sequence, which are then transferred to the full multiple sequence alignment.
4. Infer phylogeny. A maximum-likelihood phylogenetic tree is inferred with FastTree2 [48,49] under the general time reversible model. The tree is rooted on the early consensus sequence.

Ancestral sequence reconstruction.
HyPhy is used to infer ancestral sequences at the internal nodes of the phylogeny, using joint maximum likelihood reconstruction and the generalised time-reversible model (GTR) [50].
6. Multidimensional scaling. A distance matrix is computed for all HCQC sequences using the Tamura Nei 93 distance [51]. Metric multidimensional scaling [52] (implemented in scikit-learn [53]) is used to find a two-dimensional embedding of the sequences that approximates their pairwise distances. This embedding with this evolutionary distance is meant to show relationships that cannot be represented in a phylogenetic tree, because of recombination.
8. Position-specific changes. Entropy and Jensen-Shannon divergence are computed for each position in each time point.
The results of these analyses are provided to the user in an interactive web application, described next.

Web application
The FLEA web app is built using modern web design principles. It consists of two parts: a Javascript client-side app, written using the Ember.js [55] framework, and a server-side REST (REpresentational State Transfer) service for serving JSON-formatted data. There are two main benefits to using this decoupled pattern for scientific web applications. First, the clientside code only needs to be downloaded once, at the start of the session. The data are requested from the server and cached as needed. Once everything is loaded, the visualizations run entirely in the browser with no delays for page loads. Second, the REST service may be reused by other apps and third-party tools.
The web app presents the results of the FLEA analysis as a series of interactive visualizations. The report is organized into the following sections.
Multidimensional scaling. A two dimensional embedding of the HQCSs is visualized as a bubble plot, showing changes in population structure over time, as shown in Fig 4. This visualization has been especially useful for investigating populations with superinfection, or with multiple founders, where aggressive recombination between vastly different env variants precludes the use of phylogenies.
Evolutionary trajectory. The evolutionary trajectory viewer plots evolutionary and phenotypic metrics for each time point and multiple regions in the amplicon, giving a high-level Full-Length Envelope Analyzer (FLEA): A tool for longitudinal analysis of viral amplicons overview of population dynamics over time. Fig 5 shows the plot for the entire gp160 region of HIV Env, which is generated with the D3.js plotting library [56].
Sequences. The multiple sequence alignment of all the HQCS sequences is the foundation for all subsequent analyses. It is displayed in the amino acid sequences viewer, which contains a custom alignment browser and an interactive motif dynamics plot, as shown in Fig 6. Protein structure. The protein structure viewer maps evolutionary metrics to an interactive three-dimensional structure of the protein, customized from PDB ID 5FUU, a recently resolved cryo-EM structure [57], and rendered using pv [58]. Missing residues are rendered as spheres which are positioned by Bézier curve interpolation. dN/dS ratios, Jensen Shannon divergence, and entropy may all be mapped to the protein structure, as shown in Fig 7. The same metrics are also plotted in one dimension for each time point, as shown in Fig 8. The protein visualization interacts with the sequence viewer by showing alignment positions and highlighting the residues in the selected sequence motif.
Trees. The tree viewer renders a tree browser with phylotree.js [59], as shown in Fig 9. Leaf nodes are scaled to the copy number of their sequence. The tree zoom level, layout, and coloring is interactively modifiable. Motifs selected in the sequence viewer are mapped to the tree. Ancestral nodes are colored by motif, allowing inferred changes to be tracked through the entire phylogeny.

Results
The entire pipeline was run on HIV env reads from donor P018, which are available from the NCBI Sequence Read Archive under BioProject PRJNA320111, and were sequenced as part of [36] on the RS-II instrument, using the older generation P5/C3 PacBio sequencing chemistry. The full dataset contains 58,468 CCS reads. The reads are split across six time points, which

Results on simulated data
The true sequences and copy numbers are not known for the P018 data. In order to assess the accuracy of our inferred sequence population, we used the HQCSs from a previous FLEA run to simulate a gold standard dataset on which to assess the FLEA pipeline.
The simulation procedure starts with the HQCSs and copy numbers from the FLEA results on P018, then augments them with additional mutated sequences to create a gold standard set of templates. Mutated sequences were added because our clustering strategy may artificially merge similar templates. For each template, noisy reads with a SMRT-style error profile were sampled. Full details of the simulation process appear in the supporting information. These simulated reads were sent through the FLEA pipeline, both with and without frame correction.
The resulting QCS and HQCS sequences were compared to the ground truth using Earth Mover's Distance (EMD), using normalized copy numbers for the population weights and edit distance for the distance matrix. The fully constrained EMD has units that can be directly interpreted as the average change per nucleotide necessary to transform one sequence population into another. We also calculate two variants of EMD for further insight into how well the  To see the effect of sequencing runs of different depths, the experiment was repeated for 300, 1,000, 3,000, and 10,000 reads per time point. The results, which appear in Table 1 The full-length env sequencing protocol yields approximately 10,000 reads per run; the P018 data averaged 9,744 reads per time point. Therefore, these results with n = 10, 000 suggest that FLEA is capable of taking a full sequencing run of CCS reads from a diverse viral population with an average of 9.56 errors per sequence and inferring HQCSs with an average of 1.01 errors per sequence, which corresponds to an average error rate of 0.038%. Moreover, these error rates are mostly caused by low-abundance sequences in both the true population and the

Results on real data: Donor P018
FLEA was run directly on the P018 sequences, and the results are summarized here. The full results of this run are available to view at http://flea.murrell.group/view/P018. Fig 2 shows the number of sequences from the V03 time point that make it to each stage of the quality and consensus pipelines. At three months post infection, the majority amino-acid sequence variant is shared by 52.1% of the population, and the next most common variants accounts for just 8.66%. This relative lack of diversity is consistent with early infection dynamics. By 37 months post infection there is much more diversity: the most common variant accounts for only 3.96% of the population. Donor P018 shows signs of potential N332 glycan specificity, as shown by the motif trajectories in Fig 6. The glycan supersite, centered around N332 in V3, is a common target for broadly-neutralizing antibodies [60] because these sites are often conserved, so mutations in these regions are associated with escape [61]. A year into sampling (V12), mutations 328R and 330H dominate, and the majority of sequences also contain 339N from 22 months (V22) onwards.
The error rates of PacBio CCS sequences are usefully predicted by the QV scores provided by the instrument [36]. We show (Fig. G through Fig. L in S1 Text that the effective number of bases that are corrected in each CCS read (as measured by the difference between that CCS and the HQCS to which it contributes) was extremely well predicted (Spearman's rho from 0.69 to 0.76) by the QV scores. This result is especially encouraging given that our pipeline does not currently exploit these QV scores beyond the initial filtering step. Further, PacBio sequences have higher indel than substitution rates, and this was recapitulated in the number of corrected indels vs substitutions, although this ratio appeared to vary from one time point to the next.

Availability and future directions
A public instance of FLEA is hosted at http://flea.datamonkey.org. The Python source code for the FLEA pipeline can be found at https://github.com/veg/flea-pipeline. The client-side application is available at https://github.com/veg/flea-web-app. A live demo of the P018 results can be found at http://flea.murrell.group/view/P018, with an explanatory page at: http://murrell. group/FLEAexplained/.
The FLEA pipeline analyzes longitudinal full-length env sequences and provides visualizations of the results. Using simulations, we show that FLEA is capable of inferring accurate HIV env consensus sequences and population frequencies. Despite each CCS read containing an average of ten errors, our approach distinguishes variants that differ by as little as one base from an amplicon with high indel variation. It uses those high-quality consensus sequences to generate a codon-aware multiple sequence alignment of all time points, estimate ancestral sequences, infer the phylogenetic tree, and perform many other population-level analyses with high accuracy. These results are presented in a visualization suite that is highly general and applicable to many related sequencing problem.
While we provide a web application that should suffice for sequencing most standard Env samples from HIV-1 group M, we recommend that those who frequently engage in such sequencing, or who wish to sequence less straightforward samples (eg. SIV or SHIV), install FLEA locally. This provides a range of customization and tuning options, such as the filtering parameters and the set of reference sequences.
While our USEARCH-based clustering and consensus strategy for denoising long PacBio amplicons performs well when error rates are < 1%, there is a clear need for more sophisticated long-read de-noising algorithms that exploit the additional depth of lower quality reads that we currently discard. This will be especially beneficial for longer PacBio amplicons, because the CCS read quality distribution degrades with length. For example, while we can currently obtain around 15,000 CCS reads < 1% from a P6/C4 RS-II run of our 2.6kb env amplicon; this read count drops to * 1, 000 for full-length 9kb HIV genomes. Additionally, FLEA does sometimes erroneously collapse sequences from very similar templates, and more sophisticated approaches to amplicon denoising could likely improve upon this.
Both the pipeline and client-side visualizations are under development, with many improvements planned, including a novel clustering algorithm that reduces false positives and a novel consensus algorithm that uses quality scores and performs frame correction. We plan to integrate epitope prediction into the FLEA pipeline and add appropriate visualizations for the case when users have IC 50 values available for their sequences. Finally, FLEA will be expanded to support other amplicons.