QoALa: a comprehensive workflow for viral quasispecies diversity comparison using long-read sequencing data

QoALa: a comprehensive workflow for viral quasispecies diversity comparison using long-read sequencing data | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Method Article QoALa: a comprehensive workflow for viral quasispecies diversity comparison using long-read sequencing data Nakarin Pamornchainavakul, Declan Schroeder, Kimberly VanderWaal This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-4637890/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 27 Apr, 2026 Read the published version in PLOS Computational Biology → Version 1 posted You are reading this latest preprint version Abstract The concept of viral quasispecies refers to a constantly mutating viral population occurring within hosts, which is essential for grasping the micro-evolutionary patterns of viruses. Despite its high error rate, long-read sequencing holds potential for advancing viral quasispecies research by resolving coverage limitations in next-generation sequencing. We introduce a refined workflow, QoALa, implemented in the longreadvqs R package. This workflow begins with nucleotide position-wise noise minimization of read alignments and sample size standardization, and extends to viral quasispecies comparison across related samples. Raw read samples from five studies of different viruses (HCV, HBV, HIV, SARS-CoV-2, and IAV), sequenced by major long-read platforms, were used to evaluate these approaches. The comparative results provide novel insights into intra- and inter-host diversity dynamics in various scenarios and unveil rare haplotypes not reported in the original study, underscoring the versatility and practicality of our methodology. Figures Figure 1 Figure 2 Figure 3 Main The emergence of SARS-CoV-2 in 2019, alongside historical pandemics, emphasizes RNA viruses, with their characteristically high evolutionary rates, as the most concerning group of pathogens. Their remarkable adaptability stands as a primary factor behind their widespread success 1 . Their rapid mutation stemming from a lack of proofreading during genome replication 2 , which serves as a mechanism fostering genetic variability. This variability manifests as a cloud of closely related genetic variants of the virus within the viral population infecting a single individual, which is known as viral quasispecies 3 . Studies elucidating quasispecies dynamics in relation to virus adaptation have significantly enhanced our understanding of RNA virus microevolution, virus-host interactions, and the ability of viruses to adapt to evade host immunity 4–9 . The viral quasispecies concept has been referred to in nearly one thousand virology studies, notably for hepatitis C and B viruses, and HIV 10 , where some viral variants correlate with antiviral drug resistance 11–14 . The advent of deep sequencing technology, alongside improvements in de novo viral genome assembly techniques, has enhanced our ability to detect minority variants within the viral cloud that may hold clinical significance 15–18 . While next-generation sequencing (NGS) has been widely employed for viral quasispecies analysis 19,20 , one of the persisting challenges is the reconstruction of continuous viral quasispecies haplotypes—sets of identical genomic sequences. This challenge arises due to the short length of NGS reads (typically 100–400 bases), which often do not cover the entire targeted gene or genome, necessitating short-read assembly techniques that are often unable to discern between haplotypes that are closely related 17,21–24 . Long-read sequencing technology provides a solution to the coverage limitations of NGS technology in viral quasispecies research 21 , yielding read lengths typically exceeding 1,000 bp. Few analytical workflows or software exist that are tailored for viral quasispecies exploration using long-read data, and these primarily focus on single nucleotide variant (SNV) calling, haplotype reconstruction, or quasispecies profiling on an individual sample basis 25–33 . However, these packages do not fully address the unique challenges or harness the full potential of long-read sequencing for viral quasispecies analysis. One challenge is the higher error rates (10–30%) compared to NGS (~ 1%) 34 , affecting haplotype reconstruction accuracy. Efforts in post-sequencing error correction have spurred the development of numerous tools for improving sequencing reads generated by Oxford Nanopore Technologies (ONT) or Pacific Biosciences (PacBio) 34–38 . However, even with lower error rates, the rapid mutation rate of RNA viruses means that longer reads often differ by only one or two bases, and thus are called a unique haplotype. This ultimately results in a large number of haplotypes that each occur only once within the data, and obscures our ability to discern structure the viral quasispecies. Sequencing read depth and length also significantly influence viral quasispecies diversity measures 21,39,40 . Certain diversity metrics, such as Shannon entropy and mutation frequency, are dependent on sample size or read depth 40,41 . As yet, it is uncertain how such diversity metrics are influenced by sampling depth, particularly when comparing diversity between samples with different depths or down-sampled to the same depth, or even what minimum depth is needed to attain a robust measure of diversity. In addition, the longer the gene or genome length used for quantifying quasispecies, the greater the number of unique haplotypes or singletons due to the increased detection of mutations along its length. This variability may introduce bias, particularly when comparing viral quasispecies profiles across longitudinal or related samples. In response to these challenges, we developed longreadvqs , an R package designed for viral quasispecies comparison using long-read data that can be applied to any RNA virus. This package creates a customizable workflow for long-read noise-minimization and down-sampling as well as for grouping related haplotypes into operational taxonomic units (OTUs). In addition, this package can analyze multiple samples together, identifying common haplotypes or OTUs that recur in different samples. This feature is particularly useful when analyzing longitudinal samples collected from a single individual or multiple individuals that may be epidemiologically linked. Our suggested analytical steps — QoALa: Quasi-species Optimized & Adaptive Long-read Analysis — performed by the key functions of the package, were tested with publicly available long-read data sequenced by both ONT and PacBio technologies from the most studied viruses in the field of viral quasispecies. Results QoALa workflow and example datasets The goals of the QoALa workflow are to standardize and comprehensively compare viral quasispecies and OTU profiles across multiple samples based on user customization. For each comparison, equal length read alignments from the same gene or genomic region of interest (Fig. 1 a) are separately imported from any genome assembly pipeline’s output (see examples in Methods) and noise-minimized using position-wise nucleotide base replacement in the “ vqsassess ” function (Fig. 1 b). The noise-minimization step replaces potential erroneous nucleotides, defined as SNVs with a frequency lower than n% of the total read depth (cut-off percentage), with either the majority base or the dominant haplotype’s base at that position. This cut-off percentage can be determined using prior information such as the sequencing error rate, or by observing the change in the percentage of singleton haplotypes at different percentage cut-offs (see Methods and Fig. 1 a). Additionally, when comparing multiple samples, all alignments must be randomly down-sampled by the “ vqsassess ” function to achieve equivalent read depth (the shallowest depth among all samples is recommended) to reduce bias in diversity measurement caused by sample size disparity. SNV profiles of the samples visualized by the “ snvcompare ” function should be inspected to evaluate the result of noise-minimization (Fig. 1 c). Finally, prepared samples are pooled by the “ vqscompare ” function to ( 1 ) identify common haplotypes, ( 2 ) reclassify haplotypes into OTUs, ( 3 ) visualize diversity profiles and genetic relationships among samples (Fig. 1 d), and ( 4 ) summarize quantitative diversity metrics. We evaluated this workflow using read alignments assembled from publicly available raw reads of the five most studied viruses in the viral quasispecies field, namely HCV, HBV, HIV, SARS-CoV-2, and IAV. These datasets varied not only by sequencing technologies (PacBio or ONT) and sampling scenarios (short- or long-term longitudinal or cross-sectional), but also by the depth of coverage per gene of interest, ranging from 130 reads in the IAV's M gene segment sample (Supplementary Table 6) to 854,147 reads in the SARS-CoV-2's ORF3a gene sample (Supplementary Table 5). Besides, the level of sequencing error varied, estimated from the percentage of singleton haplotypes before noise-minimization. This value ranged from approximately 25% being singleton haplotypes in the SARS-CoV-2's ORF3a gene sample to 100% in the HCV and HIV samples (Supplementary Fig. 1). This considerable variation highlights the necessity for customized long-read data preparation on a case-by-case basis using our package, particularly focusing on two key parameters: the cut-off percentage and sample size (Methods). Effect of down-sampling on diversity metrics While read down-sampling is necessary for standardizing sample size before comparison, it may alter viral quasispecies diversity measures, especially at very low sample sizes. To assess this impact, we sub-sampled reads with replacement and computed nine diversity metrics, ranging from 10,000 down to 25 reads, depending on the original depth of each read alignment sample (Methods). This process was repeated 100 times per sub-sample size. This assessment was conducted twice: once for the alignments that were sampled after noise-minimization, and another for the alignments that were sampled before noise-minimization. As anticipated, diversity metrics directly reliant on read depth [number of haplotypes (H) and mutation frequency at the molecular level (Mfm)], indirectly reliant on read depth [Shannon entropy (H S ), mutation frequency at the molecular level (Mfe), functional attribute diversity (FAD), and nucleotide diversity at the entity level (π e )], or normalized by read depth [normalized Shannon entropy (H SN )] showed varying degrees of fluctuation corresponding to sample size (Supplementary Figs. 2–11). Down-sampling post noise-minimization from high-depth alignments (> 10,000 reads) minimally affected H S and Mfm at sample sizes over 1,000 reads (Supplementary Figs. 2, 4, and 8), while sampling from low-depth alignments (< 1,000 reads) resulted in abrupt changes in these metrics at any sample size (Supplementary Figs. 6 and 10). The Gini-Simpson index (H GS ) and nucleotide diversity (π m ) remained relatively stable across all sample sizes, irrespective of the original alignment depth (Fig. 2 ). Computing metrics from alignments sampled after noise-minimization generally yielded less varied values compared to those sampled before noise-minimization in most instances (Supplementary Figs. 3, 5, 7, 9, and 11). Regarding inter-sample comparisons, variations in metrics due to sample size and sampling strategy marginally affected comparisons in scenarios where the actual diversity measures of samples were either highly similar (Fig. 2 and Supplementary Figs. 6–11) or markedly different (Fig. 2 and Supplementary Figs. 2–5). However, down-sampling to very low sample sizes (< 100 reads) or down-sampling before noise-minimization likely resulted in false equivalences or differences because of the high variation in computed diversity metrics (Fig. 2 and Supplementary Figs. 2–11). Thus, our analysis suggests that diversity metrics from samples with read depths of < 100 should be interpreted with caution. Each dot indicates the mean of one target, boxes the 25th to 75th percentile, lines medians and whiskers extend from minimum to maximum values. Example scenarios of diversity comparison Since a single mutation can differentiate one haplotype from another, and such mutations may introduce noise and obscure the true frequency of different major genetic variants within the quasispecies, solely using haplotype-based diversity metrics to summarize the viral quasispecies profile might hinder our ability to observe the continuous dynamics of the virus across samples. To illustrate this, we compared related viral samples in various scenarios using both traditional haplotyping based on a strict mutational profile and a novel OTU assignment based on genetic distance (Methods). In the latter, haplotypes are clustered by relatedness into larger operational taxonomic units (OTUs). This allows the user to better comprehend the structure of the quasispecies, i.e., the evolving genetic diversity surrounding several main variants within the viral cloud. The studies on HIV and SARS-CoV-2 served as examples of short-term longitudinal sampling scenarios. Within a few hours of each other, env genes of four HIV samples infecting T cells in vitro exhibited numerous haplotypes, among which a common dominant haplotype was not found (Fig. 3 a). When OTU clustering was applied to the same samples, the members and proportions of all four samples were relatively similar (Fig. 3 b), with the 186 haplotypes being optimally clustered into 10 OTUs (Fig. 3 c/d). Examining five SARS-CoV-2’s S gene samples taken within a nine-day period from a single patient, the dominant haplotype and OTU occupied over half of the virus population in each sample (Fig. 3 a/b), with minor haplotypes appearing as outliers (Fig. 3 c/d). It is also worth noting that the sporadic appearance of many minority variants in the day 13 and 17 SARS-CoV-2 samples (Fig. 3 a/b) was likely due to noise rising from incomparable read depths in the original alignments, which were over ten times lower than the other three samples (Supplementary Table 5). As a model for long-term longitudinal sample comparison, studies on HCV and HBV were utilized. Analysis of the NS4 coding region of HCV revealed that viral haplotypes and their proportions underwent sudden changes after treatment. Interestingly, certain major haplotypes detected pre-treatment reappeared 18 months post-treatment (Fig. 3 a). Concurrently, OTU clustering provided better insights into the dynamics of HCV across three sampling points by grouping genetically closely related minor haplotypes into common clusters observed at varying proportions across all time points (Fig. 3 b/c/d). The within-host virus population dynamics were also clearly illustrated by HBV's S gene. Both haplotype and OTU classifications underscored the gradual turnover of predominant haplotypes from timepoint one to four over a 30-month period, and a completely different population emerged at timepoint five, occurring over 80 months after the initial sampling (Fig. 3 a/b). Lastly, in the scenario of the outbreak represented by IAV's M gene segment samples, viruses were collected cross-sectionally among patients in two wards (B and C) within the same hospital. At the haplotype level, only two mostly identical sets of virus populations were found in more than one patient (B1, B2 and C1, C3), without other prominent links between cases or wards (Fig. 3 a). However, upon regrouping them into new OTUs, we discovered that all patients in ward B were partly infected with the same genetically similar sub-population, while all patients in ward C were infected with a similar set of OTUs. Furthermore, the predominant OTU found in two ward B patients (B1, B2) was also present as a minority group in all ward C patients, suggesting potential inter-ward transmission of IAV (Fig. 3 b/c/d). The comparative profiles of viral quasispecies, whether analyzed at the haplotype or OTU levels, based on other genes from the same viral samples (including HBV’s P gene, HIV’s gag gene, and SARS-CoV-2’s ORF3a gene), exhibited visual similarities to the chosen genes depicted in Fig. 3 (Supplementary Fig. 12 to 14). This suggests the absence of linkage disequilibrium between the observed genes of these three viruses, indicating that they may adequately represent the population dynamics at the whole-genome level. Interestingly, one comparable result between the original and our studies is the number and frequency of haplotypes found in the day 11 and day 15 SARS-CoV-2 samples. We captured a similar quasispecies diversity in the S gene samples and discovered two SNV positions in the ORF3a gene that created two rare haplotypes not reported in the original study (Supplementary Table. 7 and Ko et al., 2021 46 ). Discussion The longreadvqs package was developed as a tool for both quantitative and qualitative analysis of viral quasispecies diversity, addressing the challenges posed by varied error rates and read depths inherent in long-read sequencing technologies as well as noise introduced by rare SNVs in the alignment. During the testing phase, we carefully considered these variations across diverse sampling scenarios and multiple viral species, ensuring the tool’s general applicability. The strengths of this tool lie in its ability to customize parameter settings based on prior information and in providing comprehensive visualizations that elucidate the dynamics of virus microevolution across multiple samples. The analytical findings derived from the QoALa workflow presented in this study provide unique perspectives and insights not previously explored in the original research. For instance, we utilized raw read data to reconstruct within-host virus haplotypes and OTUs, emphasizing an overview of population structure changes over time or among outbreak cases. In contrast, the original studies focused on specific aspects, such as mutation spectra 42 , deletions 43 , gene splicing 44 , or phylogenetic analysis of consensus sequences 45 . To build read alignments, we employed a uniform genome assembly workflow for all raw read samples, which was not tailored or reproduced from the original protocols created for each dataset. Furthermore, in the comparison step, we fixed the number of sample sizes and k-means clusters for all scenarios. Hence, discrepancies in details between the original findings and ours may have arisen, which we acknowledge as a limitation on our part. Nevertheless, it's essential to note that the objective of our study does not revolve around delivering new findings from each dataset but rather demonstrating the usage and versatility of this package. Previous studies have established a robust foundation in measuring viral quasispecies diversity and developing analytical workflows for NGS data over the past decade 20,40,41,47 . However, much of this research cannot be directly applied to long-read sequencing data. Long-read quasispecies analyses are hindered by the significantly higher level of noise present in long-read data driven by the higher likelihood that longer reads will have at least one SNV across their length. For NGS short-read data, it is recommended to employ techniques such as rarefaction, resampling, and fringe trimming based on haplotype frequency to estimate appropriate sample sizes and minimize bias from unbalanced samples before comparing diversity metrics 41 . The presence of a misrepresented large proportion of singleton haplotypes (Supplementary Fig. 1), compounded by the high sequencing error rate of current long-read technologies, hinders us from following such approaches. Instead of excluding low-frequency haplotypes as fringe trimming does, we opt to retain all reads and replace low-frequency SNVs with the majority base at their respective positions. This approach offers two primary benefits. Firstly, low-frequency haplotypes, comprising both true mutations and errors, are not eliminated from the alignment; rather, potential erroneous SNVs within the haplotype are smoothed to the mode. This advantage is particularly crucial when comparing shallow depth read alignments, where removing reads belonging to low-frequency haplotypes would severely reduce sample size. Secondly, the false haplotype diversity resulting from long-read errors are mitigated by consolidating previously low-frequency haplotypes into larger haplotype groups (OTUs) that share common SNVs, which are more likely to represent true mutations and provide more clarity on the structured diversity existing within quasispecies. Despite these advantages, many factors need to be considered in selecting the cut-off percentage for position-wise noise-minimization. Factors such as gene mutation rate, sequencing error rate, and the trade-off between retaining versus unifying rare haplotypes are critical. Especially when dealing with a lengthy alignment where singleton haplotypes are more common, one may have to decide whether to prioritize retaining most mutation details, thus making it difficult to quantify overall population complexity, or to focus on the broader picture of viral cloud diversity while potentially discarding some details that may or may not be significant. Standardizing unbalanced sample sizes by down-sampling can introduce bias in diversity comparisons, especially when integrating noise-minimization. Since sampling techniques used for NGS data 41 are not appropriate for long-read data, we redesigned sampling strategies by integrating the concept of rarefaction and repeat sampling, taking into account the order of sampling before or after noise-minimization to simulate possible options for analysis. Here, we demonstrate the sensitivity of diversity metrics across our example scenarios. Based on this analysis, down-sampling after noise-minimization is our recommended technique. However, repeated down-sampling to observe the distribution of metrics at different sample sizes, at least between the largest and smallest samples in the dataset, is encouraged to better understand the sensitivities of diversity metrics to sampling depth in a particular project. Ultimately, we suggest interpreting quantitative metrics comparisons alongside qualitative profiles that uncover shared viral variants between related samples, either at the haplotype or OTU levels, to better illuminate the dynamics of viral quasispecies, as made possible by the QoALa workflow implemented in the longreadvqs package. Methods Sequencing reads data To evaluate our analytic approaches, we gathered raw sequencing reads from five projects focusing on whole genome sequence (WGS) of five different viruses: hepatitis C virus (HCV), hepatitis B virus (HBV), human immunodeficiency virus (HIV), severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), and influenza A virus (IAV), which are among the most studied species in the field of viral quasispecies. These reads were generated using either PacBio Single Molecule, Real-Time (SMRT) sequencing or Oxford Nanopore Technologies (ONT) and were sourced from the National Center for Biotechnology Information’s Sequence Read Archive (NCBI SRA). Specifically, we selected projects that obtained either longitudinal viral samples or cross-sectional samples from a single outbreak, where the genetic relationships were suitable for comparison at the quasispecies level. PacBio read samples were obtained from three projects. First, three HCV samples from a chronic hepatitis patient, were collected before and after treatment over a span of 18 months 42 . Second, five samples were collected at various time points within a period of 100 months from an untreated HBV-infected patient 43 . Third, five longitudinal SARS-CoV-2 samples were collected between 8 to 17 days after clinical onset from a single patient 46 . Data from the remaining two projects were generated using ONT, which included four samples from HIV-infected T cells collected up to 24 hours post-infection 44 , and seven samples from different patients infected by IAV during a nosocomial outbreak in the same hospital 45 . The NCBI accession numbers of a total of 24 raw read samples and five projects were listed in Supplementary Table 1. Genome assembly and read alignment generation We conducted comparable genome assembly steps for both PacBio and ONT raw reads. Initially, we employed the default settings of HiFiAdapterFilt v3.0.1 48 and Porechop v0.2.4 49 to trim PacBio and ONT adapter sequences, respectively. Subsequently, poor-quality PacBio reads—those shorter than 1,000 bp or constituting the worst 10% of reads based on the final score—were filtered out using Filtlong v0.2.1 50 . Similarly, ONT reads with an average quality score < 7 were removed using NanoFilt v2.8.0 51 . We then aligned the trimmed and filtered reads against the NCBI reference genomic sequence of each viral species (Supplementary Table 1) using minimap2 v2.26’s presets 52 tailored for each sequencing technology ("map-hifi" for PacBio and "map-ont" for ONT). Subsequently, to prepare the resulting read alignments for downstream analyses, SAM files from read mapping were converted to FASTA format and fragmented into genes or coding regions based on the genomic annotation of each reference sequence, utilizing SeqKit v2.7.0 53 . It was ensured that the reads retained in the alignment covered the full length of the particular gene or region (Fig. 1 a). The depth and length of the final alignments were documented in Supplementary Tables 2 to 6. Sequencing read noise-minimization A main function of the longreadvqs package, “ vqsassess ,” aims to minimize noise from sequencing errors, artifacts, or rare mutations by performing position-wise nucleotide base replacement. We hypothesize that any SNV in a read alignment with a base frequency at each position less than the specified cut-off percentage is likely an artifact that should be replaced with the majority base of that position (Fig. 1 b). The cut-off percentage can be determined based on prior knowledge, such as the estimated sequencing error rate of the technology used or documented mutation or evolutionary rates of the studied virus. However, the selection of the cut-off for noise-minimization can be further guided by the “ pctopt ” function, which demonstrates the percentage of singleton haplotypes (haplotypes with only a single read member) in the alignment. The percentage of singleton haplotypes decreases as the cut-off percentage increases, as low-frequency SNVs are replaced, creating more groups of identical reads (Fig. 1 a). The normal range for the percentage of singleton haplotypes found in previous studies using more accurate technologies like Sanger or NGS sequencing can serve as a reference for this step. Alternatively, the plot showing changes in the percentage of singleton haplotypes can be used to identify the cut-off value where the percentage of singleton haplotypes ceases decreasing, indicating a plateau (Fig. 1 a). An alternative option for noise-minimization involves replacing potential erroneous low-frequency SNV bases with the base of the dominant haplotype. However, this method is not recommended, since either the base used for replacement could also be an error, or the “dominant” haplotype before noise-minimization may be relatively rare or non-existent (100% singleton haplotypes). In some cases, a single cut-off percentage cannot be generalized for the entire gene or region of interest because the amount of noise may vary throughout the sequence length. For example, low ONT base calling accuracy is often reported in homopolymer regions (continuous identical bases) 54,55 or soft clipping (unaligned regions) may persist at the ends of reads after mapping in both technologies 56,57 . The “ snvcompare ” function, which visualizes SNV distribution across sequence length between different samples (Fig. 1 c), can help identify noise-rich regions indicated by the overaccumulation of SNV positions. Subsequently, the “ vqscustompct ” function can be applied to readjust the cutoff percentage for noise-minimization at specific region(s) (Fig. 1 b). To validate the noise-minimization workflow, we extracted read alignments from specific genes or regions of five exemplary viral species. These alignments were chosen based on their sequence length, with each exceeding or closely approximating 1,000 bases, and their depth of coverage surpassing 100 reads, while ensuring that the soft-clipped region, if present, comprised less than 25% of the total length. Utilizing the " pctopt " function, we determined the optimal cut-off percentage for minimizing errors, aiming to achieve nearly 0% singleton haplotypes, or identifying the point where the percentage of singleton haplotypes in the median sample ceased decreasing, indicating a plateau phase (Fig. 1 a). The selected alignments (with their respective cut-off percentages) comprised the NS4 coding region (15%) of HCV, the P (5%) and S (5%) genes of HBV, the env (15%) and gag (22%) genes of HIV, the S (1%) and ORF3a (1%) genes of SARS-CoV-2, and the M gene segment (10%) of IAV (Supplementary Fig. 1). Down-sampling and diversity metric sensitivity For the demonstration of our package usage, we specified basic settings for the " vqsassess " function. After the noise-minimization step described above, we then down-sampled to either the size of the shallowest depth sample or to 1,000 reads (for samples with a depth over 10,000 reads), for every sample. These steps were taken to prepare the samples for the final quasispecies comparison within the same viral species or project. However, it's important to note that read depth or sample size significantly influences some quantitative metrics used for quasispecies diversity measures 40,41 . Aggressive down-sampling may result in information loss and misinterpretation. To evaluate the impact of down-sampling on between-sample diversity comparison, we conducted a sensitivity analysis of nine diversity metrics as fully described by Gregori et al. 40 , including the number of haplotypes (H), the Shannon entropy (H S ), the normalized Shannon entropy (H SN ), the Gini-Simpson index (H GS ), the functional attribute diversity (FAD), the mutation frequency at the entity level (Mfe), the nucleotide diversity at the entity level (π e ), the mutation frequency at the molecular level (Mfm), and the nucleotide diversity (π m ). We utilized two to four samples of one gene of interest per species for the analysis. Using the QSutils v1.18.0 58 embedded in our package, all nine metrics were computed from both unsampled and down-sampled read alignments. For samples with a depth over 10,000 reads (HCV, HBV, and SARS-CoV-2), down-sample sizes started from 10,000 reads and were gradually halved until the final size reached 78 reads. For samples with low depth (HIV and IAV), sample sizes were set at 300, 150, 100, 50, and 25 reads. Random down-sampling with replacement was repeated 100 times for each sample size after noise-minimization. The same repeated sampling strategy was also implemented before noise-minimization. Distributions of each metric were visually compared between different sample sizes and approaches (down-sampling after versus before noise-minimization). Viral quasispecies comparison Once the within-species samples were noise-minimized and standardized to equal depths, they were ready for comparison using the " vqscompare " function (Fig. 1 d). This function aggregates the prepared read alignments of the listed samples, initially identifying shared haplotypes between them, and then visualizes the proportion of unique haplotypes for each sample as a color-coded bar plot. In addition, it reclassifies haplotypes into new operational taxonomic units (OTUs) based on the genetic distance matrix, utilizing the " dist.dna " function in ape v5.7.1 59 , from the SNV alignment extracted from the pooled read alignment. In detail, the distance matrix is transformed into dissimilarity coordinates, which are subsequently clustered into new OTUs using classical multidimensional scaling and k-means clustering, respectively, through the " cmdscale " and " kmeans " functions in stats v4.3.1 60 . The number of OTUs or clusters of genetically closely related haplotypes must be customized by specifying the number of clusters (k). The proportions of OTUs are also depicted in a color-coded bar plot. The clustering scheme of OTUs, along with major haplotypes within them, is illustrated with corresponding colors in multidimensional scale (MDS) plots (Fig. 1 d). All resulting plots are generated using ggplot2 v3.4.4 61 . Apart from the comparative plots, the " vqscompare " function offers several other valuable outputs for in-depth investigation. These include noise-minimized down-sampled read and SNV alignments for all samples, along with classified haplotype and OTU identifications. Additionally, the nine quasispecies diversity metrics for each sample, computed using QSutils v1.18.0 58 , are tabulated. Another set of metrics are calculated based on consensus reads and read frequencies of each OTU, rather than haplotypes, to simplify the diversity quantification in a larger scale. To exemplify such comparative features for our examples, noise-minimized down-sampled alignments of each selected gene and viral species, prepared by the " vqsassess " function, are listed as input for individual " vqscompare " analyses. The number of clusters (k) for OTU classification via k-means clustering was set to 10 for every run. The resulting summary plot was used to illustrate the dynamics of viral quasispecies diversity across different example scenarios, ranging from short-term (hours to days in HIV and SARS-CoV-2 datasets) to long-term (months to years in HCV and HBV datasets) longitudinal samples, as well as samples from the same outbreak cohort (IAV dataset). Declarations Contributions N.P. and K.V. conceived, initiated, and designed the study. N.P. developed the longreadvqs package, collected raw read data, implemented the QoALa workflow, conducted all analyses, and drafted the manuscript. K.V. and D.C.S. supervised the development of the package and workflow. All authors contributed to the final version of the paper. Author Contribution N.P. and K.V. conceived, initiated, and designed the study. N.P. developed the longreadvqs package, collected raw read data, implemented the QoALa workflow, conducted all analyses, and drafted the manuscript. K.V. and D.C.S. supervised the development of the package and workflow. All authors contributed to the final version of the paper. Acknowledgement We thank J. Gregori and M. Guerrero-Murillo from VHIR Vall d’Hebron Research Institute for their support in troubleshooting the dependency of our package. We thank TGS. Williams from Guy's and St Thomas' NHS Foundation Trust for providing additional information of the IAV study. Individually, we thank D. Makau and J. Baker for testing the package. Data Availability All raw read sequencing data used as example datasets for this study are publicly available at the National Center for Biotechnology Information’s Sequence Read Archive (NCBI SRA), as well as being part of previously published publications. Accession numbers and corresponding publications were tabulated in Supplementary Table 1. Code availability All codes used in this study are available as an R package longreadvqs and accessible via the CRAN repository at https://cran.r-project.org/web/packages/longreadvqs/index.html References Carrasco-Hernandez R, Jácome R, López Vidal Y. Ponce De León, S. Are RNA Viruses Candidate Agents for the Next Global Pandemic? A Review. ILAR J. 2017;58:343–58. Steinhauer DA, Domingo E, Holland JJ. Lack of evidence for proofreading mechanisms associated with an RNA virus polymerase. Gene. 1992;122:281–8. Domingo E, Sheldon J, Perales C. Viral Quasispecies Evolution. Microbiol Mol Biol Rev. 2012;76:159–216. Domingo E. Quasispecies Structure and Persistence of RNA Viruses. Emerg Infect Dis. 1998;4:521–7. Mandary M, Poh. Impact of RNA Virus Evolution on Quasispecies Formation and Virulence. IJMS. 2019;20:4657. Vignuzzi M, Stone JK, Arnold JJ, Cameron CE, Andino R. Quasispecies diversity determines pathogenesis through cooperative interactions in a viral population. Nature. 2006;439:344–8. Woo H-J, Reifman J. A quantitative quasispecies theory-based model of virus escape mutation under immune selection. Proc. Natl. Acad. Sci. U.S.A. 109, 12980–12985 (2012). Domingo E, García-Crespo C, Perales C. Historical Perspective on the Discovery of the Quasispecies Concept. Annu Rev Virol. 2021;8:51–72. Lauring AS. Within-Host Viral Diversity: A Window into Viral Evolution. Annu Rev Virol. 2020;7:63–81. PubMed. PubMed https://pubmed.ncbi.nlm.nih.gov/ . Metzner K. The significance of minority drug-resistant quasispecies. In: Geretti AM, editor. Antiretroviral Resistance in Clinical Practice. London: Mediscript; 2006. Monaco DC, Zapata L, Hunter E, Salomon H, Dilernia DA. Resistance profile of HIV-1 quasispecies in patients under treatment failure using single molecule, real-time sequencing. AIDS. 2020;34:2201. Perales C. Quasispecies dynamics and clinical significance of hepatitis C virus (HCV) antiviral resistance. Int J Antimicrob Agents. 2020;56:105562. Kai Y, et al. Baseline quasispecies selection and novel mutations contribute to emerging resistance-associated substitutions in hepatitis C virus after direct-acting antiviral treatment. Sci Rep. 2017;7:41660. Margeridon-Thermet S, et al. Ultra‐Deep Pyrosequencing of Hepatitis B Virus Quasispecies from Nucleoside and Nucleotide Reverse‐Transcriptase Inhibitor (NRTI)–Treated Patients and NRTI‐Naive Patients. J INFECT DIS. 2009;199:1275–85. Rozera G, et al. Massively parallel pyrosequencing highlights minority variants in the HIV-1 env quasispecies deriving from lymphomonocyte sub-populations. Retrovirology. 2009;6:15. Baaijens JA, Aabidine AZE, Rivals E, Schönhuth A. De novo assembly of viral quasispecies using overlap graphs. Genome Res. 2017;27:835–48. Fritz A, et al. Haploflow: strain-resolved de novo assembly of viral genomes. Genome Biol. 2021;22:212. Houldcroft CJ, Beale MA, Breuer J. Clinical and biological insights from viral genome sequencing. Nat Rev Microbiol. 2017;15:183–92. Lu I-N, Muller CP, He FQ. Applying next-generation sequencing to unravel the mutational landscape in viral quasispecies. Virus Res. 2020;283:197963. Posada-Cespedes S, Seifert D, Beerenwinkel N. Recent advances in inferring viral diversity from high-throughput sequencing data. Virus Res. 2017;239:17–32. Huang A, Kantor R, DeLong A, Schreier L, Istrail S, QColors. An algorithm for conservative viral quasispecies reconstruction from short and non-contiguous next generation sequencing reads. In Silico Biol. 2012;11:193–201. Hong LZ, et al. BAsE-Seq: a method for obtaining long viral haplotypes from short sequence reads. Genome Biol. 2014;15:517. Mardis E. R. DNA sequencing technologies: 2006–2016. Nat Protoc. 2017;12:213–8. Dilernia DA, et al. Multiplexed highly-accurate DNA sequencing of closely-related HIV-1 variants using continuous long reads from single molecule, real-time sequencing. Nucleic Acids Res. 2015;43:e129–129. Huang DW. Towards Better Precision Medicine: PacBio Single-Molecule Long Reads Resolve the Interpretation of HIV Drug Resistant Mutation Profiles at Explicit Quasispecies (Haplotype) Level. J Data Min Genomics Proteom 07, (2016). Dudouet P, et al. SARS-CoV-2 quasi-species analysis from patients with persistent nasopharyngeal shedding. Sci Rep. 2022;12:18721. Artyomenko A, et al. Long Single-Molecule Reads Can Resolve the Complexity of the Influenza Virus Composed of Rare, Closely Related Mutant Variants. J Comput Biol. 2017;24:558–70. Jiao X, et al. QuasiSeq: profiling viral quasispecies via self-tuning spectral clustering with PacBio long sequencing reads. Bioinformatics. 2022;38:3192–9. Link RW, et al. HIV-Quasipore: A Suite of HIV-1-Specific Nanopore Basecallers Designed to Enhance Viral Quasispecies Detection. Front Virol. 2022;2:858375. Luo X, Kang X, Schönhuth A. Strainline: full-length de novo viral haplotype reconstruction from noisy long reads. Genome Biol. 2022;23:29. Ng TT-L, et al. Long-Read Sequencing with Hierarchical Clustering for Antiretroviral Resistance Profiling of Mixed Human Immunodeficiency Virus Quasispecies. Clin Chem. 2023;69:1174–85. Su J, Li S, Zheng Z, Lam T-W, Luo R. ClusterV-Web: a user-friendly tool for profiling HIV quasispecies and generating drug resistance reports from nanopore long-read data. Bioinf Adv. 2024;4:vbae006. Morisse P, Marchet C, Limasset A, Lecroq T, Lefebvre A. Scalable long read self-correction and assembly polishing with multiple sequence alignment. Sci Rep. 2021;11:761. Sahlin K, Medvedev P. Error correction enables use of Oxford Nanopore technology for reference-free transcriptome analysis. Nat Commun. 2021;12:2. Wang L, Qu L, Yang L, Wang Y, Zhu H. NanoReviser: An Error-Correction Tool for Nanopore Sequencing Based on a Deep Learning Algorithm. Front Genet. 2020. 10.3389/fgene.2020.00900 . Xiao C-L, et al. MECAT: fast mapping, error correction, and de novo assembly for single-molecule sequencing reads. Nat Methods. 2017;14:1072–4. Salmela L, Walve R, Rivals E, Ukkonen E. Accurate self-correction of errors in long reads using de Bruijn graphs. Bioinformatics. 2017;33:799–806. Zagordi O, Däumer M, Beisel C, Beerenwinkel N. Read length versus Depth of Coverage for Viral Quasispecies Reconstruction. PLoS ONE. 2012;7:e47046. Gregori J, et al. Viral quasispecies complexity measures. Virology. 2016;493:227–37. Gregori J, et al. Inference with viral quasispecies diversity indices: clonal and NGS approaches. Bioinformatics. 2014;30:1104–11. Nakamura F, et al. Mutational spectrum of hepatitis C virus in patients with chronic hepatitis C determined by single molecule real-time sequencing. Sci Rep. 2022;12:7083. Arasawa S, et al. Evolutional transition of HBV genome during the persistent infection determined by single-molecule real-time sequencing. Hepatol Commun. 2023;7:e0047–0047. Nguyen Quang N, et al. Dynamic nanopore long-read sequencing analysis of HIV-1 splicing events during the early steps of infection. Retrovirology. 2020;17:25. Williams TGS et al. Feasibility and clinical utility of local rapid Nanopore influenza A virus whole genome sequencing for integrated outbreak management, genotypic resistance detection and timely surveillance. Microb Genomics 9, (2023). Ko SH, et al. High-throughput, single-copy sequencing reveals SARS-CoV-2 spike variants coincident with mounting humoral immunity during acute COVID-19. PLoS Pathog. 2021;17:e1009431. Knyazev S, Hughes L, Skums P, Zelikovsky A. Epidemiological data analysis of viral quasispecies in the next-generation sequencing era. Brief Bioinform. 2021;22:96–108. Sim SB, Corpuz RL, Simmonds TJ, Geib SM. HiFiAdapterFilt, a memory efficient read processing pipeline, prevents occurrence of adapter sequence in PacBio HiFi reads and their negative impacts on genome assembly. BMC Genomics. 2022;23:157. Wick RR, Judd LM, Gorrie CL, Holt KE. Completing bacterial genome assemblies with multiplex MinION sequencing. Microb Genomics 3, (2017). Wick R. rrwick/Filtlong. (2024). De Coster W, D’Hert S, Schultz DT, Cruts M, Van Broeckhoven C. NanoPack: visualizing and processing long-read sequencing data. Bioinformatics. 2018;34:2666–9. Li H. New strategies to improve minimap2 alignment accuracy. Bioinformatics. 2021;37:4572–4. Shen W, Le S, Li Y, Hu F, SeqKit. A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation. PLoS ONE. 2016;11:e0163962. Sarkozy P, Jobbágy Á, Antal P. Calling Homopolymer Stretches from Raw Nanopore Reads by Analyzing k-mer Dwell Times. in EMBEC & NBC 2017 (eds. Eskola, H., Väisänen, O., Viik, J. & Hyttinen, J.) vol. 65 241–244Springer Singapore, Singapore, (2018). Huang Y-T, Liu P-Y, Shih P-W. Homopolish: a method for the removal of systematic errors in nanopore sequencing by homologous polishing. Genome Biol. 2021;22:95. Delahaye C, Nicolas J. Sequencing DNA with nanopores: Troubles and biases. PLoS ONE. 2021;16:e0257521. Zhang S-J, et al. Isoform Evolution in Primates through Independent Combination of Alternative RNA Processing Events. Mol Biol Evol. 2017;34:2453–68. Guerrero-Murillo M. QSutils. [object Object] https://doi.org/10.18129/B9.BIOC.QSUTILS (2018). Paradis E, Schliep K. Ape 5.0: An environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics. 2019;35:526–8. R Core Team. R: A language and environment for statistical computing. R Foundation Stat Comput (2019). Ginestet C. ggplot2: Elegant Graphics for Data Analysis. J Royal Stat Society: Ser (Statistics Society) 174, (2011). Additional Declarations No competing interests reported. Supplementary Files SupplementaryFigures.docx SupplementaryTables.xlsx Cite Share Download PDF Status: Published Journal Publication published 27 Apr, 2026 Read the published version in PLOS Computational Biology → Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-4637890","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Method Article","associatedPublications":[],"authors":[{"id":323420000,"identity":"a5047979-382d-471a-8fb8-eb98cb9a7ea9","order_by":0,"name":"Nakarin Pamornchainavakul","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAtUlEQVRIiWNgGAWjYNACAwkefgiLmVgtFTYykg2kaTmTZmNwgFgt5uynEx9Xth3mMb6RnfaAocI6sYGQFsue3M2GZ4FazG7kbjdgOJNOWIvBgdxtko0QLdskGNsOE6Hl/FuIFuMZIC3/iNECNFyy4Uwaj4EESEsDEVosZ7zdbNhQYcMjcebtNomEY+nGBLWY8+dufNhgIGHP3w605UONtSxhh6HwEggpx9QyCkbBKBgFowAbAAC1HT+JjLLAbgAAAABJRU5ErkJggg==","orcid":"","institution":"University of Minnesota","correspondingAuthor":true,"prefix":"","firstName":"Nakarin","middleName":"","lastName":"Pamornchainavakul","suffix":""},{"id":323420002,"identity":"c0e4e9d3-991a-432b-af5e-8773661d85fd","order_by":1,"name":"Declan Schroeder","email":"","orcid":"","institution":"University of Minnesota","correspondingAuthor":false,"prefix":"","firstName":"Declan","middleName":"","lastName":"Schroeder","suffix":""},{"id":323420004,"identity":"c0f28896-1705-456f-9293-32dbd60cce64","order_by":2,"name":"Kimberly VanderWaal","email":"","orcid":"","institution":"University of Minnesota","correspondingAuthor":false,"prefix":"","firstName":"Kimberly","middleName":"","lastName":"VanderWaal","suffix":""}],"badges":[],"createdAt":"2024-06-25 16:36:15","currentVersionCode":1,"declarations":{"humanSubjects":false,"vertebrateSubjects":false,"conflictsOfInterestStatement":false,"humanSubjectEthicalGuidelines":false,"humanSubjectConsent":false,"humanSubjectClinicalTrial":false,"humanSubjectCaseReport":false,"vertebrateSubjectEthicalGuidelines":false},"doi":"10.21203/rs.3.rs-4637890/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-4637890/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1371/journal.pcbi.1014208","type":"published","date":"2026-04-28T00:00:00+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":60573890,"identity":"097e7d7b-fff0-435c-8dbf-5745b106f73f","added_by":"auto","created_at":"2024-07-18 10:07:33","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":784821,"visible":true,"origin":"","legend":"\u003cp\u003eQoALa workflow and main functions of \u003cem\u003elongreadvqs\u003c/em\u003e package.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ea\u003c/strong\u003e, Sequencing read alignment that completely covers a gene of interest is used as input for subsequent analyses in the \u003cem\u003elongreadvqs\u003c/em\u003e package. The percentage cut-off, which is used to define low-frequency SNVs potentially as noise, can be estimated from prior knowledge or the percentage of singleton haplotypes in the alignment. \u003cstrong\u003eb\u003c/strong\u003e, Sample preparation involves noise-minimization, where nucleotide bases with frequencies lower than the percentage cut-off are replaced with the majority base (recommended) or the base of the dominant haplotype at the respective position, followed by read down-sampling for sample size standardization. The percentage cut-off can also be customized for error-rich regions. \u003cstrong\u003ec\u003c/strong\u003e, The SNV profile after sample preparation should be compared between samples to visually check for remaining potential noise in the alignment. \u003cstrong\u003ed\u003c/strong\u003e, All prepared samples are pooled to identify shared haplotypes and reclassify them into OTUs based on genetic distance. Green texts indicate key functions (italicized) and arguments in the \u003cem\u003elongreadvqs\u003c/em\u003e package. Solid boxes represent main steps in the QoALa workflow, while dashed boxes denote optional steps.\u003c/p\u003e","description":"","filename":"floatimage1.png","url":"https://assets-eu.researchsquare.com/files/rs-4637890/v1/35c40593803a6f197dee7356.png"},{"id":60573889,"identity":"fcd63fb3-f28c-4a9b-9b88-1eeeb56ab873","added_by":"auto","created_at":"2024-07-18 10:07:33","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":1264041,"visible":true,"origin":"","legend":"\u003cp\u003eEffect of down-sampling on key diversity metrics for five viruses (virus:gene).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ea\u003c/strong\u003e, Distribution of the Gini-Simpson index (H\u003csub\u003eGS\u003c/sub\u003e) for each viral species: gene samples computed at different sample sizes (100 repeated random samplings with replacement for each size) after (top) and before (bottom) noise-minimization. \u003cstrong\u003eb, \u003c/strong\u003eDistribution of nucleotide diversity (π\u003csub\u003em\u003c/sub\u003e) with a layout similar to that of panel a. Bottom legends indicate the names of samples for each virus dataset. For all plots, horizontal lines indicate medians, boxes the 25th to 75th percentile, whiskers the lowest and highest values excluding outliers, and dots the outliers.\u003c/p\u003e","description":"","filename":"floatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-4637890/v1/017ac0b9fdea2070e6288238.png"},{"id":60573080,"identity":"400c927e-f6b3-4d55-a8f3-5514a6713750","added_by":"auto","created_at":"2024-07-18 09:59:33","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":1087404,"visible":true,"origin":"","legend":"\u003cp\u003eExamples of viral quasispecies and OTU comparisons for five viruses (virus:gene).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ea\u003c/strong\u003e, Stacked bar chart presenting the proportion of each unique haplotype (same color) identified across different samples within each viral species: gene dataset. \u003cstrong\u003eb\u003c/strong\u003e, Stacked bar chart presenting the proportion of new OTUs, each consisting of genetically closely related haplotypes. \u003cstrong\u003ec\u003c/strong\u003e, Multidimensional scaling (MDS) plot of haplotypes' pairwise SNV distance (dot), showing the five largest OTUs (circle) with a color scheme corresponding to plot b. \u003cstrong\u003ed\u003c/strong\u003e, MDS plot of all OTUs.\u003c/p\u003e","description":"","filename":"floatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-4637890/v1/eb811b2f1e74765b9c0ab9e8.png"},{"id":108212092,"identity":"09579fc1-4732-4a4e-8878-64368c4a4ede","added_by":"auto","created_at":"2026-04-30 13:56:55","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":2965626,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-4637890/v1/b39f3ac2-64a2-4c20-a44d-4491213a2872.pdf"},{"id":60573083,"identity":"6ef16039-107e-4f29-bbce-28cfff0a938a","added_by":"auto","created_at":"2024-07-18 09:59:33","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":5682451,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryFigures.docx","url":"https://assets-eu.researchsquare.com/files/rs-4637890/v1/1496ce898062de101212af3d.docx"},{"id":60573079,"identity":"068e1ca2-4065-4cb2-ac3e-05d63fc270ad","added_by":"auto","created_at":"2024-07-18 09:59:33","extension":"xlsx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":20233,"visible":true,"origin":"","legend":"","description":"","filename":"SupplementaryTables.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-4637890/v1/41566bc2a99f097917ab7cd4.xlsx"}],"financialInterests":"No competing interests reported.","formattedTitle":"QoALa: a comprehensive workflow for viral quasispecies diversity comparison using long-read sequencing data","fulltext":[{"header":"Main","content":"\u003cp\u003eThe emergence of SARS-CoV-2 in 2019, alongside historical pandemics, emphasizes RNA viruses, with their characteristically high evolutionary rates, as the most concerning group of pathogens. Their remarkable adaptability stands as a primary factor behind their widespread success\u003csup\u003e1\u003c/sup\u003e. Their rapid mutation stemming from a lack of proofreading during genome replication\u003csup\u003e2\u003c/sup\u003e, which serves as a mechanism fostering genetic variability. This variability manifests as a cloud of closely related genetic variants of the virus within the viral population infecting a single individual, which is known as viral quasispecies\u003csup\u003e3\u003c/sup\u003e. Studies elucidating quasispecies dynamics in relation to virus adaptation have significantly enhanced our understanding of RNA virus microevolution, virus-host interactions, and the ability of viruses to adapt to evade host immunity\u003csup\u003e4\u0026ndash;9\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eThe viral quasispecies concept has been referred to in nearly one thousand virology studies, notably for hepatitis C and B viruses, and HIV\u003csup\u003e10\u003c/sup\u003e, where some viral variants correlate with antiviral drug resistance\u003csup\u003e11\u0026ndash;14\u003c/sup\u003e. The advent of deep sequencing technology, alongside improvements in de novo viral genome assembly techniques, has enhanced our ability to detect minority variants within the viral cloud that may hold clinical significance\u003csup\u003e15\u0026ndash;18\u003c/sup\u003e. While next-generation sequencing (NGS) has been widely employed for viral quasispecies analysis\u003csup\u003e19,20\u003c/sup\u003e, one of the persisting challenges is the reconstruction of continuous viral quasispecies haplotypes\u0026mdash;sets of identical genomic sequences. This challenge arises due to the short length of NGS reads (typically 100\u0026ndash;400 bases), which often do not cover the entire targeted gene or genome, necessitating short-read assembly techniques that are often unable to discern between haplotypes that are closely related\u003csup\u003e17,21\u0026ndash;24\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eLong-read sequencing technology provides a solution to the coverage limitations of NGS technology in viral quasispecies research\u003csup\u003e21\u003c/sup\u003e, yielding read lengths typically exceeding 1,000 bp. Few analytical workflows or software exist that are tailored for viral quasispecies exploration using long-read data, and these primarily focus on single nucleotide variant (SNV) calling, haplotype reconstruction, or quasispecies profiling on an individual sample basis\u003csup\u003e25\u0026ndash;33\u003c/sup\u003e. However, these packages do not fully address the unique challenges or harness the full potential of long-read sequencing for viral quasispecies analysis. One challenge is the higher error rates (10\u0026ndash;30%) compared to NGS (~\u0026thinsp;1%)\u003csup\u003e34\u003c/sup\u003e, affecting haplotype reconstruction accuracy. Efforts in post-sequencing error correction have spurred the development of numerous tools for improving sequencing reads generated by Oxford Nanopore Technologies (ONT) or Pacific Biosciences (PacBio)\u003csup\u003e34\u0026ndash;38\u003c/sup\u003e. However, even with lower error rates, the rapid mutation rate of RNA viruses means that longer reads often differ by only one or two bases, and thus are called a unique haplotype. This ultimately results in a large number of haplotypes that each occur only once within the data, and obscures our ability to discern structure the viral quasispecies.\u003c/p\u003e \u003cp\u003eSequencing read depth and length also significantly influence viral quasispecies diversity measures\u003csup\u003e21,39,40\u003c/sup\u003e. Certain diversity metrics, such as Shannon entropy and mutation frequency, are dependent on sample size or read depth\u003csup\u003e40,41\u003c/sup\u003e. As yet, it is uncertain how such diversity metrics are influenced by sampling depth, particularly when comparing diversity between samples with different depths or down-sampled to the same depth, or even what minimum depth is needed to attain a robust measure of diversity. In addition, the longer the gene or genome length used for quantifying quasispecies, the greater the number of unique haplotypes or singletons due to the increased detection of mutations along its length. This variability may introduce bias, particularly when comparing viral quasispecies profiles across longitudinal or related samples.\u003c/p\u003e \u003cp\u003eIn response to these challenges, we developed \u003cem\u003elongreadvqs\u003c/em\u003e, an R package designed for viral quasispecies comparison using long-read data that can be applied to any RNA virus. This package creates a customizable workflow for long-read noise-minimization and down-sampling as well as for grouping related haplotypes into operational taxonomic units (OTUs). In addition, this package can analyze multiple samples together, identifying common haplotypes or OTUs that recur in different samples. This feature is particularly useful when analyzing longitudinal samples collected from a single individual or multiple individuals that may be epidemiologically linked. Our suggested analytical steps \u0026mdash; QoALa: Quasi-species Optimized \u0026amp; Adaptive Long-read Analysis \u0026mdash; performed by the key functions of the package, were tested with publicly available long-read data sequenced by both ONT and PacBio technologies from the most studied viruses in the field of viral quasispecies.\u003c/p\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e\n \u003ch2\u003eQoALa workflow and example datasets\u003c/h2\u003e\n \u003cp\u003eThe goals of the QoALa workflow are to standardize and comprehensively compare viral quasispecies and OTU profiles across multiple samples based on user customization. For each comparison, equal length read alignments from the same gene or genomic region of interest (Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ea) are separately imported from any genome assembly pipeline\u0026rsquo;s output (see examples in Methods) and noise-minimized using position-wise nucleotide base replacement in the \u0026ldquo;\u003cem\u003evqsassess\u003c/em\u003e\u0026rdquo; function (Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003eb). The noise-minimization step replaces potential erroneous nucleotides, defined as SNVs with a frequency lower than n% of the total read depth (cut-off percentage), with either the majority base or the dominant haplotype\u0026rsquo;s base at that position. This cut-off percentage can be determined using prior information such as the sequencing error rate, or by observing the change in the percentage of singleton haplotypes at different percentage cut-offs (see Methods and Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ea). Additionally, when comparing multiple samples, all alignments must be randomly down-sampled by the \u0026ldquo;\u003cem\u003evqsassess\u003c/em\u003e\u0026rdquo; function to achieve equivalent read depth (the shallowest depth among all samples is recommended) to reduce bias in diversity measurement caused by sample size disparity. SNV profiles of the samples visualized by the \u0026ldquo;\u003cem\u003esnvcompare\u003c/em\u003e\u0026rdquo; function should be inspected to evaluate the result of noise-minimization (Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ec). Finally, prepared samples are pooled by the \u0026ldquo;\u003cem\u003evqscompare\u003c/em\u003e\u0026rdquo; function to (\u003cspan class=\"CitationRef\"\u003e1\u003c/span\u003e) identify common haplotypes, (\u003cspan class=\"CitationRef\"\u003e2\u003c/span\u003e) reclassify haplotypes into OTUs, (\u003cspan class=\"CitationRef\"\u003e3\u003c/span\u003e) visualize diversity profiles and genetic relationships among samples (Fig. \u003cspan class=\"InternalRef\"\u003e1\u003c/span\u003ed), and (\u003cspan class=\"CitationRef\"\u003e4\u003c/span\u003e) summarize quantitative diversity metrics.\u003c/p\u003e\n \u003cp\u003eWe evaluated this workflow using read alignments assembled from publicly available raw reads of the five most studied viruses in the viral quasispecies field, namely HCV, HBV, HIV, SARS-CoV-2, and IAV. These datasets varied not only by sequencing technologies (PacBio or ONT) and sampling scenarios (short- or long-term longitudinal or cross-sectional), but also by the depth of coverage per gene of interest, ranging from 130 reads in the IAV\u0026apos;s M gene segment sample (Supplementary Table 6) to 854,147 reads in the SARS-CoV-2\u0026apos;s ORF3a gene sample (Supplementary Table 5). Besides, the level of sequencing error varied, estimated from the percentage of singleton haplotypes before noise-minimization. This value ranged from approximately 25% being singleton haplotypes in the SARS-CoV-2\u0026apos;s ORF3a gene sample to 100% in the HCV and HIV samples (Supplementary Fig. 1). This considerable variation highlights the necessity for customized long-read data preparation on a case-by-case basis using our package, particularly focusing on two key parameters: the cut-off percentage and sample size (Methods).\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec4\" class=\"Section2\"\u003e\n \u003ch2\u003eEffect of down-sampling on diversity metrics\u003c/h2\u003e\n \u003cp\u003eWhile read down-sampling is necessary for standardizing sample size before comparison, it may alter viral quasispecies diversity measures, especially at very low sample sizes. To assess this impact, we sub-sampled reads with replacement and computed nine diversity metrics, ranging from 10,000 down to 25 reads, depending on the original depth of each read alignment sample (Methods). This process was repeated 100 times per sub-sample size. This assessment was conducted twice: once for the alignments that were sampled after noise-minimization, and another for the alignments that were sampled before noise-minimization.\u003c/p\u003e\n \u003cp\u003eAs anticipated, diversity metrics directly reliant on read depth [number of haplotypes (H) and mutation frequency at the molecular level (Mfm)], indirectly reliant on read depth [Shannon entropy (H\u003csub\u003eS\u003c/sub\u003e), mutation frequency at the molecular level (Mfe), functional attribute diversity (FAD), and nucleotide diversity at the entity level (\u0026pi;\u003csub\u003ee\u003c/sub\u003e)], or normalized by read depth [normalized Shannon entropy (H\u003csub\u003eSN\u003c/sub\u003e)] showed varying degrees of fluctuation corresponding to sample size (Supplementary Figs.\u0026nbsp;2\u0026ndash;11). Down-sampling post noise-minimization from high-depth alignments (\u0026gt;\u0026thinsp;10,000 reads) minimally affected H\u003csub\u003eS\u003c/sub\u003e and Mfm at sample sizes over 1,000 reads (Supplementary Figs. 2, 4, and 8), while sampling from low-depth alignments (\u0026lt;\u0026thinsp;1,000 reads) resulted in abrupt changes in these metrics at any sample size (Supplementary Figs. 6 and 10). The Gini-Simpson index (H\u003csub\u003eGS\u003c/sub\u003e) and nucleotide diversity (\u0026pi;\u003csub\u003em\u003c/sub\u003e) remained relatively stable across all sample sizes, irrespective of the original alignment depth (Fig. \u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e). Computing metrics from alignments sampled after noise-minimization generally yielded less varied values compared to those sampled before noise-minimization in most instances (Supplementary Figs.\u0026nbsp;3, 5, 7, 9, and 11).\u003c/p\u003e\n \u003cp\u003eRegarding inter-sample comparisons, variations in metrics due to sample size and sampling strategy marginally affected comparisons in scenarios where the actual diversity measures of samples were either highly similar (Fig. \u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e and Supplementary Figs. 6\u0026ndash;11) or markedly different (Fig. \u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e and Supplementary Figs. 2\u0026ndash;5). However, down-sampling to very low sample sizes (\u0026lt;\u0026thinsp;100 reads) or down-sampling before noise-minimization likely resulted in false equivalences or differences because of the high variation in computed diversity metrics (Fig. \u003cspan class=\"InternalRef\"\u003e2\u003c/span\u003e and Supplementary Figs. 2\u0026ndash;11). Thus, our analysis suggests that diversity metrics from samples with read depths of \u0026lt;\u0026thinsp;100 should be interpreted with caution.\u003c/p\u003e\n \u003cp\u003eEach dot indicates the mean of one target, boxes the 25th to 75th percentile, lines medians and whiskers extend from minimum to maximum values.\u003c/p\u003e\n\u003c/div\u003e\n\u003cdiv id=\"Sec5\" class=\"Section2\"\u003e\n \u003ch2\u003eExample scenarios of diversity comparison\u003c/h2\u003e\n \u003cp\u003eSince a single mutation can differentiate one haplotype from another, and such mutations may introduce noise and obscure the true frequency of different major genetic variants within the quasispecies, solely using haplotype-based diversity metrics to summarize the viral quasispecies profile might hinder our ability to observe the continuous dynamics of the virus across samples. To illustrate this, we compared related viral samples in various scenarios using both traditional haplotyping based on a strict mutational profile and a novel OTU assignment based on genetic distance (Methods). In the latter, haplotypes are clustered by relatedness into larger operational taxonomic units (OTUs). This allows the user to better comprehend the structure of the quasispecies, i.e., the evolving genetic diversity surrounding several main variants within the viral cloud.\u003c/p\u003e\n \u003cp\u003eThe studies on HIV and SARS-CoV-2 served as examples of short-term longitudinal sampling scenarios. Within a few hours of each other, \u003cem\u003eenv\u003c/em\u003e genes of four HIV samples infecting T cells in vitro exhibited numerous haplotypes, among which a common dominant haplotype was not found (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003ea). When OTU clustering was applied to the same samples, the members and proportions of all four samples were relatively similar (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eb), with the 186 haplotypes being optimally clustered into 10 OTUs (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003ec/d). Examining five SARS-CoV-2\u0026rsquo;s \u003cem\u003eS\u003c/em\u003e gene samples taken within a nine-day period from a single patient, the dominant haplotype and OTU occupied over half of the virus population in each sample (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003ea/b), with minor haplotypes appearing as outliers (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003ec/d). It is also worth noting that the sporadic appearance of many minority variants in the day 13 and 17 SARS-CoV-2 samples (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003ea/b) was likely due to noise rising from incomparable read depths in the original alignments, which were over ten times lower than the other three samples (Supplementary Table\u0026nbsp;5).\u003c/p\u003e\n \u003cp\u003eAs a model for long-term longitudinal sample comparison, studies on HCV and HBV were utilized. Analysis of the NS4 coding region of HCV revealed that viral haplotypes and their proportions underwent sudden changes after treatment. Interestingly, certain major haplotypes detected pre-treatment reappeared 18 months post-treatment (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003ea). Concurrently, OTU clustering provided better insights into the dynamics of HCV across three sampling points by grouping genetically closely related minor haplotypes into common clusters observed at varying proportions across all time points (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eb/c/d). The within-host virus population dynamics were also clearly illustrated by HBV\u0026apos;s \u003cem\u003eS\u003c/em\u003e gene. Both haplotype and OTU classifications underscored the gradual turnover of predominant haplotypes from timepoint one to four over a 30-month period, and a completely different population emerged at timepoint five, occurring over 80 months after the initial sampling (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003ea/b).\u003c/p\u003e\n \u003cp\u003eLastly, in the scenario of the outbreak represented by IAV\u0026apos;s M gene segment samples, viruses were collected cross-sectionally among patients in two wards (B and C) within the same hospital. At the haplotype level, only two mostly identical sets of virus populations were found in more than one patient (B1, B2 and C1, C3), without other prominent links between cases or wards (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003ea). However, upon regrouping them into new OTUs, we discovered that all patients in ward B were partly infected with the same genetically similar sub-population, while all patients in ward C were infected with a similar set of OTUs. Furthermore, the predominant OTU found in two ward B patients (B1, B2) was also present as a minority group in all ward C patients, suggesting potential inter-ward transmission of IAV (Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003eb/c/d).\u003c/p\u003e\n \u003cp\u003eThe comparative profiles of viral quasispecies, whether analyzed at the haplotype or OTU levels, based on other genes from the same viral samples (including HBV\u0026rsquo;s \u003cem\u003eP\u003c/em\u003e gene, HIV\u0026rsquo;s \u003cem\u003egag\u003c/em\u003e gene, and SARS-CoV-2\u0026rsquo;s \u003cem\u003eORF3a\u003c/em\u003e gene), exhibited visual similarities to the chosen genes depicted in Fig. \u003cspan class=\"InternalRef\"\u003e3\u003c/span\u003e (Supplementary Fig. 12 to 14). This suggests the absence of linkage disequilibrium between the observed genes of these three viruses, indicating that they may adequately represent the population dynamics at the whole-genome level. Interestingly, one comparable result between the original and our studies is the number and frequency of haplotypes found in the day 11 and day 15 SARS-CoV-2 samples. We captured a similar quasispecies diversity in the \u003cem\u003eS\u003c/em\u003e gene samples and discovered two SNV positions in the \u003cem\u003eORF3a\u003c/em\u003e gene that created two rare haplotypes not reported in the original study (Supplementary Table. 7 and Ko et al., 2021\u003csup\u003e46\u003c/sup\u003e ).\u003c/p\u003e\n\u003c/div\u003e"},{"header":"Discussion","content":"\u003cp\u003eThe \u003cem\u003elongreadvqs\u003c/em\u003e package was developed as a tool for both quantitative and qualitative analysis of viral quasispecies diversity, addressing the challenges posed by varied error rates and read depths inherent in long-read sequencing technologies as well as noise introduced by rare SNVs in the alignment. During the testing phase, we carefully considered these variations across diverse sampling scenarios and multiple viral species, ensuring the tool\u0026rsquo;s general applicability. The strengths of this tool lie in its ability to customize parameter settings based on prior information and in providing comprehensive visualizations that elucidate the dynamics of virus microevolution across multiple samples.\u003c/p\u003e \u003cp\u003eThe analytical findings derived from the QoALa workflow presented in this study provide unique perspectives and insights not previously explored in the original research. For instance, we utilized raw read data to reconstruct within-host virus haplotypes and OTUs, emphasizing an overview of population structure changes over time or among outbreak cases. In contrast, the original studies focused on specific aspects, such as mutation spectra\u003csup\u003e42\u003c/sup\u003e, deletions\u003csup\u003e43\u003c/sup\u003e, gene splicing\u003csup\u003e44\u003c/sup\u003e, or phylogenetic analysis of consensus sequences\u003csup\u003e45\u003c/sup\u003e. To build read alignments, we employed a uniform genome assembly workflow for all raw read samples, which was not tailored or reproduced from the original protocols created for each dataset. Furthermore, in the comparison step, we fixed the number of sample sizes and k-means clusters for all scenarios. Hence, discrepancies in details between the original findings and ours may have arisen, which we acknowledge as a limitation on our part. Nevertheless, it's essential to note that the objective of our study does not revolve around delivering new findings from each dataset but rather demonstrating the usage and versatility of this package.\u003c/p\u003e \u003cp\u003ePrevious studies have established a robust foundation in measuring viral quasispecies diversity and developing analytical workflows for NGS data over the past decade\u003csup\u003e20,40,41,47\u003c/sup\u003e. However, much of this research cannot be directly applied to long-read sequencing data. Long-read quasispecies analyses are hindered by the significantly higher level of noise present in long-read data driven by the higher likelihood that longer reads will have at least one SNV across their length. For NGS short-read data, it is recommended to employ techniques such as rarefaction, resampling, and fringe trimming based on haplotype frequency to estimate appropriate sample sizes and minimize bias from unbalanced samples before comparing diversity metrics\u003csup\u003e41\u003c/sup\u003e. The presence of a misrepresented large proportion of singleton haplotypes (Supplementary Fig.\u0026nbsp;1), compounded by the high sequencing error rate of current long-read technologies, hinders us from following such approaches.\u003c/p\u003e \u003cp\u003eInstead of excluding low-frequency haplotypes as fringe trimming does, we opt to retain all reads and replace low-frequency SNVs with the majority base at their respective positions. This approach offers two primary benefits. Firstly, low-frequency haplotypes, comprising both true mutations and errors, are not eliminated from the alignment; rather, potential erroneous SNVs within the haplotype are smoothed to the mode. This advantage is particularly crucial when comparing shallow depth read alignments, where removing reads belonging to low-frequency haplotypes would severely reduce sample size. Secondly, the false haplotype diversity resulting from long-read errors are mitigated by consolidating previously low-frequency haplotypes into larger haplotype groups (OTUs) that share common SNVs, which are more likely to represent true mutations and provide more clarity on the structured diversity existing within quasispecies.\u003c/p\u003e \u003cp\u003eDespite these advantages, many factors need to be considered in selecting the cut-off percentage for position-wise noise-minimization. Factors such as gene mutation rate, sequencing error rate, and the trade-off between retaining versus unifying rare haplotypes are critical. Especially when dealing with a lengthy alignment where singleton haplotypes are more common, one may have to decide whether to prioritize retaining most mutation details, thus making it difficult to quantify overall population complexity, or to focus on the broader picture of viral cloud diversity while potentially discarding some details that may or may not be significant.\u003c/p\u003e \u003cp\u003eStandardizing unbalanced sample sizes by down-sampling can introduce bias in diversity comparisons, especially when integrating noise-minimization. Since sampling techniques used for NGS data\u003csup\u003e41\u003c/sup\u003e are not appropriate for long-read data, we redesigned sampling strategies by integrating the concept of rarefaction and repeat sampling, taking into account the order of sampling before or after noise-minimization to simulate possible options for analysis. Here, we demonstrate the sensitivity of diversity metrics across our example scenarios. Based on this analysis, down-sampling after noise-minimization is our recommended technique. However, repeated down-sampling to observe the distribution of metrics at different sample sizes, at least between the largest and smallest samples in the dataset, is encouraged to better understand the sensitivities of diversity metrics to sampling depth in a particular project. Ultimately, we suggest interpreting quantitative metrics comparisons alongside qualitative profiles that uncover shared viral variants between related samples, either at the haplotype or OTU levels, to better illuminate the dynamics of viral quasispecies, as made possible by the QoALa workflow implemented in the \u003cem\u003elongreadvqs\u003c/em\u003e package.\u003c/p\u003e"},{"header":"Methods","content":"\u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eSequencing reads data\u003c/h2\u003e \u003cp\u003eTo evaluate our analytic approaches, we gathered raw sequencing reads from five projects focusing on whole genome sequence (WGS) of five different viruses: hepatitis C virus (HCV), hepatitis B virus (HBV), human immunodeficiency virus (HIV), severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), and influenza A virus (IAV), which are among the most studied species in the field of viral quasispecies. These reads were generated using either PacBio Single Molecule, Real-Time (SMRT) sequencing or Oxford Nanopore Technologies (ONT) and were sourced from the National Center for Biotechnology Information\u0026rsquo;s Sequence Read Archive (NCBI SRA). Specifically, we selected projects that obtained either longitudinal viral samples or cross-sectional samples from a single outbreak, where the genetic relationships were suitable for comparison at the quasispecies level.\u003c/p\u003e \u003cp\u003ePacBio read samples were obtained from three projects. First, three HCV samples from a chronic hepatitis patient, were collected before and after treatment over a span of 18 months\u003csup\u003e42\u003c/sup\u003e. Second, five samples were collected at various time points within a period of 100 months from an untreated HBV-infected patient\u003csup\u003e43\u003c/sup\u003e. Third, five longitudinal SARS-CoV-2 samples were collected between 8 to 17 days after clinical onset from a single patient\u003csup\u003e46\u003c/sup\u003e. Data from the remaining two projects were generated using ONT, which included four samples from HIV-infected T cells collected up to 24 hours post-infection\u003csup\u003e44\u003c/sup\u003e, and seven samples from different patients infected by IAV during a nosocomial outbreak in the same hospital\u003csup\u003e45\u003c/sup\u003e. The NCBI accession numbers of a total of 24 raw read samples and five projects were listed in Supplementary Table\u0026nbsp;1.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec9\" class=\"Section2\"\u003e \u003ch2\u003eGenome assembly and read alignment generation\u003c/h2\u003e \u003cp\u003eWe conducted comparable genome assembly steps for both PacBio and ONT raw reads. Initially, we employed the default settings of HiFiAdapterFilt v3.0.1\u003csup\u003e48\u003c/sup\u003e and Porechop v0.2.4\u003csup\u003e49\u003c/sup\u003e to trim PacBio and ONT adapter sequences, respectively. Subsequently, poor-quality PacBio reads\u0026mdash;those shorter than 1,000 bp or constituting the worst 10% of reads based on the final score\u0026mdash;were filtered out using Filtlong v0.2.1\u003csup\u003e50\u003c/sup\u003e. Similarly, ONT reads with an average quality score\u0026thinsp;\u0026lt;\u0026thinsp;7 were removed using NanoFilt v2.8.0\u003csup\u003e51\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eWe then aligned the trimmed and filtered reads against the NCBI reference genomic sequence of each viral species (Supplementary Table\u0026nbsp;1) using minimap2 v2.26\u0026rsquo;s presets\u003csup\u003e52\u003c/sup\u003e tailored for each sequencing technology (\"map-hifi\" for PacBio and \"map-ont\" for ONT). Subsequently, to prepare the resulting read alignments for downstream analyses, SAM files from read mapping were converted to FASTA format and fragmented into genes or coding regions based on the genomic annotation of each reference sequence, utilizing SeqKit v2.7.0\u003csup\u003e53\u003c/sup\u003e. It was ensured that the reads retained in the alignment covered the full length of the particular gene or region (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea). The depth and length of the final alignments were documented in Supplementary Tables\u0026nbsp;2 to 6.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec10\" class=\"Section2\"\u003e \u003ch2\u003eSequencing read noise-minimization\u003c/h2\u003e \u003cp\u003eA main function of the \u003cem\u003elongreadvqs\u003c/em\u003e package, \u0026ldquo;\u003cem\u003evqsassess\u003c/em\u003e,\u0026rdquo; aims to minimize noise from sequencing errors, artifacts, or rare mutations by performing position-wise nucleotide base replacement. We hypothesize that any SNV in a read alignment with a base frequency at each position less than the specified cut-off percentage is likely an artifact that should be replaced with the majority base of that position (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eb). The cut-off percentage can be determined based on prior knowledge, such as the estimated sequencing error rate of the technology used or documented mutation or evolutionary rates of the studied virus. However, the selection of the cut-off for noise-minimization can be further guided by the \u0026ldquo;\u003cem\u003epctopt\u003c/em\u003e\u0026rdquo; function, which demonstrates the percentage of singleton haplotypes (haplotypes with only a single read member) in the alignment. The percentage of singleton haplotypes decreases as the cut-off percentage increases, as low-frequency SNVs are replaced, creating more groups of identical reads (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea). The normal range for the percentage of singleton haplotypes found in previous studies using more accurate technologies like Sanger or NGS sequencing can serve as a reference for this step. Alternatively, the plot showing changes in the percentage of singleton haplotypes can be used to identify the cut-off value where the percentage of singleton haplotypes ceases decreasing, indicating a plateau (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea).\u003c/p\u003e \u003cp\u003eAn alternative option for noise-minimization involves replacing potential erroneous low-frequency SNV bases with the base of the dominant haplotype. However, this method is not recommended, since either the base used for replacement could also be an error, or the \u0026ldquo;dominant\u0026rdquo; haplotype before noise-minimization may be relatively rare or non-existent (100% singleton haplotypes).\u003c/p\u003e \u003cp\u003eIn some cases, a single cut-off percentage cannot be generalized for the entire gene or region of interest because the amount of noise may vary throughout the sequence length. For example, low ONT base calling accuracy is often reported in homopolymer regions (continuous identical bases)\u003csup\u003e54,55\u003c/sup\u003e or soft clipping (unaligned regions) may persist at the ends of reads after mapping in both technologies\u003csup\u003e56,57\u003c/sup\u003e. The \u0026ldquo;\u003cem\u003esnvcompare\u003c/em\u003e\u0026rdquo; function, which visualizes SNV distribution across sequence length between different samples (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ec), can help identify noise-rich regions indicated by the overaccumulation of SNV positions. Subsequently, the \u0026ldquo;\u003cem\u003evqscustompct\u003c/em\u003e\u0026rdquo; function can be applied to readjust the cutoff percentage for noise-minimization at specific region(s) (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eb).\u003c/p\u003e \u003cp\u003eTo validate the noise-minimization workflow, we extracted read alignments from specific genes or regions of five exemplary viral species. These alignments were chosen based on their sequence length, with each exceeding or closely approximating 1,000 bases, and their depth of coverage surpassing 100 reads, while ensuring that the soft-clipped region, if present, comprised less than 25% of the total length. Utilizing the \"\u003cem\u003epctopt\u003c/em\u003e\" function, we determined the optimal cut-off percentage for minimizing errors, aiming to achieve nearly 0% singleton haplotypes, or identifying the point where the percentage of singleton haplotypes in the median sample ceased decreasing, indicating a plateau phase (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ea). The selected alignments (with their respective cut-off percentages) comprised the NS4 coding region (15%) of HCV, the \u003cem\u003eP\u003c/em\u003e (5%) and \u003cem\u003eS\u003c/em\u003e (5%) genes of HBV, the \u003cem\u003eenv\u003c/em\u003e (15%) and \u003cem\u003egag\u003c/em\u003e (22%) genes of HIV, the \u003cem\u003eS\u003c/em\u003e (1%) and \u003cem\u003eORF3a\u003c/em\u003e (1%) genes of SARS-CoV-2, and the M gene segment (10%) of IAV (Supplementary Fig.\u0026nbsp;1).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eDown-sampling and diversity metric sensitivity\u003c/h2\u003e \u003cp\u003eFor the demonstration of our package usage, we specified basic settings for the \"\u003cem\u003evqsassess\u003c/em\u003e\" function. After the noise-minimization step described above, we then down-sampled to either the size of the shallowest depth sample or to 1,000 reads (for samples with a depth over 10,000 reads), for every sample. These steps were taken to prepare the samples for the final quasispecies comparison within the same viral species or project.\u003c/p\u003e \u003cp\u003eHowever, it's important to note that read depth or sample size significantly influences some quantitative metrics used for quasispecies diversity measures\u003csup\u003e40,41\u003c/sup\u003e. Aggressive down-sampling may result in information loss and misinterpretation. To evaluate the impact of down-sampling on between-sample diversity comparison, we conducted a sensitivity analysis of nine diversity metrics as fully described by Gregori et al.\u003csup\u003e40\u003c/sup\u003e, including the number of haplotypes (H), the Shannon entropy (H\u003csub\u003eS\u003c/sub\u003e), the normalized Shannon entropy (H\u003csub\u003eSN\u003c/sub\u003e), the Gini-Simpson index (H\u003csub\u003eGS\u003c/sub\u003e), the functional attribute diversity (FAD), the mutation frequency at the entity level (Mfe), the nucleotide diversity at the entity level (π\u003csub\u003ee\u003c/sub\u003e), the mutation frequency at the molecular level (Mfm), and the nucleotide diversity (π\u003csub\u003em\u003c/sub\u003e).\u003c/p\u003e \u003cp\u003eWe utilized two to four samples of one gene of interest per species for the analysis. Using the QSutils v1.18.0\u003csup\u003e58\u003c/sup\u003e embedded in our package, all nine metrics were computed from both unsampled and down-sampled read alignments. For samples with a depth over 10,000 reads (HCV, HBV, and SARS-CoV-2), down-sample sizes started from 10,000 reads and were gradually halved until the final size reached 78 reads. For samples with low depth (HIV and IAV), sample sizes were set at 300, 150, 100, 50, and 25 reads. Random down-sampling with replacement was repeated 100 times for each sample size after noise-minimization. The same repeated sampling strategy was also implemented before noise-minimization. Distributions of each metric were visually compared between different sample sizes and approaches (down-sampling after versus before noise-minimization).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003eViral quasispecies comparison\u003c/h2\u003e \u003cp\u003eOnce the within-species samples were noise-minimized and standardized to equal depths, they were ready for comparison using the \"\u003cem\u003evqscompare\u003c/em\u003e\" function (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ed). This function aggregates the prepared read alignments of the listed samples, initially identifying shared haplotypes between them, and then visualizes the proportion of unique haplotypes for each sample as a color-coded bar plot. In addition, it reclassifies haplotypes into new operational taxonomic units (OTUs) based on the genetic distance matrix, utilizing the \"\u003cem\u003edist.dna\u003c/em\u003e\" function in ape v5.7.1\u003csup\u003e59\u003c/sup\u003e, from the SNV alignment extracted from the pooled read alignment. In detail, the distance matrix is transformed into dissimilarity coordinates, which are subsequently clustered into new OTUs using classical multidimensional scaling and k-means clustering, respectively, through the \"\u003cem\u003ecmdscale\u003c/em\u003e\" and \"\u003cem\u003ekmeans\u003c/em\u003e\" functions in stats v4.3.1\u003csup\u003e60\u003c/sup\u003e. The number of OTUs or clusters of genetically closely related haplotypes must be customized by specifying the number of clusters (k). The proportions of OTUs are also depicted in a color-coded bar plot. The clustering scheme of OTUs, along with major haplotypes within them, is illustrated with corresponding colors in multidimensional scale (MDS) plots (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003ed). All resulting plots are generated using ggplot2 v3.4.4\u003csup\u003e61\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eApart from the comparative plots, the \"\u003cem\u003evqscompare\u003c/em\u003e\" function offers several other valuable outputs for in-depth investigation. These include noise-minimized down-sampled read and SNV alignments for all samples, along with classified haplotype and OTU identifications. Additionally, the nine quasispecies diversity metrics for each sample, computed using QSutils v1.18.0\u003csup\u003e58\u003c/sup\u003e, are tabulated. Another set of metrics are calculated based on consensus reads and read frequencies of each OTU, rather than haplotypes, to simplify the diversity quantification in a larger scale.\u003c/p\u003e \u003cp\u003eTo exemplify such comparative features for our examples, noise-minimized down-sampled alignments of each selected gene and viral species, prepared by the \"\u003cem\u003evqsassess\u003c/em\u003e\" function, are listed as input for individual \"\u003cem\u003evqscompare\u003c/em\u003e\" analyses. The number of clusters (k) for OTU classification via k-means clustering was set to 10 for every run. The resulting summary plot was used to illustrate the dynamics of viral quasispecies diversity across different example scenarios, ranging from short-term (hours to days in HIV and SARS-CoV-2 datasets) to long-term (months to years in HCV and HBV datasets) longitudinal samples, as well as samples from the same outbreak cohort (IAV dataset).\u003c/p\u003e \u003c/div\u003e "},{"header":"Declarations","content":"\u003cp\u003e \u003ch2\u003eContributions\u003c/h2\u003e \u003cp\u003eN.P. and K.V. conceived, initiated, and designed the study. N.P. developed the \u003cem\u003elongreadvqs\u003c/em\u003e package, collected raw read data, implemented the QoALa workflow, conducted all analyses, and drafted the manuscript. K.V. and D.C.S. supervised the development of the package and workflow. All authors contributed to the final version of the paper.\u003c/p\u003e \u003c/p\u003e\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\u003cp\u003eN.P. and K.V. conceived, initiated, and designed the study. N.P. developed the longreadvqs package, collected raw read data, implemented the QoALa workflow, conducted all analyses, and drafted the manuscript. K.V. and D.C.S. supervised the development of the package and workflow. All authors contributed to the final version of the paper.\u003c/p\u003e\u003ch2\u003eAcknowledgement\u003c/h2\u003e\u003cp\u003eWe thank J. Gregori and M. Guerrero-Murillo from VHIR Vall d\u0026rsquo;Hebron Research Institute for their support in troubleshooting the dependency of our package. We thank TGS. Williams from Guy's and St Thomas' NHS Foundation Trust for providing additional information of the IAV study. Individually, we thank D. Makau and J. Baker for testing the package.\u003c/p\u003e\u003ch2\u003eData Availability\u003c/h2\u003e\u003cp\u003eAll raw read sequencing data used as example datasets for this study are publicly available at the National Center for Biotechnology Information\u0026rsquo;s Sequence Read Archive (NCBI SRA), as well as being part of previously published publications. Accession numbers and corresponding publications were tabulated in Supplementary Table 1.\u003c/p\u003e\u003c/div\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003eCode availability\u003c/h2\u003e \u003cp\u003eAll codes used in this study are available as an R package \u003cem\u003elongreadvqs\u003c/em\u003e and accessible via the CRAN repository at \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://cran.r-project.org/web/packages/longreadvqs/index.html\u003c/span\u003e\u003cspan address=\"https://cran.r-project.org/web/packages/longreadvqs/index.html\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e \u003c/div\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eCarrasco-Hernandez R, J\u0026aacute;come R, L\u0026oacute;pez Vidal Y. Ponce De Le\u0026oacute;n, S. Are RNA Viruses Candidate Agents for the Next Global Pandemic? A Review. ILAR J. 2017;58:343\u0026ndash;58.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSteinhauer DA, Domingo E, Holland JJ. Lack of evidence for proofreading mechanisms associated with an RNA virus polymerase. Gene. 1992;122:281\u0026ndash;8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDomingo E, Sheldon J, Perales C. Viral Quasispecies Evolution. Microbiol Mol Biol Rev. 2012;76:159\u0026ndash;216.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDomingo E. Quasispecies Structure and Persistence of RNA Viruses. Emerg Infect Dis. 1998;4:521\u0026ndash;7.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMandary M, Poh. Impact of RNA Virus Evolution on Quasispecies Formation and Virulence. IJMS. 2019;20:4657.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eVignuzzi M, Stone JK, Arnold JJ, Cameron CE, Andino R. Quasispecies diversity determines pathogenesis through cooperative interactions in a viral population. Nature. 2006;439:344\u0026ndash;8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWoo H-J, Reifman J. A quantitative quasispecies theory-based model of virus escape mutation under immune selection. \u003cem\u003eProc. Natl. Acad. Sci. U.S.A.\u003c/em\u003e 109, 12980\u0026ndash;12985 (2012).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDomingo E, Garc\u0026iacute;a-Crespo C, Perales C. Historical Perspective on the Discovery of the Quasispecies Concept. Annu Rev Virol. 2021;8:51\u0026ndash;72.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLauring AS. Within-Host Viral Diversity: A Window into Viral Evolution. Annu Rev Virol. 2020;7:63\u0026ndash;81.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePubMed. \u003cem\u003ePubMed\u003c/em\u003e \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://pubmed.ncbi.nlm.nih.gov/\u003c/span\u003e\u003cspan address=\"https://pubmed.ncbi.nlm.nih.gov/\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMetzner K. The significance of minority drug-resistant quasispecies. In: Geretti AM, editor. Antiretroviral Resistance in Clinical Practice. London: Mediscript; 2006.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMonaco DC, Zapata L, Hunter E, Salomon H, Dilernia DA. Resistance profile of HIV-1 quasispecies in patients under treatment failure using single molecule, real-time sequencing. AIDS. 2020;34:2201.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePerales C. Quasispecies dynamics and clinical significance of hepatitis C virus (HCV) antiviral resistance. Int J Antimicrob Agents. 2020;56:105562.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKai Y, et al. Baseline quasispecies selection and novel mutations contribute to emerging resistance-associated substitutions in hepatitis C virus after direct-acting antiviral treatment. Sci Rep. 2017;7:41660.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMargeridon-Thermet S, et al. Ultra‐Deep Pyrosequencing of Hepatitis B Virus Quasispecies from Nucleoside and Nucleotide Reverse‐Transcriptase Inhibitor (NRTI)\u0026ndash;Treated Patients and NRTI‐Naive Patients. J INFECT DIS. 2009;199:1275\u0026ndash;85.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRozera G, et al. Massively parallel pyrosequencing highlights minority variants in the HIV-1 env quasispecies deriving from lymphomonocyte sub-populations. Retrovirology. 2009;6:15.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBaaijens JA, Aabidine AZE, Rivals E, Sch\u0026ouml;nhuth A. De novo assembly of viral quasispecies using overlap graphs. Genome Res. 2017;27:835\u0026ndash;48.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFritz A, et al. Haploflow: strain-resolved de novo assembly of viral genomes. Genome Biol. 2021;22:212.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHouldcroft CJ, Beale MA, Breuer J. Clinical and biological insights from viral genome sequencing. Nat Rev Microbiol. 2017;15:183\u0026ndash;92.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLu I-N, Muller CP, He FQ. Applying next-generation sequencing to unravel the mutational landscape in viral quasispecies. Virus Res. 2020;283:197963.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePosada-Cespedes S, Seifert D, Beerenwinkel N. Recent advances in inferring viral diversity from high-throughput sequencing data. Virus Res. 2017;239:17\u0026ndash;32.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHuang A, Kantor R, DeLong A, Schreier L, Istrail S, QColors. An algorithm for conservative viral quasispecies reconstruction from short and non-contiguous next generation sequencing reads. In Silico Biol. 2012;11:193\u0026ndash;201.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHong LZ, et al. BAsE-Seq: a method for obtaining long viral haplotypes from short sequence reads. Genome Biol. 2014;15:517.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMardis E. R. DNA sequencing technologies: 2006\u0026ndash;2016. Nat Protoc. 2017;12:213\u0026ndash;8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDilernia DA, et al. Multiplexed highly-accurate DNA sequencing of closely-related HIV-1 variants using continuous long reads from single molecule, real-time sequencing. Nucleic Acids Res. 2015;43:e129\u0026ndash;129.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHuang DW. Towards Better Precision Medicine: PacBio Single-Molecule Long Reads Resolve the Interpretation of HIV Drug Resistant Mutation Profiles at Explicit Quasispecies (Haplotype) Level. J Data Min Genomics Proteom 07, (2016).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDudouet P, et al. SARS-CoV-2 quasi-species analysis from patients with persistent nasopharyngeal shedding. Sci Rep. 2022;12:18721.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eArtyomenko A, et al. Long Single-Molecule Reads Can Resolve the Complexity of the Influenza Virus Composed of Rare, Closely Related Mutant Variants. J Comput Biol. 2017;24:558\u0026ndash;70.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJiao X, et al. QuasiSeq: profiling viral quasispecies via self-tuning spectral clustering with PacBio long sequencing reads. Bioinformatics. 2022;38:3192\u0026ndash;9.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLink RW, et al. HIV-Quasipore: A Suite of HIV-1-Specific Nanopore Basecallers Designed to Enhance Viral Quasispecies Detection. Front Virol. 2022;2:858375.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLuo X, Kang X, Sch\u0026ouml;nhuth A. Strainline: full-length de novo viral haplotype reconstruction from noisy long reads. Genome Biol. 2022;23:29.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNg TT-L, et al. Long-Read Sequencing with Hierarchical Clustering for Antiretroviral Resistance Profiling of Mixed Human Immunodeficiency Virus Quasispecies. Clin Chem. 2023;69:1174\u0026ndash;85.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSu J, Li S, Zheng Z, Lam T-W, Luo R. ClusterV-Web: a user-friendly tool for profiling HIV quasispecies and generating drug resistance reports from nanopore long-read data. Bioinf Adv. 2024;4:vbae006.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMorisse P, Marchet C, Limasset A, Lecroq T, Lefebvre A. Scalable long read self-correction and assembly polishing with multiple sequence alignment. Sci Rep. 2021;11:761.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSahlin K, Medvedev P. Error correction enables use of Oxford Nanopore technology for reference-free transcriptome analysis. Nat Commun. 2021;12:2.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWang L, Qu L, Yang L, Wang Y, Zhu H. NanoReviser: An Error-Correction Tool for Nanopore Sequencing Based on a Deep Learning Algorithm. Front Genet. 2020. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003e10.3389/fgene.2020.00900\u003c/span\u003e\u003cspan address=\"10.3389/fgene.2020.00900\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eXiao C-L, et al. MECAT: fast mapping, error correction, and de novo assembly for single-molecule sequencing reads. Nat Methods. 2017;14:1072\u0026ndash;4.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSalmela L, Walve R, Rivals E, Ukkonen E. Accurate self-correction of errors in long reads using de Bruijn graphs. Bioinformatics. 2017;33:799\u0026ndash;806.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZagordi O, D\u0026auml;umer M, Beisel C, Beerenwinkel N. Read length versus Depth of Coverage for Viral Quasispecies Reconstruction. PLoS ONE. 2012;7:e47046.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGregori J, et al. Viral quasispecies complexity measures. Virology. 2016;493:227\u0026ndash;37.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGregori J, et al. Inference with viral quasispecies diversity indices: clonal and NGS approaches. Bioinformatics. 2014;30:1104\u0026ndash;11.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNakamura F, et al. Mutational spectrum of hepatitis C virus in patients with chronic hepatitis C determined by single molecule real-time sequencing. Sci Rep. 2022;12:7083.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eArasawa S, et al. Evolutional transition of HBV genome during the persistent infection determined by single-molecule real-time sequencing. Hepatol Commun. 2023;7:e0047\u0026ndash;0047.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNguyen Quang N, et al. Dynamic nanopore long-read sequencing analysis of HIV-1 splicing events during the early steps of infection. Retrovirology. 2020;17:25.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWilliams TGS et al. Feasibility and clinical utility of local rapid Nanopore influenza A virus whole genome sequencing for integrated outbreak management, genotypic resistance detection and timely surveillance. Microb Genomics 9, (2023).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKo SH, et al. High-throughput, single-copy sequencing reveals SARS-CoV-2 spike variants coincident with mounting humoral immunity during acute COVID-19. PLoS Pathog. 2021;17:e1009431.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKnyazev S, Hughes L, Skums P, Zelikovsky A. Epidemiological data analysis of viral quasispecies in the next-generation sequencing era. Brief Bioinform. 2021;22:96\u0026ndash;108.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSim SB, Corpuz RL, Simmonds TJ, Geib SM. HiFiAdapterFilt, a memory efficient read processing pipeline, prevents occurrence of adapter sequence in PacBio HiFi reads and their negative impacts on genome assembly. BMC Genomics. 2022;23:157.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWick RR, Judd LM, Gorrie CL, Holt KE. Completing bacterial genome assemblies with multiplex MinION sequencing. Microb Genomics 3, (2017).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eWick R. rrwick/Filtlong. (2024).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDe Coster W, D\u0026rsquo;Hert S, Schultz DT, Cruts M, Van Broeckhoven C. NanoPack: visualizing and processing long-read sequencing data. Bioinformatics. 2018;34:2666\u0026ndash;9.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLi H. New strategies to improve minimap2 alignment accuracy. Bioinformatics. 2021;37:4572\u0026ndash;4.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eShen W, Le S, Li Y, Hu F, SeqKit. A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation. PLoS ONE. 2016;11:e0163962.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSarkozy P, Jobb\u0026aacute;gy \u0026Aacute;, Antal P. Calling Homopolymer Stretches from Raw Nanopore Reads by Analyzing k-mer Dwell Times. in \u003cem\u003eEMBEC\u003c/em\u003e \u0026amp; \u003cem\u003eNBC 2017\u003c/em\u003e (eds. Eskola, H., V\u0026auml;is\u0026auml;nen, O., Viik, J. \u0026amp; Hyttinen, J.) vol. 65 241\u0026ndash;244Springer Singapore, Singapore, (2018).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHuang Y-T, Liu P-Y, Shih P-W. Homopolish: a method for the removal of systematic errors in nanopore sequencing by homologous polishing. Genome Biol. 2021;22:95.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDelahaye C, Nicolas J. Sequencing DNA with nanopores: Troubles and biases. PLoS ONE. 2021;16:e0257521.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eZhang S-J, et al. Isoform Evolution in Primates through Independent Combination of Alternative RNA Processing Events. Mol Biol Evol. 2017;34:2453\u0026ndash;68.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGuerrero-Murillo M. QSutils. [object Object] \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://doi.org/10.18129/B9.BIOC.QSUTILS\u003c/span\u003e\u003cspan address=\"10.18129/B9.BIOC.QSUTILS\" targettype=\"DOI\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e (2018).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eParadis E, Schliep K. Ape 5.0: An environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics. 2019;35:526\u0026ndash;8.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eR Core Team. R: A language and environment for statistical computing. R Foundation Stat Comput (2019).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGinestet C. ggplot2: Elegant Graphics for Data Analysis. J Royal Stat Society: Ser (Statistics Society) 174, (2011).\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":true,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"","lastPublishedDoi":"10.21203/rs.3.rs-4637890/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-4637890/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eThe concept of viral quasispecies refers to a constantly mutating viral population occurring within hosts, which is essential for grasping the micro-evolutionary patterns of viruses. Despite its high error rate, long-read sequencing holds potential for advancing viral quasispecies research by resolving coverage limitations in next-generation sequencing. We introduce a refined workflow, QoALa, implemented in the \u003cem\u003elongreadvqs\u003c/em\u003e R package. This workflow begins with nucleotide position-wise noise minimization of read alignments and sample size standardization, and extends to viral quasispecies comparison across related samples. Raw read samples from five studies of different viruses (HCV, HBV, HIV, SARS-CoV-2, and IAV), sequenced by major long-read platforms, were used to evaluate these approaches. The comparative results provide novel insights into intra- and inter-host diversity dynamics in various scenarios and unveil rare haplotypes not reported in the original study, underscoring the versatility and practicality of our methodology.\u003c/p\u003e","manuscriptTitle":"QoALa: a comprehensive workflow for viral quasispecies diversity comparison using long-read sequencing data","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-07-18 09:59:29","doi":"10.21203/rs.3.rs-4637890/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"927c433f-c6c7-4091-99e2-65020bb3248b","owner":[],"postedDate":"July 18th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[],"tags":[],"updatedAt":"2026-04-30T13:56:30+00:00","versionOfRecord":{"articleIdentity":"rs-4637890","link":"https://doi.org/10.1371/journal.pcbi.1014208","journal":{"identity":"plos-computational-biology","isVorOnly":true,"title":"PLOS Computational Biology"},"publishedOn":"2026-04-28 00:00:00","publishedOnDateReadable":"April 28th, 2026"},"versionCreatedAt":"2024-07-18 09:59:29","video":"","vorDoi":"10.1371/journal.pcbi.1014208","vorDoiUrl":"https://doi.org/10.1371/journal.pcbi.1014208","workflowStages":[]},"version":"v1","identity":"rs-4637890","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-4637890","identity":"rs-4637890","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}