May 8, 2007 - evolutionary constraints gene expression binding interfaces. Proteins ...... We thank Casey Bergman, Steven Littler, and Simon Hubbard for.
Specificity in protein interactions and its relationship with sequence diversity and coevolution Luke Hakes, Simon C. Lovell, Stephen G. Oliver, and David L. Robertson* Faculty of Life Sciences, University of Manchester, Michael Smith Building, Oxford Road, Manchester, M13 9PT, United Kingdom Edited by Philip P. Green, University of Washington School of Medicine, Seattle, WA, and approved March 12, 2007 (received for review November 10, 2006)
evolutionary constraints 兩 gene expression 兩 binding interfaces
P
roteins rarely act alone; rather, they carry out their activities through a multitude of interactions with other proteins or molecules, either pairwise or in association with multiple subunits, as is the case for protein complexes (1). In yeast, the organism in which the most comprehensive studies have been done (2, 3), almost all proteins take part in some kind of interaction (4). Protein–protein interactions (PPIs) are thus fundamental to almost all biological processes and are of particular importance in relating genotype to phenotype, not only at the molecular level, but at all levels up to and including that of the organism (1). The availability of complete genome sequences and recent breakthroughs in technology to identify PPIs have allowed a shift toward a systems-level approach to biology. This focuses on interactions rather than individual components, permitting genes, their regulation, and their products to be studied as interacting components of networks and complexes rather than in isolation or as linear pathways (5). In this network-based view of function, the pleiotropic nature of the majority of proteins is central, and the roles of proteins are considered in the context of functional modules (6, 7). This interaction-centric approach is a key advance, permitting an understanding of the components that make up these networks. Individual PPIs are generally mediated through localized sites or ‘‘patches’’ situated on the exposed surface of the protein. These patches have been found to exhibit particular characteristics (8–14) in terms of their solvation potential, interface residue propensity, hydrophobicity, planarity, protrusion, and accessible surface area. These properties are common to all interaction surfaces but may also contain some of the informawww.pnas.org兾cgi兾doi兾10.1073兾pnas.0609962104
tion needed to ensure binding specificity. Because the specific binding characteristics of these patches are inherently dependent on their underlying amino acid sequence, substitutions in the interface of one member of an interacting pair may be compensated by complementary substitutions in the interface of the other. This results in coevolution of the two binding partners. By ‘‘coevolution,’’ we are referring to the coadaptation of interacting proteins. That is, substitution of an amino acid residue in the interface of one protein will select for the coordinated mutation and subsequent substitution of an amino acid in the interface of its binding partner (i.e., compensatory change). These changes are a direct result of maintaining optimal structural and functional integrity. Correlated evolution can arise through other mechanisms, such as shared evolutionary history (phylogenetic covariation) (15). However, we distinguish correlated evolution, which arises because of generic external factors acting independently on two proteins that bind, from coevolution, whereby change in one binding partner has a direct influence on change in the other binding partner. We assume that the majority of such changes will be mediated at or near the binding surface (16). Two broad approaches have been developed to identify the coevolution of binding partners with the intent of predicting probable PPIs: (i) methods that are site-specific and (ii) those based on the analysis of the entire amino acid sequence. Sitespecific methods often explicitly seek to model compensating mutations in interfaces, mapping sites with correlated mutations back to protein structures (17). Far more widely used is the ‘‘entire-sequence’’ (or ‘‘mirrortree’’) approach (18–26), in which pairwise-distance matrices derived from amino acid sequence alignments are compared, and the detection of significant correlations is used to infer probable coevolution. This approach has the advantage that phylogeny can be explicitly taken into account to reduce the detection of false positives (24, 25). Several authors have discussed the likely origins of the correlation signals detected by the mirrortree approach (and related entire-sequence methods), with the majority suggesting that such signals arise through mutations in one binding partner being compensated for by complementary mutations in the other (18, 20, 22, 23, 26). For this to be true, the surfaces of the proteins must be coevolving, and the residues involved in this process must comprise a relatively large fraction of the protein sequence. However, for the majority of proteins, the proportion of residues that constitute the binding interface is known to be relatively small (11). Moreover, the coevolutionary signal is relatively weak and may not be detectable in binding partners that make transient interactions (26). Given these problems, how can such Author contributions: L.H. and S.C.L. contributed equally to this work; L.H., S.C.L., S.G.O., and D.L.R. designed research; L.H. and S.C.L. performed research; L.H., S.C.L., and D.L.R. analyzed data; and L.H., S.C.L., S.G.O., and D.L.R. wrote the paper. The authors declare no conflict of interest. This article is a PNAS Direct Submission. Abbreviations: PPI, protein–protein interaction; ROC, receiver operator characteristic. *To whom correspondence should be addressed: E-mail: david.robertson@manchester. ac.uk. © 2007 by The National Academy of Sciences of the USA
PNAS 兩 May 8, 2007 兩 vol. 104 兩 no. 19 兩 7999 – 8004
EVOLUTION
Studies of interacting proteins have found correlated evolution of the sequences of binding partners, apparently as a result of compensating mutations to maintain specificity (i.e., molecular coevolution). Here, we analyze the coevolution of interacting proteins in yeast and demonstrate correlated evolution of binding partners in eukaryotes. Detailed investigation of this apparent coevolution, focusing on the proteins’ surface and binding interface, surprisingly leads to no improvement in the correlation. We conclude that true coevolution, as characterized by compensatory mutations between binding partners, is unlikely to be chiefly responsible for the apparent correlated evolution. We postulate that the correlation between sequence alignments is simply due to interacting proteins being subject to similar constraints on their evolutionary rate. Because gene expression has a strong influence on evolutionary rate, and interacting proteins will tend to have similar levels of expression, we investigated this particular constraint. We found that the absolute expression level outperformed correlated evolution for predicting interacting protein partners. A correlation between sequence alignments could also be identified not only between pairs of proteins that physically interact but also between those that are merely functionally related (i.e., within the same protein complex). This indicates that the observed correlated evolution of interacting proteins is due to similar constraints on evolutionary rate and not coevolution.
Results Assessment of Coevolution. We tested the ability of the mirrortree
approach (19) to detect PPIs in a subset of the yeast proteome that was composed of proteins involved in at least one physical interaction within a protein complex with a known crystal structure. Our final data set comprised 94 proteins from 32 complexes. For each pairing, we used three multiple alignments representing the whole sequence, the surface regions, or the interaction interface. The accuracy of the mirrortree prediction method using each alignment was assessed as described by Pazos et al. (24). Following calculation of interaction scores for all pairings for each protein, the resulting list was sorted in descending order, and the percentage of proteins in the sorted list that scored higher than the first true interacting pair was calculated. Interactions were defined as ‘‘true’’ if they were either identified from three-dimensional structures of complexes or were found in a subset of BioGRID (36) (see Methods). Thus, for the lower the percentage of false positives, the better the prediction. In accordance with previous observations using prokaryotic data (19), the majority of the true interacting partners are found toward the top of the list. The average percentage of false positives above the first true interacting pair for the wholesequence alignment was 20.1% (median, 10.2%). Somewhat surprisingly, we found no improvement when only surface residues were used (20.9% false positives; median, 11.6%) or when we limited the data set to those residues directly associated with the interaction interface (22.7% false positives; median, 10.3%). All three of these values are comparable with those observed in the previous work of Pazos et al. (19), who demonstrated a false-positive discovery rate of 23.4% when they applied this method using whole-sequence alignments to predict PPIs in Escherichia coli. To evaluate the performance of the mirrortree approach in each region of the alignment, a three-way comparison was 8000 兩 www.pnas.org兾cgi兾doi兾10.1073兾pnas.0609962104
Table 1. Comparison of the performance of mirrortree when different regions of the proteins are used Whole sequence
Surface Interface
Surface
P
ROC (z)
P
ROC (z)
0.576 0.078
1.51 2.88
0.101
2.47
P values are given for the comparison using the method of Pazos and coworkers (19) (one-sided sign test). Critical ratio (z) given for the comparison of ROC areas calculated using the method of Delong et al. (39).
undertaken. Relative performances were assessed by using two different methods. In the first method, the average rank of the positive cases on the ordered list for each protein in turn (24) was calculated. A one-sided sign test was used to determine statistical significance (37) with the null hypothesis being that there was no difference between the performances of each alignment type. The results (summarized in Table 1) reveal no significant differences among the three, with the greatest difference in performance observed when the performance of the whole alignment is compared with that of the interaction interface. As a second method to determine performance, the area under a receiver operator characteristic (ROC) curve was calculated (Fig. 1). From this we were able to obtain a more formal quantification of the performance of the mirrortree method with each alignment. The value of the ROC area represents the probability that a randomly chosen protein pair selected from our test set is correctly identified as either interacting or noninteracting (38). Calculation of this area by using the method of DeLong et al. (39) gives us values of 0.599 (2 ⫽ 2.8 ⫻ 10⫺4) for the method using the complete protein sequence, 0.560 (2 ⫽ 2.7 ⫻ 10⫺4) for the surface region of the protein, and 0.548 (2 ⫽ 3.4 ⫻ 10⫺4) using only the interfacial residues. Each method was then compared directly by calculating the critical ratio (z) between different modalities (sequence sets), defined as: z⫽
AUC1 ⫺ AUC2
冑Var关AUC1兴 ⫹ Var关AUC2兴 ⫺ 2Covar关AUC1, AUC2兴 , [1]
where AUC1 and Var[AUC1] refer to the observed area under the curve and its estimated variance of the ROC area associated with modality 1, AUC2 and Var[AUC2] refer to the corresponding quantities for modality 2, and where 2Covar[AUC1, AUC2] represents the covariance between variance estimates of the two areas (39). The critical ratio can be converted into a probability
Sensitivity
a localized evolutionary signal be detected by analysis of the whole-protein sequence? To better understand the processes underlying the apparent coevolution of binding partners, and in an attempt to improve the existing methodology, we restricted our analysis of correlated evolution to the surface residues of proteins of a known structure. Because interactions should be mediated through surface patches, this approach should increase the proportion of interacting residues and, thereby, improve the signal. We found that this was not the case. Furthermore, by restricting our analysis to those residues involved in binding interfaces, we found that the degree of correlated evolution actually decreases. This suggests that correlated evolution over the length of the protein sequence is separate from and unrelated to any true coevolution that occurs in the binding interface. To test this hypothesis, we determined the extent to which protein divergence at the sequence level in duplicates generated by a whole-genome duplication event (27–30) corresponds to changes in function as measured by numbers of shared interactions. We found no simple relationship between the sequence divergence of the protein products of a pair of duplicated genes and the proportion of the interacting partners they share. Correlated evolution must, therefore, be due to interacting proteins being subjected to similar evolutionary constraints. Because the level of gene expression is one of the main determinants of protein sequences’ rate of evolution (31–35), we investigated the role of gene expression on the correlation of interacting proteins. Our results indicate that the correlation between sequence alignments of interacting proteins is most likely due to members of the same functional module having similar evolutionary constraints and, consequently, similar, and thus correlated, rates of evolution.
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1-Specificity
Fig. 1. ROC plot showing the performance of the mirrortree method when different parts of the protein sequence are used in the calculation. Whole sequence is shown in blue, surface region of the protein is shown in red, and protein interaction interface is shown in yellow. Predictive performance of a method using the abundance ratio (green) between pairs of proteins is also illustrated.
Hakes et al.
sequence divergence as measured over the whole length of a protein and the protein partners bound by that protein, we examined pairs of paralogs. Paralogous pairs will have had identical protein sequences at the point of duplication and so will have had the same binding surfaces and specificity. Over time, their sequences will diverge, and their binding specificities may, as a consequence, change. The set of whole-genome duplicates previously identified (27) thus constituted an excellent data set with which to test the hypothesis that paralogous duplicates in PPI networks have a tendency to exhibit a correlation between sequence divergence and functional similarity. To do this, we determined the relationship between sequence divergence and the ratio of interactions in common between pairs of paralogs. We used this ratio rather then the raw number of shared interaction partners because it adjusts for possible errors introduced by the presence of proteins with a large number of interactions. Given the number of false-positive interactions identified by experimentation (40), highly connected proteins would, simply by chance, be more likely to share interactions with their duplicate partner. Fig. 2 shows the proportion of interactions that are in common between duplicate pairs and how this varies with sequence divergence. A value of 1 on the y axis indicates that each member of a duplicate pair has exactly the same binding specificity, whereas a value of 0 indicates that they have no interactions in common. Intermediate values indicate that some of the interactions have been gained, lost, or altered. It is evident (i) that there was a substantial change in the specificity of PPIs over evolutionary time as measured by the number of nonsynonymous substitutions per site (dN) (Fig. 2 A), and (ii) that this change in specificity was not correlated with sequence divergence, at least when this divergence was calculated over the whole length of the protein sequence. This result appears to be independent of the measure of sequence divergence, with similar plots obtained when the number of synonymous substitutions per site (dS) values are estimated (data not shown) or by consideration of branch lengths in three-way phylogenetic trees (Fig. 2B). Furthermore, limiting the analysis to the subset of duplicate gene pairs in which one member has been identified as having undergone accelerated evolution, presumably as a result of functional change (27), again indicated no correlation between the two metrics (shown in red in Fig. 2 A and B). Analysis of Gene Expression. Because expression level is known to particularly influence the rate of sequence evolution (31–35), we were interested to see whether it would be possible to predict PPIs using abundance measurements alone and how those predictions would compare with those generated by using the mirrortree method. Thus, we repeated our earlier analysis substituting the mirrortree score for our newly calculated abundance ratio and quantified the results by using both the ROC Hakes et al.
Shared int. ratio
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
0
0.2
0.4
0.6
0.8
1
1.2
dN
B
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
1
2
3
4 5 Branch Ratio
6
7
8
Fig. 2. Plot of sequence divergence versus shared interactions between paralogs. The x axis represents sequence divergence, measured by either dN (A) or branch ratio (B), and the y axis represents shared interactions, measured by the ratio of interactions in common between pairs of paralogs. Pairs in which one paralog has been identified as undergoing accelerated evolution (27) are shown in red.
method and that of Pazos et al. (24). As can be seen from the curves presented on the ROC plot (Fig. 1), the abundance ratio significantly outperformed the mirrortree method regardless of which part of the sequence was used. Moreover, when we assessed the performance of the abundance ratio measure by using the method of Pazos et al. (24), we found that the measure provides a marginal improvement, with the average percentage of false positives above the first true interacting pair being 18.7% (compared with the whole-sequence mirrortree at 20.9%). If the correlated evolution between binding partners was mainly a result of gene expression influencing the rate of sequence evolution in similar ways, this would presumably be a direct result of interacting members of protein complexes tending to exhibit similar levels of gene expression. By calculating the ratio between the mRNA abundance levels for the parent genes encoding the proteins, we were able to determine the extent to which the expression levels of the proteins in the pairings are similar. For interacting proteins, we found an average mRNA abundance level of 3.3, indicating that these pairs had parent genes whose absolute expression levels are much more similar than those of random noninteracting protein pairings (Fig. 3). This indicates that pairs of interacting proteins have mRNA abundance levels that are more similar (P ⬍ 1 ⫻ 10⫺4) than randomly selected pairs of noninteracting proteins. To confirm the link between expression levels and correlated sequence evolution, we tested whether we could detect significant correlated evolution among proteins that do not interact directly but which are found within the same protein complex. We first determined whether we could detect a significant mirrortree signal within this subset of protein pairings. The mean interaction score was calculated for both interacting proteins (239 pairs) and those noninteracting proteins present within the same protein complex (91 pairs). For negative controls, we calculated the interaction score for ‘‘unassociated’’ pairs, i.e., either 239 random noninteracting pairs or 91 random noninterPNAS 兩 May 8, 2007 兩 vol. 104 兩 no. 19 兩 8001
EVOLUTION
Relationship Between the Divergence of Protein Sequences and That of Protein Interactions. To understand the relationship between
A
Shared int. ratio
by referring to tables of the normal distribution, with values of z ⱖ 1.96 taken as evidence that the ‘‘true’’ ROC areas are different. The results of this analysis are summarized alongside those for the sign test in Table 1. In accordance with the ranking analysis comparison of the method, using either the whole sequence or surface regions reveals no significant difference between the two methods. However, when the calculation was performed using residues within the interaction interface, we observed a small but significant reduction in prediction accuracy. We conclude, therefore, that the correlated evolution between binding partners detected by the mirrortree method was not mediated primarily through binding surfaces. Consequently, correlated evolution is unlikely to have arisen from coevolving, compensating mutations between binding proteins.
Osbser vation
900 800 700 600 500 400 300 200 100 0
2
4
6 8 10 12 14 16 18 20 Average abundance ratio
Fig. 3. Average ratio between mRNA abundance level for the protein pairings. Physically interacting protein pairs are indicated by the red arrow, whereas noninteracting proteins present within the same protein complex are marked by the blue arrow. The distributions are for 10,000 random samples of noninteracting proteins (237 pairs, red) and noninteracting proteins taken from different protein complexes (91 pairs, blue).
acting pairs in which each member resides in a different complex. Permutation testing was used to determine the significance of the result. We found that the average interaction score was significantly greater (P ⬍ 1 ⫻ 10⫺4) for both directly interacting proteins and for noninteracting proteins found in the same complex (Fig. 4) than for proteins in unassociated pairs, indicating that correlated evolution was occurring in both of these sets of proteins. Interestingly, noninteracting pairings from the same complex had a higher average interaction score than interacting proteins in general (0.55 versus 0.68). When we also analyzed gene expression (Fig. 3), we found that, just like interacting proteins, noninteracting proteins present within the same protein complex had more similar absolute mRNA abundance levels than random pairs of noninteracting proteins (P ⬍ 1 ⫻ 10⫺4). Discussion It has been suggested that approaches for detecting PPIs based on the entire sequence alignments are detecting coevolution (18, 20–22, 24, 25). If this were the case, we would expect to improve the predictive power of these methods by limiting the analysis to use only those parts of the sequence that comprise the protein surface because these are enriched for interface residues. Furthermore, by using the interface region alone, we would expect to see an additional improvement. In fact, we found no improvement when the method was applied exclusively to either region. We found that the degree of correlated evolution was similar when either the whole sequence or just the surface residues were used and was less when the analysis was limited to residues within the interaction interface. Most conclusively, noninteracting pro-
Osbservation
3500 3000 2500 2000 1500 1000 500 0 0.39 0.44 0.49 0.54 0.59 0.64 0.69
Average interaction score Fig. 4. Average interaction scores. Interacting protein pairs are marked by a red arrow; noninteracting pairs present within the same protein complex are marked by a blue arrow. The distributions represent the average interaction scores for 10,000 random samples of noninteracting pairs (red) and noninteracting pairs present from different protein complexes (blue). 8002 兩 www.pnas.org兾cgi兾doi兾10.1073兾pnas.0609962104
teins found in the same complex had equivalent or higher interaction scores when compared with interacting proteins. Because the degree of correlated evolution between binding partners is higher over the entire protein sequence compared with the binding interface, we investigated the relationship between sequence divergence and changes in specificity. Other authors have used accelerated evolution of sequence as direct evidence for the evolution of function (27, 41). Therefore, because function is mediated (at least in part) through binding partners (42), some correlation between these metrics should be evident. However, we found that there is no correlation between the shared interaction ratio for a duplicate pair and the divergence of their amino acid sequences. This is true whether sequence divergence is measured by the number of synonymous or nonsynonymous substitutions per site (dS and dN, respectively), or by assessment of duplicate pair branch lengths within three-way phylogenetic trees, including a Kluyveromyces waltii ortholog (Fig. 2). Furthermore, when we limited the analysis to only those duplicates in which one member has undergone accelerated evolution, which others have assumed indicates the acquisition of a new function (27, 41), we still observed no correlation. These results indicate that approaches based on the entire protein sequence are not detecting true coevolution. Perhaps this is unsurprising given that, within our test cases, the interface region comprises only ⬇13% of the total protein sequence on average. Even taking into account residues subtending the binding interface, the proportion of sequence involved in the interaction will still only comprise a small fraction of the entire sequence length. Our result highlights the misinterpretation of the mechanisms that underlie the processes of coevolution and correlated evolution. Clearly, in some cases true coevolution must still be occurring to maintain PPIs (26). However, this will be mediated through the localized sites or patches, situated on the exposed surface of the protein and across which binding is mediated. Crucially, correlations based on entire sequences are not detecting changes that are required to maintain an interaction. Because coevolution cannot account for the evolutionary signal that mirrortree and related methods are using to predict PPIs, external factors must be constraining the evolutionary rates of interacting proteins. Potential factors include functional importance as measured by functional category, position in the PPI network, pleiotropy, expression level, essentiality or dispensability, etc. For an extensive discussion of this topic, see the recent review by Pal et al. (43). Interestingly, it has been observed that gene expression is one of the strongest determinants of the rate of sequence evolution (31–35). Thus, we postulated that approaches based on whole-protein sequences are detecting similarity of evolutionary rates in part as a result of similarity in gene-expression levels. By using quantitative transcriptome data directly to predict PPIs, we found that this approach significantly outperformed the mirrortree method (Fig. 1). Investigating this further, we found that interacting proteins have similar levels of gene expression (Fig. 3), confirming that correlated sequence evolution is most probably due to interacting proteins being constrained in similar ways and, consequently, having similar rates of evolution across their entire sequences. Furthermore, both interacting and noninteracting proteins in protein complexes have expression levels that are much more similar than random noninteracting proteins (Fig. 3). That correlated sequence divergence between interacting proteins can be extended to members of the same complex that do not come into contact [also see Chen and Dokholyan (44)], implies that gene expression is exerting its influence in a similar manner on members of the same protein complex. Such a tendency to maintain balanced levels of gene expression is presumably a result of a requirement for interacting proteins to be present in similar quantities, for Hakes et al.
Methods Assessing Coevolution Between Interacting Protein Pairs. All protein
structures deposited within the Protein Quaternary Structure database (47) were downloaded and triaged to include only complexes composed of unique chains that had experimentally determined homologs in yeast. By extracting the protein sequences of every protomer from each complex and searching for it against the yeast genome by using BLAST (48) (cut-off, E ⱕ 1 ⫻ 10⫺10, fraction of conserved residues ⱖ35%), we were able to identify and make use of protein complexes determined in other species that are structural homologs of those found in S. cerevisiae. Actual physical interactions between protomers within the complexes were identified by computational analysis of the crystal structure by using an empirically equivalent algorithm to the full atom contact method used by Gong et al. (49); an interaction between two protomers was defined as the Hakes et al.
occurrence of at least five instances in which C␣ atoms within the chains come within 7.5 Å of each other. This set was included in the ‘‘true positives’’ set of interacting proteins. A set of negative controls was created from all possible pairings between those proteins. The mirrortree approach (19) was used to computationally predict PPIs between all pairings. Each protein sequence was searched against the ‘‘nonredundant’’ database from the National Center for Biotechnology Information (cut-off, E ⫽ 1 ⫻ 10⫺5) to identify homologous sequences and a set of proteins containing only the top hit from each species was created. These sets were then reduced so that they only included sequences from species common to both homologues. Sets containing ⱖ12 sequences were then aligned by using MUSCLE (50). After imposing the limitations discussed above, our final data set comprised 4,708 protein pairings, of which 239 were identified as being true-positive interactions either through structural analysis (105 pairings) or their presence in the filtered BioGRID data set (134 pairings). In addition to analyzing the full-length alignments, columns were taken from these alignments representing both the surface and interface regions of each of the proteins. Surface regions were identified by calculating residue accessibilities by using the rolling probe algorithm (51) and were defined as including those residues with ⱖ20% of their surface area exposed (for both the main and side chains). The interaction interface was identified by calculating the solvent-accessible surface area for residues in the protomer alone and in the complete complex. If the accessible surface area of a residue was different in the isolated protomer and the complex, it was defined as being in the interface. For the negative controls, in instances where a protein contained multiple interaction interfaces, one of the interfaces was selected at random. Distance matrices were computed for each of the alignments by using PROTDIST using the Jones–Taylor–Thornton model for amino acid replacements per site (52) and for the interaction score between the two seed proteins calculated by determining the linear correlation coefficient between the two matrices. Thus, for two proteins X and Y with n species in common in their multiple sequence alignment, the interaction score rXY can be calculated as follows:
冘冘
n⫺1
rXY⫽
冑冘 冘
n
共dX ij ⫺ dX兲䡠共dY ij ⫺ dY兲
i⫽1 j⫽i⫹1
n⫺1
n
i⫽1 j⫽i⫹1
冑冘 冘 n⫺1
共dX ij ⫺ dX兲 2䡠
,
[2]
n
共dY ij ⫺ dY兲 2
i⫽1 j⫽i⫹1
where dXij is the distance between species i and j in the tree of protein X, and dYij is the distance of species i and j in the tree of protein Y. Analysis of Shared Protein Interactions in the Context of Sequence Divergence. Genome duplicates were taken as those described by
Kellis et al. (27) [similar results were obtained when duplicates, as described in the yeast gene order browser, were used (53)], and the number of synonymous and nonsynonymous substitutions per site (dS and dN, respectively) were calculated by using the method of Yang and Nielsen (54) as implemented in PAML (55). Protein sequence divergence between the duplicates was assessed by studying the evolution of each member of the pair with respect to the other and the nonduplicated K. waltii ortholog as described previously (27). The value bdiff (measured in percent) was defined as the difference between the branch lengths of the S. cerevisiae branches for each member of the duplicate pair relative to the K. waltii branch in a three-way phylogenetic tree. A larger bdiff represents a greater difference PNAS 兩 May 8, 2007 兩 vol. 104 兩 no. 19 兩 8003
EVOLUTION
example, to maintain the stoichiometry within protein complexes (45). Gene expression levels, measured indirectly by the codon adaptation index, have previously been found to be correlated between interacting proteins. Some authors have taken this to indicate coevolution for gene expression of interacting proteins (46). Others have also found that proteins within the same functional modules evolve at more similar rates and exhibit coevolution for expression levels (44). Whether gene expression is really coevolving (46) or not (44) remains to be seen. However, our results would suggest that similarity of gene expression level, although not the only determinant of evolutionary rate (43), is very probably one of the main determinants that enables whole sequence methods of coevolution detection to function. Our findings do not invalidate the use of correlated evolution to detect interacting proteins. Indeed, we have demonstrated that the mirrortree method, previously used to predict PPIs in prokaryotes, can be successfully applied to a eukaryotic system. Moreover, when the ranking method of Pazos et al. (24) is used as a metric, we obtain a similar level of accuracy to that obtained (19) in prokaryotes. This is particularly encouraging in light of the prevailing assumption that this method cannot be applied to interaction data sets from organisms such as Saccharomyces cerevisiae because of the multidomain nature and low complexity regions of eukaryotic genes (24). However, this apparent success should be treated with some caution because inspection of ROC plots for the data reveals a significant decrease in the overall predictive performance of the method compared with that observed in prokaryotic studies. ROC analysis provides us with a more formal measure of performance for this type of analysis. This result therefore indicates that (i) the ranking method previously used by Pazos et al. (24) is a poor indicator of the overall performance of the mirrortree method as a predictor of PPIs, and (ii) although it is possible to apply the method to eukaryotic data sets, it does not perform as well as when applied to prokaryotic data. In addition, absolute levels of gene expression perform better than comparison of sequence diversity (Fig. 1). It will be beneficial to integrate other predictors of rates of evolution (43) into a more extensive model for the prediction of PPIs. In conclusion, the results described in this study help to improve our understanding of the relationship between sequence divergence and changes in binding specificity and allow us to unravel the often confused effects of coevolution and correlated evolution. The finding of no simple relationship between sequence divergence and function is unsurprising because evolutionary measures derived from the entire length of the sequence are not informative about protein specificity. Therefore, it is clear that evolution of binding specificity is (to some degree) independent of evolution over the entire sequence and is, as a consequence, more subtle than previously thought.
between the branch lengths of the S. cerevisiae proteins and, thus, indicates that one member may have undergone ‘‘accelerated evolution’’ and subsequently acquired a new function. Protein interaction data were extracted from the BioGRID database (36), and all nonphysical interactions were excluded to produce a filtered BioGRID data set. Nonphysical interactions were defined as those in which the method of detection was annotated as one of the following: synthetic lethality, dosage rescue, synthetic growth defect, synthetic rescue, epistatic miniarray profile, dosage lethality, phenotypic enhancement, phenotypic suppression, or dosage growth defect. Exclusion of these classes of interactions left a data set containing 5,124 unique proteins that participate in 31,330 interactions. For duplicate pairs in which both members were identified as interacting with at least one other protein (377 pairs), the shared interaction ratio was then calculated by using the following equation: srat ⫽
s⫻2 , n1 ⫹ n2
[3]
1. Pereira-Leal JB, Levy ED, Teichmann SA (2006) Philos Trans R Soc Lond B Biol Sci 361:507–517. 2. Gavin AC, Aloy P, Grandi P, Krause R, Boesche M, Marzioch M, Rau C, Jensen LJ, Bastuck S, Dumpelfeld B, et al. (2006) Nature 440:631–636. 3. Krogan NJ, Cagney G, Yu H, Zhong G, Guo X, Ignatchenko A, Li J, Pu S, Datta N, Tikuisis AP, et al. (2006) Nature 440:637–643. 4. Grigoriev A (2003) Nucleic Acids Res 31:4157–4161. 5. Bork P, Jensen LJ, von Mering C, Ramani AK, Lee I, Marcotte EM (2004) Curr Opin Struct Biol 14:292–299. 6. Hartwell LH, Hopfield JJ, Leibler S, Murray AW (1999) Nature 402:C47–52. 7. Spirin V, Mirny LA (2003) Proc Natl Acad Sci USA 100:12123–12128. 8. Argos P (1988) Protein Eng 2:101–113. 9. Chothia C, Janin J (1975) Nature 256:705–708. 10. Janin J, Chothia C (1990) J Biol Chem 265:16027–16030. 11. Janin J, Miller S, Chothia C (1988) J Mol Biol 204:155–164. 12. Jones S, Thornton JM (1995) Prog Biophys Mol Biol 63:31–65. 13. Jones S, Thornton JM (1996) Proc Natl Acad Sci USA 93:13–20. 14. Hoskins J, Lovell S, Blundell TL (2006) Protein Sci 15:1017–1029. 15. Fares MA, Travers SA (2006) Genetics 173:9–23. 16. Gloor GB, Martin LC, Wahl LM, Dunn SD (2005) Biochemistry 44:7156–7165. 17. Pazos F, Valencia A (2002) Proteins 47:219–227. 18. Goh CS, Bogan AA, Joachimiak M, Walther D, Cohen FE (2000) J Mol Biol 299:283–293. 19. Pazos F, Valencia A (2001) Protein Eng 14:609–614. 20. Goh CS, Cohen FE (2002) J Mol Biol 324:177–192. 21. Ramani AK, Marcotte EM (2003) J Mol Biol 327:273–284. 22. Kim WK, Bolser DM, Park JH (2004) Bioinformatics 20:1138–1150. 23. Tan SH, Zhang Z, Ng SK (2004) Nucleic Acids Res 32:W69–W72. 24. Pazos F, Ranea JA, Juan D, Sternberg MJ (2005) J Mol Biol 352:1002–1015. 25. Sato T, Yamanishi Y, Kanehisa M, Toh H (2005) Bioinformatics 21:3482–3489. 26. Mintseris J, Weng Z (2005) Proc Natl Acad Sci USA 102:10930–10935. 27. Kellis M, Birren BW, Lander ES (2004) Nature 428:617–624. 28. Dietrich FS, Voegeli S, Brachat S, Lerch A, Gates K, Steiner S, Mohr C, Pohlmann R, Luedi P, Choi S, et al. (2004) Science 304:304–307. 29. Dujon B, Sherman D, Fischer G, Durrens P, Casaregola S, Lafontaine I, De Montigny J, Marck C, Neuveglise C, Talla E, et al. (2004) Nature 430:35–44. 30. Wolfe KH, Shields DC (1997) Nature 387:708–713.
8004 兩 www.pnas.org兾cgi兾doi兾10.1073兾pnas.0609962104
where srat is the shared interaction ratio, s is the number of interactions shared between the two proteins, and n1 and n2 are the number of interactions for ORF1 and ORF2, respectively. Analysis of Gene Expression. For the analysis of mRNA abundance,
data were taken from the study of Holstege et al. (56). We defined a metric mRNA abundance level for each protein pair, defined as the ratio between the mRNA abundance levels of each of the genes encoding the proteins. For each data set under investigation, we then calculated the average mRNA abundance level for the set. Permutation testing on ⬎10,000 iterations was performed to obtain a P value from the sample-specific permutation distribution of this ratio. We thank Casey Bergman, Steven Littler, and Simon Hubbard for helpful discussions. L.H. is supported by a CASE Studentship from the Biotechnology and Biological Sciences Research Council and AstraZeneca. Work on protein interactions in D.R.’s and S.G.O.’s laboratories is supported by research grants from the Biological Sciences Research Council. 31. Agrafioti I, Swire J, Abbott J, Huntley D, Butcher S, Stumpf MP (2005) BMC Evol Biol 5:23. 32. Drummond DA, Bloom JD, Adami C, Wilke CO, Arnold FH (2005) Proc Natl Acad Sci USA 102:14338–14343. 33. Drummond DA, Raval A, Wilke CO (2006) Mol Biol Evol 23:327–337. 34. Pal C, Papp B, Hurst LD (2001) Genetics 158:927–931. 35. Duret L, Mouchiroud D (2000) Mol Biol Evol 17:68–74. 36. Stark C, Breitkreutz BJ, Reguly T, Boucher L, Breitkreutz A, Tyers M (2006) Nucleic Acids Res 34:D535–539. 37. Bland M (2000) An Introduction to Medical Statistics (Oxford Univ Press, Oxford). 38. Bamber D (1975) J Math Psychol 12:387–415. 39. DeLong ER, DeLong DM, Clarke-Pearson DL (1988) Biometrics 44:837–845. 40. von Mering C, Krause R, Snel B, Cornell M, Oliver SG, Fields S, Bork P (2002) Nature 417:399–403. 41. Wolfe KH (2004) Curr Biol 14:R392–R394. 42. Baudot A, Jacq B, Brun C (2004) Genome Biol 5:R76. 43. Pal C, Papp B, Lercher MJ (2006) Nat Rev Genet 7:337–348. 44. Chen Y, Dokholyan NV (2006) Trends Genet 22:416–419. 45. Papp B, Pal C, Hurst LD (2003) Nature 424:194–197. 46. Fraser HB, Hirsh AE, Wall DP, Eisen MB (2004) Proc Natl Acad Sci USA 101:9033–9038. 47. Henrick K, Thornton JM (1998) Trends Biochem Sci 23:358–361. 48. Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ (1997) Nucleic Acids Res 25:3389–3402. 49. Gong S, Yoon G, Jang I, Bolser D, Dafas P, Schroeder M, Choi H, Cho Y, Han K, Lee S, et al. (2005) Bioinformatics 21:2541–2543. 50. Edgar RC (2004) Nucleic Acids Res 32:1792–1797. 51. Connolly ML (1983) Science 221:709–713. 52. Felsenstein J (1989) PHYLIP, Phylogeny Inference Package, Version 3.2, Cladistics 5:164–166. 53. Byrne KP, Wolfe KH (2005) Genome Res 15:1456–1461. 54. Yang Z, Nielsen R (2000) Mol Biol Evol 17:32–43. 55. Yang Z (1997) Comput Appl Biosci 13:555–556. 56. Holstege FC, Jennings EG, Wyrick JJ, Lee TI, Hengartner CJ, Green MR, Golub TR, Lander ES, Young RA (1998) Cell 95:717–728.
Hakes et al.
