This is an open access article published under an ACS AuthorChoice License, which permits copying and redistribution of the article or any adaptations for non-commercial purposes.
Article Cite This: Biochemistry XXXX, XXX, XXX-XXX
pubs.acs.org/biochemistry
Sequence-Based Prediction of Cysteine Reactivity Using Machine Learning Haobo Wang,†,‡ Xuemin Chen,†,‡ Can Li,⊥,∥ Yuan Liu,†,‡,∥ Fan Yang,†,‡ and Chu Wang*,†,‡,∥ †
Synthetic and Functional Biomolecules Center, Beijing National Laboratory for Molecular Sciences, Key Laboratory of Bioorganic Chemistry and Molecular Engineering of Ministry of Education, Peking University, Beijing 100871, China ‡ Department of Chemical Biology, College of Chemistry and Molecular Engineering, Peking University, Beijing 100871, China ⊥ Department of Chemical Engineering, Tsinghua University, Beijing 100084, China ∥ Peking-Tsinghua Center for Life Sciences, Peking University, Beijing 100871, China S Supporting Information *
ABSTRACT: As one of the most intrinsically reactive amino acids, cysteine carries a variety of important biochemical functions, including catalysis and redox regulation. Discovery and characterization of cysteines with heightened reactivity will help annotate protein functions. Chemical proteomic methods have been used to quantitatively profile cysteine reactivity in native proteomes, showing a strong correlation between the chemical reactivity of a cysteine and its functionality; however, the relationship between the cysteine reactivity and its local sequence has not yet been systematically explored. Herein, we report a machine learning method, sbPCR (sequence-based prediction of cysteine reactivity), which combines the basic local alignment search tool, truncated composition of kspaced amino acid pair analysis, and support vector machine to predict cysteines with hyper-reactivity based on only local sequence features. Using a benchmark set compiled from hyper-reactive cysteines in human proteomes, our method can achieve a prediction accuracy of 98%, a precision of 95%, and a recall ratio of 89%. We utilized these governing features of local sequence motifs to expand the prediction to potential hyper-reactive cysteines in other proteomes deposited in the UniProt database. We validated our predictions in Escherichia coli by activity-based protein profiling and discovered a hyper-reactive cysteine from a functionally uncharacterized protein, YecH. Biochemical analysis suggests that the hyper-reactive cysteine might be involved in metal binding. Our computational method provides a large inventory of potential hyper-reactive cysteines in proteomes and is highly complementary to other experimental approaches to guide systematic annotation of protein functions in the postgenome era. “isoTOP-ABPP” has been developed to identify cysteines with intrinsically heightened reactivity16 in native proteomes as well as cysteines that are sensitive to modification by endogenous and exogenous fragment electrophiles.9,17 These experimental studies have revealed a rich resource of functional cysteines that will facilitate the detailed mechanistic study of protein functions and enable inhibitor screening and drug design. Computational and bioinformatic methods have also been developed to predict functional cysteines in proteins, such as those involved in post-translational modifications,18−20 metal binding,21 disulfide bond formation,22,23 and catalytic reactions.24 They largely depend on existing databases for feature extraction, training, and prediction. With regard to the intrinsic
C
ysteine is one of the least abundant residues among the 20 protein-encoding amino acids;1 however, with a pKa of ∼8.5, its side chain thiol group could be easily deprotonated into a thiolate with greatly enhanced nucleophilicity when it is situated in a slightly perturbed local microenvironment under physiological-pH conditions.2 Because of its unique chemical properties, cysteine performs diverse functions3 such as stabilizing protein structures,4 binding metals,5 performing catalysis,6 and undergoing numerous post-translational modifications that can regulate protein functions.7 It is therefore of substantial interest to identify functional cysteines with enhanced chemical reactivity, which will contribute to the mechanistic and functional study of individual enzymes as well as accelerate the discovery of functionally unannotated proteins in light of vast arrays of genomic and proteomic data. Extensive efforts have been devoted to profiling functional cysteines in proteomes,8−11 including those with specific posttranslational modifications, involved in metal binding12 and as part of enzymes’ catalytic sites.13,14 In particular, a quantitative activity-based protein profiling (ABPP)15 method named © XXXX American Chemical Society
Special Issue: Future of Biochemistry Received: September 10, 2017 Revised: October 7, 2017
A
DOI: 10.1021/acs.biochem.7b00897 Biochemistry XXXX, XXX, XXX−XXX
Article
Biochemistry cysteine reactivity (nucleophilicity), PROPKA25 is an empirical prediction method for predicting the pKa of a cysteine based on structural features, including hydrogen bonds, desolvation effects, and charge−charge interactions. A decision tree-based classifier program named “COPA” was developed to locate cysteine sites that are prone to redox-mediated regulation in proteins.26 More recently, a novel computational algorithm “HAL-Cy” together with its web server “Cy-Preds”27 was implemented to distinguish reactive cysteines from poorly reactive ones by combining a structure-based hydrogen network energy potential with a sequence-based weighted scoring matrix derived from the isoTOP-ABPP data. The method performed well in terms of prediction accuracy and precision. However, the requirement of a protein structure to evaluate the hydrogen bonding network around the cysteine of interest precludes its general applicability to all predicted protein products obtained by genome sequencing. In the work presented here, we reported the development of a sequence-based computational algorithm named sbPCR (sequence-based prediction of cysteine reactivity) to identify reactive cysteines from any sequenced genome (Figure 1). A
ments that were performed using proteomes from human cancer cell lines.16 As defined by the authors, we grouped all cysteines with isoTOP-ABPP ratios of 5 as the negative control (“the poorly reactive set”). A window of 41 amino acids was selected to center around any given cysteine for sequence analysis. To increase the number of data points in each set and reduce the risk of overfitting, each peptide was BLASTed (E value of ≤0.1) to retrieve homologous sequences from other species. After removal of duplicated sequences, the final data set contains 565 hyper-reactive cysteines and 3880 control cysteines. Peptide Sequence Encoded by CKSAAP. CKSAAP (composition of k-spaced amino acid pairs) is a widely used encoding scheme in bioinformatics. The frequency of each kspaced (k = 0, 1, 2, 3, ...) amino acid pair in any given sequence was counted. During sbPCR, we used k values of 0, 1, 2, 3, and 4 for spacing and considered all 20 standard amino acids plus a generic placeholder for pair composition. This results in a corresponding vector with 2205 (21 × 21 × 5) units. Each unit (“feature”) represents the frequency of the k-spaced amino acid pair found in the fixed 41-residue peptide. We derived such a feature vector for each cysteine-containing peptide in the data set and fed all these vectors into a support vector machine (SVM) for training and testing. We modified the CKSAAP scheme with a variable sequence window (size from 7 to 41) for feature extraction. The window size that produces the best F1 score (see below) was selected for each feature. F(i) = n+
[xi̅ (+) − xi̅ ]2 + [xi̅ (−) − xi̅ ]2 1 n 1 n ∑k+= 1 [xk(+, i ) − xi̅ (+)]2 + n − 1 ∑k−= 1 [xk(−, i ) − xi̅ (−)]2 −1 −
(−) where F(i) represents the F1 score of feature i, x(+) i̅ , xi̅ , and xi̅ represent the average values for positive samples, negative samples, and all samples, respectively, and n+ and n− are the numbers of positive samples and negative samples, respectively. The F1 score is a measure of discrimination between positive and negative samples with respect to feature i. The larger the F1 score, the more important the feature. Support Vector Machine (SVM). The open-source library LIBSVM28 was employed for SVM analysis in sbPCR. The RBF kernel (radial basis kernel), as recommended by the library by default, was adopted. The database was randomly divided into a training set (75%) and a testing set (25%). A prediction model was built on the basis of the training set, and its performance was evaluated by running on the testing set. An average performance was obtained by repeating the process 100 times. The test precision, accuracy, recall, and F-measure were defined as below, where TP, FP, TN, and FN represent the numbers of true positives, false positives, true negatives, and false negatives, respectively.
Figure 1. Overview of the sbPCR (sequence-based prediction of cysteine reactivity) method.
knowledge database of reactive cysteines in the human proteomes was first derived from the quantitative chemoproteomic profiling data by isoTOP-ABPP and amended with homologous sequences by BLAST. The local sequences with reactive cysteines and background controls are then processed by a “skip-gram”-like algorithm to generate motif features, which can be fed into a support vector machine (SVM) to generate an effective combination of predictors. We evaluated the performance of sbPCR and applied it further to predict reactive cysteines across all kingdoms of life as collected in the UniProt database. We validated our predictions in Escherichia coli by performing additional isoTOP-ABPP experiments and identified a novel hyper-reactive cysteine from a previously uncharacterized protein, YecH. Preliminary biochemical data suggested that the cysteine might be involved in metal binding in the bacteria.
■
precision =
TP TP + FP
accuracy =
TP + TN TP + TN + FP + FN
recall =
METHODS Preparation of the Data Set of Hyper-Reactive Cysteines. We compiled the data set of hyper-reactive cysteines from the previously published isoTOP-ABPP experi-
TP TP + FN
F ‐measure = B
2 × precision × recall precision + recall DOI: 10.1021/acs.biochem.7b00897 Biochemistry XXXX, XXX, XXX−XXX
Article
Biochemistry
Figure 2. Encoding local sequence motifs around hyper-reactive cysteines for machine learning. (a) Visualization of the local sequence motifs of the 74 hyper-reactive cysteines and 612 poorly reactive cysteines (compiled from ref 16) using pLogo (http://plogo.uconn.edu). The over- and underrepresented residues are shown above and below the x axis, respectively. The horizontal red line indicates the threshold of the Bonferroni corrected p value of 0.05. (b) Encoding a peptide sequence by composition of k-spaced amino acid pairs (CKSAAP), e.g., C2C for two cysteines separated by two amino acids and Q3K for a pair of glutamine and lysine separated by three amino acids. (c) Improving the motif identification by CKSAAP with a variable selection window (tCKSAAP). In the illustrated example, a fixed window CKSAAP will identify a C2C motif in both peptides regardless of whether the cysteine of interest is hyper-reactive. Adjusting the size of the selection window will help to exclude the C2C motif that is away from the poorly reactive cysteine and improve the discriminatory power of the motif. (d) Generation of local sequence features of hyperreactive cysteines for machine learning. The local sequences (41 residues long) of hyper-reactive and poorly reactive cysteines are extracted from human isoTOP-ABPP data and BLASTed to obtain homologous sequences. The hyper-reactive (“positive”) and poorly reactive (“negative”) sequences are encoded as multiple features by tCKSAAP, which are fed to a support vector machine (SVM). The SVM will generate a prediction model that can best distinguish hyper-reactive cysteines from poorly reactive ones.
phase resin (5 μm). The sample-loaded tube was attached to a 15 cm laser pulled column (75 μm inner diameter, fused silica) packed with 10 cm of Aqua C18 reverse-phase resin (3 μm). The peptides were then eluted using an optimized gradient from 5 to 100% buffer B in buffer A (buffer A consisted of 100% water and 0.1% formic acid, and buffer B consisted of 20% water, 80% acetonitrile, and 0.1% formic acid). In positive ion mode, full-scan mass spectra were acquired from m/z 350 to 1800 using the Orbitrap mass analyzer and the mass resolution was 70000. MS/MS fragmentation was performed in a data-dependent mode, of which the 20 most intense ions were selected from each full-scan mass spectrum for high-energy collision-induced dissociation (HCD) and MS/ MS analysis. MS/MS spectra were acquired with a resolution setting of 17500 using the Orbitrap analyzer. Other parameters in the centroid format were as follows: isolation window, 2.0 m/z units; default charge, 2+; normalized collision energy, 28%; maximum IT, 50 ms; dynamic exclusion, 20.0 s. Data Processing. LC−MS/MS data are searched against the UniProt E. coli database using the ProLuCID29 methodology with a static modification of cysteine (+57.02146 Da) for iodoacetamide alkylation and variable oxidation of methionine (+15.9949 Da). For isoTOP-ABPP data, variable modifications (+464.28596 and +470.29977 Da for light and heavy TEV tags, respectively) are set on cysteine residues. The searching results are filtered by DTASelect,30 and peptides are also restricted to fully tryptic with a defined peptide false positive rate of 1%. Quantification was performed with CIMAGE16 as previously described.
Feature Ranking and Selection. For each feature, we evaluated its ability to distinguish hyper-reactive cysteines from control cysteines and computed an F1 score according to the aforementioned formula. All 2205 features were ranked in descending order by their F1 scores. We started with the top 2000 features for SVM prediction and repeatedly reduced the number of features by half. For each set of features, we tuned parameters C and γ using the grid method in SVM to obtain the best prediction model. A 5-fold cross-validation was performed to evaluate the performance of the prediction model by Fmeasure. Cysteine Profiling in E. coli Proteomes by isoTOPABPP. E. coli K-12 strain BW25113 (the Keio collection) was grown in LB medium at 37 °C overnight to an optical density at 600 nm (OD600) of 0.8 in an oxygenated shaker. Cells were harvested by centrifugation at 4000g for 10 min, resuspended in phosphate-buffered saline (PBS), and lysed by sonication. The protein concentration of the E. coli lysate was determined by the BCA protein assay (Pierce BCA Protein Assay Kit, Thermo Fisher Scientific). Proteome samples were diluted to concentrations of 2 mg of protein/mL of solution in PBS, and two aliquots (1 mL each) were subjected to the standard isoTOPABPP procedures as described previously.16 Liquid Chromatography−Tandem Mass Spectrometry (LC−MS/MS). LC−MS/MS analysis was performed on a Q Exactive Plus mass spectrometer (Thermo Fisher Scientific) coupled to a DIONEX UltiMate 3000 high-performance liquid chromatography system. Condensed samples were pressureloaded onto a pre-equilibrated fused silica capillary tube (100 μm inner diameter) packed with 2 cm of Aqua C18 reverseC
DOI: 10.1021/acs.biochem.7b00897 Biochemistry XXXX, XXX, XXX−XXX
Article
Biochemistry
Figure 3. Evaluating the performance of tCKSAAP feature selection in sbPCR. (a) Top 10 features selected by F1 score in tCKSAAP. (b) Number of tCKSAAP features used in machine learning and the resulting prediction performance as evaluated by the F-measure. The top 250 features (colored red) are used in the final prediction by sbPCR. (c) Distribution of the lengths of peptides within each group of features selected by tCKSAAP. (d) Comparison of the prediction performance by F-measure between tCKSAAP (dashed line) and CKSAAP (filled triangles) with different fixed peptide length windows. (e) Comparison of the prediction performance between tCKSAAP and CKSAAP by receiver operating characteristic (ROC) curves. The areas under the curve (AUCs) of tCKSAAP-based (black) and CKSAAP-based (red) models are 0.9938 and 0.9777, respectively, both of which are significantly larger than that from a random model (purple line).
Recombinant Expression of YecH. YecH was subcloned into a pET-21a (+) vector by amplification from the corresponding chromosomal DNA using the following primers: forward primer, 5′-GGAATTCCATATGGACTCTATTCACGGTC-3′; reverse primer, 5′-CCGCTCGAGGTGACGGCAAATCTTAC-3′ (with a C-terminal His6 tag) or 5′-CCGCTCGAGCTAGTGACGGCAAATC-3′ (without a His6 tag). The NdeI and XhoI sites were used in cloning, and the subcloned vector was verified by sequencing. The YecH(C42S) mutant was generated by QuikChange site-directed mutagenesis using the following primers and verified by sequencing: forward primer, 5′-TCGTCGGCAGAAGGGATGAC-3′; reverse primer, 5′-CTGCCGACGAGGTGTGAAAAC-3′. Wild-type YecH or the C42S mutant was transformed into BL21(DE3) and grown in 1 L of LB medium containing 1 mM ampicillin while being shaken at 37 °C to reach an OD600 of 0.8. The cells were then induced with 1 mM isopropyl β-D-1thiogalactopyranoside (IPTG) and grown at 18 °C for an additional 14 h. The cell pellet was resuspended in a buffer containing 30 mM Tris-HCl (pH 8.5) and lysed by sonication. YecH (pI = 6.17) was purified by anion-exchange chromatography (Q Sepharose, 17-0510-01) using a gradient from 0 to 500 mM NaCl. The fractions were checked by sodium dodecyl sulfate−polyacrylamide gel electrophoresis (SDS−PAGE) and combined accordingly. In-Gel Fluorescence. BL21 lysates overexpressing YecH (wild type or C42S mutant) were diluted to 0.2 mg/mL in PBS. The samples were labeled with 0, 2, 5, 10, or 20 μM IA-alkyne probe for 1 h at room temperature in the dark. Click chemistry was performed with 120 μM CY5-azide, 1 mM TCEP, 100 μM TBTA ligand, and 1 mM CuSO4. The reaction was allowed to proceed at room temperature for 1 h before being quenched
with 5× gel loading buffer. The samples were separated by 15% SDS−PAGE, and probe labeling was visualized by in-gel fluorescence using a Bio-Rad ChemiDoc MP Imaging System. The samples were also immunoblotted with an anti-His6 antibody to normalize the amount of YecH. Ultraviolet−Visible (UV−vis) Absorption Spectroscopy and Inductively Coupled Plasma-Optical Emission Spectrometry (ICP-OES). UV−vis absorption spectra of wildtype YecH and the C42S mutant (purified after anion-exchange chromatography) were scanned at 25 °C in an Agilent 8453 UV−vis spectrophotometer. The samples were also diluted 50fold using 2% nitric acid and boiled to decompose the organic compounds before ICP-OES analysis on a Prodigy 7 (PerkinElmer) spectrometer at the Analytical Instrument Center of Peking University.
■
RESULTS AND DISCUSSION Sequence Feature Analysis. We collected the data of cysteine reactivity profiling from the isoTOP-ABPP experiments published previously (details in Methods) and visualized specific sequence motifs surrounding hyper-reactive residues using pLogo31 (Figure 2a). The sequence plot was generated using 74 hyper-reactive sites (R < 2) as the foreground and 612 poorly reactive sites (R > 5) as the background, which was designed to eliminate potential sequence bias introduced by a preference for iodoacetamide probe labeling and lysine/ arginine enrichment from trypsin digestion. The logo shows that phenylalanine, tryptophan, and cysteine are overrepresented at positions −5, −1, and 2, respectively. Notably, the CXXC motif has been widely reported not only to be an active site motif of thiol disulfide oxidoreductases but also to play an important role in disulfide isomerases.32,33 D
DOI: 10.1021/acs.biochem.7b00897 Biochemistry XXXX, XXX, XXX−XXX
Article
Biochemistry
Figure 4. Prediction of hyper-reactive cysteines across all genomes by sbPCR. (a) Summary of the prediction results. (b) Analysis of the molecular function of predicted hyper-reactive cysteines by DAVID.38,39 (c) Analysis of predicted hyper-reactive cysteines by a sequence similarity network (SSN).43 The network was visualized with an edge E value cutoff of 10−12. Top-ranking clusters of the network are shown in the top inset. Eight selected clusters that contain known hyper-reactive cysteines are illustrated in the close-up at the bottom. Each cluster is highlighted with a specific color based on the functional classification of these proteins.
Feature Extraction by tCKSAAP. To formally define these sequence-based features for machine learning, we resorted to a widely used bioinformatic encoding method called CKSAAP (composition of k-spaced amino acid pair).34−36 The method represents a sequence given by multiple pairs of amino acids and encodes them in the format of “AkB” where A and B are the residues and k is the space separation between them in sequence. For example, a CXXC motif is represented as “C2C” as a “feature” and an EXXK motif is defined as “E2K” as another “feature” (Figure 2b). In the original CKSAAP
application, the length of the sequence window for feature extraction is fixed, and this will undesirably mix partial longrange interaction motifs with the local sequence motifs around the hyper-reactive cysteines. For example, a C2C feature (as if from a local disulfide) may contaminate the data set of hyperreactive cysteines as noise if it happens to be located next to a poorly reactive cysteine (Figure 2c). To address these limitations, we introduced a modified CKSAAP scheme with a variable sequence window for feature extraction. For each k-spaced amino acid pair (feature), we E
DOI: 10.1021/acs.biochem.7b00897 Biochemistry XXXX, XXX, XXX−XXX
Article
Biochemistry
Figure 5. Validating the prediction of sbPCR in E. coli by isoTOP-ABPP. (a) Workflow of isoTOP-ABPP. (b) Summary of the number of hyperreactive and poorly reactive cysteines predicted by sbPCR and/or validated by isoTOP-ABPP in E. coli proteomes. (c) Raw chromatographic peaks of the six hyper-reactive cysteines predicted by sbPCR and validated by isoTOP-ABPP (R < 2). Red and blue traces are traces of light and heavy probe-labeled peptides with quantified ratios and protein names shown above and below, respectively. (d) Summary of the six hyper-reactive cysteines predicted by sbPCR and validated by isoTOP-ABPP, including their protein names, residue positions, and peptide sequences with the hyper-reactive sites marked by asterisks.
centered around the cysteine of interest and truncated the flanking sequences by n residues (n = 3−20). For each defined window, we evaluated the F1 score (see Methods) of the feature to discriminate the hyper-reactive cysteines from poorly reactive ones in the training set and kept the window length with the best F1 score. We named this new algorithm truncated CKSAAP (tCKSAAP) and reasoned that it would tailor the peptide with an appropriate length for specific feature extraction that can be used for SVM-based machine learning. The workflow is shown in Figure 2d. Development and Evaluation of the sbPCR Algorithm. The SVM is a machine learning algorithm widely used as a binary classifier. The training set composed of input vector x and label y (1 and −1 for positive and negative samples, respectively). It aims to find the optimal hyperplane, which has the largest margin to the closest support vectors, to separate the positive and negative samples. We applied tCKSAAP to generate all possible 2205 (21 × 21 × 5) features and ranked them by F1 scores (Figure 3a and Supplementary Table 1).When the top 2000 features (the input vector) were fed into the SVM to distinguish the hyper-reactive cysteines (positive samples) from the poorly reactive cysteines (negative samples), we obtained a discriminatory power of 95.43% by F-measure (see Methods). To avoid the risk of overfitting with a highdimensional vector in SVM, we performed a series of tests with the number of features cut by half, and the results showed that an F-measure of >90% can be retained with only the top 250 features (Figure 3b).
We also evaluated the top-ranking features in terms of sequence window length for tCKSAAP analysis, and as shown in Figure 3c, the top 250 features are significantly enriched with those with a window length of ≤8 whereas the low-ranking features showed a much flatter distribution. The result suggests that the local sequences around the cysteine of interest possess more discriminating power in predicting whether it has enhanced reactivity. Combining the F-measure parameter in cross validation and the length distribution data, we chose the top 250 features to constitute the final vector. A small sacrifice in accuracy may help us not only reduce the training time but also lower the risk of overfitting. We predicted using the features extracted by the classic CKSAAP with a fixed peptide length, and the F-measure results demonstrated that the new tCKSAAP algorithm outperforms the original method at any tested length (Figure 3d). To further evaluate the predictive power of sbPCR, we plotted the true positive rates of the method versus the false positive rates in a receiver operating characteristic (ROC) curve. The test performances of the two methods were compared: one with the top 250 features extracted with the standard CKSAAP using a fixed window of 17 residues (which gives the best F-measure) and the other using the tCKSAAP as described above. The area under the curve (AUC) can be a good criterion for evaluating the prediction algorithm, and as shown in Figure 3e, both methods have strong predictive power. It is also clear that the modified tCKSAAP method (AUC = 0.9938) works better than the standard CKSAAP method (AUC = 0.9777). Using an independent benchmark test with 5-fold cross validation (see F
DOI: 10.1021/acs.biochem.7b00897 Biochemistry XXXX, XXX, XXX−XXX
Article
Biochemistry
Figure 6. Biochemical characterization of the hyper-reactive cysteine in YecH. (a) Colors of resuspended bacteria in a buffer solution (left) and purified proteins by anion-exchange chromatography (Q column) of WT YecH and the C42S mutant. (b) In-gel fluorescence (top) and anti-His6 Western blots of the wild-type (WT) and mutant (C42S) YecH that are labeled with different concentrations of IA-alkyne probes. (c) SDS−PAGE of purified YecH after anion-exchange chromatography with or without reduction of 2-mercaptoethanol (BME) in the gel loading buffer. Possible formation of YecH oligomers in the WT sample (as indicated by arrows). (d) UV−vis absorption spectra of purified YecH (WT and C42S) after anion-exchange chromatography. The WT sample has two unique peaks at 317 and 427 nm. (e) Contents of iron, calcium, and zinc ions in the YecH samples (WT and C42S) purified after anion-exchange chromatography by ICP-OES. The results show that the WT sample is more enriched with iron and zinc ions than the mutant is.
Methods), the tCKSAAP-based method finally achieved an average prediction accuracy of 98%, a precision of 95%, and a recall ratio of 89%. Prediction of Novel Hyper-Reactive Cysteines in All Proteomes. After evaluating the performance of sbPCR in the benchmark test, we applied the method to predict hyperreactive cysteines in all available sequences recorded in the UniProt database, which contain 553231 proteins from 13367 organisms. Of the 2725789 cysteines in these proteins, we predicted 90532 hyper-reactive sites in 68079 unique sequences, which constitutes 3.32% of the total population (Figure 4a and Supplementary Table 2). Because a strong correlation between the cysteine reactivity and functionality has been demonstrated by isoTOP-ABPP in human proteomes,16 we proceed to examine the functional relevance of our predicted hyper-reactive cysteines. We first employed DAVID37,38 (database for annotation, visualization, and integrated discovery) to analyze their molecular functions. We found that metal ion binding, nucleic acid binding, and aldehyde dehydrogenase (NAD) activity are among the top functionally enriched clusters (Figure 4b). It has been reported that a metal cation (acting as a Lewis acid) can assist in the deprotonation of the sulfhydryl groups in cysteines under physiological conditions,39,40 which is consistent with our prediction. The family of aldehyde dehydrogenases has also been known to have hyper-reactive cysteines in their active sites that are activated by divalent metal ions.41 To classify the sequence and function diversity of the predicted hyper-reactive cysteines, we performed a sequence similarity network (SSN) analysis of peptides harboring the 68079 unique hyper-reactive cysteines from sbPCR prediction (after sequence redundancy is removed) using the EFI-EST web tool (http://efi.igb.illinois.edu/efi-est/).42 In SSN analysis, the most related proteins are grouped together in clusters and the cluster size will allow us to assess the conservation of hyper-
reactive cysteines across all organisms, which may also imply the functional importance of these residues. As shown in Figure 4c, the largest cluster consists of thioredoxins, protein disulfide isomerases, and peptide-N(4)-(N-acetyl-β-glucosaminyl) asparagine amidases, all of which contain thioredoxin-like domains and can serve as good positive controls. Clusters 3, 4, and 13 are zinc finger proteins, RNA methyltransferases, and 7-cyano7-deazaguanine synthases, respectively, that all bind metal ions, including zinc43 and iron.44 Clusters 6 and 12 are γ-glutamyl phosphate reductases and retinal dehydrogenases/betaine aldehyde dehydrogenases, respectively, that are also known to have hyper-reactive cysteines. Clusters 5 and 8 are amino acid tRNA ligases, and though it is unclear how to explain the functions of the hyper-reactive cysteines in these enzymes, they were at least validated by isoTOP-ABPP in human proteomes. The presence of all these positive controls in our SSN analysis confirms that the sbPCR method can predict many known hyper-reactive cysteines. These encouraging results prompt us to further validate the novel reactive sites from our prediction in real biological systems. Quantitative Profiling of Reactive Cysteines in E. coli by isoTOP-ABPP. We chose a prokaryotic organism, E. coli, to validate our predictions independently. Three replicates of isoTOP ABPP experiments (Figure 5a) were performed using lysates of E. coli, and 884 cysteines were quantified with valid ratios in total (Supplementary Table 3), among which 47 cysteine sites have ratios of 5. Within these hyper-reactive and poorly reactive sets, sbPCR predicted 20 hyper-reactive sites, and six of them actually have isoTOP-ABPP ratios of
