metaFlye: scalable long-read metagenome assembly using repeat graphs

Abstract

Long-read sequencing technologies have substantially improved the assemblies of many isolate bacterial genomes as compared to fragmented short-read assemblies. However, assembling complex metagenomic datasets remains difficult even for state-of-the-art long-read assemblers. Here we present metaFlye, which addresses important long-read metagenomic assembly challenges, such as uneven bacterial composition and intra-species heterogeneity. First, we benchmarked metaFlye using simulated and mock bacterial communities and show that it consistently produces assemblies with better completeness and contiguity than state-of-the-art long-read assemblers. Second, we performed long-read sequencing of the sheep microbiome and applied metaFlye to reconstruct 63 complete or nearly complete bacterial genomes within single contigs. Finally, we show that long-read assembly of human microbiomes enables the discovery of full-length biosynthetic gene clusters that encode biomedically important natural products.

Extended Data Fig. 1. Information about metaFlye, Flye, Canu, miniasm, and wtdbg2 assemblies of the individual genomes in the SYNTH64 dataset.

NGA50 (in megabases) and reference coverage (in percentages) reported for all genomes from the SYNTH64 dataset. Genomes are ordered in the increasing mean NGA50 across all assemblers. Challenging genomes that have closely related species or strains in the metagenome are marked with (!). Grey bars on the NGA50 plot represent the length of the longest chromosome in the reference sequence for each genome (a theoretical upper bound for NGA50). NGA50 is shown in logarithmic scale (not shown for values lower than 100 kb or if the reference coverage is below 50%). The full metaQUAST report for the SYNTH64 dataset is provided in Supplementary Table 1.

Extended Data Fig. 2. NGAx plots for the mock community datasets (HMP mock, ZymoEven GridION, ZymoLog GridION).

NGA(x) is the statistic computed for contigs that are broken at their misassembly breakpoints (if any). NGA(x) is the highest possible number L such that all broken contigs that are longer than L cover at least X% of the reference. Plots were generated by metaQUAST using all available references for each dataset. Flye failed to assemble the ZymoLog datasets due to poor k-mer indexing (Methods).

Extended Data Fig. 3. Base-pair accuracy analysis for assemblies of the mock community datasets (HMP, ZymoEven GridION, and ZymoLog GridION).

Heatmaps showing the number of mismatches and short indels per 100 kbp for each species reference, computed using metaQUAST. Blue and red colors correspond to the values higher and lower than the median, respectively. Statistics were not computed for genomes with no assembled sequence (“-” symbol). Flye failed to assemble the ZymoLog datasets due to poor k-mer indexing (Methods).

Extended Data Fig. 4. The ORF lengths distribution and the GC content distribution of metaFlye and Canu assemblies of the sheep microbiome.

The ORF length distribution suggests similar base-level accuracy for both assemblies.

Extended Data Fig. 5. Taxonomic assignments of sheep microbiome assemblies.

(a) metaFlye contigs assignment at the phylum level visualized with BlobTools. (b) Length distributions of metaFlye and Canu contigs within each assigned superkingdom.

Extended Data Fig. 6. Statistics of simple bubbles for the metaFlye assemblies human gut and cow rumen.

(Left) the human gut dataset with 615 bubbles, and (right) the cow rumen dataset with 1510 bubbles. Bubble counts exclude loops, and include roundabouts with two edges.

Extended Data Fig. 7. Analysis of sequence overlap between 19 human gut samples.

Multi-way sequence alignments were computed using SiebliaZ. (left) The proportions of unique and shared sequences in each sample. An assembled segment within a sample is called unique if it has no alignments against sequence from any other samples. Otherwise, the segment is shared. (right) The total amount of sequence for each multiplicity bin. A sequence fragment belongs to the multiplicity bin X if it is shared by exactly X samples.

Figure 1.. metaFlye repeat annotation and examples of simple bubbles, superbubbles, and roundabouts.

(a) The subgraph of an assembly graph formed by four distinct genome sub-paths. Repeat and unique edges are shown in color and black, respectively. metaFlye identifies edges X, Y, and Z as repetitive by analyzing the distinct read-paths through the sub-graph. (b) A simple bubble formed by two strains. (c) A superbubble formed by three strains. (d) A roundabout formed by two strains, one of which shares a repeat with a different region of the metagenome.

Figure 2.. Information about Canu, Flye, metaFlye, miniasm, and wtdbg2 assemblies of the individual genomes in the SYNTH181 dataset.

Assembled fraction and NGA50 are reported for all 181 reference genomes from the simulated dataset. Genomes are ordered in the decreasing mean assembled fraction (left) and NGA50 (right) across five assemblers. NGA50 is the statistic computed for contigs that are broken at their misassembly breakpoints (if any). NGA50 is the highest possible number L such that all broken contigs that are longer than L cover at least 50% of the reference. NGA50 is not shown for values lower than 10 kbp or if the reference coverage is below 50%. 77 (metaFlye), 141 (Flye), 109 (Canu), 106 (miniasm) and 109 (wtdbg2) NGA50 values were filtered this way. The full metaQUAST report is provided in Supplementary Table 2.

Figure 3.. Per-species reference coverage and NGA50 statistics for the mock community datasets (HMP, ZymoEven GridION, ZymoLog GridION) computed using metaQUAST.

The read coverage for each species is given in the brackets after the species name. NGA50 values are not reported for assemblies with reference coverage below 50%. Blue and red colors correspond to the values higher and lower than the median, respectively. Flye failed to assemble the ZymoLog datasets due to poor k-mer indexing (Methods). Extended Data Figure 3 provides the base-pair quality analysis for the same datasets.

Figure 4.. Information about strains in the sheep microbiome revealed by metaFlye.

(a) An assembly graph of a single connected component in the sheep microbiome dataset before strain collapsing (visualized using Bandage). The component represents a bacterial genome of the Clostridia class with 92% conserved marker completion (computed using CheckM). There are 20 simple bubbles (shown in green) and 10 superbubbles (shown in yellow) that account for 1.2 Mbp out of 2.4 Mbp long genome. (b) Distribution of length and branch sequence identities of 1141 bubbles (excluding loops and including roundabouts with only two edges) in the sheep microbiome assembly. The length is defined as the length of the longest branch in a simple bubble.