Mapping the stabilome: a novel computational method for classifying metabolic protein stability

Features relevant to protein stability

Protein degradation via the proteasome is mediated through poly-ubiquitination [8]. However, there are a number of other post-translational modifications (PTMs) for which a role in either targeting proteins to or protecting proteins from the degradative machinery has been suggested. These modifications include phosphorylation, prolyl hydroxylation, glycosylation and small ubiquitin-related modifier (SUMO) conjugation [9]. For example, phosphorylation can either promote or inhibit ubiquitination by regulating the E3 ligases responsible for ubiquitination [10]. Glycosylation is used as a form of protein quality control in the endoplasmic reticulum, the folding location of most secreted and integral membrane proteins [11], and this modification can act as a signal for misfolded proteins to be translocated to the cytosol for degradation.

A variety of features have been proposed to influence the targeting of a protein to the ubiquitin proteasome system. The N-end rule is one of the best documented examples of sequence-based degradation signals [12]. The rule originally stated that the in vivo half-life of a protein is associated with the identity of its N-terminal residue, causing high levels of binding selectivity for the E3 ligases that target proteins for ubiquitin-mediated degradation. The rule classified N-terminal residues as stabilizing, or belonging to one of three classes of destabilizing residues: primary, secondary or tertiary. Recently, the rule was extended when it was discovered that N-terminal acetylation of most amino acids acts to create N-degrons [13]. Now it is believed that all but two amino acids (glycine and proline) can act as degradation signals upon acetylation, and are therefore considered “destabilizing”. However, the authors note that many proteins may still have a high level of metabolic stability despite the presence of an N-degron. It is possible that long-lived proteins are protected from degradation through complex formation or folding that makes the N-degron inaccessible.

It has been suggested that structural disorder is correlated with protein stability [14]. A bioinformatic study on yeast data found that structural disorder had a weak, but significant inverse correlation with protein half-life [14]. However, a separate study of protein structural disorder found that highly disordered proteins (where high disorder is defined as a protein with greater that 80% sequence disorder) had far greater metabolic stability than proteins with low structural disorder [15].

A correlation between the frequency of certain amino acids has also been reported[1, 7]. One study demonstrated that the frequencies of tryptophan, cysteine, leucine and threonine were negatively correlated with protein stability, and conversely, that glutamic acid, aspartic acid, lysine and asparagine were positively correlated with protein stability [1]. The exact biological mechanisms underlying the relationship between amino acid composition and stability are unknown, though there are likely to be a variety of factors. For example, the PEST hypothesis claims that short hydrophilic sequence segments enriched in certain residues are correlated with protein instability [16]. What mechanisms are behind this are unknown, and it has also been suggested that PEST regions are actually sequence areas enriched in amino acids that confer phosphorylation modification sites [10].

Previous work on modeling protein stability

There are a number of ways to define metabolic stability. For example, one study presented a probabilistic method for classifying the metabolic stability in vitro of chemical compounds [17]. However, to the best of our knowledge there is currently only one published predictor for intracellular protein stability [7]. The contribution of this published work is two-fold: Huang and colleagues identified optimized sets of features relevant to protein stability, and created a predictor that classifies protein stability based upon the “best” feature vectors. Using the Yen [1] data set, the authors created four classes of protein stability (based on a discrete measurement called the protein stability index, or PSI): short (PSI < 2), medium (2 ≤ PSI ≤ 3), long (3 ≤ PSI ≤ 4) and extra long (PSI ≥ 4). The authors went on to define a list of 376 possible feature components that are believed to contribute to protein stability. These included various biochemical/physiochemical attributes of proteins (such as amino acid composition, hydrophobicity and polarity), protein subcellular locations, KEGG enrichment scores, and the number of complexes a protein is involved in.

To determine what features distinguish between different classes of stability, the authors defined three 2-class problems based on their four stability classes: (1) short and medium vs long and extra long, (2) short vs medium and (3) long vs extra long. For each problem, they ordered the feature vectors according to the maximum relevance and minimum redundancy method, which ranks elements in a vector according to their relevance to a target, and redundancy against each other [18]. To optimise the elements in the feature vectors, they used incremental feature selection (IFS) with nearest neighbor (NN) (using 1-cosine distance as the metric) to classify proteins based on increasing feature vector sizes. For each vector size, they used jack-knife cross-validation to determine the accuracy of NN. For problem (1), they reported an overall accuracy of 72.8% using 62 features; for problem (2), an accuracy of 69.8% with 43 features, and for problem (3), an accuracy of 67.8% with 122 features. On a cautionary note, the automated IFS scheme is likely to invoke a “selection bias”, whereby features are optimal for that particular data set, leading to inflated performance figures [19].

The authors reported that localisation to the cell membrane was an important contributor to predicting protein stability. In our own analysis on the set of proteins used in their study, we found that proteins localised to the membrane or the cell surface were significantly over-represented in the “short” class compared to other proteins (Additional file 1: Data Set 1). Furthermore, the most significant feature for problem (1) according to the optimal feature vector is the frequency of hydrophobic residues. We found that an average of 35% of residues in the “short” class proteins were classified as hydrophobic, compared to only 29% in the “extra long” class. Considering the hydrophobic nature of transmembrane domains, this is likely a reflection of the “short” class proteins localised to the membrane. If the GPSP data is indeed biased towards membrane proteins, then it appears that the NN model shares that bias.

The NN method as employed by Huang and colleagues [7] simply classifies a query protein based on the class of the “nearest” protein. As a result, this method does not give us information on how given features influence stability – whether they correlate with stable or unstable proteins. However, the use of probabilistic machine learning methods such as Bayesian networks can give a transparent and explainable model of protein stability. Though Huang and colleagues have presented the first attempt at computationally modeling protein stability, this work can be improved through the removal of potential experimental bias from the data, and the use of a computational method that can explain its predictions.

Bayesian networks

We chose to use Bayesian networks to model metabolic stability and the relationships between relevant features for several reasons. First, a Bayesian network can represent a large number of different types of features to influence the outcome, recognizing only dependencies that we believe exist. Second, by virtue of its probabilistic nature, uncertain observations can be incorporated, missing data can be managed, and flexible queries to probe and explain predictions can be entertained [20, 21]. Finally, though most observations are strictly “Boolean” (true or false), we are able to integrate “continuous” scores produced by methods particularly suited to challenging though largely independent sub-problems, such as support-vector machines (SVMs) applied to detect protein sequence similarity and position-weight matrices (PWMs) applied to recognize PTM sites.

We discuss the selection of relevant variables and their relationships in Section “Data resources and feature identification”, including the use of latent (un-observed) variables. Parameters are set by using the expectation-maximization (EM) algorithm on a training set [22]. We also separate out data used to parameterize PWMs to score PTM sites and data used to train SVMs equipped with the 1-spectrum kernel to map a protein sequence to its stability class [23].

Identifying stability groups of proteins

where R _i is the proportion of cells in bin i for a given gene. This metric for stability classification was also used by Huang and colleagues [7], as well as other studies that have employed the GPSP data [15, 25]. However, the selection of PSI classification thresholds is non-trivial. Therefore, instead of using the PSI, we first explored how proteins grouped based on the original weight distributions across the 7 bins (using Euclidean distances between data points). Figure 1 shows a cluster dendrogram and a heat-map highlighting the R values that the proteins are enriched in. We noted in particular two distinct groups representing the extremes of protein stability. These two groups seemed a natural choice for investigating what features influence stability, as well as constructing a binary classifier. The bottom group, henceforth labelled “unstable”, contains proteins enriched in the highly unstable R1 and R2 bins. The “stable” group contains proteins that are enriched in the stable R5, R6 and R7 bins. The remainder were classified as “non-assigned”. Approximately 20% of the proteins fell into the unstable class, another 20% into the stable class, and the final 60% fell into the non-assigned class.

Identifying features correlated with stability

We investigated five types of features believed to be related to protein degradation: (1) PTMs, (2) domain and architecture types, (3) N-terminal residues, (4) structural disorder and (5) amino acid composition. For each feature (with the exception of amino acid composition) we used Fisher’s Exact test (corrected for multiple testing) to determine whether there was a significant over- or under-representation in the stable or unstable protein classes, relative all other proteins. The results from the statistical analyses on proteins in the unstable class are summarized in Table 2, and Table 3 contains the results for proteins assigned to the stable class. Phosphorylation and acetylation modifications were both over-represented in stable proteins while being under-represented in unstable proteins: 72% of stable proteins were phosphorylated compared to 30% of unstable proteins, and 36% of stable proteins were acetylated compared to only 4% of unstable proteins. The opposite trend was seen in glycosylation, where only 0.3% of stable proteins were glycosylated, compared to 6% of unstable proteins.

The most significant outcomes of the domain and architecture analysis was the strong over-representation of transmembrane domains and signal peptides in the set of unstable proteins. These were over-represented in unstable proteins while being under-represented in stable proteins: 63% of unstable proteins contained transmembrane domains compared to 4% of stable proteins. These findings were supported by an analysis of GO terms, which also showed that unstable proteins were significantly more likely to contain a transmembrane domain than were proteins from the stable or non-assigned class. Signal peptides were present in 44% of unstable proteins compared to 3% of stable proteins.

The results for the N-terminal residue analysis were highly significant, with destabilizing N-terminal residues (“destabilizing” being defined according to the original N-end rule [12]) being over-represented in unstable proteins while being under-represented in stable proteins. We found that approximately 50% of unstable proteins had a “destabilizing” N-terminal residue compared with a far lower proportion (under 10%) in the remaining proteins. The opposite was seen in the stable class, where under 2% of proteins contained destabilizing N-terminal residues.

We grouped proteins according to their structural disorder using the disEMBL software (http://dis.embl.de/html/download.html) referring to the types used previously by Linding and colleagues [26]. The disorder types are defined as “loops/coils” (secondary structures that can serve as conditions for structural disorder), “hot loops” (a subset of loops/coils with high mobility and higher disorder propensity) and “rem465” (remark465 entries in PDB indicating missing coordinates in X-Ray structure). A total of 15 classifications were used, based on the 3 disorder types as defined by Linding and colleagues [26], and 5 levels of disorder as used in an earlier study [15]. Levels were classed according to the percentage of a sequence containing disorded residues of a given type. The levels were: (1) ≥ 0% and ≤ 20% of sequence containing disordered residues, (2) > 20% and ≤ 40%, (3) > 40% and ≤ 60%, (4) > 60% and ≤ 80 %, (5) > 80% and ≤ 100%.

There were a number of disorder types found to be statistically significant. Tables 2 and 3 show that in both stable and unstable classes there was an over-representation of the very low level disorder “hot loops: 0-20%” class. This disorder type was present in 44% of unstable proteins, and present in about 50% of stable proteins. The low level 20-40% disorder class for all disorder types were over-represented in the unstable protein class. The highly disordered class of 80-100% hot loops and Loops/Coils was over-represented in the stable set of proteins.

A predictive model of protein stability

We are interested in resolving what features have a causal relationship with the classification of metabolic stability. We are also interested in establishing dependencies between features relevant to stability. Both these issues can be addressed using machine learning models. We tested the ability of a strawman SVM trained on sequence data (we refer to this model as “SVM”), a BN integrating the aforementioned features (Tables 2 and 3) excluding sequence features (we refer to this model as “BN”), and a complete model (henceforth “BN+SVM”; Figure 2), to classify proteins into stable and unstable classes. The complete model incorporates the SVM output as an additional continuous variable.

We evaluated model performance through receiver operating characteristic (ROC) analysis and calculating the area under the ROC curve (AUC) [27]. We evaluated the performance of three versions of the model: SVM, BN, and BN+SVM. Models were first evaluated on a data set made up of 743 proteins from the stable protein set and 794 proteins from the unstable protein set using ten-fold cross-validation over five different data set splits. Figure 3 shows the performance of the SVM, BN and BN+SVM models for this data set. All three models have AUCs well above random (0.5), and the BN+SVM model performs best. These results are consistent over multiple cross-validation runs with different data set splits. The mean AUC of the BN+SVM model was 0.85 with a standard deviation of 0.0026. The mean AUC of the BN model was 0.84 with a standard deviation of 0.0012, and the mean AUC of the SVM model was 0.76 with a standard deviation of 0.0043. To measure the performance of the models using a threshold, we took the threshold that occurs at the maximum F-score and calculated sensitivity, specificity and Matthews correlation coefficient (MCC). These results including the AUCs are summarised in Table 4.

We also did a performance comparison with Huang and colleagues’ [7] IFS and NN model (see Methods for an explanation of the comparison) in order to determine which method classifies protein stability most accurately. We took a subset of proteins that overlapped with their “extra long” class and our stable class, as well as a subset of proteins that overlapped with their “short” class and our unstable class. Evaluating the NN model on this subset gave an AUC of 0.899 (Figure 4). We also evaluated the performance of our models on this data, and found that the BN+SVM and BN models both had superior performance to the NN model with AUCs of 0.935 and 0.92, respectively (Figure 4)

We created a trimmed data set with transmembrane/secreted proteins removed to test the models on proteins not affected by N-terminal fluorescent tagging. Re-visiting the statistical analysis on this data set, we found that transmembrane domains, signal peptides and glycosylation were no longer over-represented in unstable proteins. We evaluated the performance of the three models with twenty five-fold cross validation due to a smaller number of observations. Figure 5 shows ROC curves for the 3 models trained and tested on the set that was cleansed from proteins that are potentially subject to experimental bias. All three models perform worse on the smaller data set, though performance is still well above random. The BN+SVM model had an average AUC of 0.75 with a standard deviation of 0.0052. The BN model performed worse with an average AUC of 0.66 and a standard deviation of 0.0046. In contrast to the results obtained using the original data set, the SVM model performed better than BN with an average AUC of 0.73 and a standard deviation of 0.0041. For each model we also calculated the sensitivity, specificity and MCC that occurs at the threshold corresponding to the maximum F-score (Table 4).

Classifying the stability of all human proteins

There are several uses for a global prediction of protein stability values. First, the predictions allow us to ascertain how well the model generalizes. Secondly, we can measure the stability of all proteins through estimating how many proteins are classified as stable or unstable. Finally, we can investigate what features are correlated with stability on a global scale, and compare the findings with experimental data.

We used the BN+SVM model to predict the stability of 30,019 protein isoforms catalogued in HPRD. Two predictions were made, one with the BN+SVM model trained on the full data set (prediction 1, or P1), and one with the BN+SVM model trained on the trimmed data set (P2). Additional file 2: Figure S1 shows a density plot of P1, and Additional file 3: Figure S2 shows a density plot of P2. As the BN+SVM model produces a probability representing the belief that the protein is stable, thresholds are required to classify protein predictions into stable and unstable groups. The thresholds can be modified according to a desired level of true positives / false positives, and are most likely be different between P1 and P2. For our analysis of P1, an unstable protein was defined as having a score less than 0.2, while a stable protein scored greater than 0.75. This resulted in 29% of proteins being assigned to an unstable class, and 26% to a stable class (Table 5). When applying these thresholds to proteins from the GPSP training set, 60% of unstable proteins (as according to the cluster analysis) are correctly classified as unstable, and 42% of stable proteins are correctly classified as stable.

Given these predicted stability groups we again tested for the presence of PTMs and domain/architecture types. A full catalogue of statistically significant features is contained in Additional file 4: Data Set 2 (domains) and Additional file 5: Data Set 3 (PTMs). Apart from the features already noted (Table 3), PTMs including methylation, S-nitrosylation and sumoylation were over-represented in stable proteins, and under-represented in unstable proteins. A number of new domain and architecture types were found as well. For example, nuclear localization signals (NLS), nuclear export signals, (NES) and other nuclear related domains were over-represented in stable proteins and under-represented in unstable proteins. A similar result was found by Yen and colleagues [1], who noted that stable proteins were enriched in nuclear GO terms.

For P2, proteins scoring less than 0.3 were assigned to the unstable class and those scoring greater than 0.7 were assigned to the stable class. 14% of proteins fell into the unstable class, and 13% were assigned to the stable class. We found that signal peptides were under-represented in both stable and unstable proteins. Transmembrane domains were under-represented in stable proteins, and though they were over-represented in unstable proteins, the significance level was greatly reduced to what was seen in P1 where we saw a highly significant over-representation of secreted and transmembrane proteins in the unstable class. It appears that when trained on the trimmed data set, the model is scoring secreted and transmembrane proteins such that they are largely falling into the “non-assigned” classification – an expected outcome given the lack of those proteins in the training data. This is supported by an under-representation of glycosylation modifications in the unstable class. These results indicate that when trained on the trimmed data set, the model is no longer making the counterintuitive decisions that may be a result of N-terminal tagging interference with protein localisation.

One useful feature of stability values is the ability to estimate the effect of global stability on specific groups of proteins. To that end, we investigated the global stability of signaling proteins using a list of approximately 6000 proteins annotated with the biological process GO term “signaling”. We found that signaling proteins were highly over-represented in the P1 unstable class (Fisher’s exact test, P = 3.08e-29). There was also a smaller, though still significant over-representation of signaling proteins in the original unstable class determined by cluster analysis (Fisher’s exact test, P = 0.019). However, the opposite was seen in P2 where signaling proteins were significantly under-represented in the unstable class (Fisher’ exact test, P = 2.54e-14) and significantly over-represented in the stable class (Fisher’s exact test, P = 1.34e-16). We also noted several domain types such as SH2 and SH3 known to be involved in signal transduction that were over-represented in both the P1 and P2 stable classes [28].

In summary, the two predicted sets of stability values have some shared, and some distinct features that govern whether a protein might be classified as stable or unstable. We have made the predicted stability values for P1 and P2 available in Additional file 6: Data Set 4, as well as SVM scores for all proteins (Additional file 7: Data Set 5). The Supplementary Material also contains the feature vectors that were used for training/testing the models and generating the predicted values (Additional file 8: Data Set 6, Additional file 9: Data Set 7 and Additional file 10: Data Set 8). While the stability values for P1 are reliable predictions for intracellular proteins without N-terminal signaling sequences, caution should be exercised when considering secreted and transmembrane proteins.

Data resources and feature identification

To identify discrete stability classes we used the normalised (values adding to one) 7-dimensional GPSP micro array data [1]. We grouped the proteins using the hierarchical clustering method provided in the standard R package. The method calculates Euclidean distances between vectors and groups the vectors with the closest distance in a pairwise manner, moving outwards to cluster pairs, and then groups, in a hierarchical manner (Figure 1). As we were interested in identifying dichotomous binary classes that could be used in training a classifier, we chose the two clusters at the top of the hierarchy that were grouped into stable and unstable proteins respectively.

To understand what specific PTMs and protein domain/architecture types correlate with measured protein stability, we used data from HPRD [30]. We also classified proteins according to the original N-end rule as described by Varshavsky and colleagues [12]. The beginning of the mature peptide is first predicted. To detect the location of signal peptide cleavage sites, we used the online predictor signalP (http://www.cbs.dtu.dk/services/SignalP). The protein set was further processed according to the N-end rule as described previously [2, 31]. For proteins not containing a signal peptide, the initiating methionine residue was removed when the second residue was one of C, G, A, S, T, V or P. After processing, the new N-terminal residue was classified as destabilizing if either R, K, H, F, L, W, I or Y. Otherwise it was classified as stabilizing. We then compared the presence of destabilizing N-terminal residues between the stable and unstable protein groups.

An issue that needs to be considered is the bias potentially caused by the N-terminal tag in the data set generated by Yen and colleagues [1]. In order to create a subset of the data that was free from secreted and transmembrane proteins, to explore the possibility of such bias, we used the online predictors SignalP and TMHMM (http://www.cbs.dtu.dk/services/TMHMM/), respectively [32]. Approximately 30% of proteins were predicted to be either secreted or integral to the membrane. This figure is consistent with those previously reported in whole proteome analysis [33].

Models for classifying protein stability

Comparison with previous work

The performance of the SVM, BN and BN+SVM models was compared with the IFS plus NN method described by [7] and provided by the authors. We used a subset of proteins contained in our “stable” class that overlapped with the “extra long” class used by [7], as well as a subset of our “unstable” class that overlapped with their “short” class. When evaluating the performance of our models on this data subset we trained using the full data set (i.e. the proteins in our stable and unstable classes), but tested only using these overlapping subsets.

To evaluate the performance of the NN predictor, we used their set of optimal feature components and values for the “short/medium vs long/extra long” classifier. To produce a continuous output from NN for analysis by ROC, we used the following function: score=D⁻−D⁺ , where for a given test sample, D⁻is equal to the minimum distance between the test sample and all negative training samples, and D⁺ is equal to the minimum distance between the test sample and all positive training samples. For our analysis, proteins in the “unstable/short” class were considered negative, while proteins in the “stable/extra long” class were considered positive.