The effect of taxonomic, host-dependent features and sample bias on virus host prediction using machine learning and short sequence k-mers

Data collection

Viruses infecting host taxa Mammalia, Viridiplantae and Insecta are the most relevant groups in terms of fecal virome studies and include enough virus taxa to explore the capacity of machine learning to predict virus host. Arboviruses (arthropod-borne viruses that infect hosts both mammalian and insect hosts) were excluded from the final dataset as their genomic features are potentially shaped by different hosts. The final data set consisted of 1363 genomes from 42 virus families represented by 160 genera (Supplementary Table S1). Of these viruses, 49.5%, 24.7%, and 25.8% infect mammals, insects, and plants, respectively. Since evolutionary clade Viridiplantae is the most variable (includes species that are the most distant from each other), the taxonomy of each host that belong to this clade was investigated. As a result, 351 out of 352 plants represented flowering plants (class Magnoliopsida) and one host species belonged to cycads (class Cycadopsida, clade Acrogymnospermae).

Predicting hosts of known and novel taxa

Metagenomic sequencing can yield genomes of viruses that belong to known taxa, to novel genera and novel families. Identification and prediction of virus hosts in the latter cases can be most challenging. Therefore, three approaches to dividing the dataset into training and test samples were used as detailed above. First, the training and test set were sampled randomly with virus families and hosts equally presented in each set. Random Forest, LightGBM, XGBoost and Support Vector Classifier models were trained to predict virus host based on the frequency of nucleotide k-mers in viral genomes (Fig. 2A, Supplementary Fig. S1). Even mononucleotide frequencies were quite good predictors of a virus host (median weighted F1-scores for ten replicates ranged between 0.78 and 0.80 for all four ML models compared to 0.32 for a dummy prediction). Di-, tri- and tetra-nucleotide frequencies were all good predictors of a virus host, with median weighted F1-score values between 0.92 and 0.96. This exceeded the prediction quality of the Host Taxon Predictor used as one reference method (median weighted F1-score equaled 0.89), but was close to the results of tBLASTx (median F1-score equaled 0.93) chosen as another reference. The best median weighted F1-score was achieved by the SVC model using 4-mers. Longer k-mers worked well for the SVC model but were inferior for other ML models. Even though not superior over tBLASTx in terms of precision in a simple task of predicting a host of known viruses across huge datasets, pre-trained ML classifiers significantly speed up the calculations: BLAST algorithms perform search through the database for every query sequence to find similar sequences with known hosts, while an ML model, once trained, can predict the host of a new sequence after a simple computation of genome features. In some cases, ML classifiers provided even more accurate prediction compared to the tBLASTx method.

A more realistic and challenging task for host prediction is the case of unknown novel genera and families. When training and test datasets did not include closely related viruses, all methods produced significantly lower prediction quality (Fig. 2B). Mononucleotide frequencies were the least informative for virus host prediction, with median weighted F1-scores between 0.60 and 0.64, which were still higher than the dummy prediction weighted F1-score of 0.29. Longer k-mers yielded scores between 0.70 and 0.79 using all four models, with k = 4 and SVC model providing slightly better results (median weighted F1-score = 0.79). Reference methods, Host Taxon Predictor and tBLASTx, produced median weighted F1-scores of 0.72 and 0.68, respectively. The robustness of prediction (variation between ten training/test sample replicates) fell sharply. Typical difference between the most and the least successful replicates was about 0.20 for non-overlapping genera, compared to < 0.05 for closely related virus datasets, which allowed (and, statistically, almost guaranteed) that viruses from a genus are present in both training and test datasets (Fig. 2A). However, even in the worst training/test dataset replicates the host prediction results were significantly better than a dummy prediction. Therefore, ML could predict well a host of novel virus genera using short k-mers. Machine learning models (methods reported here and the Host Taxon Predictor) were clearly superior over tBLASTx. Although the latter provided a fair result in half of replicates, the spread of performance scores among dataset replicates was higher (0.39) compared to the other models.

Predicting hosts of novel virus families is yet a more challenging task (Fig. 2С). Mono-nucleotide frequencies were moderately predictive of the virus host, with median weighted F1-score ranging between 0.50 and 0.53 vs. the dummy prediction value of 0.32, and with a very good reproducibility between replicates and ML models. Longer k-mers yielded essentially the same median weighted F1-scores, but much higher variation between replicates, from 0.17 (i.e. twice lower than a dummy prediction) to 0.76. Host Taxon Predictor provided essentially the same prediction quality (median weighted F1-score = 0.49) as ML-based approaches suggested here. Therefore, longer k-mers can be predictive of a host of unknown families, but this is not robust and is highly dependent on training and test sample composition. tBLASTx was apparently not applicable for identifying hosts of viruses representing unknown families.

The effect of taxonomic richness and training dataset composition on prediction quality

To better understand the reasons for variations in prediction quality across different realizations of the train-test split, we examined the composition of the training dataset for non-overlapping genera strategy (Fig. 2B) in greater detail. Since the primary objective was to maintain the ratio of susceptible hosts while ensuring previously known viruses (virus genera) did not appear in the test sample, the number of taxa (virus families or genera) in the training dataset varied across split realizations. Specifically, the number of unique families and genera in ten train-test split iterations ranged from 17 to 23 and from 32 to 48 respectively. We observed that the F1-score tended to increase with greater taxonomic richness (the number of taxa), though no statistically significant correlation was found between these variables (Supplementary Fig. S2, Pearson’s r = 0.349, p-value = 0.323).

The observed variation in prediction quality may stem from the classifier’s uneven performance across different virus families. Hosts of some families were much more likely to be predicted. The highest prediction probability was observed for Picornavirales and other viruses with a ssRNA(+) genome. The top-performing model (SVC trained on nucleotide 4-mers) correctly predicted hosts in fewer than half of the cases for plant viruses of the Bromoviridae family and misclassified hosts in nearly half of the cases for viruses of the Rhabdoviridae family – a group that includes viruses infecting all three host classes (Fig. 3). It may be tempting to speculate that Rhabdoviridae, a unique virus group that infects all three host classes, did not adapt its nucleotide content to any, but the current data do not allow such conclusion. Additionally, lower accuracy was observed for mammalian viruses of the Paramyxoviridae family and plant viruses of the Secoviridae family.

Nucleotide and amino acid feature sets

Predicting hosts of novel virus genera appeared a feasible task even though the taxonomic signal was largely lacking. Incorporating a range of complementary genome features – such as nucleotide, amino acid (or functional amino acid class) k-mer frequencies – has been shown to improve accuracy in virus host prediction tasks, as demonstrated in previous studies²⁸. To further improve host prediction quality in this task, nucleotide-based features were compared to amino acid features and their combinations (Fig. 4, Supplementary Fig. S3). Combinations of NT feature sets in the range between 1 and 7 or arbitrary combinations of NT feature sets (data not shown) did not lead to any prediction quality gain. The most informative features varied across classifiers trained on the same feature set but different train-test splits. Despite variability, certain features consistently appeared in the top 10 most important across all models trained on the same feature set (Fig. 5). K-mers with a known mechanism of variation among virus taxa (CpG, a target of ZAP, UpA, a target of RNAse L, ApC/ApT, a target and a result of APOBEC editing in SARS-CoV-2) were among the important features, but not on the leading positions. Moreover, dinucleotide features were not found at all when comparing importance of all possible k-mers with k between 1 and 7.

Classifier trained using only AA k-mers showed poor performance compared to models trained on NT feature sets, although even the least successful replicates were better than the dummy prediction (weighted F1-score = 0.29). Increasing amino acid k-mer length reduced the weighted F1-score. A combination of feature sets from different sequence representations (AA and NT) did not differ significantly from the corresponding NT feature sets.

Principal component analysis (PCA) plots revealed no clear separation between viruses infecting different host classes in the space of nucleotide or amino acid k-mers frequencies (Supplementary Fig. S4); however, those belonging to families with fewer classification errors (e.g., Peribunyaviridae, Coronaviridae) exhibited more defined clustering.

In the amino acid k-mer space, data points were generally densely clustered and viruses infecting distinct host classes were indistinguishable, although a few well-defined groups could be observed. These corresponded to plant-infecting viruses of the Potyviridae family and mammalian-infecting viruses of the Picornaviridae family.

While classifiers trained on amino acid k-mers showed overall worse prediction performance, they outperformed nucleotide-based models for 25 out of 178 virus genera (Supplementary Fig. S5). Although training on amino acid features improved prediction for certain Insecta-infecting genera (e.g. Entomobirnavirus, Culicidavirus, Gambievirus, Jonvirus, Mobuvirus, Peropuvirus), the accuracy of predicting viruses infecting insects for the model remained substantially weaker than nucleotide-based predictions. F1-score 95% confidence intervals for Insecta hosts were of [0.651, 0.788] for the best nucleotide model versus [0.521, 0.630] for amino acid models.

Notably, 21 of the 25 genera that benefited from using AA k-mers have negative-sense RNA genomes, with 13 concentrated in just two families: Paramyxoviridae (6 genera) and Rhabdoviridae (7 genera). One could speculate that their genomes are packed with nucleoprotein and this less exposed to cytoplasmic factors affecting sequence content, while they are known to contain several highly conserved amino acid motifs. A modest but consistent improvement in prediction accuracy was observed for plant-infecting positive-sense RNA viruses, particularly in the Cucumovirus and Ilarvirus genera (family Bromoviridae).

Application of host prediction methods to simulated next generation sequencing (NGS) data

The task of host prediction is relevant for the analysis of raw sequencing reads and contigs. Genetic data acquired from metagenomic analysis is represented by reads, which are usually assembled into contigs, and rarely into full genomes. Metavirome reads classification implies predicting hosts not only by full genomes but also by short genome fragments. To solve this task, classifiers trained on full genomes were tested on 400 and 800 nucleotide fragments to infer hosts of viruses that belong to genera absent in the training dataset (Fig. 6A, C, Supplementary Fig. S6, S7). The best results were obtained by using nucleotide di- and tri-mers, although the prediction quality was lower than for full genomes. Median weighted F1-scores were 0.56–0.59 for 400 nt fragments and 0.60–0.63 for 800 nt fragments compared to 0.71–0.75 range for full genomes on the same datasets for RF, LGMB and XGB models. The SVC model showed less reliable and generally worse results in this particular task. Longer k-mers were much less informative for classification of short genome fragments, especially for 400 nt, yielding lower weighted F1-score and significant variation between replicates. Median weighted F1-scores obtained here were marginally better than those of the Host Taxon Predictor that uses a very similar approach, while tBLASTx demonstrated very poor host prediction using short genome fragments for viruses of previously unseen genera.

Much better host prediction for short genomic fragments could be obtained when classifiers were trained on fragments of the same length as used for testing (Figs. 6B and D and 7, Supplementary Fig. S6, S7). Dinucleotide frequencies provided a very good result, which was only somewhat improved by using nucleotide triplets. The RF models performed worse on longer k-mers in this task, in contrast to classifiers trained on full genomes. Other three models yielded median weighted F1-scores between 0.60 and 0.63 using trimers and 400 nucleotide fragments and 0.62–0.67 using 800-nucleotide fragments. This was 0.04 better than the same models trained using full genomes (Fig. 7). Longer k-mers (5–6 nt in length) were marginally worse than 2–4 nt fragments, but not as appalling as the same 5–6 nt k-mers applied for the same task of identifying 400 nt fragments using models trained on full genomes instead (Fig. 6A). Therefore, different k-mer length might be preferred for classification of long and short fragments. Most model and feature combinations were better than the HTP, which uses essentially the same features and similar algorithms, and yielded median weighted F1-scores of 0.58 and 0.61 for 400 and 800-nucleotide long fragments, respectively. Using tBLASTx search on a fragments dataset to predict hosts by virus fragments failed completely (data not shown), therefore results of tBLASTx trained on full genomes (same as in Fig. 6A, C) are presented in Fig. 6B, D.

For both fragment sizes datasets LGBM, XGB and SVC models showed very close performance for the most feature sets. RF classifiers, on the contrary, were clearly less capable of predicting virus host using fragments datasets.

With decrease in genome fragments length classifiers drop in median weighted F1-score. Apparently, shorter fragments lose signals inherent in genomes, which could be recognized by classifiers trained on full-length genomes dataset.