
==== Front
Bioinformatics
Bioinformatics
bioinformatics
Bioinformatics
1367-4803
1367-4811
Oxford University Press

10.1093/bioinformatics/btae398
btae398
Microbiomes and Planetary Health
AcademicSubjects/SCI01060
skandiver: a divergence-based analysis tool for identifying intercellular mobile genetic elements
Zhang Xiaolei Brian Department of Biological Sciences, Carnegie Mellon University, Pittsburgh, PA 15213, United States
Department of Computational Biology, Carnegie Mellon University, Pittsburgh, PA 15213, United States

Oualline Grace Department of Biological Sciences, Carnegie Mellon University, Pittsburgh, PA 15213, United States
Department of Computational Biology, Carnegie Mellon University, Pittsburgh, PA 15213, United States

Shaw Jim Department of Mathematics, University of Toronto, Toronto, ON M5S2E4, Canada

Yu Yun William Department of Computational Biology, Carnegie Mellon University, Pittsburgh, PA 15213, United States
Department of Mathematics, University of Toronto, Toronto, ON M5S2E4, Canada
Department of Computer and Mathematical Sciences, University of Toronto at Scarborough, Toronto, ON M1C1A4, Canada

Corresponding author. E-mail: ywyu@cmu.edu (Y.W.Y.)
† Joint first authors.

9 2024
04 9 2024
04 9 2024
40 Suppl 2 Proceedings of ECCB2024 ii155ii164
© The Author(s) 2024. Published by Oxford University Press.
2024
https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

 

Motivation: Mobile genetic elements (MGEs) are as ubiquitous in nature as they are varied in type, ranging from viral insertions to transposons to incorporated plasmids. Horizontal transfer of MGEs across bacterial species may also pose a significant threat to global health due to their capability to harbor antibiotic resistance genes. However, despite cheap and rapid whole-genome sequencing, the varied nature of MGEs makes it difficult to fully characterize them, and existing methods for detecting MGEs often do not agree on what should count. In this manuscript, we first define and argue in favor of a divergence-based characterization of mobile-genetic elements.

Results: Using that paradigm, we present skandiver, a tool designed to efficiently detect MGEs from whole-genome assemblies without the need for gene annotation or markers. skandiver determines mobile elements via genome fragmentation, average nucleotide identity (ANI), and divergence time. By building on the scalable skani software for ANI computation, skandiver can query hundreds of complete assemblies against >65 000 representative genomes in a few minutes and 19 GB memory, providing scalable and efficient method for elucidating mobile element profiles in incomplete, uncharacterized genomic sequences. For isolated and integrated large plasmids (>10 kb), skandiver’s recall was 48% and 47%, MobileElementFinder was 59% and 17%, and geNomad was 86% and 32%, respectively. For isolated large plasmids, skandiver’s recall (48%) is lower than state-of-the-art reference-based methods geNomad (86%) and MobileElementFinder (59%). However, skandiver achieves higher recall on integrated plasmids and, unlike other methods, without comparing against a curated database, making skandiver suitable for discovery of novel MGEs.

Availability and implementation

https://github.com/YoukaiFromAccounting/skandiver

Natural Sciences and Engineering Research Council of Canada 10.13039/501100000038 RGPIN-2022-03074 Carnegie Mellon University 10.13039/100008047 NSERC CGS-D
==== Body
pmc1 Introduction

Mobile genetic elements (MGEs) are snippets of sequences that can move around (either intracellularly within a genome, or intercellularly across species) independent of the host genomes (Shapiro 1983). There exist many biological mechanisms associated with MGEs (viral phages, plasmids, transposons, etc.) but the end results of especially intercellular trafficking of genetic material is often an accelerated evolutionary time scale (Frost et al. 2005).

To cite just one topical example, antibiotic resistance has been a critical challenge in healthcare, due to the overuse and misuse of antibiotics (Durrant et al. 2020). Antibiotic resistance genes thus contribute substantial fitness to the bacteria, and are strongly selected for by evolutionary pressures (Davies and Davies 2010). Bacteria often share genes through horizontal gene transfer (Forster et al. 2022), and it is thought that MGEs are instrumental in the ability of so-called ESKAPE pathogens to quickly develop resistance (Partridge et al. 2018). Identifying and detecting these MGEs is essential to gain insights into their formation and evolution, which may provide a foundation for developing preventative treatments to reduce the rise of antibiotic-resistant strains (Johansson et al. 2021). Furthermore, the rapid evolution of these resistance genes highlights the need to detect MGEs without needing to reference a library of known contigs of antibiotic resistance.

Despite the importance of MGEs, there exist limitations with current bioinformatics software. One might quite reasonably hope that with the plethora of sequencing data now available to practitioners, we would be able to just throw data at the problem, but there are always trade-offs (Berger and Yu 2023). Many tools rely on gene-based classification methods, utilizing database searches and sequence homology in order to identify target sequences. One such tool, MobileElementFinder (Johansson et al. 2021), aligns assembled contiguous sequences (contigs) from a library of known mobile elements to input bacterial DNA sequences, in order to predict MGEs in those sequences. Another tool, geNomad, is a classification framework that can either use marker genes or a deep learning model to identify plasmids and viruses from gene annotations (Camargo et al. 2023). However, both methods are limited to finding MGEs that are similar to previously discovered and annotated MGEs—MobileElementFinder explicitly through alignment, and geNomad implicitly through the machine learning model training set. Furthermore, although fairly efficient, alignment and machine learning are both somewhat slower primitives compared to some of the more modern sketching-based sequence similarity techniques that we take advantage of in this manuscript (Shaw and Yu 2023).

To address the limitation of prior tools being restricted to MGEs similar to already annotated ones, we target specifically intercellular MGEs from first principles using evolutionary divergence times. Evolutionary divergence is typically measured in units of MYA (million years ago), corresponding to our current best guess based on sequence analysis of proteins and genomes as to the timeframe of the last common ancestor between two species (Doolittle et al. 1996). Our key insight is that mobile elements are by definition untethered to the background evolutionary process governing the rest of the genome—intercellular MGEs appear, mostly unchanged, in otherwise unrelated species. We can thus predict that a sequence from a genome is a mobile element if the set of genomes it appears in have high mutual evolutionary divergence times, but that sequence does not appear in more closely related genomes. Figure 2 illustrates this idea graphically through an example from Acinetobacter baumannii. (a) Random sequences from a genome typically do not match any other species. (b) Conserved elements show up in many other species, but primarily only those are related (i.e. have low divergence times from the reference genome). (c) Mobile elements show up even in many distantly-related species with high evolutionary divergence times.

In this manuscript, we introduce our method skandiver for finding intercellular MGEs. Instead of using a previously annotated library of MGEs, skandiver measures the evolutionary divergence of the species a putative MGE has jumped between, enabling detection of MGEs database-based tools like geNomad and MobileElementFinder cannot. skandiver is built around our prior tool skani (Shaw and Yu 2023) and the TimeTree of Life (Kumar et al. 2022). Skani is able to take a sequence and quickly and efficiently calculate average nucleotide identity (ANI) and the fraction of genomes aligned to one another (AF) of metagenomic sequences against a large dataset of genomes. This allows us to swiftly observe what chunks of sequences may map to multiple species of bacteria. We utilize this by looking for sequences that have strong matches to a variety of bacterial species (as measured by TimeTree of Life evolutionary divergence times), indicating that it may be mobile in nature.

Using skandiver, we can detect MGEs in genomic sequences, and produce a summary table of their location, as well as the divergence time between the species that it mapped to. Researchers input assemblies, which are then fragmented and fed through skani, to find where those fragments may occur in the GTDB database (Parks et al. 2022). Then, fragments with strong matches to more than one species are taken, and the divergence time between those species is calculated to determine if that fragment is a conserved or mobile element. If the divergence time is low, it is likely a DNA sequence shared among species due to being essential. However, large divergence times point to the possibility of the fragment of DNA traveling between different species of bacteria. We output these MGE sequences and their corresponding data in an easy-to-read format, making it convenient for researchers with minimal coding experience. The speed and efficiency of utilizing this wrapper script with skani provide a versatile tool that operates with minimal memory for researchers to keep pace with the rapid trafficking of MGEs. Through its application, researchers will be able to quickly pinpoint sequences of DNA within bacterial genomes that are putative MGEs.

2 Materials and methods

2.1 Algorithm design

The methodology behind using skandiver to identify potential MGEs is outlined in Fig. 1. The process can be divided into four main steps:

Figure 1. A composite framework for identifying mobile elements from divergence times. skandiver analyzes whole genomes by preprocessing contigs into smaller genome fragments. Using skani, average nucleotide identity (ANI) and align fraction (AF) are found for each fragment. skandiver filters genome fragments by their ANI and AF to a set of representative genomes, and then performs a phylogenetic tree traversal to retrieve divergence times between every successful mapping to a representative genome. skandiver outputs the set of genome fragments that have mapped strongly to evolutionarily distant genomes, thereby linking these fragments with potential mobile genetic elements.

Fragment genomes: Genome assemblies of interest are fragmented/chunked into smaller regions, typically regions 10 000 base pairs in length, depending on the granularity of specific genomic elements to be identified—in this manuscript, all our results are with 10-kb chunks for technical reasons related to skani’s default implementation, so smaller MGEs may be missed.

Skani search against a database of representative genomes: Fragmented genomes are searched against a database of representative genomes using skani, a robust and accurate ANI comparison tool.

Calculate divergence times: skani search results are filtered to only high ANI and AF; we then annotate significant matches with evolutionary divergence times using the TimeTree of Life (Kumar et al. 2022).

Predict genomic elements: we analyze patterns within skani search results and divergence times to output a set of potential MGEs, with genomic position and divergence time information.

2.2 Augmenting skani search for genome fragments

For whole-genome assemblies and other similar metagenome-assembled genomes (MAGs) for which genomes may be contaminated or incomplete, skani search accurately and efficiently computes ANI between a query genome and a reference genome (Shaw and Yu 2023). It does this by chaining a set of exact nonoverlapping k-mer matches between the query and reference genomes to estimate the ANI. However, skani is designed for computing global similarity, rather than local similarity. For MAGs with contigs that are large (>150 000 bp), skani search will only return a hit if there is sufficiently high ANI and alignment between the entire reference and query. For example, when searching five whole-genome contig-assembled query assemblies of Enterococcus faecium against the Genome Taxonomy Database’s (GTDB) (Parks et al. 2022) representative genomes as reference, skani search returns either zero matches between the query genome and the representative genome, or returns one match to the assembly’s original species (E.faecium).

For finding MGEs, we are more interested in local similarity, so we augmented skani search to focus on genome fragments. We thus introduced a preprocessing step in which MAGs are fragmented into nonoverlapping chunks of sequence (10 kb in length by default). This step aimed to optimize the sensitivity of the skani search by providing a segmented representation of overall genomic content within metagenomic datasets, though it does restrict skandiver to only finding MGEs that are >10 kb—also, because of edge effects related to the overlap between chunks and the MGEs, >20 kb length MGEs are needed to guarantee that at least one chunk comprises only MGE sequence. Fragmentation was accomplished using a simple Python script that utilizes the SeqIO and SeqRecord modules from the Biopython library (Cock et al. 2009); the script takes as input a directory containing genomic sequences in fasta format, a desired output directory, and the preferred fragmentation size in nucleotides. The resulting fragments are then concatenated into a singular multifasta file, preserving original genomic identification information as well as fragment location within the genome.

Following genome fragmentation, repeated skani searches of five whole-genome contig-assembled query assemblies of E.faecium against the GTDB representative genome database returned over 1500 matches of over 90% ANI between the query contig fragments and the reference database. More importantly, rather than returning only matches to a query’s original species, fragments now matched to a diverse array of bacterial species. For example, fragment 40 581–42 248 of NZ_JADVBE010000002.1 Enterococcus faecium strain VRE32783 contig 00019, which previously was undetected by plain skani search, exhibited matches to Streptococcus pasteurianus, Blautia argi, Dilemma fastidiosa, Eubacterium ramulus, and many more reference bacteria following genome fragmentation. This is of particular interest for the detection of MGEs, many of which are insertion sequences or plasmids less than 10 kb in size (Khedkar et al. 2022).

2.3 Filtering for potential MGEs

ANI has previously been used to identify insertion sequences between 70 bp to 200 kb by MGEfinder, which required query sequences to share at least 98.5% ANI with reference sequences to be classified as an MGE (Durrant et al. 2020). However, as our analysis method aims to directly identify MGEs from the query sequences themselves without the use of reference MGEs, we used a slightly more lenient requirement of at least 95% ANI and 90% align fraction (fraction of the query genome that aligns to the reference genome)—intuitively, two instantiations of an MGE that would have both mapped to the same reference sequence can be further apart from each other.

Another justification for our thresholds comes from prior work for identification of uncultivated virus segments—using the Minimum Information about any (x) Sequence (MIxS) standard, Roux et al. (2019) suggested thresholds of 95% ANI and 85% alignment fraction for identification of virus segments. Our ANI threshold was thus in line, though we use a slightly higher align fraction threshold. All skani search results were extracted by using the awk Linux command and converted to comma-separated value file format following filtering for ANI and align fraction threshold.

2.4 Estimating divergence time

Divergence time between reference and query species was determined using TimeTree of Life, a public database of evolutionary and divergence times between major clades of organisms (Kumar et al. 2022). Using the GTDB Taxonomy Browser resource, we created a list of the 20 000 most common bacteria taxons. A Newick file representing a tree of divergence times for these bacteria taxons was then generated using the TimeTree database. Importantly, taxa for which divergence time was estimated with considerable uncertainty or for which discrepancies arise from individual time estimations were labeled according to their “adjusted divergence time” rather than their “median divergence time.” This process was repeated until the complete timetree of GTDB common species was built.

Theoretically, potential MGEs should exhibit a distinct pattern of divergence times relative to their reference species, as we illustrated in Fig. 2. One might be tempted to try to characterize this pattern by considering evolutionary fitness and likelihood of mutations: some MGEs, such as those for antibiotic resistance genes, increase fitness relative to random segments of DNA, and certainly a large fraction of MGEs have coding regions within them.

Figure 2. Mobile genetic elements and conserved elements have distinct evolutionary divergence patterns. Three separate genome fragments of Acinetobacter baumannii strain 81 (DADYEO010000104.1) were run in skandiver, and their phylogenetic mappings visualized: each point is a distinct match to a representative species of known divergence time from A.baumannii. The nature of these genome fragments was verified using BLAST search and were characteristic of (a) a random section of A.baumannii genome not known to be mobile or conserved, (b) a conserved region of A.baumannii genome encoding ilvD dihydroxy-acid dehydratase protein, and (c) an A.baumannii plasmid X4-300 (CP064203.1).

However, there exists a much simpler characterization that we can perform based on recency of horizontal movement: a sequence that is shared by two species with high divergence times either must have remained largely unchanged over millions of years of evolution, or have horizontally moved more recently, which is the very definition of an intercellular MGE. This phenomenon should be represented by matches between a query species and reference species with high divergence times. One way to easily further distinguish between conserved elements and mobile elements is by whether the sequence is also present in other closely related species, since if it is not, then it likely is not a conserved element. Fortunately, this is rarely necessary, since we are operating at the genetic level—even highly conserved proteins show substantial changes at the genetic level after sufficiently many million years of evolution, so conserved element will have an average divergence time that is much lower than that for mobile elements.

This theoretical framework guided the validation of our method (Fig. 2). We ran skandiver on five whole-genome assemblies of A.baumannii to determine if MGEs can be distinguished by their pattern of divergence times. Divergence times were manually extracted from skandiver results and visualized using GraphPad Prism 10.

We then manually extracted a few representative genome fragments, showcasing the divergence pattern of random sequence, a conserved element, and a mobile element. Characteristics of query sequences were verified using BLAST (Altschul et al. 1990) search and plasmid/viral databases IMG/PR (Camargo et al. 2024) and IMG/VR (Paez-Espino et al. 2016). As shown in Fig. 2, where 2A represents a random genome fragment, 2B represents a fragment corresponding to conserved protein ilvD dihydroxy-acid dehydratase (verified via BLAST search), and 2C represents a fragment corresponding to the MGE plasmid X4-300 (CP064203.1), three distinct patterns of divergence can be observed.

In Fig. 2a, the random fragment simply matches to its original species with divergence time of 0, which is not indicative of a conserved or mobile element.

In Fig. 2b, the conserved fragment matches to many species within the same Acinetobacter clade, with divergence times ranging from 0 to 5 MYA, indicating of a conserved element residing within a single bacterial clade due to shared ancestry.

In Fig. 2c, the MGE matches to both many species within the Acinetobacter clade and also a wide range of diverse species with divergence times ranging from 1000 to 3000 MYA. The broader range of divergence times supports the notion of an MGE capable of horizontal movement across distinct and diverse bacterial clades over an extended evolutionary timeline.

These distinctive evolutionary patterns within genomic elements are what skandiver uses to differentiate between conserved and MGEs within an assembly.

2.5 Benchmarking details

We compared skandiver against MobileElementFinder (Johansson et al. 2021) and geNomad (Camargo et al. 2023)—although skandiver is not fully comparable to these other methods because skandiver does not rely on a curated database, these seemed like the closest possible comparisons. Default detection settings were used for both MobileElementFinder and geNomad. In skandiver, skani search/dist was run using c = 150, m = 1000, t = 10. The genome fragmentation length was set to 10 000 bp, but can be adjusted depending on the expected mobile element size. Benchmarks were evaluated on a Dell Poweredge 660xs 1U server with two 16-core Xeon 5416S processors and 512 GB memory using 10 threads.

3 Results and discussion

3.1 skandiver, MobileElementFinder, and geNomad find different putative mobile elements

We analyzed seven different contig-level whole-genome assemblies and highlighted regions identified as putative mobile elements (Fig. 3). skandiver and MobileElementFinder found comparable amounts of putative mobile elements—207 fragments of 1483.8 kb total for skandiver and 140 fragments of 1402.1 kb total for MobileElementFinder. geNomad found somewhat fewer—19 fragments of 625.7 kb nucleotides total, though it should be added that geNomad is trained to predict a subset of the various types of mobile elements. However, more interesting, with one major exception, the putative mobile elements found are very different, almost disjoint in some cases, with very little to no overlap between methods (Supplementary Fig. S1). The only notable exception was in Fig. 3d, where the entire contig we analyzed was itself a plasmid—geNomad marked the entire plasmid as an MGE, whereas skandiver and to a lesser extent MobileElementFinder marked overlapping chunks hits. Thus, it seems that the three different approaches represented here—divergence for skandiver, reference database alignment for MobileElementFinder, and machine learning on a database for geNomad—are orthogonal to each other, despite having similar aims. From these discordant results, we conclude that none of the methods is fully capturing the wide breadth of MGEs in nature.

Figure 3. Comparison of mobile element finding profiles of skandiver, MobileElementFinder (MEFinder), and geNomad: the three tools find different sets of putative genetic mobile elements. Eight randomly selected contig-level whole-genome assemblies were run through skandiver, MEFinder, and geNomad. For each whole-genome assembly, the regions corresponding to the potential mobile elements found by each method are labeled; (a) Pseudomonas aeruginosa strain PA2818 chromosome, (b) Stenotrophomonas maltophilia strain ISMMS3, (c) Streptomyces chartreusis strain ATCC-14922 chromosome, (d) Enterobacter hormaechei strain 14102KSN plasmid, stripes signify overlap between the methods, (e) Streptococcus anginosus strain VS152 chromosome, (f) Escherichia coli strain 19SZHZ663Rt chromosome, (g) Staphylococcus pseudintermedius strain NCTC5661. (h) The amount of nonmerged regions deemed a potential mobile element by each method. (i) The amount of nucleotides deemed part of a potential mobile element by each method.

3.2 skandiver can find novel MGEs that other methods miss

For further validation of skandiver’s hits, we do a deep dive into Fig. 3a in Fig. 4. We chose this example because skandiver only had three substantial hits, so we can do an exhaustive analysis without any cherry picking. Analyzed here were a set of MAGs from Pseudomonas aeruginosa. While there was very little overlap between the three methods, skandiver and MobileElementFinder were both able to find a potential MGE at around 1.6 Mb into a whole-genome assembly of P.aeruginosa (CP129688.1) (Fig. 4a and d). Upon performing BLAST nucleotide search of this fragment against the NCBI Database (Wheeler et al. 2007), it was revealed that this fragment shared 100% query coverage and 99.97% identity with the P.aeruginosa plasmid 2017–45-85 (CP109756.1) Similarly, there was another potential MGE at around 6.7 Mb that was also found by both MobileElementFinder and skandiver (Fig. 4c and f).

Figure 4. skandiver finding Pseudomonas aeruginosa plasmid p12CC3_GES-5 previously elusive to geNomad and MEFinder within P.aeruginosa strain PA2818 chromosome assembly. (a), (b), and (c) correspond to regions of CP129688.1 designated as potential mobile elements by skandiver. (d) Visualization of species matches for fragment 1 600 000–1 609 999 found by skandiver (overlaps with fragment 1 612 765–1 614 092 found by MEFinder). (e) Visualization of species matches for fragment 5 270 000–5 279 999 found by skandiver (does not overlap with MEFinder or geNomad). (f) Visualization of species matches for fragment 6 770 000–6 870 000 found by skandiver (overlaps with fragment 6 787 235–6 788 562 found by MEFinder).

More interestingly, skandiver was able to find a novel MGE that was previously uncharacterized by both MobileElementFinder and geNomad at approximately 5.27 Mb into the same whole-genome assembly (Fig. 4b and e). We were able to visualize this potential mobile element’s sharedness between species by plotting the adjusted divergence time in MYA against the different species the fragment mapped to following skani search. As seen in Fig. 4e, this fragment of interest as determined by skandiver mapped strongly to both Pseudomonas dentrificans and Pseudomonas yangonensis with negligible divergence times, but also to Azospira sp. and Steroidobacter dentrificans with significant divergence times. We verified that the fragment has elements of mobility by performing a BLAST Search against the NCBI Database, which revealed that this particular fragment shared 98% query coverage and 94.071% identity with the Raoultella planticola plasmid p12CC3_GES-5 DNA (LC735983.1). Interestingly, this element was not selected as a mobile element by MobileElementFinder due to it being an uncharacterized plasmid within MobileElementFinder’s database of mobile elements. Additionally, this element was not selected as a mobile element by geNomad, possibly due to containing an uncharacterized nucleotide makeup pattern within the plasmid that was unspecific to any chromosome, plasmid, or viral marker. Thus, skandiver can efficiently identify novel mobile elements within whole-genome assemblies without using gene annotation or markers.

To demonstrate the speed and scalability of skandiver, we further attempted to verify whether the p12CC3_GES-5 DNA plasmid discovered within a P.aeruginosa whole-genome assembly was present in any other assemblies or strains of P.aeruginosa. We downloaded a set of 200 P.aeruginosa complete genome assemblies from NCBI Assembly, and used skandiver to efficiently analyze all the assembly fragments. skandiver took about 20 min to process all 200 assemblies, with over 110 000 genome fragments. The assembly fragments were queried against the p12CC3_GES-5 DNA plasmid using the baseline skani algorithm, which processed 138 800 P.aeruginosa query sequences in approximately 0.19 s. From the results, it was determined that four unique P.aeruginosa strains in the 200 assemblies queried (JAPEVK010000001.1, JAOVYS010000001.1, MPBS01000001.1, LOJK01000001.1) exhibited considerable matches of 95% ANI or higher against the reference plasmid. skandiver can not only identify novel MGEs but also rapidly search large datasets of MAGs for any strains that may share an MGE.

3.3 Runtime, memory, and recall metrics

skandiver can efficiently query 100 separate whole-genome assemblies against the GTDB database of representative genomes (>65 000 bacteria species) in approximately 9 min (Fig. 5). skandiver is about 5 times faster when evaluating 100 whole-genome assemblies simultaneously compared to MobileElementFinder and geNomad. While MobileElementFinder and geNomad use significantly less memory than skandiver for smaller sets of assemblies, this gap reduces for 100 whole genomes, as skandiver is better at scaling up, with a consistent memory footprint even as the scale of the analysis increases. This demonstrates skandiver’s scalability, with consistent speed and memory usage that allows researchers to analyze extensive datasets with speed and precision.

Figure 5. skandiver scales well to larger datasets. Comparison of skandiver against MEFinder and geNomad in terms of time (left column) and memory usage (right column) required for completing a query. skandiver, MEFinder, and geNomad were run on datasets consisting of (a) 1 whole-genome assembly and (b) 100 whole-genome assemblies.

For completeness, we also extracted a list of large known plasmids (>10 kb) from NCBI (accession IDs available on Github repository), and then ran these sequences for detection of MGEs (Fig. 6). We considered skandiver to have correctly detected an MGE if one of its substring chunks was detected as mobile. Plasmids were either artificially inserted into a hypothetical genome of Escherichia coli (CP008805.1), or directly queried as independent isolates in the case of Fig. 6a and b, respectively. Due to the extremely large size of the hypothetical genome (>1 Gb), the genome was broken down into fragments of 5 million bp in size prior to analysis by geNomad. Notably, many of these plasmids are likely in the curated databases that geNomad and MobileElementFinder use for training/alignment, so it is not surprising that they have higher recalls in the case of isolated plasmids. However, skandiver still performs reasonably, especially for the largest plasmids, despite not using a curated database or a set of marker protein annotations.

Figure 6. Recall of mobile element finding software on known plasmid database. Comparison of recall across skandiver, geNomad, and MEFinder on large plasmids (>10 kb) from NCBI. Recall was assessed for either (a) plasmids embedded within a genome of Escherichia coli (CP008805.1) or (b) isolated. For embedded plasmids, a database of large plasmids from NCBI was inserted into a genome (CP008805.1 Escherichia coli O157: H7 str. SS17, complete genome). Plasmids were placed 10 kb apart from each other, creating a large hybrid genome containing all plasmids embedded within a genome, thus representing a scenario where MGEs are integrated. This also ensured that the plasmid start location was not aligned with the 10-kb chunks used by skandiver. A plasmid was considered found by a method if that method labeled a region of the hybrid genome to be an MGE, and that region contained one of the inserted plasmids.

In particular, skandiver exhibits significantly higher recalls when identifying plasmids embedded within a genome, which we believe is due to the other software using different detection methods for isolates versus embedded sequences (Fig. 6a). For metagenomic analysis, this may prove particularly useful when querying entire assembled genomes or contigs for potential mobile regions. In addition to this, skandiver’s accuracy, while lower than MobileElementFinder and geNomad for isolated plasmids smaller than 1 million bp, demonstrates increasing accuracy as MGE size increases (Fig. 6b). skandiver plateaus around near-100% accuracy for plasmids around 3 million bp in size, while the other tools plateau around 85% accuracy. This suggests that skandiver may be particularly proficient in identifying large MGEs compared to existing software.

To ensure skandiver reliably identifies mobile elements while minimizing false positives, we benchmarked skandiver against a list we curated of over 2000 conserved elements obtained from NCBI (accession IDs available on Github repository). Our results showed that skandiver had an overall false positive rate (FPR) of 2.2%, while geNomad had a lower FPR of 1.3%, and MobileElementFinder achieved an FPR of 0%. When filtering our conserved element dataset to only include elements above 10 kb in size, none of the methods produced false positives (Table 1). As such, skandiver’s ability to avoid false positives is comparable to existing mobile element detection software, especially in the >10 kb regime this manuscript works in.

Table 1. skandiver, MobileElementFinder, and geNomad correctly do not identify known conserved elements as mobile.

Protein ID	Species	Length (bp)	Genome position	Function	skandiver	geNomad	MEFinder	
UFZ14059.1	Streptomyces sp.	13 718	18 513	PKS I	N	N	N	
UFZ14060.1	Streptomyces sp.	11 240	32 331	PKS I	N	N	N	
ABW96540.1	Streptomyces spiroventicillatus	32 483	9369	Polyketide synthase	N	N	N	
ABW96541.1	Streptomyces spiroventicillatus	17 132	26 694	Polyketide synthase	N	N	N	
ABW96542.1	Streptomyces spiroventicillatus	17 312	43 846	Polyketide synthase	N	N	N	
Comparison of >10-kb conserved elements analyzed by different mobile element detection software. We curated a list of over 2000 conserved genes from NCBI was used. On the full list (including conserved genes <10 kb), skandiver had an overall false positive rate (FPR) of 2.2%, geNomad 1.3%, and MobileElementFinder 0%. When we filtered down to only the five conserved genes above 10 kb in size, no false positives were found.

4 Conclusion

In this manuscript, we argued in favor of an evolutionary-divergence based characterization of mobile-genetic elements, which bypasses some of the limitations of existing bioinformatics tools like geNomad and MobileElementFinder, which rely on reference databases (either through training or alignment, respectively). Our new tool, skandiver, is a scalable and versatile approach for the identification of potential MGEs within whole-genome assemblies that can detect novel MGEs. We have shown that skandiver is efficient enough to be run on hundreds of genomes and provides an orthogonal way of finding large putative MGEs that does not require annotation or training data. skandiver excels in scenarios involving high ANI values (>90%) and comparisons against diverse genomic datasets, making it especially suited for finding the distinct evolutionary patterns reflective of mobile elements.

However, we also showed that there is very little concordance on real data between methods for finding MGEs. All three methods are indeed finding MGEs, but none of the methods has particularly high sensitivity. Given that all of these methods have reasonable runtimes, it would not be an unexpected workflow for a practitioner to use geNomad or MobileElementFinder to find known MGEs, and then use skandiver to try to find some of the remaining uncharacterized MGEs. Ultimately, despite all the progress made in recent years, it seems that our ability as a bioinformatics community to comprehensively characterize MGEs at scale is still incomplete, and there exists substantial scope for further scientific exploration in this domain.

Future directions for skandiver in particular could include refining parameters and exploring methodologies to allow skandiver to distinguish different classes of (especially smaller) MGEs (insertion sequences, transposases, prophages, integrative conjugative elements, etc.). Furthermore, skandiver’s ability to efficiently requery large datasets for previously discovered mobile elements may allow it to construct a profile of the preferred mobile target genomes of mobile elements. Another way to leverage efficiently querying large datasets is to use a pangenome of strain-specific genomes, rather than just the GTDB. This may allow for increased sensitivity. One of the difficulties is that the TimeTree of Life may not include evolutionary divergence times between those genomes, but it may be possible to either compute those directly or use ANI as a rough proxy More generally though, we hope that future MGE detectors will incorporate our divergence-based metric with alignment and machine learning models to build a more comprehensive understanding of MGEs.

Supplementary Material

btae398_Supplementary_Data

Supplementary data

Supplementary data are available at Bioinformatics online.

Conflict of interest

None declared.

Funding

This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) [grant number RGPIN-2022-03074]; as well as startup funding from Carnegie Mellon University. J.S. was supported by an NSERC CGS-D scholarship. This paper was published as part of a supplement financially supported by ECCB2024.

Data availability

The data and code underlying this article are available in: https://github.com/YoukaiFromAccounting/skandiver.
==== Refs
References

Altschul SF , GishW, MillerW  et al  Basic local alignment search tool. J Mol Biol  1990;215 :403–10.2231712
Berger B , YuYW.  Navigating bottlenecks and trade-offs in genomic data analysis. Nat Rev Genet  2023;24 :235–50.36476810
Camargo AP , CallL, RouxS  et al  IMG/PR: a database of plasmids from genomes and metagenomes with rich annotations and metadata. Nucleic Acids Res  2024;52 :D164–73.37930866
Camargo AP , RouxS, SchulzF  et al  Identification of mobile genetic elements with geNomad. Nat Biotechnol  2023:1–10.36653493
Cock PJA , AntaoT, ChangJT  et al  Biopython: freely available python tools for computational molecular biology and bioinformatics. Bioinformatics  2009;25 :1422–3.19304878
Davies J , DaviesD.  Origins and evolution of antibiotic resistance. Microbiol Mol Biol Rev  2010;74 :417–33.20805405
Doolittle RF , FengD-F, TsangS  et al  Determining divergence times of the major kingdoms of living organisms with a protein clock. Science  1996;271 :470–7.8560259
Durrant MG , LiMM, SiranosianBA  et al  A bioinformatic analysis of integrative mobile genetic elements highlights their role in bacterial adaptation. Cell Host Microbe  2020;27 :140–53.e9.31862382
Forster SC , LiuJ, KumarN  et al  Strain-level characterization of broad host range mobile genetic elements transferring antibiotic resistance from the human microbiome. Nat Commun  2022;13 :1445.35301310
Frost LS , LeplaeR, SummersAO  et al  Mobile genetic elements: the agents of open source evolution. Nat Rev Microbiol  2005;3 :722–32.16138100
Johansson MHK , BortolaiaV, TansirichaiyaS  et al  Detection of mobile genetic elements associated with antibiotic resistance in salmonella enterica using a newly developed web tool: MobileElementFinder. J Antimicrob Chemother  2021;76 :101–9.33009809
Khedkar S , SmyshlyaevG, LetunicI  et al  Landscape of mobile genetic elements and their antibiotic resistance cargo in prokaryotic genomes. Nucleic Acids Res  2022;50 :3155–68.35323968
Kumar S , SuleskiM, CraigJM  et al  TimeTree 5: an expanded resource for species divergence times. Mol Biol Evol  2022;39 :msac174.35932227
Paez-Espino D , ChenI-MA, PalaniappanK  et al  IMG/VR: a database of cultured and uncultured DNA viruses and retroviruses. Nucleic Acids Res  2016;45 :D457–65.27799466
Parks DH , ChuvochinaM, RinkeC  et al  GTDB: an ongoing census of bacterial and archaeal diversity through a phylogenetically consistent, rank normalized and complete genome-based taxonomy. Nucleic Acids Res  2022;50 :D785–94.34520557
Partridge SR , KwongSM, FirthN  et al  Mobile genetic elements associated with antimicrobial resistance. Clin Microbiol Rev  2018;31 :10–1128.
Roux S , AdriaenssensEM, DutilhBE  et al  Minimum information about an Uncultivated Virus Genome (MIUVIG). Nat Biotechnol  2019;37 :29–37.30556814
Shapiro J, ed.  Mobile Genetic Elements. New York, NY: Academic Press Inc., 1983.
Shaw J , YuYW.  Fast and robust metagenomic sequence comparison through sparse chaining with skani. Nat Methods  2023;22 :1661–5.
Wheeler DL , BarrettT, BensonDA  et al  Database resources of the national center for biotechnology information. Nucleic Acids Res  2007;36 :D13–21.18045790
