Metagenomic assembly is the main bottleneck in the identification of mobile genetic elements

MetaMobilePicker

MetaMobilePicker, available at http://metamobilepicker.nl, is an open-source software pipeline linking detected ARGs to MGEs in metagenomes. MetaMobilePicker was written in Python 3.10 using the workflow management system Snakemake (Mölder et al., 2021) and was designed to run primarily on high-performance computation infrastructures, as MetaSPAdes (Nurk et al., 2017) is resource-intensive (Zhang et al., 2023). Versions of required software are listed in Table S1. The pipeline is installable as a Python package or by using conda from the bioconda channel. The pipeline can be divided into four parts: preprocessing metagenomic short-read sequences, genome fragment reconstruction, ARG annotation, MGE identification and output construction linking ARGs to MGEs. An overview of the overall workflow is displayed in Fig. 1. In the preprocessing steps, sequencing reads are deduplicated for PCR artifacts and filtered for read quality and (host) contamination. This is done using the QC module of Metagenome-Atlas (Kieser et al., 2020). Next, a metagenome assembly is performed using metaSPAdes (Nurk et al., 2017) on high-quality reads with default parameters as specified in the configuration file. Very short contigs of length less than 1 kbp are discarded for their limited use in downstream analyses. MGEs are identified using PlasClass (Pellow, Mizrahi & Shamir, 2020), ISEScan (Xie & Tang, 2017) and DeepVirFinder (Ren et al., 2020), respectively. ARGs are then annotated using ABRicate (Seemann, 2023) on the ResFinder (Zankari et al., 2012) database. Finally, all output files are combined, linking ARGs to MGEs if present on the same contig.

Validation using simulated data

In order to validate the results of the different tools included in MetaMobilePicker (PlasClass, ISEScan, DeepVirFinder and ABRicate), we constructed a simplified simulated metagenomic dataset. We selected the genomes of seven highly resistant bacterial species from the PATRIC database (Wattam et al., 2014 accessed on March 3rd 2021). For this dataset, we selected a representative completed genome from E. coli, S. aureus, Enterococcus faecalis, Mycobacterium tuberculosis, Salmonella enterica, Klebsiella pneumoniae and Acinetobacter baumannii. In addition to these genomes, five phage genomes specific to the selected genera were selected. Selected organisms and associated plasmids and phages can be found in Table 1. Using the selected genomes, we simulated a dataset of 20 million paired-end reads (10 M per end) of 150 bp using InSilicoSeq (Gourlé et al., 2019) using the MiSeq error profile. To test if 10 million simulated reads were sufficient to resolve the complexity of the dataset without introducing large gaps when comparing the metagenomic assembly to the reference genomes, we simulated 40 million reads per end and subsampled them at intervals of 5 million reads between 5 and 30 million. We used MetaQuast (Mikheenko, Saveliev & Gurevich, 2016) to calculate the genome completeness per genome. In order to simulate different relative abundances of different genomes, we generated abundance profiles for all species using a lognormal distribution. Additionally, a copy number for the plasmids was determined using a geometric distribution with probability $P=min\left(1,\frac{log10\left(L\right)}{r}\right)$ where L is the length of the plasmid in bp, to simulate smaller plasmids having a higher probability of having a higher copy number than larger plasmids and r is a reduction factor based on the order of magnitude of the sequence length. For plasmids, we set r = 7 and for phages we set r = 5. This copy number was multiplied by the abundance of the corresponding genome, and the abundances were normalized to sum to one. The complete genome assemblies, plasmids and phages were annotated for ARGs using ABRicate with the ResFinder. Additionally, ISEScan was used to identify the present IS elements present in the sequences. To backtrace the contigs assembled during the pipeline’s run associated with the plasmids and phages, we used Minimap2 (Li, 2018) as part of the MetaQuast workflow. Contig alignments shorter than 65 bp and with identity lower than 95% were removed, as per MetaQuast default.

To assess the performance of plasmid and phage classification, we calculated precision, defined in Eq. (1), recall, defined in Eq. (2), and the F1-score, defined as the harmonic mean between precision and recall. To assess the performance of the annotation of IS elements and ARGs, we calculate the sensitivity as defined in Eq. (3). (1)${\displaystyle Precision=\frac{Correctlypredictedcontigs}{Correctlypredictedcontigs+Falsepositives}}$ (2)${\displaystyle Recall=\frac{Correctlypredictedcontigs}{Correctlypredictedcontigs+Falsenegatives}}$ (3)${\displaystyle Sensitivity=\frac{CorrectlypredictedISelementsorARGs}{TotalnumberofISelementsorARGs}.}$

Bases not covered in simulated data

To infer the number of bases not covered by the simulated reads per reference genome, we used BWA MEM (Li, 2013) to align the simulated reads to the reference genomes. Using Samtools (Li et al., 2009), we inferred the coverage depth at each position on the reference genomes, and considered bases with a depth of 0 as not covered, as no reads mapping to these bases were found in the dataset.

Establishing a baseline to validate the performance of the detection of MGEs and mARGs in metagenomes

We used reference genomes of seven bacterial species, including nine plasmids and five genus-specific phages, to simulate metagenomic reads as a validation data set for the tools included in MetaMobilePicker (Table 1). To establish a baseline for the performances of the detection tools for IS elements, plasmids, phages and ARGs, we ran ISEScan, PlasClass, DeepVirFinder and ABRicate with the Resfinder database on the completely sequenced reference whole-genome sequences as available in RefSeq and cross-checked with the available NCBI annotation. This resulted in the detection of 417 IS elements in the reference genomes, distributed among the seven bacterial genomes with most IS elements detected in the E. coli chromosome (n = 78) and least in the smallest of the two E. coli plasmids (n = 3) (Fig. 2). PlasClass identified plasmids and chromosomes without errors and misclassified two phages as plasmids (out of a total of five) based on a cut-off score of 0.8. DeepVirfinder showed a perfect classification of phages and non-phages on the reference data. Finally, to establish a baseline for the AMR identification step on the finished reference genomes, we annotated the reference genomes using the ResFinder database with ABRicate. This identified a total of 55 predicted ARGs from 38 different families. Of these 55, 28 were only found once by ABRicate using the default parameters, five ARGs were found twice (tetB, ermA, ant(9)-Ia, aph(3′)-III and blaSHV-12), three ARGs were found three times (blaTEM, tetA and sul2), and two ARG (aph(6)-Id and aph(3″)-Ib) were found four times. Two genes (tet(B) and aph(3″)-Ib) were found with two different alleles. In the cases where multiple copies were identified, the ARGs were annotated in the same replicon in two cases, in different replicons of the same species in one case and to different replicons in two or more species in seven cases. Of these 55 ARGs, 29 genes are present on plasmids, belonging to 20 different gene families. ARGs present in the reference genomes, copy number and mobility of replicons are shown in Table S2.

Detection of MGEs and ARGs in the metagenome-assembled contigs has a moderate accuracy

From the 20 different genomic replicons from seven bacterial species, we simulated 20M metagenomic paired-end reads (10 million per end) using InSilicoSeq with different relative abundances (Table 1). The simulated dataset of 10M paired end reads was used to run MetaMobilePicker. After running the QC and assembly modules, the identification steps identify three types of MGEs: (i) plasmids, (ii) IS elements and (iii) phage sequences.

Metagenomes often have a read depth ranging from 10 million - 70 million per sample (Gweon et al., 2019), albeit with a far larger complexity than the validation data set simulated here. To test to what extent the read depth influenced the resolution of the genomes in our dataset, we simulated 40 million reads per end and subsampled them at 5, 10, 15, 20, 25 and 30 million reads and assembled them. We then used MetaQuast to calculate the completeness per genome. The average genome completion per sample depth is shown in Fig. 3. This figure shows that genome coverage completeness did not improve much after 10 million reads, with a maximum at 20 million reads, although the difference between 10 and 20 million reads was small (0.14%). We concluded that 10 million simulated reads were sufficient to resolve the complexity of this specific dataset. To maintain the defined relative abundances, and avoid this being skewed by the sub-sampling process, we used InSilicoSeq to simulate a definitive dataset with 10 million reads per end, totalling 20 million reads. Running the Metagenome Atlas (Kieser et al., 2020) QC module on these reads resulted in a total of 9,797,445 reads per end and 42,065 orphaned reads. The built-in MetaSPAdes assembly step resulted in 1,513 contigs, of which 887 contigs were larger than 1 kbp with an N50 of 82,245 and an L50 of 97. The total assembly length accumulated to 29,747,115 bp. Only contigs longer than 1 kbp length were used for the classification of MGEs and ARGs. To detect IS elements in the simulated metagenome validation data, we ran the ISEscan (Xie & Tang, 2017) module of MetaMobilePicker on the 887 contigs longer than 1 kbp. ISEScan uses a prebuilt HMM database to identify insertion sequences. In total, 241 IS elements in 144 metagenomic contigs were identified compared to 417 IS elements in the reference genomes, resulting in a sensitivity of 0.585. In order to predict plasmids in the metagenomic contigs, PlasClass as part of the MetaMobilePicker pipeline was applied to contigs larger than 1 kbp, with a cut-off classification score of 0.8. PlasClass uses a deep neural network and does not require a reference database during runtime. This resulted in 111 contigs predicted as plasmids (Fig. S1). Of these 111, 63 contigs were correctly predicted as plasmids and 756 contigs were correctly predicted as chromosomal. Additionally, 48 and 17 contigs were falsely predicted as plasmid and chromosomal, respectively, resulting in a precision of 0.57, recall of 0.79, and an F1 score of 0.66. Phages were predicted with the MetaMobilePicker module for DeepVirFinder (Ren et al., 2020), which uses a logistic regression model to predict phages and was tested as one of the best-performing phage prediction tools for metagenomics in a recent study (Ho et al., 2023). In order to measure the performance of DeepVirFinder on our metagenomic assembly, we cross-referenced the contigs predicted by DeepVirFinder with a classification score greater than 0.95, with the contigs that originated from the five phages in the dataset. This analysis showed 27 predictions, of which five originated from the phages added to our community. These five phage contigs were the full-length assemblies of these phages. Of the 22 contigs not originating from our added phages, functional annotation using blastx (Camacho et al., 2009) showed 14 having a direct link to phage DNA, most likely originating from prophages. Another four contigs were putatively linked to phages, containing hits not exclusively associated with phages. The remaining four contigs showed no clear link to phage DNA. The false positive contigs were notably shorter than the true positive phage contigs. Of the 22 false positives, three were larger than 10 kbp, with a median of 2,568 bp, indicating limited room for full-length genes and genomic context on many of the contigs. As we were unable to determine the number of prophage- or phage-related genes not identified by DeepVirFinder, we did not take the 18 phage-related genes into account when calculating the classification metrics. This resulted in a recall of 1.0, a precision of 0.555 and an F1-score of 0.713. To benchmark the performance of ARG annotation on our validation set, we cross-referenced the metagenomic hits with the genomic hits (Fig. 4). Of the 55 ARGs in the reference genomes, 38 were identified on the correct reference chromosome or plasmid in the metagenome-assembled contigs, resulting in a sensitivity of 0.697. For 16 of the remaining 17 genes, we identified the gene once, while they were present in multiple copies in the reference genomes. These 16 hits comprise 10 unique ARGs. The remaining hit was not found in the metagenomic contigs. Additionally, for the genes for which more than one allele was present (tet(B) and aph(3″)-Ib), only one allele was found. Only one sul2 gene was found more than once in the metagenomic assembly, which was found twice compared to an expected copy number of three in the reference genomes. The identified genes are displayed in Table S2. An overview of the classification metrics for plasmids and phages can be found in Table 2. An overview of the annotation metrics for IS and ARGs can be found in Table 3.

To find potential causes for the observed percentage of IS elements that were not present in the contigs, for the low precision of the plasmid identification, and for the low number of annotated ARGs, we further investigated read coverage, the influence of contig length cutoffs and assembly of the different MGEs and ARGs.

Lacking read coverage explains a minority of missing MGEs in the assembled metagenome

To investigate whether missing or lower read coverage of the simulated data could explain the low precision observed for the identification of MGEs and mARGs in our metagenome validation data, we mapped the simulated reads against the reference genomes and validated the depth and breadth of coverage of the MGE annotated positions. This showed a varying range of percentage of bases not covered between 0.001% (E. coli plasmid) and 3.67% (E. faecalis plasmid), with an average of 0.48% missing bases with a standard deviation of 1.16% (Table S3). This indicates that only a small portion of reference bases was not covered by the simulated reads. None of these bases with missing coverage aligned to ARGs. The E. faecalis chromosome contains a total of 24,900 bases not covered, possibly caused by the low simulated abundance of this species. These are spread out over the genome, and cause eight IS elements to have a contiguous region of bases not covered of 100 bp. Additionally, two more IS elements have contiguous regions not covered larger than 50 bp (Table S4), both in the M. tuberculosis genome. The bases not covered predominantly align to the chromosome sequences, with only one phage (containing a singular base not covered), and two of the plasmids having bases not covered. These plasmids are the aforementioned E. faecalis plasmid, and one E. coli plasmid having a singular base not covered. In contrast, each of the seven chromosome sequences have bases not covered ranging from 178 (0.006% of the genome) to 24,900 (0.96% of the genome) bases. Lacking read coverage was therefore, at most, responsible for 3.67% of the overall recall for plasmids and for less than 0.01% of the overall recall for phages. Especially for IS in the E. faecalis chromosome, the missing read coverage caused difficulty in assembly and subsequent annotation. The IS elements containing regions not covered of 100 bp were not assembled into contigs in both the per-species assembly and the complete assembly.

Smaller contig length cutoffs do not improve recall of MGEs

We investigated if the contig length cutoff of 1 kbp, as applied by MetaMobilePicker could influence the recovery of MGEs. We analyzed the benchmark data set with a contig, length cutoff of 500 bp instead of MetaMobilePicker’s default of 1 kbp, annotated them as described previously and compared all outcomes. This resulted in 954 contigs (>500 bp), instead of 887 (>1 kbp). Of the additional contigs ranging between 500 and 1,000 bp, none contained an ARG. Furthermore, the contigs with a length between 500 and 1,000 bp that were predicted as plasmids were predominantly of chromosomal origin (11 of 15 contigs), increasing the number of false positives more than the number of true positives and lowering the precision for plasmid prediction from 0.57 (>1 kbp) to 0.54 (>500 bp). Likewise, four contigs with a length between 0.5 and 1 kbp were predicted as phages but were of chromosomal origin. Of these sequences, two were homologs of phage proteins, and the other two had a potential association with phage DNA. Therefore, all four of these sequences were subsequently not counted in the calculation of the phage metrics. This suggested that lowering the sequence length threshold increased the number of false positive predictions, which was not compensated by the increase of true positive predictions.

Complexity of metagenome assembly is not causal for assembly performance

When split by species, the prediction of plasmids by PlasClass exhibited a false-positive bias, especially for contigs from E. coli (Fig. 5). Of the 48 false positive plasmid predictions, 30 originated from the E. coli genome. The false positive predictions of phages did not show this bias, and no false positive annotations were found for ARGs and IS elements. To investigate if this overrepresentation was caused by a highly fragmented genome assembly, we separated the simulated reads per genome and used each set of reads as input for MetaMobilePicker separately. The resulting set of classified plasmid contigs was compared to the reference genome sequences and the metagenomic contigs using metaQuast. The classification scores for each species are displayed in Supplemental Table S5. Most notably, the number of false positive plasmid predictions originating from the E. coli genome for the single-species-only assembly was as high as in the metagenomic assembly with 30 false positive predictions. Comparing the false positive-predicted contigs from the single-species-only assembly to the metagenomic assembly showed 29 identical pairs, and one pair where the contig originating from the metagenomic assembly was 274 bp shorter (total length = 1,622 bp) than the contig from the reference sequences. For the other species, the single-species-only assembly also showed largely the same outcome as the metagenomic assembly with regard to plasmid-predicted contigs, with the exception of S. enterica. In S. enterica, the assembly of the plasmid improved notably in the single-species-only assembly (metagenome: 16 mapping contigs, average length 6.3 kbp, single-species-only: four mapping contigs, average length 28.1 kbp). However, as the S. enterica plasmid represented only one false positive prediction in the metagenome assembly, the S. enterica F1 score for metagenome or single-species only differed because of the higher degree of fragmentation of the metagenomic assembly (F1_meta = 0.96, F1_{singlespecies} = 0.8). Comparison of the wrongly predicted sequence lengths shows similar results (bp_meta = 9.8 kbp, bp_{singlespecies} = 11.0 kbp). Based on this, we note that the increased complexity of metagenomic assembly does not substantially influence the recovery of MGEs in the benchmark dataset.

Assembly collapses identical insertion sequences and ARGs

To test if the low sensitivity of the identification of IS elements and ARGs was caused by difficulty assembling these elements, we first tested if all IS elements in the reference genomes were represented in the metagenomic contigs (Fig. S2). We mapped the metagenomic contigs larger than 1 kbp to the reference genomes. This showed that for 324/417 IS elements (77.70%) annotated in the reference genomes, a corresponding metagenomic contig was mapped, whereas a matching metagenomic contig was missing (less than 20% of the IS elements covered) for 93/417 IS elements present in the reference genomes (22.30%). This includes the 10 IS elements not assembled due to missing read coverage. Of the 324 reference IS elements with a corresponding metagenomic contig, 90 (27.80%) corresponded to 16 metagenomic contigs which mapped to multiple IS elements in the reference genomes. These 16/417 metagenomic contigs are therefore likely to represent assembly collapses. Furthermore, 16 IS elements annotated in the reference genomes (4.90%) corresponded only partially (between 20% and 80% of the IS elements covered) to a metagenomic contig. The remaining 221 IS elements in the reference genome (69.70%) corresponded to a metagenomic contig uniquely mapped to that IS element. To contrast the detection sensitivity of ISEscan before and after metagenome assembly, we also calculated the maximum achievable amount of IS elements detected in the metagenome-assembled contigs. To this end, we discounted all IS elements that were lost or only partially assembled during metagenome assembly. When taking only the IS elements into account that were present in the metagenomic contigs, and subtracting the duplications of the ambiguously mapped contigs, the maximal number of IS contigs ISEscan could potentially identify was 248. ISEscan detected 231, thereby achieving a sensitivity of 93.15%. From this, we conclude that collapse of IS elements during the assembly process was a major contributor to the sub-optimal identification sensitivity of IS elements.

As mentioned above, 16 of the ARGs in the reference genomes were not identified with the correct copy number. Further investigation of the metagenomic contigs mapped to the reference genomes showed that 13 of these genes did not have a corresponding metagenomic contig. The remaining three genes mapped partially to a contig but lacked annotation. This suggested that similar to the IS, ARGs were collapsed during the metagenomic assembly step.