Characteristics of predicted network ($G_{predicted}$) and newly discovered viruses

The binary model performed well in predicting the presence of links formed due to sharing of hosts between two virus nodes in the network (mean positive predictive value=0.99, sensitivity=0.96, F-score 0.97, Fig. S6). The distribution of predicted probability for all links using the binary model showed clear bimodal distribution (Fig. S7a). The accuracy scores as a function of precision and recall indicated good model performance beyond 0.15 predicted probability for the binary model (Fig. S8). Hence, as a more conservative approach and to give more weight to the precision, we decided to use 0.7 as an optimum threshold for detecting a positive link between two nodes (viruses). The performance of the multilabel model varied for taxonomical orders, with a high to moderate performance for predicting taxonomical group and order of ‘humans’ and Cetartiodactyla (Figs. S7, S9). For 531 novel viruses, we identified 184,055 possible links to other viruses formed due to sharing of hosts (based on the optimum probability threshold of 0.7 identified for the binary model) generating the predicted network ($G_{predicted}$, Fig. 1, Fig. S7a). For these predicted links, between two viruses, the multiclass model was able to estimate the potential taxonomic order of the shared species for 175,113 links. For the remaining links, the model was not able to confidently predict a specific taxonomic order. Jaccard coefficient consistently showed high importance with all three importance metrics tested (gain, cover, and weight), indicating predictive ability in identifying missing links between unipartite viruses formed due to sharing of hosts (Fig. S10). Although, we suspect that these neighborhood-based topological features will be more predictive in a bipartite network setting due to more information flow between nodes.

Empirical biological networks are rarely scale-free (network with large hubs and showing a power-law distribution for degree)¹⁴ but a recently published study with host-host projected networks where links are represented by sharing of pathogens between hosts, has shown scale-free nature where models with power-law distributions showed the best fit for host-parasite networks¹⁵. Similarly, both observed ($G_c$) and predicted ($G_{predicted}$) networks provided evidence that some viruses shared significantly larger numbers of hosts, creating hubs of preferential attachment and showed weak evidence of scale-free nature (loglikelihood ratio test p>0.05). The degree (Fig. 2a, e, f) and betweenness centrality (Fig. 2b, i, j) distributions for predicted network ($G_{predicted}$) showed longer tails (Kolmogorov-Smirnov, p<0.05) than the observed network ($G_c$) both at network and virus family level. Mean network degree for all virus families reduced significantly with the addition of newly discovered viruses that were predicted to have fewer links than known viruses, indicating lower host plasticity for novel viruses than known viruses or insufficient adjustment of reporting bias (Fig. S11).

Based on a linear regression model with node-level permutations (10,000 permutations), our adjustment for search effort (PubMed hits) was found to have no effect on the degree (p=0.39, Fig. S12) and betweenness centrality (p=0.22, Fig. S13), but did significantly affect the eigenvector (p<0.05, Fig. S14) and clustering coefficient (p<0.05, Fig. S15) of novel viruses. These results indicate that sampling and reporting efforts affect our understanding of the predilection towards certain species as illustrated by clustering in the network, but do not affect the prediction of missing host links quantified by degree centrality within the network. Many of the newly discovered viruses were mostly detected in only one species (mean=1.32, SD±0.99, n=944). Long tails of centrality distributions generated for the predicted network ($G_{predicted}$) and comparatively lower centrality measures for novel viruses, when compared with known viruses, support a tendency for newly discovered viruses to be more host-specific than previously recognized viruses, a pattern that should be further evaluated with additional sampling effort to identify the full host range for novel viruses.

Importantly, a comparison between virus families of novel viruses showed that novel coronaviruses had a higher degree (p<0.001, Fig. 2c, Fig. S12), betweenness (p=0.02, Fig. 2d, Fig. S13), and eigenvector (p<0.001, Fig. S14) centralities in the predicted network compared to newly discovered viruses in all other virus families (Fig. 2c, d, g). In addition, the raw detection data showed significantly higher host diversity for novel coronaviruses with a mean of 2.02 (SD±2.03, n=114) unique host species (maximum of 15 species) compared to 1.22 (SD±0.70, n=834) for other novel viruses detected in this study. This finding raises concern about the ability of novel coronaviruses to infect a greater number of species than viruses from other families. The recently emerged SARS-CoV-2 and the previously emerged SARS-CoV-1, have shown a wide host breadth¹⁶. These predictions for novel coronaviruses highlight their key ecological properties that can influence spillover into humans. Following coronaviruses, novel flaviviruses showed significantly higher betweenness centrality (p<0.001). Host taxonomic order for novel viruses had no significant association with the degree centrality of the virus in the predicted network. Predicted network characteristics not only differentiate virus families based on network characteristics but also predict network characteristics that are key in understanding the ecology of a novel virus and its behavior within the network community of hosts, including the expected breadth of host species most likely to be infected by that novel virus.

Prioritizing novel viruses for further characterization

For the 531 novel viruses, we developed prioritization metrics that inform on the ecological and evolutionary tendencies for spillover based on number of human links with known viruses predicted by the multiclass model. Novel viruses from Herpesviridae, Rhabdoviridae, Coronaviridae, Adenoviridae, Astroviridae, and Paramyxoviridae families not only showed a high median probability of sharing human links with known viruses (Fig. S16) but also were predicted to have large numbers of human links in the predicted network ($G_{predicted}$). Novel viruses of the Picobirnaviridae and Rhabdoviridae families detected here have been speculated to be hyper-parasites infecting bacteria and insects and were identified in mammalian host samples. Hence the predicted associations for these virus families should not be inferred as infection but only as detection in host samples (e.g., potentially insect viruses detected in oral swab samples from bats). Based on generalized linear mixed models, search effort (PubMed hits) was not associated with the predicted number of human links (p=0.24, Table S1) nor the mean probability of sharing human links for novel viruses (p=0.778, Table S2).

For a relative comparison of zoonotic risk for novel viruses, a prioritization score was developed based on the predicted probability of links being human and the number of shared human links in the predicted network for a given virus. To understand the performance of the prioritization score, we compared scores for known zoonotic and non-zoonotic viruses generated by the ensemble of both binary and multi-class models. Results indicated significantly higher prioritization scores for known zoonotic viruses (Fig. S17, p<0.001) compared to known non-zoonotic viruses. Prioritization scores were derived essentially from the prediction of new/yet unobserved network links generated by the virus with another virus formed due to sharing of hosts. However, models were unable to predict new links for well-recognized viruses that have numerous hosts, such as Rabies virus and West Nile virus, and consequently resulted in a prioritization score of zero. Figure 3a–d shows the top ten and bottom five novel viruses from four virus families for relative comparison based on the prioritization score (Figs. S18–S24). PREDICT_CoV-15 found in two Phyllostomidae bats from South America (Artibeus lituratus, Sturnira lilium) scored the highest prioritization score in all novel viruses. Other top ten novel coronaviruses based on the prioritization score included viruses detected in Phyllostomidae bats (PREDICT_CoV-4, PREDICT_CoV-13, PREDICT_CoV-11, PREDICT_CoV-5). Out of these, PREDICT_CoV-11 was also detected in Mormoopidae species (Pteronotus personatus) and PREDICT_CoV-5 was found in Vespertilionidae species (Bauerus dubiaquercus) during the surveillance. These also included coronaviruses detected in Southeast Asian Pteropodidae bat species such as PREDICT_CoV-16 and PREDICT_CoV-22. PREDICT_CoV-22 was also detected in Hipposideridae bat species (Hipposideros lekaguli). PREDICT_CoV-78 detected in multiple bat and rodent species of Southeast Asia also showed a high prioritization score. These model outcomes, especially the prioritization score, provide a data-driven tool to quantify zoonotic risk for novel viruses. Even though the model is trained on numerous data points for known zoonotic and non-zoonotic viruses, individual predictions for newly discovered viruses would only require data on hosts and virus family if used within our modeling framework.

Open in a separate window

Fig. 3

Prioritization metrics for novel viruses to understand zoonotic risk.

Top ten and bottom five newly discovered viruses from six virus families (a–d) with the virus prioritization scores based on multiclass model predictions. Annotations show the score and support represented by the number of human links predicted.

Development of ${\mathit{G}}_{\mathit{c}}$

a. Centrality measures of observed network (${\mathit{G}}_{\mathit{c}}$)

To test if centrality measures (degree centrality, betweenness centrality, eigenvector centrality, clustering coefficient) for viral nodes in the observed network ($G_c$) vary significantly between viral families, we firstly used the Kolmogorov-Smirnov (KS) test. KS test is routinely used to identify distances between cumulative distribution functions of two probability distributions and is largely used to compare degree distributions of networks^22,23. For each viral family, distributions of centrality measures (degree centrality, betweenness centrality, and eigenvector centrality) and clustering coefficient within the observed network ($G_c$) were compared with the distribution of all nodes in the network using the two-tailed KS test. Secondly, a linear regression model with virus family as a categorical variable and the number of PubMed hits as a covariate to adjust for sampling bias were fitted to understand associations of viral families with centrality measures.

After fitting the model, node-level permutations were implemented. For each random permutation, the output variable was randomly assigned to covariate values and the model was re-fitted. Finally, a p-value was calculated by comparing the distribution of coefficients from permutations with the original model coefficient.

Network topology feature selection

Using the observed network (${\mathit{G}}_{\mathit{c}}$), multiple network topological features for all node (virus) pairs were calculated. The following are topological network features calculated. Features data type, definition and methods to calculate these features are presented in Table S3.

1. The Jaccard coefficient: a commonly used similarity metric between nodes in information retrieval, is also called an intersection of over the union for two nodes in the network. In the unipartite network generated here, it represents the proportion of common neighbor viruses from the union of neighbor viruses for two nodes. Neighbor viruses are defined as viruses with which the virus shares at least a single host.

2. Adamic/Adar (Frequency-Weighted Common Neighbors): Is the sum of inverse logarithmic degree centrality of the neighbors shared by two nodes in the network²⁴. The concept of Adamic Adar index is a weighted common neighbors for viruses in the network. Within network prediction, the index assumes that viruses with large neighborhoods have a less significant impact while predicting a connection between two viruses compared with smaller neighborhoods.

Both Jaccard and Adamic Adar coefficients have been routinely used for generalized network prediction and have shown high accuracy in predicting missing links in networks, specifically bipartite networks²⁵, the information flowing through neighborhoods formed by two nodes might not always be enough to have similar predictive power in an unipartite network. This warrants use of other topology features along with neighborhood-based features.

3. Resource allocation: Similarity score of two nodes defined by the weights of common neighbors of two nodes. Resource allocation is another measure to quantify the closeness of two nodes in the network and hence to understand the similarity of hosts they infect.

4. Preferential attachment coefficients: The mechanism of preferential attachment can be used to generate evolving scale-free networks, where the probability that a new link is connected to node x is proportional to k²⁶.

5. Betweenness centrality: For a node in the network betweenness centrality is the sum of the fraction of all-pairs shortest paths that pass through it. The feature that we used for training the supervised learning model was the absolute difference between of betweenness centralities of two nodes. The difference between the betweenness centrality represents the difference in the sharing observed by two viruses in the pair.

6. Degree centrality: The degree centrality for a node v is the fraction of nodes it is connected to. The feature that we used for training the supervised learning model was the absolute difference between degree centralities of two nodes. Unlike the difference in the betweenness centrality, the difference in degree centrality only looks at the difference in the number of observed host sharing.

7. Network clustering: All nodes were classified into community clusters using Louvain methods²⁷. A binary feature variable was generated to describe if both the nodes in the pair were part of the same cluster or not. If both viruses are from the same cluster, it represents a similar host predilection than when both viruses are not from the same cluster hence accounting for the evolutionary predilection of viruses (or virus families) to infect a certain type of host.

These topological network characteristics come with certain limitations when it comes to the unipartite network of viruses with links formed due to shared hosts and might not truly represent the flow of information between nodes as compared to a bipartite network. Therefore, to account for these limitations, we use multiple network features as weak learners in our model building characteristics summarizing the network through the use of several quantitative metrics. In addition to this, we estimated the feature importance of these metrics in predicting missing links between viruses to quantify the information pasting through these links.

Pearson’s correlation coefficients were calculated to identify highly correlated features and for choosing features for model training (Fig. S33). Virological features included in model training were categorical variables describing the virus family of both the nodes in the pair, followed by a binary variable if both the viruses belong to the same virus family. During the model development, PubMed hits generated three predictive features for each pair of viruses on which model training and predictions were conducted. These included two features representing PubMed hits for the two viruses in the pair (PubMed_V1, PubMed_V2) and the absolute difference between PubMed_V1 and PubMed_V2 to account for differences in sampling bias between the two viruses.

Cross-validation and fitting generalized boosting machine (GBMs) models

A nested-cross-validation was implemented for the binary model while simple cross-validation was implemented for the multiclass model (multiple output categories). The parameters of the binary model were first hyper-tuned using a cross-validated grid-search method. Values were tested using a grid search to find the best-performing model parameters that showed the highest sensitivity (recall). The parameters tested for hypertuning and their performance are provided in the supplementary material (supplementary results and Table S5). For further cross-validation of the overall binary model, all the viruses were randomly assigned to five groups. For each fold, the viruses assigned to a group were dropped from the data, and a temporary training network (${\mathit{G}}_{\mathit{t}}\mathit{)}$ was constructed, assuming that this represented the current observed status of the virus-host community. For all possible pairs in ${\mathit{G}}_{\mathit{t}}$ (both that sharing and not sharing any hosts) ten topological and viral characteristics were calculated as training features (Table S4). Categorical features were one-hot-encoded and numeric features were scaled. An XGBClassifier model with binary: logistic family was trained using the feature dataset to predict if virus pairs share hosts (1,0 encoded output). The cross-validation was also used to determine the optimum decision threshold for determining binary classification (Fig. S6) and a precision-recall curve was used to identify positive predictive value and sensitivity at the optimum threshold (Fig. S8).

The multiclass model was implemented in the same way, creating an observed network ($G_c$) based on species-level sharing of hosts and randomly dropping viruses to generate a training network ($G_t$) to train the XGboost model. The output variables were generated based on the taxonomical orders of shared hosts. A pair of viruses can share multiple hosts, hence we trained a multioutput-multiclass model. Humans were considered an independent category of taxonomical order (label) and were given a separate label from primates. For fine-tuning the multiclass model, we started with the best performing parameters of the binary model and manually tested 5 combinations of model parameters by adjusting values of the learning rate, number of estimators, maximum depth, and minimum child weight (Supplementary code and results).

We used three methods to estimate the importance of features for our binary model. Specifically, improvement in accuracy brought by branching based on the feature (gain), the percentage of times the feature appears in the XGboost tree model (weight), and the relative number of observations related to the specific feature (cover). Results for feature importance are shown in supplementary results (Fig. S10).

Missing links for novel viruses, binary and multiclass prediction

The wildlife surveillance data represented a sampling of 99,379 animals (94,723 wildlife, 4656 domesticated animals) conducted in 34 countries around the world between 2009–2019 (Table S6)¹. Specimens were tested using conventional Rt-PCR, Quantitative PCR, Sanger sequencing, and Next Generation Sequencing protocols to detect viruses from 28 virus families or taxonomic groups (Table S7). Testing resulted in 951 novel monophyletic clusters of virus sequences (referred to as novel viruses henceforth). Within 951 novel viruses, 944 novel viruses had vertebrate hosts that were identified with certainty based on barcoding methods and field identification. Host species identification was confirmed by cytochrome b (cytb) DNA barcoding using DNA extracted from the samples²⁸. We predicted the shared host links between novel viruses and known viruses using binary and multiclass models in the following steps. Out of 944 novel viruses discovered in the last ten years, we were able to generate predictions for 531 novel viruses that were detected in species already classified as hosts within the network. The remaining 413 viruses were the first detection of any virus in that species and thus host associations could not be informed by the observed network (${\mathit{G}}_{\mathit{C}}$) data.

1. A new node representing the novel virus was inserted in the observed network (${\mathit{G}}_{\mathit{c}}$). Using the list of species in which the novel virus was detected, new edges were created with known viruses that are also known to be found in those hosts. This generated a temporary network for the novel virus (${\mathit{G}}_{temp}$). If the novel virus was not able to generate any edges with known viruses, meaning the host in which they have been found was never found positive for any known virus, predictions were not performed.

2. Using ${\mathit{G}}_{temp}$ feature values were calculated for the novel virus (betweenness centrality, clustering, and degree). For all possible pairs of the novel virus with known viruses that are not yet connected with each other through an edge in ${\mathit{G}}_{temp}$ a feature dataset was generated (Jaccard coefficient_{(novel virus, known virus)}, the difference in betweenness centrality of the novel virus and known virus, if the novel virus and known virus were in the same cluster, the difference in degree centrality_{(novel virus, known virus)}, if the novel virus and known virus were from same virus family, the difference in PubMed hits_{(novel virus, known virus)}, PubMed hits for the novel virus, PubMed hits for the known virus). Studies and nucleotide sequences for novel viruses are expected to be published and shared on PubMed’s Nucleotide database and in various peer-reviewed publications. Data associated with GenBank accession numbers and nucleotide sequences for novel viruses are presented in Supplementary Data 3 and Supplementary Data 4 respectively. At the time of development of the model, data for all viruses was not shared in a format that would reflect on PubMed’s database, we decided to use the number of unique species the virus was detected in the last ten years of wildlife surveillance conducted by the USAID PREDICT project. These detections will be reflected in PubMed’s Nucleotide database and search term eventually, hence we considered them as a proxy for search terms conducted for known viruses. Currently, evaluation of the effects of this substitution of PubMed hits with the number of detections for novel viruses is not possible with limited data on novel viruses but needs to be reevaluated as more studies are published on these novel viruses. To further evaluate the association between PubMed hits through search term and Genbank hits, we ran a generalized linear regression model with PubMed hits as dependent variable and Genbank hits as intendent variable, accounting for virus families.

The results indicated that Genbank hits had statistically significant predictive value in predicting PubMed hits (β= 0.72, p<0.005) even after accounting for various virus families. Multiple virus families showed statistically different estimates than the reference virus family (Adenoviridae) indicating a significantly different association than other virus families. Results of the generalized linear regression model are presented in Table S8.

3. Using this dataset for the novel virus, a binary presence of a link between the novel virus and known viruses was predicted using the trained binary model. The taxonomic order of the host link was predicted using the trained multiclass model.

4. For each possible link, the binary model predicted the probability of sharing a link, and the multiclass model predicted multivariate outcomes of taxonomic orders and associated probabilities. A threshold of 0.70 for the binary prediction model was used to classify if the link is present or not and only those links were explored for their corresponding multiclass model outputs.

5. The multiclass model showed higher performance for correctly classifying links as “human” hosts than other numerous avian and mammalian taxonomic orders. Hence, the multiclass model outputs were summarized into either humans or other taxonomic groups. For the novel virus, a list of known viruses with the predicted link was generated. Using the hosts of these known viruses and the taxonomic order in which the novel virus was detected, a list of most likely species was generated based on the overall frequency of the host species. For understanding the likelihood of infecting humans two factors were considered to be of importance. Firstly, the number of links where humans are predicted as shared hosts with known viruses ($n$) and the average model-predicted probability of those links. A representation was generated incorporating the probability and available model support in terms of number links to reflect the likelihood and compare viruses relative to each other.

To test if virus family, the taxonomic order of hosts in which novel viruses were detected, and the number of times the viruses were detected (equivalent to PubMed hits for known viruses) influenced node (virus) level network centrality measures in the predicted network (${\mathit{G}}_{\mathit{p}}$) a linear regression model was fitted with centrality measures.

For each of the random 10,000 node-level permutations, the output variable (centrality measure) was randomly assigned to covariate values and the model was re-fitted. A p-value was calculated by comparing the distributions of coefficients with the original model coefficient. These models were fitted for degree centrality, betweenness centrality, eigenvector centrality, and clustering coefficient of novel viruses in the predicted network.

Prioritization score for novel viruses

Generalized Linear Mixed Models were used to understand the association effects of virus family, taxonomic order of the host and PubMed hits on the number of predicted human links and mean probability of the predicted links. The models were fit using glmmTMB and glm packages in R. For relative comparison of zoonotic risk and for prioritizing novel viruses for further characterization, a prioritization metric was developed based on the predicted probability of sharing the humans as hosts with known viruses ($p_{sharinghumans}$) and the number of predicted shared human links ($n_{humans}$) in the predicted network for the given virus ($G_{predicted}$). Distributions for both $p_{sharinghumans}$ and $n_{humans}$ were normalized and multiplied to generate a single score for a virus and for appropriate relative comparisons between novel viruses. To understand the behavior of the prioritization score when predicting the zoonotic risk of novel viruses, we also compared prioritization scores of known zoonotic and non-zoonotic viruses using the Kolmogorov-Smirnov test.