High-resolution phage-host assignment through key proteins using large language models

Domain: biology
Published: 2026-03-20
Editorial outcome: major_revision

Abstract
Viral sequences in diverse environments remain largely uncharacterized, impeding our comprehension of their genetic makeup, biological interactions, and potential applications. This underscores an urgent need for innovative analytical methods. Here, we present the VirHost Hunter framework, which employs phage tails and lysins, bypassing the requirement for full genomes, for efficient and high-resolution host assignment. By harnessing Protein Language Models and Vision Transformers, VirHost Hunter captures protein functional homology despite sequence dissimilarity, significantly boosting prediction accuracy. In the scenario of disease-associated gut bacteria, the calibrated VirHost Hunter surpasses existing methods, doubling phage host assignments, expanding taxonomic reach, and revealing previously uncharacterized phages targeting gut bacteria, including Akkermansia and Prevotella. Therefore, we establish a gut phage lysin database, enabling the synthesis of a lysin that effectively and specifically targets an obesity-promoting bacterium. VirHost Hunter's precision and scalability mark a significant leap forward in virome research and present a promising avenue for microbiome therapies.

Introduction
Virome is a significant component of Earth’s ecosystems and has a profound impact on ecological and human health. In various environments, uncharacterized viral genomes and sequences widely exist due to limitations in current analytical techniques and are referred to as viral dark matter. This concept highlights the need for innovative approaches to uncover and understand these hidden viral entities1. The intricate interplay between bacteria and their viruses - bacteriophages (phages) - has garnered significant attention in recent years, fueled by advances in predictive modeling and therapeutic applications. Identifying the host range of phages is essential in studying phage resistance of bacteria, coevolution of phage-bacteria2, 3, the influence of community context on phage-bacteria systems4, and the role of phages in human health and diseases5, 6. In clinical settings, phages have already been adopted to treat infections caused by drug-resistant bacteria, offering advantages in precision medicine due to their host specificity and minimal disturbance of normal gut flora7, 8. Whereas phages also hold great promise in modulating gut microbiota, their efficacy hinges on the availability of phages targeting gut bacteria, particularly those associated with chronic diseases. Only a handful of phages have been reported targeting gut anaerobes, and it has been implicated that isolating gut phages is arduous9, 10.

Phage lysins have demonstrated effective antimicrobial effects in animal models, food industry, and clinical therapies11-15, presenting broad industrial and medical applications. Lysins often exhibit a broader range of specificity than phages, which typically target at the genus or species level, and are easy to synthesize, making them especially suitable under scenarios where phages are unavailable or where available phages are too host-specific to apply16. However, existing lysin databases mainly focus on clinical pathogens that cause infectious diseases and have not been developed to target the broader gut microbial community17, limiting their application in gut microbiota modulation. Predicting phage hosts and establishing a phage lysin database by leveraging gut phage databases to specifically target gut bacteria, serve as an alternative solution.

Various computational approaches have been developed for predicting phage hosts, falling into two categories: alignment-dependent and alignment-free methods (Table S1). Alignment-dependent methods rely on phage marker genes18-20, phage-host relatedness21, and CRISPR spacers22, 23, but they have limitations such as database size, strong data source dependency, alignment parameters (e.g., sequence similarity thresholds, allowed mismatches/gaps), and applicability only to phages with specific marker genes or CRISPR signals24, 25. Alignment-free approaches utilize phage-bacterium interaction matrices26-37, phage whole genomes29, 38 or sequences of receptor-binding proteins (RBPs)39, 40 to predict phage hosts. Notably, Gonzales et al.40 utilized protein language models (PLMs) for feature extraction from RBPs, highlighting the potential of computational techniques in this area. Protein language models (PLMs) are computational models inspired by techniques developed in the field of natural language processing (NLP). This approach significantly improves contextual understanding and enables the identification of complex patterns that were previously difficult to discern.

Current viral databases predominantly use alignment-dependent methods and CRISPR spacers for host assignment, resulting in incomplete coverage and limited recall. For instance, Paez-Espino et al.41 identified 9,992 putative virus-host associations covering only 7.7% of metagenomic viral contigs (mVCs) in their study of Earth’s virome41. In the past three years, several human gut virome databases have also been released. The Metagenomic Gut Virus (MGV, 2021) database assigned hosts to 81% (n = 153,892) of the phages, followed by 69% (n = 31,259) in the Cenote Human Virome Database (CHVD, 2021), 42% (n = 13,954) in the Gut Virome Database (GVD, 2020), and 29% (n = 40,932) in the Gut Phage Database (GPD, 2021)42-45. The GPD had the most stringent criteria resulting in the lowest recall of the four databases, i.e., it only utilized CRISPR spacers from 2,898 high-quality genomes of cultured human gut bacteria and tolerated zero mismatches across the whole length of the spacers. Therefore, a high-quality alignment-free method, with improved machine learning models and input features, can be complementary to CRISPR spacers method to increase the sensitivity of host prediction without compromising precision.

Indeed, a recent tool, iPHoP46, integrates alignment-dependent and alignment-free methods for host prediction, including Blast47, CRISPRS 22, 23, WIsH26, VHM31, 48, and PHP32. Although iPHoP is a comprehensive framework integrating six approaches for phage host prediction, the authors discussed its limitations, including slow running time and the fact that it only achieves genus-level resolution, which may impact its practical applications. Alignment-free computational methods based on host-specific proteins such as tails and lysins, instead of the whole genomes of phages, offer a complementary strategy. In some contexts, these methods can be as effective or even more practical than whole-genome-based methods. There are several advantages: (1) they require minimal data input, avoiding vast redundant information and overuse of computing resources; (2) they can make predictions based on incomplete genome assemblies resulting from metagenome or metavirome sequencing, as long as key host-interacting proteins, such as tails or lysins, are present; (3) they can achieve high-resolution host prediction, potentially at the species or strain level, for phage therapy applications49-54; and (4) they facilitate applications in synthetic biology, including host range modulation by swapping or engineering phage RBPs50, 55, delivery vehicles based on proteins recognizing and attach to host surfaces56, and therapeutic agents based on lytic proteins breaking down bacterial cell walls57.

In this study, we develop VirHost Hunter, a framework for phage host prediction that combines protein and DNA feature representations, constructs datasets from specific proteins, and leverages Protein Language Models together with a Vision Transformer to enhance prediction performance. We verify the roles of each key component of our design by conducting control analyses, followed by a comprehensive comparison to other methods across family to species levels. We calibrate the model to facilitate its application towards disease-associated gut bacteria and validate its robustness under targeted scenarios. We apply the calibrated model to the GPD and identify a great number of phages targeting disease-associated gut bacteria, including previously uncharacterized phages targeting renowned bacterial species whose phages have hardly been characterized before. To further promote the application of the resource, we extract lysins from the GPD with expanded host assignments to establish a repository. As a proof of concept, we select a lysin from the repository and synthesize it to verify its function against an obesity-promoting bacterium. This work elucidates the design of a predictive framework for phage host prediction and provides insights in how to utilize machine learning to advance genomics data mining and protein function prediction. Deciphering gut phages using this tool not only enhances our understanding in phage diversity and phage-bacteria interactions, but also facilitates the downstream application of gut phage resources into disease intervention.

Results
Designing a phage host prediction framework
To predict the host of phages, full genome sequences of phages and bacteria are usually used. However, whole-genome-based methods introduce a significant amount of nonessential data, including proteins unrelated to host recognition or infection, which can create noise and interfere with the prediction accuracy, resolution, and efficiency. Concentrating on phage proteins conferring specificity—those directly involved in the infection cycle—offers a more targeted approach. In our previous study58, we showed that models using only tail protein sequences achieved performance comparable to those using full viral genomes. Some methods have utilized receptor-binding proteins (RBPs) to predict hosts39, 40, but it can be challenging to annotate RBPs for many phages. We initially counted the number of RBPs and tail proteins in 7,598 phage genomes from NCBI (December 29, 2021), revealing an average of 1.33 RBPs and 15.24 tail proteins per phage (Figure S1). Additionally, other tail components also mediate specific phage-host interactions and are structurally associated with RBPs, collectively influencing host recognition specificity59. Therefore, we expanded the dataset to include specific proteins beyond tail fibers and tail spikes: non-RBPs of phage tails, such as tail sheath, tail tube, baseplate, and tail collar proteins; and lysins, which are enzymes highly active against bacterial cell wall60. These proteins are key to the infection cycle while more widely annotated and are thus included for prediction as well.

Proteins may share low sequence similarity while still performing similar functions across diverse species, rendering traditional sequence alignment methods less effective in capturing these functional similarities. To overcome the challenge of predicting host specificity using these proteins, particularly when sequence similarity is low, we employed protein language models (PLMs). ProtT5 was chosen because of its superior performance in functional prediction tasks and exceptional ability in handling low-homology sequences, as demonstrated by Gonzales et al. who compared multiple PLMs40. PLMs provide a powerful solution by learning deep contextual and functional patterns within protein sequences, enabling them to capture viral protein function and viral biology even in cases of minimal sequence homology61.

Because the same protein sequence can be encoded by different DNA sequences, we incorporate DNA sequence features of tail proteins and lysins into the framework. DNA sequences provide additional insights into phages’ genomic context, such as codon usage bias, GC content, and nucleotide frequency, which can further refine predictions by accounting for genomic stability and evolutionary constraints39. To uncover long-range dependencies and global patterns in DNA sequence data, we utilized a Vision Transformer (ViT)62, 63, a language model capable of handling variable-length sequences. As a language model, ViT63, through its self-attention mechanism, captures intricate patterns and contextual relationships within DNA sequences, such as nucleotide interactions that influence genetic structures and functions. These patterns can reveal insights into genetic structures, functions, or relationships that are not easily discerned by examining individual sequences alone.

As a result, we present the VirHost Hunter framework with the above characteristics (Figure 1A). It consists of two primary components: a feature extractor and a classifier. The feature extractor integrates protein sequence embeddings from the ProtT5 model64, along with physicochemical and k-mer features derived from DNA sequences via a deep neural network (DNN)65 comprising a three-path convolutional neural network (CNN; Figure S2) and ViT for multi-scale feature extraction. By jointly leveraging protein and DNA representations, the VirHost Hunter implements a multimodal approach for robust host prediction. The final classification step uses a multi-layer perceptron (MLP) and a Random Forest (RF) classifier66, with RF refining high-confidence predictions from the MLP to improve accuracy. We next investigate whether combining protein and DNA features, constructing datasets from specific proteins, and using language models such as PLMs and ViT enhances host prediction as expected, respectively.

Using both protein and DNA features improves learning over either alone
To evaluate whether integrating protein and DNA features offers superior performance compared to using either individually, we conducted ablation experiments using two benchmark datasets: the Bacteriophage RBP (Drug-Resistant receptor-binding proteins, DRRBP) dataset (n=4,845)39 and the Bacteriophage Tail Proteins (Drug-Resistant tail, DRTail) dataset (n=12,509). We measured performance using accuracy (ACC), precision, and F1 scores under three experimental conditions: using only protein features, using only DNA features, and using a combination of both. It is shown that relying on a single type of feature led to inconsistent model performance across different datasets (Figure 2A). Particularly, in the DRRBP dataset, models that used only protein features outperformed those that used only DNA features. In contrast, for the DRTail dataset, DNA features alone provided better performance than protein features. This inconsistency reveals the limitations of using only one feature type, as neither approach fully captures the complexity of phage-host interactions.

On the other hand, integrating both protein and DNA features consistently improved model performance across all datasets and metrics. For instance, using both feature sets together resulted in the highest performance, with an accuracy of 0.9081 and 0.8927, precision of 0.9090 and 0.8930, and F1 scores of 0.9077 and 0.8925 on the DRRBP and DRTail datasets, respectively, significantly outperforming models that used either feature set alone. This demonstrates that integrating protein and DNA features not only enhances predictive accuracy but also provides greater consistency and stability across datasets, particularly in the context of bacteriophage host prediction.

Phage tail components and lysins drive host prediction without full-genome data
To confirm that using all tail components - RBPs, tail sheath, tail tube, baseplate, tail collar, etc - for host prediction is feasible, we conducted a 10-fold cross-validation on the DRRBP and DRTail datasets using our method, DeepHost (CNN-based DNA feature processing)29, Boeckaerts et al.'s method (RF with handcrafted features)39, and M. Gonzales et al.'s method (ProtT5 protein embeddings with RF)40 (Table S2). The results suggest that host prediction accuracy through all tail proteins is comparable to that via RBPs, emphasizing the utility of incorporating all tail proteins since they exist 10 times more than RBPs. Additionally, VirHost Hunter outperformed the other methods, achieving an accuracy of 0.9081 and 0.8927, precision of 0.9090 and 0.8930, and F1 scores of 0.9077 and 0.8925 on the DRRBP and DRTail datasets, respectively (Table S2).

To compare the efficacy of phage tails and lysins with that of non-specific proteins for host prediction, we conducted evalution using head proteins and terminases from the same phage datasets. Phage tails and lysins consistently outperformed head proteins and terminases across family, genus, and species levels (Figure 2B and Figure S3). Across all tested sequence similarity thresholds, phage tails achieved the highest accuracy, followed closely by lysins. In contrast, head proteins and terminases reached significantly lower accuracy, with a notable decline in performance at lower sequence similarity thresholds, especially at the species level (Figure S3). This highlights that phage tails and lysins maintain their predictive power, even at reduced sequence similarity, unlike the non-specific control proteins. We further demonstrated that VirHost Hunter can reach species-level resolution and outperformed existing methods in accuracy, precision, and F1 score on a multi-taxonomic dataset of 7,598 phage genomes (Figure 1B, Supplementary Results, Figure S4A and Table S3-S4). As shown, the training loss decreases steadily and converges to near zero, while the validation loss stabilizes without divergence, indicating that the model achieves good generalization and does not suffer from overfitting (Figure S5).

Functional homology is captured even in low-similarity sequences
To demonstrate the ability of VirHost Hunter to capture functional homology using protein language model (PLM) embeddings and DNA features extracted via a Vision Transformer (ViT), we evaluated its performance across datasets with varying sequence similarities. Using CD-HIT67, we partitioned the multi-taxonomic dataset into subsets with sequence similarity thresholds of 50%, 60%, 70%, 80%, and 90%, enabling us to assess VirHost Hunter’s capability of predicting phage host based on functional relationships rather than strict sequence homology.

We benchmarked VirHost Hunter against existing models, including Boeckaerts et al.39, DeepHost29, and M. Gonzales et al.40, across different similarity thresholds. VirHost Hunter consistently outperformed the other methods across taxonomic ranks—family, genus, and species—highlighting its superior ability to leverage functional homology in low-similarity datasets (Figure 2C and Figure S6).

Crucially, as sequence similarity decreased, the performance gap between VirHost Hunter and the other methods widened, particularly at the family and genus levels. This underscores the increasing importance of capturing functional homology in low-similarity regions, where conventional sequence similarity-based methods typically fail. VirHost Hunter’s integration of protein language models (PLMs), such as ProtT5, and DNA sequence features enable it to move beyond reliance on sequence similarity alone. Instead, it identifies deeper functional relationships, resulting in robust and accurate predictions, even under low similarity conditions.

Robust phage host prediction for application under targeted scenarios
Given the substantial impact of gut bacteria on human health, such as inflammatory bowel disease (IBD)68, 69, colorectal cancer70-72, and metabolic diseases73-77, obtaining more information of phages targeting these bacteria is advantageous. We can expand our knowledge in gut phage-bacteria interactions, gut phage diversity, and facilitate their therapeutic application. Phage information can be obtained either through co-culturing with bacterial host or mining data from high-throughput sequencing. However, gut phages, especially those targeting obligate anaerobes, are hard to culture and isolate. Investigating gut phages by analyzing sequencing data is, therefore, usually considered more efficient. Having validated the superior performance of VirHost Hunter, including accuracy, precision, and resolution, we next evaluate its effectiveness in identifying phages targeting disease-associated gut bacteria. We compiled a dataset consisting of 60 gut bacterial species associated with various diseases, including carotid atherosclerosis, inflammatory bowel disease (IBD), and obesity (Figure 1B, Supplementary Data 1). We annotated prophage tails and lysins from the dataset, yielding a total of unique 328,701 tail proteins and 312,565 lysins. We calibrated VirHost Hunter model using these sequences across 29 families, 40 genera, and 60 species.

Consistent with previous evaluations, VirHost Hunter outperformed the other three tested methods when applied to this dataset. At the family, genus, and species levels, VirHost Hunter-tail (based on gut phage tails) yielded ACC scores of 0.9516, 0.935, and 0.9132, Precision scores of 0.9513, 0.9341, and 0.9112, and F1 scores of 0.9512, 0.9342, and 0.9115, respectively (Figure S4B , Table S5). VirHost Hunter-lysin (based on gut phage lysins) exhibited ACC scores of 0.9817, 0.9756, and 0.9590, Precision scores of 0.9817, 0.9755, and 0.958, and F1 scores of 0.9817, 0.9755, and 0.9582, respectively (Figure S4B, Table S6). We also examined how sample sequence similarity would affect model performance. VirHost Hunter consistently outperformed other methods across all similarity thresholds and taxonomic levels (Figure S7), further highlighting its robustness in predicting gut phage hosts associated with chronic diseases.

To further validate VirHost Hunter’s performance on isolated gut phages, we used a previously reported collection of cultivated gut phages78 targeting Bacteroides, Bifidobacterium, Clostridium, Dorea, Eggerthella, Enterococcus, Phocaeicola, Streptococcus, and Lactococcus (Figure 1B). 702 tail proteins and 373 lysins were extracted from 156 gut phages, all with experimentally verified host data. To ensure comparability with existing CRISPR-based host prediction benchmarks, we adopted the same precision thresholds reported by Dion et al.22. In their study, allowing different numbers of mismatches led to trade-offs between precision and recall: E-value of 10-9 led to the highest precision of 95% and the lowest recall of 2.5%, zero mismatch resulted in 84% precision and 31% recall, while allowing two mismatches yielded 69% precision and 49% recall. These thresholds were then chosen for evaluation and comparison. At a 95% precision cutoff, VirHost Hunter correctly identified hosts for 73/156 phages at the family level and 58/156 at the genus level, while CRISPR-based method yielded no assignments due to its low recall rate (Table 1). At 84% and 69% cutoffs, VirHost Hunter performed comparably with the CRISPR-based method, and combining both methods further improved the accuracy to 101/156 (84% cutoff) ,113/156 (69% cutoff) at the family level and 107/156 (84% cutoff), 117/156 (69% cutoff) at the genus level (Table 1). At the species level, VirHost Hunter correctly identified 9/63, outperforming CRISPR, which yielded no valid assignments (95% cutoff). At a 84% cutoff, VirHost Hunter and CRISPR achieved 20/63 and 28/63 correct predictions, respectively, showing similar performance. At a 69% cutoff, CRISPR outperformed VirHost Hunter with 32/63 versus 23/63 correct predictions. Combining both methods improved prediction across cutoffs, with 36/63 correct predictions at a 84% cutoff and 43/63 at a 69% cutoff compared to either method alone (Table 1, Table S11), highlighting their complementary strengths.

To sum up, VirHost Hunter demonstrated superior performance in comparison to the other three alignment-free methods tested. Furthermore, it significantly outperformed the CRISPR-based method in an independent isolated gut phages dataset at a 95% precision cutoff and achieved comparable performance at 84% and 69% precision cutoffs. Additionally, our experiment revealed that the combination of VirHost Hunter and the CRISPR-based method significantly enhances the proportion of true positive predictions, particularly for high-resolution phage-host predictions in the gut microbiota at the species level. Overall, these results highlight the scalability of VirHost Hunter across different environments.

The repository of phages targeting disease-associated gut bacteria is vastly expanded
The four most recently published gut virus databases typically adopted commensal bacteria as their CRISPR libraries. Among them, the GVD42, the MGV44, and the CHVD45 employed loose cutoffs compared to the GPD43, which allowed zero mismatch and resulted in the lowest host assignment rate. Although these databases provide an unbiased representation of gut viral diversity, a tailored approach is necessary for specific application scenarios, such as for intestinal pathobionts. Considering that the GPD had the lowest host assignment recall of 28.66% among the four databases, and as evaluated by Dion et al.22 the precision was 84% at the genus level, we used VirHost Hunter to assign hosts for GPD at both 95% and 84% precision thresholds, aiming to explore the dark matter in the human gut associated with chronic diseases (Figure 1C).

Using our customized annotation pipeline, we identified 163,590 lysins and 388,894 tail proteins from 142,809 assembled gut phages in the GPD. We applied precision filters of 84% and 95% to predict hosts at different taxonomic levels. Based on phylogenetic composition analysis, the prediction results covered 8 phyla, 13 classes, 21 orders, 29 families, 40 genera, and 58 species, including 42 species of obligate anaerobes (Figure 3, Supplementary Data 4). The host assignment for each phage integrated the predictions from both tail and lysin proteins. Notably, 7 families could only be assigned by VirHost Hunter-lysin but not VirHost Hunter-tail, including Eubacteriaceae, Atopobiaceae, Leuconostocaceae, Prevotellaceae, Peptoniphilaceae, Gemellaceae, and Aerococcaceae (Figure 3). We evaluated the host assignment results of VirHost Hunter using 95% and 84% precision and compared that with the previous results of the GPD. We found that both VirHost Hunter-tail or VirHost Hunter-lysin can enhance the host assignment of gut phages. At 95% precision, VirHost Hunter newly assigned host to 15.91% (22,724/142,809) of the GPD phages, with 10.98% (15,677/142,809) by VirHost Hunter-tail and 9.41% (13,432/142,809) by VirHost Hunter-lysin (Figure 4A). At 84% precision, VirHost Hunter newly assigned host to 33.99% (48,545/142,809) of the GPD phages, with VirHost Hunter-tail contributing 20.16% (28,790/142,809) and VirHost Hunter-lysin contributing 25.37% (36,236/142,809), boosting the final host assignment ratio to 62.66% (89,478/142,809) (Figure 4A). These data illustrate that excelling VirHost Hunter on either tails or lysins can enhance the host assignment of gut phages, while combing the results of VirHost Hunter based on different key proteins and that of the CRISPR method could optimize the outcome.

By integrating VirHost Hunter and the CRISPR-based method, we assessed the improvement and refinement of host assignment results in the GPD. Both VirHost Hunter-tail and VirHost Hunter-lysin significantly enhanced host taxonomic classification compared to the previous results. The host assignment results of VirHost Hunter-tail newly covered 3 families, 8 genera, and 20 species and those of VirHost Hunter-lysin covered additional 5 families, 12 genera, and 25 species (Figure 4B-D). Overall, at the family level, VirHost Hunter identified phages targeting 5 previously uncharacterized families accounting for 1.38% of total assignments under 84% precision, while phages targeting Aerococcaceae were not detected at the 95% cutoff. Lachnospiraceae and Bacteroidaceae, recognized as the two most prevalent host families by both VirHost Hunter and the CRISPR-based method, collectively accounted for over 50% of total assignments at both the 84% and 95% cutoffs (Figure 4B). At the genus level, VirHost Hunter identified phages targeting 12 previously uncharacterized genera accounting for 21.58% of total assignments at the 84% cutoff, while three of the previously uncharacterized genera were not detected at the 95% cutoff. Bacteroides is the most abundant host genus identified by both VirHost Hunter and the CRISPR-based method (Figure 4C). At the species level, VirHost Hunter identified phages targeting 25 previously uncharacterized species accounting for 0.14% of total assignments at the 84% cutoff, while four of the previously uncharacterized species were not predicted at the 95% cutoff. Notably, VirHost Hunter identified phages targeting Cronobacter sakazakii as predominant, which was not detected by the CRISPR-based method, likely due to differences in training datasets (Figure 4D).

In the refined database, there are five newly annotated host families, including Aerococcaceae, Akkermansiaceae, Gemellaceae, Prevotellaceae, and Xanthomonadaceae (Figure 1D, Figure 3 and Figure 4B). Phages have potential in diverse applications, including lysing specific bacteria to inhibit pathogens, and integrating into host genomes for in situ microbiota manipulation. Akkermansia muciniphila is commonly regarded as a beneficial gut bacterium, while has also showed potential associations with certain autoimmune diseases and neurological disorders, such as schizophrenia79. Identifying phages that target this bacterium may be utilized for developing undescribed gene-editing and phage display tools, and for modulating the gut microbiota to intervene in disease development. However, phages targeting Akkermansia muciniphila have never been characterized by any previous publications. We successfully identified 36 phages targeting Akkermansia muciniphila at 95% precision cutoff and 95 phages at 84% precision cutoff, and we examined the 36 phages at the more stringent cutoff. The genome sizes of the Akkermansia muciniphila phages range from 11,830 to 92,135 bp and the GC content ranged from 49.11% to 60% (Figure S8). The number of CDS is between 21 and 127 and the annotation rate is between 23.81% and 52.17% (Figure S8). Prevotella copri within the Prevotellaceae family, is another renowned species associated with rheumatoid arthritis80-82 and type 2 and type 1 diabetes mellitus83-85. Megaphages were the only phages reported to target the Prevotella copri, but previous attempts to isolate them failed86. We successfully identified 15 phages targeting Prevotella copri at a 95% cutoff and 22 phages at an 84% cutoff, and we examined the 15 phages at the most stringent cutoff. It was shown that the genome sizes of the Prevotella copri phages range from 12,114 to 127,100 bp and the GC content range from 39.17% to 48% (Figure S8). The number of CDS is between 16 and 166 and the annotation rate is between 23.17% and 56.25% (Figure S8). We selected representative phages targeting Akkermansia muciniphila and Prevotella copri using CD-hit with thresholds of coverage of 0.6 and identity of 0.6 and annotated their genomes using our refined pipeline. It was shown that the functional elements of phages mainly include lysis, lysogenic-related, structure, DNA maintenance, packaging and assembly, replication and transcription, transport, and regulation (Figure 4E).

Diversity and geographic distribution of gut phages
To gain further insights from the expanded host assignments, we first analyzed the phylogenetic lineages of phages in the refined database (Figure 1D). Out of the 89,478 phages, 11.42% were classified under six viral families, including Siphoviridae, Myoviridae, Podoviridae, Herelleviridae, Tectiviridae, and Microviridae, covering all taxonomic classifications identified in the GPD (Figure 5A). The remaining 88.57% of assigned phages were unclassified (Figure 5A). Compared to the previous results of the GPD, VirHost Hunter increased host assignments by nearly 1-fold to Siphoviridae and Herelleviridae phages, approximately 2-fold to Podoviridae and Myoviridae phages, and 33.3% to Microviridae phages, significantly enhancing host assignments across multiple taxonomic levels. As a result, we assigned hosts to 79,332 unclassified phages, 4,566 Siphoviridae phages, 2,902 Myoviridae phages, 2,598 Podoviridae phages, 75 Herelleviridae phages, 4 Microviridae phages and 1 Tectiviridae phage (Figure 5A). It is noteworthy that Microviridae, a class of tailless phages, were assigned hosts by VirHost Hunter-lysin instead of VirHost Hunter-tail as expected. Therefore, it is important to combine the results of VirHost Hunter-tail and VirHost Hunter-lysin for downstream analyses. These findings demonstrate the broad applicability of VirHost Hunter for host prediction across diverse phage lineages, regardless that the phages are with or without tails.

Given the large number of phages predicted to target identical hosts, we evaluated phage diversity within bacterial families across diverse phyla by calculating the ratio of viral cluster numbers to phage counts sharing the same host family, with a higher value indicating greater viral diversity. We observed a wide distribution of phage diversity across bacterial families, especially within Bacteroidetes. Notably, 23 bacterial families exhibited the highest viral diversity (with a value of 1), indicating that each viral cluster within these host families contained only one viral sequence. Among these, 15 families belong to the phylum Firmicutes (Figure 5B), which is consistent with the GPD database showing that Firmicutes harbor significantly higher viral diversity. Additionally, eight other families, including Pseudomonadaceae, Neisseriaceae, Muribaculaceae, Moraxellaceae, Dermabacteriaceae, Corynebacteriaceae, Coprobacteriaceae, and Cellulomonadaceae, which belong to the phyla Actinobacteria, Bacteroidetes, and Proteobacteria, also demonstrated the highest viral diversity (Figure 5B). In contrast, the top 5 families with the lowest viral diversity all belong to the phylum Bacteroidetes, including Bacteroidaceae, DTU089, Marinifilaceae, Rikenellaceae, and Tannerellaceae (Figure 5B).

To gain insights into the relationship between the viral cluster number of host families and their geographic distribution, we analyzed the dominant families and performed principal coordinate analysis (PCoA) (Figure 1D and Figure 5C-D). The results showed that Asia and Europe have a higher total phage count compared to other regions, which may be attributed to the greater number of human metagenomic sequencing studies conducted in Asia and Europe (Figure 5C). Host families show similar geographic distribution patterns across continents, with Lachnospiraceae and Bacteroidaceae dominating across all continents, indicating their role as hosts for globally prevalent gut phages (Figure 5C). Similarly, Asia, Europe and Africa have more overlapping regions, suggesting similar phage compositions, while Oceania, North America and South America are more distinct, indicating different phage communities (Figure 5C). The PCoA result reveals a similar pattern, showing that while some continents do not completely overlap with others in host bacterial community composition, the slight differences observed are not statistically significant (Pr(>F) = 0.065), indicating similar phage compositions (Figure 5D).

An expansive lysin repository countering a broad array of gut bacteria
Considering VirHost Hunter's precision in predicting gut phage and lysin hosts and its complementarity with the CRISPR-based method (Figure 4A, Supplementary Results), we established the Gut Phage Lysin Database (GPLD), which encompasses 117,698 lysins precisely targeting 29 disease-related gut bacterial families (Figure 1E, Supplementary Data 3). Of these, 35.20% (n=41,429) can be identified through both VirHost Hunter and the CRISPR-based method, 13.27% (n=15,617) were identified using the CRISPR-based method, and 51.53% (n=60,652) were exclusively identified by VirHost Hunter. Hydrolases, holins, and endolysins were the predominant functional categories (Figure 6A). To better understand the functionality, stability, and potential applications of lysin proteins, we conducted various analyses focusing on their physicochemical properties. Their secondary structures were composed predominantly of turn, sheet, and helix (Figure 6B), varying in length from 30 bp to 5,811 bp, with a mean length of 195 bp. The molecular weight ranged from 2.8 kDa to 65 kDa, with a mean of 21.6 kDa (Figure 6C), suggesting favorable attributes for efficient synthesis and manipulation. A majority (81.98%) were stable with an instability index below 40 (Figure 6C). Amino acid frequency analysis indicated a prevalence of hydrophobic alanine, leucine, and isoleucine, potentially enhancing protein stability and function.

To analyze the functional diversity and sequence-function relationships within the gut lysin protein family, we employed the sequence similarity network (SSN) tool87, 88, which provided insights into their sequence and functional divergence. Lysins were grouped into 603 clusters based on sequence similarity, with nodes colored according to their host taxonomical phylum (Figure 7A). The SSN was predominantly populated by proteins from Firmicutes (n=4,559), followed by Proteobacteria (n=996), Bacteroidetes (n=906), and Actinobacteria (n=600) (Figure 7A). The protein clusters in the SSN were categorized according to their respective protein types, indicating that proteins with similar functions might possess conserved domains that confer these functions (Figure 7A). Notably, holins, differing from other protein types in the network, exhibit high diversity in their sequences and structures (Figure 7A). Furthermore, proteins targeting the same host phylum tend to cluster together, implying that lytic proteins, much like phages, exhibit a high degree of host specificity. This host specificity is likely mediated by conserved regions within the lytic proteins that are essential for identifying and binding to host cells.

To uncover the conserved functional motifs and the underlying mechanisms, we generated sequence logos89 for three representative clusters. Cluster 1, the largest cluster, contained 1,735 protein sequences targeting Firmicutes, Actinobacteria, Bacteroidetes, and Proteobacteria, including holins (n=3,170), hydrolases (n=3,098), endolysins (n=1,125) and lysis proteins (n=796). The sequence logo analysis revealed three conserved motifs (RHTKAPAVLIECCFVDNKDD, NVTVHRDFANKSCPG, and RSWCSSSAANDNRAITIEVA), all located in the N-acetylmuramoyl-L-alanine amidase domain, which is crucial for phage-mediated bacterial lysis (Figure 7B). Cluster 2 comprised 1,213 representative holins sequences, and two conserved motifs were identified in the toxin secretion domain, which could facilitate the release of lytic enzymes to lyse bacterial cells (Figure 7B). Cluster 3 consisted of 296 hydrolase representative sequences, and its motif was mainly associated with the N-(deoxy)ribosyltransferase-like domain, which functions in degrading bacterial cell walls during phage infection (Figure 7B). Functionally important residues were conserved within putative isofunctional clusters, while motif and domain analyses revealed structural and functional differences between different types of lysin proteins. The findings have valuable implications for the design and engineering of lysins and their application in lysin therapy.

Lysin Ply491_6 effectively and specifically inhibits an obesity-promoting bacterium
Obesity has emerged as a significant global health concern, with the gut microbiome implicated in its onset and progression90 . Comparative analyses have revealed distinct microbiome profiles between obese and non-obese individuals, suggesting association between certain bacterial genera and obesity, including Bacteroides, Megamonas, Ruminococcus, Dorea, Coprococcus, Fusobacterium, Blautia, and Eubacterium73, 91, 92. While phage therapy holds promise for modulating the gut microbiota, the lack of reported phages targeting Megamonas and our failure in repetitive attempts to isolate Megamonas phages prompted our investigation into the therapeutic potential of lysins from the Gut Phage Lysin Database (GPLD) against this bacterial genus (Figure 1F).

We identified 505 unique lysin sequences specific to Megamonas from GPLD, which were clustered into 167 distinct clusters (Figure 8A). Among these, Ply491_6 (ivig_491_6) is the representative sequence of the protein cluster with the highest number of proteins (Figure 8A and Figure S9). The cDNA sequence encoding Ply491_6 spans 561 base pairs. Ply491_6 comprises 187 amino acids, with a molecular weight of 20.8 kDa and a theoretical isoelectric point (pI) of 5.37. Ply491_6 exhibits hydrophilicity, with a grand average of hydropathicity (GRAVY) value of -0.207. The instability index is 25.93, indicating that Ply491_6 is a stable protein. Additionally, Ply491_6 is devoid of signal peptides and transmembrane regions and is structurally characterized by four predominant α-helices alongside multiple β-sheets (Figure 8B). Ply491_6 shares high sequence identity (99.46%) with QIW89318.1, a cell wall hydrolase autolysin from Caudoviricetes sp., and contains a conserved N-acetylmuramoyl-L-alanine amidase domain.

Therefore, we synthesized and purified the Ply491_6 protein for in vitro assays to verify its lytic activity against Megamonas. We incubated Ply491_6 with Megamonas rupellensis and monitored the bacterial turbidity over time. It was shown that Ply491_6 effectively lysed bacterial cells at concentrations as low as 20 µg/mL, with a significant reduction in bacterial turbidity observed within 150 minutes (Figure 8C). To further assess the specificity of Ply491_6, we measured its lytic activity against other high-abundance gut bacteria and common probiotics, including Bacteroides fragilis, Clostridium perfringerns, Ruminococcus gnavus, Bifidobacterium longum, Lacticaseibacillus paracasei, and Lactiplantibacillus plantarum. Ply491_6 demonstrated minimal impact on the viability of these bacteria (Figure 8D-E). These results underscore the efficacy and specificity of Ply491_6 against Megamonas, positioning it as a promising candidate for targeted bactericidal therapy against obesity-associated dysbiosis. These findings provide valuable insights into phage-bacteria interactions in the gut and offer essential data to support the development of precision therapies targeting intestinal pathobionts.

Methods
Establishment of the VirHost Hunter framework
VirHost Hunter consists of two primary components: a feature extractor and a classifier (Figure 1A). When extracting features for VirHost Hunter, we utilize three distinct tools to process phage specific proteins and their corresponding DNA sequences, resulting in three types of features: protein sequence embeddings from the pre-trained ProtT5 model, physical-chemical characteristics of DNA sequences, and k-mer features of DNA sequences extracted via a DNN network. These features will be elaborated upon in detail below.

For protein sequence representation, we leveraged the capabilities of the pre-trained protein language model ProtT5 to generate dense vector representations (embeddings) of protein sequences. Specifically, we utilized only the encoder portion of the ProtT5 model. The encoder integrates essential components such as a multi-head attention mechanism and feedforward layers, enabling it to capture intricate relationships between amino acid residues in the input protein sequence. This process yields rich embedding vectors containing valuable information regarding protein structure and functionality. We extract the average embedding vector from the last layer of the pre-trained model to generate the embedded feature vector, resulting in a 1024-dimensional feature vector for each protein sequence.

The physical-chemical features employed to represent DNA sequences align with the methodology proposed by Boeckaerts et al39. These features encompass nucleotide frequency, GC content, codon frequency, and codon usage bias, amounting to a total of 133 dimensions for the representation of DNA sequences.

To preserve the intrinsic sequence information of DNA sequences, we encoded them following the approach outlined by Wang et al. in their study DeepHost29. This method represents DNA sequences through K-mer frequency. Subsequently, we construct a deep neural network (DNN). The DNN incorporates a convolutional neural network with three paths, each outputting a different number of channels, facilitating the capture of feature information at varying scales (Figure 1B). By using multiple channels in parallel and fusing their outputs, the model can simultaneously learn abstract features at different levels. Subsequently, we leverage the Vision Transformer (ViT)62, 100, utilizing the self-attention mechanism of the ViT model to capture global relationships and multi-channel feature representations, yielding richer original sequence feature embeddings.

Subsequently, we concatenate the three types of features learned by the model into a unified vector, which serves as input for the classifier ensemble. This ensemble includes an MLP neural network, an autoencoder, and a random forest. The training process follows four sequential phases: (1) The MLP and feature extractor are trained end-to-end using Adam optimization with cross-entropy loss; (2) The feature extractor parameters are then frozen ("stabilized") to maintain consistent feature representations; (3) The autoencoder is trained on these fixed features using reconstruction loss; (4) Finally, the random forest is trained on the processed features through an exhaustive grid search that optimizes tree depth, node splitting criteria, and ensemble size to maximize prediction accuracy. Throughout this process, the MLP neural network updates parameters while the feature extractor remains fixed after initial training. We employ the softmax function as the activation function, cross-entropy as the loss function, and utilize the Adam algorithm to optimize the loss function.

Due to the inherent characteristics of the softmax function, which may produce overconfident predictions in incorrect categories, we implement a blended prediction approach combining MLP and random forest outputs. This integration occurs through three steps: First, we train the MLP and feature extractor, then stabilize the feature extractor parameters before training the autoencoder and RF. During testing, the model executes a blending operation represented by the “×” symbol in Figure 1A. The MLP generates softmax-normalized class probabilities while the RF produces independent probability estimates. These outputs are linearly combined with equal weighting (0.5) to yield the final prediction according to the equation:
𝐹𝑖𝑛𝑎𝑙 𝑃𝑟𝑒𝑑𝑖𝑐𝑡𝑖𝑜𝑛 = 0.5×𝑀𝐿𝑃 𝑂𝑢𝑡𝑝𝑢𝑡+0.5×𝑅𝐹 𝑂𝑢𝑡𝑝𝑢𝑡
This design enables the RF to calibrate potential MLP overconfidence while preserving accurate predictions. Implementation details are provided in Supplementary Table S7.

Bioinformatics pipeline for phage genome annotations
A bioinformatics pipeline was developed to enable the rapid and efficient annotation of phage tails and lytic proteins. The pipeline involved several steps. Firstly, proteins predicted from phage genomes using Prodigal v2.6.3 (-f gff -c -p meta)101. Secondly, the predicted proteins were aligned against multiple databases, including 1) the NR phage protein database using Blastp v2.3.0 (-evalue 1e-5 -max_target_seqs 1 -outfmt ‘6 qseqid sseqid stitle pident length mismatch gapopen qstart qend sstart send evalue bitscore’)47, 2) Uniref phage protein database using phmmer v3.1b2 (-E 1e-5), 3) Uniprotkb phage protein database using phmmer v3.1b2 (-E 1e-5) and 4) TIGRFAM, SMART, CDD, ProSiteProfiles, SUPERFAMILY, PRINTS, PANTHER, Gene3D, PIRSF, Pfam, Coils, and MobiDBLite database using hmmscan v3.1b2 (-E 1e-5). BLASTp results were prioritized; if no confident match was found, we used a consensus from phmmer and hmmscan.

To ensure annotation reliability, alignment results were further filtered to retain only those with ≥50% sequence coverage and ≥50% identity, and hits were considered valid only when supported by one or more curated domain databases102-106. The final annotation was generated by integrating results across all tools and databases, emphasizing high-confidence matches107, 108. This conservative strategy prioritizes specificity over sensitivity, which is particularly important when constructing datasets for downstream analyses that rely on accurate functional labels.

Phage tail proteins were identified using the keyword "tail" and lysin proteins were extracted using keywords such as "lysis," "lyase," "lysin," "holin," "hydrolase," "spanin," and "endolysin" from final annotations via Seqkit v0.16.0109.

Construction of benchmark datasets
Complete phage genomes were collected from NCBI using specific keywords related to bacterial hosts, including ‘Staphylococcus’, ’Acinetobacter’, ‘Escherichia’, ‘Clostridium’, ‘Klebsiella’, ‘Pseudomonas’, and ‘Salmonella’. A total of 3,116 phage genomes were collected. Protein annotation was performed using the bioinformatics pipeline, resulting in 22,151 phage tail proteins (21,264 from the pipeline and 887 from a published paper by Boeckaerts et al.39). From these, 7,493 RBPs were screened out using specific keywords related to the tail protein functions, including ‘fiber’, ‘fibre’, and ‘spike’.

Three filters were applied to clean the tail proteins and RBP datasets: 1) sequences with lengths shorter than 50 amino acids or longer than 1,500 amino acids were removed, 2) sequences containing undetermined amino acid ‘X ’in protein sequences or undetermined nucleotides ‘N ’in CDS were excluded, and 3) identical protein sequences with different hosts were discarded to remove redundancy. The final benchmark datasets consisted of 4,845 RBPs in DRRBP and 12,509 tail proteins in DRTail, respectively.

Construction of tail protein and lysin datasets at multi-taxonomic levels
Phage genomes from the viral category were screened in the NCBI Refseq database as of December 29, 2021 (https://www.ncbi.nlm.nih.gov/genome/browse/#!/viruses). Those contain partial genomes and coding sequences were excluded from the dataset. Genomic sequences in FASTA format and annotation files in GBFF format were downloaded from the corresponding table on the NCBI FTP site. This screening process resulted in a total of 7,598 phage genomes for further analysis. Next, information related to the host organism was extracted from the annotation files (GBFF format) using a custom script. If the ‘host=’ filed was empty, the species information mentioned in front of the phage in the GenBank tile (ORGANISM) was selected as the host information. For instance, if the host information of phage AF234172 was empty, we selected ‘Escherichia’ as the host based on the record ‘ORGANISM: Escherichia virus P1’. Then, NCBI taxonomy toolkit, TaxnoKit, was used to obtain the taxonomy ID and taxonomy level of the host organism (taxnokit name2taxid –show-rank). The host taxonomic information was transformed into a standard format including phylum, class, order, family, genus, species, and strain (taxonkit lineage | taxnokit reformat | cut -f 1,3). This process resulted in the compilation of phage-host taxonomic rank information.

We counted the number of RBPs and tail proteins in the datasets and observed an average of 1.33 RBPs and 15.24 tail proteins per phage (Figure S1). We also found that 53.10% of phages lacked RBPs, prompting us to construct multi-taxonomic levels dataset at different taxonomic ranks, enabling the establishment of a tail proteins-based VirHost Hunter (VirHost Hunter-tail) for broader applications. We filtered the phage data and created a phage tail protein dataset, including 37 families, 54 genera, and 57 species. Additionally, we trained VirHost Hunter on lysins – another type of host-specific protein – using 37,469 lysin protein sequences from the same 7,598 phages to construct lysin-based VirHost Hunter (VirHost Hunter-lysin). The lysin dataset comprised 37 families, 42 genera and 47 species.

Family, genus, and species datasets for tail proteins were constructed based on the taxonomic ranks obtained in the previous step. The three datasets were filled using the same three filters as used in constructing the benchmark datasets. Category with fewer than 50 counts in each taxonomic dataset were discarded. After filtering, there were 47 families, 72 genera, and 120 species remaining in the tail protein datasets. These three datasets were used to train the VirHost Hunter-tail model. To address bias issues observed in certain taxa, taxa with precision lower than 0.7 were eliminated from the datasets. For example, Enterobacteriaceae has a recall of 0.9098 and precision of 0.6941 at family level, Escherichia has a recall 0.7125 and precision of 0.3526 at genus level, and Escherichia coli has a recall of 0.7373 and precision of 0.3199 at species level. As a result, the final set taxa include 37 families, 54 genera, and 57 species left. Each dataset was randomly split into training, validation, and testing sets with a proportion of 6:2:2.

The original lysin dataset contained a total of 37,469 protein sequences. The same three filters as previously applied to the benchmark datasets were used to clean the lysin dataset. However, the maximum allowed sequence length was set to 1,000 amino acids since protein sequences with a length of over 1,000 amino acids accounted for less than 2%. The screening process and building procedures for the VirHost Hunter-lysin dataset followed a similar approach to VirHost Hunter-tail. Taxa with precision lower than 0.62 were eliminated from the family, genus, and species taxonomy datasets based on the training results. This step ensured reliable predictions for the remaining taxa. After eliminating low-precision taxa, the final lysin datasets consisted of 37 families, 42 genera, 47 species. Each dataset was randomly split into training, validation, and testing sets with a proportion of 6:2:2.

Gut prophage dataset construction for model calibration
To calibrate the VirHost Hunter model for improved identification of gut phages, particularly those targeting disease-associated gut bacteria, we first identified 60 gut bacterial species associated with diseases. These associations were based on positive correlations from large-scale cohort metagenomic analyses, validation of pathogenicity through animal experiments, or studies on pathogenic mechanisms (Table S5). We then downloaded 13,483 high-quality genomes of these bacterial species from NCBI RefSeq. Prophages were identified within these bacterial genomes using ProphageHunter110, a reliable tool designed to detect prophages with high precision. Subsequently, tail and lysin proteins were identified through stringent sequence alignment, resulting in a dataset of 328,701 unique tail proteins and 312,565 lysins. Since these prophages integrate into the bacterial genome rather than causing immediate lysis, their phage-host relationships are inherently accurate, providing definitive host information. We used this information as the host data for the tail and lysin proteins of prophages to calibrate our model.

Additional filter for higher precision at multi-taxonomic ranks
Since the range of categories that our model can cover is limited, an additional filter was implemented to VirHost Hunter trained in the multi-taxonomic levels’ dataset and the gut prophages dataset, to generate an ‘Unknown ’output for any given input that exceeded the prediction range. To determine the appropriate cutoff for this filter, two datasets were constructed: Positive Control, which comprised samples from the test dataset, and Negative Control containing samples not belonging to any predefined classes in the training dataset. The recall and precision on the Positive Control and the specificity on the Negative Control were illustrated in Figure S10-11. These figures showed that more stringent cutoffs resulted in higher precision and lower recall. This phenomenon occurred because as the cutoff increased, more data were classified as ‘Unknown’, and the remaining data was considered more reliable by VirHost Hunter.

To benchmark VirHost Hunter’s performance against other methods, we considered the work of Dion et al.22 where they evaluated the precision and recall of a CRISPR spacer-based method under different cutoffs of mismatch numbers or e-value at the genus level. They found that with an e-value of 10-9, the method achieved the highest precision of 95% but the lowest recall of 2.5%. With zero mismatches, the method achieved 84% precision and 31% recall. By tolerating two mismatches, the method obtained a balanced performance of 69% precision and 49% recall. Accordingly, several probability cutoffs were selected at the family, genus, and species levels to achieve the same precision values of 95%, 84%, and 69%, respectively (Table S8). Consequently, when the precision on the Positive Control and the specificity on the Negative Control surpassed 95%, VirHost Hunter demonstrated a precision of 95%.

Application of synthetic phage lysin data for host prediction
The latest SQL file (v2021_04) was downloaded from the largest available Phage Lytic Protein database (PhaLP)17. The SQL file provided lysin IDs, corresponding phage genome IDs, lysin annotation information, host taxonomy information, and experimental support information. Phage genome IDs that were not used in VirHost Hunter construction were marked in the dataset collected from NCBI, resulting in 3,448 phage genomes. Lysins that were synthesized and experimentally validated and their corresponding 31 phages were screened out from the dataset. A total of 138 tail proteins were annotated using the custom bioinformatics pipeline. Three phages could not be annotated with tail proteins, leaving final real-world evidence of 31 phage genomes and 138 phage tail protein sequences.

Application of VirHost Hunter in GPD
Gut Phage Database (GPD) and the corresponding taxonomy information table by Camarillo-Guerrero et al. were downloaded43. Based on the ‘Host_range_taxo ’field in the information table, phages were categorized into two groups: those with host information and those without host information. Phage tail and lysin proteins were annotated using a custom annotation pipeline. A total of 163,590 lysin sequences and 388,894 tail protein sequences were obtained from 111,355 phages, which accounted for 77.97% of the total 142,809 phage genomes. A comprehensive dataset of 42,586 proteins was downloaded from NCBI (as of February 22, 2024) using keywords(lysis protein, lysin, lyase, holin, hydrolase, and endolysin AND phage). Lysins encoded by gut phages were identified by comparing with the dataset using BLASTP with a threshold of 60% identity and 50% coverage. VirHost Hunter-lysin model (95%, 84%), VirHost Hunter-tail model (95%, 84%), and CRISPR-based method (84%) were used to construct the Gut Phage Lysin Database (GPLD) targeting human gut commensal bacteria.

Sequence similarity network for the Gut Phage Lysin Database (GPLD)
Biopython was employed to conduct descriptive statistical analysis of the GPLD database, which included aspects such as protein categories, secondary structure proportions, length, amino acid composition, molecular weight, isoelectric point, and stability index. The results were visually represented using ggplot2111. A sequence similarity network was established using a tool developed by Miguel M. Sandin (available at https://github.com/MiguelMSandin/SSNetworks), which was based on lysin sequences clustered using CD-HIT with a 70% similarity threshold. The network construction parameters were set at an identity level of 35% and a coverage level of 50%. The resulting networks were visualized using Cytoscape. Furthermore, MEME was utilized to identify conserved motif sites within the three primary clusters of the network, employing default parameters.

Identification of Megamonas-targeting lysin from GPLD
A total of 505 unique lysin sequences specific to the genus Megamonas were identified from the GPLD. These sequences were subsequently clustered into 167 distinct groups using the CD-HIT with a sequence similarity threshold of 95% and a coverage threshold of 90% (-c 0.95 -aL 0.9). Ply491_6 (ivig_491_6), representing the largest cluster among these groups, was chosen for in-depth characterization and experimental validation of lytic activity. This process involved the prediction of signal peptides using SignalP (https://services.healthtech.dtu.dk/services/SignalP-6.0/), identification of transmembrane regions with HMMTOP (https://services.healthtech.dtu.dk/services/TMHMM-2.0/), and assessment of physicochemical properties via ProtParam (https://web.expasy.org/protparam/).

Synthesis and purification of Ply491_6
To synthesize and purify Ply491_6, the gene encoding Ply491_6 was synthesized and subcloned into the pET-30a(+) plasmid using NdeI and XhoI restriction sites, which was obtained from General Biosystems (China) company. The plasmids were constructed and transformed into BL21 (DE3) competent cells. These transformed cells were cultured on agar plates containing kanamycin at a final concentration of 50 µg/mL at 37°C. Colonies were picked from the plates and cultured until the optical density at 600 nm (OD600) reached 0.6-0.8. Protein expression was induced by adding IPTG to a final concentration of 0.5 mM, followed by incubation of the cultures for an additional 4 hours at 37°C. Cells were then harvested, lysed, and the lysates were subjected by SDS-PAGE to verify protein expression. Then the proteins were purified using Ni-NTA affinity chromatography. The purified proteins were dialyzed into phosphate-buffered saline (PBS) containing 300 mM NaCl, 10% glycerol, and adjusted to pH 7.4, followed by filter sterilization.

Bacterial strains
Megamonas rupellensis strain 150922 was used for lysin activity assay of Ply491_6 in vitro. Bacteroides fragilis bf2 (BF1), B. fragilis bf5 (BF2), Clostridium perfringerns 0840 (CP1), C. perfringerns 0812 (CP1), Ruminococcus gnavus 1177 (RG1), R. gnavus 1186 (RG2), Bifidobacterium longum 4486 (BL1), B. longum 2366 (BL2), Lacticaseibacillus paracasei LAC-F (LP1), L. paracasei LAC-J (LP2) and Lactiplantibacillus plantarum SZHD0015 (L. plantarum) were used for comparing lysin activity of Ply491_6 in vitro. All bacterial strains were isolated from human feces. All bacterial strains were grown overnight in BHI-YH (Brain Heart Infusion medium supplemented with 5 g/L yeast extract, 5 mg/L hemin). To maintain anaerobic conditions, all media and buffers were additionally supplemented with 0.5 g/L L-cysteine hydrochloride and 0.25 g/L anhydrous sodium sulfide, serving as reducing agents.

Lytic activity and specificity of Ply491_6
M. rupellensis strain 150922 was grown overnight, diluted 1:100, and grown to the midlogarithmic phase. The bacterial cells centrifuged, washed, and resuspended in phosphate buffered saline (PBS, pH 7.4) to an OD of 0.9. Phage lysin Ply491_6 was added to bacterial suspension with a final concentration of 20 μg/mL. Each concentration was plated in a U-bottomed 96-well plate in triplicate technical replicates. Ply491_6 dilution plates were then incubated at 37°C in a BioTek Epoch2 Microplate Spectrophotometer(BioTek Instruments, Inc., USA) for 240 minutes. The OD600 was measured every 10 min.

To verify the specificity of Ply491_6, B. fragilis (n=2), C. perfringerns (n=2), R. gnavus (n=2), B. longum (n=2), L. paracasei (n=2) and L. plantarum (n=1) strains were each grown overnight. The bacterial cells were then centrifuged, washed, and resuspended in PBS, and were then incubated with 20 µg/mL Ply491_6 or PBS at 37°C for 240 min. The OD600 was measured every 10 min.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this paper, Du et al. propose VirHost Hunter, which is a machine learning tool that combines language models to predict virus-host associations. They apply their model to the Gut Phage Database and uncover many new phage-bacteria associations. Additionally, using VirHost Hunter, the authors establish a database of gut phage lysin with host assignments. Finally, based on predictions made by VirHost Hunter, they synthesized a lysin and experimentally validated its effectivity in lysing host cells.\nThe work is interesting as it covers a lot of ground from building a machine learning model and a database to providing an experimental proof-of-concept validation of one of their computational predictions.\nMajor comments:\nGiven that the proposed machine learning model is central to this work, it warrants a more detailed description. There are several questions about it that are not touched upon.\n- Why was ProtT5 chosen among other many possible choices.\n- ViT seems to be a transformer model meant for images, and it seems to not have been peer-reviewed yet. Why did the authors choose to use it instead of directly applying a DNA language model.\n- How exactly are the output of the MLP and Random Forest \"integrated\" or \"blended\" ? (related to this, in Fig 1A, what does the x operator do to the outputs of the MLP and random forest)?\n- How is the classifier trained? The text mention updating parameter of the MLP, how about the parameters of the random forest? In L697, what does it mean to \"stabilize\" the parameters of the feature extractor? How are the autoencoder and random forest trained?\n- The input to the VirHost Hunter relies on identifying the key proteins (RBPs, tail proteins, lysins, etc.) based on existing annotation or by using some annotation tools. It is likely that these annotation tools produce false positive annotations or suffer from low recall. Have the authors investigated this matter?\nMinor comments:\n- The title is a bit vague and can be made more explicit. For example, instead of \"NLP-enhanced framework\", \"(biological) language models\" is more informative, and instead of \"decrypting viral dark matter\", the authors can be more detailed.\n-L80 : \"... data source dependency ...\": what does it mean? Isn't it true for any tool?\n-L80 : \"alignment parameters\": how is this a limitation?\n-L85 : \"PLMs are a subset of NLP\" : the former is a mathematical model, the latter refers to a field/domain. please rephrase.\n-L96 : specify date of accession of the databases.\n-L108 : Please provide references to the tools mentioned\n-L109 : \"is the most comprehensive tool to date\" : how so? in what sense?\n-L122 : \"...integrating highlights...\" : English sounds off, please rephrase.\n-L126 and elsewhere: It is not clear what \"calibrate\" means. Does it mean training or evaluation/validation or both?\n-L305: thresholds are mostly arbitrary, but 84% and 69% are curious choices.\n-L329: in what sense are the databases \"comprehensive\"?\n-L337 : \"using our optimized pipeline\": in what sense it it \"optimized\"?\n-L433 : first occurrence of VC. what is it an abbreviation of ?\n-Figure 1B: n = 393,19 and n= 172,00 : are digits missings or the commas misplaced?\n-Fig 3 A: There is too much information in this figure. There are so many items in the Host Family/Genus/Species legend, that it is difficult to match the colors in the legend to the plot. Maybe it might be easier on the eye if the number of rings is reduced to just one precision level? The plot does not have labels 1 through 6, even though the caption refers to it. Also, please fix \"Newly Taxonomy\".\n- Figs 3C -- E : y axis label missing.", "reviewer_2": "Summary:\nDu et al. present an approach for finding phage targets for human gut-associated prokaryotes. Their approach is unique in that it uses only phage tail and lysin sequences to predict whether a phage will infect a host or not. It also uses modern AI and machine learning techniques to optimize these predictions, deploying recent innovations in protein language modeling. They demonstrate their model out-performs CRISPR-based approaches and were able to link many more phages to potential hosts in the human gut compared to previous approaches. They also used their method to identify a lysin for a culture gut microbe and were able to synthesize and effectively kill the host, and that it did not have the same effect on other common gut bacteria.\nStrengths:\n1. The authors do a good job outlining the approach and leveraging modern frameworks to build a tool for phage-host\n2. The authors explain very well how much VirHost Hunter adds to the prediction of hosts for gut phages, and how with this tool we may be able to make better inferences about the phage ecology of the gut, and potential phage therapies for pathogenic hosts from these predictions.\n3. There are many strong phage lysin and tail comparisons in this study, from genomic to physiochemical.\nWeaknesses:\n1. Not demonstrated that all the parts of the methodology are essential/necessary for good predictions. That is, does each additional layer improve the model, or would a simpler CNN + random forest have the same accuracy?\n2. It was not demonstrated that using tails and lysins is as good or better than using whole viral genomes.\n3. Only one lysin used as a proof of concept. Given the high accuracy reported, I’d like to see at least one more example if possible, or any reports where it did NOT work despite the prediction algorithm saying a lysin would work on a host. For people to use this tool, they need to trust it’s going to give them useful information for experimentation.\n4. The viral database to train the model is outdated. The viruses used were identified with VirSorter and VirFinder for the Gut-phage database. I believe using modern tools such as GeNomad and VirSorter2 will vastly improve the quality of the training dataset.\n5. The reasoning for choosing Megamonas doesn’t make a lot of sense. I’m not convinced on the biology that there is an obesity-inducing bacteria, or that killing this microbe in the gut will allow an individual to revert to the non-obese state. I understand if the motivation is just that it’s the only microbe that’s isolated in the lab, but this reasoning doesn’t make sense to me given there are many other pathogens with lots of clinical evidence and motivation behind them this tool could help find therapies for.\nMajor Revisions:\nLine 112-120: The current framing may overstate the power of using lysins and tail proteins for host prediction, as it overlooks the dynamic nature of phage-host co-evolution. Phages and hosts are constantly evolving—which means their interaction isn't static. The genomes of both phages and hosts carry rich, co-evolutionary signals (e.g., shared genes, codon usage, k-mer patterns) that are lost when relying solely on lysins.\nWhile lysins might identify phages capable of attachment at a given point in time, this approach risks missing others that could become infectious as receptors evolve. Therefore, rather than claiming lysins are superior, it would be more accurate to say they hold promise to be as effective—or in some cases more effective—than using the full viral genome, depending on the context.\nLine 225-227: The approach of comparing phage tails and lysins to head proteins and terminases does not prove these components should be used over whole genomes. It would be most appropriate to see phage tails and lysins compared against whole phage genomes, or even include phage tails + lysins as it’s own feature so that per-feature contributions to prediction accuracy are assessed. It’s concluded in 234-235 that whole genomes may be redundant, but it’s only shown that phage tails and lysins are better than heads and terminases, which make up an incredibly small part of the phage genome and is not a fair representation of the whole genome.\nLine 308-312: These precision rates at the species level seem very low compared to the scores reported for the test set in the previous paragraph. Can you explain more why that is? Also why were no other phage-host prediction algorithms compared for the isolate example, since they were compared for the test dataset? iPHoP is not confident at the species level, but can still make predictions from each vOTU to each MAG, and I’m curious how well a tool like this performs on these datasets. CRISPR methods do not feel like a fair comparison since in my experience very few host predictions are based on CRISPR, few hosts carry CRISPR, and the chances of having the phage in CRISPR memory are potentially low, and only indicative of historical infection.\nLine 525: I am not convinced phage therapy will do anything for obesity. I was convinced for pathogens, whereby a microbe colonizes the gut causing a diseased state such as C. difficile or Pseudomonas aeruginosa infection of a wound. These states prevent the return to a healthy microbiome. However, it’s not clear to me that these hosts cause or prevent an obese individual to return to a non-obese state. My general assumption is that everything is everywhere, the environment selects, and by eating certain foods, it recruits certain types of microbes. However, has it been shown that these microbes persist once a healthy diet has been assumed? There just doesn’t seem like enough evidence to support that killing these microbes would return an individual to a non-obese state, or that these microbes are impeding an obese individual from returning to a non-obese state. This may be due to my lack of knowledge, and please correct me here if there are studies with this kind of evidence. I can just imagine many other good motivating diseases associated with the gut that are not obesity, and think that would be a better motivating example if it’s proven the diseased state is caused by the microbes, and killing those microbes would allow a diseased individual to return to the non-diseased state.\nMinor Revisions:\nLine 68-69: Need citations for this statement: “Lysins possess moderate host specificity, are easy to synthesize and are 68 especially suitable under scenarios where phages are unavailable, or the available 69 phages are too host-specific to apply.”\nLine 70-72: Is there a way you can strengthen this gap? The definition of a pathogen and commensal is somewhat subjective, and can vary person-to-person. I’d get rid of “rather than gut commensals” and just say they’re focused on this subset of clinical pathogens, and have not been developed to target the broader gut microbial community.\nLine 100: “resulting the lowest” to “resulting in the lowest”\nLine 113 “the” to “of the”\nLine 115 “they can handle incomplete genome assemblies resulting from virome sequencing” Only if they capture the phage tail\nLine 116-117: Do you have enough supporting evidence to make this claim? What previous studies have shown that phage tails can predict species/strain level? I’d like to see a citation here.\nLine 131: “As a proof of concept, we select a lysin from the repository and synthesize it to verify its function against an obesity-inducing bacterium.” How hard was this to do? I’d like to see a few, especially of cases that you would not have expected.\nLine 152-153: “Therefore, we expanded the dataset to include specific proteins beyond tail fibers and tail spikes: non RBPs of phage tails, such as tail sheath, tail tube, baseplate, and tail collar proteins;”\nYou just explained in the previous sentence that many genes are non-essential. How are these tail sheath, tube, baseplate, and collar proteins essential? Do you have a reason to believe they help the phage confer specificity to the host? Is there evidence to suggest certain tail fibers and tail spikes have associated other tail proteins? Potentially showing some kind of clustering of these proteins that certain tail proteins do indeed cluster/associate with specific tail fibers and spikes would support their inclusion in the model.\nLine 233-236: This is expanding/reiterating my major revision suggestion: You compare head and terminase to phage tail and lysin, which an A compared to B, but a better comparison would be A + B, meaning if you use both these features, is there an improvement in the model, similar to combining DNA and protein level features? Furthermore, I’d like to see this analysis expanded to comparing whole genome (A + B + C), so we could have the base model, which is just tail proteins (A), then tail proteins + capsid and terminases (A + B) and then the whole genome, which already includes A + B, and C is other genes (A + B + C). If we don’t see much improvement over using A alone, then I’d agree with your claim that we don’t need to use whole genomes.\nLine 296-312: Is it ensured that these hosts/tail proteins/lysins are NOT seen by the model already?\nLine 296: These scores seem pretty low compared to the model scores--this makes me think the model is overfitting. What is the sequence space of these phages and it is quite different from the model? Is the model not generalizing well to these phages?\nHow much of the training data for this model is experimentally verified? How are the phage-host associations made? I’m wondering if a smaller training dataset of experimentally verified phage-host pairs would be stronger. There may be a lot of noise or certain type of bias in non-experimentally verified phage-host pairs. CRISPR libraries, for example, may not actually be good indicators of phage-host pairs at the phage tail / lysin level, because they’re indicative of historical infection, but not of what can infect the host right now based on its current makeup. So there’s a lot of noise there. Even if some of the phages can still infect based on CRISPR knowledge, there may also be a lot of noise of phages that can no longer infect. I think i personally need to understand more about how good CRISPRs are at predicting if a phage can infect. Just need to double check that the training set is not composed of these.\nHow is CRISPR not species specific? These CRISPRs probably come from specific strains, and is variable between strains of the same species, so I would have thought CRISPR is species-specific.\nLine 384-386: The claim is Akkermansiaceae modulates multiple diseases, which yes these studies cited do show an association between the microbe and disease, but the papers cited point to higher abundance of Akkermansiaceae decreasing obesity and stabilizing glucose levels in mice (citation 71), and that it is a probiotic in the gut (citation 72). Furthermore, it is indicated in this review, https://www.nature.com/articles/s41575-022-00631-9, that this is a beneficial microbe, with increased abundance associated with decreases in obesity and diabetes. Therefore, it doesn’t make sense why we’d want to phage target this microbe. I’d suggest changing the motivation for why this microbe was chosen.\nLine 417: How does this compare to IMG/VR4? GPD predictions I would consider to be “old” and very conservative since these phages were identified using VirFinder and VirSorter, which are shown to have extremely low recall and precision. I’m curious how this tool stacks up against other host prediction tools like iPHoP. If you give iPHoP only the phage tail and lysin seqeuences, how well does it predict hosts for these? Or if you have the full genomes of these phages, how do iphop predictions compare?\nLin 436: should be Firmicutes, not firminutes.\nLine 437-445: Is more viral diversity a function of host richness? That is, if you have more host genome representatives, is that correlated with having more viral diversity within that taxonomic group?\n“work we only verified the scenario of searching for phages and lysins targeting disease628 associated gut bacteria. Future research should aim to refine VirHost Hunter by 629 incorporating a broader range of datasets, including diverse protein datasets and 630 environmental contexts. A great advantage is that VirHost Hunter only requires input 631 of key proteins, which can be extracted from prophages integrated within bacterial 632 genomes and fragmented phage genomes from metagenomic sequencing, vastly 633 expanding the scale of datasets.” (Du et al., p. 14)\nMy personal opinion is training this model on a more diverse set of phages would vastly improve the model and it’s applicability. I understand there is a large body of work done and eagerness to publish, but I can’t see how there’s anything specific about gut phage-microbe interactions or tail proteins, and that these predictions would likely benefit from including high quality phage-host training data from other sources, such as the curated dataset built in VHIP. I think this is a diverse system just like any other, and teaching the model the “rules” of phage host interactions from any system can only improve it.\nLines 600-700: Was the model built in an iterative fashion? Was each layer justified? That is, does each additional complexity of the model actually improve the model’s ability to make a prediction at the classifier step?\nLine 747-748: What is the integrity of the validation set where this phage infects this host? Are these from experimentally verified experiments? I want to know how good this phage-host pairing method is. So of the 7,598 phages, are these all isolates from experimentally tested infections in these hosts? Or are they from binned metagenome-assembled genomes, or other phage-host prediction tools?\nLine 775: If I understand correctly, E. coli was NOT included in the training set? This is strange to me because E. coli is the most well studied microbe. I’d assume it also has some very well-studied phages. I wonder if this is due to “junk” from NCBI where phages that do not actually infect E. coli are included. Are you able to include only the experimentally verified phages? This may improve the model in other ways, since many host assignments in NCBI are incorrect and not validated. Training on a curated dataset may improve the model, such as VHIP’s phage-host dataset (Bastien 2024, PLOS Computational Biology).\nLine 828: The MGV database is much bigger than the GPD database and with more complete phages. The GPD database used VirFinder and VirSorter to identify phages from metagenomes and isolates. These tools are quite old, and have been shown to be much less accurate than modern tools like VirSorter2 and GeNomad. I’d like to see these phage genomes re- classified with GeNomad, and filtered with an updated CheckV database, to remove any potential non-phage genomes that would impact the model’s precision.\nThe host assignments being made using CRISPR arrays works well, so the host prediction method I’m not worried about the GPD. They also did rigorous stats for prophages binned with host genomes in GPD.\nLine 1194: Why is CRISPR-based zero and N.D. for species?\nFigure 2B: Why would species-level accuracy be higher than family or genus level?\nFigure 4B: Choose different colors perhaps? Hard to read the brown, red, burnt orange.\nFig 6B: can go into the supplemental. This takes up a lot of space and the text is not readable. I see the reason is to point out the placing of the lysin, but it’s not necessary.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Authors’ Response to Reviewers’ Comments\nPaper title: Decrypting viral dark matter through key proteins using large language models\nAuthors:\nZhihua Du"}]

Decision letter