Protein-protein interactions (PPIs) play a central role in virtually all cellular processes. Proteins interact with other proteins to accomplish specific biological functions, such as DNA replication or RNA transcription, gene translation, gene regulation and protein transport, as well as signal transduction. Identification of interaction sites between two binding proteins is essential to understand complex formation and investigate their function (e.g., [13, 20, 53]). In particular, information about specific amino acid residues that play essential roles in protein interactions usually has a wide range of applications such as design of the targets of drugs and antimicrobials (e.g., [40]).

Despite continual improvement, certainly over the last decade, experimental techniques for large scale determination of PPIs are not yet able to provide comprehensive coverage over all PPIs in the detail needed to allow better understanding of the evolutionary and physical forces that govern them (e.g., [13, 23, 24]).

During the past decades, several types of computational methods have been developed for protein interaction prediction. Docking and modeling approaches that rely mainly on surface complementarity and electrostatics to predict structural complexes. These approaches fit together two known structures through their interacting surfaces, or predict protein–protein interaction sites from known monomer structures [13, 23, 41, 47]. However, these methods require structure information of proteins, which remains relatively scarce and expensive. Therefore with the increasing amount of sequence data from sequencing initiatives, a method that only uses sequence information without known structures to predict protein-protein interaction sites is becoming increasingly attractive.

Several such computational methods aim to predict the possibility of interaction between two proteins [15, 35, 52]. Perhaps the most well-known technique for predicting PPIs from sequence data is the ‘mirror tree’ method (e.g., [38]). This method infers interactions from the correlation of evolutionary patterns, as seen in phylogenetic trees representing each of the interaction partners. However, this correlation may instead arise from functional relatedness as well as a number of other general evolutionary mechanisms [8, 29].

Predicting intra-protein and inter-protein residue-residue contacts from sequence covariation has recently revived [19, 22, 32–34, 51, 55]. This is directly due to the availability of large amounts of sequence data and the recent development of so-called direct-coupling methods (e.g., [55]. The idea has been studied in the eighties (e.g., [1]) and nineties (e.g., [16, 27]). The main limitation of these methods is that, typically, five times more sequences than the alignment length (the ‘5 L’ rule) are required [21, 22, 32, 33, 51]. For most proteins, this is not available. In addition, the application to inter-protein residue contacts is hampered by the need to construct large correlated alignments. Here, for each sequence an ortholog in the other alignment must be included, so that positional variations of the alignment of one interaction partner may be correlated with those of the other protein.

For identifying protein-protein interaction (PPI) sites, often conservation measures on sequence features are used [45]. For example, ISIS by Ofran and Rost combine PSI-blast profiles and predicted solvent accessibility and secondary structure to predict interface sites [35, 36]. SPPIDER [43] uses in addition several structure-derived features in an elaborate Machine learning approach. In addition, sequence and network features [12, 48], as well as conservation in combination with specificity [31] are also used to predict interaction sites. Several findings indicate that the interface rim tends to be more conserved than the interface core (e.g., [5, 18, 44]), while localized conservation of single residues can indicate interaction hot spots [9, 35, 50]. At the level of PPI networks, mixed results are being reported. Some conserved PPI network motifs appear related to conserved sequence motifs [12, 48]. Overall conservation patterns, however, are found to be weak and mostly not significant (e.g., [28, 42]).

Although progress has been made in predicting binding sites from sequence information, the problem remains far from solved and several limitations persist. First, extracting evolutionary information from sequence data critically depends on sequence alignments containing large numbers of sequences. Second, most methods rely on a combination of structural and sequence features (e.g., [52, 54]). While combined methods can achieve high prediction performance, the performance of sequence-only methods remains modest [35, 37, 42].

Specificity of interaction, i.e. differences between groups of homologs that display different interactions has previously been reported. Pirovano et al. [39] identified interface residues by comparing homologs with different binding partners. Manning et al. [31] predicted positions which define sequence subfamily specificity, where some of these positions were binding sites. Based on a dataset of yeast interaction data and fungal ortholog groups, it has been suggested that, in addition, specificity between non-interacting and heteromeric interacting protein pairs might be used to detect the interaction sites [14]. Interestingly, here only up to one hundred sequences were needed to detect the specificity signal between binding and non-binding groups, far fewer than the ‘5 L’ needed for covariation-based methods. However, the performance of their predictions is only just above random, indicating a need for a cleaner dataset for obtaining proof of principle.

In this paper, we investigate whether specificity between interacting and non-interacting subgroups can be used to predict interaction sites. To address this question, we chose homodimers as a use case to construct interacting subgroups and monomers to constitute non-interacting subgroups. In this way, we can confirm that all sequences in the interacting subgroup physically interact, and that we have a sub-group of monomers known not to (self) interact. Furthermore, the specificity signal is from compositional differences of one chain rather than multiple chains as would be the case when comparing heteromeric interacting groups with non-interacting groups. All homodimers and monomers were obtained from PISA which is a resource for exploring marcromolecular interfaces [26]. We compiled a new database derived entirely from crystallized proteins in the PDB [4], and compared homodimers with homologous non-interacting monomers in a multiple sequence alignment. Starting with 9152 homodimers and 13,355 monomers, we constructed 1,592 pair groups for which we predicted putative interface residues. We found that the compositional differences between interacting and non-interacting subgroups pinpoint interface positions. We also found that various filters on the input sequences yielded a stronger specificity signal and a better prediction performance. Finally, we relate our method with a sequence-only method, SPPIDER [43].

We applied our analysis on the 10 pair-group sets created based on sequence identity (i.e., the complete set, <99 %, <98 %, <95 %, <90 %, <80 %, <70 %, <60 %, <50 % and <40 %) for Test-set 1. The number of pair-group alignments in each of the sets can be found in Additional file 1: Table S1. With lower sequence identity cut-off values, the numbers of homodimer and monomer sequences within each of the groups and the total number of groups are reduced. The complete set (Test-set 1) consists of 1593 pair-group alignments, each containing on average 25 homodimer and 14 monomer sequences per group.

Conservation differs for buried sites

From the alignments created, for each position the entropy was calculated to quantify the evolutionary conservation pattern. Figure 2 shows the conservation difference between three types of positions: interface, other surface and buried residues in the homodimer subgroup. For each alignment, the overall average entropies for each of the three types of positions were calculated for the homodimer subgroup.

As described in methods, interface sites stand for the residues for which accessible surface area becomes buried during association; surface residues are solvent-accessible, but not interface; and buried residues are inaccessible. When filtering on maximum sequence identity from <40 % to <100 % (<100 % yielding the complete dataset), average entropies went down because more similar sequences were included in the datasets. It was not surprising that the buried (inaccessible) sites were more conserved than the other two (interface and other surface sites). During evolution, proteins usually conserve their hydrophobic core to keep structural stability. However, the conservation differences between interface and the rest of the protein surface were small over the whole range of % ID cut-offs.

In summary, it is easy to separate the buried sites from surface residues because the conservation pattern differs. However, it is virtually impossible to distinguish the interaction sites from other surface sites using sequence conservation information only, because differences in conservation are generally negligible.

Specificity differs between interface, surface and buried positions

Next, we calculated the specificity between interacting and non-interacting groups for the three aforementioned types of positions using Sequence Harmony. In our hypothesis, the lower SH scores should be located at the interface positions. We took the complete dataset to show the specificity. Sequence Harmony detected compositional differences between interacting and non-interacting sub-groups for the three different types of positions. The overall average SH value for interface is 0.358, for other surface 0.375 and buried position 0.380 (p-value ≤ 0.05 between interface and other surface using two-sided Student's t-test). Indeed, the interface positions show the lowest SH scores. In other words, there is signal present in the specificity.

We investigated whether there was a subset of our pair-groups which could score much stronger SH signals. Two parameters were used to make different dataset selections. First is the length of the High-Scoring Segment Pair (HSP) between query homodimer and its first monomer hit as given by BLAST. The HSP length cut-off was used to exclude hits which only consisted of small local alignments from the BLAST results. Our pair groups should reflect homologous relationships between two groups of homologs: the homodimer group which conserves interaction, and the monomer group which conserves non-interaction. Since longer HSPs reflect stronger evidence for a homologous relation, we consider longer HSP hits reported by BLAST as more relevant. With increasing HSP-length, the difference between the average SH scores of the different positions grows steadily, as can be seen in Fig. 3. This is caused by the SH scores of the interaction sites decreasing more rapidly than the scores for the two other types of positions. The differences of average SH values between interface and surface become even more significant from 100 aa to 300 aa minimal HSP length (p-value ≤ 10⁻⁵ using two-sided Student's t-test). In other words, the specificity signal appears to be able to separate interaction sites out of both the other surface sites and buried residues, and better for the more significant, longer HSP pair groups. The numbers of pair groups retained after filtering on HSP-length and sequence length are reported in Additional file 1: Table S2 and S3.

The second selection parameter was the interacting (homodimer) sequence length. Additional file 1: Figure S1 shows the changes of average SH score with increasing sequence length (of the interacting protein). The figure shows that with increasing sequence length, the differences between the SH scores of the three groups remained rather stable, unlike the diverging trend observed for HSP-length selection (Fig. 3.).

Interface prediction depends on HSP length and % ID filtering

In light of the observed trends in differences of SH scores, we tested if the specificity signal could be used to predict interaction sites. Predictions were validated against interaction site annotations obtained from the PISA database (see Methods for details). For all pair groups, the Area Under the Curve (AUC) of Receiver Operator Characteristic (ROC) plots was calculated as a measure of performance.

Figure 4 presents the average AUC versus increasing HSP-length. As can be seen, the performance is rising with increasing HSP length, as expected based on Fig. 3. Moreover, the performance also increases with lower sequence identity cut-offs. For the dataset filtered at 50-99 % ID, the average AUC reaches a maximum of around 0.64-0.65, a feasible range for prediction purposes. Figure 4B shows that the AUC values may be increased further by combining the sequence identity cut-off with a minimum sequence length. Interestingly, the <100 % ID set performs worse than the others at higher HSP length cut-offs. This set contains many single-residue mutants that appear to bias the scoring. The highest average AUC increased to about 0.7 with the dataset selected by sequence identity <70-98 %, sequence length >400 aa and HSP length between interacting and non-interacting group >200 aa. The AUC with the datasets filtered at a 40 % identity cut-off, or at HSP-length >250 aa are not shown because too few pair-groups remain (see also Additional file 1: Table S3).

Thus sequence length and HSP length were used in combination to select a final dataset. In this part we select the dataset filtered at 80 % ID for performance reasons. Figure 5 shows the ROC plots using different HSP cut-off values and a minimal length of 400 aa. All ROC plots show a better-than-random performance of prediction. Also with minimal sequence length 400 aa, the average AUC was improving when the HSP length was increasing. Here, we could achieve an AUC of 0.7 using an HSP length of 175 aa or 200 aa.

Predicting the interface that stabilizes the ligand binding regions in a phosphatase and a kinase family

To illustrate the impact of accurate interface prediction, we here show details of two protein families with the highest AUC in ROC performance in our dataset, indicating the specificity signal strongly correlates with the interface region.

The first example is chronophin, a cofilin-activating phosphatase of the haloacid dehalogenase (HAD) super-families, directly dephosphorylates cofilin with high specificity [17]. C2a-type HAD Hydrolases, usually show interaction between two chains, froming a homooligomer chronophin. This phosphatase subfamily is only active as oligomer, because the interactions are crucial for substrate specificity loop positioning, while other subfamilies may not all require oligomerization [25]. We constructed pair groups using 3QGM as query sequence following the scheme described in methods. The interacting group contains all four dimeric C2a-type HAD Hydrolases (3QGM, 1ZJJ, 1YV9, 2OYC), and the non-interacting groups contains all 11 monomeric C1-type HAD Hydrolase group (2NYV, 2HSZ, 2HI0, 2AH5, 4EX6, 3MC1, 3D6J, 3KBB, 3KZX, 2HDO, 3SD7). The AUC of ROC for prediction is 0.756 (Additional file 1: Figure S4). With the SH value cut-off <0.2, our method identifies 85 positions, including 21 binding sites out of all 31 interface sites: Recall 67.7 %, FPR 27.7 %, Precision 24.7 %. Fig. 6 shows the interface and predicted interface sites in the structure. The binding site is clearly delineated by the amino acids correctly predicted to be interacting.

The second example is Amino-imidazole riboside kinase (1TZ6), which has a homodimeric structure with one active site per monomer. The active site is covered by a lid which is supposed to be a morphological marker for evolution within the ribokinase superfamily [3, 56]. The homodimeric structure is formed through lid-to-lid interactions [3, 56]. The query monomer sequence found is 2ABS, adenosine kinase, another member of the ribokinase superfamily which can be active as a monomer. We identify 53 positions and detect 13 interacting sites out of all 24 binding sites below the SH value cut-off 0.2: Recall 54.2 %, FPR 14.7 %, Precision 24.5 %. The AUC of ROC reaches 0.746 (Additional file 1: Figure S4.)

Interaction site prediction using sequence information alone remains a challenging problem. It is particularly important in the context of increasing protein sequence data and given the relative paucity of structure information which always needs expensive and time-consuming experiments. We demonstrate that sequence specificity information from interacting proteins and their non-interacting homologs is able to detect interaction sites. To the best of our knowledge, this is the first time that predicting interaction sites using subfamily specificity by including non-interacting information is performed at a scale beyond a few protein families.

Our results show that prediction is well beyond random: The SH signal is able to obtain ROC values greater than 0.6 AUC. The AUC increases to 0.7 if the dataset is filtered on sequence length and HSP-length between interacting and non-interacting groups. HSPs are formed by BLAST local alignments based upon residue similarity, thus residues should be similar. However, this goes against looking for differences between two subgroups, which is the hallmark of SH scoring. Since shorter HSPs should have more similar alignment positions to get a high-enough BLAST score, longer HSPs may contain less similar positions; i.e., that differ more between the subgroups, here implying interface residues that confer interaction specificity. Additionally, alignments comprising shorter HSPs contain larger regions of higher divergence, which will likely lead to more false positive predictions.

For our analysis, on average only 25 homodimer and 14 monomer sequences were used in a group. Furthermore, we do not lose prediction performance with low numbers of sequences, above the minimum of 5 required in our analysis. This is a vast improvement to covariance-based methods that require an estimated five times more sequences than the length of the alignment to detect interacting sites. Thus, this opens up the possibility of deriving interaction signals from genomes with little sequence data available, or from sparsely sampled protein families.

We observe that conservation of interaction sites is indistinguishable from other surface sites (Fig. 2), which corroborates observations by others (e.g., [7, 28]). This also helps explain why predicting interaction sites only using sequence conservation information still remains a very difficult problem. We also observe low correlation (R = 0.22) between AUC of ROC plots and sequence length of query homodimers, which suggests that the size of protein and its interface are unimportant factors in interaction. That is similar to what Dhole et al. [10] reported recently.

The definition of 'non-interacting' in our manuscript is not that the monomer can not interact with all the other proteins but the monomer loses the interaction with another monomer (itself). It is reported that, the same protein involved in different interaction might have different binding sites [35]. If the monomer in our non-interacting group also binds to other proteins, the interaction sites might be different. Then, the sequence specificity between the interacting sub-group (homodimer) and monomer sub-group can still be used to pinpoint the homodimer interaction sites.

Our dataset is obtained from PISA and PDB and might include any bias which the PDB has. To test this, we map the homodimer query sequences from each group in Test-set 1 (both Complete and Selected) onto CATH superfamilies and calculate the overall performance of our method (Average area Under Curve, AUC) for each superfamily. Our Test-set 1 covers all four main classes, 62.5 % (25/40) architectures, 11.5 % (158/1375) topologies, 8 % (219/2738) superfamilies in CATH while only comprising 1593 out of 69058 (2.3 %) total proteins in CATH. Our method does not only predict superfamiles which are enriched in PDB but also those for which few (homodimer) structures are solved. For the selected dataset, the sequences map to 7 superfamilies. Interestingly, we also selected one superfamily which our method can not predict well. The results of this can be found in Additional file 1: Figure S7 and S8.

In order to show that we exploit useful information from sequence (sequence specificity between interacting and non-interacting homologs) to predict interaction sites, we compared our results with other sequence-based methods. A clean comparison is difficult, since methods may rely on a particular way in which the data is prepared (e.g., particular database, or % identity filtering), and, in our case the inclusion of a non-interacting group. Figure 7 shows the prediction performance of other interacting site prediction methods using their own dataset from a recent benchmark [10]. Note that, due to the lack of a non-interacting group, we cannot apply our method to the Dhole et al. benchmark dataset. For comparison, performance of our method and SPPIDER using our dataset are also shown (cf. Table 1). On our dataset we achieve similar precision, but considerable higher recall than SPPIDER, and we get similar precision and recall compared to the best ones reported by Dhole et al. [10]. Overall, we conclude that our method performs well compared to counterparts, although alternative methods may perform better on individual datasets.

We tested the robustness of our interface prediction to changes in the selection of homodimer and monomer sequences. For selecting sequences into either subgroup, homodimer or monomer, we use an effective e-value cut-off which is based on the first occurrence of a sequence of the other class in the BLAST results; i.e. the first monomer found for a homodimer query determines the e-value cut-off for creating this homodimer group (and vice versa). By changing this e-value cut-off we changed the selection of sequences and thereby test the robustness. We created three variations: 1) By taking the lowest of the two e-values for the homodimer and monomer groups in a pair-group, we obtain a stricter selection; this tests the effect of missing data. 2) By relaxing the threshold up to the point where 20 % of the other class was mixed in, we obtain a looser selection; this tests the effects of polluting the data. 3) By taking a fixed, strict but arbitrary e-value cut-off of 10E-10, we test the case where no complete annotation of interaction status of the sequences is available. Figure 8 summarizes the results for these three variations. For comparison, also the original selections are shown, and p-values for the comparisons are listed in Additional file 1: Table S7. A stricter threshold (variation 1) does not affect results, and neither does a fixed 10E-10 cut-off (variation 3). The polluted or 'mixed' data (variation 2) does yield a lower AUC, but only slightly so. This means our method for predicting interfaces is stable with respect to missing data, incorrect annotations, as well as to a particular choice of e-value cut-off. This opens up the possibility of extending the method to use a single e-value cut-off to delineate both groups, which means a full annotation of interaction status of all proteins may not be required.

Although our results demonstrate an advance in predicting interaction sites from sequence information alone, it is clear that, like other approaches, our approach also holds some limitations. First, detection of functionally conserved homologs is still a difficult problem. Thus, care should be taken in selecting homologous proteins that are likely to conserve interaction. Initially, our procedure was based on ortholog clustering databases, such as COG and OrthoMCL, but unfortunately, interaction appears to be an ill conserved property within these orthologous groups [28]. This introduces considerable noise into the specificity signal [14]. Second, for our prediction we need homologous sets of interacting and non-interacting groups. Not for all interacting proteins will one be able to identify non-interacting homologs. Proteins interacting with different partners might be used instead [14, 39], but that remains to be investigated further.

We have shown it is possible to predict interaction sites out of all residues by combining sequence and group specificity information. When used as a prediction method in its current form, on homodimer versus monomer data the Sequence Harmony specificity signal yields similar precision as other signals but may obtain higher coverage.