Proteome-wide prediction of acetylation substrates

Training Set and Key Assumptions.

We used histones as a training set because of the wealth of information known about their PTM patterns and well-developed purification and analytical detection methods, and focused on the major human core histones bearing a total of 56 lysines (H2A: 13; H2B: 19; H3: 13; H4: 11) (Fig. 1A). To date, MS and antibody data suggest that there are 23 “validated” acetylated lysines and 33 lysines that have not yet been observed as acetylated in human histones based on literature [supporting information (SI) Table S1]. We sought to uncover additional acetylation sites within the “not observed” class of lysines in a systematic, rigorous manner via our computational method. We selected parameters that could influence our ability to predict acetylation sites on histones by making a series of assumptions. First, we focused our attention on short stretches of amino acids N- and C-terminal of all 56 lysines. Because structural studies of published KAT domains coupled with peptide substrates typically do not exceed 14–20 aa in length (10, 11), a sliding window of a maximum number of 12 residues flanking each lysine was chosen (Fig. 1B). Residues most proximal to the lysine were given the highest weight (Fig. 1B), assuming that these residues are most important for enzyme recognition, as several studies have shown (10, 11). Second, we varied standard BLAST sequence alignment parameters, including gap penalty, extension, insertion, and deletion scores (Fig. 1C). For lysines in the extreme N- and C-terminal region, such as H3K4 or H2AK129, we normalized the raw alignment score based on the length of the sequence. Additionally, both orientations of the protein sequence (N-terminal to C-terminal or vice versa) were weighted equally. For sequences with lysines that are located in close proximity to each other, such as H3K36 and H3K37, we restricted our alignment matrix so that these sequences did not receive an alignment score. This restriction prevented our training set to be overrepresented with sequences from overlapping fragments of the same protein. Finally, we compensated for structural accessibility by penalizing buried lysines and improving the score of accessible lysines (12). This, however, did not influence our ability to predict acetylation sites on histones, and therefore was not included in our further computations.

We performed a hierarchical clustering of core histone lysines based on the sequences surrounding each of these given lysines. All 56 histone core sequences were aligned to one another, creating a matrix of pairwise alignment scores, generating a hierarchical tree of histone sequences (Fig. 1D). We next classified each lysine into 1 of 2 categories based on its acetylation status reported in literature: “validated” (23 lysines) or “not observed” (33 lysines) (Table S1). Finally, we visually categorized each of the 56 lysines by color coding our tree based on the acetylation status of each lysine (Fig. 1E).

To assess how robust our clustering was and how well it could actually predict lysine acetylation, we took all 56 lysines and performed a leave-one-out cross-validation (LOV) (13) by iteratively excluding one lysine from our training set. Next, we reconstructed the hierarchical tree with the remaining 55 lysines and incorporated the excluded single lysine observation as test data. For each set and combination of predefined parameters (stated above) and in a single run, we performed a LOV analysis to examine the predictive power on all 56 lysines to discover which set of parameters best optimized classification power. If 2 lysines were in overlapping fragments of the same protein, we excluded both of these lysines from our training set when either lysine was a test case. We took each test lysine (total of 56) and traversed through our training tree to find which subgroup of sequences our target sequence formed the tightest cluster with.

A receiving operating curve (ROC) analysis was performed on our test dataset (Fig. S1), where the statistics measure used was the area under curve (AUC). An AUC of 1 represents a perfect prediction, and an AUC of 0.5 random predictions. Each point on a single curve of the ROC plot was calculated by measuring the false positive versus true positive rate of the performance on all 56 lysines for a given parameter(s) under a cutoff alignment score. If the test lysine clustered within a group of validated acetylated lysines (Fig. 2A, red) above the cutoff score, the lysine was predicted to be acetylated. Conversely, if the test lysine clustered within a group of not-observed lysines (Fig. 2A, green) above the alignment score, the lysine was predicted as not acetylated. The default status of the lysine when it did not fall into the above criteria was not acetylated. The best ROC plot achieved an AUC of 0.80, and the parameters in this case included 6 weighted residues to both the left and right of the tested lysine (Fig. S1). A threshold for prediction was also determined based on this plot. To test the significance of this score, we applied the previous procedure to 1,000 random permutations of the labels of the observed and not-observed lysines. The median AUC in these permutations was 0.64, and the maximum score was 0.79; thus, our AUC was statistically significant (P < 0.001).

Computational Prediction of Novel Human Histone Acetylation Marks and in Vivo Validation by Mass Spectrometry.

After hierarchical clustering of all our lysine-embedded histone sequences, we next sought to predict novel acetylation sites in the human core histones. As our tree illustrates in Fig. 2, not-observed lysines that clustered tightly with validated acetylated lysines (green sequences in gray captions) were potential acetylation targets because of their similar sequence constitution. Based on the threshold, determined by the ROC plot, we selected these as candidate sites. The previous method predicted 7 unique acetylation sites in the human core histones; 4 in H2A (K9, K13, K125, K127), 1 in H2B (K116), 1 in H4 (K44), and 1 in H3 (K37) (Fig. 2). This large number of predictive sites was unexpected because histones have been intensely investigated for PTMs in recent years. To test whether these predicted lysines are acetylated in vivo, we used an MS-based approach to examine histone peptides from human cell lines that were asynchronously growing and treated without any HDAC inhibitors (see SI Text). All peptides containing the predicted lysines were identified, and importantly, 4 of our 7 predicted acetyl-lysines were experimentally validated: H2AK9, H2AK13, H2AK125, and H2AK127 (Figs. S2 and S3). Histones H3 and H2B from sodium butyrate-treated human cells also showed H3K37 and H2BK116 acetylation, but because these marks were observed only under these special conditions (see Discussion), we did not count them as validated.

In summary, we correctly predicted 4 of the 7 acetyl-lysine sites, suggesting that our algorithm is capable of identifying acetylation sites in human histone proteins.

Nonhistone Sequence-Based Dataset Prediction and Validation.

Because our computational analysis revealed a high level of sequence homogeneity among acetylated lysines within histone proteins, leading to the successful prediction of unique modified residues, we next wondered if our approach might also enable us to predict nonhistone acetylation sites.

In our first approach, we included a dataset that contained both nuclear and cytosolic proteins from HeLa cells, which were immunoprecipitated with a pan-acetyl antibody (Fig. S4A and Table S2) and identified by MS (14). The precipitate contained peptides with a total of 1,413 lysines, and 51 previously validated acetylation sites. With PredMod, we were able to predict 34 (67%) of these sites correctly (Fig. S4A) when they were surrounded by 6 residues to the left and right (AUC = 0.75, sensitivity S_n = 0.66, specificity S_p = 0.94) (Fig. 3A, orange curve). In total, 6% (85) of the total number of lysines were predicted that were not validated as acetylated (Fp < 6%). Fp is a maximum false positive rate; a true negative count cannot be accurately determined because many of these lysines could potentially be acetylated, but not detected under the experimental procedures used.

In our second dataset, we compiled a list of 32 proteins containing 1,378 lysines with 73 of these reported in literature to be acetylated in vivo and/or in vitro (Fig. S4B and Table S3). With PredMod, we predicted 39 of 73 (53%) lysine marks accurately with Fp < 6.5% (AUC = 0.74, S_n = 0.58, S_p = 0.93) when these were surrounded by six residues to the left and right (Fig. 3B, orange).

Both test datasets exhibited a decrease in performance when larger numbers of residues N- and C-terminal to the target lysine were used (Fig. 3, blue line), suggesting that KATs may recognize a smaller and defined set of residues. Overall, findings from both approaches revealed that our selected parameters for histones were also valid for the prediction of acetylated nonhistone substrates using an ROC analysis approach.

Analysis of Acetylation Motifs.

We next sought to understand which amino acids play a critical role in acetylation site selection, and asked whether there were preferences for certain amino acids near the target acetylated lysines in our datasets. Notably, when we examined the surrounding residues (six residues to the left and right) of a validated acetylated lysine versus a not-observed one in human histone and nonhistone proteins we discovered an enrichment for small residues (G/A in pink), lysines (K in green), and phosphorylatable residues (S/T in blue) (Fig. 4). To test whether the observed enrichment of G, K, S was statistically significant, we determined the frequency of these residues flanking a lysine in the entire human proteome. We noticed that on average, these residues were of significantly higher frequency in our datasets than in the human proteome. We used the hypergeometric test to measure the statistical relevance of this observation (Table S4). Our findings show that the most significant P values were found in the category of small residues (P < 0.01 in multiple flanking positions; Fig. 4, tick marks), suggesting that small amino acids, perhaps due to their sterically undemanding side chains, could accommodate the flexibility of the substrate, thus allowing protein docking and catalysis. This observation was in agreement with a previous study (9), which revealed that glycine preceding lysine was common among acetylated lysines. In conclusion, we were able to identify a significant enrichment of mainly small amino acids and lysines surrounding validated acetylated lysines in comparison with not-observed ones, suggesting that KAT enzymes have a general need for specific residues for recognition and/or activity. These observations are in agreement with studies of several KATs with test substrates (10, 11).

S. cerevisiae Proteome-Wide Prediction and in Vivo Validation.

The previous predictions were performed with human proteins, and we therefore wondered whether our algorithm would also be able to predict acetylation sites in proteins from other organisms. Because histone acetylation has been studied extensively in budding yeast, we assessed the performance of our model on a proteome-wide dataset that included acetylated peptides in S. cerevisiae (15) (see SI Text). In addition, we experimentally validated our predicted acetylation sites in candidate yeast nonhistone proteins in vivo.

In our first approach, we examined in vitro a proteome-wide dataset of acetylated peptides of S. cerevisiae that contained 356 peptides, including acetylated histone peptides (see SI Text). This dataset allowed us to approximate the number of yeast acetylation events on a global level (0.6%; see SI Text), and the substrates themselves allowed us to further validate our prediction algorithm. We filtered these protein-derived peptides according to their cellular compartment (nuclear vs. cytoplasmic) (16), and correctly predicted 43% of acetylation events on nuclear proteins (79 lysines total; AUC = 0.71, S_n = 0.41, S_p = 0.92, Fp < 4%) and 30% on the cytoplasmic proteins (248 lysines total; AUC = 0.70, S_n = 0.31, S_p = 0.90, Fp < 5%). We also noted that nuclear yeast proteins showed a similar enrichment for small residues surrounding the target lysine, as found in the human substrates (Fig. S5).

In our second approach, we validated our predictions on 3 yeast candidate proteins that had previously not been published to contain acetylated sites: Spt6 (17), Sir3 (18), and Eaf7 (19). We expressed and purified our tap-tagged candidate proteins in S. cerevisiae (Fig. 5A) and subsequently subjected them to MS. With PredMod, we predicted 15 sites to be acetlyated of 416 total lysines in our 3 candidate proteins combined. Four of these, within our top 6 ranked predicted sites (Fig. 5B), were validated as acetylated by MS and therefore predicted correctly (Fig. S6). The total number of acetylated lysines in the yeast proteome is ≈0.6%; therefore, our in vivo hit rate of ≈25% is of reasonable accuracy.

These findings show the power of PredMod for identifying bona fide acetylation sites in nonhistone proteins, and further display the strength of using histone sequences as a useful guide for nonhistone acetylation prediction.

Cell Lines.

Mammalian cell lines were grown in Iscove's DMEM supplemented with 10% FCS and penicillin/streptomycin at 37 °C and 5% CO₂.

Histone Isolation.

Nuclei were isolated and histones acid-extracted from asynchronously growing, untreated cells as previously described (23). See SI Text for further details.

MS Analysis of Histones.

Experimental details are described in SI Text.

MS Analysis of Yeast Nonhistone Acetylation Sites.

Tagged cells of our nonhistone proteins were lysed under cryogenic conditions. Tandem TAP-tag purification was performed on candidate yeast proteins as described (24), and eluates run on SDS-PAGE gels and stained with Coomassie. Protein bands were in-gel digested with trypsin or chymotrypsin, and peptides extracted. Details of these methods are provided in SI Text.

Datasets.

Training set: 56 human and S. cerevisiae core histone lysine sequences were collected from the Swiss-Prot database (http://ca.expasy.org/sprot/). Test set: source of nuclear protein and pan-acetyl antibody datasets are described in Results. For information on the budding yeast proteome-wide dataset, see SI Text.

Hierarchical Clustering Analysis.

We performed hierarchical clustering on the sequences surrounding each of the 56 histone lysines. All 56 sequences were aligned to one another, creating a matrix of pairwise alignment scores; our metric was based on these pairwise scores. Sequence alignment scores were computed by performing BLAST local alignments using the NCBI BLAST 2.0 server. A standard BLOSUM62 evolutionary substitution matrix was applied (25).

Statistical Analysis.

ROC calculations are described in the main text. Hypergeometric probability calculation: Pr = (_K^m)(_{n − K}^{N − m})(_Nⁿ) (N, all lysines in human proteome; K, number of times the particular residue is seen flanking in each position in human proteome; n, total number of lysines in each independent validation dataset; m, number of times the particular residue is seen flanking in each position in validation dataset). Sensitivity (S_n) was calculated as the total number of correctly identified acetylation sites from the positive dataset divided by the total positive dataset. Specificity (S_p) was calculated as the total number of negative sites that were not predicted to be acetylated divided by the total negative dataset size. For additional information, please see SI Text.

Sequence Logos.

Sequence logos for displaying the flanking residue distribution of all lysines in our training and test datasets were created according to ref. 26.

Software URL.

Our acetylation prediction software, PredMod, can be found at www.cs.cornell.edu/w8/∼amrita/predmod.html (see SI Text).