The benchmarking pipeline ( Figure 1) started with an assessment of the quality of all long-read FASTQ inputs using NanoPlot ⁵¹ to ensure data conformity, particularly regarding the median read length. This study evaluates the performance and accuracy of seven assemblers on three different computational systems. Subsequently, the quality and accuracy of each assembler’s output are assessed using QUAST. Additionally, sequence similarity and HIV-1 subtype of all assembled contigs are investigated. The benchmarking pipeline is containerized, enabling its execution on any cloud or non-cloud environment ⁵² supporting either Docker or Singularity.

This study aimed to assess the accuracy, performance, and computational efficacy of seven candidate long-read assemblers, categorized into two approaches for haplotype reconstruction: (1) de novo long-read assemblers, which include Canu, ⁴⁷ Goldrush, ⁴⁸ MetaFlye, ⁴² and Strainline, ²⁸ and (2) reference-based long-read assemblers, which comprise HaploDMF, ⁴³ iGDA, ⁴⁴ and RVHaplo. ⁴⁵

The assembly software, including Canu, MetaFlye, and iGDA, was installed using the Micromamba package manager, which utilized recipes from the conda-forge ( https://anaconda.org/conda-forge) and bioconda channels ( https://anaconda.org/bioconda). An additional channel, zhixingfeng ( https://anaconda.org/zhixingfeng), was used specifically for iGDA as per the installation guidelines. Each tool was set up in separate environments to adhere to the developer’s specified version of dependencies.

For HaploDMF, RVHaplo, and Strainline, specifications were not available in Micromamba, necessitating manual installation. The source codes for HaploDMF ( https://github.com/dhcai21/HaploDMF) and RVHaplo ( https://github.com/dhcai21/RVHaplo) were downloaded from their respective GitHub repositories, with dependencies installed via Micromamba using environment specification files provided by the developers. For Strainline, the source code was sourced from its GitHub repository ( https://github.com/HaploKit/Strainline). Dependencies were managed and installed with Micromamba, except for Daccord ( https://github.com/gt1/daccord) and Metabat2 ( https://bitbucket.org/berkeleylab/metabat/src/master/), which were obtained by downloading the latest releases from the developers’ websites. Licensing for these tools varies: Strainline, GoldRush, HaploDMF, and RVHaplo are licensed under GPL3.0, Canu and iGDA under GPL2.0, and MetaFlye under the BSD-3-Clause license. Table 1 provides an overview of all candidate long-read assemblers.

Benchmarking analyses were conducted on three different computational systems: a server, a workstation PC, and a standard home PC. Specifications of each system is indicated in Table 2. The performance parameters of each assembler were measured as follows: 1) total CPU utilization, indicating the overall workload on the CPU, 2) memory usage, and 3) total runtime, defined as the elapsed wall-clock time from the start of assembly to the generation of output contigs. A dedicated Python script was developed to track the resource usage of each process of the long-read assemblers, capturing CPU and memory utilization every 100 milliseconds at each stage of the assembly process. As the pipeline operates within a containerized environment facilitated by Docker, unrelated processes such as operating system operations were filtered out during the measurement of resource utilization.

Long-read FASTQ files were simulated from four HIV-1 genome mixtures: (1) 2 group M HIV-1 subtypes (2M), (2) 2 CRF subtypes (2C), (3) 1 group M HIV-1 subtype and 1 CRF subtype (1M1C), and (4) 2 group M HIV-1 subtypes and 1 CRF subtype (2M1C), see Table 3. Each genome mixture comprised 100 sets of FASTA files containing corresponding HIV-1 genomes randomly selected from the Los Alamos HIV sequence databases ⁵³ ( https://www.hiv.lanl.gov/). All 400 HIV-1 FASTA sets underwent data simulation to generate 400 individual long-read FASTQ files using NanoSim V3.1.0. A list of mixtures within the 400 FASTA sets is available in Extended data, Table S6. ⁵⁴ Additionally, 20 long-read FASTQ files, comprising 5 samples from each mixture ( Table 3), were generated with varying amplicon sizes as described in the assembly quality assessment.

The data simulations were conducted using NanoSim V3.1.0 ( https://github.com/bcgsc/NanoSim/releases/tag/v3.1.0) in metagenomic mode, ^{55 , 56} employing a pretrained read profile model named human NA12878 DNA FAB49712 guppy. Additional simulation conditions included a maximum read length of 12,000 nt, a minimum read length of 7,500 nt, a median read length of 9,000 nt, and a standard deviation of read length set to 0.75. The parameter “--perfect” was configured as True to simulate error-free reads. A depth of coverage of 2,000 was selected, deemed sufficient for variant calling and quasispecies detection. ^{57 , 58} The resulting FASTQ files ( Table 3) were similar to those generated by an Oxford-ONT R9.4 chemistry using a Guppy 3.1.5 basecaller. Furthermore, this study benchmarked the assemblers using experimental data sourced from published studies (Extended data, Table S1). ⁵⁴ All FASTA templates for data simulation, ⁵⁹ simulated FASTQ files, ⁶⁰ and QC results of the simulated FASTQ files ⁶¹ are available as Underlying data.

Sequence similarity was assessed using local NCBI BLASTN ⁶² against a customized database of 12,000 HIV-1 genomes from the Los Alamos HIV-1 sequence database, ⁵³ with contigs showing >95% sequence similarity considered. QUAST ⁶³ was utilized to compare contigs from different assemblers, with regular QUAST for single viral isolate datasets and MetaQUAST ⁶³ for multiple isolate inputs. Primary QUAST parameters assessed were number of contigs, sizes, N50, % genome fraction, and total aligned bases, along with assembly correctness measured by average mismatches and indels per 100,000 aligned bases. Additionally, assembly quality was evaluated by average completeness of major HIV-1 ORFs (e.g., gag, gag-pol/pol, and env). The results generated from the pipeline are available as Underlying data. ⁶⁴

All seven assemblers were executed on three computational systems to measure wall-clock time usage, memory usage, and CPU utilization. Each assembler processed 100 simulated 2M FASTQ files. Figure 2 illustrates the time usage of all FASTQ files processed by different assembly pipelines across the three computational systems. As shown in Extended data, Table S2, ⁵⁴ the workstation server completed all 700 assemblies in 3 days, 9 hours, 50 minutes, and 28.31 seconds, while the workstation PC required 3 days, 23 hours, 24 minutes, and 2.64 seconds. Due to a failed Strainline-mediated assembly with the generic home PC, the remaining 600 assemblies took 2 days, 14 hours, 49 minutes, and 16.79 seconds. Notably, GoldRush emerged as the fastest assembler, with approximately 30 seconds or less per assembly under all three computational systems, while MetaFlye was the slowest, taking approximately 1,370 seconds or more per assembly.

For memory usage, Strainline could not be tested with the generic home PC due to a segmentation fault error from memory limitation. From the 189-Gb workstation server, the order of maximum memory usage, from highest to lowest, was Strainline, MetaFlye, GoldRush, HaploDMF, iGDA, Canu, and RVHaplo ( Figure 3A and B). Similarly, from the 64-Gb workstation PC, the order of maximum memory usage was Strainline, Canu, MetaFlye, HaploDMF, iGDA, RVHaplo, and GoldRush ( Figure 3C and D). From the 8-Gb generic home PC, the order of maximum memory usage was MetaFlye, iGDA, Canu, GoldRush, RVHaplo, and HaploDMF ( Figure 3E and F). Extended data, Table S3 ⁵⁴ presents the maximum memory usages of the six assemblers.

Figure 4 demonstrates the CPU utilizations of all assemblers, with the percentages of maximum CPU usages shown in Extended data, Table S3. ⁵⁴ On the workstation server ( Figure 4A), the assemblers with the highest to lowest averaged CPU usage were HaploDMF, Strainline, GoldRush, iGDA, RVHaplo, Canu, and MetaFlye. Both GoldRush (3,625.9%) and Canu (647.4%) did not overutilize the CPU and remained within the acceptable limit for 36 cores. On the workstation PC ( Figure 4B), the order of CPU usage was Strainline, HaploDMF, iGDA, Canu, RVHaplo, GoldRush, and MetaFlye. Notably, MetaFlye overutilized CPU at 1,437.4%, whereas Canu demonstrated only 637.5% peak utilization. Under the generic home PC system ( Figure 4C), the CPU usage order was GoldRush, iGDA, RVHaplo, MetaFlye, Canu, and HaploDMF. Both HaploDMF (726.3%) and Canu (628.9%) did not overutilize the CPU. Finally, Correlations between computational factors (e.g., CPUs, maximum memory, and assemblers) and runtime were analyzed using Jamovi statistical software with a non-parametric one-way ANOVA. ^{65 , 66} As demonstrated in Table 4, assembler selection significantly impacted runtime with a large effect size (ε ² = 0.903***), while no statistically significant differences in runtime were observed across different levels of CPUs and maximum memory.

The minimum read length was determined by testing the assemblers with FASTQ files from four different median read lengths: 1,000-nt, 2,000-nt, 3,000-nt, and 4,000-nt ( Table 3). This assessment aimed to ensure that the assembled contigs exhibited high contiguity, genome completeness, and recovered intact major HIV-1 open reading frames (ORFs) such as gag, gag-pol/pol, and env. The evaluation was conducted using the server system.

The first evaluation of the minimum read lengths was based on the assembled contig size. Increasing median read lengths resulted in longer contigs assembled by Canu ( Figure 5A). MetaFlye ( Figure 5B) and Strainline ( Figure 5C) could not process the 1,000-nt inputs, but Strainline assembled significantly longer contigs from longer reads. However, GoldRush failed the assessment. The remaining de novo assemblers, Canu, MetaFlye, and Strainline, processed the 4,000-nt median read length inputs and yielded contig sizes relatively close to a 9-kb HIV-1 genome. All reference-based assemblers performed well with all four inputs, generating contig sizes close to that of an HIV-1 genome ( Figure 5D-F). Complete statistical evaluation details can be found in Extended data, Table S9. ⁵⁴

Sequence similarity and genome fraction of assembled contigs were evaluated across different median read lengths. For the 1,000-nt reads (Extended data, Figure S1A), ⁵⁴ Canu yielded >5,000-nt contigs with >99.5% similarity, while iGDA and RVHaplo generated >8,000-nt contigs with 97% and 99% similarity, respectively. RVHaplo showed superior recovery of major HIV-1 ORFs (Extended data, Figure S2A). ⁵⁴ At the 2,000-nt reads (Extended data, Figure S1B), ⁵⁴ Canu and MetaFlye produced >7,000-nt contigs with >99% similarity, whereas Strainline’s contigs were ~3,000-nt with >98% similarity. All reference-based assemblers yielded longer contigs with greater similarity and improved ORF recovery (Extended data, Figure S2B). ⁵⁴ At the 3,000-nt reads (Extended data, Figure S1C), ⁵⁴ Strainline only produced 4,000-nt contigs with >99.50% similarity, while others generated longer contigs with higher similarity and improved ORF recovery (Extended data, Figure S2C). ⁵⁴ Finally, at the 4,000-nt reads (Extended data, Figure S1D), ⁵⁴ all assemblers produced longer contigs with higher similarity. All, except Canu and HaploDMF, demonstrated satisfactory % recovery of the major HIV-1 ORFs (Extended data, Figure S2D). ⁵⁴ RVHaplo yielded the highest averaged genome fraction (Extended data, Figure S1E). ⁵⁴

A further analysis ( Table 5) revealed a statistically significant positive correlation between the median read lengths of the inputs and the lengths of the assembled contigs (ρ=0.185***), as well as between the median read lengths and sequence similarity (ρ=0.163***). Additionally, assembler selection demonstrated a positive correlation with the size of assembled contigs (ρ=0.131***), but not with sequence similarity (ρ=-0.220). Thus, a key finding from this study suggested that a minimum read length of 2,000-nt is necessary for all assemblers, with optimal outcomes achieved at 4,000-nt. In summary, longer median read sizes are associated with higher quality contigs, characterized by longer contig length and greater sequence similarity.

This experiment evaluated the performance of assemblers in generating contigs from simulated FASTQ inputs containing heterogeneous HIV-1 sequences. These sequences were simulated from four HIV-1 genome mixtures: 2M (2 group M subtypes), 2C (2 circulating recombinant forms, CRFs), 1M1C (1 group M subtype and 1 CRF), and 2M1C (2 group M subtypes and 1 CRF), as shown in Table 3. The evaluation criteria included contiguity, genome completeness, the recovery of major HIV-1 open reading frames (ORFs), and the error rate. Initially, contig sizes were assessed, revealing that most fell within a range of 8,500–9,200 nt, close to the size of the HIV-1 genome (Extended data, Figure S3). ⁵⁴ Notably, MetaFlye produced shorter contigs, ranging from 7,800 to 8,800 nt. Detailed statistical analysis of this assessment is available in Extended data, Table S10. ⁵⁴

Considering both contig length and averaged genome fraction, all assembler pipelines processed the four datasets to yield contigs exceeding 8,000 nt, close to the size of the HIV-1 genome ( Figure 6A–D). However, the MetaFlye pipeline generated contigs with a relatively low (<90%) averaged genome fraction ( Figure 6E). Except for MetaFlye, all assembler pipelines demonstrated satisfactory recovery of major HIV-1 ORFs ( Figure 7A–D). De novo and reference-based assemblers showed a clear distinction in 2M, 2C, and 2M1C HIV-1 mixtures ( Figure 6A, B, and D), with de novo assemblers producing contigs with >95.5% sequence similarity. Reference-based assemblers HaploDMF and RVHaplo yielded contigs with relatively high sequence similarity (98–99.5%), while iGDA produced 95–97% sequence similarity. However, iGDA generated fewer contigs from the 2M mixture inputs ( Figure 6A) due to its reliance on overlapping SNVs for contig extension, which may fail with highly similar sequences. The presence of CRFs allowed iGDA to recover more contigs ( Figure 6B–D). Individual assembly evaluation statistics for each sample are available in Extended data, Table S10. ⁵⁴

Since all input data consisted of error-free reads, the assembly correctness ( Figure 8) was assessed by comparing the contigs with their respective HIV-1 genomes. The correctness of each assembler was evaluated based on the average number of mismatches per 100,000 aligned bases ( Figure 8A) and the average number of indels per 100,000 aligned bases ( Figure 8B). Among the four de novo assemblers, all introduced a few errors to the contigs, with Strainline exhibiting a greater degree of accuracy compared to the rest. MetaFlye, however, showed the least assembly correctness. Regarding the reference-based assemblers, HaploDMF and RVHaplo demonstrated similar degrees of assembly correctness, while iGDA exhibited the least. Interestingly, indels were predominantly associated with all reference-based assemblers in this study.

The final assessment investigated the HIV-1 subtype recall rates of the assemblers processing the 4 HIV-1 mixture datasets. In the 2M dataset ( Figure 9A), the reference-based assemblers, such as HaploDMF, RVHaplo, and iGDA, correctly recalled most of the subtypes, while the de novo assemblers exhibited lower recall rates. Strainline showed the best performance with >80% recall rates on all common subtypes, except for an overestimation on subtype B (118.52%). Canu, GoldRush, and MetaFlye exhibited low recall rates due to their collapsed assembly approach poorly handling sequences of relatively similar HIV-1 subtypes. Results from the 2C ( Figure 9B) and the 1M1C data ( Figure 9C) were similar to those observed in the 2M data. Interestingly, results from the 2M1C data ( Figure 9D) showed much lower recall rates (<60%) by all assemblers on all HIV-1 group M subtypes due to either a reduction in coverage or an error correction process. For the CRFs, the recall rates were similar to those shown in the 2C and 1M1C datasets.

Based on the results from the previous assessments, Strainline, MetaFlye, and HaploDMF were selected for further investigation using publicly available experimental sequencing data of HIV-1 and other viruses. Strainline, despite its high memory usage, offered rapid de novo viral haplotype reconstruction. MetaFlye, chosen for its compatibility with lower-end devices and acceptable memory usage and runtime, served as a de novo assembler. Despite its long runtime, HaploDMF was included as a representative reference-based assembler with memory efficiency comparable to MetaFlye.

This experiment assessed the performance of selected assemblers, including Strainline, MetaFlye, and HaploDMF, using experimental data of HIV-1 and other pathogenic viruses (Extended data, Table S1). ⁵⁴ The quality of assembled contigs was evaluated through QUAST, which provided metrics such as the number of contigs, sizes, N50, % genome fraction, and total aligned bases. The experiment was conducted using a workstation server.

HIV-1 recombinant forms and NL4-3

The dataset comprised plasma HIV-1 genomes (TRN01, TRN08, TRN09) and NL4-3 laboratory isolate from the MinION sequencer (BioProject: PRJDB13369). ⁶⁷ All assemblers generated NL4-3 contigs >9,000 nt ( Figure 10A). MetaFlye and HaploDMF gave contigs with 95.27% and 91.22% sequence similarity, respectively. Strainline produced nine contigs, the best at 94.34% similarity (Extended data, Table S4). ⁵⁴ For TRN01, all assemblers produced contigs identified by REGA (version 3.46) ⁶⁸ as subtype B and F recombinants, similar to a previous finding. ⁶⁷ Similarly, TRN08 yielded subtype B contigs. However, no recombinant forms with CRF01_AE were reconstructed. For TRN09, Strainline and HaploDMF replicated original findings, while HaploDMF also identified an A-C recombinant. In contrast, MetaFlye generated subtype A and CRF07_BC gagpol contigs, and a C-CRF01_AE recombinant. REGA subtyping results are in Extended data, Table S12. ⁵⁴ For the 10 5’ half genomes of NL4-3 (BioProject: PRJNA938445), ⁶⁹ all assemblers produced contigs meeting both length and alignment criteria. MetaFlye achieved the highest performance with an 89.26% average genome fraction, followed by Strainline at 85.93% and HaploDMF at 84.64%. Detailed alignment statistics are provided in Extended data, Table S13. ⁵⁴

Monkeypox (Mpox)

The first sequencing data, derived from MinION sequencing of a Mpox propagated in a CV-1 cell line (SRA: ERR10963128), ⁷⁰ contained 2,588 reads post-filtering, with an average length of 4,633 nt and a maximum length of 9,070 nt. HaploDMF generated 3 contigs, the longest being 188,872 nt, with 99.35% genome fraction. MetaFlye produced 7 contigs, the longest at 129,209 nt, with 95.98% genome fraction. Strainline yielded 12 contigs, the longest being 33,081 nt, with a genome fraction of 78.07% ( Figure 10B and Extended data, Table S4 ⁵⁴ ). The second dataset, derived from MinION sequencing of a Mpox from a pustular lesion (SRA: SRR15830920), ⁷¹ comprised 30,990 reads after filtering, with an average length of 9,048 nt and a maximum length of 59,093 nt. MetaFlye and HaploDMF exhibited similar performance, with MetaFlye assembling 122 contigs with an N50 of 13,335 nt and the longest having a 98.58% genome fraction. HaploDMF produced 5 contigs with an N50 of 173,134 nt and the longest contig having a genome fraction of 98.62%. In contrast, Strainline generated 7 contigs with an N50 of 30,255 nt and the longest contig of 101,047 nt. Only the 65,261-nt contig of those 7 contigs demonstrated 51.31% genome fraction ( Figure 10C and Extended data, Table S4 ⁵⁴ ).

Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2)

The dataset comprises GridION sequencing data of a SARS-CoV-2 clinical sample retrieved from the SRA database (SRP250446). ⁷² Strainline generated 3 contigs: 1 full-length and 2 duplicated half-genome contigs, with a 99.93% genome fraction. MetaFlye and HaploDMF each produced a single contig, with genome fractions of 96.77% and 93.11%, respectively. However, MetaFlye’s contig lacked the first 500 nucleotides of the genome, while HaploDMF’s contig had a 2,000-nt gap when aligned to the reference genome ( Figure 10D and Extended data, Table S4 ⁵⁴ ).

Polio virus mixture

This simulated ONT dataset consisted of 6 poliovirus type 2 sequences. ²⁸ In Figure 11A and Extended data, Table S5, ⁵⁴ Strainline and HaploDMF generated contigs with N50 values of 7,453-nt and 7,436-nt, respectively, closely matching the 7,400-nt genome of Polio virus. Strainline produced 6 contigs with an average genome fraction of 99.34%, all sharing >95% sequence similarity with 6 viral haplotypes. HaploDMF yielded 3 contigs with >95% sequence similarity and an average genome fraction of 99.80%. MetaFlye, however, produced 3 shorter contigs with an N50 of 4,830-nt. These contigs shared >95% sequence similarity with 5 haplotypes, but with an average genome fraction of 70.13%.

Hepatitis C virus (HCV) mixture

This simulated ONT dataset consisted of 10 strains of hepatitis C virus (HCV), Subtype 1a. ²⁸ In Figure 11B and Extended data, Table S5, ⁵⁴ Strainline and HaploDMF yielded contigs with N50 values of 9,286-nt and 9,309-nt, respectively. Strainline produced 12 contigs with an average genome fraction of 99.85%, all of them sharing >95% sequence similarity with all 10 haplotypes with 2 duplicated contigs. HaploDMF generated 3 contigs with a 99.99% average genome fraction, all exhibiting >95% sequence similarity with HCV haplotypes. MetaFlye produced 5 shorter contigs with an N50 of 5,670-nt, and the longest contig (8,006-nt) shared >95% sequence similarity with 5 HCV haplotypes. However, MetaFlye contigs had an average genome fraction of 77.40%.

ZIKA Virus (ZIKV) mixture

This simulated ONT dataset comprised 15 strains of Zika virus (ZIKV). ²⁸ In Figure 11C and Extended data, Table S5, ⁵⁴ Strainline, MetaFlye, and HaploDMF generated contigs with N50 values of 10,264-nt, 5,842-nt, and 10,264-nt, respectively, closely matching the size of the ZIKV genome (11,000-nt). Strainline produced 13 contigs that recovered 14 ZIKV haplotypes, with an average genome fraction of 99.84%. MetaFlye yielded 15 contigs with a lower average genome fraction of 78.56%, recovering 10 ZIKV haplotypes. Interestingly, HaploDMF produced 10 contigs, each with an average genome fraction of 99.90%, and recovered 11 ZIKV haplotypes.

In summary, the selected assemblers showed strong performance across different experimental ONT datasets. MetaFlye exhibited superior performance with experimental ONT data of full-length or half-length HIV-1 genomes. Strainline demonstrated excellent haplotype recovery with ONT sequencing data of various viral mixtures, including SARS-CoV-2, Poliovirus, HCV, and ZIKV. Additionally, MetaFlye showcased better performance on the metagenomic experimental ONT data of the 200-kb Mpox.

The HIV-64148 benchmarking pipeline is designed for portability, ensuring it functions across different computational environments with minimal setup. As a result, the containerized HIV-64148 can be executed on either Docker or Singularity platforms. The pipeline offers the flexibility to select from a range of assemblers—Canu, GoldRush, MetaFlye, Strainline, HaploDMF, iGDA, and RVHaplo—to best match the computational system or meet specific study requirements. It generates a comprehensive report in HTML format, detailing contigs/haplotype identities, and abundances. Additionally, the report includes protein sequence alignments of drug-targeted HIV-1 genes, alongside profiles of drug resistance and susceptibility associated with these genes. This is achieved by querying the assembled sequence against the Stanford University HIV Drug Resistance Database ^{73 – 78} via the provided Sierra web service 2 API, as illustrated in Figure 12. The downstream analysis pipeline could also serve as a model for implementation with other viruses. The source code for building the containerized pipeline is available in a public code repository: https://github.com/STTLab/HIV-64148.

Not applicable.