Benchmarking of long-read assemblers for prokaryote whole genome sequencing

Simulated read sets

Simulated read sets (read sequences generated in silico from reference genomes) offer some advantages over real read sets when assessing assemblers. They allow for a confident ground truth – i.e. the true underlying genome is known with certainty. They allow for large sample sizes, in practice limited only by computational resources. Also, a variety of genomes and read set parameters can be used to examine assembler performance over a wide range of scenarios. For this study, we simulated 500 read sets to test the assemblers, each using different parameters and a different prokaryote genome.

To select reference genomes for the simulated read sets, we first downloaded all bacterial and archaeal RefSeq genomes using ncbi-genome-download v0.2.10 (14333 genomes at the time of download)⁸. We then performed some quality control steps: excluding genomes with a >10 Mbp chromosome, a <500 kbp chromosome, any >300 kbp plasmid, any plasmid >25% of the chromosome size or more than 9 plasmids (Extended data, Figure S1)⁹. We then ran Assembly Dereplicator v0.1.0 with a threshold of 0.1, resulting in 3153 unique genomes¹⁰.

To produce a final set of 500 genomes with 500 plasmids, we randomly selected 250 genomes from those containing plasmids, repeating this selection until the genomes contained exactly 500 plasmids. We then added 250 genomes randomly selected from those without plasmids. Any ambiguous bases in the assemblies were replaced with ‘A’ to ensure that sequences contained only the four canonical DNA bases.

We then used Badread v0.1.5 to generate one read set for each input genome¹¹. The parameters for each set (controlling read depth, length, identity and errors) were randomly chosen to ensure a large amount of variability (Extended data, Figure S2)⁹. Note that not all of these read sets were sufficient to reconstruct the original genome (due to low depth or short read length), so even an ideal assembler would be incapable of completing an assembly for all 500 test sets.

For genomes containing plasmids, the read depth of plasmids relative to the chromosome was also set randomly, with limits based on the plasmid size (Extended data, Figure S3)⁹. Large plasmids were simulated at depths close to that of the chromosome while small plasmids spanned a wider range of depth. This was done to model the observed pattern that small plasmids often have a high per-cell copy number (i.e. may be high read depth) but can be biased against in library preparations (i.e. may be low read depth)¹². All replicons (chromosomes and plasmids) were treated as circular sequences in Badread, so the simulated read sets do not test assembler performance on linear sequences.

Real read sets

Despite the advantages of simulated read sets, they can be unrealistic because read simulation tools (such as Badread) may not accurately model all relevant features: error profiles, read lengths, quality scores, etc. Real read sets are therefore also valuable when assessing assemblers. The challenge with real read sets is obtaining a ground truth genome against which assemblies can be checked. Since many reference genome sequences are produced using long-read assemblies, there is the risk of circular reasoning – if we use an assembly as our ground truth reference, our results will be biased in favour of whichever assembler produced the reference.

To avoid this issue, we used the datasets produced in a recent study comparing ONT (MinION R9.4) and PacBio (RSII CLR) data which also included Illumina reads for each isolate¹³. For each of the 20 bacterial isolates in that study, we conducted two hybrid assemblies using Unicycler v0.4.7: Illumina+ONT and Illumina+PacBio¹⁴. Unicycler works by first generating an assembly graph using the Illumina reads, then using long-read alignments to scaffold the graph’s contigs into a completed genome – a distinct approach from any of the long-read assemblers tested in this study. We ran the assemblies using Unicycler’s --no_miniasm option so it skipped its Miniasm-based step which could bias the results in favour of Miniasm/Minipolish. We then excluded any isolate where either hybrid assembly failed to reach completion or where there were >50 nucleotide differences between the two assemblies as determined by a Minimap2 alignment¹⁵. I.e. the Illumina+ONT and Illumina+PacBio hybrid assemblies needed to be in near-perfect agreement with each other. This left six isolates for inclusion. The above process may have biased these isolates in favour of easier-to-assemble genomes, as more complex genomes would be more likely to encounter inconsistencies between the two Unicycler assemblies.

The ONT and PacBio read sets for these isolates were quite deep (156× to 535×) so to increase the number of assembly tests, we produced ten random read subsets of each, ranging from 40× to 100× read depth. This resulted in 120 total read sets for testing the assemblers (6 genomes × 2 platforms × 10 read subsets). The Illumina+ONT hybrid assembly was used as ground truth for each isolate.

All real and simulated read sets¹⁶ and reference genomes¹⁷ are available as Underlying data.

Assemblers tested

We assembled each of the read sets using the current versions of eight long-read assemblers: Canu v2.1, Flye v2.8, Miniasm/Minipolish v0.3/v0.1.3, NECAT v20200803, NextDenovo/NextPolish v2.3.1/v1.3.1, Raven v1.3.0, Redbean v2.5 and Shasta v0.7.0. Default parameters were used except where stated, and exact commands for each tool are given in the Extended data, Figure S4⁹. Assemblers that only work on PacBio reads (i.e. not on ONT reads) were excluded (HGAP¹⁸, FALCON¹⁹, HINGE²⁰ and Dazzler²¹), as were hybrid assemblers which also require short read input (Unicycler¹⁴ and MaSuRCA²²).

Canu has the longest history of all the assemblers tested, with its first release dating back to 2015. It performs assembly by first correcting reads, then trimming reads (removing adapters and breaking chimeras) and finally assembling reads into contigs²³. Its assembly strategy uses a modified version of the string graph algorithm²⁴, sometimes referred to as the overlap-layout-consensus (OLC) approach.

Flye takes a different approach to assembly: first combining reads into error-prone disjointigs, then collapsing repetitive sequences to make a repeat graph and finally resolving the graph’s repeats to make the final contigs²⁵. Of particular note to prokaryote assemblies, Flye has options for recovery of small plasmids (--plasmids) and uneven depth of coverage (--meta), both of which we used in this analysis.

Miniasm builds a string graph from a set of read overlaps – i.e. it performs only the layout step of OLC. It does not perform read overlapping which must be done separately with Minimap2, and it does not have a consensus step, so its assembly error rates are comparable to raw read error rates. A separate polishing tool such as Racon is therefore required to achieve high sequence identity²⁶. For this study, we developed a tool called Minipolish to simplify this process by conducting Racon polishing (two rounds by default) on a Miniasm assembly graph²⁷. To ensure clean circularisation of prokaryote replicons, circular contigs are ‘rotated’ (have their starting position adjusted) between polishing rounds. Minipolish also comes with a script (miniasm_and_minipolish.sh) which carries out all assembly steps (Minimap2 overlapping, Miniasm assembly and Minipolish consensus) in a single command, and subsequent references to ‘Miniasm/Minipolish’ refer to this entire pipeline.

NECAT follows an approach similar to Canu: first correcting the input reads, then building an assembly from the corrected reads²⁸. Both the correction and assembly steps are progressive, using multiple processing steps to achieve better accuracy/completeness.

NextDenovo is a performance-oriented assembler, which like Canu and NECAT performs read-correction at the start of its pipeline. It performs the first two steps of OLC (overlap and layout), leaving the final step (consensus) to a separate tool: NextPolish²⁹. We used both tools in conjunction in this study, referred to as ‘NextDenovo/NextPolish’.

Raven (previously known as Ra) is another tool which takes an OLC approach to assembly³⁰. Its overlapping step shares algorithms with Minimap2, and its consensus step is based on Racon, making it similar to Miniasm/Minipolish. It differs in its layout step which includes novel approaches to remove spurious overlaps from the graph, helping to improve assembly contiguity.

Redbean (previously known as Wtdbg2) uses an approach to long-read assembly called a fuzzy Bruijn graph³¹. This is modelled on the De Bruijn graph concept widely used for short-read assembly³² but modified to work with the inexact sequence matches present in noisy long reads.

Shasta is an assembler designed for computational efficiency³³. To achieve this, much of its assembly pipeline is performed not directly on read sequences but rather on a reduced representation of marker k-mers. These markers are used to find overlaps and build an assembly graph from which a consensus sequence is derived.

Computational environment

All assemblies were run on Ubuntu 18.04 instances of Australia’s Nectar Research Cloud which contained 32 vCPUs and 128 GB of RAM (r3.xxlarge flavour). To guard against performance variation caused by vCPU overcommit, the assemblers were limited to 16 threads (half the number of available vCPUs) in their options. Any assembly which exceeded 24 hours of runtime or 128 GB of memory usage was terminated.

Assembly assessment

Our primary metric of assembly quality was contiguity, defined here as the longest single Minimap2 alignment between the assembly and the reference replicon, relative to the reference replicon length. This provides a simpler picture of assembly quality than is created by QUAST (which quantifies misassemblies and other metrics such as NG50) but is appropriate for cases where complete assembly is likely². Contiguity of exactly 100% indicates that the replicon was assembled completely with no missing or extra sequence (Extended data, Figure S5A)⁹. Contiguity of slightly less than 100% (e.g. 99.9%) indicates that the assembly was complete, but some bases were lost at the start/end of the contig (Extended data, Figure S5B)⁹. Contiguity of more than 100% (e.g. 101%) indicates that the contig contains duplicated sequence via start-end overlap (Extended data, Figure S5C)⁹. Much lower contiguity (e.g. 70%) indicates that the assembly was not complete due to fragmentation (Extended data, Figure S5D)⁹, missing sequence (Extended data, Figure S5E)⁹ or misassembly (Extended data, Figure S5F)⁹. Contiguity values were determined by aligning the contigs to a tripled version of the reference replicon, necessary to ensure that contigs can fully align even with start-end overlap and regardless of their starting position relative to that of the linearised reference sequence (Extended data, Figure S6)⁹. To encourage longer alignments, Minimap2 was run with the asm20 preset, chain elongation threshold of 10 kbp, banding threshold of 10 kbp, Z-drop score of 1000 and inversion Z-drop score of 500. The script for conducting this analysis (assess_assembly.py) is available in Extended data.

Contiguity values were determined for each replicon in the assemblies – e.g. if a genome contained two plasmids, then the assemblies of that genome have three contiguity values: one for the chromosome and one for each plasmid. A status of ‘fully complete’ was assigned to assemblies where all replicons (the chromosome and any plasmids if present) achieved a contiguity of ≥99%. If an assembly had a chromosome with a contiguity of ≥99% but incomplete plasmids, it was given a status of ‘complete chromosome’. If the chromosome had a contiguity of <99%, the assembly was deemed ‘incomplete’. If the assembly was empty or missing (possibly due to the assembler prematurely terminating with an error), it was given a status of ‘empty’. Computational metrics were also observed for each assembly: time to complete and maximum RAM usage.

Reliability

Reliability was assessed using each assembler’s performance on the real read sets (Figure 2A). When considering only the chromosome, NextDenovo/NextPolish was the most reliable assembler, followed by Raven, NECAT, Canu and Flye – all were able to complete the chromosome in over three-quarters of their assemblies. If plasmids are also considered, then Canu was the most reliable assembler followed by Flye. Miniasm/Minipolish and Shasta were moderately reliable, completing over half of the chromosomes. Redbean was the least reliable and completed less than half of the chromosomes.

Robustness

Robustness was assessed using each assembler’s performance on the simulated read sets (Figure 1A) which contained a large amount of variation on many metrics (Extended data, Figure S7). NextDenovo/NextPolish and Raven were the most robust assemblers, able to complete the chromosome in over three-quarters of their assemblies. Flye, Redbean and Canu performed best in cases of low read depth, able to complete assemblies down to ~10× depth (Extended data, Figure S7A)⁹. Raven, NextDenovo/NextPolish and NECAT performed best with low-identity read sets (Extended data, Figure S7B)⁹. The assemblers performed similarly with regards to read length, except for Shasta which required longer reads (Extended data, Figure S7C)⁹. The assemblers were similarly unaffected by random reads, junk reads, chimeric reads or adapter sequences (Extended data, Figure S7D–F)⁹. Read glitches (local breaks in continuity) were more likely to cause assembly problems for Canu, NECAT and Shasta (Extended data, Figure S7G)⁹.

Identity

In our real read tests, Flye achieved the highest overall assembled sequence identity (Figure 2D). Canu achieved high sequence identity on PacBio reads. Miniasm/Minipolish, NextDenovo/NextPolish and Raven did well on ONT reads. For each assembler, real PacBio reads resulted in higher identities than real ONT reads. For the simulated reads (which contain artificial error profiles), results were more erratic, with Canu and Flye performing best (Figure 1D).

Regarding the maximum indel error size in the assemblies, Flye performed best, usually producing assemblies with errors no larger than 10 bp (Figure 1E/Figure 2E). NECAT and Shasta performed poorly, usually producing errors larger than 10 bp. The other assemblers had a large variance in this metric, sometimes producing assemblies with small errors and sometimes with large errors.

The nature of read errors depends on the sequencing platform and basecalling software used, so these results may not hold true for all read sets. Platform-specific post-assembly polishing tools (including Nanopolish⁷, Medaka³⁶ and Arrow³⁷) are routinely used to improve the accuracy of long-read assemblies³⁸, and these can often achieve assembly identities of >99.9% for ONT read sets and >99.999% for PacBio read sets (i.e. better than any of the assemblers were able to achieve on their own). Identity can be further increased by polishing with Illumina reads where available (e.g. with Pilon³⁹). Therefore, the sequence identity produced by the assembler itself is potentially unimportant for many users. However, large-scale indel errors may be less easily fixed using polishing tools and therefore could be of greater relevance.

Resource usage

Canu was the slowest assembler tested on both real (Figure 2F) and simulated (Figure 1F) read sets, sometimes taking hours to complete. Its runtime was correlated with read accuracy and read set size, with low-accuracy and large read sets being more likely to result in a long runtime.

Flye was typically faster than Canu, taking less than 15 minutes for the real read sets and usually less than an hour for the simulated read sets. It sometimes took multiple hours to assemble simulated read sets, and this was correlated with the amount of junk (low-complexity) reads, suggesting that removal of such reads via pre-assembly QC may be beneficial. Flye had the highest RAM usage of the tested assemblers and its RAM usage was correlated with read N50 and read set size, with long and large read sets being more likely to result in high RAM usage.

Shasta, Redbean and Raven were the fastest assemblers, typically completing assemblies in less than 5 minutes. While not tested in this study, Racon (which is used in Minipolish) and Raven can be run with GPU acceleration to further improve speed performance. Raven and Shasta had the lowest memory usage, typically requiring less than 4 GB of RAM.

Circularisation

Of all assemblers tested, Miniasm/Minipolish and NECAT most regularly achieved exact circularisation (contiguity=100%) (Figure 1B/Figure 2B). Flye often excluded a small amount of sequence (tens of bases) from the start/end of circular contigs (contiguity <100%), and Raven typically excluded moderate amounts of sequence (hundreds of bases). Contiguities for Canu and NextDenovo/NextPolish usually exceeded 100%, indicating a large amount (thousands of bases) of start/end overlap. The amount of overlap in a Canu or NextDenovo/NextPolish assembly was correlated with the read N50 length (Extended data, Figure S7C)⁹. Redbean and Shasta were both erratic in their circularisation, often producing some sequence duplication (contiguity >100%) but occasionally dropping sequence (contiguity <100%).

In addition to cleanly circularising contig sequences, it is valuable for a prokaryote genome assembler to clearly distinguish between circular and linear contigs. This can provide users with a clue as to whether or not the genome was assembled to completion. Flye, Miniasm/Minipolish, Raven and Shasta produce graph files of their final assembly which can indicate circularity. Canu indicates circularity via the ‘suggestCircular’ text in its contig headers. NECAT, NextDenovo/NextPolish and Redbean do not signal to users whether a contig is circular.

Plasmids

Canu and Flye were the two assemblers most able to assemble plasmids at a broad range of size and depth (Extended data, Figures S8, S9)⁹. Miniasm/Minipolish also performed well, though it failed to assemble plasmids if they were very small or had a very high read depth. Raven was able to assemble most large plasmids but not small plasmids. NECAT, NextDenovo/NextPolish, Redbean and Shasta were least successful at plasmid assembly.

Circularisation of plasmids followed the same pattern as for chromosomes, with Miniasm/Minipolish, Flye and NECAT most consistently achieving clean circularisation (Figure 1C/Figure 2C)⁹. For smaller plasmids, start/end overlap could sometimes result in contiguities of ∼200% – i.e. the plasmid sequence was duplicated in a single contig. This was most common with Canu and NextDenovo/NextPolish, though it occurred with other assemblers as well.

Ease of use

Most assemblers tested were relatively easy to use, either running with a single command (Canu, Flye, Raven and Shasta) or providing a convenience script to bundle the commands together (Miniasm/Minipolish and Redbean). NECAT requires a configuration file be prepared, making it somewhat cumbersome to run. NextDenovo/NextPolish was the most difficult to run, requiring multiple commands and multiple configuration files. All were able to take long reads in FASTQ format as input (Extended data, Figure S4)⁹. We encountered no difficulty installing any of the tools by following the instructions provided.

Some of the assemblers needed a predicted genome size as input (Canu, NECAT, NextDenovo/NextPolish and Redbean) while others (Flye, Miniasm/Minipolish, Raven and Shasta) did not. This requirement could be a nuisance when assembling unknown isolates, as it may be hard to specify a genome size before the species is known.

Configurability

While we ran our assemblies using default and/or recommended commands (Extended data, Figure S4)⁹, some of the assemblers have parameters which can be used to alter their behaviour. Raven was the least configurable assembler tested, with few options available to users. Flye offers some parameters, including overlap and coverage thresholds. Miniasm/Minipolish, NECAT, NextDenovo/NextPolish, Redbean and Shasta all offer more options, and Canu is the most configurable with hundreds of adjustable parameters. Many of the available parameters are arcane (e.g. Miniasm’s ‘max and min overlap drop ratio’ or Shasta’s ‘pruneIterationCount’), and only experienced power users are likely to adjust them – most will likely stick with default settings or only adjust easier-to-understand options. However, the presence of low-level parameters provides an opportunity to experiment and gain greater control over assemblies and are therefore appreciated even when unlikely to be used.

Another aspect worth noting is whether an assembler produces useful files other than its final assembly. Canu and NECAT stand out in this respect, as they create corrected and trimmed reads in their pipelines which have low error rates and are mostly free of adapters and chimeric sequences. Canu and NECAT can therefore be considered not just assemblers but also long-read correction tools suitable for use in other analyses.

Assembler summaries

Canu v2.1 was the slowest assembler and suffered from large circularisation problems. However, it was quite reliable and did well with plasmids. Its main strength is in its configurability, so power users who are willing to learn Canu’s nuances may find that they can tune it to fit their needs. However, it is probably not the best choice for users wanting a quick and simple prokaryote genome assembly.

Flye v2.8 was a strong and well-balanced performer in our tests: reliable, robust and good with plasmids. It also produced the fewest large-scale indel errors in its assemblies. However, it often deleted some sequence (usually on the order of tens of bases) when circularising contigs and had the highest RAM usage of assemblers tested.

Miniasm/Minipolish v0.3/v0.1.3 was not the most reliable assembler but was fairly robust to read set parameters. Its main strength is that it was the most likely to consistently achieve perfect contig circularisation (as this is a specific goal of its polishing step). It was also one of the better assemblers for plasmids, especially regarding clean circularisation of plasmid sequences.

NECAT v20200803 performed reliably with chromosome assembly in the real read sets and was second only to Miniasm/Minipolish for contig circularisation. However, it failed to assemble many plasmids and was cumbersome to run.

NextDenovo/NextPolish v2.3.1/v1.3.1 was resource-efficient and very good at completing chromosomes in both simulated and real read sets, but it performed poorly on plasmid assembly. It was also the most cumbersome assembler to run, requiring multiple commands.

Raven v1.3.0 was reliable and robust for chromosome assembly and used very little RAM. However, it suffered from worse circularisation problems than Flye (often deleting hundreds of bases) and wasn’t good with small plasmids.

Redbean v2.5 assemblies tended to have glitches in the sequence which caused breaks in contiguity, making it perform poorly in both reliability and robustness. This makes it a less-than-ideal choice for long-read prokaryote read sets.

Shasta v0.7.0 was the fastest assembler tested and had low RAM usage, but it had the worst robustness and second-worst reliability. It is therefore more suited to assembly of large genomes in resource-limited settings (the use case for which it was designed) than it is for prokaryote genome assembly.

Underlying data

Figshare: Read sets. https://doi.org/10.26180/5df6f5d06cf04¹⁶.

These files contain the input read sets (both simulated and real) for assembly.

Figshare: Reference genomes. https://doi.org/10.26180/5df6e99ff3eed¹⁷.

This file contains the reference genomes against which the long-read assemblies were compared. For the simulated read sets, these genomes were the source sequence from which the reads were generated.

Figshare: Assemblies. https://doi.org/10.26180/5df6e2864a658³⁴.

These files contain assemblies (in FASTA format), times and terminal outputs for each of the assemblers.

Extended data

Zenodo: Long-read-assembler-comparison. https://doi.org/10.5281/zenodo.2702442⁹.

This project contains the following extended data:

Extended data are also available on GitHub.

Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).