Multi-level machine learning prediction of protein–protein interactions in Saccharomyces cerevisiae

Residue level positives and negatives

The first step of our procedure was to build a training dataset for level-I predictor. We employed the sliding window technique to extract fragments of protein sequence. In this work we refer to it as extraction window.

Positive examples in the training set were formed from pairs of fragments in which central residues interact, and additionally certain number of other residues within a specified distance from the central one interact. We refer to the required number of interacting residues as interaction threshold and to the maximal distance from the central residue as the maximal neighbour interaction distance. By introducing this restriction we deliberately focused on strong interactions, filtering out the weaker ones, which could be just noise in the crystallisation process.

In the data preparation phase, we fixed the maximal neighbour interaction distance at 10 residues. We have chosen this value following the studies of Youn et al. (2007) and Kauffman & Karypis (2010), who predicted binding residues from protein sequences and used the window of size 21 (1 central residue, 10 residues to the right, 10 to the left). As for interaction threshold, we tested the values of 0, 5, 10, 15, and 20 interacting residues within the maximal interaction distance.

Let us observe that the values of the maximal neighbour interaction distance and extraction window size are not necessary the same. One can imagine identifying residue level positives using larger interaction distance, and then encoding features of only a few central residues using small extraction window, and vice versa. Indeed, in our work we tested different sizes of extraction window ranging from 3 to 31 independent of the fixed interaction distance.

In our data extracted from PDB we had no natural source of negative residue interactions. Pairing sequence fragments fully randomly could result in a lot of noise and false negatives. Therefore, we decided to extract non-interacting pairs of fragments only from interacting protein pairs. Fragment pairs without any interacting residues or fragment pairs in which one fragment has some interactions but the other has none were considered negatives. This way it was guaranteed that at least one of the fragments did not come from the interface region. The number of potential negatives was much larger than the number of positives, therefore we decided to keep the imbalance ratio at 3:1. The required amount of negatives was therefore sampled at random. This is a common practice in machine learning since most of the algorithms perform poorly on datasets with large class imbalance (see Chawla (2005) for review).

Protein level positives and negatives

The second step of our methodology was preparing data for training level-II predictor. We used the same dataset of interacting protein pairs from PDB database as positives, and generated negative examples. The construction of high quality negative examples is very difficult. Common methods for generating negatives include drawing random pairs of biomolecules from all known proteins found in a specific organism (Saha et al., 2014), or only from the selected subset of the whole proteome, namely from the proteins occurring in positive examples (Chang, Syu & Lin, 2010). We strongly believe that such methods have their inherent drawbacks because they ignore network properties of the underlying protein interactome. We used the following procedure instead:

• •

  Let G₁ be a graph representing positive examples. Denote V = v₁, …, v_n as the set of its vertices. Each vertex in V represents a protein and each edge v_i, v_j represents an interaction. Let [Deg(v₁), …, Deg(v_n)] be a vector containing degrees of vertices from V. Let G₂ be a graph of negative interactions. At first it has vertices identical to G₁ and no edges.

• •

  While there exist v such that Deg(v) > 0:

  1. Find vertex v with the largest Deg(v).

  2. Find vertex u if exist such that:

     1. There is no edge (v, u) in G₁.

     2. u has as large Deg(v) as possible.

     3. Distance d(u, v) in G₁ is as large as possible.

  3. If u exist:

     1. Add edge (u, v) to G₂.

     2. Deg(v)←Deg(v) − 1

     3. Deg(u)←Deg(u) − 1

  4. else: Deg(v)←0

Such schema of constructing negative protein pairs set is unbiased, i.e., the protein composition of the positives and the negatives remains identical. Every single protein has the same number of positive and negative interactions. This forces the trained classifier to predict meaningful biophysical interactions rather than predicting general reactivity (the relative number of interactions) of a single protein. Otherwise, the best results would be achieved by a predictor, which predicts that the two proteins interact if each of them has a lot of interactions in general, regardless of their actual compatibility. It’s also important that our algorithm favours protein pairs which are remote to each other in the interaction network which reduces the risk of introducing false negatives.

Train-test split

The last step of data generation procedure is splitting samples between training and testing sets. In order to truly evaluate our method in a realistic setup, we split the benchmarking dataset at the protein level, not at the residue level. This made our goal more difficult as compared to previous works that often used residue-level splitting of benchmarking dataset. The schematic depiction of the train-test split is given by Fig. 2. The Level I classifier was evaluated through train-test experiment with a relatively large datasets. The Level-II classifier was evaluated through cross-validation experiment which used the testing set only—information from the training set came only in the form of the trained level-I classifier. Such schema eliminated the risk of overoptimistic performance estimates caused by the same data appearing during training and testing phases.

Level-I features

We trained the level-I predictor on interacting sequence fragments of proteins from the training set and tested it on the testing set. Input for level-I predictor consisted of pairs of sequence fragments of the length of extraction window.

The following sets of features were considered:

• Raw sequence—raw sequence of amino acids encoded numerically.

• HQI8—sequence of amino acids encoded with High Quality Indices (Saha, Maulik & Plewczynski, 2012). These features were constructed by clustering AAindex database (Kawashima et al., 2007). Each amino acid is described by 8 values representing its physicochemical and biochemical properties.

• DSSP structure—secondary structure of the protein extracted from PDB complex with DSSP software. It was limited to the three basic symbols: E—β-sheet, H—α-helix, C—coil.

• PSIPRED structure—secondary structure of the protein predicted from sequence with PSIPRED software. It was limited to the three basic symbols: E—β-sheet, H—α-helix, C—coil.

As the core classifier, we evaluated two popular machine learning methods: Random Forest and Support Vector Machine. Both algorithms are commonly used in bioinformatics and are considered the best off-the-shelf classifiers (Yang et al., 2010). Their parameters values were chosen through a grid search. As the performance measure we have chosen ROC AUC (area under Receiver Operating Characteristic curve).

Level-II features

To infer a binary interaction between two proteins, we consider all possible interactions between their sequence segments as predicted by level-I predictor. An output of this phase is a matrix of likelihoods with the dimension equal to the multiplied proteins’ lengths. Each prediction score is a real number between 0 and 1. Sample matrices of scores for a positive and negative case are presented in Fig. 3.

To transform the 2D matrix into an input vector suitable for level-II predictor, we extracted the following features (numbers in parentheses denote the number of values in the final feature vector):

• the mean and variance of values over the matrix (2),

• the sums of values in 10 best rows and 10 best columns (20),

• the sums of values in 5 best diagonals of the original and the transposed matrix (10),

• the sum of values on intersections of 10 best rows and 10 best columns (1),

• the histogram of scores distributed over 10 bins (10),

• graph features: fraction of nodes in the 3 largest connected components (3).

The graph features require further explanation. Predicted contacts between residues were represented as a bipartite graph. Nodes in the graph represented residues and edges represented predicted contact. To make the graph more realistic biologically, for each node we left only 3 strongest outgoing edges. We set the value of this threshold (3) following the observation that in our PDB structures the mean number of interactions of a single interacting residue is between 2 and 3. In the trimmed graph we calculated fractions of nodes contained in 3 largest connected components. Those values were also appended to the feature vector.

Evaluation of the level-II predictor

The performance of the level-II predictor was evaluated through a variant of stratified 30-fold cross-validation performed on the protein level. Each fold contained $\frac{1}{30}$ of positive protein pairs and $\frac{1}{30}$ of negative protein pairs from the testing set. There was no overlap between splits on the pair level, but there was still an overlap on the level of single proteins which constitute pairs. We observed that this introduced a huge bias into evaluation results (similar observation was previously made by Park & Marcotte (2012)). To fix our cross-validation scheme we applied the following procedure:

1. Let $O=\left\{\left(p_1^1,p_1^2\right),\left(p_2^1,p_2^2\right),\dots ,\left(p_n^1,p_n^2\right)\right\}$ be a set of all protein pairs.

2. For each fold F ⊂ O:

   1. Build a set P composed of all proteins occurring in F: ${\displaystyle P=\left\{x:{\exists }_y\left(x,y\right)\in F\vee \left(y,x\right)\in F\right\}.}$

   2. Build a set A ⊂ O composed only of pairs consisted of proteins occurring in P: ${\displaystyle A=\left\{\left(x,y\right):x\in P\wedge y\in P\right\}.}$

   3. Build a set B ⊂ O composed only of pairs consisted of protein not occurring in F but occurring in the testing set: ${\displaystyle B=\left\{\left(x,y\right):x\notin P\wedge y\notin P\wedge \left(x,y\right)\in O\right\}.}$

   4. Train the classifier on B set, and test it on A.

3. Collect all the predictions for A-sets, and calculate performance metrics.

The above described procedure differs from the standard cross-validation, since the number of observations in constructed test sets vary slightly, but this variance is small, and does not influence the estimated performance. Such evaluation schema does not allow for any information leak: the datasets are always balanced, and the classifier is tested on previously unseen proteins.

As classification methods for level-II predictor we used Random Forest and Support Vector Machine with parameters tuned through a grid search.

Protein sequence feature aggregation

We compared our ensemble method with various sequence feature aggregation schemas that are commonly applied in machine learning predictors of proteins interactions. To make the benchmarking results comparable between different algorithms, we used the same classification method (Random Forest) and evaluation procedure (modified 30-fold cross-validation on the testing set) as for level-II predictor. We benchmarked the following feature aggregation schemas:

1. AAC—Amino Acid Composition (Nanni, Lumini & Brahnam, 2014). Feature set is the set of frequencies of all amino acids in the sequence.

2. PseAAC—Pseudo Amino Acid Composition (Chou, 2001). Feature set consists of the standard AAC features with k-th tier correlation factors added. We calculate those correlations on HQI8 indices.

3. 2-grams (Nanni, Lumini & Brahnam, 2014). Feature set comprises of frequencies of all 400 ordered pairs of amino acids in the sequence.

4. QRC—Quasiresidue Couples (Guo & Lin, 2005). A set of AAIndices is chosen. For each index d combined values of this property d for a given amino acid pair are summed up for all the pair’s occurrences over the full protein sequence. Occurrences for pairs of residues separated from each other by 0, 1, 2…m residues. In effect, one obtains QRC^d vectors of length 400 × m. In this model we also use HQI8 indices.

5. Variation of Liu’s protein pair features (Liu, 2009). The method starts from encoding each amino acid in a protein sequence with 7 chosen physicochemical properties, thus obtaining 7 feature vectors for each sequence. For each feature vector its “deviation” is calculated: ${\displaystyle {\gamma }_{dj}=\frac{1}{n-d}\sum _{i=1}^{n-d}{x}_{ij}\times x_{\left(d+i\right),j}j=1,\dots ,7d=1,\dots ,L}$ where x_ij is the value of descriptor j for amino acid at position i in sequence P, n is the length of protein sequence P, and d is the distance between residues in the sequence. For the purpose of the comparison, we tested this method with the original 7 amino acid indices used by Liu, as well as with HQI8 features. We tested different values of L from 5 to 30 in a quick cross-validation experiment on our data and chose L = 9 as yielding the best results.

Experimental results

We evaluated carefully all subsequent steps of our method to choose optimal features and parameter values. Then we compared performance of level-II predictor with popular sequence encoding schemas. In our experiments we stuck to the rule that during classifier training we can use all the information available in PDB complexes, but in the evaluation phase only information derived from the sequence is allowed. This was to demonstrate that our method can be employed successfully in a situation when only protein sequences are known.

The first task was to decide on the set of optimal features for the level-I predictor. We have chosen ROC AUC as the performance metric. To have a complete picture, we measured both the performance of level-I predictor and the performance of level-II based on the trained level-I predictor. We tested two sources of secondary structure information—DSSP and PSIPRED—separately, but the evaluation was done using PSIPRED secondary structure (results on DSSP are given in parentheses for completeness). Initial experimentation was done with Random Forest, then we tested whether it is possible to gain anything by replacing it with a carefully tuned SVM.

Evaluation results of the predictor on different features for 5 interaction thresholds are presented in Table 1. The DSSP secondary structure performed better than sequence features. Including both secondary structure and HQI8 in the feature vector provided minor improvement on level-I, but it resulted in performance degradation on level-II. While the quality of level-II prediction is based on the quality of level-I prediction, our results show that overfocus on optimizing level-I prediction may result in overfitting visible on level-II. SVM did not perform better than Random Forest. In all further experiments we used secondary structure as the only source of features for level-I predictor, together with Random Forest as the learning algorithm.

Table 2 reports ROC AUC values for different interaction thresholds. Each value changes the definition of positive and negative example for level-I predictor. The larger the threshold, the more positive and negative examples differ. We can see that increasing interaction threshold indeed made level-I prediction task easier, but the impact on level-II prediction was non-linear. We chose 15 as the optimal threshold value.

The final decision was selecting the optimal extraction window size. Figure 4 presents level-I ROC AUC score against window size. After analysing the plot we decided to keep size 21 as it provided good performance and it was previously used in other publications. In that way we fixed extraction window size to the same value as the maximal interaction distance, used to define positive residue interactions.

After fixing the parameters, we wanted to choose classification algorithm for level-II predictor and compare performance of our multi-level representation with representations based on the aggregated protein sequence. Results are presented in Table 3. Once again Random Forest proved to be more suitable for the task than Support Vector Machine. On this kind of data level-II predictor outperformed other methods significantly.

Role of secondary structure

Our results stress the link between secondary structure and residue contacts. Compatibility of structural motives is important for contact forming. To get more insights from our results, we extracted feature importances from the trained Random Forest model of level-I predictor. Relative importance of an attribute in a decision tree is expected fraction of examples split by nodes based on this attribute. In Random Forest this value is averaged over all trees. Relative feature importances for the input vector constituted by secondary structure annotation are given by Fig. 5. It is visible that the most important are central and border residues of the fragments.

To understand what is happening on these positions, we calculated frequency of particular secondary structure patterns occurring between the residues of interacting and non-interacting fragments. Calculated frequencies are presented in Table 4. It is clear that pairs H–H occur more frequently among interacting fragments than non-interacting for all three positions. E–E pair is more likely to occur in interacting fragments on the central position, but not so much on the border. C–C pair is generally more prevalent among non-interacting fragments, but this effect is also stronger for central residues. Those results are not surprising: since secondary structure motives are stabilised by hydrogen bonds the relation to both intra-molecular and inter-molecular contacts is expected (Hubbard & Kamran Haider, 2010).

From our experiments it is clear that using the real secondary structure for classifier training is preferable over the predicted secondary structure, even if later only predicted secondary structure will be available. This limits the noise introduced during the training phase allowing the predictor to focus on the important patterns.

What is surprising at first glance is that encoding both secondary structure and protein sequence in a single feature vector improves the prediction on level-I but leads to a worse prediction quality on level-II. It is true that protein sequence contains more information than just secondary structure, but, on the other hand, including more attributes increases the risk of overfitting. This is the case of so-called bias–variance trade-off, where increasing the model complexity potentially decreases classifier bias but at the same time increases classifier variance (Geman, Bienenstock & Doursat, 1992). Here, the situation is even more complex because overfitting may occur at different levels—improving the prediction quality at level-I may in some cases make level-II task more difficult. From various feature combinations considered in our study secondary structure alone proved to be optimal, conveying important information while keeping the model complexity low.

We may speculate that this phenomena comes directly from the hierarchical nature of biological systems. Biological processes have to be robust and predictable, harnessing the dynamic of physical molecules in a constructive way. It is easier to build the upper layers of a complex system using constrained and standardized building blocks. Secondary structure motives are such stable building blocks, more constrained than the protein sequence. In that sense secondary structure may be seen as an useful data compression optimized evolutionary to perform certain functions. Our classifier performs better with a more compact, specialized representation than with the full sequence information, which requires too many free parameters to estimate.

Quality of residue contacts prediction

In our work we focused on predicting binary protein interaction and the constructed residue-level predictor served only to generate data for level-II predictor. We optimised level-II prediction quality and demonstrated that optimising level-I only could lead to overfitting on level-II. In fact, for the threshold values greater than 0, our level-I predictor was not trained to identify single residue contacts anymore; it focused on discovering only strongly interacting fragments important for inter-protein interactions. To assess the importance of this difference in task formulation, we calculated ROC AUC score for level-I predictor not on our filtered segment pairs data but on real contact matrices of the test proteins. The predictor obtained 0.65 AUC—a score lower than on the filtered data and lower than level-II predictor score. This means that level-I predictor produces useful information for predicting protein–protein interactions but is not necessarily a good predictor of specific residue contacts.

After examining predicted contact maps visually, we realised that they consist of characteristic horizontal and vertical lines. This means that certain sequence fragments are predicted to be very active and likely to interact while the others are classified as inactive. We considered the possibility that our predictor focused on predicting solvent exposure of the fragments. To verify this hypothesis we calculated relative solvent accessibility for each residue (Shrake & Rupley, 1973). We chose threshold value 0.2 to distinguish between buried and accessible residues (Chen & Zhou, 2005). All residues of proteins from our test set were divided into these two classes. Then we tested whether mean values of rows and columns of the predicted contact matrices are useful for predicting residue state (buried or accessible). Surprisingly, we obtained AUC ROC score 0.46, which means that characteristic rows and columns in our matrices are not directly correlated with just solvent accessibility.

It is possible that our level-I predictor reflects more specific protein properties such as the location of active sites. Improving our method would require an in-depth analysis of these matters. Since solvent accessibility is not correlated with the current predictions perhaps it is possible to include this information to further improve accuracy.

Quality of protein interactions prediction

We draw the reader’s attention to the fact that the performance of popular protein representation strategies evaluated on our data was generally much lower than results reported in the literature. One of the reasons may be a relatively small size of our dataset—it might not contain enough examples for a classifier to learn complex patterns. The other explanation is the way we constructed positives and negatives. In our case every protein occurred in the same number of positive and negative pairs. Moreover, we performed cross-validation with splits on the protein level, which means that no single protein occurred simultaneously in training and testing set. In such conditions any method which makes a good prediction of general proteins’ reactivity but does not consider their actual compatibility performs poorly. This observation is consistent with results obtained by Park & Marcotte (2012), who analysed the impact of performing splits on the component level instead of on the pair level on cross-validation results. They reported AUC scores of popular protein interaction prediction methods dropping from 0.7–0.8 to 0.5–0.6. We strongly believe that performing splits on the component level, i.e., keeping training and testing protein sets disjoint, results in a better estimation of predictor performance in a realistic scenario of predicting interaction between two previously unseen proteins.

Even though our evaluation procedure was carefully designed to reduce certain biases, we still do not expect it to reflect the reality of protein–protein interaction prediction perfectly. The most important issue is the proportion of positives to negatives in our datasets. On the protein level, we used balanced sets with the same number of positives and negatives. In reality the number of negatives, i.e., pairs of non-interacting proteins is much larger than the number of positives. Grigoriev (2003) estimated that there are around 26,000 hetero interactions between 6,300 proteins in yeast proteome. This would make the ratio of negatives to positives equal to: ${\displaystyle \frac{\left(\genfrac{}{}{0ex}{}{\text{6,300}}{2}\right)-\text{26,000}}{\text{26,000}}=762.1481.}$ It is not feasible to construct a dataset with such a ratio in practice. Since the negatives in our experiment are constructed through pairing proteins randomly, there might be false negatives present. The more negatives we draw, the more false negatives we will introduce. In the most pathological case all the unknown interactions we would like to predict will be labelled as negatives. One of the possible solutions to this problem is to calculate performance metrics using balanced data and then introduce corrections based on the expected negative to positive ratio.

Classifier performance is calculated from a confusion matrix consisting of the following terms: TP—true positives (positives classified as positives), TN—true negatives (negatives classified as negatives), FP—false positives (negatives classified as positives), FN—false negatives (positives classified as negatives). The correction involves multiplying terms based on negative examples in the data set—TN and FP—by constant d = 762. Let us examine the definitions of the specific performance metrics: ${\displaystyle \mathrm{Accuracy}=\frac{\mathrm{TP}+d\mathrm{TN}}{\mathrm{TP}+d\mathrm{TN}+d\mathrm{FP}+\mathrm{FN}}}$ ${\displaystyle \mathrm{Recall}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}}$ ${\displaystyle \mathrm{Precision}=\frac{\mathrm{TP}}{\mathrm{TP}+d\mathrm{FP}}.}$ AUC is the area under receiver operating characteristic, which is a plot of recall (sensitivity) vs 1-specificity, where specificity is defined as: ${\displaystyle \mathrm{Specificity}=\frac{d\mathrm{TN}}{d\mathrm{TN}+d\mathrm{FP}}=\frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}}.}$ From this equation it is clear that class imbalance does not affect recall and AUC. In the case of the final level-II predictor TP = 179, FP = 78, TN = 160, FN = 83, the corrected precision is 0.003, the corrected accuracy is 0.67. Precision score 0.003 is hardly impressive but it is still more than 2 times bigger than the score 0.0013 expected by chance when labelling protein pairs blindly. This means that the use of our predictor introduces valuable information.

In the light of these results, we believe that current methods for predicting protein–protein interactions from sequences, including ours, are not mature enough for large-scale application. The goal of reconstructing full protein interaction network is right now beyond reach. However, the developed predictors could still be applied in a variety of contexts. One of the possible applications is the initial screening of candidates for more costly in vitro experiments—choosing protein pairs for experimental verification in a more informed way would make the procedure more effective. Another possibility is using the predictor as a filter to improve quality of the interaction data collected from the existing databases. This is also an important issue because the quality of results from high-throughput experiments is still limited.