Exploring functionally related enzymes using radially distributed properties of active sites around the reacting points of bound ligands

Background Structural genomics approaches, particularly those solving the 3D structures of many proteins with unknown functions, have increased the desire for structure-based function predictions. However, prediction of enzyme function is difficult because one member of a superfamily may catalyze a different reaction than other members, whereas members of different superfamilies can catalyze the same reaction. In addition, conformational changes, mutations or the absence of a particular catalytic residue can prevent inference of the mechanism by which catalytic residues stabilize and promote the elementary reaction. A major hurdle for alignment-based methods for prediction of function is the absence (despite its importance) of a measure of similarity of the physicochemical properties of catalytic sites. To solve this problem, the physicochemical features radially distributed around catalytic sites should be considered in addition to structural and sequence similarities. Results We showed that radial distribution functions (RDFs), which are associated with the local structural and physicochemical properties of catalytic active sites, are capable of clustering oxidoreductases and transferases by function. The catalytic sites of these enzymes were also characterized using the RDFs. The RDFs provided a measure of the similarity among the catalytic sites, detecting conformational changes caused by mutation of catalytic residues. Furthermore, the RDFs reinforced the classification of enzyme functions based on conventional sequence and structural alignments. Conclusions Our results demonstrate that the application of RDFs provides advantages in the functional classification of enzymes by providing information about catalytic sites.


Background
High-throughput methods for structural genomics have produced an increasing number of protein structures to be solved by X-ray crystallography. The abundance of protein structure information in the Protein Data Bank (PDB) has increased the need and desire for structurebased function prediction [1] and has contributed to structure-based drug design [2]. However, two problems remain regarding the prediction of enzyme function. First, proteins within a superfamily, which are usually expected to share the same catalytic properties, can catalyze different reactions. There are reports that enzymes with 98% sequence identity, such as melamine deaminase and atrazine chlorohydrolase, may catalyze different reactions [3]. Second, two enzymes belonging to different superfamilies or fold classes can catalyze almost identical reactions [4].
The function of a protein can be affected by a small number of residues in a localized region of its three-dimensional structure [5]. Moreover, the specific arrangement and conformation of these residues can be crucial to a protein's function and may be strongly conserved during its evolution, even when the protein sequence and structure change significantly [5]. For example, it was reported that the positioning of the reactive region of a substrate with respect to a cofactor is generally conserved in flavoenzymes [6].
Two methods for the description of local structures have been developed for predicting enzymatic functions. First, in the element-based description of catalytic residues, the catalytic roles in an enzymatic reaction are defined as acid-base, stabilizer or modulator roles [7]. Some insight into enzymatic reactions can be gained using this method, but manual annotation is inherently required. In addition, it is often difficult to differentiate between the acid-base and stabilizer roles because most structures solved by X-ray crystallography provide no information about hydrogen atoms. The second method is based on descriptions of substructures within the local structures of enzymes [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. Many approaches to analyze and compare local structures have been proposed. One group of algorithms, which includes the PINTS [8], ETA [9][10][11] and FLORA [12] algorithms, scans protein structural databases using pre-calculated or automatically generated templates. Another group includes algorithms that compare the substructural epitopes of proteins using geometric hashing [13][14][15]. Similarly, SiteEngine [16] uses the concept of pseudocenters [17] to define the properties of the corresponding surface. None of these approaches can characterize catalytic sites and create feature vectors, even though they assess the similarity between catalytic sites.
In this study, we examine the structures of oxidoreductases and transferases using radial distribution functions (RDFs) that encode radially distributed properties of active sites centered around the reacting points of bound ligands. Thus, element-based and substructure descriptions are integrated into the RDF, assuming that catalytic roles are restricted by distances and that different catalytic residues can play identical roles. Although the topological correlation vector method of Stahl et al. [18] and WaveGeoMap, developed by Kupas et al. [19], provide feature vectors related to enzyme cavities, these descriptions use patches of active sites, regardless of the orientation of the catalytic residues. Therefore, it is still unclear whether the orientation of active sites around a reacting point is related to enzymatic function and how much of the orientation is conserved. Our method provides a different view of enzymatic function by focusing on the physicochemical properties surrounding a reacting point found in enzyme cofactors.

Characteristic physicochemical pattern of active sites
To examine how catalytic residues contribute to the radially distributed properties of active sites, we decomposed the RDF into the total charge for each residue. Figure 1 shows the contributions made to the peaks and minima of the RDFs by the various catalytic residues surrounding the carbon atom (C4N) of the nicotinamide adenine dinucleotide (NAD) molecule in 1dc6 and the iron atom (FE) of the heme (HEM) molecule in 1sog (PDB). The first local minimum in 1dc6 corresponded to the nucleophilic cysteine residue and the asparagine residue that binds NAD ( Figure 1A). The second local minimum and the last two minima were affected by the threonine residue that binds the substrate ( Figure 1A). The peak at 6 Å corresponded to the histidine residue that activates the thiol group of the cysteine residue and to the cysteine residue itself ( Figure 1A). The last two peaks were derived from the arginine residue that binds the substrate ( Figure 1A). All of these five residues in 1dc6 are known to be critical for the enzymatic reaction. The first peak in 1sog corresponded to the proximal histidine residue (the heme axial ligand) and the distal histidine residue (the proton acceptor) ( Figure 1B). The subsequent minimum at 7.5 Å was slightly affected by the histidine residues and the arginine residue (a transition state stabilizer) ( Figure 1B). The small peak at 8.5 Å shown in Figure 1B was derived from the tryptophan residue (a radical intermediate). We can show that all of the four residues described in 1sog also play an important role in the catalytic reaction. Moreover, the degenerated total charge of the catalytic residues corresponded to the RDF with a range from 0 to 5 Å for both enzymes ( Figure 1C, D). The bias of the RDFs toward a negative charge may be due to ignoring hydrogen atoms. These results show that catalytic residues are primarily responsible for the physicochemical properties of active sites.
Subsequently, to investigate whether the RDFs adequately discriminate between active sites, we selected pairs including the wild-type and a mutant form of the enzymes from the dataset. We then calculated the averaged Euclidean distances or cosine distances (1 minus the cosine similarity) of matched pairs (wild/wild or mutant/mutant) and mismatched pairs (wild/mutant) for each enzyme compared to the distance measure (the complement 100 minus the match score) obtained using SiteEngine ( Table 1). As shown in Table 1, most of the pairs were agglomerated within a Euclidean distance of 222, and the RDFs were successful in revealing their similarity. The pairs from MDH_ECOLI were located approximately at a Euclidean distance of 322 from each other; however, two identical mutants were reported to have different conformations (PDB code: 1ib6 and 1ie3) [24]. Long distances were observed between the mismatched pairs compared to those of the matched pairs (Table 1). Of the mismatched pairs, the CCPR_YEAST proteins were particularly close to each other (< a Euclidean distance of 199), and the RDFs failed to identify their dissimilarity. However, the conformation of the catalytic site in the CCPR_YEAST protein is not altered by the mutation included in this analysis (PDB code: 3ccp) [25]. These results show that the Euclidean distance between the RDFs reflects the conformational changes in the active sites. The match scores from SiteEngine were similar to the distances between the RDFs. Although the distances between the RDFs were slightly poorer than the match scores in terms of the number of successful discriminations, the feature vector of the RDF is almost equivalent to the measure of SiteEngine in discriminating between the active sites. Thus, the active sites are characterized based on the physicochemical patterns of the RDFs.

Active site properties as the critical determinants of enzyme function
To investigate whether the RDFs account for a major part of the enzyme function, clustering of the RDFs was performed using a self-organizing map (SOM) approach. Figure 2 shows the results for glyceraldehyde-3-phosphate dehydrogenase (GAPDH) and cytochrome c peroxidase (CCP). The GAPDH residues were mainly distributed in the area around node [39,6], including the two different catalytic sites ( Figure 2). Within the GAPDH distribution, 1 dc6 from Escherichia coli and 1nq5 from Bacillus stearothermophilus (PDB) were closely positioned at nodes [38,5] and [38,9], respectively. The only difference between the catalytic sites in this orthologous pair is the replacement of cysteine 149 with serine leading to a 10 4 -fold reduction in dehydrogenase activity [26] ( Figure 3A).
Similarly, the CCP residues were mainly localized in the area around node [33,10], including the two different  The w/w and m/m columns show wild-type/wild-type or mutant/mutant pairs. The w/m columns show wild-type/mutant pairs. The results with statistically significant differences between the match and mismatch are shown in bold font. The statistical significance was assessed by Wilcoxon rank sum tests with a 5% significance level.
catalytic sites ( Figure 2). Within the CCP distribution, 1sog and 1dso from Saccharomyces cerevisiae (PDB) were positioned at nodes [36,8] and [34,13], respectively. In the active site of 1dso, histidine 175 is replaced by glycine ( Figure 3B). Thus, the results show that the obtained clusters of enzymes consist of clusters of their catalytic sites, suggesting that the RDFs of active sites account for a major part of the enzyme function.

Prediction of enzyme functions based on the physicochemical properties of active sites
In this study, we sought to identify functionally related enzymes by clustering active sites. First, we utilized the EC number for assignment of RDFs to analyze the SOM clustering. An overview of the SOM is summarized in Additional files 1-4, for Additional file 1: Tables S1, Additional file 2: Table S2, Additional file 3: Figure S1, Figure 2 Nonlinear projection of RDFs. The SOM was run using an RDF with an Epanechnikov neighborhood function in a [46, 28]-sized rectangular lattice (left) and a magnified section (right). Following training, each node was colored according to the enzymes or catalytic residues in the RDFs that were mapped onto it. The size of the squares or triangles indicates the relative frequency of the mapped RDFs.
Additional file 4: Figure S2. Tables 2 and 3 show the division of the RDFs into nodes, each of which is labeled with its EC number. The partial RDFs labeled with the EC numbers indicated a well-defined segregation, discriminating among 76% of the EC numbers of oxidoreductases and among 55% of the EC numbers of transferases.
Then, to evaluate how many of the active sites are associated with enzyme functions, we performed a statistical analysis of the results of the SOM clustering. The averaged F-measure of all of the assigned EC numbers of oxidoreductases was 0.87, ranging from 0.22 to 1.00. Over 88% of the active sites of oxidoreductases were assigned to an EC number (see Additional file 5: Table S3). Similarly, the averaged F-measure of all of the assigned EC numbers of transferases was 0.88, ranging from 0.33 to 1.00. Over 88% of the active sites of transferases were assigned to an EC number (see Additional file 6: Table S4).

Prediction performance in comparison with sequence and structural alignment-based annotation
To clarify the contribution of the RDFs to the functional annotation of the enzymes, we examined the relationship of the RDFs with different measures, such as sequence and structural alignment. First, we performed statistical analyses of these measures. Tables 4 and 5 show the partial correlation coefficients between the SOM distance, active site distance, local and global sequence similarities and structural similarity. The SOM distance among the RDFs was distinct from the other methods.
Next, the SOM distances among the RDFs were evaluated for their ability to annotate enzyme function in datasets that had not been correlated with known functions by either structural or sequence alignments due to pairwise identities below 25%. Tables 6 and 7 show the area under curve (AUC) values of the SOM distances for these datasets. In oxidoreductases, these values, ranging from 0.729 to 0.746, represented higher performance compared to the values obtained using sequence and structural alignments (Table 6). In transferases, the AUC values of 0.800 and 0.790 for the datasets with pairwise identities below 15% also represented higher performance compared to the values obtained using sequence and structural alignments (Table 7). These results showed that the SOM distance predicts enzyme function, even for enzymes with weak conventional similarities. Moreover, the SOM distance outperformed the match score of the SiteEngine based on substructure.
We then confirmed the ability of detecting enzymes with pairwise identities below 25%. While the ETA detected 63 oxidoreductases and 65 transferases, the numbers of enzymes assigned to the nodes within the SOM distance of 5 nodes were 454 of oxidoreductases and 387 of transferases, suggesting that the coverage of the SOM detection was higher than that of the ETA (Table 8).

Structural genomics prediction
To perform a blind validation for proteins with unknown function, we used the SOMs trained by oxidoreductases and transferases to predict enzyme functions of 102 proteins in structural genomics. While the coverage of the ETA predictions was 31%, the SOM predictions covered 57% of the query structures (Table 9). Of the predicted EC numbers, the rates of validated prediction that the EC number is compatible with the bound ligands were 59% of the ETA predictions and 72% of the SOM predictions, suggesting the SOM predictions provide a clue to annotate these functions (Table 9).

Discussion
Without using any templates, the RDFs centered around active sites are capable of clustering oxidoreductases and transferases based on their function. In this study, we applied our method only to the oxidoreductase and The numbers indicate the RDF counts assigned to the nodes, and the number of classes is shown in parentheses. The SOM was performed by the RDF with an Epanechnikov neighborhood function in a [46, 28]-sized rectangular lattice *One class is more than 80% of the total. } The nodes were labeled using SCOP [44]. The numbers indicate the RDF counts assigned to the nodes, and the number of classes is shown in parentheses. The SOM was performed by the RDF with a cut-Gaussian neighborhood function in a [40,19]-sized rectangular lattice. *One class is more than 80% of the total. } The nodes were labeled using SCOP [44]. transferases classes of proteins. We focused on these classes for the following reasons. First, oxidoreductases exhibit a great variety of catalytic sites compared to other known classes, possibly because the redox potential is modulated by oxidoreductases. Second, the reaction centers are well-defined in oxidoreductases and transferases, consisting of a substrate and cofactor that mechanistically exchange electrons and protons. The catalytic residues are generally capable of assisting in the migration of protons from the reaction center, a role that strongly resembles the roles of other enzyme classes. For example, caspase-1 is a hydrolase that catalyzes the hydrolytic reaction of peptides; the cysteine residue nucleophilically attacks the substrate, which is followed by protonation of the histidine [27]. This mechanism closely resembles the catalytic behavior observed for GAPDH, an oxidoreductase [28]. Figure 4 shows the pattern of the Cys-His catalytic diad in 1bmq was similar to that in 1dc6. The peak shift may be due to the different position between substrate and cofactor. These similarities suggest that our method can be applied to other enzymes to predict additional protein functions. To apply our method to other enzymes, the reaction centers will require manual annotation. Reaction pairs published by KEGG RPAIR [29] that include candidate reaction centers are available for other enzymes and can be used for this purpose. Understanding the orientation of catalytic sites is important for drug design. For a given G protein-coupled     [30]. The RDFs describe the orientation of catalytic sites, detecting conformational changes as well as enzyme function (Table 1). In addition, the description of the microenvironment produced by the RDF is better than simple superposition of catalytic sites when a particular functional group is not present (Figure 3). In structural genomics, the RDFs would be advantageous for finding remote orthologs, especially when evolutionary pressure has enhanced sequence/structural divergence. Although sequence-based methods are the first choice for functional annotation, proteins with sequence identities of < 20-35% are problematic [31]. Measuring structural similarity is more informative for enzyme functions exhibiting distant relationships and/or convergent evolution. However, proteins within wellknown superfamilies sharing the same structural topology, such as TIM barrels, do not always have the same functions [32]. In these cases, the measure of structural similarity alone does not correspond to functional similarity. Therefore, a specific measure representing functionality is desirable. We focused specifically on the local features around the catalytic site. Compared to the structural alignment, the functional annotation was reinforced by focusing on the reaction center (Tables 6 and 7). It is also likely that convergent evolution of an enzyme function depend less on evolutional process than on physicochemical properties of active sites (Tables 8 and 9). For proteins with unknown function, 41% of query structures were newly classified into the EC numbers (Table 9). However, the true performance of our method will be evaluated by revealing the actual function of those proteins. The combination of results obtained using different approaches will also improve the accuracy of function predictions.

Conclusion
We propose a novel classification method for the prediction of enzymatic function based on the physicochemical properties of catalytic sites. The RDFs for predicting enzymatic functions are thus far limited to enzymes with bound ligands. For ligand-unbound structures, either homology modeling or superposition based on ligandbound structures can be applied to our method. Our results suggest that the RDF provides a different perspective compared to structural and sequence alignments by focusing on a local feature because catalytic sites are thought to be more highly conserved than the overall sequences or structures of enzymes.

Dataset of active sites
Two sets of 1,880 oxidoreductase (EC1) and 789 transferase (EC2) protein structures were initially obtained from the PDB. In the case of NMR data, we used the first model in the PDB file. To simplify the filtering of the candidate active sites, structures including at least one cofactor or analogous compound were manually selected based on the annotation of PDBsum [33]. In this study, we used the substructures within 10 Å from the reaction centers of these cofactors as active site data. The reaction centers [34] of the cofactors are extensionally defined as follows: (1) atoms associated with bond formation and cleavage; (2) atoms exhibiting a change in charge; and (3) corresponding atoms in analogous compounds (see Additional files 1 and 2, Additional file 1: Tables S1, Additional file 2: Table S2). In oxidoreductases, a cofactor generally forms a part of the reaction center, acting as a donor and acceptor. Finally, based on this definition, 4,092 oxidoreductase and 1,444 transferase active sites corresponding to reaction centers were obtained. The subsequent encoding for comparison of active sites also used the Cartesian coordinates of these reaction centers as a starting point. In addition, a set of 102 protein structures with the key words of "structural genomics" and "unknown function" in the PDB was used for a blind validation of function prediction.

Characterization of physicochemical properties of active sites
The values of physicochemical atomic properties, including the main chain of the amino acid residues, were empirically calculated by the PETRA server [35,36]. The atomic properties included were the total charge for electrostatic interactions and σ-electronegativity, π-electronegativity and effective atom polarizability for van der Waals interactions. These properties are based on the Partial Equalization of Orbital Electronegativities (PEOE) [35], which is independent of 3D structures. Because the side chains of proteins show various conformations, PEOE is suitable for describing their properties.

Physicochemical encoding of active sites for the RDFs
The RDFs integrate the Gaussian distributions proportional to a physicochemical property at a distance from a starting point. Encoding of the RDF was performed by the method of Aires-de-Sousa et al. [37] with slight modification, as described below. The RDF as a function of the distance, r, is given by the following equation: where N is the number of atoms in the active site residues; r i is a constant for the inter-atomic distance between atom i and the reaction center atom (see Additional file 1: Table  S1); σ 2 is the fluctuation of the atoms around their averaged positions; and p is an atomic property (see Additional file 7: Figure S3). Thus, the RDFs naturally combine active site structures and their physicochemical properties, which exhibit an isotropic and rotationally invariant nature. In addition, we tested the effect of large σ 2 in the RDFs to investigate the robustness to conformational change, suggesting that the RDFs were robust over a large range of B-factor (= 8π 2 σ 2 /3) in the PDB (see Additional file 8: Figure S4).

SOM clustering and SOM distance
SOMs provide a topology-preserving map using a nonlinear projection of high-dimensional data onto a lowdimensional grid [38]. The low-dimensional grid is composed of nodes that represent data clusters. The neighboring nodes are connected to each other in the sense that they receive similar updates. Hence, SOMs provide information on the similarity between nodes. The SOM was run using a batch algorithm with an Epanechnikov or cut-Gaussian neighborhood function and an initial update radius of 5 or 10 nodes via implementation in the SOM Toolbox for Matlab (Mathworks, Inc.), which was developed in the Laboratory of Computer and Information Science of the Helsinki University of Technology.
In addition to the clustering, we also defined the SOM distance, which is the Euclidean distance between the SOM locations of the nodes on the grid, to obtain the distance measure between the active sites encoded by the RDFs.
Software for the alignment of sequences, structures and active sites for comparative experiments The sequences and structures were aligned using the Smith-Waterman algorithm [39] or the Needleman-Wunsch algorithm [40], both of which are implemented in the EMBOSS program package [41], or the structurebased alignment algorithms in the MAMMOTH program package [42]. All of the pairwise alignments were performed with the default parameters. The active sites were compared using a geometric hashing algorithm implemented in SiteEngine [16].

Evaluation of SOM clustering
The F-measure is defined as a harmonic mean of both precision and recall that measures the extent to which a cluster contains only enzymes of particular EC classes and all enzymes of that EC class. A cluster was defined as all nodes labeled by an identical EC class. For a particular node in the SOM, we can calculate the centroid by finding the arithmetic mean of all of the RDFs. If an RDF in the centroid vector has a high value, then the corresponding EC class occurs frequently within the node. These EC classes can be used as labels for the node. The F-measure of a cluster with respect to an EC class was defined by the following equation: The precision of a cluster with respect to an EC class was defined as follows: where M is the number of enzymes in a specific cluster, and m is the number of enzymes of the specified EC class in the cluster. Recall is the extent to which a cluster contains all of the enzymes of a specified EC class. The recall of a cluster with respect to an EC class was defined as follows: where N is the number of enzymes in the EC class. The averaged F-measure for the validation of the classification performance was obtained by calculating the average of all of the EC classes, with 1 being the best value and 0 being the worst value.

Evaluation of the measures for predicting enzyme functions
To estimate the degree of separation between two different functions when using a certain pairwise measure, analysis of a receiver operating characteristic (ROC) curve for the SOM distance was performed as well as local, global and structural alignments. Based on the cutoffs that determine whether the protein-protein pairs are predicted to be involved in the same function, i.e., true (match) or false (mismatch), the data are divided into true positives (TP), false positives (FP), false negatives (FN), and true negatives (TN). The true positive rate (TPR) and false positive rate (FPR) are defined as follows: and The ROC curve is a graphical plot of TPR versus FPR, showing the fidelity of discrimination at varying thresholds. The AUC is defined as the area under the ROC curve, representing the overall performance of discrimination. In this study, the SOM distances represented the dissimilarities among the RDFs. In the alignments, the similarities were the percentages of the number of