Prediction of Protein Lysine Acylation by Integrating Primary

Subscriber access provided by UNIV NEW ORLEANS

Article

Prediction of protein lysine acylation by integrating primary sequence information with multiple functional features Yipeng Du, Zichao Zhai, Ying Li, Ming Lu, Tanxi Cai, Bo Zhou, Lei Huang, Taotao Wei, and Tingting Li J. Proteome Res., Just Accepted Manuscript • DOI: 10.1021/acs.jproteome.6b00240 • Publication Date (Web): 24 Oct 2016 Downloaded from http://pubs.acs.org on November 1, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Proteome Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Proteome Research

Prediction of protein lysine acylation by integrating primary sequence information with multiple functional features Yipeng Du1,7, Zichao Zhai2,7, Ying Li2, Ming Lu2, Tanxi Cai3, Bo Zhou1, 4, Lei Huang5, Taotao Wei1,*, and Tingting Li2,6*

1

National Laboratory of Biomacromolecules, Institute of Biophysics, Chinese

Academy of Sciences, Beijing 100101, China 2

Department of Biomedical Informatics, School of Basic Medical Sciences,

Peking University Health Science Center, Beijing 100191, China 3

Laboratory of Protein and Peptide Pharmaceuticals & Laboratory of

Proteomics, Institute of Biophysics, Chinese Academy of Sciences, Beijing 100101, China 4

University of Chinese Academy of Sciences, Beijing 100049, China

5

College of Information Science and Engineering, Ocean University of China,

Qingdao, China 6

Institute of Systems Biomedicine, School of Basic Medical Sciences, Peking

University Health Science Center, Beijing 100191, China 7

These authors are equal contributors.

* Correspondence should be addressed to Tingting Li ([email protected]) and Taotao Wei ([email protected]). Tel: +86-10-8280-1585 Fax: +86-10-8280-1001 1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 49

Authorship Y. D., T. W. and T. L. designed the study, collected the experimental data, and wrote the manuscript. Z. Z., Y. L., M. L., T. C., L. H., and T. L. programed the computational method. B. Z. helped to collect the experimental data. All authors contributed to the data interpretation and approved the final manuscript. T. W. and T. L. are the guarantors of this work and, as such, had full access to all of the data in the study, and take responsibility for the integrity of the data and the accuracy of the data analysis.

2


Page 3 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Abstract: Liquid chromatography-tandem mass spectrometry (LC-MS/MS) based proteomic methods have been widely used to identify lysine acylation proteins. However, these experimental approaches often fail to detect proteins in low abundance or absent in specific biological samples. To circumvent these problems, we developed a computational method to predict lysine acylation, including acetylation, malonylation, succinylation, and glutarylation. The prediction algorithm integrated flanking primary sequence determinants and evolutionary conservation of acylated lysine, as well as multiple protein functional annotation features including gene ontology, conserved domains and protein-protein interactions. The inclusion of functional annotation features increases predictive power over simple sequence considerations for four of the acylation species evaluated. For example, the Matthews correlation coefficient (MCC) for the prediction of malonylation increased from 0.26 to 0.73. Performance of prediction was validated against independent dataset for malonylation. Likewise, when tested with independent datasets, the algorithm displayed improved sensitivity and specificity over existing methods. Experimental validation by western blot experiments and LC-MS/MS detection further attested to the performance of prediction. We then applied our algorithm on to the mouse proteome and reported global-scale prediction of lysine acetylation, malonylation, succinylation and glutarylation which should 3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 49

serve as a valuable resource for future functional studies.

Keywords: acylation; acetylation; malonylation; succinylation; glutarylation; prediction; proteomics; posttranslational modification (PTM).

4


Page 5 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Introduction Protein lysine acylation is an evolutionarily conserved protein posttranslational modification (PTM) which uses Acyl-CoA in the reaction as group donor

1-2

.

Acetylation is a well-known and most studied protein lysine acylation which plays important roles in various physiological and pathological processes

3-6

.

Malonylation, succinylation and glutarylation are recently-discovered types of protein lysine acylation, the biological functions of which remain to be elucidated

7-10

.

Liquid

chromatography-tandem

mass

spectrometry

(LC-MS/MS) based proteomic methods play a key role in the identification of the four types of acylation. In 2006, Kim et al. reported the first proteomic survey of acetylation and identified 388 acetylation sites on 195 proteins from the HeLa cell line and mouse liver tissue 3. This study demonstrated that acetylation

distributes

widely

in

multiple

cellular

compartments

and

participates in diverse cellular pathways. Three years later in 2009, Choudhary et al. enlarged that number to 3,600 acetylation sites on 1,750 proteins from MV4-11, A549 and Jurkat cells, which clearly shows that acetylation is one of the major PTMs in the mammalian cellular system 4. In 2010, Zhao et al. reported the identification of 1,300 acetylation sites on 1,047 proteins from human liver tissue 5. Then, in 2012, Lundby et al. quantified 15,474 acetylation sites on 4,541 proteins from 16 rat tissue samples 6. Analysis of these datasets revealed that acetylation is conserved in the same tissue from different species, but varies from tissue to tissue. For example, the overlap of acetylation 5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 49

proteins between mouse and human liver tissue is over 70%, whereas the overlap between human liver tissue and human leukemia cell line is only 14% 5. In addition to acetylation, thousands of malonylation proteins have recently been identified in three independent proteomic studies

11-13

. Thousands of

succinylation proteins and hundreds of glutarylation proteins have also been identified using a similar proteomic approach

9, 14-16

. The identification of

malonylation, succinylation and glutarylation sites were mainly achieved in mouse liver tissue. Similar to acetylation, the modification patterns of these three types of PTM in other tissues might be different from these in the liver. Additional proteomic studies are anticipated to identify the three types of PTM in various tissues. Besides proteomic method, an alternative way is to predict modifications by the computational method. Computational prediction of acetylation has been widely achieved and proven to be useful. For example, Basu et al. predicted histone acetylation by hierarchical clustering analysis of protein sequences and found four novel histone acetylation sites

17

. Gnad et al. performed

acetylation prediction using the support vector machine (SVM) method and obtained results with high sensitivity and specificity 18. We used Acetylation Set Enrichment Based (ASEB) method to predict lysine acetyl transferase specific acetylation sites and successfully identified p300 or PCAF regulated acetylation proteins

19

. In addition, succinylation had been predicted by the

online tool, - SuccFind

20

. However, most prediction methods are based only 6


Page 7 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


on sequence features and ignore other protein attributes which may inform its modification status. Moreover, some types of PTM do not always show significant sequence features. For example, sequence characterization of malonylation is not as distinct as that of serine phosphorylation. On the other hand, other features, such as protein-protein interaction (PPI) and gene ontology (GO) annotation, are valuable for the characterizations of modified proteins. Our previous work on phosphorylation substrates prediction demonstrated that the prediction performance was significantly improved after integrating sequence features with such functional features 21. In an attempt to establish a relatively complete catalogue of acylation proteins, we developed a computational acyl-lysine prediction method. Compared with previous prediction methods which are based only on sequence information, our method greatly improved the prediction performance by integrating the functional features of proteins. We predicted the mouse proteome, and found more acetylation, malonylation, succinylation and glutarylation proteins, which have never been identified previously, and validated part of them through experimental approaches. More importantly, our prediction tool can be easily expanded to protein propionylation, butyrylation, crotonylation and other modifications, and thus should facilitate the functional studies of various PTMs in this field.

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 49

Experimental procedures Positive datasets of malonylation, succinylation, acetylation and glutarylation The datasets used as positive samples were obtained as follows. Acetylation data were downloaded from PhosphoSitePlus database (version 070515). Datasets of malonylation, succinylation and glutarylation were obtained from literatures

9, 11-12, 14-15

. Considering that most proteomic identification of

acylation was performed in mouse tissue, we only collected data from mouse tissue. Then four sets of 17 amino acid (17-aa) long sequences centered by each type of acylation sites were extracted. Sequences of less than 17-aa were ignored in order to provide a uniform sequence for calculation. In total, we obtained 10,120 acetylation sites centered 17-aa peptide sequences from 3,608 proteins, 1,292 malonylation sites centered sequences from 486 proteins, 2,554 succinylation sites centered sequences from 670 proteins and 596 glutarylation sites centered sequences from 161 proteins. These acylation sequences were used as positive samples (Table S2). The relative abundance of residues flanking acylation sites was calculated according to a previous method 3 based on the positive samples. Generally, for every amino acid, the frequency equals to the number of occurrence of each amino acid divided by the number of all amino acids observed. The relative abundance of the residue flanking acylation sites is represented by the ratio equal to the frequency of each amino acid in the acylation peptides divided by 8


Page 9 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the background frequency. The background frequency for every amino acid based on the mouse proteome (Swiss-Prot database, version 2015-03-23) is as follows: A 0.070

C 0.022 D 0.049 E 0.070 F 0.037 G 0.065 H 0.026 I 0.043

K 0.056

L 0.100

M 0.022 N 0.036 P 0.062 Q 0.048

T 0.053 V 0.061 W 0.012

R 0.056 S 0.084

Y 0.027

Background set and negative datasets The background set was constructed by all lysine centered 17-aa peptides extracted

from

the

mouse

proteome

(Swiss-Prot

database,

version

2015-03-23). To build a negative dataset, the same number of lysine centered 17-aa peptides as the positive set for each type of acylation were randomly selected from the background set. During the random selection, those peptides included in the positive dataset of each type of acylation were excluded.

Functional features Data for functional features were obtained from the following sources: 1) GO database for Biological Process (BP), 2) Cellular Component (CC) and 3) Molecular Function (MF) annotation (version 2015-06-20); 4) STRING database for PPI (version 10); and 5) Pfam database for Protein Functional Domain (release 28.0). Then for each protein, we could extract a huge number 9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 49

of functional annotations as functional features from these five types of annotations. We only used the over- or under-represented functional annotations of each type of acylation proteins in order to reduce the dimensionality of the feature space. Hypergeometric tests were utilized to detect over- or under-represented functional annotations based on each type of acylation proteins, with the mouse proteome in the Swiss-Prot database (version 2015-03-23) as a background. Annotations with Bonferroni corrected p-values of less than 1e-2 were considered as significant and used as functional features. For each of the five types of annotations, if there were too many significant functional features, the top 100 functional features with minimum p-values were used.

Conservation features To obtain the conservation information of lysine sites, we downloaded the multiple sequence alignment data of vertebrates from the eggNOG database (version 4.5)

22

. Then for each lysine site on a mouse protein, we tried to find

out whether it was conserved in human, cow, rat and zebrafish. Usually more than one homolog proteins can be found for the mouse protein in certain species. In these cases, we defined the lysine site to be conserved when it was found on at least one homolog protein.

Converting three kinds of features into numeric vectors 10


Page 11 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


To input flanking primary sequences, evolutionary conservation of acylated lysine, as well as multiple protein functional annotations as feature vectors for SVMs，we should transform them into numeric vectors. To transform primary sequence into numeric vectors, each amino acid was represented by a 20-bit binary vector (each bit is an indicator for one of twenty amino acids). For instance, lysine (K) is represented by a 20-dimensional vector [0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0]T and cysteine (C) by vector [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0]T. Therefore, if the window size of a candidate sequence is 17, the dimension of the sequence feature vector representing it is 20*17. To transform the functional features into numeric vectors, for all the over-represented or under-represented functional features detected by the Hypergeometric test, if the protein has the feature, it will be represented by “1”, otherwise, it will be assigned “0”. Therefore, for a type of acylation proteins, if totally N over/under-represented functional annotations were obtained, the dimension of functional features of all the peptides should be N. To transform the conservation information into numeric vectors, the conserved site was represented by “1”, otherwise by “0”. Then the conservation information of each lysine site is transformed to a 4-dimensional vector which indicates its conservation state in human, cow, rat and zebrafish. Fox example, if the lysine site on one mouse protein is found in all these four species, the conservation features of this site are denoted as “1111”. 11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 49

Finally, for each peptide in the positive and negative datasets of a type of acylation, the input vector for SVM should include 20*17 dimensions for sequence features, N dimensions for functional features and 4 dimensions for conservation features.

Training and testing with Support Vector Machines The support vector machines (SVMs) were utilized to predict each type of the four acylation. The machine learning package scikit-learn (scikit-learn-0.16.1) was used to construct SVM models with the radial basis function (RBF) kernel23. High sequence similarity between the training and testing sequence might cause a bias in performance evaluation. So for the positive dataset of each of the four acylation species, sequences with more than 70% sequence similarity were treated as redundant sequences and these data were deduplicated. For the negative dataset, during the process of random selection from the background set, if the selected sequence was over 70% identical with the previous selected ones, it was removed. To assess the prediction performance, for both positive and negative samples, 4/5 of each dataset were used as the training set and the remaining 1/5 were used as the testing set. Some acylation sites were on the same proteins. If two sites on the same protein are divided into the training and testing sets, there might be an overestimate of the effect of functional features. Therefore, we divided the positive samples at the protein level. The overall process was performed 1,000 12


Page 13 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


times, and the final evaluation was based on the average of these 1,000 performances. To avoid evaluation bias, the set of significant functional features was re-selected only based on samples in the positive training set in each round of test. The schematic view of the prediction strategy can be found in Figure 1. To estimate the effect of each type of features, we constructed different feature groups as follows: 1) sequence only; 2) sequence with lysine conservation; 3) sequence with GO BP; 4) sequence with GO CC; 5) sequence with GO MF; 6) sequence with PPI; 7) sequence with Protein Functional Domain; 8) sequence with lysine conservation and all functional features together.

Performance Metrics Four statistical measures were used to describe the overall performance of the prediction. They were Sn, Sp, MCC, and Area Under the Curve (AUC). The positive set and negative set were represented by P and N, respectively. We abbreviated true positive, true negative, false positive and false negative as TP, TN, FP, FN for convenience. Thus: P = TP + FN N = TN + FP

Sn measures the proportion of actual positives that are correctly identified, and Sp measures the proportion of actual negatives that are correctly identified. 13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 49

MCC evaluates true and false positives and negatives simultaneously, and constitutes a more balanced parameter. The formulas are as follows: TP TP + FN TN Sp = TN + FP

Sn =

MCC =

(TP × TN ) − ( FN × FP ) (TP + FN ) × (TN + FP ) × (TP + FP) × (TN + FN )

Building the SVM classifier for each of the four acylation The SVM prediction model of each of the four acylation was built by the corresponding positive dataset and negative dataset respectively. The samples in the negative dataset of each of the four acylation were randomly selected from the background set. The background set was built by all mouse lysine centered 17-aa peptides. So the size of the background set is very huge. The number of 17-aa sequences in the background set was about half a million. We found that for those acylation with a small positive sample size, the prediction results of two SVM models built with two different sets of negative samples were quite different. Take glutarylation which have 596 positive peptides as an example, with two randomly selected negative datasets, two glutarylation SVM prediction models can be generated. Using these two models to screen glutarylation substrate sites from the background set, the overlap of two putative glutarylation substrate site lists was 82.52% (Table 2). These results indicate that the selection of different negative samples will influence the final prediction results. To make the prediction results stable, we 14


Page 15 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


used the strategy of multiple SVM models. Specifically, for each kind of four acylation, multiple SVM prediction models were built. Each SVM model was built with the same positive set and different negative datasets. The final prediction result was decided by the vote of multiple models. Take the vote number as cutoff, candidates with vote numbers higher than or equal to the cutoff will be predicted as putative substrates. Then for each kind of acylation, using the SVM classifiers with multiple SVM models to screen substrate sites from the background set, if the overlapped putative substrate sites of two classifiers reached more than 90% with the same cutoff, we believed that the prediction results were stable. For SVM classifiers with N models, the vote number (N+1)/2 was used as the cutoff. Since the positive sample sizes were quite different for these four types of acylation, the number of models which can generate stable predictions for four types of acylation were determined, respectively.

Protein purification, in-gel digestion and peptides enrichment Several proteins in the prediction results were randomly selected for experimental validation. The FLAG-tagged proteins were expressed in 293T cells, immunoprecipated with anti-FLAG M2 affinity gel (A2220, Sigma), and eluted by FLAG peptide (F4799, Sigma). Purified proteins were performed by SDS-PAGE and blotted by anti-flag (F1804, Sigma) and anti-malonyllysine, succinyllysine and glutaryllysine (PTM-901, PTM-401 and PTM-1151, PTM 15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 49

Biolabs) antibodies as described previously12. Two immunoprecipated proteins ENO2 and MLYCD were selected for malonylation and succinylation site detection by LC-MS/MS. Briefly, the proteins were excised from the SDS-PAGE gels, minced into small pieces, reduced by 10 mM dithiothreitol (DTT), alkylated by 10 mM iodoacetamide (IAM) (45 min at RT in the dark) and digested with sequencing-grade modified trypsin (Promega) with a substrate ratio of 1:50 (w/w) overnight at 37 °C. Peptides were extracted by 60% acetonitrile and 5% formic acid (FA). Final peptides were dried and re-solubilized in NETN buffer (50 mM Tris-Cl pH 8.0, 0.5% NP-40, 100 mM NaCl, 1 mM EDTA). anti-malonyllysine or anti-succinyllysine antibody conjugated resin (PTM-904, PTM-402) were added and incubated at 4 ℃ for 12 hours with gentle shaking. The beads were washed two times with NETN buffer, once with ETN buffer (50 mM Tris-Cl pH 8.0, 100 mM NaCl, 1 mM EDTA), and once with purified water. The bound peptides were eluted by 0.1% trifluoroacetic acid (TFA) and dried using the SpeedVac. The dried peptides were re-dissolved in 0.1% FA and desalted with C18 ZipTips.

Malonylation and succinylation identification by LC-MS/MS LC-MS/MS

analyses

were performed

on

a LTQ-Orbitrap

XL mass

spectrometer (Thermo, San Jose, CA) equipped with an Eksigent nanoLC system (Eksigent Technologies, LLC, Dublin, CA). The peptide mixtures were loaded onto 20-mm ReproSil-Pur C18-AQ (Dr. Maisch GmbH, Ammerbuch) 16


Page 17 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


trapping columns (packed in-house, i.d., 150 µm; resin, 5 µm) and eluted to 150 mm ReproSil-Pur C18-AQ (Dr. Maisch GmbH, Ammerbuch) analytical columns (packed in-house, i.d., 75 µm; resin, 3 µm) for MS analysis. Elution was achieved with a 0%–35% gradient buffer B (Buffer A, 0.1% formic acid and 5% acetonitrile; Buffer B, 0.1% formic acid and 95% acetonitrile) at a flow rate of about 300 nL/min over 90 min using a data-dependent method. Full-scan MS spectra (m/z 300–1600) were acquired in the orbitrap with a resolution of 60,000 at m/z 400. Tandem-MS spectra were acquired in the ion-trap mass analyzer. Charge state screening and rejection were enabled. Unassigned charge states and charge states 1 were rejected. Dynamic exclusion was set to 90s. All MS/MS spectra were collected using normalized collision energy (setting, 35%), an isolation window of 3 m/z, and one micro-scan. The acquired MS/MS spectra were searched using the Sequest HT search engine in the Proteome Discoverer™ Software (v1.3, Thermo Fisher Scientific, Bremen, Germany) against the mouse Uniprot database (version 2.3, 50807 sequences) with the target-decoy database searching strategy. The following search criteria were employed: three missed cleavages were allowed; mass tolerance for parent ion was set to 20 ppm; mass tolerance for fragment ion was set to 0.6 Da; minimum peptide length was set at 6; Cys (+57.021 Da， carbamidomethylation) was set as a fixed modification; Met (+15.995 Da, oxidation) and Lys (+86.000 Da, malonylation) or Lys (+100.016 Da,

17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 49

succinylation) were considered as variable modifications. The global false discovery rate (FDR) for both peptides and proteins was set to 0.01.

Results Sequence preference of acylation sites and acyl-lysine prediction performance with sequence features Collection of the large set of acylation sites made the analysis of sequence preferences possible for the four types of acylation. The relative abundance of the 16 residues flanking the modified lysine were calculated and shown by heatmaps (Figure 2). For all of the four types of acylation, lysine (K) residues were abundant, while cysteine (C) residues were scarce at almost all flanking positions. This is consistent with previous analysis and probably reflects limitations of the experimental methods instead of biological meaning (see Discussion)3. Another general feature is the low frequencies of histidine (H), proline (P) and serine (S) surrounding the acylation sites (Figure 2A-D). In addition, acetylation peptides have a preference for negatively charged residues (aspartate (D) and glutamate (E)) in the immediate vicinity of the acetylation site (Figure 2A), and a preference for glycine (G) and methionine (M) at the -1 position, and arginine (R) and tyrosine (Y) at the +1 position. For malonylation peptides, there is a tendency of more D, G and isoleucine (I) in the surrounding of the modified sites (Figure 2B). The most significant feature of succinylation peptides is that a group of short chain hydrophobic amino 18


Page 19 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


acids (alanine (A), valine (V), phenylalanine (F), G and I) are enriched at almost all surrounding positions (Figure 2C). Glutarylation peptides are also characterized by a group of enriched hydrophobic amino acids but less significantly than that of succinylation peptides (Figure 2D). Abundant analysis of residues flanking acylation sites shows the sequence features of different types of acylation, which might contribute to acyl-lysine prediction. To test the prediction performance with different lengths of sequences, acylation site-centered peptides with 7, 9, 11, 13, 15, and 17 amino acids were extracted. The prediction performance was estimate by Sn, Sp, and MCC. As summarized in Table S1, the longer sequence would bring more information and the best prediction performance came from peptides with 17 amino acids. Furthermore, the prediction performance of acetylation, malonylation, and succinylation was comparable (Sn, and Sp around 60%, MCC around 0.25, Table S1), whereas the prediction performance of glutarylation was less robust (MCC around 0.15, Table S1). These results reveal poor prediction performance of lysine acylation with only sequence information.

Impact of integrating functional features and lysine conservation on acyl-lysine prediction We believe that if a given protein takes part in similar pathways and interacts with similar protein partners with certain kind of acylation proteins, it could 19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 49

probably be modified by that type of acylation. So we obtained the over- and under-represented functional annotations of each type of acylation proteins based on annotations from the GO, Pfam and STRING databases. The GO database is comprised of three sections, including CC, MF and BP, which describes the cellular location, molecular function and biological process of a protein, respectively. The Pfam database features a protein by manually curated entries of protein families. The STRING database contains information of PPI. The over- and under-represented terms from each database for the four types of acylation were used as functional features (Table S3-6), and the top two significant annotations are shown (Figure 3A-D). For example, the most significant terms of acetylation proteins in PPI are “Acaa2” and “Etfa”, which indicates that a significant number of acetylation proteins interact with these two proteins (Figure 3A). The most significant functional terms in CC are “mito matrix” and “mito inner membrane”, which indicates that most of acetylation proteins were enriched in mitochondrial matrix or inner membrane (Figure 3A). Besides the above functional features, it has been reported that acetylated lysine sites are more conserved than non-acetylated counterparts24. So the conservation of lysine might be an important feature for acyl-lysine prediction. Here we added each type of functional features and the conservation features to the sequence features separately and estimated the prediction ability for each combination. The Receiver Operating Characteristic (ROC) curves 20


Page 21 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


showed that prediction performance can be improved for all these four types of acylation after adding functional features and the conservation features (Figure 4). But the effect of the conservation features was much less than that of functional features. Among all the functional features, the PPI and CC information contribute to the prediction performance significantly, while the Pfam information contribute less than the other features. The best performance was obtained by combining all the functional features and the conservation features together (Table 1 and Figure 4). For instance, with only sequence, the Sn, Sp, MCC and AUC of malonylation prediction are 56%, 70%, 0.26 and 0.67; after combining the functional features and the conservation features, those can improve up to 79%, 93%, 0.73 and 0.93 respectively (Table 1).

Multi-model stabilization of Acyl-lysine prediction and application to the mouse proteome The above analysis indicated that after integrating functional features and conservation features, our method could predict acylated substrates effectively. Then for each kind of four acylation, a SVM classifier with multiple SVM prediction models was built. Each SVM model was built with the same positive dataset and different negative datasets. Here the multi-model strategy was employed to avoid the influence of different negative datasets on final prediction results (Details were described in the methods section). Using the 21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 49

SVM classifiers of each acylation to screen acylation substrate sites from the background set, if the overlapped putative substrate sites of two classifiers reached more than 90%, we believed that the prediction results were stable. We found that for the acetylation which includes more than ten thousands positive samples, the SVM classifiers with one model is enough to generate stable predictions (Table 2). However, for the glutarylation which there are about six hundred positive samples, the SVM classifiers with at least nine models can generate more than 90% overlapped predictions (Table 2). In the final SVM classifiers, we adopted 1, 5, 9, and 9 SVM models for acetylation, malonylation, glutarylation, and succinylation, respectively. As shown in table 2, for each of the four acylation, overlapped proportions of predicted substrate sites of two classifiers are higher than 90%. There results indicated that our multi-model strategy can avoid the influence of different negative datasets and ensure the stabilization of Acyl-lysine prediction. For those three acylation with more than one model. The number of models with positive predictions can be used as cutoff. With different number of positive predictions as cutoff, we obtained ROC curves for different types of acylation (Figure S1). On these ROC curves, we can see that when set the cutoffs as 1, 3, 5, and 5 for acetylation, malonylation, glutarylation and succinylation respectively, the sensitivity and specificity were balanced. So these cutoffs were used to predict putative acylation sites from the whole mouse proteome containing 513,789 lysine sites on 16,696 proteins 22


Page 23 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(Swiss-Prot database, version 2015-03-23). The final prediction results for each type of acylation were provided as supplementary tables (Table 3 and Table S7, S8). In addition, other cutoffs can also be used according to user’s preferences regarding more comprehensive coverage or fewer false positives.

Evaluating performance of the prediction method In a previous study, Colak et al. reported another proteomic study of lysine malonylation which identified 4,042 malonylation sites on 1,426 proteins in mouse liver13. These data were not included in our positive dataset, which was suitable for independently testing of the prediction performance of our method. Thus, we extracted the malonylation sites with sequence length of 17 amino acids and obtained 3,492 malonylation sites on 1,158 proteins. In the above whole mouse proteome prediction, we obtained 35,959 potential malonylation lysine sites on 2,226 proteins (Table S9). The Venn diagram in Figure 5 displays the overlap among our prediction results, the positive dataset, and the dataset from Colak et al. After excluding those malonylation sites used in our positive dataset from the dataset of Colak et al., there were 2,695 malonylation sites left. We found that our method could predict 1,511 of these 2,695 malonylation sites, which Sn is 56%. We also found that there was less overlap of our positive dataset and the dataset from Colak et al., which might explain the low sensitivity calculated by the dataset from Colak et al. This suggests that more malonylation proteins remains to be identified, and that additional 23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 49

malonylation datasets should be included in the positive dataset to enhance the prediction sensitivity. The results indicate that our method can predict malonylation of proteins, even though the prediction sensitivity for the new dataset is lower than that for the positive dataset. To evaluate the prediction performance further, we randomly selected several proteins which were predicted to be modified by our method but have never been identified previously. The proteins were immunoprecipitated from transiently transfected HEK293T cells, and acylation status was detected by each type of pan antibody (Figure 6A-C). The efficiency of immunoprecipitation and the band position of each protein were detected by flag antibody. The results showed that six out of nine proteins can be malonylated, four out of six proteins can be succinylated, and three out of five proteins can be glutarylated (Figure 6A-C, marked by an asterisk). We also detected the malonylation status of proteins which were not predicted to be malonylated. None of the selected proteins were detected to be malonylated (Figure S2). Furthermore, to identify the acylated lysine sites, the immunoprecipitated ENO2

and

MLYCD

proteins

were

digested

and

enriched

with

anti-malonyllysine or anti-succinyllysine antibody conjugated resin respectively. The enriched peptides were detected by LC-MS/MS. Five malonylation lysine sites in ENO2 and six succinylation sites in MLYCD, which have been predicted by our method, were detected (Figure 7A, B). The representative MS/MS spectra of one succinylated peptides was shown in Figure 7C. All the 24


Page 25 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


spectra and sequence for the identified peptides can be found at Figure S3 and Table S10.

Comparison with other existing tools As for the performance of acetylation prediction, we compared our method with SSPKA (Species-SPecific lysine Acetylation prediction)25, which performs better than most of the existing acetylation prediction tools. We downloaded the training dataset of SSPKA including 747 positive sites and 1,478 negative sites. We then tested the SSPKA algorithm by independent positive datasets. The independent positive dataset includes 881 positive lysine sites, which were randomly selected from our positive acetylation dataset after excluding those ones in the SSPKA training set. The independent negative dataset includes 7,450 sites randomly selected from the background set after excluding sites in proteins of our or SSPKA negative dataset. For the 881 positive lysine sites in the independent positive dataset, 456 (Sn=52%) were predicted to be acetylated at the threshold 0.35 by their online tool (Table 4). For the 7450 negative lysine sites in the independent negative dataset, 1,314 (Sp=82%) were predicted to be acetylated at the threshold 0.35. Then, we trained our method by the same positive training samples of SSPKA and negative samples which were randomly selected from the background set. Our negative dataset building strategy was different from that of SSPKA, since the functional features are at protein level, if we adopt the SSPKA negative training 25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 49

dataset, which is comprised of non-acetylated sites on the same proteins of the positive samples, these functional features would not contribute to the prediction. After building of the classifier, we tested the performance of our method on the same independent positive and negative datasets. 565 out of 881 (Sn=64%) positive sites and 299 out of 7450 (Sp=96%) negative sites were predicted to be acetylated (Table 4). As for the performance of succinylation prediction, we compared our method with SuccFind26, which was the only succinylation prediction method available. Similar to the process of the above acetylation comparison, we trained our method with the same positive training dataset of SuccFind, a randomly chose negative training dataset from the background set, and tested the performance with the same positive and negative testing sets randomly selected from the independent dataset. We found that 310 out of 552 (Sn=56%) positive sites and 890 out of 3,235 (Sp=72%) negative sites were predicted to be succinylated by SuccFind, while 493 out of 552 (Sn=89%) positive sites and 291 out of 3,235 (Sp=91%) negative sites were predicted to be succinylated by our method (Table 4). These results indicate that our method represents an improvement over the existing prediction algorithms.

26


Page 27 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Discussion In the present study, we performed an integrative prediction algorithm to identify potential protein lysine acylation events by utilizing sequence information, conservation information and functional features of protein. SVM was used here since it works well when dealing with high-dimensional and low-sample size data. Sequence information has been widely used to predict unknown PTM sites. However, some types of PTM do not show significant sequence features. Thus, solely relying on sequence features is likely to be insufficient for the reliable prediction of these types of PTM. Our method solved this problem by integrating sequence information with enriched functional annotation features, which greatly improved the prediction performance. Our previous work on the prediction of phosphorylation has shown the advantage of integrating functional features21. Here, we developed this method and extended it to the prediction of various types of acylation. Motif analysis of the four types of acylation revealed that lysine residues are abundant, while cysteine residues are scarce at almost all positions surrounding the modified lysine sites. These features might reflect the bias of the identification experiments. Proteomic identification of acylation sites involves carbamidomethylation of cysteine residues, followed by enrichment of acylation peptides by antibody. The carbamidomethylation of cysteines around acylation lysine residues probably disturbed the binding between antibodies and acylation lysine residues, which caused the low identification of peptides 27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 49

with cysteines. For the acylation lysine residues, if the adjacent lysine residues were modified by the same type of acylation, lysine residues would accumulate after extracting sequences centered by the acylation lysine residues. This bias in sequence preference may obstruct the detection of real acylation lysine residues. Functional features contribute to the prediction results at the whole proteome level, which can compensate for the shortcomings of sequence features. However, those protein level functional features do not improve the prediction performance at the site level. As shown in Table S7, for most of the predicted acylation proteins, multiple sites were predicted to be acylated. We calculated the proportion of predicted acylated lysine sites on each protein (the number of acylated lysine sites/the number of all lysine sites on the same protein) and drew the distribution maps for each type of acylation. There are a percentage of proteins predicted to be acylated on all of the lysine sites (Figure S4). Comparing the distributions of proportion of acylated lysine sites on each protein between our predicted list and known acylation site list, we found that, for Succinylation and Glutarylation, the distribution patterns between predicted sites and known sites are similar (Figure S4). While for Acetylation and Malonylation, the differences were even bigger. This differences is probably due to the detection limit of MS. Another possible reason is that these sites might be acylated under other conditions. In our method, the information of sequence features contribute less than other features, which induces the 28


Page 29 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


unsatisfied performance at the site level. We recommend the usage of our method for predicting acylation at the protein level. Interestingly, functional features of the four types of acylation proteins are all involved in metabolism associated terms, such as “mitochondrial”, “ATP binding”, “fatty acid beta-oxidation”, etc. We also analyzed the proteomic data of another PTM, methylation, which contains 201 methylation sites on 144 proteins downloaded from the PhosphoSitePlus database (Table S11). In contrast to acylation, functional features of methylation proteins are involved in certain terms, including “nucleosome”, “calmodulin binding”, and “histone” (Table S11). This suggests that acylation, including acetylation, malonylation, succinylation, and glutarylation mainly functions in metabolic processes. Various types of acyl-CoA, including acetyl-CoA, malonyl-CoA, succinyl-CoA, and glutaryl-CoA are taken as group donors in the reaction of acylation. These acyl-CoAs are important intermediates in metabolic pathways, such as glycolysis, TCA cycle, fatty acid synthesis, etc. Thus, acylation might constitute one type of feed-back mechanism in the regulation of metabolic processes. Our prediction method using functional information is reasonable, because the acylation proteins would probably be metabolism-associated ones. Altogether, our method presented here provides novel and useful insights for biologists to investigate the function of various types of acylation, and should facilitate the elucidation of discovered PTMs in health and disease.

29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 49

Acknowledgements The authors thank Dr. Torsten Juelich (Peking University, China) for the critical reading and editing of the manuscript, and thank Mengmeng Zhang (Institute of Biophysics, CAS, China) for assistance with the LC-MS/MS experiment. This work was supported by the National Natural Science Foundation of China (No. 31400666, 31371337, 31671175), the National Basic Research Program of China (No. 2012CB934003), the National High-Tech Research and Development Program of China (No. 2012AA020401), the Major Equipment Program of China (No. 2011YQ030134) and the National Laboratory of Biomacromolecules.

30


Page 31 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Reference 1.

Wagner, G. R.; Hirschey, M. D., Nonenzymatic protein acylation as a carbon stress regulated by

sirtuin deacylases. Mol Cell 2014, 54, (1), 5-16. 2.

Lin, H.; Su, X.; He, B., Protein lysine acylation and cysteine succination by intermediates of energy

metabolism. ACS Chem Biol 2012, 7, (6), 947-60. 3.

Kim, S. C.; Sprung, R.; Chen, Y.; Xu, Y.; Ball, H.; Pei, J.; Cheng, T.; Kho, Y.; Xiao, H.; Xiao, L.; Grishin,

N. V.; White, M.; Yang, X. J.; Zhao, Y., Substrate and functional diversity of lysine acetylation revealed by a proteomics survey. Mol Cell 2006, 23, (4), 607-18. 4.

Choudhary, C.; Kumar, C.; Gnad, F.; Nielsen, M. L.; Rehman, M.; Walther, T. C.; Olsen, J. V.;

Mann, M., Lysine acetylation targets protein complexes and co-regulates major cellular functions. Science 2009, 325, (5942), 834-40. 5.

Zhao, S.; Xu, W.; Jiang, W.; Yu, W.; Lin, Y.; Zhang, T.; Yao, J.; Zhou, L.; Zeng, Y.; Li, H.; Li, Y.; Shi, J.;

An, W.; Hancock, S. M.; He, F.; Qin, L.; Chin, J.; Yang, P.; Chen, X.; Lei, Q.; Xiong, Y.; Guan, K. L., Regulation of cellular metabolism by protein lysine acetylation. Science 2010, 327, (5968), 1000-4. 6.

Lundby, A.; Lage, K.; Weinert, B. T.; Bekker-Jensen, D. B.; Secher, A.; Skovgaard, T.; Kelstrup, C.

D.; Dmytriyev, A.; Choudhary, C.; Lundby, C.; Olsen, J. V., Proteomic analysis of lysine acetylation sites in rat tissues reveals organ specificity and subcellular patterns. Cell Rep 2012, 2, (2), 419-31. 7.

Peng, C.; Lu, Z.; Xie, Z.; Cheng, Z.; Chen, Y.; Tan, M.; Luo, H.; Zhang, Y.; He, W.; Yang, K.; Zwaans,

B. M.; Tishkoff, D.; Ho, L.; Lombard, D.; He, T. C.; Dai, J.; Verdin, E.; Ye, Y.; Zhao, Y., The first identification of lysine malonylation substrates and its regulatory enzyme. Mol Cell Proteomics 2011, 10, (12), M111 012658. 8.

Zhang, Z.; Tan, M.; Xie, Z.; Dai, L.; Chen, Y.; Zhao, Y., Identification of lysine succinylation as a new

post-translational modification. Nat Chem Biol 2011, 7, (1), 58-63. 9.

Tan, M.; Peng, C.; Anderson, K. A.; Chhoy, P.; Xie, Z.; Dai, L.; Park, J.; Chen, Y.; Huang, H.; Zhang,

Y.; Ro, J.; Wagner, G. R.; Green, M. F.; Madsen, A. S.; Schmiesing, J.; Peterson, B. S.; Xu, G.; Ilkayeva, O. R.; Muehlbauer, M. J.; Braulke, T.; Muhlhausen, C.; Backos, D. S.; Olsen, C. A.; McGuire, P. J.; Pletcher, S. D.; Lombard, D. B.; Hirschey, M. D.; Zhao, Y., Lysine glutarylation is a protein posttranslational modification regulated by SIRT5. Cell Metab 2014, 19, (4), 605-17. 10. Hirschey, M. D.; Zhao, Y., Metabolic regulation by lysine malonylation, succinylation and glutarylation. Mol Cell Proteomics 2015. 11. Nishida, Y.; Rardin, M. J.; Carrico, C.; He, W.; Sahu, A. K.; Gut, P.; Najjar, R.; Fitch, M.; Hellerstein, M.; Gibson, B. W.; Verdin, E., SIRT5 Regulates both Cytosolic and Mitochondrial Protein Malonylation with Glycolysis as a Major Target. Mol Cell 2015, 59, (2), 321-32. 12. Du, Y.; Cai, T.; Li, T.; Xue, P.; Zhou, B.; He, X.; Wei, P.; Liu, P.; Yang, F.; Wei, T., Lysine malonylation is elevated in type 2 diabetic mouse models and enriched in metabolic associated proteins. Mol Cell Proteomics 2015, 14, (1), 227-36. 13. Colak, G.; Pougovkina, O.; Dai, L.; Tan, M.; Te Brinke, H.; Huang, H.; Cheng, Z.; Park, J.; Wan, X.; Liu, X.; Yue, W. W.; Wanders, R. J.; Locasale, J. W.; Lombard, D. B.; de Boer, V. C.; Zhao, Y., Proteomic and Biochemical Studies of Lysine Malonylation Suggest Its Malonic Aciduria-associated Regulatory Role in Mitochondrial Function and Fatty Acid Oxidation. Mol Cell Proteomics 2015, 14, (11), 3056-71. 14. Rardin, M. J.; He, W.; Nishida, Y.; Newman, J. C.; Carrico, C.; Danielson, S. R.; Guo, A.; Gut, P.; Sahu, A. K.; Li, B.; Uppala, R.; Fitch, M.; Riiff, T.; Zhu, L.; Zhou, J.; Mulhern, D.; Stevens, R. D.; Ilkayeva, O. R.; Newgard, C. B.; Jacobson, M. P.; Hellerstein, M.; Goetzman, E. S.; Gibson, B. W.; Verdin, E., 31



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 49

SIRT5 regulates the mitochondrial lysine succinylome and metabolic networks. Cell Metab 2013, 18, (6), 920-33. 15. Park, J.; Chen, Y.; Tishkoff, D. X.; Peng, C.; Tan, M.; Dai, L.; Xie, Z.; Zhang, Y.; Zwaans, B. M.; Skinner, M. E.; Lombard, D. B.; Zhao, Y., SIRT5-mediated lysine desuccinylation impacts diverse metabolic pathways. Mol Cell 2013, 50, (6), 919-30. 16. Colak, G.; Xie, Z.; Zhu, A. Y.; Dai, L.; Lu, Z.; Zhang, Y.; Wan, X.; Chen, Y.; Cha, Y. H.; Lin, H.; Zhao, Y.; Tan, M., Identification of lysine succinylation substrates and the succinylation regulatory enzyme CobB in Escherichia coli. Mol Cell Proteomics 2013, 12, (12), 3509-20. 17. Basu, A.; Rose, K. L.; Zhang, J.; Beavis, R. C.; Ueberheide, B.; Garcia, B. A.; Chait, B.; Zhao, Y.; Hunt, D. F.; Segal, E.; Allis, C. D.; Hake, S. B., Proteome-wide prediction of acetylation substrates. Proc Natl Acad Sci U S A 2009, 106, (33), 13785-90. 18. Gnad, F.; Ren, S.; Choudhary, C.; Cox, J.; Mann, M., Predicting post-translational lysine acetylation using support vector machines. Bioinformatics 2010, 26, (13), 1666-8. 19. Li, T.; Du, Y.; Wang, L.; Huang, L.; Li, W.; Lu, M.; Zhang, X.; Zhu, W. G., Characterization and prediction of lysine (K)-acetyl-transferase specific acetylation sites. Mol Cell Proteomics 2012, 11, (1), M111 011080. 20. Xu, H. D.; Shi, S. P.; Wen, P. P.; Qiu, J. D., SuccFind: a novel succinylation sites online prediction tool via enhanced characteristic strategy. Bioinformatics 2015. 21. Li, T.; Du, P.; Xu, N., Identifying human kinase-specific protein phosphorylation sites by integrating heterogeneous information from various sources. PLoS One 2010, 5, (11), e15411. 22. Huerta-Cepas, J.; Szklarczyk, D.; Forslund, K.; Cook, H.; Heller, D.; Walter, M. C.; Rattei, T.; Mende, D. R.; Sunagawa, S.; Kuhn, M.; Jensen, L. J.; von Mering, C.; Bork, P., eggNOG 4.5: a hierarchical orthology framework with improved functional annotations for eukaryotic, prokaryotic and viral sequences. Nucleic Acids Res 2016, 44, (D1), D286-93. 23. Pedregosa, F.; Varoquaux, G.; Gramfort, A.; Michel, V.; Thirion, B.; Grisel, O.; Blondel, M.; Prettenhofer, P.; Weiss, R.; Dubourg, V.; Vanderplas, J.; Passos, A.; Cournapeau, D.; Brucher, M.; Perrot, M.; Duchesnay, E., Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research 2011, 12, 2825-2830. 24. Weinert, B. T.; Wagner, S. A.; Horn, H.; Henriksen, P.; Liu, W. R.; Olsen, J. V.; Jensen, L. J.; Choudhary, C., Proteome-wide mapping of the Drosophila acetylome demonstrates a high degree of conservation of lysine acetylation. Sci Signal 2011, 4, (183), ra48. 25. Li, Y.; Wang, M.; Wang, H.; Tan, H.; Zhang, Z.; Webb, G. I.; Song, J., Accurate in silico identification of species-specific acetylation sites by integrating protein sequence-derived and functional features. Sci Rep 2014, 4, 5765. 26. Xu, H. D.; Shi, S. P.; Wen, P. P.; Qiu, J. D., SuccFind: a novel succinylation sites online prediction tool via enhanced characteristic strategy. Bioinformatics 2015, 31, (23), 3748-50.

32


Page 33 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Tables Table 1. Prediction performance of integrating various features evaluated by Sensitivity (Sn), Specificity (Sp), Matthews correlation coefficient (MCC) and AUC. seq, sequence features; Conser, conservation features; Pfam, protein domain features; MF, molecular function features; BP, biological process features; CC, cellular component features; PPI, STRING PPI features; all, an integration of all of the above features; Ace, acetylation; Mal, malonylation; Succ, succinylation; Glut, glutarylation. Sn seq 60% +Conser 62% +Pfam 58% +CC 65% +MF 54% +BP 55% +PPI 55% +all 61%

Ace Sp 64% 63% 68% 83% 77% 79% 86% 86%

MCC 0.24 0.26 0.27 0.47 0.32 0.35 0.44 0.48

AUC 0.67 0.67 0.69 0.75 0.69 0.70 0.77 0.78

Sn 56% 57% 60% 68% 64% 61% 78% 79%

Mal Sp 70% 70% 67% 86% 74% 81% 94% 93%

MCC 0.26 0.27 0.27 0.55 0.38 0.43 0.72 0.73

AUC 0.67 0.68 0.71 0.87 0.76 0.80 0.90 0.93

Sn 61% 67% 63% 85% 64% 64% 77% 83%

Succ Sp 65% 65% 67% 91% 79% 90% 92% 92%

MCC 0.26 0.32 0.30 0.77 0.44 0.57 0.70 0.76

AUC 0.66 0.67 0.69 0.93 0.76 0.83 0.89 0.94

Sn 52% 72% 57% 88% 58% 68% 84% 92%

Glut Sp 63% 61% 63% 94% 71% 86% 96% 95%

MCC 0.15 0.33 0.21 0.83 0.30 0.54 0.80 0.87

AUC 0.63 0.66 0.66 0.96 0.67 0.85 0.96 0.97

33



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 49

Table 2. Overlapped proportions of predicted substrate sites of two classifiers with different number of models. During this exploration process, we selected candidates with no less than half models positive prediction as putative substrates uniformly.

Acetylation Malonylation Succinylation

Glutarylation

Model numbers 1 1 5 1 5 9 1 5 9

Overlapped proportion 94.51% 85.62% 92.25% 77.15% 88.34% 91.29% 82.52% 88.52% 91.87%

34


Page 35 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3. Number of predicted acylation lysine sites and proteins result from the whole mouse proteome screen

Proteins

Sites

Acetylation

2,923

59,714

Malonylation

2,226

35,959

Succinylation

3,521

40,235

Glutarylation

1,795

21,574

35



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 49

Table 4. Performance comparison of our algorithm with other published tools of acetylation and succinylation. The performance was evaluated by Sensitivity (Sn) and Specificity (Sp).

Sn SSPKA 52% Acetylation Our method 64% SuccFind 56% Succinylation Our method 89%

Sp 82% 96% 72% 91%

36


Page 37 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure legends Figure 1. Schematic view of the prediction strategy. The

collected

acylation

proteins,

including

acetylation,

malonylation,

succinylation, and glutarylation were used as positive samples. Negative samples were randomly selected from the background set. Both positive and negative samples were divided into a 4/5 training set and a 1/5 testing set. The SVM models were then trained by sequence features, conservation features and different types of functional features. The functional features were extracted based only on the training set. Performances were tested by the 1/5 testing set. “Conser” represents conservation features; “Pfam” represents significant Pfam domain features; “GO MF” represents significant GO molecular function features; “GO BP” represents significant GO biological process features; “GO CC” represents significant GO cellular component features; “PPI” represents significant STRING PPI features; “ALL” represents an integration of all of the above features. Sn: sensitivity; Sp: specificity; MCC: Matthews correlation coefficient; AUC: area under receiver operating characteristic curve.

Figure 2. Heatmap showing the relative abundance of amino acid residues flanking acylation sites.

37



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 49

Seventeen amino acid long sequences centered by each type of acylation site were extracted. The frequency of every types of amino acids was calculated at each position, and normalized by the background amino acid frequency. The relative abundance was represented by the normalized frequency for every amino acid. Hierarchical clustering was based on the relative abundance of each amino acid. The heatmap was drawn according to the relative abundance of every amino acid at each position for acetylation (A), malonylation (B), succinylation (C), and glutarylation (D). The color key from red to green indicates a frequency from high to low.

Figure 3. The most significant over-represented terms of functional features. The first two significant over-represented terms from five types of functional features were shown for acetylation (A), malonylation (B), succinylation (C) and glutarylation (D). Significance was represented by adjusted p-value. See also Supplemental Table 4-7.

Figure 4. ROC curve of integration of sequence features with other features for acyl-lysine prediction. ROC curves for acetylation (A), malonylation (B), succinylation (C), and glutarylation (D) prediction were drawn. Each curve represented one integration of sequence features with other features. FPR: false positive rate; TPR: true positive rate. 38


Page 39 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. Validation of prediction results via the independent malonylation dataset. Overlapped lysine sites among the positive dataset, the dataset in Colak et al., and the prediction result of the mouse proteome.

Figure 6. Experimental validation of the prediction results. Flag-tag conjugated mouse genes were overexpressed in HEK293T cells and purified by the immunoaffinity method.

Status

of malonylation (A),

succinylation (B), and glutarylation (C) were detected by pan-antibodies respectively. Kmal, malonylation. Ksucc, succinylation. Kglut, glutarylation. GCKR: Glucokinase regulatory protein. GYK: Glycerol kinase. MLYCD: Malonyl-CoA decarboxylase. ENO2: Enolase 2. Galk1: Galactose kinase. AHSG: Alpha-2-HS-glycoprotein. BTF3: RNA polymerase B transcription factor 3. SIRT5: Regulatory protein SIR2 homolog 5. CDO1: Cysteine dioxygenase type 1. *: expected bands. Acylated lysine position was marked in red.

Figure 7. LC-MS/MS validation of the prediction results. Flag-tag conjugated mouse genes ENO2 and MLYCD were overexpressed in HEK293T cells and purified by immunoaffinity. After digestion and peptides enrichment, malonylated or succinylated peptides were identified by LC-MS/MS. (A), malonylation lysine sites of ENO2. (B), succinylation lysine 39



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 49

sites of MLYCD. Acylated lysine position was marked red. (C), A representative MS/MS spectra of succinylated peptides.

40


Page 41 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1

41



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 49

Figure 2

42


Page 43 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3

43



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 49

Figure 4

44


Page 45 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5

45



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 49

Figure 6

46


Page 47 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7

47



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 48 of 49

Abstract graphic/For TOC only

48


Page 49 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Supporting Information： Table S1. Prediction preformance of different length of amino acids evaluated by Sensitivity (Sn), Specificity (Sp) and Matthews correlation coefficient (MCC). Table S2. Acylation proteins and sites used as positive dataset for the prediction of acetylation, malonylation, succinylation and glutarylation. Table S3. Significant functional features of acetylation proteins. Table S4. Significant functional features of malonylation proteins. Table S5. Significant functional features of succinylation proteins. Table S6. Significant functional features of glutarylation proteins. Table S7. Prediction results of lysine sites of acetylation, malonylation, succinylation and glutarylation from the mouse proteome. Table S8. Prediction results of proteins of acetylation, malonylation, succinylation and glutarylation from the mouse proteome. Table S9. Overlap calculation among the malonylation prediction results, the positive dataset and an independent malonylation dataset. Table S10. Primary LC-MS/MS data for ENO2 and MLYCD. Table S11. Collected methylation proteins and functional enrichment annotations of methylation proteins. Figure S1. ROC curve of different number of models for each type of acylation. Figure S2. Malonylation status of proteins which are not predicted to be malonylated by our algorithm. Figure S3. MS/MS spectra of malonylated or succinylated peptides. Figure S4. Distribution of proteins with different proportion of acylated lysine sites on proteins for each type of the four acylation.

49


Prediction of Protein Lysine Acylation by Integrating Primary

Recommend Documents