Prediction of influenza A virus-human protein-protein interactions using XGBoost with continuous and discontinuous amino acids information

Positive dataset

In this study, a total of 10,485 pairs of IAV and human PPIs were acquired by searching the HPIDB3.0 (Host-Pathogen Interaction Database) (Ammari et al., 2016). After applying several filtering criteria, including removing duplicate values, limiting the host to humans, restricting the pathogen to the IAV, setting a minimum protein sequence length of 50 amino acids, and removing non-standard amino acids, a total of 7,011 pairs of IAV and human PPIs were obtained as positive samples for this study. The process is shown in Fig. 2.

Negative dataset

Due to the difficulty in obtaining empirical evidence of non-interaction between two proteins, there is no gold standard for constructing negative samples so far. However, the quality of negative samples affects the accuracy of PPI prediction. In response to this problem, many scholars have proposed methods for constructing negative samples. The mainstream methods for constructing negative samples are random sampling, subcellular localization based sampling, and dissimilarity-based sampling. Random sampling is the earliest mainstream method (Zhang et al., 2018). In this method, human proteins and viral proteins from the real PPI are randomly paired together, the original positive combinations are then removed from these pairs, leaving just the remaining pairs as negative samples. This method may take a significant number of positive samples as negative samples. The approach of constructing negative samples based on subcellular localization has been acknowledged and adopted by numerous researchers since its introduction (Jansen & Gerstein, 2004). Subcellular localization data provides insight into the functional position and biological role of proteins within a cell, making it biologically significant for the creation of negative samples. Nevertheless, this approach has its limitations. Ben-Hur & Noble (2006) experimentally demonstrated that although choosing negative examples as pairs of proteins that are localized to different cellular compartments creates high-quality negative examples, it also makes them easier to distinguish from interacting proteins. Dissimilarity-based sampling is proposed by Eid, ElHefnawi & Heath (2015). They hypothesized that viral proteins with high sequence similarity could theoretically interact with a large number of similar host proteins. That is, if there exist human proteins H1, H2 that interact with viral proteins V1, V2 respectively, and if the sequence similarity between V1 and V2 is greater than 80%, then it can be assumed that H1 interacts with V2, and H2 interacts with V2. Conversely, if the sequence similarity between V1 and V2 is less than 20%, then H1-V2, H2-V1 can be considered as non-interacting and can be used as negative samples. This method greatly reduces false negatives and is more biologically meaningful. This method is currently the most recognized method for constructing negative samples and is still used by many researchers today. In addition, Dey, Chakraborty & Mukhopadhyay (2020) proposed a degree-based negative sampling method, which is based on the principle that human proteins with higher degrees are more likely to be attacked by viral proteins. Here, “degree” refers to the number of connections a protein has with other proteins in the human protein-protein interaction network. The degrees of all human proteins are calculated and sorted in ascending order, and proteins with low degrees are selected as negative samples. In Dey, Chakraborty & Mukhopadhyay (2020) the dataset is single human proteins, and this degree-based method of constructing negative samples has not been applied to protein pairs. The prevalent methods for generating negative PPI datasets continue to be subcellular localization-based sampling and similarity-based sampling.

Inspired by these previous works, here we propose a new method, namely, degree and dissimilarity-based negative sampling. The negative dataset obtained by the dissimilarity-based sampling is then further selected according to the degree, and the human proteins involved in it with lower degrees are selected as the final negative samples, thereby significantly enhancing the stringency of the selection process. This approach enhances the overall rigor of the methodology.

The following is how the degree and dissimilarity-based negative sampling is carried out. First, global alignments of pathogen proteins were carried out using the Needleman-Wunsch algorithm on an all-versus-all basis. These alignments employed a linear gap penalty of 10 and utilized the BLOSUM62 scoring matrix (Peris & Marzal, 2011). The obtained comparison scores were then normalized between 0 and 1. The sequence dissimilarity distance between viral proteins V1 and V2 is obtained by subtracting the normalized scores from 1 and the dissimilarity threshold was set at 0.8. For example, for a certain viral protein, we found all other viral proteins that meet the dissimilarity threshold, obtaining all viral proteins with less than 20% similarity to that viral protein. Then the human proteins that originally interact with these viral proteins in the positive sample were selected and randomly paired with the obtained set of viral proteins. After that, we removed the protein pairs that already exist in the positive sample. The remaining pairs were considered as candidate negative samples. Once all candidate negative samples were acquired, the degrees of the human proteins within those samples were calculated. Finally, the protein pairs were sorted in ascending order according to the degree of human proteins involved. An equal number of protein pairs as in the positive samples were selected to determine the negative sample. As a result, a total of 7,011 protein pairs were designated as negative samples. Figure 3 shows the workflow of degree and dissimilarity-based negative sampling.

Independent dataset

An independent dataset was prepared to predict potential interactions between human proteins and influenza A virus proteins to get potential human target protein. In this study, we downloaded the full set of IAV and human protein sequences from the UniProt/Swiss-Prot (Boutet et al., 2007) database. After excluding non-standard amino acid sequences and setting a minimum protein sequence length of 50 amino acids, which aligns with the positive sample filtering criteria, we eventually obtained a total of 1,384 IAV proteins and 20,304 human proteins. The IAV and human proteins were randomly combined, as V (1,384) × H (20,304), generating a total of 28,100,736 pairs as full-set samples.

Conjoint triad

Protein-protein interactions are mainly driven by hydrophobic effects, with hydrogen bonding and electrostatic interactions playing a crucial role (Pontremoli et al., 2018). Electrostatic and hydrophobic interactions are affected by the dipole and volume of the amino acid side chain. The conjoint triad is a method proposed by Shen et al. (2007) to classify amino acids based on dipole and side chain volume in order to extract useful information from protein sequences and convert the sequences into feature vectors. It takes into account the properties of an amino acid and its neighbors and regards any three consecutive amino acids as a single unit. Therefore, it is possible to distinguish triads based on the class of amino acids, and by calculating the frequency of each triad type, the PPI information of the protein sequence can be converted into a numerical vector. The process of generating the descriptor vector is described as follows: the 20 amino acids are categorized into seven categories depending on dipole and side chain volumes. The seven categories are {AGV}, {DE}, {FILP}, {KR}, {MSTY}, {HNQW}, and {C}. Considering three consecutive amino acids as one unit, we can get 7 * 7 * 7 = 343 triplet types, then calculate the frequency of each triplet occurrence in the given protein sequence, and finally obtain the 343-dimensional feature vector according to the following formula.

$$d_i={\displaystyle \frac{f_i-\mathrm{m}\mathrm{i}\mathrm{n}\left\{f_1,f_2,\dots ,f_{343}\right\}}{\mathrm{m}\mathrm{a}\mathrm{x}\left\{f_1,f_2,\dots ,f_{343}\right\}},i=1,2,\dots ,343}$$

Moran autocorrelation

Protein-protein interactions can be categorized into four distinct modes of interaction: electrostatic interactions, hydrophobic interactions, spatial interactions, and hydrogen bond interactions (Arenas et al., 2015). The conjoint triad approach considers the local environment of residues, but it only considers the properties of an amino acid and its two neighboring amino acids. However, information between discontinuous amino acids is equally vital, as these discrete amino acids may be spatially close for protein folding (You et al., 2014; Wang et al., 2017). By considering this information, a more comprehensive understanding of PPIs can be achieved (Guo et al., 2008), so we also use Moran autocorrelation to represent protein sequences. The autocorrelation-based descriptors capture the spatial distribution of local and global properties in proteins by calculating correlations between the physicochemical properties of amino acids at different locations in the protein (Hosseinzadeh et al., 2012). The autocorrelation descriptors also provide access to the physicochemical information contained in the protein sequence, and here eight physicochemical properties of amino acids (Xiao et al., 2015) have been chosen to reflect these interaction patterns as much as possible: normalized average hydrophobicity scales (AccNo. CIDH920105; Cid et al., 1992), average flexibility indices (AccNo. BHAR88010; Bhaskaran & Ponnuswamy, 1988), polarizability parameter (AccNo. CHAM820101; Charton & Charton, 1982), free energy of solution in water (AccNo. CHAM820102), residue accessible surface area in tripeptide (AccNo. CHOC760101; Chothia, 1976), residue volume (AccNo. BIGC670101; Bigelow, 1967), steric parameter (AccNo. CHAM810101; Charton, 1981), and relative mutability (AccNo. DAYM780201; Dayhoff, Schwartz & Orcutt, 1978). The Moran autocorrelation descriptor is defined as follows:

$$I\left(d\right)={\displaystyle \frac{{\displaystyle \frac{1}{L-d}{\sum }_{i=1}^{L-d}\left(P_i-{\overline{P}}^{\mathrm{\prime }}\right)\left(P_{i+d}-{\overline{P}}^{\mathrm{\prime }}\right)}}{{\displaystyle \frac{1}{L}{\sum }_{i=1}^L{\left(P_i-{\overline{P}}^{\mathrm{\prime }}\right)}^2}}d=1,2,\dots ,30}$$where $P_i$ and $P_{i+d}$ denote the physicochemical properties of the $i-th$ and $i+d-th$ amino acids, and d denotes the interval between two amino acids.

Feature selection

Due to the complexity of protein structures and functions, feature extraction for proteins is more intricate than for DNA and RNA sequences (Tsubaki, Tomii & Sese, 2018). Protein sequence-based features were categorized into five groups: amino acid composition descriptors, such as amino acid composition (AAC) and conjoint triad (CT), Autocorrelation descriptors, such as Geary autocorrelation (Geary) and Moran autocorrelation (Moran); pseudo amino acid composition descriptors, such as pseudo amino acid composition (PseAAC) and amphiphilic pseudo amino acid composition (APseAAC); quasi-sequence-order descriptors, such as quasi-sequence-order (QSO) and sequence-order-coupling number (SOCN); and CTD descriptors, including composition (CTDC), transition (CTDT), and distribution (CTDD) (Liu, 2017). We evaluated features from these categories on the constructed dataset using a wrapper-based feature selection method, identifying the best-performing features within each group. Subsequently, feature combination and ablation experiments were conducted to determine the optimal feature set.

Gene ontology and pathway enrichment analysis

The functional enrichment analysis of genes by gene ontology (GO) and pathway analysis were performed on target genes. The GO and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway enrichment analysis was performed via the DAVID (Database for Annotation, Visualization and Integrated Discovery) online website (Sherman et al., 2022), and the top 20 items were selected for visualization via the ggplot2 package in R (Ito & Murphy, 2013). Significant GO and KEGG pathways were selected using a statistical threshold criterion of adjusted P-value < 0.05.

Identification of hub genes and PPI network analysis

The target genes obtained through our model prediction were mapped on the STRING database (Szklarczyk et al., 2022), and the PPI was constructed by selecting a threshold of confidence >0.90. Network analysis was performed using Cytoscape (Smoot et al., 2010) to filter out the genes with degree ≥8, which were identified as potential hub genes (Tian et al., 2019). To better identify hub genes, the top 20 hub genes were screened from this network based on connectivity by CytoHubba plugin for Cytoscape. These high confidence hub genes are considered as high confidence biomarkers of the human pathway of Influenza A virus infection.

Applying our method to other datasets

To provide a more objective evaluation of the predictive performance of our constructed model, we conducted a series of comparative experiments. Firstly, we obtained the datasets used by existing PPI prediction methods proposed in the literature. Subsequently, we applied our method to these datasets, utilizing the CT+Moran features for protein sequence extraction and employing the XGBoost model for prediction. The results obtained using our method were then compared with those reported in the original article.

The two datasets published by Zhou et al. (2018) have been widely adopted as benchmarks for evaluating the performance of state-of-the-art models in viral-human PPI prediction tasks. We refer to these datasets as Zhou’s H1N1 and Zhou’s Ebola, with each dataset named after the virus represented in the respective test set.

In Zhou’s H1N1 dataset, the training set contains 10,955 true PPIs between humans and any virus other than H1N1, along with an equal number (10,955) of negative interaction samples. The test set consists of 381 real PPIs between humans and the H1N1 virus and 381 negative interactions. It can be seen in Table 3 that the sensitivity, specificity, accuracy, and MCC of our method on Zhou’s H1N1 are 91.60%, 70.07%, 80.83%, and 0.631. Their original method’s performance on sensitivity, specificity, accuracy, and MCC are 66.39%, 65.98%, 66.19%, and 0.324. The results presented in Table 3 demonstrate the superior performance of our method, particularly in terms of sensitivity and MCC, when compared to the original research approach.

In Zhou’s Ebola dataset, the training set contains 11,341 true PPIs between humans and any virus other than Ebola and an equal number (11,341) of negative interaction samples. The test set contains 150 true PPIs between humans and Ebola viruses and 150 negative interactions. Table 4 displays the performance metrics of our method on Zhou’s Ebola dataset, including sensitivity, specificity, accuracy, F1-score, and MCC. Our method achieves a sensitivity of 96.66%, specificity of 66.00%, accuracy of 81.33%, F1-score of 83.81%, and MCC of 0.658. In comparison, the original method outlined by Zhou et al. (2018) reports a sensitivity of 90.67%, specificity of 65.33%, accuracy of 78.00%, and MCC of 0.579. The results indicate that our method is better than the original method in all aspects.

We also assessed the performance of our method on the dataset developed by Prasasty et al. (2021), referred to as Prasasty’s bacteria. Prasasty’ bacteria contain three human pathogens, Bacillus anthracis, Yersinia pestis, and Francisella tularensis. They utilized two of the bacteria as the training set, while the remaining one was used as the test set for their analysis. In our evaluation, we adopted Bacillus anthracis and Yersinia pestis as training sets, Francisella tularensis as the test set. Therefore, the training set contains 6,354 positive interactions between humans and bacteria, and the corresponding negative interactions. The test set contains 1,187 positive interactions between humans and Francisella tularensis, as well as 1,187 negative interactions. The outcomes are shown in Table 5. The method utilized in the original article by Prasasty et al. (2021) achieved the following performance on this dataset: sensitivity of 74.56%, specificity of 97.83%, accuracy of 95.84%, and precision of 76.36%. However, when applied to our dataset, our method achieves a sensitivity of 93.11%, specificity of 95.28%, accuracy of 93.25%, and precision of 95.07%. Although our model exhibits slightly lower specificity and accuracy, it displays higher sensitivity and precision. As a result, our method is superior overall.

Additionally, it is worth mentioning that the three datasets we evaluated were not based on a single virus-human PPI dataset. In these datasets, the training and test sets involve different viruses, leading to low sequence similarity between their proteins. In contrast, our dataset consists entirely of influenza A virus-human PPI data, which might suggest potential dataset bias. However, our method demonstrated strong generalization ability, consistently outperforming others across diverse datasets. This robustness highlights its potential as a foundational model for PPI prediction across viruses.

Comparison with other methods on our dataset

In order to evaluate our method more comprehensively, we also applied the proposed approach by others to our dataset and conducted a comparative analysis with our own method’s results. Specifically, we employed the method presented in Denovo (Eid, ElHefnawi & Heath, 2015), which utilized CT for protein feature extraction and SVM for prediction. The comparative outcomes are summarized in Table 6. Our original results exhibit superior performance compared to Denovo. In particular, our results demonstrate an accuracy of 96.89%, precision of 98.79%, recall of 94.85%, F1-score of 96.78%, and MCC of 0.9386. These values are 6.41%, 7.29%, 6.39%, 6.83%, and 0.1291 higher than Denovo in terms of accuracy, precision, recall, F1-score, and MCC, respectively.