article: A new monotonicity-constrained maximum likelihood approach, called Partial Order Optimum Likelihood (POOL), is presented and applied to the problem of functional site prediction in protein 3D structures, an important current challenge in genomics. The input consists of electrostatic and geometric properties derived from the 3D structure of the query protein alone. Sequence-based conservation information, where available, may also be incorporated. Electrostatics features from THEMATICS are combined with multidimensional isotonic regression to form maximum likelihood estimates of probabilities that specific residues belong to an active site. This allows likelihood ranking of all ionizable residues in a given protein based on THEMATICS features. The corresponding ROC curves and statistical significance tests demonstrate that this method outperforms prior THEMATICS-based methods, which in turn have been shown previously to outperform other 3D-structure-based methods for identifying active site residues. Then it is shown that the addition of one simple geometric property, the size rank of the cleft in which a given residue is contained, yields improved performance. Extension of the method to include predictions of non-ionizable residues is achieved through the introduction of environment variables. This extension results in even better performance than THEMATICS alone and constitutes to date the best functional site predictor based on 3D structure only, achieving nearly the same level of performance as methods that use both 3D structure and sequence alignment data. Finally, the method also easily incorporates such sequence alignment data, and when this information is included, the resulting method is shown to outperform the best current methods using any combination of sequence alignments and 3D structures. Included is an analysis demonstrating that when THEMATICS features, cleft size rank, and alignment-based conservation scores are used individually or in combination THEMATICS features represent the single most important component of such classifiers.
  In prior implementations of THEMATICS for the identification of active-site residues from the 3D structure of the query protein [3]–[5], titration curve shapes were described by the moments of their first derivative functions. These first derivative functions are essentially probability density functions and give unity when integrated over all space. In Ko et al. [3], the third and fourth central moments μ3 and μ4 of these probability functions were used. These moments measure asymmetry (skewness) and, roughly, the area under the tails relative to the area near the mean (kurtosis), respectively. In Tong et al. [5], the first moment and second central moment were also used. In each of these approaches, the moments measure deviations from normal curve shape and the analyses identify the outliers, the residues that deviate most from the normal proton occupation behavior. These prior approaches all use spatial clustering, so that outlier residues are reported as positive by the method if and only if they are in sufficiently close spatial proximity to at least one other outlier. Thus the previous THEMATICS identifications involve two stages, where the first stage makes a binary (outlier / not an outlier) decision on each residue and the second stage finds spatial clusters of the outliers. In the new approach reported here, every residue is assigned a probability that it is an active-site residue. Here, as an alternative to the clustering approach, we introduce features that describe a residue' s neighbors; we call these environment features. For a given scalar feature x, we define the value of the environment feature xenv (r) for a given residue r to be: (1) where r′ is an ionizable residue whose distance d (r′, r) to residue r is less than 9Å, and the weight w (r′) is given by 1/d (r′, r) 2. In this study, we use the same features μ3 and μ4 used in the Ko [3] approach, along with the additional features μ3env and μ4env. Thus every ionizable residue in any protein structure is assigned the 4-dimensional feature vector (μ3, μ4, μ3env, μ4env). The present approach has a number of advantages. Specifically, active residues may be selected in one step and they can be rank-ordered according to the probability of involvement in an active site. Furthermore, while THEMATICS previously has been applied to ionizable residues only, the present approach opens the door to direct prediction of non-ionizable active site residues, because the environment features μ3env and μ4env are well defined for all residues, including the non-ionizable ones. Finally, additional geometric features that are obtainable from the 3D structure only may be readily combined with the four THEMATICS features in order to enhance performance. Geometric features, such as the relative sizes of the clefts on the surface of the protein structure, have been shown to correlate with active site location [6], [7]. For instance, for the majority of single-chain proteins, the catalytic residues are in the largest cleft. However geometric features alone do not perform comparatively well for active residue prediction, particularly because they are not very selective. It is shown here that cleft size information combined with THEMATICS electrostatic features yields high performance in purely 3D structure based functional site predictions. The monotonicity-constrained maximum-likelihood approach underlying the POOL method described below is built on certain assumptions relating features used for classification to the probability that an instance having those features belongs to the positive class. Here we describe in detail the form these assumptions take when relating the THEMATICS features listed above to the probability that a residue is an active-site residue. Later we will also note that similar assumptions are reasonable when considering cleft rank and sequence conservation scores and apply them to those features as well. These THEMATICS feature-based monotonicity assumptions are as follows: More precisely, for the first assumption, we treat μ3 and μ4 as measures of degree of perturbation, and for the second we treat μ3env and μ4env as measures of overall perturbation within the spatial vicinity. These assumptions then become: Given two residues in the same protein, let their corresponding 4-dimensional feature vectors be x = (x1, x2, x3, x4) and y = (y1, y2, y3, y4). If xi≤yi for each i, the probability that the first residue is an active-site residue is less than or equal to the probability that the second residue is an active-site residue. A more elegant formulation arises from the definition of a coordinate-wise partial order on the 4-dimensional feature space by x≤y iff xi≤yi for all i, and the above monotonicity assumptions then take the simple form x≤y implies P (active|x) ≤P (active|y). Finally, there is one additional subtlety that all implementations of THEMATICS have had to address, and the current approach is no exception: the need for some kind of normalization across proteins. In Ko' s approach [3], the raw features were individually transformed into Z-scores, the deviations from the mean in units of the standard deviation, as calculated for the set of all ionizable residues within a given protein. Similarly, in Tong' s SVM approach [5], the raw features were likewise transformed into robust Z-scores, defined as the deviations from the median in units of the interquartile distance for the set of all ionizable residues within a given protein. Here a very different type of transformation is applied to each feature across the population of residues within a given protein. We call this transformation rank normalization. Within each protein, each feature value is ranked from lowest to highest in that protein, and each data point is then assigned a number uniformly across the interval [0,1] based on the rank of that feature in that protein. The highest value for that feature is thus transformed to 1, and the lowest value is transformed to 0. Note that unlike the use of Z-scores or the robust Z-scores of Tong [5], this is a nonlinear transformation of the raw feature values. For each scalar feature x, denote its within-protein rank-normalized value as, which by definition lies in [0,1]. The use of this notation is extended to feature vectors in the obvious way, i. e. . Note that this rank normalization transformation does not affect the within-protein partial order used in the assumptions. That is, x≤y is true for raw feature vectors x and y in the same protein if and only if. However, when data from multiple proteins is combined for training and the results are used to make predictions for new proteins, as is described in detail below, this actually implies an even stronger monotonicity assumption across proteins in which the within-protein rankings replace the raw feature values. This assumption is harder to justify intuitively, but such an approach is required to be able to train on multiple proteins and make predictions for novel proteins, and, as shown below, it appears to give good results. After the normalization is performed, the labeled dataset may be regarded as a collection of pairs, one for each ionizable residue in the protein, where is the 4-dimensional rank-normalized feature vector for the ith residue and the label ci is either 1 (identified as an active-site residue) or 0 (not so identified). Given such a set of training data, the mathematical problem we wish to solve is to find a maximum likelihood estimator for as a function of in [0,1]4 based on this training dataset and satisfying the constraint that whenever in the coordinate-wise partial order described above. Letting n represent the number of training examples and pi the estimate of for each i from 1 to n, we seek to maximize: (2) subject to the constraints: for each (i, j) such that This is a convex optimization problem with linear constraints. We have shown [8] that the solution to this convex optimization problem is the same as the solution to the quadratic programming problem of the minimization of: (3) subject to these same constraints. Note that while the equivalence of these solutions is well-known if there are no constraints, not every constrained maximum-likelihood problem is equivalent to the corresponding minimum squared-error problem with the same constraints. However, with these particular constraints, the two solutions are indeed identical. This latter optimization problem is a special case of the general isotonic regression problem [9], [10] and this special form lends itself to a particularly straightforward solution technique. First, at an arbitrary point in the feasible region, the set of active constraints is determined by solving the corresponding dual quadratic programming problem of finding {λi} minimizingsubject to the constraints λi≥0 for all i, where is the negative gradient and are the normal vectors to the m constraint surfaces. The kth constraint in the primal problem is active iff λk>0. Furthermore, by rescaling coordinates in the primal problem, its contours become circular and the negative gradient at any point points toward a single point, the unconstrained optimum. (Thus another formulation of the primal problem is to find the point in the feasible region closest to the unconstrained optimum.) As a consequence, the active set so determined at any feasible point is exactly the same as the active set at the solution point. But the active set simply represents equivalence classes of data points for which equality of the estimates must hold. Since equality-constrained maximum-likelihood estimates have the form (number of positives) / (total number of points), identifying which constraints are active at the solution leads immediately to the solution itself. Full details of this algorithm as well as the proof that the minimum sum-of-squared-errors solution is also the maximum-likelihood solution can be found in Tong' s dissertation [8]. We call our algorithm for solving this maximum-likelihood problem the POOL algorithm. POOL is both an acronym for Partial Order Optimal Likelihood as well as an accurate characterization of the way the method first identifies the active constraints and then simply combines the corresponding data values into “pools” to be assigned probability estimates according to the proportion of positives in that pool. The use of the POOL method with this 4-dimensional THEMATICS feature vector is denoted POOL (T4) in the Results section, where its performance is compared with other methods. Previous studies have shown that active site residues tend to be located in one of the largest clefts in a protein structure [6], [7], [11]. Indeed it has been reported that in 83% of single-chain enzymes, the active site is located in the largest cleft [11]. Nearly all active sites are principally located in one of the five largest clefts of a protein structure, with the largest cleft containing the active site for the highest fraction of enzymes and with the fractions decreasing as the size rank progresses to smaller clefts [12]. If such purely geometric analyses were to be used for active site prediction, the result would be low precision and a high false positive rate, since the active site typically constitutes just a fraction of the area of the cleft. Since such geometric analyses are purely 3D-structure based, may be performed rapidly, and constitute a very different type of information from that of THEMATICS, it makes sense to combine these data in order to enhance overall performance. One straightforward way to do this is to combine the features from both THEMATICS and cleft size rank into a single vector of input to any appropriate classifier or probability estimator. Cleft size rankings may be readily incorporated into POOL, since there is an implicit monotonicity assumption that applies to this feature as well: The probability that a cleft contains an interaction site is highest for the largest cleft in a protein structure and decreases for clefts of smaller size rank. In this study we used CASTp [13], which uses computational geometry to define and measure pockets on the protein surface, to calculate cleft information for each residue in the protein. For present purposes, every residue in a given protein is assigned an integer number corresponding to the rank of the size of the cleft to which it belongs, where 1 is the largest, 2 is the second-largest, and so on. If the atoms of a residue belong to more than one cleft, the residue is assigned the rank of the largest of these clefts. When combined with THEMATICS features, the result is a 5-dimensional input vector to which coordinate-wise monotonicity constraints are applied on all five coordinates. POOL (T4, G) denotes the estimator resulting from applying the POOL method to this five-dimensional concatenation of the four THEMATICS features and one cleft size rank. An interesting alternative to simply concatenating all features into a single vector and applying a single classifier or probability estimator to such vectors is to compute two separate probability estimates and then combine them. Consider the general problem of estimating the class probability P (c|x) for a feature vector x = (x1, x2, …, xk) formed as the concatenation of feature vectors x1, x2, …, xk. It is straightforward to show that if the naïve Bayes conditional independence assumption (4) holds for each class c, then (5) where α is a normalizing constant. This gives a computationally attractive way to consider combining probability estimates for combinations of feature sets when separate estimates are available for the individual feature sets. As with other applications of naïve Bayes, it is not necessary that the conditional independence assumption be strictly true for the results of this computation to give useful results, especially when it comes to relative rankings [14]. This then gives another approach, which we have dubbed chaining, to obtain active-site probability estimates using both THEMATICS features and cleft size rank. In this case, we use equation (5) to combine the POOL estimates based on THEMATICS with the POOL estimates for the one-dimensional cleft size rank feature. POOL estimates based on the four-dimensional THEMATICS input and those based on the one-dimensional cleft size rank are labeled POOL (T4) and POOL (G), respectively, where G stands for geometry. POOL (G) gives a simple set of active-site probabilities for each ranking. The probability estimator computed using equation (5) with POOL (T4) and POOL (G) we then call POOL (T4) xPOOL (G). Later we also incorporate a conservation score feature, based on sequence alignment, using this same technique. Non-ionizable residues do not have titration curves and thus THEMATICS does not predict them directly. Nevertheless, the non-ionizable residues in interaction sites tend to have ionizable residues in their immediate vicinity and these ionizable residues generally have perturbed titration curves [1], [5]. This was the basis for the attempt by Tong et al. [5] to identify non-ionizable active site candidate residues by their proximity to the ionizable residues selected by THEMATICS. That approach, based on SVM results and called SVM-region, yields an unacceptably high false positive rate. Here we adopt a related strategy based on POOL and demonstrate substantially improved results. Note that every non-ionizable residue has the environment features μ3env and μ4env; these serve as measures of the overall amount of titration curve perturbation in their spatial neighborhood. Thus we posit an extension to the THEMATICS monotonicity assumptions, namely: All other things being equal, a non-ionizable residue having more titration curve perturbation in its neighborhood is more likely to be an active-site residue. Thus we can apply the POOL method to non-ionizable residues separately by applying coordinate-wise monotonicity constraints to the probability estimates for the 2-dimensional feature vectors, once again using the transformations, rank-normalized within each protein, of these features. In this case, the rank normalization is performed separately on just the set of non-ionizable residues in a given protein. Furthermore, we have the same options described above for incorporating cleft or other information for these non-ionizable residues. Finally, for any given protein, we can start with separate ordered lists of probability estimates for the ionizables and the non-ionizables, however computed, and then merge these into a single ordered list. This list then gives an estimated probability, and hence a ranking, for all residues. Yet another feature that is generally taken to be predictive of functional activity in a monotonic fashion is the extent to which a given residue is found to be conserved across sequence homologues: The more conserved the residue, the more likely that residue is to be functionally important in the protein. Here we also examine combining such conservation information with THEMATICS and cleft size information. In particular, we use ConSurf [15], a sequence comparison based method that identifies functionally important regions on the surface of a protein of known three-dimensional structure, based on the phylogenetic relations between its close sequence homologues. If there are a sufficient number of sufficiently diverse homologues to the query protein, Consurf assigns a score between 1 and 9 to each residue in the query sequence based on how conserved this residue is among those homologues. The more conserved a residue is, the higher its score. We call this the conservation feature and in the results below we use C to denote its inclusion. Note that taken by itself, it is also a simple one-dimensional feature, like the cleft size rank. In this study, if Consurf returns more than 10 homologues, we use the score ConSurf assigns to each residue as its conservation feature value. For any proteins with 10 or fewer homologues, all residues in that protein are assigned a single common value for that feature. The effect of this is that all residues in such proteins have a tie for that feature, so it contributes nothing to the individual probability estimates or the residue rankings within that protein. In the Results section it is shown that using information derived solely from the 3D structure, the present method outperforms all other 3D structure based methods. When sequence conservation information is included, the resulting classifier outperforms all other methods. Especially noteworthy is that in the absence of sequence conservation information, performance is nearly as good as that with such conservation information. This is particularly significant for structural genomics proteins, for which the present method is expected to perform well, even for novel folds and orphan sequences.
  The results presented here are based on several standard measures of performance. For a standard classification problem, performance is typically measured by recall (or true positive rate) and false-positive rate. Within a specific system with tunable parameters, recall and false positive rate typically involve a tradeoff: adjusting the parameters to lower the false-positive rate also lowers recall, while raising the latter also raises the former. So to judge the performance of such a system, it is important to know the tradeoff between these two, and thus ROC curves, which plot recall against false positive rate, are presented here. In the latter subsections, it is sometimes necessary to use other performance measures in order to compare our results against those reported by others. Since our method outputs a ranked list (actually a list of probabilities) for all residues within a given protein and not a binary classification, considering the ROC curves is an especially useful way to characterize the behavior of any binary classification scheme derived from it. Among the many possibilities for creating a binary classification from such a list would be to select the top n or the top p percent or use a probability threshold. One advantage of ROC curves is that they are independent of the selection scheme. One disadvantage of using ROC curves alone, however, is that unless the curve for one method dominates (i. e. , lies completely above and to the left of) that of another, there may be no simple metric to compare these two methods. For this reason, a single number that is sometimes used as a reliable measure for comparing systems in the machine learning literature is the area under the ROC curve (AUC) [18]. We make use of this as a single numerical measure to which we can apply statistical significance tests to corroborate the apparent superiority of one method over another. In order to generate ROC curves, we need to be able to calculate recall and false-positive rate values, which come from classification problems. In the POOL method, the result for each protein is a ranked list based on the probability of a residue being in the active site. A natural way to draw a ROC curve for every protein is to move the cutoff one residue at a time from the top to the bottom of the list. The resulting ROC curve has a staircase shape: only recall increases when an active site residue is encountered and only false positive rate increases when a non-active-site residue is encountered. We define average specificity (AveS) for each protein in the set: (6) where r is the rank, N is the number of residues in a protein, pos (r) is a binary function that indicates whether the residue of a given rank r is annotated in the reference database in the active site (pos (r) = 1) or not (pos (r) = 0), and S (r) is the specificity at a given cut-off rank r. (Specificity is defined to be 1 – false positive rate.) It is not hard to see that AveS represents the area under the ROC curve (AUC) for that protein. We also compute the across-protein mean of AveS over a given set of proteins, which we call the mean average specificity (MAS) for that set. To visually compare the performance from different methods, we also generate the averaged ROC curve for each method by computing the recall and false-positive rate after truncating the list after each of the positive residues in turn, followed by linearly interpolating the value at each recall value and computing the mean of the interpolated false-positive rate values across all proteins in the dataset. From these average ROC curves we can get a strong sense of the apparent relative performance of different systems, but it is also important to be able to verify that such apparent differences are in fact statistically significant. To test the significance of the observed differences, we also perform the Wilcoxon signed-rank test [19] on AveS from these methods to estimate the probability of observing such a difference under the null hypothesis that the seemingly better-performing method is actually not better than the other. This test essentially determines which method is consistently better on a protein-by-protein basis (as measured by AUC or AveS), while the curves we display essentially demonstrate which methods perform better on average.
  In this paper, we presented the application of the POOL method using THEMATICS plus some other features for protein active site prediction. We started with the application of the POOL method just on THEMATICS features, with features similar to those used before in the SVM method [5], as well as those used in Ko and Wei' s statistical analysis [3], [4]. These results show that the POOL method outperforms all of the earlier THEMATICS methods with no cleaning of the training data and no clustering after the classification. This suggests that by relying solely on the underlying THEMATICS monotonicity assumptions, the POOL method makes better use of the training data. We also tested different ways of incorporating additional features into the learning system. Not surprisingly, the results show that in order to improve performance, we have to incorporate the right features in the right way. Even with features that were found to be helpful in improving the performance, how they are incorporated matters. The data show that chaining the results from separate POOL estimates is better than simply combining all the available features into a big POOL estimator over a higher-dimensional feature space. As mentioned earlier, the reason behind this might be overfitting, since combining features into a POOL table with high dimension causes the number of probabilities needed for estimation to grow exponentially, while the training data can only increase linearly in most cases. In other words, the high dimensionality makes the table too sparse and less accurate for probability estimates. We also extended the application of THEMATICS to all residues, not just ionizable residues, in a natural way and showed that it is effective. Although the performance for non-ionizable residues is not as good as the performance for ionizable ones, this extension does provide a way to combine features from THEMATICS, which by itself can only be applied to ionizable residues directly, with some other features. The inclusion of the non-ionizable residues results in better overall performance and also makes performance comparison with other methods more accurate and fair. The incorporation of sequence conservation information does improve the predictions when there are enough homologues with appropriate diversity. The POOL method gives us a means for easily utilizing this information when it is available, while not affecting the training and classification when it is not. When comparing with other methods, especially if the other methods use binary classification instead of a ranked list, we have to commit to a specific cutoff value and turn our system into a binary classification system. The results in this paper clearly show that the POOL method using THEMATICS and geometric features achieves equivalent or better performance than the other methods in comparison, even in cases where their methods are tested on very special groups of proteins. This makes this method more widely applicable to proteins with few or no sequence homologues, such as some Structural Genomics proteins, than other methods that use sequence alignments from homologues. Performances of the previous best methods, those of Youn and of Petrova, will degrade significantly when sequence conservation information is not available. However with THEMATICS data the approach developed here is still robust in the absence of sequence conservation information. In effect, for those proteins having an insufficient number of sequence homologues, the POOL (T) xPOOL (G) xPOOL (C) method reduces to the still highly effective structure-only POOL (T) xPOOL (G) method. Interestingly enough, when comparing the performance of POOL (T) xPOOL (G) and POOL (T) xPOOL (G) xPOOL (C) in Figure 4, it is apparent that the addition of the conservation information does improve the performance a little, but not to the extent observed previously for sequence-structure methods. Typically the conservation information is the most important input feature, and without it performance is substantially worse [24]. This suggests that the 3D structure based THEMATICS features are quite powerful compared with other 3D structure based features and can take the place of conservation information for purposes of active site prediction. This is also borne out by the analysis at the end of the Results section. When looking at the recall and false positive rates of the results from all the protein active site prediction methods, one must keep in mind that the annotation of the catalytic residues in the protein dataset is never perfect. Since most of the labeling comes from experimental evidence, some active site residues are not labeled as positive simply because no experiment was ever carried out to verify the role of that specific residue. Since we have used the CatRes/CSA annotations as the sole criteria to evaluate the performance in order to keep the comparisons consistent, the reported false positive rate is probably higher than in reality. There is evidence available to support the functional importance of some residues that are not labeled as active in the CatRes/CSA database [3], [4] and these residues have high ranks in the list from the POOL method and are classified as positive by THEMATICS-SVM and THEMATICS-statistical analysis as well. Although we evaluated the POOL method performance using filtration ratio values as a cutoff, it is just for the purpose of comparing with other protein active site prediction methods that use a binary classification scheme. The ranked list of residues based on their probability of being in the active site contains much more information than traditional binary classification. The rank of the first annotated positive residue analysis in this paper shows just one application of the extra information contained in a ranked list rather than a traditional binary label. There are many possible measurements of performance depending on the actual application by users, and in turn many possible applications that can benefit from using a ranked list form. It is noteworthy that P-Cats [25] also estimates the probability that a residue belongs to a protein active site, using a k-nearest neighbor method. The P-Cats server uses the probability estimates as the basis to assign binary labels; residues with probability larger than 0. 50 are labeled as positive and the others as negative. The method of Cheng [26] also generates a rank-ordered list based on a scoring system; these scores could in principle be translated into probability estimates. The POOL approach is amenable to the addition of other properties for the prediction of active sites [27]–[32]. We also note that the POOL methodology is applicable to other types of problems in a variety of different areas where probability depends monotonically on the input feature variables. In conclusion, we have established that applying the POOL method, with THEMATICS and other features, appears to yield the best protein active site prediction system yet found and it provides more information than other active site prediction methods.
  The three-dimensional coordinate files for the protein structures used for training and testing were downloaded from the Protein Data Bank (http: //www. rcsb. org/pdb/). In order to predict the theoretical titration curve of each ionizable residue in the structure, finite-difference Poisson-Boltzmann calculations were performed using UHBD [33] on each protein followed by the program HYBRID [34], which calculates a corresponding titration curve of the form average net charge as a function of pH. These titration curves were obtained for each ionizable residue: Arg, Asp, Cys, Glu, His, Lys, Tyr, and the N- and C- termini. The pH range we simulated for all curves is from −15. 0 to 30. 0, in increments of 0. 2 pH units. This wide theoretical pH range is necessary for proper numerical integration of the first derivative functions. The structures were processed and analyzed to obtain the central moments μ3 and μ4, as described earlier. These individual features were then rank-normalized within each protein, and thus assigned values in the interval [0,1], also as described above. This four-dimensional representation constitutes what we designate the THEMATICS features for each residue. The monotonicity assumptions for this multidimensional feature set are as described earlier. For the geometric feature, we used CASTp [13], which uses a pocket algorithm for shape measurements to calculate the cleft information for each residue in the protein. The clefts were ranked based on their sizes in decreasing order and each residue having atoms located in any cleft is assigned the rank number of the largest of the clefts where its atoms are located. One special value is assigned to every residue not on the protein surface, and another is assigned to every residue on the surface but not within any cleft. Ignoring these special values, the monotonicity assumption is that the larger the cleft to which a residue belongs, the more likely that residue is to belong to the active site. For the conservation feature we used ConSurf [15] to calculate a sequence conservation score for the residues in each protein. ConSurf takes a protein sequence and finds its closest sequence homologues using MUSCLE [35], a multiple-sequence alignment algorithm. Two sequences with similarity higher than a preset threshold are treated as homologues. ConSurf analyzes the homologues of the query sequence and determines how conserved each residue is in the query protein among these homologues. In order to normalize the result and make it comparable between different proteins with different numbers of homologues and with different degrees of overall conservation, the program labels each residue with a conservation score between 1 and 9, with 9 being the most conserved and 1 being the most variable. If there exist more than 50 homologues for the query sequence, the 50 homologues closest to the query sequence are analyzed. In this study, we only used the conservation score reported by ConSurf when there are at least 11 homologues for a protein. The monotonicity assumption applied to this feature is that the larger the conservation score for a residue, the more likely that residue is to belong to the active site. The results reported here are based on eight-fold cross-validation on a set of 64 proteins or 10-fold cross-validation on a set of 160 proteins, both taken from the Catalytic Site Atlas (CSA) database [16], [17]. The labels were taken directly from the CSA database; if a residue is identified there as active in catalysis, it was labeled as positive in our dataset. If not so identified in the CSA, we labeled it as negative. The CSA annotations, although incomplete, constitute the best source of active residue labels for enzymes. In anticipation that the POOL method would not be overly sensitive to mislabeled data, no hand tuning of the labels was performed and no residues were omitted during training, in contrast to the SVM study reported by Tong [5]. For the eight-fold cross-validation procedure, the 64-protein set was randomly divided into eight folds of eight proteins each, with seven of the eight folds (56 proteins) used for training and the remaining fold (8 proteins) used for testing. This was repeated eight times, once for each of the eight folds. Likewise, for the ten-fold cross-validation procedure, the 160-protein set was randomly divided into ten folds of sixteen proteins each, with nine of these (144 proteins) used for training the remaining fold (16 proteins) used for testing, and this was repeated a total of ten times, once for each fold. Training was performed by applying the POOL method to obtain a function for each rank-normalized feature vector in the appropriate feature space [0,1]k. Note that: k = 4 for the POOL method applied on the four THEMATICS features of ionizable residues as stated earlier, denoted by POOL (T4); k = 5 for the POOL method applied on the four THEMATICS features of ionizable residues plus the geometric feature of the cleft size, denoted as POOL (T4, G); k = 1 for the POOL method applied to the geometric feature of cleft size, denoted by POOL (G), as well as the POOL method applied to the conservation feature, denoted by POOL (C); and k = 2 for POOL applied to the environmental features for non-ionizable residues, denoted as POOL (T2). An additional detail is that for training we quantize the multi-dimensional data points. For example, for POOL (T4), each rank-normalized feature falls into one of 20 bins whose sizes vary depending on their distance from 0. 0. In particular, the lowest ranked bins cover the half-open intervals, [0. 0,0. 2), [0. 2,0. 4), [0. 4,0. 6), [0. 6,0. 7), and there are 16 more bins of width 0. 02 above that, with one special bin for 1. 0. Thus the lowest-ranking data are quantized more coarsely than the remaining data. This is appropriate since these data tend to have very low average probability of being in the active site anyway, because the vast majority of residues are negatives. Thus the inability to make fine distinctions among these low-probability candidates does not degrade the overall quality of the results. It does, however, improve the efficiency of the training procedure significantly, so this is an important component of the analysis. This is especially helpful in the 10-fold cross-validation on the 160-protein set. The typical training set of 144 proteins from this set contains about 14500 ionizable residues, which fall into more than 6000 quantized bins in the 4-dimensional space used for POOL (T4). The number of corresponding inequality constraints is about 35,000–40,000. One final detail is that the probability estimates generated by the POOL method as described here tend to have numerous ties as well as some places where there is no well-defined value. The latter places occur because the method only assigns values to existing data points (or bins containing data in the case of our use of quantization). The locally constant regions occur both because of the quantization applied to the training data at the outset and because the data pools created by the algorithm acquire a single value. In cells where no value is defined, the interpolation scheme used is to simply assign a value linearly interpolated based on the Manhattan distance between the least upper bound and the greatest lower bound for that cell based on the monotonicity constraint. Finally, since both the data pooling performed by the algorithm and this interpolation scheme tend to lead to ties, the Manhattan distance from the origin of the four THEMATICS features is used as a tie-breaker for any residues whose probability estimates are identical. This simply imposes a slight bias toward strict monotonicity even though the mathematical formulation used to determine these probabilities is based on a non-strict monotonicity assumption, making it possible to obtain well-defined rankings for all the residues in a protein.
summary: Genome sequencing has revealed the codes for thousands of previously unknown proteins for humans and for hundreds of other species. Many of these proteins are of unknown or unclear function. The information contained in the genome sequences holds tremendous potential benefit to humankind, including new approaches to the diagnosis and treatment of disease. In order to realize these benefits, a key step is to understand the functions of the proteins for which these genes hold the code. A first step in understanding the function of a protein is to identify the functional site, the local area on the surface of a protein where it affects its functional activity. This paper reports on a new computational methodology to predict protein functional sites from protein 3D structures. A new machine learning approach called Partial Order Optimum Likelihood (POOL) is introduced here. It is shown that POOL outperforms previous methods for the prediction of protein functional sites from 3D structures.
keywords: biochemistry/bioinformatics
  biochemistry/theory and simulation
  computational biology/genomics
  biochemistry/structural genomics
  mathematics/statistics
title: Partial Order Optimum Likelihood (POOL): Maximum Likelihood Prediction of Protein Active Site Residues Using 3D Structure and Sequence Properties