Functional classification of protein structures by local structure matching in graph representation

Abstract As a result of high‐throughput protein structure initiatives, over 14,400 protein structures have been solved by Structural Genomics (SG) centers and participating research groups. While the totality of SG data represents a tremendous contribution to genomics and structural biology, reliable functional information for these proteins is generally lacking. Better functional predictions for SG proteins will add substantial value to the structural information already obtained. Our method described herein, Graph Representation of Active Sites for Prediction of Function (GRASP‐Func), predicts quickly and accurately the biochemical function of proteins by representing residues at the predicted local active site as graphs rather than in Cartesian coordinates. We compare the GRASP‐Func method to our previously reported method, Structurally Aligned Local Sites of Activity (SALSA), using the Ribulose Phosphate Binding Barrel (RPBB), 6‐Hairpin Glycosidase (6‐HG), and Concanavalin A‐like Lectins/Glucanase (CAL/G) superfamilies as test cases. In each of the superfamilies, SALSA and the much faster method GRASP‐Func yield similar correct classification of previously characterized proteins, providing a validated benchmark for the new method. In addition, we analyzed SG proteins using our SALSA and GRASP‐Func methods to predict function. Forty‐one SG proteins in the RPBB superfamily, nine SG proteins in the 6‐HG superfamily, and one SG protein in the CAL/G superfamily were successfully classified into one of the functional families in their respective superfamily by both methods. This improved, faster, validated computational method can yield more reliable predictions of function that can be used for a wide variety of applications by the community.


Introduction
A wealth of new protein structures has been reported by structural genomics (SG) initiatives since 2000, but determination of the biochemical function of these structures has proved to be much more difficult than originally envisioned. Reliable methods for prediction of the function of proteins from their three-dimensional (3D) structures constitute a critical current need; such capability will add tremendous value to SG data and advance significantly our understanding of protein function at the atomic level. While structural genomics holds tremendous promise for future applications of great benefit to society, a key step toward the realization of its (still largely untapped) full potential is the ability to determine the function of the thousands of protein structures for which the biochemical function is currently unknown or uncertain.
Current methods for assigning biochemical function are generally informatics based; sequence and structure comparisons are made between the query protein and other proteins in large databases, and functional assignments are transferred based on sequence or structure similarity with previously annotated proteins. Such methods have been described in recent reviews and compilations. [1][2][3][4][5][6][7][8][9] Simple transfer of function based on global sequence or structure similarity can lead to misannotations. 10,11 Automated methods for functional annotation can cause misannotation errors to propagate through databases. Although important efforts are underway to assign correct functions to proteins, 12 there are still thousands of protein structures without functional annotations and many more are misannotated. 13 A local-structure based function prediction method, Structurally Aligned Local Sites of Activity (SALSA), has been described recently. 4,9,14,15 SALSA establishes local spatial arrays of predicted functionally active residues for sets of proteins of known, experimentally determined biochemical function. A distinctive feature of the SALSA approach is that functionally active residues for each protein structure are predicted from computed chemical and electrostatic properties using Partial Order Optimum Likelihood (POOL), [16][17][18] a machine learning method that predicts catalytically important residues using the structure of the query protein as the input. Predicted residues of common type in aligned spatial positions across a set of proteins of known, common function defines a Chemical Signature for that functional type. SALSA then matches the predicted functionally active residues for a protein of unknown function to the Chemical Signatures; a strong match of residue types in aligned spatial positions suggests that function may be transferred reliably.
In this work, a new approach to the local structure matching, Graph Representation of Active Sites for Prediction of Function (GRASP-Func), is introduced; instead of using a Cartesian coordinate representation of the active site residues and relying on global multiple structure alignments as was done previously, 14,15,19 the predicted sets of active residues are expressed in a topological graph representation. This enables much faster alignment and matching of the local active site structures. The Ribulose Phosphate Binding Barrel (RPBB), 6-Hairpin Glycosidase (6-HG), and Concanavalin Alike Lectin/Glucanase (CAL/G) superfamilies are analyzed to illustrate application of the method and to make function predictions for some of the SG proteins predicted to be members of these superfamilies. Each superfamily was chosen for this study because it is medium-sized with functional diversity and with generally good structural coverage and experimental functional characterization within each of the known functional families.
The RPBB superfamily (SCOP 20 ID 51366) has a (b/a)-barrel fold consisting of an eight-stranded parallel b barrel surrounded by eight a helices. 21 RPBB enzymes play essential roles in a variety of different metabolic pathways, including amino acid biosynthesis, pyrimidine biosynthesis, carbon fixation in plants, the nonoxidative phase of the pentose phosphate pathway (which generates ribose 5phosphate, a precursor for the biosynthesis of nucleotides), L-ascorbate metabolism, and the ribulose-monophosphate cycle. Some members of this superfamily also represent potential novel therapeutic targets for antibacterial or antifungal agents. [22][23][24] The 6-HG superfamily (SCOP ID 48208) contains all-a structures sharing a common (a/a) 6 -barrel fold. These enzymes share a similar catalytic mechanism, catalyzing the hydrolysis of glycosidic linkages in poly-or oligo-saccharides. The CAL/G superfamily (SCOP ID 49899) contains all-b proteins sharing a common antiparallel b-strand sandwich core. These enzymes are involved in biosynthesis, cellular development, and localization, and other metabolic processes. Members of both the 6-HG and CAL/G superfamilies have potential applications in biomass degradation and biofuel production. These two superfamilies have previously been analyzed by the SALSA method. 9 In this work, two approaches, SALSA and GRASP-Func, are used to predict the biochemical function of RPBB proteins of unknown function. Additionally, the second approach GRASP-Func is applied to the 6-HG and CAL/G superfamilies. First, the RPBB proteins of known function are used to generate Chemical Signatures for each of the functional families. Then the original SALSA method is applied, with alignments performed by conventional Cartesian-coordinate-based alignment programs on the entire protein structures, from which locally aligned sets of predicted active residues are generated. The 6-HG and CAL/G superfamilies have been sorted previously with SALSA. 9 We then present analysis of the three superfamilies with a new approach, wherein predicted sets of residues are expressed as graphs and local alignments are generated based on the graph representation. This new approach produces locally aligned signatures much faster and allows for more rapid, facile, larger-scale functional classification of protein structures.

Results and Discussion
Chemical signatures based on Cartesian alignment of predicted residues using SALSA The structures of proteins of known function in each superfamily were used to generate the Chemical Signatures for their respective superfamily and were chosen such that sequence homology between any two members within each family is as low as possible (Tables S3-S5, Supporting Information). For most families, at least two experimental structures are available within each family to establish the Chemical Signatures. For families with only one crystal structure available, homology models were generated using protein sequences in these functional families when available (Table S1, Supporting Information). The sequence identity matrix for the previously characterized protein structures in each superfamily was obtained using Clustal Omega 25 and is given in Tables S3, S4, and S5. For each protein, the top 9% of POOL-ranked residues were taken to be the predicted set of functional residues. Since the 6-HG and CAL/G superfamilies have been analyzed previously, 9 only the RPBB superfamily is analyzed by the SALSA method here.
Each superfamily is divided up into its respective functional families. Upon structural alignment of 31 selected RPBB proteins of known function (Table S2, Supporting Information), POOL-predicted residues were found in 24 of the aligned spatial positions and are divided into nine functional families: indole-3-glycerol phosphate synthase (IGPS), tryptophan synthase (TrpA), phosphoribosyl anthranilate isomerase (PRAI), phosphoribosylformimino-5aminoimidazole carboxamide ribotide isomerase (HisA), imidazole glycerol phosphate synthase (HisF), ribulose-phosphate 3-epimerase (RPE), orotidine 5 0 -monophosphate decarboxylase (OMPDC), keto-3-gulonate-phosphate decarboxylase (KGPDC), and hexulose phosphate synthase (HPS). Additionally, the structure of E. coli TrpC (PDB 1pii) in RPBB is bifunctional, where the N-terminal domain (1-255) catalyzes the IGPS reaction and the Cterminal domain (256-452) catalyzes the PRAI reaction. 26 The alignment of the predicted residues for these 31 previously characterized proteins is shown in Table I, in which each row represents a protein structure, with proteins of common biochemical function grouped together. The vertical columns represent spatially aligned positions, obtained from Cartesian-based alignment of the complete structures. POOL-predicted residues are shown in uppercase; aligned residues not predicted are in lowercase. The Chemical Signature residues are highlighted in yellow. Amino acids previously identified as important for catalysis, either from experimental evidence [27][28][29][30][31][32][33][34][35][36][37][38] or by sequence homology with an experimentally characterized protein, 39 are shown in boldface. The normalized SALSA scores for the known members of this superfamily are given in Table S6, Supporting Information. Table I shows that each functional family within RPBB has a unique set of predicted residue types in aligned spatial positions; these local sets of structurally aligned, predicted residues that are common to a particular biochemical function constitute the Chemical Signature for that functional family, with a unique Chemical Signature for each functional family. For example, the Chemical Signature for the IGPS Table I Each row represents a protein structure, with proteins of common function grouped together. The vertical columns represent spatially aligned positions, obtained from Cartesian-based alignment of the complete structures. POOL-predicted residues are shown in uppercase; aligned residues not predicted are in lowercase. Previously reported catalytic residues are shown in boldface. The Chemical Signature residues are shaded in yellow.

Application of SALSA to the SG members of the RPBB superfamily
The SG members of each superfamily were found from searches for proteins with a sequence or keyword match, or structural similarity to previously characterized proteins in each respective superfamily. These SG proteins, with the sources of their structures, are listed in the Table S12, Supporting Information. In the RPBB superfamily, the SG proteins are aligned with previously characterized proteins (Table I), and the aligned, POOL-predicted residues for the SG proteins are scored against the Chemical Signatures for the nine functional families.
The match score MS for SG protein j with the Chemical Signature CS for family k, calculated using scoring matrix M, is obtained as: Normalized match scores S are calculated as: so that a perfect match of aligned residues of the SG protein with those of the Chemical Signature for family k yields a score S of 1. For present purposes, the BLOSUM62 40,41 scoring matrix was used in Eqs.
(1) and (2).  Information) shows the normalized match scores S for 44 SG proteins against the Chemical Signatures for the nine functional families in the RPBB superfamily. For each functional family, the number of aligned positions N in the Chemical Signature is given in the first row. In the next row, for functional families with more than two previously characterized proteins, the range of S values within the set of previously characterized members is given (Table S6, Supporting Information). Table S7 (Supporting Information) reveals that 41 of the 44 SG proteins have high scores with one functional family and substantially lower scores with the other eight functional families. In some instances, a protein exhibiting a strong match with one function and a moderate match with another function (i.e., putative hexulose-6-phosphate synthase SgbH from Vibrio cholerae, PDB 3ieb) may exhibit some promiscuity, as has been observed for previously characterized KGPDC and HPS enzymes. 36,37 The last two proteins shown in Table  S7 (two putative N-acetylmannosamine-6-phosphate 2-epimerases, PDBs 1y0e and 1yxy) have scores below 10.10 with all nine functional families. These two proteins have similar structures to the members of the RPBB superfamily but have predicted function different from those of the RPBB proteins. For one of the superfamily members from Saccharomyces cerevisiae, originally annotated as a HisA/HisF protein (PDB 2agk), its highest score of 10.20 with the HisF family is too low to assign function and therefore it is unlikely to have any of the nine RPBB functions.
The highest match score is used to guide the SALSA functional assignment. Based on the ranges of normalized match scores obtained for the previously characterized proteins, a measure can be derived of the strength of the match to a given functional family. For each SG protein, if the highest normalized match score is greater than or equal to 0.90 or is within the range of scores obtained for the previously characterized proteins in a given functional family, then that highest score is labeled as a strong match (designated s). For normalized match scores less than the strong match threshold but greater than or equal to 0.70, the match strength is labeled moderate (m). Scores between 0.50 and 0.69 are labeled weak matches (w). Scores less than 0.50 are labeled "no match". The top SALSA annotations for each SG protein, labeled (s), (m), or (w), are listed in Table S12, Supporting Information.
Application of SALSA to the SG members of the 6-HG and CAL/G superfamilies Previously, several SG proteins in the 6-HG and CAL/G superfamilies were analyzed using the SALSA method 9 ; additional SG proteins are analyzed here. Aligning and scoring as described above, each SG protein was scored against each functional family in their respective superfamily. Table S9 (Supporting Information) shows the normalized match scores S for 11 SG proteins against the Chemical Signatures for 13 functional families in the 6-HG superfamily. For each functional family, the number of aligned positions N in the Chemical Signature is given in the first row. In the next row, for functional families with more than two previously characterized proteins, the range of S values within the set of previously characterized members is given (Table S8, Supporting Information). Table S9 (Supporting Information) reveals that fewer than half of the SG proteins can be sorted into a functional family reliably. Only uncharacterized protein BT_3781 from Bacteroides thetaiotaomicron (PDB 2p0v), uncharacterized protein BACOVA_03626 from Bacteroides ovatus (PDB 3on6), putative arhamnosidase from B. thetaiotaomicron (PDB 3cih), and putative glycoside hydrolase protein BH0842 from Bacillus halodurans (PDB 2rdy) show strong matches with one functional family (AMAN, AMAN, ALR, and ALF/ALG, respectively). Interestingly, the two SG proteins showing a strong match with the AMAN functional family (PDB 2p0v and 3on6) also show weak matching with the AGG and TRE functional families, suggesting that these two SG proteins might display some promiscuity. In this superfamily, there are a few SG proteins that show weak matching with one functional family; putative alkaline invertase from Nostoc sp. (PDB 5goo) with AGG, two putative GH105 family proteins from Klebsiella pneumoniae (PDB 3pmm) and Salmonella paratyphi (PDB 3qwt) with UGH, and two putative N-acetylglucosamine 2-epimerases from Salmonella typhimurium (PDB 2afa) and Xylella fastidiosa (PDB 3gt5) with NAE. Two SG proteins, lin0763 protein from Listeria innocua (PDB 3k7x) and putative glycosyl hydrolase from B. thetaiotaomicron (PDB 4mu9) do not show significant normalized scores with any of the functional families. The top SALSA annotations for each SG protein, labeled (s), (m), or (w), are listed in Table S12, Supporting Information.
For the CAL/G superfamily, Table S11 (Supporting Information) shows the normalized match scores S for eight SG proteins against the Chemical Signatures for the six CAL/G functional families. Similar to Table S9 (Supporting Information), the number of aligned positions N in the Chemical Signature is given in the first row, followed by the range of S values within the set of previously characterized members (Table S10, Supporting Information). Table S11 (Supporting Information) reveals that one protein, putative GH16 family protein from Mycobacterium smegmatis (PDB 3rq0), has a score of 10.40. Normally, this would be considered "no match" according to our criteria; however, since the range of scores between the previously characterized members of the family is low (0.60-0.72) due to their different substrate specificities, we have assigned a weak functional annotation to this SG protein. Table S12 (Supporting Information) lists the SALSA results and shows that the other seven SG proteins have no match with any functional family we have analyzed. These SG proteins may be in functional families that lack structural coverage or are novel functional families.
Function prediction with a graph theory approach (GRASP-Func) Here we introduce a computationally faster approach to sorting superfamilies according to biochemical function. For each protein structure in each superfamily, the set of highly-ranked POOL residues is represented as a set of points in 3D space to form a graph representation, generated by Delaunay triangulation, of the active site. These graph representations can match rapidly one active site to another. The topological graph descriptors represent each predicted residue as a single point in space, using the coordinates of the a carbon atoms. This generates a set of tetrahedra, where the residues are represented by the vertices and the edges indicate that the two joined residues are neighbors. Delaunay triangulation has been used previously for protein structural alignment by common volume superposition 42 ; here it is applied to identify similar spatially localized regions of structures.
The sets of tetrahedra that contain POOLpredicted residues for a pair of proteins are then compared using a pairwise matching algorithm, described in the Methods section. Sets of proteins with matched tetrahedra are then grouped together by this algorithm. Matches between sets of proteins of known function with a query protein of unknown function thus enable function prediction for the query protein. One of the main advantages of GRASP-Func over SALSA is that GRASP-Func does not rely on global structural alignments, which can be very time consuming and labor intensive. Additionally, when analyzing function similarity across folds, SALSA requires a manual alignment process 4 while GRASP-Func can analyze function without the need for global alignments. While SALSA makes function predictions using a table of spatially aligned, functionally important residues for protein structures within a superfamily (as illustrated in Table I), GRASP-Func uses similarity between sets of four-membered graphs and generates a figure showing the proteins of similar function grouped together; individual proteins are represented as nodes and the thickness of each edge shows the degree of similarity between the two connected proteins (as illustrated in Figs. 1-3). GRASP-Func was optimized with the RPBB superfamily; 6-HG and CAL/G superfamilies were then used to test the method.
In the RPBB superfamily, the previously characterized proteins listed in Table S2 (Supporting Information) are sorted correctly into nine groups by GRASP-Func (Fig. S3, Supporting Information). This correct classification into nine functional families is the same as the SALSA classification shown in Table I. In the 6-HG superfamily, the previously characterized proteins are sorted into 13 groups by GRASP-Func (Fig. S4, Supporting Information). This functional classification is similar to the SALSA classification, with the exception of the Phosphorylase II family (Group 12). The maltose phosphorylase from Lactobacillus brevis (PDB 1h54) and the nigerose phosphorylase from Clostridium phytofermentans (homology model NGP1) do not show a correlation using this method. This is attributed to the homology model generated for nigerose phosphorylase, which was built from the maltose phosphorylase crystal structure (PDB 1h54) template but has a low model quality score 9 (Table S1, Supporting Information). The model structure was analyzed by PROCHECK, 43 and the results showed only 88.2% of the nonglycine/proline residues (605 residues) are in the most favored regions, 10.1% (69 residues) in additionally allowed regions, 1.2% (8 residues) in generously allowed regions, and 0.6% (4 residues) in disallowed regions. A good quality model is expected to show 90% or more of the nonglycine/proline residues in favored regions. The residues in the generously and disallowed regions are located distal from the active site and may disrupt the network within the protein structure. Similarly, the 19 previously characterized proteins in the CAL/G superfamily are sorted into six biochemical functional groups by GRASP-Func (Fig. S5, Supporting Information), with the same classification as that of SALSA. The GH family 16 functional family (Group 4) shows some separation due to the different substrate specificities of the proteins of known function.
The 6-HG superfamily proteins were sorted by GRASP-Func (Fig. 2), and the results show that for seven of the 11 SG proteins, GRASP-Func is able to assign the same function as SALSA (Table S12, Supporting Information). The two putative GH105 family proteins from K. pneumoniae (PDB 3pmm, H4) and S. paratyphi (PDB 3qwt, H5) are assigned a weak (10.51) UGH function by SALSA but are assigned a URH function by GRASP-Func. Both families function by hydrolyzing their respective substrates and have a number of similar residues in their active sites. 9 However, SALSA can only obtain a reliable Chemical Signature if the family has two or more protein structures and/or sequences of known function. In this case, the URH functional family has only one known representative. It is possible that SALSA assigned UGH function over URH function because a reliable Chemical Signature for URH is unavailable. In contrast, GRASP-Func does not rely on the Chemical Signatures and global structural alignments and is able to provide functional annotations with only one known representative.
Putative a-L-fucosidase from Bacillus halodurans (PDB 2rdy, H7 in Fig. 2) is predicted to be in the ALF/ALG functional family. Upon further analysis with individual members of the functional family, SALSA predicts galactosidase function. In GRASP-Func, there is a strong match between this SG protein and the galactosidase function, as illustrated in Figure 2 by the darker edge connecting it to a-L-galactosidase from Bacteroides ovatus (PDB 4ufc, 7a in Fig. 2). Two SG proteins, putative GH76 family protein from Listeria innocua serovar 6a (PDB 3k7x, H10) and putative glycosylhydrolase from Bacteroides thetaiotaomicron (PDB 4mu9, H11) are unable to be annotated by either method. It is possible they are members of new functional families.
The CAL/G superfamily proteins were also sorted by GRASP-Func (Fig. 3). In this instance, only one SG protein, putative GH family 16 from Mycobacterium smegmatis (PDB 3rq0, C1 in Fig. 3) is able to be assigned function by both SALSA and GRASP-Func, in this case as having GH family 16 function (Table S12, Supporting Information). Specifically, Figure 3 shows that this protein likely has endo-b-1,3-glucanase activity. While neither SALSA nor GRASP-Func can assign function to the other seven SG proteins, GRASP-Func shows that the three putative b-xylosidase (PDBs 1y7b, 1yif, and 1yrz, C224 in Fig. 3, respectively) cluster together away from the other families and have a strong connection to each other as shown by the thick edges. Similarly, the two putative sugar hydrolases (PDBs 3h3l and 3nmb, C5 and C7 in Fig. 3, respectively) and the two putative glycosyl hydrolases (PDBs 3hbk and 3osd, C6 and C8 in Fig. 3, respectively) form a four-membered, well-connected cluster. These two clusters could represent new functional families in the superfamily.
The amount of time it takes to sort a set of proteins with GRASP-Func varies, depending on the degree of similarity between pairs; sets with higher variability discard larger numbers of pairs early and therefore the sorting proceeds faster. In a typical run on an Intel Xeon E3-1220 v3 CPU running at 3.10 GHz, with 16 GB of RAM, it took 15 min of clock time to obtain 2240 results. This is at least several orders of magnitude faster than the full structural alignment employed in the original SALSA method, which can take hours to run depending on the size of the superfamily being analyzed. In addition, SALSA often requires manual adjustments, or unification of multiple, smaller alignments, to obtain the best local alignments, particularly for large sets of structures. GRASP-Func also enables matching of functional types across folds; while this is possible in the original SALSA method, 9 it is slow and labor intensive because manual alignments are required. The thickness of each edge shows the degree of similarity between the two connected proteins. PDB IDs for proteins of known function: 1m4w, 1h4g, 1bcx (1a-c); 1uu4, 1h8v, 2nlr (2a-c); 1z3t, 1dy4, 2rfw (3a-c); 2ayh, 1dyp, 3ilf, 2vy0, 1mve (4a-e); 1uai, 1j1t, 1vav (5a-c); 2fir, 1y43 (6a-b). Each SG protein is numbered based on its label in Table S12, Supporting Information.
SALSA and GRASP-Func both incorporate computed chemical properties from the POOL method to predict protein function from 3D structure. Both methods are based on structure similarity at the local site of biochemical activity and both have successfully sorted members of the three superfamilies into families according to predicted biochemical function. The graph representations of GRASP-Func obviate global Cartesian alignments and therefore yield local-structure-based function assignments substantially faster and can be fully automated. Faster protein function annotation methods like GRASP-Func will help correct function misannotations in databases and provide the scientific community with correct information. This will add a substantial amount of information to the already extensive amount of work done through SG efforts.

POOL predictions
POOL predictions were made as described by Somarowthu et al. 18 SALSA predictions based on Cartesian alignments SALSA predictions were made as described by Wang et al. 15 The top 9% of the residues in the POOL rankings were taken to be the predicted, functionally active residues that are marked in the structural alignments. When more than half of the proteins in a functional family have POOL-predicted residues of common type in an aligned position, that residue becomes part of the Chemical Signature.

GRASP-Func Analysis
The protein structures were preprocessed to convert the coordinates into a set of tetrahedra and to identify the tetrahedra near the active site, based on the POOL rankings. To achieve this, first Delaunay triangulation was performed on the protein structure using Qhull. 44 The vicinity of the active site is determined by the top 50 residues in the POOL rankings. All tetrahedra that contain a POOL-predicted residue, or have a vertex connected to a POOLpredicted residue, are collected for matching analysis. In a pair of proteins P 1 and P 2 , the tetrahedra in the active site vicinity that have been identified in the preprocessing step are compared and seed pairs are sought. Seed pairs are ranked using POOL rank, residue similarity as measured by the BLO-SUM62 40,41 matrix, and lengths of the edges of the tetrahedra. If tetrahedron t j,1 in protein P 1 and tetrahedron t k,2 in protein P 2 have residues with high POOL rankings and chemical similarity, then the pair t j,1 and t k,2 is a seed pair. Then seed pairs of tetrahedra are compared according to the edge lengths, that is the distances between a carbon atoms. Additional features of a tetrahedron used in the matching algorithm are the volume, the sum of the lengths of the edges, and the relative orientation. The average volume for a tetrahedron in the RPBB superfamily is 14.4 Å 3 , so pairs of tetrahedra with a volume difference greater than 14.4 Å 3 are rejected. The average sum of edge lengths is 9.6 Å , so pairs are rejected if total edge length difference exceeds 9.6 Å . Then the vertices, which represent the individual amino acids, are analyzed further. With the set of surviving pairs, the vertex pairs v j,m,1 in t j,1 from P 1 and v k,n,2 in t k,2 from P 2 , where m and n are indices for the individual vertices in the tetrahedron, are further filtered according to the following sequential steps: