SourceFinder: a Machine-Learning-Based Tool for Identification of Chromosomal, Plasmid, and Bacteriophage Sequences from Assemblies

In this work, we built classifiers to predict the chromosomal, plasmid, or bacteriophage sources for complete and partial sequences that are commonly found in assembled short read genomic and metagenomic data sets. To train the models, we started by building a training data set from sequences that were downloaded from the National Center for Biotechnology Information (NCBI) Reference Sequence (RefSeq) and PATRIC databases (24, 25). The sequences originated from a taxonomically diverse collection of bacterial hosts: 1,213 genera and 67 families for chromosomes, 502 genera and 220 families for plasmids, and 1,071 genera and 110 families for the bacteriophages (Fig. S1).

Because some taxonomic groups were overrepresented, e.g., sequences from the Enterobacteriaceae amount to nearly 20% of the available data, we chose to down sample the collection using a k-mer-based clustering approach. The number of plasmid and bacteriophage clusters converged around the 87% k-mer similarity threshold where the three sequence classes reached approximately the equal number of clusters. Therefore, we chose to use this convergence point at the 87% k-mer similarity threshold which reduced the total number of sequences by 49%, 29%, and 31% for chromosomes, plasmids, and bacteriophages, respectively (Fig. S2). Once the similar sequences were clustered, one sequence was kept per cluster (Table S1), and the remaining sequences were excluded from the full data set (Table S2). This resulted in 9,163 chromosome, 7,750 plasmid, and 8,877 bacteriophage sequences. This homology reduction preserved nearly all of the taxonomic diversity in the data set (Fig. S1). Before clustering, the average lengths of the chromosomes were approximately 5,000,000 bases (5 mb) while the plasmids and bacteriophages were 50,000 bases (50 kb); these average lengths remained the same after the clustering (Fig. S3).

Random forest classifiers were trained using a single representative sequence from each cluster. The features were oligonucleotide k-mers that were generated by stepping by one nucleotide at a time over each chromosome, plasmid, and bacteriophage sequences (Fig. 1). The k-mer frequencies were calculated and collected into a matrix where the rows were the sequences and the columns were the normalized k-mer counts. K-mers of different lengths may provide different information to the models. Short k-mers are informative for the sequence structure, codon usage, and GC content, while longer k-mers may be valuable for unique patterns such as the presence or absence of genes or mutations. Therefore, the k-mer length was optimized to improve the model performance. Sequences were classified using k-mers from five to eight nucleotides. Overall, we found that the models built from whole sequences with different k-mer sizes did perform similarly. The highest performance was obtained with the 6-mers where the hold-out performance reached 0.992 ± 0.001 area under the curve (AUC) (0.928 ± 0.002 Matthews correlation coefficient [MCC]). Per sequence class, the model accuracy reached 0.979, 0.944, and 0.930 for chromosomes, plasmids, and bacteriophages, respectively (Table 1; Table S5).

To evaluate the accuracy and generalizability of the models for classifying incomplete sequences (a likely scenario when analyzing metagenomic assemblies), the models that were trained with the whole sequences were tested against randomly selected fragments of chromosomes, plasmids, and bacteriophages that were 5 kb in length. K-mers from five to eight were tested for the fragment classification. The highest average prediction performances reached 0.876 ± 0.001 AUC (0.419 ± 0.006 MCC) with 5-mers. Per sequence class, the model accuracy dropped to 0.018, 0.839, and 0.893 for chromosomes, plasmids, and bacteriophages, respectively. Even though the overall model performance was high with the fragments, the robustness of the model decreased compared with the evaluation made with the whole sequences. In particular, the model failed to distinguish the chromosomes from plasmids and bacteriophages (Fig. 2; Table 1; Table S5). This could suggest that the original contig lengths of the chromosomes were influencing the prediction, and would result in a negative outcome for predicting against smaller sequence fragments.

To improve the model generalizability on sequence fragments, the random forest models were regenerated by training and testing using randomly sampled 5-kb fragments from full-length chromosome, plasmid, and bacteriophage sequences with various k-mer sizes. When the model was trained and tested using a single random 5-kb fragment sampled from each sequence, the highest average model performance reached 0.911 ± 0.000 AUC (0.655 ± 0.005 MCC) with 5-mers (Table 1; Table S5). As the number of randomly sampled fragments increased from one to 10, the prediction performances increased to 0.939 ± 0.000 AUC (0.713 ± 0.002 MCC) with 5-mers. Per sequence class, the model accuracy reached 0.895, 0.652, and 0.844 for chromosomes, plasmids, and bacteriophages, respectively (Fig. 3; Table 1; Table S5). Although training with sequence fragments came with a notable improvement in the accuracy for the chromosome class, it also corresponded with a loss in accuracy for the plasmid class compared with the models trained from whole sequences.

In this analysis, there were two chromosomes, 911 plasmids, and 51 bacteriophage sequences that were shorter than 5 kb in the test and hold-out data sets. The three-class classifier trained with 5-mers identified the source of the <5-kb fragments with 1.0, 0.731, and 0.810 accuracy for chromosomes, plasmids, and bacteriophages, respectively. These accuracies for these shorter fragments match the overall accuracies obtained from the hold-out set (Fig. 3; Table 1; Table S5).

The taxon representation in the training data set was unbalanced. While we observed high overall performances, there is a risk that the model is only predictive for the well represented taxa. In order to assess the generalizability of the best fragment-based model (defined as the 5-mer model built from 10 randomly subsampled, 5-kb fragments per element) for the underrepresented taxa, performances were inspected for each host taxon at the family level (Fig. S4; Table S6). In addition, Spearman’s correlation coefficient was calculated between the hold-out performances of every host taxon and the corresponding sequence counts in the training data. A correlation was detected in case of plasmids (Spearman’s correlation coefficient = 0.403, P value < 0.05) and bacteriophages (Spearman’s correlation coefficient = 0.480, P value < 0.05). However, correlation was not significant for chromosomes (Spearman’s correlation coefficient = −0.02, P value > 0.05). These results suggest that the host representations in the training data have an impact on the accuracy of the plasmid and bacteriophage identification. Additionally, this correlation could be the consequence of the overlap between the sequence structures of plasmids and bacteriophages, e.g., GC content; therefore, more examples of plasmids and bacteriophages are required for the correct sequence classification (Fig. S5).

The presence of prophages, which are chromosomally integrated bacteriophage sequences, may cause chromosomes to be misclassified. Prophages might be a major challenge in the fragment-based model because the model relies on a small fraction of the sequences for the classification. Of the hold-out data set, 410 chromosome and 62 plasmid sequences were predicted to contain prophages using PhiSpy (26). When the whole-sequence-based model was tested with the prophage containing hold-out sequences, 1% of the chromosomes and 4.5% of the plasmids were predicted to be bacteriophages (Table S7). When the whole-sequence-based model was tested with the hold-out sequences lacking prophages, 2% of the chromosomes and 3% of the plasmids were predicted to be bacteriophages (Table S7). Overall, for the whole-sequence-based model, the presence of bacteriophage elements appears to have a negligible effect on the overall accuracy.

When the same experiment was repeated with using the model that was generated from 5-kb fragments (10 random subsequences, 5-mer features), 50 chromosome fragments and 90 plasmid fragments were predicted to contain at least 500 bp of prophage DNA. When the fragment-based model was tested with the prophage containing fragments, 20% of the chromosomal fragments and 6% of the plasmid fragments were predicted to be bacteriophages (Table S7). When the fragment-based model was tested on fragments lacking prophage elements, 5% of the chromosomal fragments and 8.1% of the plasmid fragments were classified as bacteriophages (Table S7). The presence of prophages was found to have a significant effect on the misclassification of the 5-kb chromosome fragments (P value < 0.05, Chi-square contingency test), but not for the 5-kb plasmid fragments (P value = 0.108, Chi-square contingency test). Overall, these results indicate that in the models built from short 5-kb sequences, some of the misclassification of chromosomal sequences is due to the presence of prophage elements.

To test the fragment model’s (10 random subsequences, 5-mer features) generalizability to metagenomic scaffolds, simulated metagenomic reads were generated using BBMap from the complete chromosome, plasmid, and bacteriophage sequences. The simulated reads were assembled into scaffolds where the majority of the sequences were fragmented into the multiple contigs (N50 = 40,936 bp). A total of 27,634 scaffolds (corresponding to 29% of all scaffolds) were at least 1 kb. These were used for the validation of the fragment-based model. After removing chimeric sequences, the remaining 26,484 scaffolds (chromosome: 25,888, plasmid: 397, bacteriophage: 199) were tested by the fragment model which reached the accuracy of 0.914, 0.698, and 0.794 for chromosomes, plasmids, and bacteriophages, respectively (Fig. S6). These performances are comparable with the performances for the partial sequences. Overall, these results suggest that the fragment model can classify the metagenomic scaffolds minimum 1-kb length with at least 70% accuracy which corresponds to the plasmid identification performance.

We tested our three-class classifiers (whole-sequence-based and fragment-based) against the PPR-Meta tool. First, the PPR-Meta tool was tested with the whole sequences from the hold-out data set without fragmenting the sequences. The tool predicted the chromosome, plasmid, and bacteriophage sequences with 0.868, 0.819, and 0.986 accuracy, respectively (Table S8; Fig. S7). Compared with our whole model (5-mer features), PPR-Meta performed approximately 8.3% and 11.9% poorer with the chromosome and plasmid sequences, and 6.3% better with the bacteriophage sequences. Additionally, the PPR-Meta tool was tested using the hold-out data set that was fragmented into 5-kb fragments and sampled 10 times. The tool predicted chromosome, plasmid, and bacteriophage sequences with 0.780, 0.749, and 0.977 accuracy, respectively (Table S8; Fig. S7). Compared with our fragment model (10 random subsequences, 5-mer features), PPR-Meta performed approximately 12.8% poorer with the chromosome sequences, but 14.9% and 15.8% better with the plasmid and bacteriophage sequences, respectively. Finally, the source of the simulated metagenomic scaffolds was detected using the PPR-Meta tool. The chromosome, plasmid, and bacteriophage derived contigs were identified with 0.820, 0.877, and 0.935 accuracy, respectively (Table S8; Fig. S7). Likewise, PPR-Meta performed approximately 10.3% poorer with the chromosome sequences, but 25.6% and 17.8% better with the plasmid and bacteriophage sequences, respectively. These results suggest that PPR-Meta is highly precise for the bacteriophage identification, yet our whole and fragment models are more accurate at identifying the chromosome sourced sequences. In terms of total runtime, PPR-Meta’s runtime increased with the contig length (Table S8) in contrast to our tool, as more models were used for large contigs. Overall, our tool could be run efficiently and as accurately as PPR-Meta regardless of contig length.

The random forest classifier that was trained with 5-mers and 5-kb fragments that were subsampled 10 rounds from each element is implemented as a command line tool well as a web service called SourceFinder. The command line version of SourceFinder is available through Bitbucket (https://bitbucket.org/vlagri/sourcefinder/src/master/). The command line tool requires Anaconda packages (version: 4.4.0) and KMC (version: 3.0.0) and can run on a Linux machine with 1 CPU, 10 GB RAM memory (for 1,000 input elements), and 22 GB storage (for downloading trained random forest models). It makes predictions for each contig based on the average class probabilities of five different random forest models. The web service is available through the Center for Genomic Epidemiology (https://cge.food.dtu.dk/services/SourceFinder/). It has a user option for selecting the number of times a contig is sampled. The web service runs on backend hardware and thus, besides the web browser, has no other computing requirements. Runtimes are generally 1 minute per genome. An example job submission and output are shown for the online SourceFinder tool in Fig. 4.

A set of 18,022 complete representative bacterial chromosomal sequences were downloaded from the National Center for Biotechnology Information (NCBI) Reference Sequence (RefSeq) database; 10,852 plasmid and 12,930 bacteriophage sequences were downloaded from the Pathosystems Resource Integration Center (PATRIC) databases in August 2020 (30, 31). The taxonomic information was extracted from the sequence headers. The extracted genus level information was translated to the family level using the NCBI Taxonomy information by the Python ete2 package (version:2.3.10) (32, 33). To remove the redundant sequences while retaining the diversity of the data, all sequences were clustered based on the k-mer similarities using the KMA tool (version: 1.1.0) for a given similarity threshold (34). KMA applies the Hobohm-1 algorithm to cluster the similar sequences, and it was performed per type of data with the following parameters: -k16, -NI, -Sparse and for the following query and template k-mer similarity thresholds: 81%, 84%, 87%, 90%, 93%, 95%, 96%, and 99%. For the machine-learning training, testing, and validation, one representative sequence was kept from each cluster for a given k-mer similarity threshold, and the other sequences were excluded.

In order to predict the three sequence classes (chromosome, plasmid, and bacteriophage), we generated multiclass machine-learning models. These models were trained on matrices containing k-mer frequencies from whole or fragmented sequences (Fig. 1). All of the input sequences were divided into nucleotide k-mers and the frequencies of the k-mers were counted using the KMC tool (version: 3.0.0) (35). KMC considers each k-mer and its complement, counting only the k-mer that is highest in a lexicographical sorting. KMC was performed for k-mer sizes 5, 6, 7, and 8 with the following parameters: “-fm” and “-ci1.” These optional parameters describe the type of input as a multi-fasta file and the minimum k-mer frequency of 1 that should be reported, respectively.

For each k-mer length, frequencies were collected into a matrix where the rows represented the isolate inputs and the columns represented unique k-mers (Fig. 1). Each matrix contained all three classes of sequences: chromosome, plasmid, and bacteriophage. To eliminate the bias introduced by differences in sequence lengths, the k-mer frequencies were normalized by dividing the counts with the total number of used k-mers of the given sequence. Each normalized matrix was utilized for training, testing, and validation of the machine-learning model.

Random forest was chosen for the machine learning due to its ease of interpretability, ability to handle multi-class problems, and robustness against overfitting (36). Models were generated using ensemble.RandomForestClassifier from the Python Scikit-learn library (version: 0.23.2) (37). To assess the model performances and overfitting, 10% of the input data were held out for validation, and the remaining was used for training and testing the models using a 5-fold cross-validation, which was applied through model_selection.GroupShuffleSplit from the Python Scikit-learn library (version: 0.23.2) (37).

Selected random forest parameters were tuned using the random grid search method through model_selection.GridSearchCV from the Scikit-learn library (Table S3) on the training and testing data (version: 0.23.2) (37). The grid search was completed only for the normalized whole sequences and all k-lengths, repeated 100 times. The optimal random forest parameters were: tree depth of 400, 800 trees, and data bootstrap as false. The remaining parameters were kept as default. Unless otherwise stated, these parameters optimal for the whole sequences were used for all random forest models described in the study. Each tuned random forest model was trained and tested five times and the ensemble of the trained models was applied to the hold-out data set for the final validation. Each classifier provides class probabilities that sum to 1.0. For each tested sequence, the category with the highest class probability was accepted as the final prediction.

Machine-learning model performances were assessed using multiple performance metrics: the receiver operating characteristic (ROC) area under the curve (AUC), Matthews correlation coefficient (MCC), macro F-score, macro precision, macro recall, accuracy, and confusion matrix using Scikit-learn (version: 0.23.2) (37). AUC was calculated using the one versus rest option to handle the multi-class outputs and the performances were macro-averaged. Statistical tests, including Spearman’s correlation coefficient and the chi-square contingency test, were calculated using stats.spearmanr stats.chi2_contingency functions from SciPy (version: 1.2.2) (38). Additionally, the significance threshold for P value was set to 0.05.

In order to avoid potential biases due to sequence length, models were built from subsampled sequence fragments. Fragmentation was performed by randomly sampling complete sequences using the subsequence length of 5,000 bases (5 kb). When sequences were shorter than 5 kb, the entire sequence was selected. Sequences larger than 5 kb were assigned a random starting point between the first position and 5 kb before the last position. We also tested sampling multiple sequence fragments from a single sequence. In theory, multiple fragmentation could improve the model performances by increasing the coverage of the sequences in the training data. To assess the impact of multiple fragmentation on the model performances, sampling rounds of one, five, or 10 were tested. In other words, if the sampling round was five, five different 5-kb fragments were sampled from each given sequence. When a sequence was shorter than 5 kb, the entire sequence was used multiple times. Additionally, to test the robustness of the method, the fragmentation process was repeated five times to confirm the model performances with different sets of random fragments.

The trained models were named after the type of input that was used in model training. In the “whole”-sequence models, the model was trained and tested with the matrix composed of k-mer frequencies of whole sequences where each sequence was considered a row in the matrix. In the sequence “fragment” models, the model was trained and tested with the matrix composed of k-mer frequencies of fragments where each fragment was considered a row in the matrix. Furthermore, if the fragment sampling round was greater than one (five or 10), fragments originating from the same sequence were kept in the same training, testing, or hold-out data set to prevent possible overfitting.

In order to assess the generalizability of the fragment model against simulated metagenomic reads, 100 sequences per class were randomly selected from the hold-out data set (Table S4). In total, metagenomic reads were simulated from 300 full-length sequences with random coverage (from 10× to 50×) to simulate coverage similar to metagenomic sequencing using the randomreads.sh tool in the BBMap tool suite (version: 38.90). The following parameters were used: the read length was set to 150 bp, minimum insertion size to 300 bp, maximum insertion size to 600 bp, and phred quality score to 35; also, the tool was used with the metagenome and paired end options (39). The simulated metagenomic reads were concatenated into a single file and assembled into metagenomic contigs using the spades.py tool from the SPAdes software (version: 3.15.2) with the –meta, −12, and –only-assembler options (40). The metagenomic scaffolds shorter than 1 kb were discarded. The original source of scaffolds was determined by realigning the assemblies to the 300 reference sequences using the nucmer function of the Mummer tool (version: 4.0.0beta2) with the default parameters (41). The metagenomic scaffolds that aligned to multiple sequence classes were discarded due to likely chimerism.

In order to understand bacteriophage predictions occurring within chromosomal and plasmid sequences, prophages were identified from the chromosome and plasmid sequences using the PhiSpy tool (version: 4.2.17) (26). As required by the PhiSpy tool prior to prophage detection, sequences were annotated with genomic features in the GenBank format (.gbk) using Prokka (version: 1.14.5) (42) with the default parameters. Afterwards, the PhiSpy tool was executed with the annotated files by setting the minimum contig size parameter to 1 kb.

PPR-Meta uses a bi-path convolutional neural network and is the state-of-the-art three-class classifier for identifying chromosome, plasmid, and bacteriophage derived sequences at once. The PPR-Meta tool (version: 1.1) was tested with the hold-out data set, including 908 chromosome, 758 plasmid, and 913 bacteriophage sequences using the default parameters. The comparison was performed for the whole sequences, fragmented sequences of 10 rounds of sampling, and simulated metagenomic scaffolds. The tool performances were reported using the accuracy metric.

The Python 3.6.10 scripts used for this project are available on Bitbucket (https://bitbucket.org/vlagri/sourcefinder/src/master/) All genomes are available through the NCBI RefSeq database (https://ftp.ncbi.nlm.nih.gov/refseq/release/) and the PATRIC resource (https://www.patricbrc.org) using the genome identifiers in Tables S1. The fragment-based model is available as a web service at https://cge.food.dtu.dk/services/SourceFinder/.