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 and PacBio 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 structural differences between the two assemblies as determined by a Minimap2 alignment¹⁵. This left six isolates for inclusion.

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 six long-read assemblers: Canu v1.9, Flye v2.6, Miniasm/Minipolish v0.3, Raven v0.0.5, Redbean v2.5 and Shasta v0.3.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 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.

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 64 GB of RAM (m3.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 64 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. 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)⁹.

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’. If the assembly terminated due to exhausting the available RAM, it was given a status of ‘out of memory’. Computational metrics were also observed for each assembly: time to complete and maximum RAM usage.

Reliability

When considering only the chromosome, Raven was the most reliable assembler, closely followed by Flye – both were able to complete the chromosome in over three-quarters of the real read sets (Figure 2A). If plasmids are also considered, then Flye was the most reliable assembler. Miniasm/Minipolish and Canu were moderately reliable, completing over half of the real read set chromosomes. Redbean and Shasta were the least reliable and completed less than half of the chromosomes.

Robustness

Flye, Miniasm/Minipolish and Raven were the most robust assemblers, able to complete over half of the assemblies attempted with the simulated read sets (Figure 1A). Flye and Redbean performed best in cases of low read depth, able to complete assemblies down to ~10× depth (Extended data, Figure S7A)⁹. Raven performed the 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 well-tolerated by the assemblers except for Redbean and Shasta (Extended data, Figure S7G)⁹.

Identity

In our real read tests, Canu achieved high sequence identity on PacBio reads, Miniasm/Minipolish and Raven did well on ONT reads, and Flye did well on both platforms (Figure 2C). 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, Miniasm/Minipolish and Raven performing best (Figure 1C).

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. Post-assembly polishing tools (including Racon²⁶, Nanopolish⁷, Medaka³³ and Arrow³⁴) are routinely used to improve the accuracy of long-read assemblies³⁵, and 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.

Resource usage

Canu was the slowest assembler tested on both real (Figure 2D) and simulated (Figure 1D) 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 occasionally hit our 64 GB limit for simulated read sets. 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.

Miniasm/Minipolish, Raven and Redbean were comparable in performance, typically completing assemblies in less than 10 minutes and with less than 16 GB of RAM. 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. Shasta was the fastest assembler and had the lowest memory usage.

Circularisation

Of all assemblers tested, Miniasm/Minipolish was the only one to regularly achieve exact circularisation (contiguity=100%), due to Minipolish’s polishing pipeline (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). Canu’s contiguities usually exceeded 100%, indicating a large amount (thousands of bases) of start/end overlap. The amount of overlap in a Canu 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 and Shasta produce graph files of their final assembly which can indicate circularity. Canu indicates circularity via the ‘suggestCircular’ text in its contig headers. Raven 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. Redbean and Shasta were least successful at plasmid assembly.

Circularisation of plasmids followed the same pattern as for chromosomes, with only Miniasm/Minipolish consistently achieving clean circularisation (Extended data, Figure S10)⁹. 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, though it occurred with other assemblers as well.

Ease of use

All 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). All were able to take long reads in FASTQ format as input, with the exception of Shasta which required reads to first be converted to FASTA format (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, Flye and Redbean) while others (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, 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 stands out in this respect, as it creates corrected and trimmed reads in its pipeline which have low error rates and are mostly free of adapters and chimeric sequences. Canu can therefore be considered not just an assembler but also a long-read correction tool suitable for use in other analyses.

Assembler summaries

Canu v1.9 was the slowest assembler and not the most reliable or robust. Its 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.6 was an overall strong performer in our tests: reliable, robust and good with plasmids. However, it requires a genome size parameter, tended to delete some sequence (usually on the order of tens of bases) when circularising contigs and could be excessive in its RAM usage when assembling simulated read sets.

Miniasm/Minipolish v0.3 was not the most reliable assembler but was fairly robust to read set parameters. Its main strength is that it was the only assembler to consistently achieve perfect contig circularisation (as this is a specific goal of its polishing step). Also, it does not require a genome size parameter to run, which makes it easier to run than Canu, Flye or Redbean for unknown genomes.

Raven v0.0.5 was the most reliable and robust assembler for chromosome assembly. However, it suffered from worse circularisation problems than Flye (often deleting hundreds of bases) and wasn’t good with small plasmids. Like Miniasm/Minipolish, it does not require a genome size parameter.

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, combined with its erratic circularisation performance and requirement to specify genome size, make it a less-than ideal choice for long-read prokaryote read sets.

Shasta v0.3.0 was the fastest assembler tested and used the least RAM, but it had the worst reliability and robustness. 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).