Here, a scalable, accurate, reliable, and robust protein functional site comparison algorithm is presented. The key components of the algorithm consist of a reduced representation of the protein structure and a sequence order-independent profile-profile alignment (SOIPPA). We show that SOIPPA is able to detect distant evolutionary relationships in cases where both a global sequence and structure relationship remains obscure. Results suggest evolutionary relationships across several previously evolutionary distinct protein structure superfamilies. SOIPPA, along with an increased coverage of protein fold space afforded by the structural genomics initiative, can be used to further test the notion that fold space is continuous rather than discrete. functional site ͉ structure T he evolutionary relationship between protein sequences, protein structures, and their associated function(s) remains a central topic of molecular biology and one resulting in the development of many computational methods (1-3). A central question is: What were the early protein folds and how did these folds change over long evolutionary time scales (4-7)? Comparative genomics studies and structural and phylogenetic analyses (8-10) have established that a subset of proteins, dominated by the structure classification of proteins (SCOP) (11) ␣/␤ class, were likely present in the last universal common ancestor (12, 13). Concurrently, growing evidence suggests that recurring substructures, that is, 3D fragments of noncontiguous sequence shared between different folds, may be clues that protein fold space is more continuous than discreet (14, 15) . The sequence/ structure similarity of such substructures correlates well with the similarity of function found between the different folds containing these substructures (16). The notion that protein fold space is a continuum is further supported by recent studies that show that protein domains can adopt different topologies through combination, swapping, deletion (4, 17, 18) , and cyclic permutation (19, 20) of subdomains. Likewise, new folds can emerge from accretion (21) or embellishment (22) of substructures around a core of conserved secondary structures. Given these findings concerning the dynamic nature of protein structure and the possible continuous nature of protein fold space, it is important to distinguish proteins that share a common ancestor (divergent evolution) from those that have adopted common structural constraints (convergent evolution). Typically, evolutionary relationships between protein sequence, structure, and function are deduced from the respective comparisons among known genes and their products. These comparisons are made at various levels, from genome sequences to protein domains and motifs to biochemical pathways. Such comparisons may miss important relationships because sequence relationships may be too weak to detect, and/or fail to identify complex evolutionary events such as domain swapping and cyclic permutation. Likewise, differences in global protein structure may disguise a true evolutionary relationship that exists between substructures. One approach, which involves the comparative analysis of substructures, including functional sites between proteins (1, 23-26), has been successful in detecting evolutionary relationships between different fold superfamilies and has been applied mostly to enzyme families. One study of 31 diverse enzyme superfamilies revealed that functional diversity during evolution is achieved by local sequence variation and domain shuffling (24). Such functional diversity can also be observed within a single SCOP superfamily. For example, within the protein kinase-like superfamily, it has been suggested that atypical kinases diverged early in evolution from protein kinases (26). In doing so, the overall catalytic mechanism is retained through a high level of conservation associated with the ATP binding cassette, thus preserving phosphorylation, yet the substrate binding motif exhibits significant diversity. In the case of mechanistically diverse enzymes, whose members catalyze different overall reactions but share a partial reaction, it has been found that these enzymes use a similar active site to generate a common intermediate, then direct the intermediate to different products in different active sites (25). Beyond these case studies, the global evolutionary relationship of functional sites across fold space has not been systematically studied and remains elusive. Global functional site comparison has been thwarted by the lack of efficient and accurate computational tools to undertake such a large scale comparison and a lack of rigorous statistics to test their similarity. The work described herein is a step toward accurate and efficient functional site comparison and analysis and is subsequently applied to seek out new evolutionary relationships. Although the concept of functional site matching is not new, and a variety of approaches have been attempted (27-47), it has not proven an easy task to design and implement a practical software solution with performance that is close to that of routine sequence comparison. These site comparison algorithms usually consist of three interrelated components; the representation of the functional site, an algorithm to superimpose two sites and a method to score their similarity. The functional site is usually represented either by a coordinate set with certain physicochemical or evolutionary properties, or by 3D shape descriptors that define pockets within the protein (44). The coordinate set can consist of atoms (28) We use this algorithm here to confirm or propose hitherto unseen evolutionary relationships and note that a larger scale study is underway. The algorithm proposed is different from previous algorithms in several respects. First, our algorithm searches for similar functional sites by scanning whole proteins in the spirit of local sequence alignments, i.e., it is not necessary to predefine functional sites during the search. Second, most of the available algorithms score similarities between whole functional sites. However, the complete ATP and NAD binding sites may look different overall, but they may share a similar subpocket for binding chemically and conformationally similar fragments, adenine in this case. This subsite similarity would be missed by most current algorithms. Our algorithm is designed to detect similar subsites across all of fold space. Third, we proposed a sequence order-independent profile-profile alignment (SOIPPA) algorithm, which is more general and sensitive than others (29, 38, 48) , because the sequence order of amino acid residues in the functional site is not necessarily conserved. Fourth, with the aim of applying our approach to homology models and low resolution structures, the algorithm was purposely designed and tested with a C␣ only representation of the protein structure. Finally, the statistical significance of the similarity is evaluated by using a nonparametric statistical method based on a fold distribution model. The algorithm is evaluated on a benchmark set and a control group that includes 247 and 101 nonredundant protein chains of diverse folds with and without adenine binding pockets, respectively. Our results show that SOIPPA is suitable for application on a large-scale. For example, in studies reported elsewhere, we have successfully applied SOIPPA to identify off-targets for pharmaceuticals for which the primary drug-receptor complex is present in the Protein Data Bank (PDB). Searching by using the primary site information, we have identified off-targets for selective estrogen receptor modulators (49) and have been able to repurpose an existing drug for treatment of drug resistant tuberculosis (S. Kinnings, L.X., and P.E.B., unpublished data). In both cases the primary and secondary targets are in different gene families. Here, we focus on a few examples that support the notion that fold and functional space may be continuous rather than discrete. Comparison of Functional Site Superposition Algorithms. Two datasets were used in the comparison study. The 247-benchmark consisted of 247 nonredundant protein chains known to bind an adenine containing ligand. The 101-control consisted of 101 nonredundant protein chains believed not to bind an adenine containing moiety. From the 247-benchmark, 30,381 benchmark pairs were generated, and, from the 101-control, 24,947 control pairs were generated in an all-by-all pairwise comparison (see Methods). All benchmark and control pairs were used to evaluate the performance of the sequence order-independent profileprofile alignment (SOIPPA) algorithm [see Methods and supporting information (SI) Text]. We further implemented several existing functional site comparison algorithms that employ maximum cliques (28, 30, 40) and similarity measures, using substitution matrices (29, 38), physiochemical properties, and amino acid grouping (40, 43) . Using the benchmark and control data, their performance was compared with SOIPPA with respect to alignment quality and sensitivity/specificity during site searches. We also include results from PSI-BLAST (50) and CE (51) in our comparisons to relate functional site similarity to more traditional sequence and structural similarity, respectively. To measure the alignment quality, we use the root mean square deviation (RMSD) between ligands bound to the aligned functional sites. In our 247-benchmark, all proteins bind to ligands containing a chemically identical and conformationally rigid adenine plus a chemical and/or conformational variable component. Although the complete functional sites may be different as a result of binding to ligands with different compositions or conformations, they may share similar subsites that bind to similar ligand fragments. Our purpose is to detect such similar subsites. More specifically, the RMSD between common molecular moieties from diverse ligands is used to measure the alignment quality of the functional site comparison algorithm. SOIPPA provides the best performance when compared with other methods. As shown in Fig. 1 , using SOIPPA, 6.5% and 25.9% of pairs are aligned with RMSD values of Ͻ2.0 and Ͻ5.0 Å, respectively. Lists of these pairs are given in Dataset S1 and Dataset S2. Using a BLOSUM45 substitution matrix (52), the percentage of aligned pairs drops to 5.3% and 22.3%, respectively. Using a chemical similarity-based method (53), only 4.2% and 18.1% of pairs are aligned with RMSD values of Ͻ2.0 and Ͻ5.0 Å, respectively. Using a maximum size clique method based on amino acid grouping (40), only 3.1% and 15.7% pairs are aligned with RMSD values of Ͻ2.0 and Ͻ5.0 Å, respectively. Moreover, the frequency of aligned ligands when performing global structure similarity, using CE, is significantly lower than that achieved with SOIPPA. Specifically, for pairs of proteins where the RMSD values of aligned ligands are Ͻ2.0 and Ͻ5.0 Å, the percentage of dissimilar structures (defined by a CE Z score Ͻ3.5) are Ϸ40% and Ϸ60%, respectively (Fig. S1 ). These results indicate that the functional site is more structurally conserved than the global protein structure. Although similar functional sites may arise from convergent evolution (69), the high occurrence among dissimilar structures raises the possibility of very significant divergent evolution. It is noteworthy that SOIPPA uses the whole-protein structure when performing the ligand site search and comparison and does not require that the ligand or functional site be known a priori. SOIPPA detects local similarity between two proteins by aligning two global structures in a manner that is independent of sequence order, important because sequence order is not always conserved (4) (Fig. S2, Fig.  S3, Fig. S4, and Fig. S5 ). Beyond identifying similar subsites that bind to common molecular fragments with rigid conformations from different ligands, SOIPPA also detects similar subsites that bind flexible molecular fragments that adopt similar conformations ( Fig. S3 and Fig. S5 ). To evaluate the performance of SOIPPA in detecting evolutionary relationships, we first see whether the algorithm can identify known sequence and structural homologous within the same SCOP superfamily. Among the 247-benchmark pairs, Ͻ5% of them (1,230 pairs) are from the same SCOP superfamily. If only these 1,230 pairs are considered as true positives, Fig. 2a illustrates the performance of PSI-BLAST (50), CE (51), and SOIPPA in detecting remote homologous that belong to the same SCOP superfamily. For a false-positive ratio of 0.05, the coverage of PSI-BLAST, CE, and SOIPPA is 0.55, 0.60, and 0.75, respectively. If SCOP superfamilies are taken as the gold standard for defining remote evolutionary relationships, these results illustrate the well known fact that structure is more conserved than sequence. However, global structure comparison falls significantly short of SOIPPA, which takes both evolutionary profiles and structural constraints within the functional site into account. Consequently, it is more sensitive in detecting remote evolutional relationships than either PSI-BLAST or CE. The question then becomes, can SOIPPA detect functional similarities missed by SCOP, that is, relationships across superfamilies? In the 247-benchmark, there are 15,058 pairs aligned from different known SCOP superfamilies, and Ϸ30% of them are identified by SOIPPA with a false-positive ratio of 0.05 (Fig. 2b) . Fig. 2b also shows that the sequence and structural similarity of these cross-superfamily pairs is not significant. To illustrate that the sensitivity and specificity of functional site comparison, using SOIPPA, is improved over other methods, we define all of the 30,381 247-benchmark pairs to be true positives (Fig. 3) . The false-positive rate for SOIPPA is Ϸ30% lower than that for a weighted method using an amino acid substitution matrix or chemical similarity at a 15% true-positive ratio. In the low false-positive region, the weighted methods (SOIPPA, substitution matrix, and chemical similarity) perform better than nonweighted amino acid grouping methods. However, if high false positives are tolerated, amino acid grouping performs best. As a baseline, the performance of CE is also included in the comparison. Clearly, a global structural comparison algorithm, such as CE, cannot detect the similarity between proteins when the overall structures have changed. For the 5,917 structurally related pairs (Ϸ20% of 247-benchmark pairs) that have the same CATH topology, the coverage by CE is lower than SOIPPA at a false-positive ratio Ͻ0.05 (Fig. S6a) . Moreover, the frequency of well aligned ligands from global alignment by CE is also lower than from SOIPPA (Fig. S7 ). For those cross-CATH topology pairs, it is not surprising that the coverage by CE is reduced to Ͻ10%. However, the coverage of SOIPPA is only slightly lower than that of the same CATH topology pairs (Fig. S6b) . These results suggest that SOIPPA is not sensitive to overall structural changes, implying that local functional sites are more conserved than global sequence or structure. This then leads us to the question: Can evolutionary relationships be detected that are not defined by SCOP? Evolutionary Linkage Across Fold and Functional Space. Using SOIPPA to scan a nonredundant set of protein structures from the complete PDB (7,644 chains including both apo and holo structures) against several functional sites, we have detected statistically significant similarities between fold families whose global structures are very different yet whose local functional sites appear invariant over evolutionary time scales. Due to lack of structural and/or functional annotations on the whole set of structures, the true and false positives cannot be defined conventionally. However, the false-positive rate of the predication can be approximated by the P value from the fold distribution score (see SI Text). One example is the relationship identified between the protein kinase-like, the SAICAR-synthesase-like, and the ATP grasp superfamilies after scanning the 7,644 proteins against ATP (PDB ID code 1ODB, SAICAR synthase-like) and ADP (PDB ID code 2HGS, ATP-grasp) binding sites. These three superfamilies show significant binding site similarities with each other (Fig. S8 and Fig. S9 ). Their evolutionary relationship has been proposed by others through manual analysis of structure (26, 54), even though their global structural similarities cannot be detected by structural comparison methods. For example, the CE similarity Z score and RMSD are 2.3 and 5.17 Å between PDB structure 1WBP (protein kinase-like) and 2HGS (ATP-grasp), respectively. Given that such remote evolutionary relationships  can be established automatically by SOIPPA with high statistical significance we have a tool with the potential to discover unseen evolutionary relationships throughout protein fold space. We explore this notion further with some specific findings. By searching against the 7,644 proteins, using the ATP binding site of PEP carboxykinase (PCK) (PDB ID code 1AYL), besides two top ranked PCK structures, a protein annotated as a P loop nucleotide triphosphate hydrolase (NTH) is ranked as the second most significant cross-superfamily hit (PDB ID code 1LS1; P ϭ 8.4 ϫ 10 Ϫ5 ). Although there is no structural annotation for the most significant cross-superfamily hit (PDB ID code 1YRB; P ϭ 8.4 ϫ 10 Ϫ5 ), its structure is similar to PDB 1LS1 (CE z score, 4.9; RMSD, 2.99 Å). Other high ranked hits with P Ͻ 2.0 ϫ 10 Ϫ3 are also similar to the P loop NTH fold (Table  S1 and Table S2 ). The functional site alignment revealed by SOIPPA further suggests that these two superfamilies are evolutionarily related even though there is no observable global structure similarity. Besides the conserved Walker motif A (55) (GXXXXGKT/S) that is the common site for binding nucleotides in many proteins, several residues important for ATP binding and catalytic activity are also conserved. The N-terminal histidine residue (His-232) in PCK contacts the ␥-phosphate of ATP and contributes to either activation of the substrate or stabilization of the transition state (56, 57). H232 of PCK can be superimposed onto H248 of NTPase, although the sequence order is switched to the C terminus. Likewise, residues D268 and D269 in PCK play roles in stabilizing a Mg 2ϩ -Mn 2ϩ bimetal cluster at the active site (56, 58). An evolutionary linkage tree and multiple functional site alignment (Fig. S4 and Fig. S10 ) indicate that these two residues are strictly conserved in several other enzymes containing a P loop NTH fold. It is noteworthy that the overall conformation of the ATP molecule bound to these two classes of proteins is different, with an almost 180°flip in the adenine fragment. However, the conformations of the triphosphate tails are quite similar (Fig. S5) . Our predication based on only a subsite similarity is consistent with previous speculation based on detailed structure analysis (59). However, both SCOP (11) and CATH (60) classify them differently missing an apparent evolutionary relationship. The NAD-binding Rossmann fold is one of the most common protein folds and observed in a large number of enzyme families (61) and believed to have evolved early (8). Evidence suggests that significant structural changes were accumulated from the original Rossmann fold (11, 61). By scanning structures against a NAD-binding site (PDB ID code 2C5A), SOIPPA detected similar functional sites from different SCOP superfamilies that have not been observed before. Except for the nicotinamide adenine dinucleotide (FAD) binding site from the FAD/NAD binding fold (P ϭ 5.8 ϫ 10 Ϫ5 ), all other hits with a P Ͻ 2.0 ϫ 10 Ϫ3 contain the Rossmann fold, but with versatile domain insertions and/or sequence permutations (Table 1) . Although their cofactor binding sites may bind different ligands, they share conserved sequence motifs (Fig. 4) . Structural superimposition of the ligands indicates that the common adenine moiety from these cofactors is well aligned (Fig. S11, Fig. S12 , and Fig. S13 ). Given that these proteins all adopt the same Rossmann fold topology, divergent evolution would seem the most likely scenario (62). Among the hits (Table 1) , structures with SIS and glycerate kinase I domains do not have cocrystallized ligands. Their aligned residues are both located in deep pockets, indicating that they are potential ligand binding sites (Fig. S14) . Moreover, their functions suggest that they can bind to nucleotides. For example, glycerate kinase catalyses the phosphorylation of (R)-glycerate to 3-phospho-(R)-glycerate in the presence of ATP (63). Functional Sites as a Linkage to Trace Distant Evolutionary Relationships Across Fold Space. It has been found that convergent evolution to achieve similar enzyme active sites is a common event (64). However, divergent evolution retaining only a ligand binding site cannot be discounted in some case given that the percentage of similar ligand binding motifs is high among dissimilar structures. Unlike enzyme active sites that typically only involve three or less spatially distinct and noncontiguous residues, many of the similar ligand binding motifs found here contain a number of conserved consecutive residues. We hypothesize that it is more efficient for nature to reuse such local features to bind the same or similar molecular moiety, but to modify other structural components to achieve different functions than to reinvent the whole structure from scratch. This is particularly true for adenine-binding folds studied here. They are found throughout fold space and have been predicated among the most ancient folds by several studies (12, 13, (65) (66) (67) (68) . It is possible that a small molecule containing an adenine fragment was the first ligand recognized by a protein (69). Gross topology of protein can change dramatically in cases already attributed to divergent evolution, using mechanisms of structure drift, segment swapping, the insertion of additional structures, and fusion or permutation by duplication events (4, 70, 71) . Such events imply that the sequence order of conserved motifs may change. Thus, it is difficult to detect evolutionary relationships between different fold and functional families with overall sequence or structural similarity, although several studies have revealed that evolutionary relationships can be established within superfolds, using profile-profile comparison (72), or between different folds with similar short sequence motifs (73). Here, we assume that the functional site carries the evolutionary fingerprint, because it is one of the most important parts of the structure that must be conserved. Indeed, our results suggest that functional site residues and their spatial arrangements are conserved even when the global structure has changed. As such, the functional site can be used as a common linkage from which to potentially trace the evolutionary transition of the protein topology and discover underlying genetic events (15) and physical rules (17, 74) that govern structural changes. A number of in silico and in vitro evolutionary studies have suggested that the sequence or profile similarity between proteins is a strong indication of divergent evolution (5, 6, 72, 73, (75) (76) (77) (78) . Moreover, it is necessary to establish the structural relationships based on a rigorous statistical test (79). By combining sensitive profile-profile alignment, structural constraints independent on the sequence order, and a rigorous statistics test, it is possible for us to propose divergent evolutionary relationships across fold space. Our results are consistent with existing observations, namely that convergent evolution of domain architectures is rare (62). However, to date, only a small percentage of such divergent structural relationships have been identified across fold space. Beyond new methods, structural genomics, because it strives to fill in gaps in the coverage of protein fold space (80), will also provide new clues as to the evolution and continuity of this space. A new algorithm is introduced for functional site comparison that is based on a reduced protein structure representation and a sequence order independent profile-profile alignment. Using a well defined adenine-binding pocket from various folds as a benchmark, it has been shown that the proposed algorithm outperforms both nonweighted and sequence-weighted methods. Although there is still room for significant improvement in the performance of the current implementation, it provides a framework for developing a robust, reliable, accurate, and scalable functional site comparison algorithm. We show that the algorithm has the ability to detect new evolutionary relationships across existing current discreet descriptions of fold and functional space opening the door to tracing fold changes during evolution. Benchmark and Control Data. From the Research Collaboratory for Structural Bioinformatics Protein Data Bank (81), a set of 247 protein monomer chains, which are bound to the following ligands, ATP, ADP, NAD, FAD, Sadenosylmethionine (SAM), and SAH, were selected to use as a benchmark (247-benchmark). All of these ligands include adenine as a common molecular fragment. Structures with multiple heteromeric chains were not considered. The sequence identity between any pair of chains was Ͻ30%. These chains cover 106 SCOP superfamilies and 152 enzyme classifications (EC). In addition, 81 and 70 chains were yet to be given SCOP and EC assignments. As a control, a set of proteins not bound to a ligand containing ribose, adenine, flavin, and nicotinamide were extracted from the PDB. Subsequently, redundant chains were removed by using a 30% sequence identity cutoff against each other and with the 247-benchmark. The final control contained 101 protein chains (101-control). The chains from the 247-benchmark were aligned against each other in a pairwise fashion to generate 247 ϫ (247 Ϫ 1)/2 ϭ 30,381 benchmark pairs. Among them, Ϸ5% of the pairs shared the same SCOP superfamily. The 247 ϫ 101 ϭ 24,847 control pairs were obtained by aligning the 247-benchmark against the 101 chains in the 101-control set. Sequence Order Independent Profile-Profile Alignment of Functional Sites. Protein structures are represented by Delaunay tessellation of C␣ atoms and characterized with geometric potentials as described fully in ref. 82. Each C␣ atom is assigned a probability distribution and position specific score matrix of 20 aa. The regular tessellation of the protein structure can be considered a graph representation in 3D space. The protein is scanned and aligned to the functional site by finding the maximum-weight common subgraph between two encoded protein graphs (83). Full details of the algorithm are provided in the SI Text. Performance Evaluation. The alignment quality between the functional sites of two protein structures is evaluated by the RMSD between the common molecular moieties associated with the ligands. In addition, search performance is evaluated by a true and false-positive rate defined as follows: True positive rate ϭ true positives/(true positives ϩ false negatives), False positive rate ϭ false positives/(false positives ϩ true negatives). We sincerely thank the anonymous reviewers and the editor for their constructive suggestions in improving the manuscript.
SOIP-PA represents functional sites by using a delaunay tesselation of C -atoms REF .
Detecting evolutionary relationships across existing fold space, using sequence order-independent profile-profile alignments
Proceedings of the National Academy of Sciences of the United States of America
Proceedings of the National Academy of Sciences of the United States of America
Biology Medicine