We incorporate three interpretable, domain-knowledge-driven statistical descriptors to complement deep embeddings:

(1) Amino acid composition (AAC): Measures the frequency of each amino acid A_i in sequence S:

where n_i is the count of A_i and L is the sequence length. AAC provides coarse but robust information for classification tasks.

(2) Atomic composition (AC): Represents the proportion of each element E_j (e.g., C, H, N, O, S) in a protein:

where a_j is the number of atoms of E_j in S. AC captures the protein's fundamental chemical makeup.

(3) Molecular weight (MW): Summarizes compositional information into a physicochemical metric:

where M(A_i) denotes the molecular weight of amino acid A_i. MW can discriminate proteins based on mass profiles.

These statistical features enhance model explainability by providing intuitive biochemical insights into protein sequences.

We leverage ProtBERT and ProT5, two pre-trained protein LLMs, to generate context-rich sequence embeddings.

(1) ProtBERT: ProtBERT is a BERT-based protein language model that employs a bidirectional Transformer encoder with a masked language modeling (MLM) objective [20]. During pre-training, a subset M = {m₁, …, m_K} of positions in the input sequence S is randomly masked, and the model learns to predict the masked residues using context. The MLM loss is defined as:

where S_\m denotes S with the residue at position m replaced by a mask token. ProtBERT stacks multiple Transformer layers to iteratively refine the sequence representation:

with H⁽⁰⁾ being the input token embeddings. By leveraging contextual cues, ProtBERT learns both local motifs and long-range dependencies, yielding robust sequence features for downstream tasks.

(2) ProT5: ProT5 is a T5-based encoder–decoder model that frames protein modeling tasks in a sequence-to-sequence format. It can be pre-trained to reconstruct corrupted sequences, predict functional or structural annotations, or even generate novel protein sequences. Given an input sequence S and a target sequence Y = {y₁, …, y_T}, ProT5 models the conditional probability:

It learns by minimizing the negative log-likelihood:

The encoder first produces a hidden representation:

The decoder then uses H_enc (together with previously generated tokens y_<t) to predict the next token y_t:

Through this process, ProT5 captures rich long-range dependencies and excels in generative and multi-task settings. Training on diverse objectives endows ProT5 with a broad latent space of protein sequences, complementing ProtBERT's embeddings.

The datasets originate from a BGI research project titled Research on Artificial Phage Model Construction Based on Deep Generative Adversarial Network Learning (see text footnote ¹), providing a reliable empirical foundation for our prediction tasks. This study is part of an ongoing initiative at BGI-Shenzhen to explore phage-host prediction under real-world data constraints. The equipment is based on the NVIDIA A100 cloud platform. To evaluate the model under multi-source conditions, we utilize two distinct datasets and a third integrated dataset. To evaluate the model under multi-source conditions, we utilize two distinct datasets and a third integrated dataset.

To provide a clearer picture of the dataset composition, we include a heatmap of the host-phage interaction matrix (Figure 1) in the data description. This figure offers an overview of which phages infect which host strains in Dataset 1, highlighting the data's structure before any modeling.

Dataset 1: This dataset includes 101 host bacterial strains and 129 phages, collected under well-defined conditions, ensuring clear information on species composition and environmental context.

Dataset 2: This dataset comprises 38 host strains and 176 phages. Compared to Dataset 1, its sampling conditions (e.g., environment and strain selection) differ, yielding distinct phage–host interaction patterns.

Dataset 3: Dataset 3 combines Dataset 1 and Dataset 2 to enable performance assessment under more diverse conditions, providing a robust test of the model's adaptability to heterogeneous data sources in real-world settings.

(1) Immersion experiment and labeling strategy: The phage–host interactions were measured via immersion experiments, where phages were exposed to host cultures and infection outcomes recorded in a host–phage matrix (Figure 1). We binarized these outcomes by labeling interactions with infection values above 1.5 as 1 (significant infection) and below 1.5 as 0 (no significant infection), ensuring a clear separation between positive and negative samples.

(2) Sequence information and feature extraction: For each host and phage, we obtained the protein amino acid sequence and extracted features including traditional descriptors (AAC, AC, MW) and LLM-based embeddings (ProtBERT, ProT5), yielding a rich feature set for learning.

In summary, these curated multi-source datasets and comprehensive feature sets provide a solid foundation to assess the model's generalization and robustness. Their diversity underscores the importance of MoEPH's dynamic fusion mechanism for reliable predictions across heterogeneous conditions.

The MoEPH model proposed in this study is designed to integrate multi-source features both statistical descriptors and deep sequence embeddings extracted by large language models (LLMs) and to adaptively weight these features through a Mixture-of-Experts (MoE) mechanism. By doing so, the model produces more robust and expressive representations for subsequent classification tasks. Figure 2 provides an overview of the framework.

In A1 (Statistical feature extraction), we derive fundamental statistical descriptors such as amino acid composition (AAC), atomic composition (AC), and molecular weight (MW) from each protein sequence, denoted as Xstats∈ℝN×26 in the example illustration. These features reflect basic physicochemical properties of the protein and serve as an initial numeric reference for subsequent integration.

Concurrently, in A2 and A3 (LLM feature extraction), we employ two pre-trained protein language models: ProtBERT and ProT5 to capture deep semantic representations of the sequences, resulting in XProtBERT∈ℝN×1024 and XProT5∈ℝN×1024. These high-dimensional embeddings encode contextual, local structural, and functional information within the sequences, thereby enriching the predictive power for phage-host interaction tasks.

Next, these three kinds of features are fed into B: the MoE layer, which consists of a Router and two expert modules (Expert 1 and Expert 2). A gating network processes the statistical features X_stats to generate weights α, dictating how X_ProtBERT and X_ProT5 should be combined. This adaptive weighting ensures that the model can dynamically select the most discriminative or best-matched expert features for each sample under varying data distributions and feature patterns. In addition, the resulting weighting coefficients can provide interpretability regarding the relative importance of each LLM for different types of samples.

As depicted in C (Concatenate of features), the fused features obtained from the MoE layer are concatenated with the statistical features, producing Xcombined∈ℝN×1050 (illustrative dimensionality). This consolidated vector representation fully integrates both the “protein statistical attributes” and the “LLM-based deep features.”

Finally, D (Model) processes X_combined using a customizable prediction network (e.g., a CNN, MLP, or RNN) to output the phage-host interaction results. In this study, we employ a CNN classifier to systematically examine how the fused features improve performance (e.g., accuracy, F1-score, AUC) under a fixed network structure.

In summary, MoEPH incorporates statistical features from traditional analyses and advanced embeddings from ProtBERT/ProT5, then uses an adaptive MoE layer to effectively and interpretably combine multi-source information. This integrated solution provides a robust yet flexible approach to phage-host interaction prediction.

In protein-related prediction tasks, relying solely on a single pre-trained model (e.g., ProtBERT or ProT5) often fails to fully capture the diverse sequence patterns and structural information inherent in biological data. To address this limitation, we introduce a Mixture-of-Experts (MoE) mechanism into our model. By combining multiple pre-trained experts and adaptively assigning their importance based on sample-specific statistical attributes, the MoE layer flexibly merges various feature advantages, thereby enhancing predictive performance and overall generalization (Shazeer et al., 2017; Pearce and Zhang, 2021). The entire fusion procedure is summarized in Algorithm 1.

In our Mixture-of-Experts (MoE) layer, a gating network dynamically computes sample-specific weights to fuse the outputs of two expert models (ProtBERT and ProT5). Formally, for each sample i, let XProtBERT,i∈ℝd and XProT5,i∈ℝd denote the pre-trained feature vectors from ProtBERT and ProT5 (with d = 1024 in our implementation). Each sample also has an associated statistical feature vector Xstats,i∈ℝs (with s = 26 descriptors such as physicochemical properties and sequence composition). The gating function g(·) is realized as a two-layer feed-forward network that transforms the statistical features into a pair of importance weights (α_{i, 1}, α_{i, 2}) for the two experts. Concretely, the gating network first applies a fully-connected layer to X_{stats, i} (shared across samples) to produce a hidden representation, then a second fully-connected layer produces two logit scores which are normalized by a softmax activation:

where W1∈ℝs×64,b1∈ℝ64,W2∈ℝ64×2,b2∈ℝ2 are trainable parameters. The ReLU activation in Equation 7 introduces nonlinearity into the gating function, and the softmax in Equation 9 ensures each sample's two gating coefficients (α_{i, 1}, α_{i, 2}) form a probability distribution (i.e., α_{i, 1}, α_{i, 2}≥0 and α_{i, 1}+α_{i, 2} = 1 for each i). Importantly, this gating is dynamic and sample-specific: for each sample i, the statistical attribute vector X_{stats, i} yields its own gating weights α_{i, 1} and α_{i, 2}. This design allows the model to adaptively decide how much to rely on each expert's features based on the characteristics of that sample (rather than using fixed static fusion weights).

Using the gating weights α_{i, 1} and α_{i, 2}, the MoE layer modulates and fuses the expert outputs for each sample. Let α_i = [α_{i, 1}, α_{i, 2}] be the weight vector for sample i. We obtain the fused feature for sample i by an element-wise weighted sum of the two expert feature vectors:

In other words, the ProtBERT embedding is scaled by α_{i, 1} and the ProT5 embedding by α_{i, 2}, and then they are added together to produce a single fused representation for sample i. Stacking these results for all N samples yields the fused feature matrix:

where α:,1,α:,2∈ℝN×1 denote the two columns of α (broadcasted across the d-dimensional feature vectors), and ⊙ denotes element-wise (Hadamard) product. Through Equations 10–12, the gating weights effectively modulate the contribution of each expert: if α_{i, 1}≫α_{i, 2} for a given sample, the fused representation X_{fused, i} will be dominated by ProtBERT's features, whereas if α_{i, 2} is larger, ProT5's features are emphasized. This adaptive fusion flexibly leverages the strengths of both experts, allowing the model to favor the expert that is more informative for each particular sample's attributes.

After obtaining the fused LLM-based features X_fused, we integrate them with the original statistical features. Specifically, we concatenate each sample's fused vector with its statistical descriptor vector to form the final combined feature:

In our implementation d = 1024 and s = 26, so Xcombined∈ℝN×1050. This concatenation preserves the original handcrafted features alongside the fused deep features, ensuring that downstream classifiers receive a comprehensive feature set. The entire gated fusion procedure is summarized in Algorithm 2.

After extracting and preprocessing the features, we employ a convolutional neural network (CNN) as the classification backbone for phage-host interaction prediction, integrating the MoE module for feature fusion throughout the training process. While we fix this CNN architecture for consistency, our main objective is to demonstrate how the proposed MoEPH framework leverages multi-source embeddings to boost predictive performance, the advantage would extend similarly to other model architectures. The overall training and evaluation pipeline is summarized in Algorithm 3.

This study evaluates the proposed model on three datasets (Dataset1, Dataset2, and a merged Dataset3 which combines the former two). Dataset 1 contains phage–host pairs collected under relatively consistent experimental conditions (with well-defined species compositions and environmental factors), whereas Dataset 2 comes from a more complex ecological background, yielding greater heterogeneity in phage–host interaction patterns. By merging these two sources into Dataset 3, we impose more stringent demands on the model's adaptability to heterogeneous, multi-source inputs. This strategy was chosen to preserve each dataset's unique characteristics and to assess whether the model generalizes across different sources. Had we merged the datasets from the beginning, any source-specific patterns or performance differences would be hidden. By first testing on each dataset individually, we can demonstrate MoEPH's robust performance under each condition, and then confirm its adaptability on the merged Dataset 3.

To further challenge the model's robustness, each dataset is examined under three class imbalance settings: the original imbalanced distribution (Raw), an Over-sampling variant, and an Under-sampling variant. In Over-sampling, instances from the minority class are replicated to balance the number of positive and negative samples (e.g., if N_pos and N_neg denote the counts of positive and negative samples with N_pos<N_neg, additional positive instances are randomly duplicated until N′pos≈Nneg). In Under-sampling, the opposite approach is applied: majority-class instances are randomly removed until N′neg≈Npos, thereby equalizing class counts. These three sampling methods simulate varying degrees of class imbalance encountered in real-world scenarios, enabling a comprehensive evaluation of the model's robustness across different data distributions.

In this study, we adopt several widely recognized metrics to evaluate the predictive capability of our model under class-imbalanced conditions (Fawcett, 2006; Davis and Goadrich, 2006; Saito and Rehmsmeier, 2015; Huang and Ling, 2005). Below, we provide their mathematical definitions and briefly discuss their relevance in the context of phage–host interaction prediction.

The proportion of all predictions that are correct.

The proportion of predicted positive cases that are truly positive.

The proportion of actual negative cases correctly identified (true negative rate).

The harmonic mean of precision and recall.

The area under the ROC curve (TPR vs. FPR), summarizing performance across all thresholds.

The area under the precision–recall curve, reflecting the precision–recall trade-off.

Figure 4 provides a comparative visualization of the final feature representations produced by different methods across the three datasets. Each vertical column corresponds to a specific feature channel, and the color intensity indicates the magnitude of the feature values (with darker shades representing higher values and lighter shades indicating lower values). These heatmaps are derived from features output by our MoEPH model after the expert fusion layer, just before classification, under various configurations: using only statistical features, using ProtBERT embeddings, using ProT5 embeddings, and using the fused MoE approach (ProtBERT+ProT5). Compared to the traditional statistical baseline, the MoEPH-based models display more distinct block and stripe patterns, indicating that MoEPH captures finer-grained, context-dependent sequence features than conventional methods. Each subfigure corresponds to one method on one dataset (Figures 4a–l), allowing side-by-side comparison of how feature distributions vary by method and data complexity. The key observations are as follows:

We conducted an in-depth evaluation involving 36 different experimental settings, spanning three datasets, four algorithms, and three data sampling strategies. As shown in Table 1, the MoEPH series—whether incorporating ProtBERT, ProT5, or both—achieved high performance in most scenarios, covering 239 out of 252 metrics (i.e., 94.9%). Its maximum accuracy reached 99.6%, significantly surpassing the current state-of-the-art methods (see Table 2). Moreover, in more complex datasets, the accuracy (ACC) was improved by as much as 31%, rising from 0.535 to 0.845 (Figure 5). In some over-sampled scenarios, certain cross-validation folds with very few positive instances yielded a Sensitivity of 1.0. We note that this perfect recall is due to the extremely low number of positives in those folds, reflecting class imbalance.

Figure 6 illustrates the comparison between our proposed MoEPH variants (ProtBERT, ProT5, and fused) and the Statistical approach, as well as PredPHI (Li et al., 2020), under the Dataset 2 (Under-sampling) scenario. The metrics presented include Accuracy (ACC), Sensitivity (Sens), and Specificity (Spec). We observe that MoEPH (ProtBERT) and MoEPH (ProT5) attain ACC values of 0.84 and 0.845, respectively, outperforming both the Statistical method (0.53) and PredPHI (0.78) by a notable margin. They also exhibit superior Sensitivity, indicating that even under severe under-sampling, the model can still capture a larger portion of positive samples, thus reducing the rate of missed detections. Furthermore, to comprehensively evaluate performance within the same CNN classification architecture, we expand our discussion to include additional metrics, analyzing the results from three key perspectives: varying protein representation algorithms, different sampling strategies, and multiple datasets.

MoEPH holds promising potential to innovate patient care in the context of antibiotic-resistant infections. By rapidly and accurately predicting phage–host interactions, MoEPH could assist clinicians in selecting effective phage therapies tailored to a patient's drug-resistant bacterial infection, exemplifying precision medicine in infectious disease treatment. This approach could be integrated into the clinical workflow as a decision-support tool, where its robust and interpretable predictions provide physicians with high-confidence recommendations for alternative treatments when antibiotics fail. The model's emphasis on interpretability and reliability builds the trust necessary for clinical adoption, ensuring that healthcare providers can understand and rely on its suggestions. Ultimately, a trustworthy AI system like MoEPH could streamline the management of AMR cases—improving treatment outcomes by offering timely, personalized therapeutic options and potentially integrating into hospital infection control and antibiotic stewardship programs.