Prediction and analysis of redox-sensitive cysteines using machine learning and statistical methods

Statistics

We studied the differences between composition of residues in the 3D neighborhood of Cys+ and Cys−. We considered the 10 spatially nearest residues. Significance values are indicated as follows: *p < 0.05; **p < 0.01; ***p < 0.001.

Figure 1 shows the deviation of the composition of residues in the spatial neighborhood of Cys+ from the composition of residues in the spatial neighborhood of Cys− of protein set 1. The frequencies of lysine, glycine, proline, serine, and threonine are significantly elevated. It has been suggested (Tanner et al. 2011) that serine and threonine may be involved in the formation of sulfenyl esters together with sulfenic acids formed by the oxidation of cysteines, while lysines and histidines may be involved in the formation of sulfenamides, such as in the case of Protein Tyrosine Phosphatase 1B (PTP1B) (Salmeen et al. 2003; Sarma and Mugesh 2007), which may be a possible explanation for the relative abundance of these amino acids. The frequency of leucine is significantly reduced.

Figure 1:

Relative frequencies of residues in the 3D neighborhood of Cys+. Red bars indicate significantly elevated frequencies and green bars significantly reduced frequencies. The error bars are 3 × σ with σ being the Poisson standard deviation for the total number of counts. The frequencies of the amino acids glycine (G), lysine (K), proline (P), serine (S), and threonine (T) are significantly elevated. The relative frequency of leucine (L) is significantly reduced.

When grouped for their features, significant differences could be found between the two groups of cysteines for all features, see Supplementary Figure S2. We found the strongest difference for positive charge as well as the presence of aliphatic side chains. While there is much overlap between the data for modifiable and unmodifiable cysteines, these significant differences may be one important piece of the puzzle when trying to predict modifications.

We created a sequence logo (Crooks et al. 2004) of the different frequencies of the closest amino acids in the sequence neighborhood of the cysteines in the two groups of protein set 2, as seen in Figure 2. The hydrophobic leucine (L) and cysteine (C) are reduced in the sequence logo around a Cys+, especially and significantly at position −1 and positions −3, 3 and 6, respectively.

Figure 2:

Sequence logo (Crooks et al. 2004) of differences between the residues in the neighborhood of Cys+ and Cys−. Sequence position in relation to cysteine are shown on the x-axis, percentage difference of Cys+ in relation to Cys− on the y-axis. Enriched residues around Cys+ are shown at the top, depleted residues at the bottom. Size of symbols is proportional to the difference between the two samples. Only differences with a p-value <0.05 according to the t-test are shown.

Positively charged amino acids like lysine (K) and arginine (R) are unaffected in the direct sequence neighborhood of Cys+, but are enriched in the further sequence neighborhood, i.e., for more than four positions away from the cysteine, confirming the results of Chen et al. (2015), who also found an abundance of positively charged residues around S-nitrosylation sites and a reduced occurrence of C. Note that positively charged amino acids are dominant also in the 3D neighborhood, see Figure 1 and Supplementary Figure S2. Aromatic amino acids like phenylalanine (F) and tryptophan (W) were also reduced in the sequence neighborhood of Cys+, see Supplementary Figure S1.

We were unable to confirm the presence of an acid-base motif which was found in some previous studies on smaller datasets (Greco et al. 2006; Hess et al. 2005). Marino and Gladyshev (2010) proposed the presence of a modified acid-base motif, consisting of a positively charged residue in close proximity to the cysteine and a negatively charged amino acid up to eight Å away.

We detected a higher abundance of phosphorylated residues in the sequence neighborhood of Cys+ than Cys− in protein set two according to Uniprot (Consortium 2019) annotations. On average, each Cys+ had around 0.097 phosphorylated serines, threonines and tyrosines in its neighborhood of 10 residues upstream and downstream, while we observed only 0.055 phosphorylated residues around Cys−, making phosphorylation about 1.76 times more common around Cys+ (data not shown). We found similar results for other post-translational modifications, like acetylation. We discovered 0.063 acetylated residues in the neighborhood of the average Cys+, with only 0.045 modified lysines present around Cys−, a factor of 1.48 (data not shown). Ubiquitination also appeared more commonly around Cys+, where we located an average of 0.016 ubiquitination sites around modified cysteines, while we only detected an average of 0.002 ubiquitination sites around unmodified cysteines. These results were statistically significant even when corrected for the different abundances of serines and lysines around Cys+ in the case of phosphorylation and ubiquitination with p-value <0.001 according to the Mann-Whitney U test. Results for acetylation were not statistically significant.

We explored the accessible surface area as predicted by the algorithm DSSP (Define Secondary Structure of Proteins) (Kabsch and Sander 1983) as well as HSE in the 3D neighborhood of cysteines in protein set 1. Cys+ showed a higher accessible surface area than Cys− with distributions that differed with high statistical significance (p-value <4 × 10⁻⁶), despite much overlap in the range of values, see Supplementary Figure S3. This appears reasonable, as an exposed cysteine should be assumed to be more easily accessible to ROS. HSE showed similar results, see Supplementary Figure S4 and Supplementary Figure S5. We also tested the accessibility and HSE of residues in the 3D neighborhood of Cys+. Again, both DSSP and HSE showed that residues around Cys+ were significantly more accessible.

We examined the relative frequencies of secondary structure elements (SSEs) in the neighborhood of Cys+ in protein set 1. In the 3D neighborhood of Cys+, we found the frequencies of β-strands as well as of bends and unstructured loop regions to be significantly elevated, see Supplementary Figure S6. Statistical significance could not be found for the other SSEs, often due to their low frequency of occurrence in general. This differs from the results of Marino and Gladyshev (2009), who found a marked preference for both α-helical and loop geometries around thiol oxidoreductases, testing a more limited dataset of 75 structures. We found that Cys+ themselves had a much higher chance than Cys− to be present in loop structures, and found a higher incidence of β-strands upstream and α-helices downstream from Cys+ than the reverse, while the ratios for Cys− were more balanced, confirming the findings of Fomenko et al. (2007).

We created a heatmap of the Pearson correlation coefficient between the features in protein set 1, see Supplementary Figure S7. The heatmap showed a high negative correlation between accessibility and HSE, as both features display related physical properties, whereby HSE shows a stronger correlation to cysteine redox-sensitivity. Features of residues closer to the investigated cysteine tended to have a stronger correlation to its redox-sensitivity than those of more distant residues. The most highly correlated features were SSE and HSE of the cysteine.

Machine learning

We applied all three machine learning algorithms on the datasets to be able to infer redox-sensitive cysteine sites in proteins. We used both a dataset consisting of proteins from mammals, plants and fungi as well as a dataset consisting of proteins strictly from mammals. We used 5-fold cross-validation to guard against overfitting. We also implemented feature selection to remove superfluous data.

The best Area Under the Curve (AUC) value for the Receiver Operating Characteristic (ROC) curve was 0.72 for both the imputed mammal and the full imputed datasets using the ET algorithm, see Figure 3, Table 1 and Table 2. The ROC curves depict average values over the different cross-validation folds.

Figure 3:

ROC curves of the ET (blue), RF (green) and SVM (orange) algorithms for the full structural (left), sequence (middle) and imputed sequence (right) datasets. Red line shows a completely random prediction. The ROC curve shows average values over the different cross-validation folds.

To compare the performance of the algorithms with the statistical results, we ran training and testing with cross-validation using only a subset of our features in the same configuration, see Supplementary Table S1. In general, we achieved the highest AUC value when only using the residue identity in the neighborhood of cysteines. We also received favorable results by using only HSE of close residues, which was superior to using only accessibility. This is in agreement with the results of the feature heatmap, which indicated a higher correlation of cysteine redox-sensitivity with HSE. The AUC value using secondary structure was on a similar level as for accessibility.

We evaluated the variance of our features using the ANOVA F-test to compare the feature importances based on the not imputed sequence dataset for all species, see Supplementary Table S2. A higher score signifies that the feature is more dependent in the target variable, i.e. the redox activity of the cysteine. We find that in general, features concerning the cysteine itself, like its secondary structure or exposure, show a high impact on the modifiability of the cysteine. Residue identity also showed high importance, followed by HSE and accessibility.

We trained the ET algorithm, which showed the highest AUC value during model training, on the full dataset 2, using imputation to make novel predictions for cysteine sensitivity in proteins. These proteins were also collected from the RedoxDB and are suspected to contain redox-sensitive cysteines, but their sequence positions are not known. The list included 121 proteins with 575 cysteines. We ran the predictions 10 times for more robust results. Out of the 121 proteins, 44 included cysteines with a prediction score larger than 0.45, which we applied as a threshold for a cysteine likely being a Cys+. We chose this value, instead of the usual value of 0.50, since the algorithm tends to strongly underestimate the number of Cys+ due to the large number of Cys− in our dataset. 11% of the cysteines were predicted as Cys+. For the full list of positively predicted cysteines, see Table 3 and Table 4. We compared the predictions with known data of post-translational modifications from the UniProt (Consortium 2019) database. We found that two of our predictions, CYS37 and CYS40 in Protein Data Bank (PDB) file 1X5D (available at https://www.rcsb.org/structure/1x5d) of protein disulfide-isomerase A6, are known to form a redox-active disulfide bond according to UniProt. None of the other examined cysteines were known locations of redox activity according to UniProt. These findings reflect positively on the reliability of our methods. We predicted redox activity in disulfide bridges not known to be redox-active in two more proteins, CYS 120 and CYS 127 in PDB file 1B56 (Hohoff et al. 1999) for fatty acid-binding protein 5, where both cysteines were predicted as CYS+, as well as CYS105 and CYS137 in PDB files 3O0W (Kozlov et al. 2010) and 3POW (Chouquet et al. 2011) for calreticulin, where only CYS137 was predicted to be CYS+. These results do not seem unexpected and may even be a positive outcome, as disulfide bridges are one of the redox modifications we aim to predict, and not all disulfide bridges that are known to be redox-active according to redoxDB are marked as such in other databases like UniProt. Several other proteins contained disulfide bridges, which were not predicted as redox-active, like prothrombin (3HK3 (Gandhi et al. 2009)), the cytochrome b6-f complex iron-sulfur subunit (1RFS (Carrell et al. 1997)) and the non-specific lipid transfer protein (1BWO (Charvolin et al. 1999)).

Use case

After using different algorithms to create models for the prediction of redox-sensitive cysteines, we applied our models to generate predictions for the NDUFS1, MT-ND3, NDUFA2 and NDUFC2 subunits of mammalian respiratory complex I, based on the structural data from PDB entries 6G2J (Agip et al. 2018), 6G72 (Agip et al. 2018), 5LC5 (Zhu et al. 2016), 5LNK (Fiedorczuk et al. 2016) and 5XTD (Guo et al. 2017).

Mitochondrial complex I is with its mass of around 1 MDa the largest protein complex in the respiratory chain. It is responsible for the oxidization of NADH produced mainly by the Krebs cycle, transferring electrons to the ubiquinone pool in the process (Wirth et al. 2016). The respiratory chain continuously reduces O₂ into H₂O, whereby a small amount of ROS in the form of O 2 − is generated. It has also been proposed that ROS formed at complex III could have a direct feedback on complex I cysteines (Larosa and Remacle 2018).

Which cysteines are likely targets of such activity could be predicted by our models. We applied the models trained with the random forest (RF), support vector machine (SVM) and extra trees (ET) algorithms to the complex, using the parameters and training data detailed in the previous sections. Highest AUC was received for the mammalian structural and imputed datasets. For illustration, we used the PDB entry with the accession number 5XTD of mitochondrial complex I (Guo et al. 2017), as it is the most complete structure of human origin.

We compared the predictions of all algorithms based on the mammalian structural and imputed datasets to the experimental results from Chouchani et al. (2013). 25% of the cysteines of the relevant subunits were predicted correctly by all algorithms in both datasets, see Supplementary Table S5. Additionally, 30% were predicted correctly by the majority of algorithms. Overall, 68% of the individual predictions of the algorithms for the investigated cysteines agreed with experimental data. For the 3D structure of complex I with cysteines experimentally validated to be modified, see Figure 4. For a closer look at the structural features of CYS367, CYS554 and CYS564 in subunit NDUFS1, see Supplementary Figure S8. Supplementary Figure S9 shows all investigated cysteines as well as their structural features, with features that are more commonly observed in Cys+ in green, in red features observed in Cys−, and in black those observed in both groups at about the same frequency. The machine learning models tend to show positive predictions for cysteines with a high number of features that are statistically associated with redox modification.

Figure 4:

3D structure of human complex I using PDB entry 5XTD (Guo et al. 2017). The subunits NDUFS1, MT-ND3, NDUFA2 and NDUFC2 are colored yellow, cyan, blue and purple, respectively. Green dots indicate the positions of experimentally verified (Chouchani et al. 2013) Cys+. Red dots indicate the positions of Cys−.

Data sets

We generated two data sets, the first one contains structurally resolved proteins, whereas the second one was extended by proteins with unknown structure. For the protein set 1, we collected 439 redox-active proteins with 644 Cys+ and 1692 Cys− from the latest update of RedoxDB (Sun et al. 2012) up to May 2020. We selected only proteins with known structure stored in the PDB (Berman et al. 2000) and at least one Cys+. To assign SSEs and accessibility to amino acids, we applied the algorithm DSSP (Kabsch and Sander 1983). Protein set one included 369 mammal proteins, 25 plant proteins and 45 fungal proteins.

For the second data set, we considered 1097 additional redox-active proteins listed in the RedoxDB for which no structures were available. For these proteins, we imputed structural properties of amino acid residuess in proteins, such as, e.g., secondary structure and accessibility, using the values from protein set 1. We applied the function IterativeImputer with Bayesian Ridge estimator (Pedregosa et al. 2011). Combined with model training, we performed imputation separately within each cross-validation fold. By imputation, we could extend the data set by 834 mammal proteins, 127 plant proteins, and 136 fungal proteins, such that it consisted of 1536 redox-active proteins in total. The proteins contained a total number of 2250 Cys+ and 13,373 Cys−. During model training, we removed 55% of randomly chosen Cys− from the dataset to reduce bias towards Cys− in our model. We additionally removed 80% of Cys− from the incomplete, imputed dataset for the same reason.

Statistical methods

We performed statistical analyses of the environment of both Cys+ as well as Cys− to show that the chosen features should theoretically enable our models to make useful predictions. We performed all statistical tests on the not imputed dataset before pre-processing. The neighborhood of residues is one of the main predictors of post-translational modifications and catalytic activity of residues (Blom et al. 2004). Good results have been obtained for the prediction of disulfide bridges (Passerini and Frasconi 2004). No specific amino acid sequence motifs around Cys+ have been found (Greco et al. 2006; Marino and Gladyshev 2010). We applied the Mann–Whitney U test (Mann and Whitney 1947) to test for significance and corrected for multiple testing, using the Benjamini-Hochberg Procedure (Benjamini and Hochberg 1995) as Bonferroni correction.

Feature extraction and pre-processing

We pre-processed protein set one by two methods. The first method processed amino acid residues in the sequence, while the second method utilized the 3D neighborhood in the structure.

The first method extracted 414 features for each cysteine. We considered the 20 closest amino acids in the sequence, according to a threshold adapted from Passerini and Frasconi (2004). The second method extracted 70 features for each cysteine. The features were intended to describe physical properties and were computed for each cysteine and each of the 13 amino acid residues in its closest 3D neighborhood (Euclidean distance). Preliminary statistical tests showed that statistical significance decreased with larger neighborhood (data not shown).

Only the first method was used for the imputed dataset, since the 3D neighborhood of cysteines was unknown for the proteins with unknown structure.

For each amino acid, we computed 24 features:

• – 17 values of physicochemical properties of the side chain (abbreviations in parentheses): molecular mass (Mass), volume (Vol), surface area (Area), three values for the likelihood of either large, regular or small solvent exposed area (SEA1, SEA2, SEA3), three values of propensities for the SSEs, α-helix, β-strand, or turn (Alpha, Beta, Turn), and eight binary values for the chemical classification, which are polar, non-polar, charged, positive, tiny, small, aromatic and aliphatic (Polar, Non-P, Charge, Positive, Tiny, Small, Aromatic, Aliphatic). We adopted standard values for the physicochemical properties for the amino acids from Lin et al. (2011).

• – Half Sphere Exposure values HSE1 and HSE2: We applied the function HSExposure of BioPython in version 1.74. The HSE is defined here as the number of amino acid neighbors within two half-spheres with a radius of 12 Å. The sphere is divided into two halves by a plane perpendicular to the C _β -C _α vector.

• – Relative accessibility: To determine the accessibility from the 3D structure, we applied DSSP (Kabsch and Sander 1983). The accessible surface area is defined as the residue water-exposed surface in Ų. Relative accessibility is defined here as the accessible surface area divided by the maximum accessible surface area as defined by Tien et al. (2013).

• – Four values to assign an SSE to the residue according to DSSP:

  • – α-helix (H), 3-helix (G), 5-helix (I), which we counted as “helix”

  • – residue in isolated β-bridge (B), extended strand (E), hydrogen-bonded turn (T), which we counted as “strand”

  • – bend (S),

  • – loop/irregular structure (−)

The secondary structure was encoded in a four-dimensional binary vector of length one.

For the sequence based method, the residues were treated separately and in sequence order. The 17 physicochemical properties plus the HSE1, HSE2 and relative accessibility for each of the 20 closest residues sum up to 400 values. The vectors for the SSE assignments were added together into four values. The number of known PTMs (phosphorylation, acetylation or ubiquination) of residues among the 20 closest residues of the investigated cysteine were also used as three separate features. We completed the set of values by the seven features of the cysteine itself, not including the redundant 17 physicochemical properties of the cysteine. For method 2, we added all values for any of the physicochemical properties together into one value for each of the 17 properties to avoid the introduction of a sequential arrangement into the data, as amino acids in a spatial neighborhood are not ordered consecutively. The vectors for the SSE assignments were again also added together into four values. We treated HSE and relative accessibility separately for each residue, resulting in 39 additional features. Together with the values for PTMs and the secondary structure, HSE and relative accessibility for the cysteine itself, this sums up to 70 features. We collected 3D positions of residues from the PDB files of highest resolution. For the imputed dataset, all data, except for the 17 physicochemical properties of the 20 neighboring residues and the PTMs, were imputed. All values were scaled using the StandardScaler function of scikit-learn (Pedregosa et al. 2011). We applied feature selection using the SelectKBest function of scikit-learn (Pedregosa et al. 2011) with the f_classif score function to compute the ANOVA F-value of the features and reduce their number to 30. ANOVA is a statistical test for determining whether the means of two or more samples of data came from the same distribution or not. An F-test is used here to calculate the ratio between the explained and unexplained variance by a statistical test like ANOVA. Feature selection was performed independently within each cross-validation fold.

Machine learning methods and tools

For modeling and prediction of Cys+, we utilized the following machine learning techniques: SVM (Cortes and Vapnik 1995), RF (Ho 1995) and ET (Geurts et al. 2006). We applied the package scikit-learn 0.22 (Pedregosa et al. 2011) in Python version 3.7. We used the support vector classification (SVC) function with an RBF (radial basis function) kernel for SVM prediction. For all algorithms, 5-fold cross-validation was performed to avoid overfitting.

To evaluate the algorithms, we used the AUC value for the ROC curve, which displays the True Positive Rate (TPR) against the False Positive Rate (FPR) at different thresholds for a positive prediction. TPR is the ratio between the number of Cys+, which were accurately predicted divided by the full number of Cys+ in the dataset. FPR is the ratio between the number of Cys− which were falsely predicted as Cys+ divided by the full number of Cys− in the dataset. The AUC is the probability that the algorithm will rank a randomly chosen Cys+ higher than a randomly chosen Cys−. A value of one would be a perfect score, while a value of 0.5 signifies a completely random classification. For detailed information on the machine learning methods, see the section Methods in the supplemental material.