As shown in Figure 1B, the DeepPPISP proposed by Zeng et al. (Zeng et al., 2020) for PPIS prediction had three types of input: position-specific scoring matrix (PSSM), secondary structure, and raw protein sequences. The PSSM is an excellent feature extractor for protein sequences and thus have widely been applied to problems in the field of computational biology, such as predicting protein post-translational modification (Huang et al., 2013; Huang et al., 2014; Dehzangi et al., 2017), membrane type (Wang et al., 2019), protein-RNA binding site (Liu et al., 2021), and structure (Guo et al., 2021). The quality of PSSM features is closely associated with the underlying multiple sequence alignments. Although there are many multiple sequence alignment algorithms including HIMMER (Eddy, 2011; Wheeler and Eddy, 2013) (Johnson et al., 2010) and Hhbilits (Remmert et al., 2012), PSI-BLAST (Altschul et al., 1997) is still a popular multiple sequence alignment and homology search algorithm. Here, PSI-BLAST was used to search NCBI’s non-redundant (NR) sequence database with three iterations and an E-value threshold of 0.001.

Many protein-protein interfaces are related to secondary structures (Taechalertpaisarn et al., 2019). Information about protein secondary structure is helpful to predict PPIS. The DSSP program (Touw et al., 2015) was used to generate nine state secondary structures: α-helix, 3₁₀- helix, π-helix, β-bridge, β-strand, β-turn, bend, loop or irregular, and no secondary structure. Therefore, each amino acid residue corresponded to a 9-dimensional vector. The primary protein sequence is valuable information and thus is essential to predict protein properties. One-hot encoding was used to encode the protein sequences. There are 20 kinds of common amino acids in the protein sequences, so each amino acid residue corresponds to a 20-dimensional 0/1 vector. The protein-protein interaction is closely associated with neighboring residues of interaction sites. The local feature of interaction sites contributes to the identification of PPIS. The sliding window method was used to collect the neighboring residues of the interaction sites. The size of the sliding window was seven. For example, if the interaction site was at position i, residues at position i-3, i-2, i-1, i, i+1, i+2, and i+3 were separated. Because each residue corresponds to a 20-dimensional PSSM feature, a 9-dimensional secondary structure feature, and a 20-dimensional one-hot feature vector, a window of seven amino acid residues was encoded into a 343-dimensional vector which was called the local feature.

Protein-protein interaction is not only linked to the local information of interacting sites, but also to global information. Zeng et al. (2020) demonstrated that the inclusion of global information improved the performance of predicting PPIS. A 500-residue peptide was used to represent the global feature of PPIS. If the number of amino acid residues in the protein sequence was less than 500, it was padded with a 0. Each peptide corresponds to a 500*49-dimensional vector called a global feature.

The local and the global features were fed into the DeepPPISP (Zeng et al., 2020). The DeepPPISP was made up of one embedding layer, three different scale convolutions, two fully connected layers, and an output layer (Figure 1B). For more detail, readers can refer to the reference (Zeng et al., 2020).

Both the local features or global features would contain a certain degree of noise. The dimension is large, especially for global features. The DeepPPISP was used to extract a more informative representation. The DeepPPISP was trained on the training data in a supervised manner. The local and global features were fed into the trained DeepPPISP, and the input to the first fully connected layer was the abstract representation of the raw features. Compared with the raw features, the abstract representation was of low dimension and had low noise.

The XGBoost proposed by Chen and Guestrin, 2016 belongs to Gradient Boosting Decision Tree (GBDT) (Ke et al., 2017), and both are tree boosting algorithms. Compared with traditional tree boosting, the XGBoost used a theoretically justified weighted quantile sketch for approximate learning, a novel sparsity aware algorithm for handling sparse data, and an effective cache-aware block structure for out-of-core tree learning (Chen and Guestrin, 2016). In addition, the XGBoost performed faster as it exploited parallel and distributed computing. The XGBoost has such a significant superiority that it has widely been used in many areas including machine learning and data mining challenges.

The XGBoost is an addition model. At each iteration, the XGBoost learns a new tree that fits the residual between the predicted result of the previous trees and the true values of the training samples.

Assume that D = { ( x i , y i ) | | D | = n , x i ∈ R m , y i ∈ R } denotes a training set, where m and n represented the numbers of features and samples, respectively. At the t-th iteration, the aim of the XGBoost is to learn a function f t so that y ^ i t = y ^ i t − 1 + f t ( x i ) (1) where y ^ i t − 1 is the fitting value of the previous t−1 trees for the i-th sample. To search for f t , the loss function with the regularization was used as the objective function: obj =   ∑ i = 1 n l ( y i , y ^ i t ) + ∑ i = 1 n Ω ( f t ( x i ) ) = ∑ i = 1 n l ( y i , y ^ i t − 1 + f t ( x i ) ) + Ω ( f t ) + constant , (2) where l was the loss function which was generally defined as l ( y i , y ^ i t ) = ( y i − y ^ i t ) 2 . (3) ∑ i = 1 t Ω ( f i ) denotes the regularization. The loss function l was approximated by the second-order Taylor series, namely l ( y i , y ^ i t − 1 + f t ( x i ) ) ≈ l ( y i , y ^ i t − 1 ) + g i f t ( x i ) + 1 2 h i f t 2 ( x i ) , (4) where g i = ∂ l ( ( y i , y ^ i t − 1 ) ∂ y ^ i t − 1 and h i = ∂ 2 l ( ( y i , y ^ i t − 1 ) ∂ y ^ i t − 1 ∂ y ^ i t − 1 were the first- and the second-order gradients of the loss function with respect to y ^ i t − 1 respectively. Ω ( f t ) was defined by Ω ( f t ) = γT + 1 2 ∑ j = 1 T ω j 2 , (5) where T was the number of leaf nodes and ω j was the weight of the j-th leaf node. The objective function was equivalently rewritten as  obj = ∑ i = 1 n [ g i f t ( x i ) + 1 2 h i f t 2 ( x i ) ] + γT + 1 2 λ ∑ j = 1 T ω j 2 . (6)

The set of instances of the leaf node j was defined by I j = { x i | q ( x i ) = j } . (7)

The objective function was further represented as obj = ∑ j = 1 T ( ∑ i ∈ I j g i ) ω j + 1 2 ∑ j = 1 T ( ∑ i ∈ I j h i +  λ ) ω j 2 + γT (8)

Given a fixed tree q ( x ) , the optimal value of each leaf node was calculated by ω j * = ∑ i ∈ I j g i ∑ i ∈ I j h i +  λ , (9) and the optimal value of the whole tree was calculated by  ob j * = − 1 2 ∑ j = 1 T ( ∑ i ∈ I j g i ) 2 ∑ i ∈ I j h i +  λ + γT . (10)

It was expensive and impossible to exhaust all the possible trees for the training data. In practice, the greedy algorithm was used, which started from one node and iteratively split the node. Assume that before the node was split, the objective function of the tree was ob j 1 = − 1 2 ∑ j = 1 T − 1 ( ∑ i ∈ I j g i ) 2 ∑ i ∈ I j h i +  λ + γT − 1 2 ( ∑ i ∈ I k g i ) 2 ∑ i ∈ I k h i +  λ . (11)

After the node k was split into the left tree   I L and the right tree I R , the objective function was  ob j 2 = − 1 2 ∑ j = 1 T − 1 ( ∑ i ∈ I j g i ) 2 ∑ i ∈ I j h i +  λ + γ ( T + 1 )   − 1 2 ( ∑ i ∈ I L g i ) 2 ∑ i ∈ I L h i +  λ — 1 2 ( ∑ i ∈ I R g i ) 2 ∑ i ∈ I R h i +  λ . (12)

The gain of node splitting was calculated by  gain = ob j 1 − ob j 2 = 1 2 ( ∑ i ∈ I L g i ) 2 ∑ i ∈ I L h i +  λ + 1 2 ( ∑ i ∈ I R g i ) 2 ∑ i ∈ I R h i +  λ − 1 2 ( ∑ i ∈ I k g i ) 2 ∑ i ∈ I k h i +  λ − γ . (13)

The gain was used to assess the split candidates.

To investigate the ability of the features to discriminate protein-protein interaction sites from non-interaction sites, we used the Uniform Manifold Approximation and Projection (UMAP) (McInnes et al., 2020) to depict the first two principal components. The UMAP is a powerful tool for dimension reduction and visualization. As shown in Figure 2, the features processed by the DeepPPISP demonstrated a tighter cluster than the raw features, indicating that features generated by the DeepPPISP were more discriminative. To further evaluate the performance of the preprocessed features, we performed 5-fold cross-validation and independent tests. Figure 3A showed the ROC curves of the 5-fold cross-validation over both the preprocessing features and raw features, while Figure 3B depicted the ROC curves of the independent tests. The performance of preprocessed features is equivalent to or better than those of raw features. It must be pointed out that the user-defined parameters were identical in the XGBoost classifiers. Comparison with other methods

Due to its versatile roles in the cellular process, the identification of protein-protein interaction sites is increasingly becoming a hot topic and is also a challenging task. Over the past decades, more than 10 methods have been proposed to predict protein-protein interaction sites (Patel et al., 2006; Du et al., 2009; Murakami and Mizuguchi, 2010; Wang et al., 2014; Zhang and Kurgan, 2019; Northey et al., 2018; Zeng et al., 2020; Chen et al., 2012; Šikić et al., 2009; Fiorucci and Zacharias, 2010; Dosztányi et al., 2009; La and Kihara, 2012; Bradford and Westhead, 2005; Chen and Jeong, 2009; Chung et al., 2006; Fernandez-Recio et al., 2004; Shoemaker et al., 2010; Ofran and Rost, 2007; Qin and Zhou, 2007; Liang et al., 2006; Li et al., 2007; Zhou and Shan, 2001; Neuvirth et al., 2004; Porollo and Meller, 2007; Segura et al., 2011; Qiu and Wang, 2012; Wei et al., 2016; Zhu et al., 2020; Guo et al., 2018; Kuo and Li, 2016; Wang et al., 2021; Maheshwari and Brylinski, 2015; Li, 2020; Dick and Green, 2016; Wang et al., 2019; Zhao et al., 2017; Jia et al., 2016; Deng et al., 2020; Singh et al., 2014; Hou et al., 2017; Li et al., 2012; Wang et al., 2019; Bagchi, 2015 #412; Zhang and Kurgan, 2019). We compared the proposed method with six other state-of-the-art methods. These six competing methods were DeepPPISP (Zeng et al., 2020), SCRIBER (Zhang et al., 2019), IntPred (Northey et al., 2018), RF_PPI (Hou et al., 2017), SPRINGS (Singh et al., 2014), PSIVER (Murakami and Mizuguchi, 2010), ISIS (Ofran and Rost, 2007), and SPPIDER (Porollo and Meller, 2007). PSIVER was a Naïve Bayes-based classifier that used features from PSSM and accessibility, while SPPIDER combined fingerprints with information from the sequences and structures for PPIS predictio. Both SPRINGS and ISIS were neural network-based methods. The former used evolutionary information, averaged cumulative hydropathy, and predicted relative solvent accessibility, while the latter used structural features and evolutionary information. RF_PPI was a random forest-based classifier for PPIS prediction, while the DeepPPISP was a deep learning-based classifier. The performances of these seven methods over the independent test were listed in Table 1.

The DeepPPISP-XGB method achieved the highest value in terms of Recall, F1-score, AUROC, AUPRC, and MCC, and it reached the second-highest performance in terms of Precision. Although ISIS got the best ACC, its performance in other respects was lower than those of DeepPPISP-XGB. The DeepPPISP-XGB method improved the Recall by 4.7%, 5.5%, 11.6%, 11.2%, 2.6%, 15.6%, 26.2%, and 16.5%, in comparison with DeepPPISP, SCRIBER, IntPred, RF.PPI, SPRINGS, PSIVER, ISIS, and SPPIDER, respectively. The DeepPPISP-XGB method increased F1-score and MCC by 0.5% and 0.3%, and the AUROC by 1%, in comparison with DeepPPISP.

K-fold cross-validation is a common method in regression or classification questions. In the k-fold cross-validation, the training set was split into k parts. One part was tested and other k−1 parts were trained. The procedure was performed k times. We carried out 10-fold cross-validations, and the principle was shown (Supplementary Figure S1). Figure 4 showed ROC curves for the 10-fold cross-validations. The mean and the standard deviation of the AUROCs were 0.741 and 0.006, respectively. Supplementary Table S1 lists the ACC, Precision, Recall, F1-score, AUROC, AUPRC, and MCC for each cross-validation.

To further evaluate the predictive performance of the DeepPPISP-XGB method, four machine learning algorithms were used for PPIS prediction. Decision tree (Safavian and Landgrebe, 1991) is a widely utilized classification algorithm, which is made up of the root node, internal nodes, and leaf node. Random forest (RF) (Breiman, 2001) is an ensemble learning algorithm. It consists of many weak classifiers which determine the sample category. Extremely randomized tree (ERT) (Geurts et al., 2006) is similar to RF but the decision tree of ERT is randomly divided. Support vector machine (SVM) is a statistical algorithm proposed by Boser et al. (Boser et al., 1992). These classifiers were implemented in the Scikit-Learn package (v0.24.2) which has been widely utilized in computational biology. The ROC curves and the precision-recall curves are shown in Figure 5. The XGBoost classifier obtained an AUROC value of 0.681 and an AUPRC value of 0.339 on the independent test, significantly better than four classifiers.

After removing global features, we trained DeepPPISP-XGB. The user-defined parameters of the DeepPPISP-XGB were the same as the previous. Table 2 shows the performance of predicting PPIS by using local features alone. The ROC and the precision-recall curves were displayed in Figure 6. The experimental results showed that the inclusion of the global features was beneficial to improve PPIS prediction, which was in agreement with the findings of Zeng et al. (2020).