Predicting the animal hosts of coronaviruses from compositional biases of spike protein and whole genome sequences through machine learning.

The COVID-19 pandemic has demonstrated the serious potential for novel zoonotic coronaviruses to emerge and cause major outbreaks. The immediate animal origin of the causative virus, SARS-CoV-2, remains unknown, a notoriously challenging task for emerging disease investigations. Coevolution with hosts leads to specific evolutionary signatures within viral genomes that can inform likely animal origins. We obtained a set of 650 spike protein and 511 whole genome nucleotide sequences from 222 and 185 viruses belonging to the family Coronaviridae, respectively. We then trained random forest models independently on genome composition biases of spike protein and whole genome sequences, including dinucleotide and codon usage biases in order to predict animal host (of nine possible categories, including human). In hold-one-out cross-validation, predictive accuracy on unseen coronaviruses consistently reached ~73%, indicating evolutionary signal in spike proteins to be just as informative as whole genome sequences. However, different composition biases were informative in each case. Applying optimised random forest models to classify human sequences of MERS-CoV and SARS-CoV revealed evolutionary signatures consistent with their recognised intermediate hosts (camelids, carnivores), while human sequences of SARS-CoV-2 were predicted as having bat hosts (suborder Yinpterochiroptera), supporting bats as the suspected origins of the current pandemic. In addition to phylogeny, variation in genome composition can act as an informative approach to predict emerging virus traits as soon as sequences are available. More widely, this work demonstrates the potential in combining genetic resources with machine learning algorithms to address long-standing challenges in emerging infectious diseases.

Introduction

The ongoing COVID-19 pandemic remains a significant public health emergency. Since the first identified cases in China in December 2019, this outbreak of respiratory disease has developed into a global crisis, with over 100 million cases worldwide to date [1]. The causative virus was termed ‘severe acute respiratory syndrome-related coronavirus 2’ (SARS-CoV-2) [2] and is a previously unknown betacoronavirus that likely emerged through zoonotic transmission from contact with non-human animals [3,4]. However, the precise origins of the current pandemic remain inconclusive at present [5].

Two other betacoronaviruses have zoonotically emerged to cause significant human epidemics. Severe acute respiratory syndrome-related coronavirus (SARS-CoV) emerged in China in 2002 via an intermediate host of masked palm civets (Paguma larvata) in live animal markets [6,7], and Middle East respiratory syndrome related coronavirus (MERS-CoV) emerged in Saudi Arabia in 2012 via an intermediate host of dromedary camels [8,9], with considerable evidence that both originated in bats [10-13]. Four additional coronaviruses are known to be endemic within humans, causing mild common cold-like illness (Alphacoronavirus: Human coronaviruses 229E and NL63; Betacoronavirus: Human coronaviruses HKU1 and OC43).

All viruses in the family Coronaviridae feature similar structural proteins, including a spike glycoprotein on the outer viral surface. This protein allows entry to host cells via attachment to cell receptors at the receptor binding domain (RBD) region and subsequent fusion with cell membranes. As such, this protein exhibits high variation even between closely related coronaviruses that may correspond to profound differences in receptor usage and/or different tissue and host preferences [14]. While SARS-CoV-2 shows high genetic similarity to bat coronaviruses, particularly bat coronavirus RaTG13 (matching 96% sequence identity) [4], its spike protein instead exhibits differences among key amino acid residues of the RBD [3,15]. Based on this, SARS-CoV-2 is predicted through structural [15,16] and in vitro experimental models [17,18] to have highly efficient binding to the human angiotensin-converting enzyme 2 (ACE2) receptor, a feature that has likely contributed to its efficient human-to-human transmissibility. Several methods, e.g., ENC analysis of codon usage for SARS-CoV-2 [19] and positive selection analysis for MERS-CoV [20,21] have demonstrated that the spike proteins of these coronaviruses experienced substantial selection pressure during emergence compared to other genomic regions. Therefore, as a key molecular determinant of host range [14,22], evolutionary adaptation of coronavirus spike proteins represent a key opportunity to further understand their host range constraints.

Beyond selection acting at specific loci, viral adaptation can also manifest through broad-scale genomic signatures. Viruses exhibit biased genome composition, for example, in non-uniform use of synonymous codons [23]. Furthermore, coevolution within different hosts may indirectly lead to selection for particular compositions, as reported for nucleotide and dinucleotide usage within avian and human influenzaviruses [24,25] and codon pair usage of arboviruses within their insect vectors and mammalian hosts [26]. The Coronaviridae are no exception—different coronaviruses (including SARS-CoV-2) vary in their genome composition, with particularly complex patterns of codon usage in spike protein coding sequences [27], which could potentially contain important evolutionary signal regarding host origin.

Machine learning has recently gained substantial attention as a methodology in comparative modelling of emerging diseases. These methods are capable of decomposing signal in high-dimensional genomic information (a limitation of regression frameworks) without the need for sequence alignment. Genomic machine learning analyses have demonstrated the ability to not only classify viruses from recurring viral genome motifs [28], but also classify their broad host origins [29-32]. Specifically considering coronaviruses, support vector machines and random forests have been trained on various genomic features to predict host group, including nucleotide and dinucleotide biases [33], amino acid composition [34] or sequence k-mers [35]. However, previous model predictions are mostly concentrated upon bats or humans, and few analyses explicitly address the spike protein (but see [35]). The exact potential of genome composition to predict host origin therefore remains unclear.

We aimed to use machine learning to understand how the complex genomic signatures of coronaviruses might predict their hosts and determine the importance of such signature in the spike protein. Specifically, we trained random forest models on compositional biases for spike protein and whole genome nucleotide sequences and compared their performance. A limitation of these approaches is that model predictions can be a) strongly influenced by the inherent viral sampling biases in available data [36] and b) reflect evolutionary relatedness of viral lineages rather than true host-associated signal [37]. Therefore, we conduct data thinning of overrepresented viruses before training models and control for similarity between related sequences by holding out entire species during cross-validation. We demonstrate the use of machine learning as a reliable method to estimate host origins of future novel coronaviruses in humans and livestock.

Materials and methods

Data extraction and processing

Spike protein or whole genome sequence data for coronaviruses were identified within GenBank, using search terms where successive searches were conducted replacing ### with taxonomic identifiers for each species and unranked sub-species belonging to the family Coronaviridae within the NCBI taxonomy database [38] (n = 1585 taxonomic ids total). Matching sequences were then extracted and further filtered to exclude incomplete or truncated sequences based on a) metadata labels and b) length restrictions, discarding any spike protein sequences < 2 kilobases (kb) and any whole genome sequences outside a range of 20kb– 32kb. We accepted both spike protein coding sequences within whole genome sequences and standalone complete spike protein sequences, excluding those only covering individual S1 or S2 subunits. All sequence data searching, filtering and extraction was conducted with R package `rentrez`v1.2.2 [39] (see also [40]).

Host classification

For each spike protein or whole genome sequence, host names were extracted from the host organism metadata field before being resolved to the standard NCBI taxonomy using the R package `taxizedb`v0.1.9.93 [41] (see also [40]). Host names were automatically resolved to the highest taxonomic resolution possible and any unmatched host names were resolved manually, discarding sequences with missing/unresolvable names.

We then constructed a new variable broadly describing host category of each sequence, defined at various taxonomic levels: human (species Homo sapiens), camelid (family Camelidae), swine (family Suidae), carnivore (order Carnivora), rodent (order Rodentia), and bird (class Aves). Following a previous analysis [29], we included two categories to represent bats (order Chiroptera): suborder Yinpterochiroptera (families Craseonycteridae, Hipposideridae, Megadermatidae, Pteropodidae, Rhinolophidae, and Rhinopomatidae) and suborder Yangochiroptera (all other families), based on their evolutionary divergence [42] and differences in ecology and host-virus relationships [43,44]. Sequences not conforming to any of the above host categories were excluded from further analysis.

Genomic feature calculation

We then calculated several features describing genome composition biases of each spike protein and whole genome coding sequence at nucleotide, dinucleotide or codon level. Firstly, nucleotide biases were calculated as simple proportion of A, C, G or U content. Dinucleotide biases were calculated as the ratio of observed dinucleotide frequency to expected based on nucleotide frequency, following [29]: where d denotes frequency of dinucleotide xy, n and n denote frequency of individual nucleotides x and y, and D denotes total dinucleotides and N total nucleotides for length of the given sequence. Biases were calculated separately for each dinucleotide at each position within codon reading frames (i.e., at positions 1–2, 2–3 or 3–1) as dinucleotides spanning adjacent codons may be subject to more extreme biases [45,46]. Finally, Relative Synonymous Codon Usage (RSCU) was also calculated for each codon including stop codons, following [47]: where n denotes number of codons synonymous for amino acid i and c denotes frequency of j codon encoding for such amino acid. In total, this gave 116 genomic features for use in predictive models (4 nucleotide biases, 16 * 3 dinucleotide biases, 64 codon biases). The Effective Number of codons (ENc) [48] was also calculated for each sequence as a summative metric of magnitude of codon bias. All calculations for whole genomes considered nucleotide sequences as-is rather than as-read; sequence strings duplicated by frameshifting among the ORF1a and ORF1b replicase protein were discarded to avoid disproportionate weighting in modelling analyses. Nucleotide and dinucleotide frequencies were obtained using R package `Biostrings`, v2.56.0 [49], and codon frequencies and ENC values were obtained using R package `cordon`, v.1.6.0 [50].

To assess over- or underrepresentation of codons among spike protein and whole genome datasets, we conducted one-sample t-tests comparing RSCU for each codon to a null value of 1 (excluding stop codons and codons with no other synonymous codons), applying Bonferroni corrections for multiple testing. Clustering of RSCU within-genus was assessed by constructing heatmaps using package `gplots`, v3.0.3 [51] before extracting associated dendrogram clusters and calculating Normalised Mutual Information (NMI) with genus, a measure of co-occurrence scaled between 0 (no mutual information) and 1 (perfect information) (see [52]). NMI was calculated using function `NMI()`in package `aricode`, v.1.0.0 [53], cutting dendrograms into k = 10 clusters and excluding viruses unclassified into a genus.

Machine learning analysis

To quantify the potential for genome composition biases to predict coronavirus host category, we used random forests, an ensemble machine learning approach that aggregates over a large number of individual classification tree models [54]. We selected random forests as a method that natively handles multiclass classification problems, intuitively handles interactions between predictor features, and can provide interpretability in feature-outcome explanatory relationships (see below) [55].

These predictive modelling methods are often sensitive to training data composition, a potential problem given the heavily biased patterns of viral sampling [36]. Several viruses appeared highly overrepresented in extracted data, for example, 1802 spike protein and 555 whole genome sequences were available for porcine epidemic diarrhea virus (S1 Table), compared to the overall median of 1 sequence per virus in both cases. As viral sequences (and therefore, genome composition) are expected to be highly similar within-species or subspecies, we therefore conducted a heuristic data thinning procedure prior to modelling, sampling a maximum of twenty sequences per host category per virus, thinning for approximately 5% of viruses (S1 and S2 Tables).

Machine learning methods can also be sensitive to class imbalances in outcomes [56] and slight imbalances between host categories remained after data thinning (S2 Table). Preliminary models suggested predictive performance was robust to this imbalance, being comparable or decreasing when alternative resampling methods were applied (S1 Methods and S3 and S4 Tables).

Zoonotic and epizootic coronaviruses are unlikely to have experienced substantial coevolution within the novel host following cross-species transmission. Instead, their evolutionary signature in genome composition is much more likely to reflect the original or donor host. Therefore, zoonotic or epizootic coronavirus sequences sampled from novel hosts (i.e., SARS-CoV, SARS-CoV-2, MERS-CoV in humans, totaling m = 9 taxonomic identifiers, full list in S1–S4 Data; swine acute diarrhea syndrome coronavirus in swine) were held out from model training. Of these, we retained the zoonotic sequences as a prediction set of interest. By generating predictions for zoonotic viruses from sequences alone, we aimed to investigate model utility when presented with the scenario of a newly identified emerging coronavirus. We excluded human enteric coronavirus from all analyses as the zoonotic potential of this virus remains unclear.

As genome composition is highly homogenous for closely related virus sequences, we established training data for random forest models using an outer loop of hold-one-out validation applied to coronaviruses (S1A Fig), i.e., rather than exclude randomly sampled sequences as a test set, we constructed and validated random forests excluding all sequences from a given virus species or unranked subspecies in each instance. This allows host predictions for unseen (i.e., novel) viruses based on values of compositional features, rather than indirectly predicting host by the proxy of viral similarity. For each outer loop, model parameters were then optimized within an inner loop of 10-fold cross-validation, dividing the data into ten folds, training models using a grid search of parameter combinations on nine and identifying the best parameters when applied to the tenth fold in each case (S1B Fig). The parameter set yielding the highest prediction accuracy on validation folds was retained (S2 Fig).

Model performance was then assessed by applying each random forest to its respective held-out coronavirus sequences as a test set. Probabilities of host categories were obtained by dropping sequences down each individual tree model within the random forest and averaging host category prevalence of resulting terminal nodes [57]. Similar predictions were generated for each zoonotic coronavirus sequence by averaging probabilities across all random forests to generate grand mean probabilities (S1A Fig). To investigate explanatory relationships between genomic biases and hosts of coronaviruses, variable importance (calculated as relative mean decrease in Gini impurity) and partial dependence (calculated as marginal probability of each host category) associated with each genomic feature were averaged across all random forests.

All random forests were constructed using 1000 trees and implemented using R package `ranger`, v0.12.1 [58]. All analyses were initiated from a fixed random seed chosen a priori. Model performance was insensitive to choice of random seed (S5 Table). All data processing and modelling were conducted within R v4.0.0 [59]. Supporting data and code are available at https://github.com/lbrierley/cov_genome.

Results

Genome composition across the Coronaviridae

In total, we identified n = 3595 nucleotide sequences for coronavirus spike proteins and n = 1815 whole genome sequences that met inclusion criteria. These were thinned to n = 650 spike protein sequences from m = 222 coronaviruses and n = 511 whole genome sequences from m = 185 coronaviruses for use in further analysis (S1 Fig and S2 Table), spanning 40 identified host genera and 58 identified host species (S1 and S2 Data).

Broadly consistent genome composition biases were observed across the diversity of all coronavirus sequences examined. For both spike protein and whole genome sequences, all A- and U-ending codons were significantly overrepresented (one-sample t-tests of RSCU, Bonferroni-corrected p < 0.001) except AUA, GCA, CUA, GGA, and GUA, and all C- and G-ending codons were underrepresented (Bonferroni-corrected p < 0.001) except UUG (and AGG in the case of spike proteins only) (Figs 1 and 2). Hierarchical clustering based on RSCU values suggested codon usage was less distinct between genera within spike protein sequences (Normalised Mutual Information value = 0.268; Fig 1) than within whole genomes (NMI = 0.428; Fig 2). However, clear separation of deltacoronaviruses was observed for both cases, as these appeared to have less extreme biases in codon usage. This was confirmed by ENc calculation; deltacoronaviruses had higher ENc values than other coronavirus genera (Table 1). Considering dinucleotide biases, compositional bias was typically more extreme for dinucleotides spanning adjacent codons, i.e., position 3–1 (S3 Fig), and the characteristic coronavirus CpG suppression was also observed.

Fig 1

Codon biases (RSCU) across coronavirus spike protein sequences examined.

Heatmaps of coronavirus codon usage bias (RSCU) associated with each codon in each spike protein sequence (n = 650). Main colour scale denotes RSCU value, a null value of 1 (black) indicating no difference in codon usage from expectation, with blue and red representing under- or overrepresentation respectively. Dendrogram colour bar denotes taxonomic genus.

Fig 2

Codon biases (RSCU) across coronavirus whole genome sequences examined.

Heatmaps of coronavirus codon usage bias (RSCU) associated with each codon in each whole genome sequence (n = 511). Main colour scale denotes RSCU value, a null value of 1 (black) indicating no difference in codon usage from expectation, with blue and red representing under- or overrepresentation respectively. Dendrogram colour bar denotes taxonomic genus.

Table 1

ENc values across coronavirus genera.

Effective Number of Codons (ENc) for coronaviruses, stratified by genus and genome sequence type. ENc values are calculated as grand means, i.e., mean ENc was calculated per coronavirus by averaging sequences before means of means were calculated per genus by averaging coronaviruses. SD denotes standard deviation.

	Spike protein sequences		Whole genome sequences
Genus	Mean	SD	Mean	SD
Alphacoronavirus	48.65	3.12	45.02	3.00
Betacoronavirus	47.77	5.04	47.03	4.79
Gammacoronavirus	46.36	1.83	46.11	0.74
Deltacoronavirus	54.79	3.43	51.75	3.22
(unclassified)	47.89	3.41	48.56	3.04
total	48.77	4.36	46.81	4.22

Host predictions of random forest models

Random forest models trained on nucleotide, dinucleotide and codon bias features of spike protein sequences predicted coronavirus hosts with 73.5% accuracy during hold-one-out cross-validation (Table 2). Genome composition of spike proteins appeared just as informative as whole genomes despite being much smaller in sequence length, as both models achieved very similar performance in all diagnostic measures (Table 2).

Table 2

Predictive performance of random forest models.

Model diagnostics describing overall performance when applied to predict host category of held-out coronaviruses not used for model training. CI denotes confidence interval, Kappa denotes Cohen’s Kappa statistic, mAUC denotes multiclass area-under-curve statistic, and F1macro denotes F1 score calculated using macro-averaging (performance on each host category weighted equally).

Predictor features	Accuracy (95% CI)	Kappa	mAUC	F1macro
Spike protein	0.735 (0.700, 0.769)	0.696	0.898	0.757
Whole genome	0.728 (0.687, 0.766)	0.688	0.902	0.758

Patterns of host-specific predictive performance were evident during hold-one-out cross-validation. Random forests trained on both spike protein and whole genome sequence compositional features most easily distinguished bird, carnivore and rodent host categories (Fig 3 and S6 and S7 Tables). Less powerful predictive performance was obtained for livestock (i.e., swine and camelid) host categories with these sequences often predicted as having bat (suborder Yangochiroptera) hosts, including all MERS-CoV sequences sampled from camels.

Fig 3

Random forest host predictions based on coronavirus genome composition.

Stacked bar plots of predicted probabilities of each host category for coronavirus sequences. Predictions were obtained from ensemble random forest models trained on A) spike protein and B) whole genome composition features. Panels depict sequences from each metadata-derived host category and colour coding denotes model-predicted host category. Stacks represent individual coronavirus sequences, ordered from largest to smallest probability of the correct host, i.e., greater panel area matching the correct host category indicates better overall model performance. Non-zoonotic coronavirus sequences originating from humans (human coronaviruses HKU1, NL63, OC43, 229E) are labelled for clarity. Versions stratified by genera and species are provided as S4 and S5 Figs.

Human host origins appeared particularly difficult to characterise, with model-predicted hosts appearing more uncertain using spike protein features than whole genome features (Fig 3); while human coronaviruses HKU1 and NL63 were more confidently correctly classified based on whole genomes, human coronaviruses OC43 and 229E were more confidently misclassified as having camelid or Yinpterochiroptera hosts. The reciprocal also only occurred using whole genomes, i.e., several camel coronaviruses were predicted to have human hosts.

We then applied random forests to those sequences of zoonotic viruses sampled from humans and excluded from model training: SARS-CoV, SARS-CoV-2, and MERS-CoV (S1 Fig and S3 and S4 Data). As these viruses have experienced little coevolution following zoonotic spillover, their genome composition signal likely gives indications about their ultimate or proximate animal host origins. MERS-CoV was overwhelmingly predicted to have camelid hosts (Fig 4) and SARS-CoV was predicted with less certainty as having carnivore hosts, consistent with the respective known intermediate hosts of camels and palm civets (order Carnivora). Contrastingly, SARS-CoV-2 was predicted mostly strongly to have a bat (suborder Yinpterochiroptera) host. Host predictions for zoonotic viruses were consistent between models using spike protein and whole genome features (Fig 4).

Fig 4

Random forest predictions based on zoonotic coronavirus genome composition.

Stacked bar plots of predicted probabilities of each host category for zoonotic coronavirus sequences sampled from humans. Predictions were obtained from ensemble random forest models trained on A) spike protein and B) whole genome composition features. Colour coding denotes model-predicted host category. Stacks represent individual coronavirus sequences.

Variable importance of random forest models

The most informative genomic features towards predicting coronavirus hosts were a mixture of dinucleotide and codon biases (Fig 5), with dinucleotide biases appearing slightly more informative for spike protein sequences and codon biases appearing slightly more informative for whole genome sequences. However, predictive power of individual genomic features did not hold between spike protein and whole genome sequences; only weak correlation was observed between ranked variable importance from both analyses (Spearman’s rank, ρ = 0.191, p = 0.042) (Figs 5, S6 and S7). Partial dependence plots suggested the strongest individual discriminating feature to be GG dinucleotides at positions 1–2; an overrepresentation of this dinucleotide within the spike protein sequence clearly distinguished bird hosts from mammalian hosts (S6 Fig), consistent with the greatest predictive performance for bird coronaviruses (S6 Table).

Fig 5

Variable importance of genomic features.

Variable importance of genome composition features in ensemble random forest models predicting coronavirus host category from whole genome sequences (x axis) and spike protein sequences (y axis), with labelling of top ten most informative features from both analyses. Points denote mean values of relative decrease in Gini impurity associated with each feature across A) m = 222 and B) m = 185 random forests during hold-one-out cross-validation. Colour key denotes genomic feature type.

Discussion

We observe biased dinucleotide and nucleotide usage across the family Coronaviridae, and demonstrate that these genome composition biases contain sufficient evolutionary signal such that they can predict animal host origin. We show that training random forests on these features of spike proteins is equally as informative as using whole genome sequences in predicting hosts of novel (i.e., unseen during model training) coronaviruses, with bird, carnivore and rodent viruses having the highest prediction accuracy. When applied to human coronavirus sequences from previous epidemics (SARS-CoV, MERS-CoV), random forest model predictions consistently represented the intermediate hosts. In the case of SARS-CoV-2, where the exact transmission pathway remains unknown, models predicted sequences to have a bat host (suborder Yinpterochiroptera).

Variability in genome composition

Among our dataset of 222 coronaviruses, we observed A- and U-ending codons to be overrepresented and C- and G-ending codons to be underrepresented (Figs 1 and 2), a commonly noted trait in other studies [19,60]. Elsewhere, CpG dinucleotide bias has been proposed as a specific determinant of host (and tissue) range of coronaviruses [61], on the basis that CpG dinucleotides are targeted by zinc finger antiviral proteins and their suppression is therefore linked to immune evasion. We observed consistent CpG suppression (S3 Fig) and CG dinucleotides (positions 1–2) ranked 6th and 8th in feature importance for spike proteins and whole genomes (Fig 5), respectively.

However, spike proteins display different patterns of codon usage from other viral proteins [27], reflected in the lack of correlation between genome composition feature importance in random forests trained on spike proteins versus whole genomes (Figs 5, S6 and S7). This indicates spike proteins contain complex evolutionary signatures reflecting their distinct role in host-virus coevolution, supporting our approach in using many features beyond single dinucleotides [62]. The same evolutionary signatures would not be expected of proteins not involved in direct molecular interactions with hosts; in a comparative model trained on genome composition of the envelope protein (which has codon usage consistent with mutational bias rather than selective pressure [19]) predictive performance was reduced (S8 Table). Further work is needed to address the individual products of the ORF1ab polyprotein as an additional area of the genome under selective pressure for codon usage [19].

Model predictions of human coronaviruses

During model validation, human hosts appeared more challenging to correctly predict than other host types. Although the endemic human coronaviruses are common, they are also thought to have their ultimate evolutionary origins within non-human animals [11], which may explain this difficulty. In particular, human coronaviruses NL63 and HKU1 were more consistently predicted as having human hosts than human coronaviruses OC43 and 229E, especially when using whole genome sequences (Fig 3). Human coronavirus NL63 is estimated to have a more ancient common ancestor with bat coronaviruses than 229E [63,64], implying longer coevolutionary history within human hosts has resulted in a more consistently identifiable genomic signature.

Although several mutations of SARS-CoV-2 are becoming fixed in the population, e.g. D614G in the spike protein [65], the virus has experienced only weak purifying selection [66] and sequences remain extremely similar over the course of the pandemic. As such, our approach cannot identify host adaptation “in real-time”; rather, we examine variation generated over much longer macroevolutionary histories.

Instead, we would expect viruses that have transmitted cross-species more recently to retain the genome composition signature of their previous hosts, having experienced little coevolution within the novel host. Applying our finalised models to zoonotic human virus sequences may therefore give an indication of their proximate or ultimate animal host origin (Fig 4).

Human sequences of SARS-CoV were predicted to have a carnivore host, consistent with the known intermediate host of palm civets (Paguma larvata). Much previous work has shown human and civet SARS-CoV sequences to have high similarity, with adaptive mutations concentrated within the spike protein (specifically, the receptor binding domain) [6,7,22], which may explain the stronger prediction of carnivore hosts when using spike proteins than whole genome sequences (mean (SD) predicted probabilities: 0.902 (0.007) vs 0.688 (0.024)). Similarly, human sequences of MERS-CoV were strongly predicted as having camelid hosts, consistent with the intermediate host of dromedary camels. The detection of evolutionary signatures corresponding to these intermediate hosts implies that these coronaviruses circulated in those hosts for sufficient time for coevolution to shape genome composition before cross-species transmission to humans. For MERS-CoV, camel infections have been recognised as far back as at least the 1980s [67,68].

The origins of SARS-CoV-2 have been heavily speculated upon since its discovery, though there remains no compelling evidence towards the animal source of the first human infections. Our random forest models trained on genome composition of both spike protein and whole genome sequences predicted SARS-CoV-2 as having a bat host (suborder Yinpterochiroptera). Alignment-based and phylogenetic approaches showed the most closely related virus to be bat coronavirus RaTG13, a virus sampled from a horseshoe bat (Rhinolophus affinis) belonging to this suborder [4], and more widely, the Rhinolophidae family are the most likely ancestral hosts of the Sarbecovirus genus [69].

While our model predictions support bats as the ultimate origin of SARS-CoV-2, the involvement of intermediate hosts remains unclear. Although the Malayan pangolin (Manis javanica) was proposed early in the pandemic [70,71], recent analyses have argued there is absence of evidence for this [72,73]. The methods used here are unable to identify intermediate hosts without sufficient sequence availability, and lack of such from pangolins (order Pholidota) preclude us from directly testing this hypothesis. However, selection analyses indicate SARS-CoV-2 could reasonably have exhibited efficient human infectivity and human-to-human transmissibility following direct transmission from bats [66,73], i.e., without strict need for prolonged selection within an intermediate host.

Future directions

A natural comparison to these methods is phylogenetic analyses, which can estimate traits such as host type from reconstructing viral ancestry based on sequence similarity. There is challenging confounding between molecular characteristics and sequence similarity, i.e., variation in genome composition may actually be predictive of viral lineage rather than host type [37], essentially acting as a proxy for phylogenetic relatedness. To separate these signals, viral similarity needs to be considered in model construction [31]; by using a cross-validation procedure holding out entire coronavirus species or sub-species, we attempt to distinguish genomic signatures arising from convergent evolution within specific hosts, rather than from viral similarity. A more generalised scope of study across multiple viral families would allow holdout of entire families during cross-validation [31], removing any further phylogenetic proxy effects.

Additional challenges are created by the unavoidable, systematic gaps in sampling coverage. For example, disproportionate sampling to identify viruses in wildlife similar to those already known to affect humans or livestock may introduce bias to predictive models [36]. Although we address this by data thinning, our model predictions, particularly for zoonotic coronaviruses, are likely influenced by the range of known viruses with available sequence data. These issues highlight the need for a careful choice of training dataset in modelling studies, but also for wider sampling and surveillance of coronaviruses among the wild virome [74], especially considering their high public health risks.

Although we focus on compositional features, predictive approaches using other sequence properties may improve more mechanistic understanding of host range. For example, amino acid composition and physicochemical similarity between contiguous amino acid residues can predict human origin of coronavirus spike sequences [34]. Similarly, Young et al. have recently demonstrated the use of multiple types of genomic features in combination to predict infected hosts and found physicochemical classification of amino acid k-mers to achieve similar predictive power to nucleotide k-mers [31]. More widely, hydrophobic and hydrophilic composition of host receptors shows some predictive signal towards virus sharing [75], hydropathy being of mechanistic importance during virus-receptor binding, e.g. for murine coronavirus [76]. These properties could be used as additional features and improve machine learning model-derived predictions.

This emphasises an additional key question for future modelling studies distinct from host origin—whether genomic traits can predict the zoonotic potential of newly discovered animal coronaviruses [77]. As this is strongly determined by molecular mechanisms of virus-host receptor interactions [22], these predictions may be best inferred by model frameworks combining genomic features of both spike proteins and host receptors.

Conclusion

By training machine learning models on genome composition across the Coronaviridae, we demonstrate a detectable evolutionary signature predictive of host type rooted in a region of the genome that is key to host shifts. Our random forest predictions add to the growing evidence COVID-19 ultimately originated within bats, though further work is needed to understand the potential for intermediate hosts. Characterising spike proteins (and by extension, their interaction with host receptors) may provide a fruitful path to further understanding zoonosis risk among coronaviruses.

Alternative sampling methodologies for random forest training datasets.

(DOCX)

Click here for additional data file.

Spike protein sequences used in analysis with metadata and model-predicted host category.

List of non-zoonotic coronavirus genome sequences used in spike protein machine learning analysis along with genus, GenBank accession ID, metadata-derived host category and predicted host category from random forest when respective coronavirus was held-out during validation. Coronavirus names refer to taxonomic entities as defined by NCBI taxonomy (whether species or unranked sub-species). Model predictions are given as single category with highest probability, along with probabilities associated with each individual category.

(CSV)

Click here for additional data file.

Whole genome sequences used in analysis with metadata and model-predicted host category.

List of non-zoonotic coronavirus genome sequences used in whole genome machine learning analysis along with genus, GenBank accession ID, metadata-derived host category and predicted host category from random forest when respective coronavirus was held-out during validation. Coronavirus names refer to taxonomic entities as defined by NCBI taxonomy (whether species or unranked sub-species). Model predictions are given as single category with highest probability, along with probabilities associated with each individual category.

(CSV)

Click here for additional data file.

Spike protein sequences of zoonotic human coronaviruses with metadata and model-predicted host category.

List of zoonotic coronavirus genome sequences sampled from humans used in spike protein machine learning prediction along with genus, GenBank accession ID, and predicted host category, aggregating over all random forests. Coronavirus names refer to taxonomic entities as defined by NCBI taxonomy (whether species or unranked sub-species). Model predictions are given as probabilities associated with each individual category.

(CSV)

Click here for additional data file.

Whole genome sequences of zoonotic human coronaviruses with metadata and model-predicted host category.

List of zoonotic coronavirus genome sequences sampled from humans used in whole genome machine learning prediction along with genus, GenBank accession ID, and predicted host category, aggregating over all random forests. Coronavirus names refer to taxonomic entities as defined by NCBI taxonomy (whether species or unranked sub-species). Model predictions are given as probabilities associated with each individual category.

(CSV)

Click here for additional data file.

Distribution of sequences represented per coronavirus.

Descriptive statistics of number of sequences per coronavirus species or unranked subspecies (i.e., unique taxonomic ids) pre- and post-data thinning procedure to a maximum of 20 sequences per host-virus combination. Data shown separately for coronavirus spike proteins and whole genome sequences.

(DOCX)

Click here for additional data file.

Number of coronavirus sequences represented per host category.

Number of genome sequences pre and post-data thinning procedure (to a maximum of 20 sequences per host-species combination) and number of coronavirus species or unranked subspecies (i.e., unique taxonomic ids) sourced from each host category. Data shown separately for coronavirus spike proteins and whole genome sequences. Note that although most coronaviruses were only known to infect a single host category, several coronaviruses infected multiple host categories and are represented across multiple counts of taxonomic ids.

(DOCX)

Click here for additional data file.

Predictive performance of random forest models using spike protein predictor features applying alternative sampling methodologies.

Model diagnostics describing overall performance using spike protein predictor features as in Table 2 when applying alternative sampling methodologies for class imbalance. Each methodology was applied to training data, retaining the same outer validation sets in all cases for comparability. CI denotes confidence interval, Kappa denotes Cohen’s Kappa statistic, mAUC denotes multiclass area-under-curve statistic, and F1macro denotes F1 score calculated using macro-averaging (performance on each host category weighted equally).

(DOCX)

Click here for additional data file.

Predictive performance of random forest models using whole genome predictor features applying alternative sampling methodologies.

Model diagnostics describing overall performance using whole genome predictor features as in Table 2 when applying alternative sampling methodologies for class imbalance. Each methodology was applied to training data, retaining the same outer validation sets in all cases for comparability. CI denotes confidence interval, Kappa denotes Cohen’s Kappa statistic, mAUC denotes multiclass area-under-curve statistic, and F1macro denotes F1 score calculated using macro-averaging (performance on each host category weighted equally).

(DOCX)

Click here for additional data file.

Variance of predictive performance of random forest models over different random seeds.

Model diagnostics describing overall performance as in Table 2, repeating analyses for an additional ten random seeds. All metrics given represent the mean diagnostic with standard deviation in brackets. Kappa denotes Cohen’s Kappa statistic, mAUC denotes multiclass area-under-curve statistic, and F1macro denotes F1 score calculated using macro-averaging (performance on each host category weighted equally).

(DOCX)

Click here for additional data file.

Predictive performance of random forest models using spike protein predictor features, stratified by host category.

Model diagnostics describing overall performance when applied to predict host category of held-out coronaviruses. Balanced accuracy denotes 0.5*(sensitivity + specificity).

(DOCX)

Click here for additional data file.

Predictive performance of random forest models using whole genome predictor features, stratified by host category.

Model diagnostics describing overall performance when applied to predict host category of held-out coronaviruses. Balanced accuracy denotes 0.5*(sensitivity + specificity).

(DOCX)

Click here for additional data file.

Predictive performance of random forest models using envelope protein predictor features.

Model diagnostics describing overall performance using envelope protein predictor features as in Table 2. The same data thinning and cross-validation methodologies were applied as in Materials and Methods, thinning 1781 down to 458 envelope protein sequences from 144 coronaviruses. CI denotes confidence interval, Kappa denotes Cohen’s Kappa statistic, mAUC denotes multiclass area-under-curve statistic, and F1macro denotes F1 score calculated using macro-averaging (performance on each host category weighted equally).

(DOCX)

Click here for additional data file.

Structured data partitioning and machine learning procedure.

Data partitioning diagram indicating the machine learning procedure used, distinguishing A) the outer loop using hold-one-out cross-validation applied to coronavirus species or unranked subspecies (with the aim of validating model performance) from B) the inner loop using 10-fold cross-validation (with the aim of optimising model parameters). Distribution of outcome classes (host category) were preserved when sampling data folds during each inner loop. Zoonotic coronavirus sequences sampled from humans that were not used for model training are also distinguished. RF denotes a random forest model, while ‘spike’ and ‘wgs’ refer to spike protein feature dataset and whole genome feature dataset, respectively.

(TIF)

Click here for additional data file.

Performance of random forests during parameter optimisation.

Performance of random forest models during inner loop of 10-fold cross-validation for parameter optimisation during A) spike protein and B) whole genome analysis. Y axis denotes prediction accuracy on test fold, X axis denotes minimum number of genome sequences in nodes at which tree algorithm continues to split, and colour denotes number of genomic features randomly considered at each splitting point. Boxes denote median and upper/lower quartiles, with whiskers extending to 1.5*IQR.

(TIF)

Click here for additional data file.

Dinucleotide biases by position within codon reading frames.

Calculated dinucleotide biases for coronavirus A) spike protein and B) whole genome sequences. Boxes denote median and upper/lower quartiles, with whiskers extending to 1.5*IQR. Colour denotes position of dinucleotide within codon reading frames, i.e., blue boxes denote dinucleotides spanning adjacent codons (position 3–1). Dashed grey line denotes null value of 1, indicating no difference in dinucleotide usage from expectation.

(TIF)

Click here for additional data file.

Random forest host predictions based on coronavirus genome composition, blocked by genus.

Stacked bar plots of predicted probabilities of each host category obtained from random forest models trained on A) spike protein and B) whole genome composition features as in Fig 3 when separated and ordered by genus (either Alphacoronavirus, Betacoronavirus, Gammacoronavirus, Deltacoronavirus, or “U” to indicated unassigned), with secondary ordering of largest to smallest probability of the correct host.

(TIF)

Click here for additional data file.

Random forest host predictions based on coronavirus genome composition, blocked by species or subspecies.

Stacked bar plots of predicted probabilities of each host category obtained from random forest models trained on A) spike protein and B) whole genome composition features as in Fig 3 when separated and ordered by species or unranked subspecies (i.e., unique taxonomic ids), with secondary ordering of largest to smallest probability of the correct host. See S1 Text for species key.

(TIF)

Click here for additional data file.

Partial dependence plots for most informative genomic features of spike proteins.

Model-predicted marginal probability of each coronavirus host category as functions of the four most informative genome composition bias features of spike protein sequences. A)–D) depict random forests trained on spike protein sequences and E)–H) depict probabilities as functions of the same features within random forests trained on whole genome sequences for comparison. Lines denote median values across A)–D) m = 225 and E)–H) m = 187 random forests during hold-one-out cross-validation. Shaded areas denote 2.5th and 97.5th percentiles. Colour key denotes host category.

(TIF)

Click here for additional data file.

Partial dependence plots for most informative genomic features of whole genome sequences.

Model-predicted marginal probability of each coronavirus host category as functions of the four most informative genome composition bias features of whole genome sequences. A)–D) depict random forests trained on whole genome sequences and E)–H) depict probabilities as functions of the same features within random forests trained on spike protein sequences for comparison. Lines denote median values across A)–D) m = 187 and E)–H) m = 225 random forests during hold-one-out cross-validation. Shaded areas denote 2.5th and 97.5th percentiles. Colour key denotes host category.

(TIF)

Click here for additional data file.

Number key for coronavirus species or unranked subspecies as depicted in S5 Fig.

(DOCX)

Click here for additional data file.

6 Jan 2021

Dear Brierley,

Thank you very much for submitting your manuscript "Predicting the animal hosts of coronaviruses from compositional biases of spike protein and whole genome sequences through machine learning" for consideration at PLOS Pathogens. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

Thank you for your submission. Both reviewers appreciated your attention to an interesting topic. We would be happy to consider a revised manuscript. However, both reviewers raise a number of comments, which are not entirely overlapping. A revision would need to adequately address these comments and would likely be sent out for re-review.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Adam S. Lauring

Section Editor

PLOS Pathogens

Adam Lauring

Section Editor

PLOS Pathogens

Kasturi Haldar

Editor-in-Chief

PLOS Pathogens

​orcid.org/0000-0001-5065-158X

Michael Malim

Editor-in-Chief

PLOS Pathogens

orcid.org/0000-0002-7699-2064

***********************

Thank you for your submission. Both reviewers appreciated your attention to an interesting topic. We would be happy to consider a revised manuscript. However, both reviewers raise a number of comments, which are not entirely overlapping. A revision would need to adequately address these comments and would likely be sent out for re-review.

Reviewer's Responses to Questions

Part I - Summary

Please use this section to discuss strengths/weaknesses of study, novelty/significance, general execution and scholarship.

Reviewer #1: The authors analyze spike protein and whole genomes of various coronaviruses using measures such as dinucleotide and codon usage biases to predict animal host using random forests. The approach appears to be surprisingly successful given the relatively simple measures used. Such an approach may be useful in determining genomic biases that are preferred by various coronaviruses as they adapt to their respective hosts.

Reviewer #2: The paper uses the compositional signatures of various coronaviruses to predict hosts with the aim of identifying the animal host of SARS-CoV-2. To do so, the authors use the random-forest machine-learning algorithm for classification, followed by analysis of the importance of each genomic feature.

While the development of algorithms for host prediction are extremely important (the current coronavirus emergence scenario is a good example), my general opinion is that such methods have to be transparent in their predictive mechanism. A black box approach, as used in the current study is potentially susceptible to a range of artefacts, of which over-fitting of compositional features through the underlying clustering of sequences by phylogeny relationships represents a substantial underlying compounding factor.

As I see it, (although I would be interested to hear the authors’ opinion on this), the basic problem is the compositional similarities used in the algorithm, such as codon choice and dinucleotides that specify amino acids, can also originate simply from virus sequences being similar to each other through genetic relatedness. Assignment of host based on compositional features might then simply represent a match to the phylogenetic clusters of viruses infecting that particular host. Indeed, although the figure is unlabelled, genus assignments seem to correlate with their host predictions in Fig. 2. Predicted hosts of pig coronaviruses (delta, alfa and beta), might correspond to the three blocks in host prediction. The alfa and bertacoronavirus blocks (on the right) furthermore look remarkably similar to the patterns of the human alfa and beta human seasonal human coronavirus. Similarly for camel viruses, where the first block might be 229E (Alfacoronavirus) and the second MERS-CoV-2-related (Betacoronavirus) based on the relative numbers of camel viruses provided in Table S1. At the very least the figure should be annotated to provide the genus assignments to demonstrate the extent to which these contribute to host prediction patterns.

In relation to that, I think it is a bit disingenuous to claim that MERS-CoV-2 was excluded from the training (lines 246-248) since the camel betacoronaviruses are very similar to MERS-CoV found in humans and it is therefore no surprise that the latter are then found to have a predicted camel origin (Fig. 3).

In terms of informative features, the observation that biases in codon positions 1 and 2 were the most predictive of host (with GG being the highest rank) might simply follow differences in amino acid composition between sequences that co-associated with genus.

**********

Part II – Major Issues: Key Experiments Required for Acceptance

Please use this section to detail the key new experiments or modifications of existing experiments that should be absolutely required to validate study conclusions.

Generally, there should be no more than 3 such required experiments or major modifications for a "Major Revision" recommendation. If more than 3 experiments are necessary to validate the study conclusions, then you are encouraged to recommend "Reject".

Reviewer #1: -The authors consider the spike protein as a special case but it may be that other proteins such as orf1ab or the RdRp may yield better results. At the very least some discussion should be added as to what would be expected if other gene(s) were analyzed. Also, a more detailed description of spike protein variation in the introduction would improve the justification for focusing on spike. The current description is fairly minimal.

-There are 116 genomic features used for the ML dataset. Can these be explicitly listed somewhere, e.g. in a supplementary table?

-The under-sampling strategy leads to a dramatic reduction in data utilized. Might it be possible to avoid this by instead introducing class weights during training?

-As far as host categories were concerned, did the model allow for the possibility of same virus but multiple hosts? Multiple bats, for example, may be able to harbor the same virus.

-The inner/outer loop cross validation strategy is unusual and makes it unclear whether the train and test sets are really independent. If they are not independent, the results may be marred by overfitting. How different are the results if only one or the other of these strategies (leave-one-out or 10-fold cross-validation) were used?

-The authors compare their "spike protein" and "whole genome" models at least twice: in the second paragraph of results ("Hierarchical clustering based on RSCU 182values suggested codon usage was less distinct between genera within spike protein sequences") and then in the discussion ("which may explain the stronger prediction of carnivore hosts when using spike protein sequences''). However, in neither case they provide any statistical support for their statements. Are these differences statistically significant? Similarly, the statement “Among our dataset of 225 coronaviruses, we observed A- and U-ending codons to be overrepresented and G- and C- ending codons to be underrepresented” in the discussion is not analyzed statistically.

Reviewer #2: 1. There is no published source code for this analysis, rendering it unable to be completely reproduced. At the minimum, all code used to generate the results and figures should be published.

2. The paper does not make clear the rationale behind using random forests for modelling the hosts of coronaviruses. Given the class imbalance was resolved with under-sampling, other statistical classification approaches such as support vector machines, logistic regression, and k-nearest neighbors may also be useful. Indeed, some of these approaches (such as k-NN) are robust to class imbalance in the first place and no explanation for their exclusion is provided in the paper. A fuller comparison of the methods and their relative performances would strengthen the methodological advances claimed in the manuscript.

3. The description of the training method is not fully clear and would be served by further explanation. For example, the reason the parameters were optimized within the inner loop is not fully explained within the main manuscript.

4. The paper is unclear whether dinucleotide frequencies, RSCU, and ENc were computed using a package and, if so, which package.

5. In relation to reproducibility, the manuscript does not describe whether a fixed random seed was used for the analysis and chosen beforehand. Doing so with the code supplied (see comment 1 above) would allow reproduction of the exact results presented as well as examination of the method over a variety of random seeds to ascertain its general performance.

6. Similarly, it is not explained in the paper how many times the analyses were run. In the absence of such information and reproducible source code, multiple testing can be a concern for the validity of the results.

7. The paper does not fully describe on the rationale for selecting twenty random samples per host category per virus for under-sampling. There are methods for under-sampling based on selecting nearest neighbors that can maintain or improve performance while also reducing class imbalance. Examples of such methods are near miss and condensed nearest neighbor under-sampling.

**********

Part III – Minor Issues: Editorial and Data Presentation Modifications

Please use this section for editorial suggestions as well as relatively minor modifications of existing data that would enhance clarity.

Reviewer #1: -Typo on line 157: should read “rather than”.

Reviewer #2: Lines 36, 39. Why are these called intermediate hosts?

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example see here on PLOS Biology: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see

5 Mar 2021

Submitted filename: Response_to_review_letter.docx

Click here for additional data file.

9 Apr 2021

Dear Brierley,

We are pleased to inform you that your manuscript 'Predicting the animal hosts of coronaviruses from compositional biases of spike protein and whole genome sequences through machine learning' has been provisionally accepted for publication in PLOS Pathogens.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Pathogens.

Best regards,

Adam S. Lauring

Section Editor

PLOS Pathogens

Adam Lauring

Section Editor

PLOS Pathogens

Kasturi Haldar

Editor-in-Chief

PLOS Pathogens

​orcid.org/0000-0001-5065-158X

Michael Malim

Editor-in-Chief

PLOS Pathogens

orcid.org/0000-0002-7699-2064

***********************************************************

Reviewer Comments (if any, and for reference):

Reviewer's Responses to Questions

Part I - Summary

Please use this section to discuss strengths/weaknesses of study, novelty/significance, general execution and scholarship.

Reviewer #1: The revised version is satisfactory, and the manuscript is ready for publication.

Reviewer #2: The authors have performed an excellent job of modifying the manuscript to address the reviewers' comments, in particular in their demonstration of effective host prediction independently of phylogenetic clustering. It's up to the authors, but in many ways, Fig. S4 is more informative in this respect than Fig. 2 (ie. showing moderately effective host prediction in each genus) and therefore removing one obvious source of compounding error with the approach. I would be happy see Fig. 2 updated if the authors wished.

**********

Part II – Major Issues: Key Experiments Required for Acceptance

Please use this section to detail the key new experiments or modifications of existing experiments that should be absolutely required to validate study conclusions.

Generally, there should be no more than 3 such required experiments or major modifications for a "Major Revision" recommendation. If more than 3 experiments are necessary to validate the study conclusions, then you are encouraged to recommend "Reject".

Reviewer #1: (No Response)

Reviewer #2: None

**********

Part III – Minor Issues: Editorial and Data Presentation Modifications

Please use this section for editorial suggestions as well as relatively minor modifications of existing data that would enhance clarity.

Reviewer #1: (No Response)

Reviewer #2: See above - Fig. 2

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

15 Apr 2021

Dear Brierley,

We are delighted to inform you that your manuscript, "Predicting the animal hosts of coronaviruses from compositional biases of spike protein and whole genome sequences through machine learning," has been formally accepted for publication in PLOS Pathogens.

We have now passed your article onto the PLOS Production Department who will complete the rest of the pre-publication process. All authors will receive a confirmation email upon publication.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any scientific or type-setting errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript. Note: Proofs for Front Matter articles (Pearls, Reviews, Opinions, etc...) are generated on a different schedule and may not be made available as quickly.

Soon after your final files are uploaded, the early version of your manuscript, if you opted to have an early version of your article, will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting open-access publishing; we are looking forward to publishing your work in PLOS Pathogens.

Best regards,

Kasturi Haldar

Editor-in-Chief

PLOS Pathogens

​orcid.org/0000-0001-5065-158X

Michael Malim

Editor-in-Chief

PLOS Pathogens

orcid.org/0000-0002-7699-2064