Protein protein interactions, evolutionary rate, abundance and age

Data

The following Protein-Protein interaction datasets were used in this study: DIP, MIPS, BIND, GRID, MINT, INTACT [19–24]. A self interaction was counted as one interaction against the interacting protein. Duplicates and duplicates by virtue of inversion were removed from the interactions sets.

Protein sequences for interaction sets, where possible were downloaded from the dataset's corresponding website. In all other cases they were obtained from either SGD or UniProt depending on the annotation of proteins [25, 26]. Not all interacting proteins could be assigned protein sequences.

The MIPS dataset contained two types of interactions, physical and genetic. Physical interactions are those ascertained from Y2H studies and purified complexes while genetic interactions were obtained from suppression mutation and synthetic lethality tests. These two types were treated as two different sets.

The Intact dataset housed a small set of interactions, that contained results from small to medium scale experiments as well as the results from the four large scale studies [27–30]. These were treated as individual sets.

Evolutionary distance

Best Reciprocal Hit (BRH) orthologs in each interaction set were found using the BLASTP program [31] and comparing to the entire proteome of either Mus musculus [1] or Saccharomyces paradoxus [2]. The evolutionary rate was estimated using two methods. The first [32] required us to numerically solve the equation q = [ln(1 + 2d)]/2d, where q is the proportion of identical sites between aligned sequences and d is the evolutionary distance. The second method used the gamma distance correction [7]d = α [(1 - p)^-1/α - 1] where d is the evolutionary distance between two protein sequences, p is the number of different residues, and α is the estimated gamma shape parameter, α = 1.53 [33].

Randomisation

The randomisation test was conducted by systematically selecting a protein in a dataset, and assigning it an evolutionary rate by sampling the distribution of evolutionary distances. These random rates were then plotted against the number of interactions. 100 sets were generated, their correlation coefficients were calculated and compared to the correlation coefficient of the original experimental set.

Overlap

We universalised the labelling of all protein interactors in order to overcome the use of different notations to mark proteins. This was done by matching the sequence of each protein interactor to every sequence in the yeast genome (entire genome downloaded from SGD). Only 100% sequence matches were reannotated with the GenBank id, (GI Code). Protein attrition was no more then 5% in all datasets.

Abundance

Gene expression level data was taken from the Young lab [34].

Error rates

We used three methods to assess the accuracy of the different interaction data sets.

The first was the expression profile reliability (EPR) index [14]. The EPR index is calculated using an online server [19]. An expression based distance score is calculated for all interacting protein pairs in a set.

The resulting distribution of distance scores is compared to the distance score distributions of standard interacting and noninteracting sets. The comparison yields the approximate percentage of true interactions in the set.

The second error rate indicator, the Reference Index, involves comparing each interaction against a reference set, following work done by Von Mering et al (2002). The reference set used was the DIP_Core dataset which is considered to be a good interaction set [19]. This dataset contains protein interactions that have been computationally verified or observed in more than one large-scale experiment or those that come from small scale experiments. The percentage of interactions, from the dataset of interest, present in the DIP_Core dataset is taken as an indicator of the reliability of that set.

The third estimator of error rate was the percentage of interacting proteins that shared the same subcellular localisation. In an extension to the logic that interacting proteins would share similar functional roles, it is also possible to say that they would share similar subcellular compartments [35]. The number of interactions in which both interacting partners share the same compartment is used to give a measure of error within an interaction dataset. The subcellular localisations of yeast proteins into 19 compartmental categories is available [36]. The Localisation Similarity (LS) index was the fraction of interactions in which both proteins were from the same compartmental category.

In order to quantify the effect of error rates on each dataset, we ranked the dataset according to each measure of error. 1 being the highest rank and 0 if no error measure existed for the dataset. By calculating the mean rank of each dataset we obtain a consensus measure of error.

Evolutionary excess retention

To estimate the age of a protein we calculated the Evolutionary Excess Retention (ER), previously used as a measure for conservation of a protein [18]. The ER is a value that depicts the propensity of a protein to have orthologs in other fully sequenced genomes. It should be noted that ER does not estimate the exact age of a protein and is not necessarily correct for all proteins, as it does not identify gene loss or consider horizontal gene transfer. Therefore proteins that have a high ER value are most likely to be old but proteins with a low ER value may not necessarily be new.

We measured for orthologs in S. cerevisiae and H. sapiens, D. melanogaster, C. elegans, M. musculus and A. thaliana. Orthologs were taken from the InParanoid database [37] and we only used those core pairs of each cluster that had a confidence of 100%.

Evolutionary rate

Due to nomenclature and curation errors a few proteins without sequence data remained (Table 1). This resulted in the loss of some interactions from our final analysis. The BIND dataset was particularly affected.

Figure 1A shows the correlation between evolutionary rate and number of interactions for the DIP_Core dataset. The Spearman's rank correlation coefficient returned by the DIP_Core data was -0.1682 with a P-value of 1.29e^-11 indicating the statistical significance of this weak correlation. This negative correlation suggests that proteins with a larger number of interactions tend to evolve slower.

Varying the model used to estimate evolutionary rate between Grishin's and Ota's method had little effect. The Spearman's rank correlations for all the datasets, using Grishin's method for evolutionary rate, are shown in Table 1. None of the sets returned very strong correlations. The strongest correlation was observed in the INTACT_Ho dataset, a fully experimental set obtained by complex purification. INTACT_Ito, a purely Y2H dataset shows no correlation at all (Spearman's ρ: 0.03, P-value: 0.8905). The INTACT_SMALL dataset has a very small number of interactions and it returned the second strongest correlation, however its statistical significance was low (P-value: 0.0508). Figure 2 shows all the datasets ranked by their correlation values. The five datasets with the worst correlations, (Spearman's ρ > -0.1) returned high P-Values (> 0.01).

We also estimated the evolutionary rate of proteins by finding orthologs in the species S. paradoxus. Using a more closely related species to estimate evolutionary rate resulted in a greater number of orthologs. In all but the BIND dataset, we found orthologs for over 90% of the interacting proteins. Figure 2 shows the Spearman's rank corrolations for evolutionary rate estimated from S. paradoxus orthologs using Grishin's method, against number of interactions.

Error rates and overlap

Not all error rate measures agree for particular datasets. For example, the MIPS_Genetic dataset has a high EPR Index. This would indicate a large number of true positives, yet when comparing to a reference set the overlap is only 4.08% which is a low value compared to other datasets. The corresponding localisation similarity is at 27%, which is on the lower end of the LS index spectrum.

The INTACT_SMALL dataset possessed a strong correlation between connectivity when cerevisiae-musculus orthologs were used yet the statistical significance of this correlation was very low. Furthermore it returns an abnormally high EPR Index. This is because of the statistical nature of both tests. The small size of the INTACT_SMALL dataset makes the significance of any results highly dubious. The DIP_Core dataset, our reference dataset, has one of the lowest error rates. It possesses the second highest EPR Index and also a high LS index. The ITO and UETZ datasets (both large scale Y2H experiments) have high error rates, corroborating previous error rate analysis [13].

To obtain a consensus measure of the error rates, we calculated the mean accuracy rank for each dataset, based upon its rank in each measure. Figure 3 shows a graph of the consensus measure for all the datasets with the exception of the BIND, INTACT_SMALL and MIPS_GENETIC datasets. These three datasets are excluded as accurate error rates cannot be calculated for them.

Randomisation

A randomisation test was carried out to demonstrate that the observed negative correlation was not just a simple combination of the distribution of evolutionary rate and number of interactions. The randomisation study was conducted on those datasets that returned a significant negative correlation, to clarify that it was in fact a biological factor that was responsible for the correlation.

Figure 4 shows the results from a randomisation study conducted on the DIP_Core dataset. It can be seen that the correlation of the DIP_Core dataset stands out amongst the correlations of the randomised datasets. This was also the case for all datasets that showed a significant correlation.

Abundance

We consistently found a negative correlation between evolutionary rate and the expression level of a protein (Spearman's ρ ~ -0.4). In other words proteins that were more abundant in the cell tended to evolve at a slower rate (Table 1).

We also found that datasets that displayed a significant correlation between interactions and evolutionary rate, also displayed similarly strong correlations between abundance and numbers of interactions. The correlation between abundance and interactions was positive, i.e. proteins with high expression levels would participate in more interactions (Table 1).

Partial correlations for evolutionary rate and connectivity when protein abundance is controlled for were also calculated (Table 1). Controlling for protein abundance reduced the magnitude of any significant correlation between evolutionary rate and connectivity.

Evolutionary excess retention

Figure 5 shows the relationship between ER and protein connectivity. The proteins have been logarithmically binned by their number of interactions. It can be seen that proteins with a high number of interactions tend to be older (have a higher ER) than proteins with fewer interactions.

Figure 6 shows the relationship between the expression levels of proteins and ER. The relationship is similar to that observed between ER and interactions, however the correlation is not as strong in this case.

Connectivity vs evolutionary rate

In our analysis of several different datasets we found that a correlation between protein interactions and evolutionary rate exists in some datasets and not in others (Figure 2). Where correlations were observed, they were weak, (between Spearman's ρ: -0.1 and -0.25)) but statistically significant (P-value: < 10^-3). The weakness of observed correlations could be attributed to the incomplete nature of the known yeast interactome. It has been estimated that the ~6000 proteins in yeast participate in at most 40,000 interactions [38–40]. GRID, the largest dataset contains 4907 proteins and a total of 17,598 interactions. Considering this also includes a percentage of false positives, it is clear that only a fraction of yeast interactions have been measured. There is also very little overlap between interactions returned from different experimental methods. Results returned by Gavin et al (2002) measure 3957 interactions using the TAP method. Only 63 of these interactions can be found in the Uetz et al's (2000) dataset. Although this could be due to different experimental methods favouring different types of interactions, the lack of agreement denotes the partial nature of the picture to date.

The large fraction of missing data on evolutionary distance is another factor that may explain the weakness of correlations. When searching for orthologs in the M musculus we found BRH orthologs for only half of the proteins in the GRID dataset, (Table 1). This resulted in 50% of the nodes from the interactome missing from our final correlation graph (Figure 1). A highly interacting node which is missing one interaction would simply move its position on the graph, however if the node itself is missing the graph will be missing a point. In this case 50% of the points are missing and it is entirely possible that the 50% that are present may have the wrong number of interactions.

To address the issue of missing orthologs, we searched for orthologous proteins in the more closely related species S paradoxus. This resulted in far more orthologs being found (over 80%), however the strength and general pattern of the correlations remained the same as before (Figure 2). The slight disparity between the cerevisiae-musculus and cerevisiae-paradoxus based correlations is probably due to the relatively small amount of evolutionary change that occurred between the cerevisiae and paradoxus. However the correlation still remains weak even though we located more orthologs. It therefore seems that the primary reason for the weak magnitude of the correlation may be the incomplete nature of the network. Our error rate analysis showed that the accuracy of sets varied (Table 2). When considering the accuracy of the sets it is important that we consider the three error rate indicators collectively. Considering any error rate indicator on its own could be misleading as the experimental methods used could bias a particular error rate indicator. For example the MIPS_Genetic dataset has a high EPR index (74.9%). This indicator on its own suggests that this dataset is highly accurate. However the Reference index is 4.08% and the LS index is 27.52%. The EPR index value is explained by the experimental methods used to obtain the interactions in this dataset. The experimental methods used (synthetic lethality and suppression analysis) check for functional interactions. Functionally related genes tend to be expressed in a similar manner and so a high EPR index would be expected.

In order to translate error rates into a meaningful consensus based representation, we calculated the average rank over the three independent measures to give us an Average Accuracy Rank (Figure 3). Three datasets have been omitted. The MIPS_GENETIC dataset as it is a set of functional interactions rather than physical interactions. The BIND dataset as we could only calculate one of the three error rates for it, and the INTACT_SMALL dataset as it will give uncertain error rates due to its small size. From Table 2 and Figure 3, it is clear that the MIPS_PHYSICAL dataset shows no correlation and its accuracy level is amongst the lowest. In general we find that datasets that demonstrate stronger correlations between connectivity and evolutionary rate are more accurate whereas datasets that show no correlation are found to be less accurate.

The UETZ dataset has relatively high accuracy levels for all three of our error measures. Its consensus accuracy is higher than that shown by GRID and MINT (Figure 3). However it shows no statistical correlation between connectivity and rate of change. The UETZ dataset is obtained via the Y2H method. This has previously been shown to be an inaccurate experimental process [13]. A plausible explanation for the lack of correlation is its lack of representation of highly interacting proteins. The UETZ dataset contains a comparable number of proteins to the GAVIN and HO sets, yet it contains a significantly lower number of interactions (Table 1). The UETZ dataset contains 1438 interactions for 1328 proteins which averages to no more than 1.08 interactions per protein, whereas the GAVIN and HO sets average over 2.3 interactions per protein. It is fair to say that with just 1.08 interactions per protein, proteins with more than one interactions are highly under-represented in the UETZ dataset.

The HO dataset has a low accuracy according to our consensus measure (Figure 3), yet it returns a strong-correlation. It possessed the lowest EPR Index and LS index from amongst all the datasets. Von Mering's analysis estimated the experimental method used in this set to have an accuracy of only 2% [13]. The strength of the correlation between connectivity and evolutionary rate in this dataset could be due to a previously discussed artificatually generated association between connectivity and expression level [9]. Specifically the artificial correlation shown by the HO dataset could arise from the experimental method used to generate the interactions. Ho et al used the HMS_PCI protocol, where the bait proteins are transiently overexpressed. This overexpression may have led to the detection of a large number of false interactions for highly expressed genes. However our partial correlation analysis does not support this conclusion, as we find that if we control for expression level, the correlation observed between connectivity and evolutionary rate still exists.

An analysis of the overlap between the accumulative datasets (sets containing data from many sources) and single experimental method datasets (HO, GAVIN, ITO, UETZ) further corroborates findings from the error rate analysis (Table 3). The DIP_Core dataset, has very little overlap with the inaccurate ITO dataset. This gives further support to our initial assertion that the DIP_Core dataset contains a large fraction of good interactions. Interestingly the three accumulative datasets, BIND, MIPS_Genetic and MIPS_Physical, which showed no correlation between interaction and connectivity, had very little overlap with the GAVIN dataset. The GAVIN dataset is obtained by the affinity purification method and our error rate analysis considers it to be quite accurate. The BIND and MIPS databases are missing a very large fraction of affinity-purification data.

It was also noted that the MIPS_Genetic dataset has very little overlap with any of the single experimental method datasets. To a certain extent this is to be expected as the MIPS_Genetic set contains functional interactions as opposed to physical interactions. The lack of congruence between the MIPS_Genetic set and physical interaction datasets highlights the stark differences between functional interactions and physical interactions.

Abundance vs evolutionary rate

We used the mRNA expression levels in yeast as a measure of abundance. Table 1 shows the correlations observed in all the datasets, between the three factors, evolutionary rate, abundance and interactions. Strong and significant correlations between abundance and evolutionary rate are detected for all the datasets bar the SMALL dataset. A possible explanation for the absence of a correlation in the SMALL dataset is that within this set only 60 proteins had both evolutionary rate and expression information. As has previously been the case, such a small set of nodes may be lacking enough information to display a significant correlation.

Datasets with high accuracy also demonstrated a relationship between abundance and interactions. A positive correlation (where abundant proteins tend to possess more interactions), of a similar magnitude to that observed between interactions and evolutionary rate was seen. This correlation was not observed in sets which were considered to be inaccurate. A simple explanation for this could be that proteins which are broadly present in the cell will have a greater functional role and therefore will participate in many interactions.

When comparing the two correlations, the strength of the abundance vs. evolutionary rate correlation is far stronger than the correlation of interactions vs. evolutionary rate. An explanation for this could be that the interactome is far from complete as discussed earlier. Expression data on the other hand, is far more exhaustive, with expression levels known for 6172 proteins [34]. Expression data is also thought to be of a better quality [15].

Previously it was suggested that affinity-purification methods were biased in that they measured more interactions for highly expressed proteins [9]. This assertion was based on the observation, only in affinity purification sets, of a positive correlation between number of interactions and expression levels, i.e. highly expressed proteins had more interactions. It is a questionable claim, as it can be said that highly expressed proteins are more abundant because of their important functional role, and such a role may require it to interact with many proteins.

Our findings throw further doubt on the claim as we also observed positive correlations between expression levels and number of interactions in accumulative datasets as well as our "golden standard" dataset DIP_Core. Accumulative datasets contain interaction information from different sources, small scale, Y2H and TAP methods. The DIP_Core dataset, our dataset of true interactions, compiled and verified from different sources shows a significant positive correlation between expression and interactions (Spearman's ρ:0.1755).

Nevertheless Bloom et al went on to conclude that because of the stronger positive correlation between expression levels and evolutionary rate, the expression levels were responsible for any correlation between the number of interactions and the evolutionary rate.

To judge what effect, if any, expression levels had on the relationship between connectivity and evolutionary rate we calculated the partial correlations based on the Spearman's rank correlation. When controlling for expression, we found that the strength of the correlation between interactions and evolutionary rate did decrease slightly, in cases where it was observed in the first instance (Table 1). This however does not imply that expression levels are the reason why we observe a correlation between interactions and rate of change. Judging by the strength of the correlation between expression and evolutionary rate, we believe that expression is simply a better predictor of evolutionary rate than connectivity.

Interactions vs age

We found proteins with a high ER value tend to participate in more interactions. Conducting our analysis on the DIP_Core dataset, we found a very strong correlation between ER and interactions (Figure 5). This supports the belief that hub proteins are more likely to be older than non-hub proteins and corroborates previous work [41, 42]. The ER differential can be explained by the theory of preferential attachment [43]. New protein's once having entered the interactome by a growth process, are more likely to form connections with proteins that are already highly interacting. As a result of this process proteins that are present in the interactome for a longer period will accrue more interactions, i.e. hub proteins are older. The scale-free nature of the interaction network could also be explained by such a growth process [44]. We also examined the relationship between ER and expression. This correlation was slightly weaker than the correlation between ER and interactions, yet still significant (Figure 6). The correlation suggests that proteins that tend to be older are more abundant. Collectively this suggests that older proteins are not only highly expressed but also participate in more interactions.

In order to ensure that expression was not causing a bias in the relationship between ER and interactions, we took all the proteins from several expression bins and checked for an association between ER and interactions. This stratification analysis is an effective way of checking if expression has any affect on the relationship between ER and connectivity (Figure 7). For the DIP_Core dataset a correlation, between ER and number of interactions, of strength ρ: 0.98301 was detected, this correlation remained when we examined an expression based bin in which a large number of proteins (> 50) were present (rho : 0.91762). This indicates that the abundance levels of proteins had very little effect on the relationship between ER and interactions.