FINDSITE-metal: Integrating evolutionary information and machine learning for structure-based metal binding site prediction at the proteome level

Dataset

The metal binding proteins used in this study were obtained from the Metalloprotein site Database and Browser (MDB) ⁴⁹, which provides quantitative information on all metal-containing sites available from structures in the PDB ⁵⁰. Only proteins bound to the following eight metal ions were included in the dataset: Ca, Co, Cu, Fe, Mg, Mn, Ni and Zn. For each binding metal, a non-redundant set was compiled using PISCES ⁵¹. For proteins 50-600 residues in length, redundancy was removed at the 35% pairwise sequence identity level. The final non-redundant dataset comprises 860 proteins, of which 201, 29, 35, 117, 152, 87, 21, and 251 bind to Ca, Co, Cu, Fe, Mg, Mn, Ni and Zn, respectively. For each metal binding site, a set of binding residues was identified using the interatomic contacts provided by the LPCsoftware ⁵² with the remaining residues classified as non-binding. The list of proteins and associated metal binding ions can be found at http://cssb.biology.gatech.edu/findsite-metal/.

Protein structure modeling

For each target protein, we have constructed several models of different quality in terms of their RMSD ⁵³ from the native structure. In addition to the crystal structures, we use three sets of uniformly distorted structures with an average RMSD of 2, 4 and 6 Å from native. The distorted structures were generated starting from the crystal structures by a simple Monte Carlo procedure that deforms protein structures to a desired resolution ⁵⁴. Moreover, we apply a state-of-the-art template-based structure prediction algorithm ⁵⁵ to construct a set of weakly homologous protein models. First, for each target protein, distantly related template structures (<35% sequence identity to the target) were identified in a non-redundant PDB library by our meta-threading procedure that employs the SP3 ⁵⁶, SPARKS2 ⁵⁷ and PROSPECTOR_3 ^58,59 algorithms. Subsequently, full-length models were assembled by chunk-TASSER ⁵⁵. Finally, all-atom models from the top ranked chunk-TASSER structures were constructed by Pulchra ⁶⁰ and energy minimized in the CHARMM22 force field ⁶¹ using the Jackal modeling package ⁶². The benchmark dataset in terms of the average quality of protein structures as assessed by RMSD ⁵³, TM-score ⁶³ and MaxSub ⁶⁴ is summarized in Table I.

Template identification

FINDSITE-metal is a template-based procedure for metal binding site prediction. Here, template proteins are identified in a non-redundant PDB library using meta-threading that employs three threading procedures: SP3 ⁵⁶, SPARKS2 ⁵⁷ and PROSPECTOR_3 ^58,59. Only weakly homologous (<35% sequence identity to the target) template structures that have a Z-score of ≥4 reported by at least one threading method are included. The initial set of templates provided by threading is used to retrieve all metal binding protein structures from the PDB that are homologous to at least one threading-identified template. Multiple instances of a metal binding template protein (>90% sequence identity) are only retained if they bind either to a different metal ion or to the same metal but in a different location, with a distance of >4 Å. Otherwise, only one PDB structure is included. Again, we remove all PDB templates with >35% sequence identity to the target. Finally, only those template structures that have a TM-score to the provided target structure of ≥0.4 are retained. If distorted or modeled proteins are used as the targets for metal binding site prediction, the TM-score is calculated vs. these structures. This structure similarity threshold ensures that the template-to-target structure alignments generated by fr-TM-align ^65,66 are statistically significant.

Metal binding site prediction

Similar to the original FINDSITE approach ^7,48, FINDSITE-metal employs structure alignments provided by fr-TM-align ^65,66 to superimpose metal-binding templates detected by threading onto the target structure (either crystal or modeled). Subsequently, upon global superposition of the template structures, template-bound metal ions are clustered using an average linkage clustering procedure and the resulting clusters are ranked by the number of binding metals. Each cluster represents a putative metal binding site with the predicted metal location at the cluster geometrical center (averaged coordinates of all templatebound metals).

Binding residue prediction

For each metal binding cluster, the initial set of binding residues is calculated as follows: Each target residue is assigned a probability that corresponds to the fraction of templates that have a residue in equivalent position in contact with a metal, including pseudo counts:

where $p_i^B$ is the metal binding probability, c is the number of templates that have the equivalent residue in contact with a metal, n is the total number of templates and f_i is the frequency of occurrence of residue i in UniProtKB/Swiss-Prot ⁶⁷ (see Release notes for UniProtKB/Swiss-Prot release 56.0). Residue equivalences are calculated from structure alignments generated by fr-TM-align.

Prediction of binding metal preferences

Similarly, for each putative binding site, we calculate the preferences toward different metal ions using a fraction of templates that bind particular type of metal:

where $p_j^M$ is the probability of binding metal j (we use eight different metal types: Ca, Co, Cu, Fe, Mg, Mn, Ni and Zn), c is the number of templates that bind metal j, n is the total number of templates and f_j is the frequency of occurrence of metal j in a non-redundant subset of the MDB ⁴⁹.

The uncertainty of metal binding (ME) is quantified by the Shannon’s entropy ⁶⁸:

Low ME values are indicative of relatively homogenous metal binding sites, i.e. similar locations in evolutionarily related proteins tend to bind the same type of metal. We use ME to construct a reliable confidence index for binding metal prediction.

Machine learning

The accuracy of metal binding residue prediction is further improved by machine learning using classification-based Support Vector Machines (SVM). Here, we use libSVM 2.9 ⁶⁹ to build a C-SVC model with a radial basis function. To avoid memorization of the dataset, we use a 2-fold cross validation protocol. The complete dataset of the target complexes was randomly divided into two subsets with < 35% sequence identity between any two proteins that belong to the different subsets. Subsequently, each subset was used to train the model, and then predictions were made for the remaining targets, excluded from the training procedure. The constructed SVM model employs the set of 25 features summarized in Table II. We use the trained SVM classifier to assign each residue with a probability to bind a metal ion (probability being a positive). One of the features used in SVM is the TM-score to native estimated for the target structure. Using crystal structures, this value is 1.0. For structures distorted to 2, 4 and 6Å RMSD, we train the model on the TM-score values calculated vs. native structures; however, in validation, we use random TM-score values sampled from a normal distribution calculated using the mean and standard deviation for a given dataset of distorted structures (see Table I). For example, the estimated TM-score for each target from the set of structures distorted to 4 Å RMSD is selected randomly from a normal distribution with a mean of 0.75 and a standard deviation of 0.07. For TASSER models, the TM-score is estimated from the C-score, the confidence score calculated from TASSER simulations ⁷⁰. This is described in the following section.

TM-score estimation

Previously, the C-score was introduced as a structure prediction confidence index ⁷⁰. A positive C-score indicates that the modeled structure is very likely to be topologically similar to native at the structurally significant level. Here, C-score values are used to estimate the TM-score vs. native for protein structures modeled by TASSER. To find the correlation between the C-score and the TM-score, we use the results of large-scale benchmark simulations carried out for a standard, non-redundant dataset of proteins that cover the PDB at the 35% sequence identity. The dataset consists of 1,489 single protein chains up to 100 residues in length, 2,494 proteins between 100-200 residues and 1,203 larger proteins between 200-300 residues. Structure models are constructed by TASSER from weakly homologous template structures (<35% sequence identity) identified by threading ^58,59. The regression analysis is performed using the C-score values calculated from TASSER simulations and the TM-score values calculated vs. crystal structures of the targets.

Binding site re-ranking

The identity of metal-binding residues is typically highly conserved across a set of evolutionarily related proteins. We use this observation to re-rank binding sites by a sequence profile score, which involves the summation over all twenty amino acid types of the product of the probabilities that a given amino acid occupies the equivalent position in the target and template and is derived from structure alignments generated by fr-TM-align for a set of weakly homologous metal-binding templates. Putative binding sites, initially ranked by the number of binding metals, are re-ranked by the total sequence profile score calculated over binding residues predicted by the SVM.

Prediction confidence

FINDSITE-metal uses three confidence indexes that estimate the chances of i) the metal position to be predicted within a distance of 4 Å, ii) the Matthew’s correlation coefficient for the binding residues to be at least 0.5 and iii) the binding metal type to be correctly predicted. The prediction confidence is estimated by a Naïve Bayes classifier ⁷¹ from a set of features listed in Table II. A separate classifier is trained for each confidence index. Similar to the binding residue prediction by SVM, we use a 2-fold cross validation protocol. We find that the raw instance scores returned by the Bayesian classifier make rather poor confidence estimates, because they are grouped around the extreme values (0 and 1). Therefore, we apply a calibration procedure to normalize confidence estimates generated by the classifier. Here, we use the Pool Adjacent Violators (PAV) algorithm ⁷² that transforms the raw scores into well-calibrated posterior probabilities, which are further used as the confidence estimates.

Comparison to a sequence-based method

We compare the performance of FINDSITE-metal to SVM-Prot, an SVM-based method that predicts the functional class of metal-binding proteins from sequence derived physicochemical properties ⁷³. Here, we use the same non-redundant MDB dataset, which is described in the previous sections. To make results comparable between SVM-Prot and FINDSITE-metal, for the latter, we employ weakly homologous TASSER models as the target structures. Moreover, we use only distantly related (<35% sequence identity to target) metal-bound template structures identified by meta-threading. In this manner, the predictions are made by both approaches solely using sequence information as the input. The SVM-Prot server at http://jing.cz3.nus.edu.sg/cgi-bin/svmprot.cgi was queried automatically and the results were parsed for metal-binding assignments. The class of binding metal is selected based on the p-value reported by SVM-Prot. For FINDSITE-metal, the metal type is selected based on the highest preference calculated using Eq. 2. The results are assessed separately for each metal type (Ca, Co, Cu, Fe, Mg, Mn, Ni and Zn) using a standard receiver operating characteristic (ROC) analysis.

Prediction of metal-binding sites in the human proteome

Amino acid sequences of all gene products identified in the human proteome were obtained from the Ensembl genome database ⁷⁴. Here, we use 56,376 protein sequences 50-600 residues in length selected from the human assembly GRCh37, release 55. For each sequence, a Cα backbone model was built by TASSER ^70,75 from template structures identified by SP3 ⁵⁶, SPARKS2 ⁵⁷ and PROSPECTOR_3 ^58,59. Subsequently, all-atom models reconstructed by Pulchra ⁶⁰ from low-resolution TASSER structures were subject to short energy minimization using Jackal ⁶². The set of meta-threading identified templates and the modeled structures were then used by FINDSITE-metal to detect putative metal-binding sites in the human proteome. Each predicted binding site was assigned a confidence and further characterized by the prediction of metal-binding residues, the class of binding metal and the molecular function.

Metal-binding templates

FINDSITE-metal employs a set of evolutionarily related metal-binding templates selected by threading. The average number of templates per target is 18. Upon the global superposition onto the target structure, putative metal-binding sites are detected by a clustering procedure. First, for a given set of templates, we analyze what is the optimal clustering cutoff in terms of the average distance from the metal position in the target structure and the ranking ability by the cluster multiplicity. In Figure 1, we show that the average distance of the best binding site increases with the clustering threshold. Small cutoff values result in many puny clusters, one of which is typically close to the native metal binding site; however, this is at the expense of a poor ranking capability (Figure 1, inset). On the other hand, if a large cutoff is used, the ranking becomes more efficient, but the average distance from the natively bound metal position increases. As a trade-off between the accuracy of the site prediction and the ranking ability, we use a clustering cutoff of 8 Å in further calculations.

Analysis of metal binding residues

Metal binding in proteins typically requires a specific geometrical arrangement of relatively few residues, whose identity strongly depends on the type of binding metal ^13,76-78. This is shown in Figure 2, where we compare the amino acid preferences to bind different metals. The electron donors in the side chains are mainly the carboxyl oxygen atoms of Asp and Glu, the imidazole nitrogen of His and the thiol group of Cys. Moreover, the amide nitrogen and oxygen of Asn and the thioether group of Met also coordinate divalent metal ions. Calcium, magnesium and manganese ions preferentially bind to the acidic chain of Asp. Zinc, with a lower coordination number preference, is typically chelated by Cys and His. Histidine residues also have strong preference toward binding of cobalt, copper, iron, nickel and zinc atoms. These results correlate very well with the recent statistical analysis of the chemical environment of metal binding in proteins ⁷⁹. Binding patterns are potentially important for the prediction of metal binding residues. If the type of metal is correctly predicted, the differential metal binding preferences of amino acid side chains can be used to increase the accuracy of binding residue prediction.

Binding metal preferences

In this section, we show that evolutionarily remotely related proteins tend to bind similar metals in equivalent locations. For each target protein and the corresponding set of threading identified templates, we calculate the preferences toward binding native and non-native metal ions. The native metal preference is equivalent to the fraction of templates that bind the same metal as the target structure. Likewise, the non-native preference is an averaged fraction of templates that bind different metal types in similar locations. Figure 3 shows that the native metal preference is strongly correlated with the distance between the target- and template-bound metal ions upon the global superposition of their structures. For site distances less than 2 Å (4 Å), more than 70% (50%) of the templates bind the same metal type as the target structure. This fraction drops dramatically for sites >4 Å away from each other. Encouragingly, this tendency is observed not only for crystal structures (Figure 3A), but also for the distorted target structures (Figure 3 B-D). As we demonstrate in the following sections, combining information on the predicted binding metal with the differential metal binding preferences of amino acids improves the accuracy of the prediction of binding residues particularly against modeled protein structures.

TM-score estimation

Before we discuss the performance of FINDSITE-metal in metal binding site prediction, we shortly describe the results of an analysis that focuses on the relationship between the C-score and the TM-score. The TM-score provides a length-independent measure of the structural similarity between two proteins ⁶³. A significant similarity is indicated by a TM-score of >0.4. FINDSITE-metal uses the estimated TM-score to native as one of the SVM features to accurately predict metalbinding residues and to estimate the prediction confidence. For a given protein model, its TM-score can be directly calculated against the native crystal structure, if known. However, in a real prediction scenario when the experimental structure of a target is unavailable, the TM-score needs to be estimated. Most contemporary structure prediction algorithms estimate the reliability of structure modeling using some score ^80-82. TASSER, a template-based structure assembly/refinement approach ^70,75, calculates a confidence score, called the C-score. To estimate the TM-score for a target structure modeled by TASSER, we carried out a regression analysis using the results of large-scale benchmarks on a non-redundant and representative dataset of 5,186 proteins. For each modeled target structure, we calculate the TM-score vs. its crystal structure and plot it against the C-score obtained from TASSER simulations; this is shown in Figure 4. Next, we calculate the regression line, which in this case is:

with a Pearson correlation coefficient of 0.81. We use this equation to estimate the TM-score for target structures modeled by TASSER. A high correlation between the TM-score and a modified version of the C-score was also reported for I-TASSER ⁸³. We note that FINDSITE-metal does not require the exact TM-score; all benchmark results reported here were obtained using TM-score estimates rather than the exact values (see Materials and Methods). In principle, any other model quality assessment that correlates with the TM-score can also be used.

Accuracy of the predicted metal position

FINDSITE-metal predicts putative metal binding locations by the clustering of template-bound metal ions upon the global superposition of the template structures. In Figure 5, we assess how far is the predicted site center from a metal position in the crystal structure, when the crystal structures themselves as well as structures distorted to a desired RMSD are used as the targets. Using a distance of 4 Å as the hit criterion, the fraction of correctly predicted sites at the top five ranks is 69.5%, 67.2%, 58.9% and 50.8% for the crystal structures and structures distorted to 2, 4 and 6 Å RMSD, respectively. The difference between crystal and distorted structures diminishes with more promiscuous distance thresholds; for 6 Å (8 Å), the difference between the fraction of correctly predicted metal binding sites using crystal and the most distorted structures is only 8.1% (3.9%). This is a common feature of many template-based predictors that employ structure alignments, which are sensitive to the global topology rather than to the local structural features ^7,31,48.

The ability of FINDSITE-metal to correctly rank the predicted binding sites is assessed in Figure 5 (inset) for two ranking protocols. Similarly to the original FINDSITE algorithm, predicted binding locations could be ranked by the fraction of templates that share a common site. Using ranking by fraction, in 67.8%, 67.0%, 66.8% and 66.5% of the cases, the best predicted site is at rank 1 for the crystal structures and structures distorted to 2, 4 and 6 Å RMSD, respectively. As demonstrated above, binding metal types as well as residue profiles are remarkably strongly conserved across a set of evolutionarily related proteins. Here, we use this observation to re-rank the predicted sites by a sequence profile score calculated vs. metal-bound templates. Figure 5 (inset) demonstrates that ranking by a sequence profile score calculated for binding residues predicted by SVM (see Materials and Methods) yields 3-4% improvement with respect to ranking by fraction. Furthermore, ranking ability is hardly affected by the distortion of the target structure; using the crystal structures and structures distorted to 2, 4 and 6Å RMSD, the correct ranking by the sequence score is found in 70.7%, 70.6%, 70.3% and 69.3% of the cases, respectively.

Binding residue prediction

Many algorithms for metal binding residue prediction use the fact that metal binding typically involves a limited set of strongly conserved binding patterns 30,33,84. Here, we analyze the performance of FINDSITE-metal in the prediction of binding residues using two classifiers: a probability-based residue selection and machine learning that employs a set of sequence and structure derived features. A probability-based approach simply assigns a residue with a binding probability that corresponds to the fraction of templates that have a residue harboring a metal ion in the equivalent position. Figure 6 demonstrates that the overall accuracy of binding residue prediction is significantly improved when the SVM classifier is applied. For example, at the cost of a 1% false positive rate, machine learning increases the true positive rate by 7% (from 84% to 91%), 9% (from 81% to 90%), 11% (from 75% to 86%) and 12% (from 69% to 81%) for the crystal structures and structures distorted to 2, 4 and 6 Å RMSD, respectively. Similar improvement is observed in the Recall-Precision graphs shown as the inset plots in Figure 6.

Prediction of binding metal

As demonstrated in the previous sections, similar binding sites across a set of evolutionarily related proteins bind similar types of metal ions. We use this observation to predict a metal that likely binds to the detected binding sites. The fraction of dataset targets for which the binding metal was correctly predicted is shown in Figure 7, separately for each metal class. The highest accuracy, 70-90%, is observed for Fe, Cu, Zn, Ca and Mg binding proteins. The fraction of proteins correctly predicted to bind Mn, Co and particularly Ni is significantly lower; this is caused by the underrepresentation of these proteins in the PDB ^50,79. Moreover, only very distant homologues are used in this study that might have evolved to bind a different class of metal ions. It is noteworthy that the accuracy of binding metal prediction is highly insensitive to the structural distortions of the target structures.

Predicted amino acid/metal composition

Analyzing all predicted binding residues and metal ions across the non-redundant dataset, we can estimate the rate of over- and under-prediction of certain amino acids and metal types. Overall, we find that the predicted residue as well as metal composition is in an excellent agreement with these calculated directly from the PDB complexes; this is shown in Table III. Typically, the amino acid composition difference is less than a couple of percent, with the exception of Cys that is slightly under-predicted by 2.6-3.7%. For low quality protein structures, particularly these that are distorted to 4 Å and 6 Å RMSD, we observe a moderate composition excess of Asp, Glu and His residues of 1.0-5.5%. Considering the metal composition across the dataset, most binding metal types are correctly predicted with two exceptions: magnesium-binding sites that are over-predicted by 7.9-8.6% and manganese-binding sites are under-predicted by 6.3-7.0%. Nevertheless, the residue/metal composition of metal binding sites is fairly well reproduced by FINDSITE-metal predictions; this feature is important for large-scale functional assignments at the level of entire proteomes.

TASSER models as targets for FINDSITE-metal

Artificially distorted structures provide some notion about the performance of a method using different quality target structures. However, from the point of view of real applications, the most interesting results are these obtained using protein structures modeled by a state-of-the-art protein structure prediction approach. In this study, using TASSER, we constructed protein models from only weakly related template structures. In Figure 8, the performance of FINDSITE-metal using TASSER models is compared in terms of the site location accuracy and ranking capability to that using the crystal structures of the targets. Unlike other structure-based approaches to metal binding site prediction that strongly rely on the quality of the target structure ^85,86, FINDSITE-metal is insensitive to some extent to the structural distortions in modeled proteins structures. Using a 4 Å (8 Å) distance as a hit criterion, the accuracy drops by 10% (5%) if weakly homologous TASSER models are used instead of the crystal structures, with the ranking ability reduced by 4.5%. Furthermore, the high accuracy of binding residue identification is retained; this is shown in Figure 9 as a ROC plot and a Recall-Precision graph (inset). We also find a good agreement between the number of binding residues per site in the crystal structures of the complexes (3.13 ±0.99) and the number of residues predicted for TASSER models (3.69 ±1.93), see Table I.

Confidence indexes

Due to the inherent limitations of many template-based approaches to functional annotation, such as the unavailability of suitable templates or the possibility of severe topological inaccuracies in the structures modeled using remote homology, confidence indexes are required to assess the prediction reliability. FINDSITE-metal employs three confidence estimates for the predicted site distance, the set of identified binding residues and the class of binding metals. As described in the Materials and Methods, these indexes are calculated using calibrated Bayesian classifiers. In Figure 10, we present the performance of FINDSITE-metal, as assessed by the fraction of templates for which a correct prediction was made, for targets assigned with different confidence values: the binding site distance ≤4Å, the Matthew’s correlation coefficient, MCC, ≥0.5 and the confidence that the binding metal is correctly identified. Our confidence indexes correlate well with the prediction accuracy not only for the target crystal structures, but also when distorted structures and TASSER models are used as the targets. We note that the confidence is estimated for a given target without using any information on the experimental structure of the metal-protein complex.

Comparison to SVM-Prot

We compare the performance of FINDSITE-metal to SVM-Prot, a sequence-based predictor of the functional class of metal-binding proteins. In this analysis, the focus is on the prediction of a metal-binding type. The results for Ca, Co, Cu, Fe, Mg, Mn, Ni and Zn binding proteins are presented in Figure 11 and Table IV. FINDSITE-metal achieves a relatively high sensitivity of >50% for five out of eight functional classes (Fe, Cu, Zn, Ca and Mg) at a moderate false positive rate below 20%. For the remaining metal-binding proteins, the assignments are marginally better than random. The accuracy of SVM-Prot is notably better than random only for targets that bind to Fe, Zn and Mn. The sensitivity of SVM-Prot for zinc is higher than using FINDSITE-metal; however, with a much higher false positive rate, which indicates a significant over-prediction of zinc-binding proteins. SVM-Prot is also more sensitive in detecting manganese-binding targets (29% sensitivity at 5% false positive rate). Nevertheless, considering the overall accuracy of the functional class assignment, FINDSITE-metal represents a significant improvement over SVM-Prot. Since the benchmarking results reported here seem encouraging, below we describe the application of FINDSITE-metal to the entire human proteome.

Metal binding complement of the human proteome

In this study, we constructed structural models for 56,376 gene products in the human proteome, 50-600 residues in length. 34,808 of these were assigned with at least one metal binding site. The distribution of the estimated TM-score for putative metal-binding proteins in the human genome is shown in Figure 12. A TMscore of ≥0.4 was assigned to 70.7% of the targets (24,617 gene products). Structural models for these sequences are very likely to be correct, at least at the topological level. Since FINDSITE-metal tolerates to some extent structural inaccuracies in modeled structures, these results suggest that reliable predictions can be made for the majority of proteins.

The prediction confidence for the metal binding site prediction in the human proteome is presented in Figure 13. For roughly one third of the gene products, the estimated confidence that the distance of the top-ranked site is predicted within 4 Å and the binding residues are identified with a Matthew’s correlation coefficient of at least 0.5 is higher than 40-50%. For three quarters of the targets, the confidence of binding metal prediction is >50%. Considering the large number of targets, FINDSITE-metal provides confident metal binding information for thousands of gene products in the human proteome.

Figure 14 shows the statistics on the assigned metal binding class to the human metalloproteome. Nearly one third of putative metal binding proteins in the human proteome were predicted as being calcium-binding and another 30% as magnesium-binding. The third largest class (one-quarter of putative metalloproteins) consists of proteins that bind to zinc. Fe, Mn, Co, Cu and Ni were assigned to 5.4%, 3.6%, 1.9%, 1.4% and 1.2% of the targets, respectively. Benchmark results reported in the previous sections suggest that the number of proteins that bind to Mn, Co and Ni might be somewhat underestimated. Nevertheless, the composition of the metal binding complement of the human proteome identified by FINDSITE-metal is consistent with other studies ^23,79.