Abstract
Background
Modernday proteins were selected during long evolutionary history as descendants of ancient life forms. In silico reconstruction of such ancestral protein sequences facilitates our understanding of evolutionary processes, protein classification and biological function. Additionally, reconstructed ancestral protein sequences could serve to fill in sequence space thus aiding remote homology inference.
Results
We developed ANCESCON, a package for distancebased phylogenetic inference and reconstruction of ancestral protein sequences that takes into account the observed variation of evolutionary rates between positions that more precisely describes the evolution of protein families. To improve the accuracy of evolutionary distance estimation and ancestral sequence reconstruction, two approaches are proposed to estimate positionspecific evolutionary rates. Comparisons show that at large evolutionary distances our method gives more accurate ancestral sequence reconstruction than PAML, PHYLIP and PAUP*. We apply the reconstructed ancestral sequences to homology inference and functional site prediction. We show that the usage of hypothetical ancestors together with the present day sequences improves profilebased sequence similarity searches; and that ancestral sequence reconstruction methods can be used to predict positions with functional specificity.
Conclusions
As a computational tool to reconstruct ancestral protein sequences from a given multiple sequence alignment, ANCESCON shows high accuracy in tests and helps detection of remote homologs and prediction of functional sites. ANCESCON is freely available for noncommercial use. Precompiled versions for several platforms can be downloaded from ftp://iole.swmed.edu/pub/ANCESCON/ webcite.
Background
Presentday protein sequences can be used to reconstruct ancestral sequences based on a model of sequence evolution. Such knowledge about ancestral sequences is helpful for understanding the evolutionary processes as well as the functional aspects of a protein family. Existing methods of ancestral sequence reconstruction can be divided into two main categories: Maximum Parsimony (MP) methods [1,2] and Maximum Likelihood (ML) methods [35]. MP methods do not take into account biased substitution patterns between amino acids or different tree branch lengths, and cannot distinguish those equally parsimonious reconstructions [3]. ML methods do not have these limitations and generally give more reliable results than the MP methods [6]. Yang et al. [3] first developed a ML method for ancestral sequence reconstruction. Yang [7] also made a distinction between "joint" reconstruction and "marginal" reconstruction. Joint reconstruction methods intend to find the most likely set of amino acids for all internal nodes at a site, which yields the maximum joint likelihood of the tree [5]. Marginal reconstruction compares the probabilities of different amino acids at an internal node at a site and selects the amino acid that yields the maximum likelihood for the tree at that site. Marginal reconstruction can also compute probabilities of all other amino acids for that node [4]. Koshi and Goldstein [4] developed a fast dynamic programming algorithm for marginal reconstruction in the framework of Bayesian statistics, while Pupko et al. [5] proposed a fast algorithm for joint reconstruction. The computational complexities for both algorithms scale linearly with the number of sequences. Both marginal and joint reconstruction algorithms are implemented in our program.
All reconstruction methods require a phylogenetic tree inferred from a given alignment. The quality of the tree is crucial for the reliability of reconstruction. A number of methods exist for phylogenetic inference, such as maximum likelihood [8], distancebased [9] and parsimony [1]. Distancebased methods have the advantage of being simple and are able to handle a large set of sequences. They require evolutionary distances estimated for all the sequence pairs. The most common method to infer phylogeny from distances is based on the neighborjoining algorithm [9]. Bruno et al. [10] introduced a distancebased phylogeny reconstruction method called "Weighbor", i.e. "weighted neighbor joining", which takes into account the fact that errors in distance estimates are larger for longer distances. Giving similar results, Weighbor is much faster than ML phylogeny reconstruction. It is also better than other methods such as BIONJ [11] and parsimony [1], in aspects of "long branches attract" and "long branch distracts" problems [10]. Weighbor is used in our program for phylogenetic inference.
Overwhelming evidence exists for substitution rate variation across sites [1215]. For a protein family, rate heterogeneity reflects the selective pressure imposed by folding, stability and function. Gamma distribution is widely used to model the rate variation among sites [13,16,17] because of its simplicity. Nielsen [18] suggested a method for sitebysite estimation of rate factors by a Maximum Likelihood approach. Rate variation among sites has not been taken into account in the early work of ML reconstruction of ancestral sequences [4,5]. Recently, Pupko et al. [19] introduced rate variation into joint reconstruction by a branchandbound algorithm, assuming a gamma distribution of rates among sites. In our package, two methods are proposed to estimate a rate factor for each site. The first one is based on our observation that the substitution rate at a site is correlated with the conservation of the site. The more conserved the site is in a multiple sequence alignment, the smaller its substitution rate is. This empirical method, the result of which we call AlignmentBased rate factors or α_{AB}, relies only on a multiple sequence alignment and a general model of amino acid exchange. The other one is a maximum likelihood method (α_{ML}), which requires a tree. In our implementation, we incorporate α_{AB }or α_{ML }in the joint and marginal reconstruction algorithms [4,5]. α_{AB }is also used in the Maximum Likelihood estimation of evolutionary distances [20] for tree inference.
We implement a method of evolutionary simulation that introduces sitespecific rate variations in a natural way by imposing structural and functional constraints [21]. We show by simulations that the reconstruction methods can give reasonable results and that the problem of evolutionary distance underestimation [22] is alleviated by considering rate variation across sites.
Background (or equilibrium) amino acid frequencies (π) are usually estimated from the target set of sequences or from large databases of protein families. Background amino acid frequencies estimated from a small dataset tend to have bias, while amino acid frequencies from large databases may not be suitable for the specific protein family under analysis. Here, we propose a ML method to optimize the amino acid frequency vector π. The optimized π vector can give significant improvement over the likelihood of a alignment.
Information obtained from ancestral sequence reconstruction is used for two applications: homology detection and prediction of functional sites. For homology detection, ancestral sequences represent an enlargement of the sequence space around native sequences. We demonstrate that adding reconstructed ancestral sequences to a native alignment improves the detection of homologs in database searches.
A number of methods have been developed to predict functional sites from amino acid sequences [23,24]. One simple way to infer functional sites is by positional conservation of a multiple sequence alignment [25]. Lichtarge et al. [26] proposed a method called evolutionary trace to predict functional sites by analyzing the conservation of sequence subgroups. Functional divergence during the evolutionary process can be reflected in the variation of amino acid usage across different functional subgroups. We propose a new approach that uses information from ancestral sequence reconstruction to identify sites that are well conserved within individual subtrees but exhibit variability among different subtrees. By several examples, we show that these sites frequently contribute to the functional specificity of a protein family.
Results and discussion
We developed a package (ANCESCON) to reconstruct ancestral protein sequences considering rate variation among sites. Rate factors can be estimated either by an empirical method or by a maximum likelihood method. Consideration of rate variation among sites not only improves evolutionary distance estimation, but also gives more accurate ancestral sequence reconstruction. Ancestral sequences are used to improve profilebased sequence similarity searches. We also propose a new approach to predict positions with functional specificity based on the reconstruction of ancestral sequences.
Observed α, Alignment Based Rate Factor α (α_{AB}) and Rate Factor α estimated by Maximum Likelihood (α_{ML})
Evolutionary simulations based on a Zscore model introduce rate variation across sites in a natural way by incorporating structural and functional constraints specific for a protein family [21]. The simulation procedure is a Monte Carlo simulation of the amino acid substitution process. The fixation of substitutions is dictated by a simple scoring function, which is derived from the template structure and an alignment of its homologs. The number of substitutions occurring at each site can be recorded during the simulation process and the observed α at a site equals the number of recorded substitutions at that site divided by the average substitution number for all sites. To reduce sampling variance, an average observed α vector is calculated from 100 simulations.
For the alignment consisting of all the leaf node sequences generated by the simulation process, an α_{AB }vector was calculated according to equation (11) (for details see Methods). An average α_{AB }vector was derived from 100 simulations. Correlation coefficient between the average α_{AB }vector and the average observed α vector was high (data not shown). However, we found that for large observed α values, the corresponding α_{AB }values were smaller. A constant β was introduced to correct this underestimation in equation (11).
Here, α_{i }is AlignmentBased rate factor at site i. K is the number of sites in a given alignment. C_{i }is the value assigned to site i (for details see Methods).
We optimized the β value by fitting the average α_{AB }vector and average observed α vector to y = x line. Alignments for three different protein families (trypsin, carboxypeptidase and pdz domain) gave a good empirical estimation for β of about 1.3. The relation between this corrected average α_{AB }vector and average observed α vector is shown in Figure 1a for a typical example, the pdz domain (correlation coefficient 0.973).
Figure 1. a) Correlation between average α_{AB }and average observed α. b) Correlation between average α_{ML }and average observed α. α_{AB }is AlignmentBased rate factor solely depending on the given alignment. α_{ML }is rate factor estimated by maximum likelihood method, which requires an alignment and evolutionary tree inferred from the alignment. The protein family used here is the PDZ domain.
We also estimated an α_{ML }vector for each alignment generated from the simulation (for details see Methods). The average α_{ML }vector shows good correlation with the average observed α vector (Figure 1b) (correlation coefficient 0.945). α_{AB }or α_{ML }can be incorporated in likelihood calculation in marginal or joint reconstruction. Table 1 shows that improvement of logarithm likelihood of the alignment is significant when α_{AB }or α_{ML }is used.
Table 1. Difference of logarithm likelihood and CPU time when using different α vectors
Rate variation across sites can be modeled by assuming that the rate factors follow a certain type of statistical distribution. Gamma distribution [13,27] and its discrete approximations [28] are frequently used for DNA or protein sequences. Rate variation for a protein family reflects different selective pressure at different sites to maintain structure and function. Fewer substitutions are expected to occur in more conserved sites. This hypothesis has prompted us to estimate rate factors (α_{AB}) based on sequence conservation in an empirical way. The α_{AB }is compared and calibrated using the observed α as standards. Our method of estimating α_{ML }is similar to the one proposed by Nielson [18]. One problem with sitebysite rate factor estimation is the small sample size at each site, especially with a small alignment. We have used α_{AB }to eliminate outliers with very large α_{ML }estimates (for details see Methods).
Sitespecific rate factors improve distance estimation
Evolutionary distances tend to be underestimated when rate homogeneity among sites is assumed [22]. This was tested using the simulation with structural and functional constraints. For the arbitrarily selected tree shown in Figure 2, we obtained leaf node sequences in the simulation and estimated an evolutionary distance for each sequence pair by Maximum Likelihood, either incorporating α_{AB }or setting α equal to 1.0 (equation (16)). Evolutionary distances were severely underestimated (average underestimation: 0.894) without considering rate variation among sites (Figure 3a). Introducing α_{AB }in the maximum likelihood method gave more accurate distance estimation (Figure 3b), although the distances were still underestimated, especially for small distances (average underestimation: 0.286). We believe that more accurate distances will give more accurate phylogeny reconstruction using "Weighbor" [10]. Since a tree is required to estimate α_{ML}, α_{ML }is not incorporated in estimating evolutionary distance.
Figure 2. The tree used to test ancestral sequence reconstruction. This is an arbitrarily selected evolutionary tree. Evolutionary distances are shown to scale.
Figure 3. Comparison of pairwise distances between the rebuilt tree and original tree. a) distance estimation assuming no rate variation among sites; b) distance estimation with α_{AB}. The rebuilt tree is inferred from the alignment that is generated by evolutionary simulation performed on the original tree. The original tree is arbitrarily selected.
Optimization of equilibrium frequencies
A continuous minimization method by simulated annealing was used to optimize the equilibrium frequency vector π, with the objective function being the logarithm likelihood of the alignment. Our π vector optimization program was tested on four alignments, which were taken from the SH2 and SH3 superfamilies in Pfam database (version 7.3) [29]. Two alignments from the SH2 superfamily have 44 and 87 sequences respectively and both alignment lengths are 83 amino acids (including gaps). The other two alignments from SH3 superfamily have 39 and 94 sequences respectively and both alignment lengths are 57 amino acids (including gaps). For each alignment, we ran optimization 3 times starting from different random initial points. The optimized π vectors did not converge to exactly the same point, but they had a high correlation with each other (always > 0.95) and the difference of logarithm likelihood function values was small (less than 0.1%). The logarithm likelihood of the alignment, using optimized π vector, increased slightly, but significantly (Table 2), compared with the logarithm likelihood using the π vector calculated from the alignment.
Table 2. Difference of logarithm likelihood and CPU time with and without optimization of π vector
Optimization of the π vector is time consuming. The running time for reconstruction with or without optimizing π vector is 14,902 seconds and 213 seconds for SH2 alignment (44 sequences), respectively, on a Dell PowerEdge 8450 server (CPU 700MHz, RAM 8G) (Table 2). In our program, the default π vector is calculated from the alignment while the user has the option to optimize the π vector for ancestral sequence reconstruction.
Testing reconstruction
Two different methods for simulations of the evolutionary process were used, as described in Methods, to test the reliability of the reconstruction results. In the first simulation method, starting from a randomly generated root sequence, we simulated the evolutionary process to obtain leaf node sequences based on a tree and a rate matrix. This process was repeated 100 times for a given root sequence R to produce 100 alignments consisting of all leaf node sequences. For each of the 100 alignments, we used the marginal reconstruction method to obtain an amino acid probability vector for each site at the root. To reduce sampling variance, the amino acid probability vector was averaged over the 100 simulation trials. At each site, the amino acid with the highest average probability was chosen as our result of the "reconstructed amino acid" at that site. All "reconstructed amino acids" formed the reconstructed sequences R'. There is no difference between R and R', that is, the accuracy of reconstruction is 100% for the tree shown in Figure 2. For each individual simulation and its reconstruction, we checked the amino acid with the highest probability in the reconstructed probability vector of the root. If it is indeed the "reconstructed amino acid", the prediction for that simulation is correct according to the average reconstructed results. The fraction of individual predictions that are correct according to the average reconstructed results is almost always higher than the average probability of the "reconstructed amino acid", suggesting that the average probability of the "reconstructed amino acid" gives a lower estimation of the reconstruction reliability (Figure 4a).
Figure 4. a) Correlation between the average probability of "the reconstructed amino acid" and the fraction of correct predictions. b) Correlation between the fraction of correct predictions and average α_{AB }at each site. The protein family used here is the PDZ domain. Red filled points are sites with incorrect reconstruction.
For the second simulation method, we introduced rate heterogeneity across sites with structural and functional constraints [21]. For the same tree, the accuracy of reconstruction was about 90%. Sites with larger substitution rates are expected to have less reliable reconstructions. Figure 4b shows the relationship between the average α_{AB }and the fraction of individual predictions that are correct according to the "reconstructed amino acid". Sites with incorrect "reconstructed amino acids" all have large α_{AB }values. These values reflect the difficulty of reconstructing sites with large numbers of substitutions. The probabilities of the "reconstructed amino acids" are all small for sites with incorrect reconstructions (less than 0.15), suggesting that the information content of the reconstruction is low.
The second simulation method was also used to test ANCESCON along with the reconstruction programs from PAML [30], PHYLIP [31] and PAUP* [32]. All tree topologies used in reconstruction tests were inferred from real alignments. All original root sequences were taken from PDB database [33]. We had three different types of alignment testing sets. The first testing set used the same tree topology but different root sequences to generate 100 alignments (for details see Methods). The second testing set used the same root sequence but different tree topologies. The third testing set randomly selected a root sequence and a tree topology to generate 100 alignments. After 100 alignments were generated, we reconstructed the root sequence for each alignment and found the consensus root sequence for the 100 reconstructed root sequences. Finally, the consensus root sequence was compared with the original root sequence to calculate the reconstruction accuracy, i.e. the fraction of correctly reconstructed sites for the root sequence. In addition, for the third test, the paired ttest was used to calculate the onetail probability between ANCESCON and other three methods. In order to make different tree topologies comparable, those trees were scaled to make the average distance from root to all leaf nodes (d_{a}) the same for all trees and equal to the tree of pii1 (a signal transduction protein) (d_{a }= 4.23). If d_{a }was too small (e.g. 0.5), the reconstruction accuracy was always close to 1 for all reconstruction methods used. The value d_{a }= 4.23 was large enough to generate diverse sequences to differentiate 4 different ancestral sequence reconstruction methods.
For ANCESCON we had 3 different parameter settings, which included sitespecific rate factors estimated by maximum likelihood method (α_{ML}), AlignmentBased rate factors (α_{AB}) and no rate factors (equal rates among sites). Different parameters were also used for the reconstruction programs from PAML and PHYLIP to find their best reconstructions. For PAML, reconstruction was tested with parameter α (rate factor) estimated from alignment and without α. For PHYLIP, 4 different parameter settings were tried, which were combinations of with/without α estimated from alignment by PAML and with/without branch length dwelling in input tree topology. For PAUP*, default settings were used.
Table 3 shows a comparison of the reconstruction accuracy for these 4 methods. The reconstruction accuracy of ANCESCON with α_{ML }is higher than the other three methods in almost every test. Also the reconstruction accuracy of ANCESCON with α_{AB }and without α is comparable with PAML and PHYLIP methods and is much better than PAUP*. For the first testing set, the best average accuracy for ANCESCON is about 0.5, while the best average reconstruction accuracies for PAML, PHYLIP and PAUP* are 0.45, 0.39 and 0.32 respectively. Testing set 2 and 3 produce similar results. Using the paired ttest in the third testing set, we show that ANCESCON method with α_{ML }gives significantly better reconstruction than the other 3 methods. Because the sitespecific α_{ML }is very close to the true mutation rate at a site (Figure 1b), using the sitespecific α_{ML }can improve our ability to reconstruct the amino acids for ancestral sequences correctly. These reconstruction tests suggest that ANSCESCON may be a better tool to reconstruct ancestral sequences compared to PAML, PHYLIP and PAUP* if the given alignment contains more diverse sequences.
Table 3. Ancestral sequence reconstruction accuracy by different programs
Ancestral sequences used in homology detection
Thirtyeight OB (Oligonucleotide/oligosaccharide binding)fold [34] proteins and ten other alignments (adenylyl kinase, gef, globin, pdz, ph, ptb, ras, sh2, sh3 and subtilase) from the Pfam database (version 7.3) [29] were chosen to perform homology detection tests.
Given an alignment with N sequences, we had four different methods, "BEST", "SECOND BEST", "SHUFFLE" and "RANDOM", to generate another N1 sequences (for details see Methods). For each combined alignment (2N1 sequences), PSIBLAST [35] searches were performed starting from each sequence and seeded with the combined alignment (B option in the program BLASTPGP, evalue cutoff 0.01), and all found hits were pooled together.
The benchmark experiment was PSIBLAST seeded with the native alignment (N sequences). For each type of the four combined alignments, we checked hits not found by the native alignments. New hits were verified to be true positives or false positives by running PSIBLAST or HMMER [36], followed by manual inspections.
Using the 48 native alignments, a total of 13973 hits were found by the benchmark. Compared to the benchmark, the "BEST" method detected 120 new homologs and the other three methods found 69, 74 and 9 new homologs, respectively (Figure 5). Among those new homologs, "BEST", "SECOND BEST", "SHUFFLE" and "RANDOM" methods had 3, 2, 6 and 3 false positives, respectively (Figure 5). Also, "BEST", "SECOND BEST", "SHUFFLE" and "RANDOM" methods missed 61, 1070, 60 and 7811 homologs as compared to the benchmark.
Figure 5. Comparison of "BEST", "SECOND BEST", "SHUFFLE" and "RANDOM" methods in the number of new homologs detected when compared with the benchmark experiment. The methods are defined in "Methods" section. The blue portion of the bar shows the number of true positives. The red portion of the bar shows the number of the false positives.
Adding nonnative sequences to the native alignment results in a change of sequence profile for PSIBLAST searches. Random sequences can dilute the positionspecific amino acid exchange characteristics of native alignments. This effect should not improve the profile. Indeed, few new homologs are found by the "RANDOM" method. However, sequences generated by shuffling each position of the native alignment have the same conservation properties as the native alignment, and the "SHUFFLE" method detects a total of 74 new homologs. Two effects may account for this finding. First, addition of shuffled sequences to the native alignment can slightly change the estimates of pseudocount frequencies of amino acids and thus the position specific scoring matrix [35]. Second, the new version of PSIBLAST program uses compositionbased statistics with evalue estimation related to the composition of the query sequence [37]. Each shuffled sequence has its own amino acid composition that is different from the native sequences. This difference can affect the evalues of hits. The "BEST" method detects the most number of new homologs, suggesting that the reconstructed ancestral sequences resemble the native sequences. Ancestral sequences may therefore be more similar to some remote homologs than to the native sequences. The "SECOND BEST" method detects less new homologs than the "BEST" method but more than the "RANDOM" method, suggesting that the second most probable amino acids in reconstruction can still reflect some properties of native sequences. Table 4 shows homology detection results of OBfold structures using reconstructed ancestral sequences.
Table 4. Homology detection results of OBfold structures using reconstructed ancestral sequences
Prediction of functional sites
Ten wellstudied protein families (adenylyl kinase, gef, globin, pdz, ph, ptb, ras, sh2, sh3 and subtilase) from the Pfam database (version 7.3) [29] were selected to test the prediction of functional sites. To define functional sites, we considered residues falling within 5Å of any ligand to be functionally important (i.e. AP5 for adenylyl kinase). As a simple quantification of prediction accuracy, we counted the number of predictions that lie within 5Å from the ligands and consider these sites to be true positives.
Our method intends to identify those sites with high similarity within individual subtrees and high variation among subtrees. These sites are likely to contribute to functional specificity. Based on a tree partition and the reconstructions at the cutting nodes (details see Methods), we have developed a measure called specificity score (equation (27)). We expect that both highly variable sites and highly conserved sites tend to score low in our method. Ten topranking sites were selected as our predicted functional sites for each family. For comparison, we also implemented a simple conservation (SC) method [25], the evolutionary trace (ET) method [26] and the conservation difference (CD) method [21] on the 10 protein families. The results are shown in Table 5. Here, the results from these three methods tend to include invariant or highly conserved sites while the result from our method scores those sites low. Still, the number of true positives of our method is comparable to other methods for several families. For some protein families, such as gef, pdz and subtilase, our method predicts no fewer functional residues than the other three methods.
Table 5. Comparison of the true hits among the top 10 predicted sites for ANCESCON, evolutionary trace (ET), simple conservation (SC), and conservation difference (CD) methods
Figure 6 shows the mapping of our predictions on the structure for the PDZ domain family. In green color is the ligand and in red color are the functional residues predicted by our method. Six of the predicted residues are within 5Å to the peptide ligand. Nine of the predicted residues are around the ligand binding area. Only one is distant from the ligand (Figure 6).
Figure 6. Mapping top 10 predictions by ANCESCON to PDZ domain (PDB ID: 1be9) [50]. The color code scheme: ligand is shown in green and the predicted functional residues are shown in red.
Another example is the family of adenylyl kinases. Our method identified 3 residues within 5 Å to the ligand while the other 3 methods identified more such residues, most of which are in highly conserved positions such as the catalytic residues. Highly conserved residues, however, are not selected by our method since our measure is designed to emphasize on sites contributing to specificity. Figure 7 shows the Nterminal part of the alignment of adenylyl kinases, with our predictions highlighted in red and orange colors. The evolutionary tree for the alignment is shown in Figure 8. The first cutting layer (for details see Methods) results in two wellseparated subtrees. Functional annotations suggest that they contain enzymes with different substrate preferences: adenylyl kinases and uridylate kinases, respectively. Three residues (27, 54 and 89) from our predictions (red colored in Figure 9) contribute to substratebinding specificity, as have been noted before in the structural studies of the UMP kinases [38]. Figure 9 highlights our predicted functional residues on the adenylyl kinase protein structure. Most of our predictions fall within the specificity pocket.
Figure 7. A partial alignment of the Nterminal part of adenylyl kinases. Sites colored in red are our predictions that are within 5Å from the ligand. Sites colored in orange are our predictions more than 5Å apart from the ligand.
Figure 8. The evolutionary tree for the adenylyl kinase family generated by "Weighbor". The first cutting layer is shown. Evolutionary distances are shown to scale.
Figure 9. Mapping top 10 predictions by ANCESCON to adenylyl kinase domain (PDB ID: 1aky) [47]. The color code scheme: ligand is shown in green and the predicted functional residues are shown in red.
Conclusions
We developed a package (ANCESCON) to reconstruct ancestral protein sequences that takes into account the variation of substitution rates among sites. Two methods were proposed to estimate sitespecific evolutionary rates (α), namely AlignmentBased rate factor (α_{AB}) and rate factor α estimated by maximum likelihood (α_{ML}). Consideration of rate variation among sites can alleviate the underestimation of evolutionary distances. Accuracy of ancestral sequence reconstruction by our method is higher than that of PAML, PHYLIP and PAUP* when the given alignment contains more diverse sequences. We show that reconstructed ancestral sequences help to improve detection of distant homologs and prediction of functional sites with specificity.
Methods
Transition probability and likelihood calculations
For all models discussed in this paper, we assume all sites in an alignment evolve independently and according to a homogeneous, stationary and time reversible Markov process. The probability of an amino acid i to be replaced by amino acid j after a time interval t is P_{ij}(t). The transition probability matrix of 20 amino acids is written as P(t), which can be calculated as
P(t) = exp(Qt) (2)
Here, Q is the rate matrix. The nondiagonal elements q_{ij }are the instantaneous rates of change from amino acid i to amino acid j and diagonal elements q_{ii }are such that each matrix row sums up to 0. Q can be calculated by:
Q = S* diag(π) (3)
S is the matrix of amino acid exchangeability parameters [39]. π_{i }is the equilibrium frequency for amino acid i. Time reversibility implies that S is a symmetric matrix. In our program, the S matrix is taken from Whelan and Goldman [39] and the default π vector is estimated from the given alignment.
Q can be decomposed into eigenvalues (λ_{i}) and eigenvectors (u_{i}).
U = (u_{1}, ..., u_{20}) (5)
P_{ij}(t) can be calculated using the following equation,
The likelihood function [8] for an evolutionary tree T shown in Figure 10 is:
Figure 10. An evolutionary tree topology. Nodes C, D, E and F represent given protein sequences, while nodes A and B represent ancestral protein sequences, i.e. unknown sequences. d_{YZ }represents the evolutionary distance between nodes Y and Z.
Here, is the equilibrium frequency of the amino acid at the node A. is the transition probability from the amino acid at node A to the amino acid at node B after an evolutionary distance d_{AB}.
Considering that each site i has a rate factor α_{i }[13,18], we have:
t in equation (6) can be expressed as:
t = α·d (9)
d is the evolutionary distance and α is rate factor. The following restriction on the vector α holds:
Here, K is the number of sites.
AlignmentBased Rate Factor α (α_{AB}) and Rate factor α estimated by Maximum Likelihood (α_{ML})
Our program supports two methods to estimate a rate factor for each site: AlignmentBased rate factor α (α_{AB}) and Maximum Likelihoodestimated rate factor α (α_{ML}).
The estimation of α_{AB }is empirical and based on the observation that the substitution rate at a site is correlated with the conservation of the site, which, in turn, is correlated with the average transition probability among the amino acids at that site. Conserved sites are dominated by highly similar amino acids and thus have high average transition probabilities among the amino acids. The algorithm to calculate α_{AB }is as follows:
1. Set t equal to 1.0 and use equation (6) to calculate a transition probability matrix P for 20 amino acids. Equation, , is used to compute a symmetric matrix P'.
2. Calculate the average transition probability for each site and take the reciprocal: , where is the number of nongapped amino acid pairs in site i and the denominator is the sum over the transition probabilities between all amino acid pairs (j,k) at a site i.
3. For invariant sites, C_{i }is set to 0 to make it consistent with the Maximum Likelihood estimation.
4. Equation (11) is used to calculate α_{AB}, so that equation (10) holds.
If an evolutionary tree is assumed for the alignment, we can estimate the α_{ML }factors by maximizing the likelihood (equation (8)) for each site:
If some sites are highly variable, the α_{ML }at those sites can be very large, as has been previously noticed [18]. We consider these rate factors to be outliers. For these sites, we have observed that likelihood changes very little over a wide range of the α values. An empirical method is used to reduce the values of α_{ML }outliers, guided by the α_{AB }values. a Zscore of the ratio of α_{ML }to α_{AB }is calculated for each site except invariant sites:
Here, is the ratio of to for site i; is the number of sites excluding the invariant sites. If Z_{i }is greater than 3, it is reduced to 3 by decreasing the value of . We repeat this procedure until no Z_{i }for any site i is greater than 3. After removing the outliers, we scale the values so that equation (10) holds.
Amino acid frequency vector π optimization
Two methods are implemented to estimate the equilibrium frequency vector π, one derived directly from the given alignment (AlignmentBased π or π_{AB}) and the other estimated by Maximum Likelihood (π_{ML}). The likelihood for the entire alignment is a function of π with 19 variables. A continuous minimization method by simulated annealing [40] is used to optimize π, with the objective function being the logarithm likelihood of the alignment. The simulated annealing is computationally intensive and is the major reason for the long CPU time given in Table 2.
Distance matrix calculation and tree inference
A Maximum Likelihood approach is used to estimate the evolutionary distances among sequences, either considering rate variation across sites or not. The logarithm likelihood for replacing one protein sequence (A) with another protein sequence (B) after an evolutionary distance d can be written as:
Here, is the equilibrium frequency for the amino acid at site j in sequence A. is the transition probability from amino acid at site j in sequence A to amino acid at site j in sequence B after an evolutionary distance α _{j}·d. α_{j }is 1 if all sites are assumed to evolve at the same rate; otherwise the α_{AB }at site j is used for α_{j}.
An estimate of the evolutionary distance between two sequences is obtained by maximizing the likelihood function of equation (15):
Equation (16) can be solved by the bisection rootfinding method [40].
After the distance matrix is calculated, the "Weighbor" method, i.e. weighted neighbor joining, is used to infer an evolutionary tree [10].
Ancestral sequence reconstruction
Two methods are implemented to reconstruct ancestral sequences. One is a marginal reconstruction method [4], and the other is a joint reconstruction method [5]. Below are their brief descriptions.
The marginal reconstruction method [4]
We calculate P(A_{r}{A_{l}}T), which is the conditional probability of amino acid A_{r }at the root, given leaf node amino acid set {A_{l}} and a tree T. Since time reversibility is assumed, any internal node can serve as a root. Using Bayes' theorem, we have:
Here, P(A_{r}) is used here instead of P(A_{r}T) because the frequency of the root amino acid A_{r}, i.e. π_{r}, does not depend on tree T. P({A_{l}}A_{r}T) is the conditional probability of the known amino acids at the leaf nodes, given T and A_{r}. P({A_{l}}T) does not depend on A_{r}, so it is calculated as a normalization constant for P(A_{r}{A_{l}},T) terms over all 20 possible values of A_{r }to make the sum equal to 1.
For Figure 10, P({A_{l}}A_{r}T) can be expanded as:
Here, is the transition probability from amino acid at node A to amino acid at node B after an evolutionary distance d_{AB}. Equation (18) can be calculated using a recursive method suggested by Felsenstein [8].
If rate factors are used in the reconstruction of the root sequence, we have:
Here, α_{i }could be either α_{AB }or α_{ML }at site i. P(A_{C},A_{D},A_{E},A_{F } A_{A},T)_{i }is the conditional probability P(A_{C},A_{D},A_{E},A_{F } A_{A},T) at site i.
The joint reconstruction method [5]
The objective of a joint reconstruction method is to find the combination of amino acids for an internal node set {A_{i}} that maximize the conditional probability of this amino acid combination, given the leaf node amino acid set {A_{l}} and a tree T, P({A_{i}}{A_{l}},T). Using the Bayes' theorem, we have:
Because P({A_{l}}T) is the same for all amino acid combination at internal node set {A_{i}} this problem becomes finding the maximum of P({A_{l}}{A_{i}},T) *P({A_{i}}).
The details of a fast algorithm to solve equation (20) can be found in Pupko et al. [5]. We also incorporated sitespecific rate factors in this algorithm, in a similar way as equation (19)
Gaps
Due to difficulties with the probabilistic models of gaps, a simplified empirical approach is used to alleviate the problem. We assume that gaps are "supersede" letters. Gaps are considered for each site independently. If a leaf node has a gap instead of an amino acid at a site, this node will be deleted from the tree for this site. After dealing with leaves, we check all internal nodes for children. If an internal node has no children or only one child due to the leaf removal because of gaps, it will be removed from the tree and a gap will be assumed as its reconstructed state.
Simulations of evolutionary process
Two methods of simulating amino acid substitution process were used to test the reliability of reconstruction, rate factors and evolutionary distance estimation. The first simulation method was based on a homogeneous time reversible Markov model. The parameters from Whelan and Goldman [39] were chosen for our model, including the equilibrium frequency vector π and the S matrix. Given the length of a branch from a parent node to one of its child nodes and the amino acid for the parent node, we simulated the substitution process to generate an amino acid for the child node based on the transition probabilities that were calculated using equation (6). For the arbitrarily selected tree shown in Figure 2, we first generated a random sequence of 100 amino acids as the root sequence based on the amino acid frequencies from Whelan and Goldman [39]. We then simulated the random substitution process to obtain all leaf node sequences. This simulation was repeated 100 times. The resulting 100 alignments were used to test the reliability of the reconstruction result. In this simulation, each site evolved independently according to the same tree topology and branch lengths, thus there was no rate heterogeneity across sites.
The second simulation method, based on a Zscore model, introduced rate variation across sites by using structural and functional information for a specific protein family [21]. We selected three protein families for the Zscore simulations under structural and functional constraints: pdz domain (Protein DataBank (PDB) ID: 1g9o) [41], trypsin (PDB ID: 1sgt) [42] and carboxypeptidase A (PDB ID: 2ctb) [43]. Given a rooted tree, the native sequence with known structure was used as the root sequence. Simulations were made along the tree to generate sequences at any internal node or leaf node. If the evolutionary distance from a parent node to a child node was d, the child sequence was obtained after l*d accepted substitutions starting from the parent sequence, where l is protein sequence length. Simulations of the substitution process were repeated 100 times. For each site, the number of accepted substitutions was recorded and averaged over 100 simulations. Rate factors (observed α), representing site mutability, were calculated from these average substitution numbers, such that the average of rate factors is 1 (equation (10)). 100 simulated alignments were used to test the rate factor estimators (α_{AB }and α_{ML}), distance calculation methods and ancestral sequence reconstruction.
Homology detection
Testing dataset
38 OB (Oligonucleotide/oligosaccharide binding)fold [34] proteins with known structures were selected for homology detection test. OBfold has a 5stranded βbarrel structure. In the SCOP (Structure Classification of Proteins) database (version 1.55) [44], there are 7 OBfold superfamilies. The superfamily of nucleic acid binding proteins is the most populated. Diversity of many OBfold homologs extends beyond detection by automatic PSIBLAST searches. Multiple sequence alignments of native sequences were obtained from PSIBLAST searches starting from the 38 OBfold sequences with known structures. We also selected 10 alignments (adenylyl kinase, gef, globin, pdz, ph, ptb, ras, sh2, sh3 and subtilase) from the Pfam database (version 7.3) [29] for homology detection test.
Four different methods
For each alignment with N sequences, ancestral sequences for the N1 internal nodes were reconstructed. The idea is to test whether adding more sequences to a native alignment can help homology detection. Four types of combined alignments were generated, adding different sets of N1 sequences to the native alignment. In the first case, the added sequence at each internal node consisted of amino acids with the largest probability at each position. In the second case, the added sequences were made up of amino acids with the second largest probability. In the third case, we shuffled the native alignment at each position while keeping the gap pattern as in the native alignment. After shuffling, we added N1 sequences resulted from the shuffling to the native alignment. In the fourth case, N1 random sequences were generated with the overall amino acid frequencies of the native alignment. These four methods are named "BEST", "SECOND BEST", "SHUFFLE" and "RANDOM", respectively.
Prediction of functional sites
Our objective is to find sites that are well conserved within each subtree, but show high variability between different subtrees. These sites are likely to contribute to functional specificity [26,45,46].
Sequence datasets
Multiple sequence alignments of ten protein families were chosen from the Pfam database (version 7.3) [29]. These families are: adenylyl kinase (adkinase) (representing structure PDB ID: 1aky; its ligand or substrate: AP5) [47], guanine nucleotide exchange factor (gef) (1bkd; HRas) [48], globin (1a6g; HEM) [49], pdz domain (1be9; Cterminal peptide of protein CRIPT) [50], ph domain (1mai; I3P) [51], ptb domain (1shc; PTR) [52], ras (821p; GTN) [53], sh2 domain (1a09; ACE) [54], sh3 domain (1nlo; ACE) [55] and subtilase (1av7; SBL) [56]. Most of these alignments contain many sequences. We pruned and clustered the sequences in each alignment according to the length and diversity. Representative sequences were kept and used for tree inference and ancestral sequence reconstruction. This procedure was done in three steps: 1) removing fragments, 2) singlelinkage clustering and 3) completelinkage clustering, as described below.
1. For each family, there is a template sequence with known structure. The sequences, which cover less than 75% of the nongapped positions in the template sequence with amino acids, were considered to be fragments and discarded.
2. A sequence identity matrix was calculated for the remaining sequences. A single linkage clustering was done to form sequence groups at sequence identity threshold 0.8. For each group, we chose the longest sequence as a representative, discarding other members. This step reduced redundancy in the dataset.
3. An average sequence identity was calculated for the remaining sequences. We used this average identity as a threshold for complete linkage clustering to form new sequence groups. Four groups with the largest sequence numbers were chosen to form our new alignment. Any group with the same number of sequences as the fourth group was also included in the new alignment. The purpose of this step is to keep the major sequence subgroups of a family while leaving out highly divergent sequences that might be deleterious for tree inference.
Rooting
The "Weighbor" method gives an unrooted tree. For our purpose of predicting functional sites, we need to find a point on the tree that serves as the root. We used a leastsquares modification of the midpoint rooting procedure to define the root [57].
Tree partitioning
The tree was partitioned into subtrees at several levels and compared the amino acid usages within each subtree and among the subtrees. For this partitioning, we "cut" the tree into a fixed number of equaldistanced layers, using the midpoint as the root (Figure 11). Several criteria were tried for selecting the distance between adjacent layers. Empirically we found that a simple partition of the tree into 5 layers usually gave the best results. If the average distance from the root to all leaf nodes is d_{r}, then the distance between adjacent layers is d_{r}/5 (Figure 11). Each place of a "cut" between the layers corresponds to a certain ancestral sequence. We term the location of a "cut" as a "cutting" node. The marginal reconstruction method was used to reconstruct amino acid probability vectors for all the cutting nodes (Figure 11). The reconstructed probability vector of a cutting node reflects the amino acid usages of the subtree under it.
Figure 11. An example showing the different cutting layers in a rooted tree. d_{r }is the average distance from the root to all leaf nodes. Nodes i and j are neighboring cutting nodes.
Calculating specificity score for each site
We use {L_{K}} to represent the set of cutting nodes for layer L_{K}, K = 0,1,5. {L_{0}} is the root and L_{1 }is the closest layer to the root, etc.
A dissimilarity score between any neighboring cutting node pair is calculated. The definition of a neighboring cutting node pair (i, j) (Figure 11) is:
1. i ∈ {L_{K}}
2. j ∈ {L_{K+1}}
3. Node i is an ancestor of node j (all points on the path from j to root node are ancestors of node j), so that the distance between i and j is exactly d_{r}/5. Each cutting node has only one ancestral cutting node neighbor.
The dissimilarity score for cutting node j and its ancestral cutting node neighbor i, i.e. anc(j), at site m is defined as:
and are the reconstructed probabilities of amino acid A at cutting node j and its ancestral cutting node neighbor i(anc(j)), respectively.
Let , K = 1,...,5 (22)
Here, is the average dissimilarity score for layer K . N_{K }is the number of cutting nodes in layer K.
The specificity score is defined as:
reflects the difference of amino acid compositions among the major subtrees defined by the first layer. to reflect the average difference of amino acid compositions within each subtree. If the amino acids are highly conserved within each subtree but show variability among the subtrees, to are small and is large, leading to a large value of S_{m}. We set S_{m }to 0 for invariant sites. We sort the sites by their specificity scores and choose the 10 top scoring sites as our predicted functional sites. Those predicted functional sites that lie within 5 Å from the ligand(s) are considered to be true positives.
Comparison with other methods
We compared our method with three other methods for prediction of functional sites. The first method (Simple Conservation or SC) is based on sequence conservation. Highly conserved sites are considered to be functional. For each family, we sorted the sites by positional conservation [25] and chose the 10 topranking sites as the predictions. There might be ties for sites. For example, if there were 5 sites tied at the tenth conservation value and only one of them was within 5Å from the ligand(s), then its contribution to the total number of "correct predictions" was 1/5. The second method is the evolutionary trace (ET) method [26], which partitions a sequence identity dendrogram into subtrees at varying sequence identity thresholds. Sites that are invariant within each individual subtree are picked as functional sites. A higher identity threshold gives rise to more subtrees and, since conserved sites are more frequent in the subtrees with smaller sizes, lead to more predicted sites. ET analysis was performed from a low identity threshold to higher thresholds until the number of predicted sites was 10 or just above 10 (in the cases of ties). Ties were resolved similarly to the simple conservation method. The third method (conservation difference or CD) is based on the conservation differences between a native alignment and an alignment derived from the Zscore sequence design [21]. The basic idea was to differentiate sites conserved due to structural stability and sites conserved due to function. Since the pairwise potential in the Zscore design tends to weaken the conservation caused by function, functionally conserved sites tend to have a large conservation difference between the native alignment and the alignment of designed sequences. We chose 10 top ranking sites sorted by conservation difference as predictions by CD.
Authors' contributions
NVG conceived and initiated the study. All authors took part in developing methods and designing experiments. WC wrote the source code and JP analyzed the data. All authors read and approved the final manuscripts.
Acknowledgements
We thank Lisa Kinch, James Wrabl and Hua Cheng for their useful comments. This work was supported by the NIH grant GM67165 to NVG.
References

Fitch WM: Toward defining the course of evolution: minimum change for a specific tree topology.

Yang Z, Kumar S, Nei M: A new method of inference of ancestral nucleotide and amino acid sequences.
Genetics 1995, 141:16411650. PubMed Abstract  Publisher Full Text

Koshi JM, Goldstein RA: Probabilistic reconstruction of ancestral protein sequences.
J Mol Evol 1996, 42:313320. PubMed Abstract  Publisher Full Text

Pupko T, Pe'er I, Shamir R, Graur D: A fast algorithm for joint reconstruction of ancestral amino acid sequences.
Mol Biol Evol 2000, 17:890896. PubMed Abstract  Publisher Full Text

Zhang J, Nei M: Accuracies of ancestral amino acid sequences inferred by the parsimony, likelihood, and distance methods.
J Mol Evol 1997, 44 Suppl 1:S139146. PubMed Abstract  Publisher Full Text

Yang Z: PAML: a phylogenetic analysis by maximum likelihood. Version 2.0e.

Felsenstein J: Evolutionary trees from DNA sequences: a maximum likelihood approach.
J Mol Evol 1981, 17:368376. PubMed Abstract

Saitou N, Nei M: The neighborjoining method: a new method for reconstructing phylogenetic trees.
Mol Biol Evol 1987, 4:406425. PubMed Abstract  Publisher Full Text

Bruno WJ, Socci ND, Halpern AL: Weighted neighbor joining: a likelihoodbased approach to distancebased phylogeny reconstruction.
Mol Biol Evol 2000, 17:189197. PubMed Abstract  Publisher Full Text

Gascuel O: BIONJ: an improved version of the NJ algorithm based on a simple model of sequence data.
Mol Biol Evol 1997, 14:685695. PubMed Abstract  Publisher Full Text

Uzzell T, Corbin KW: Fitting discrete probability distributions to evolutionary events.
Science 1971, 172:10891896. PubMed Abstract

Yang Z: Maximumlikelihood estimation of phylogeny from DNA sequences when substitution rates differ over sites.
Mol Biol Evol 1993, 10:13961401. PubMed Abstract  Publisher Full Text

Fitch WM, Margoliash E: Construction of phylogenetic trees.
Science 1967, 155:279284. PubMed Abstract

Fitch WM: The estimate of total nucleotide substitutions from pairwise differences is biased.
Philos Trans R Soc Lond B Biol Sci 1986, 312:317324. PubMed Abstract

Gu X, Zhang J: A simple method for estimating the parameter of substitution rate variation among sites.
Mol Biol Evol 1997, 14:11061113. PubMed Abstract  Publisher Full Text

Felsenstein J: Taking variation of evolutionary rates between sites into account in inferring phylogenies.
J Mol Evol 2001, 53:447455. PubMed Abstract  Publisher Full Text

Nielsen R: Sitebysite estimation of the rate of substitution and the correlation of rates in mitochondrial DNA.
Syst Biol 1997, 46:346353. PubMed Abstract

Pupko T, Pe'er I, Hasegawa M, Graur D, Friedman N: A branchandbound algorithm for the inference of ancestral aminoacid sequences when the replacement rate varies among sites: Application to the evolution of five gene families.
Bioinformatics 2002, 18:11161123. PubMed Abstract  Publisher Full Text

Muller T, Vingron M: Modeling amino acid replacement.
J Comput Biol 2000, 7:761776. PubMed Abstract  Publisher Full Text

Pei J, Dokholyan NV, Shakhnovich EI, Grishin NV: Using protein design for homology detection and active site searches.
Proc Natl Acad Sci U S A 2003, 100:1136111366. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Yang Z: Evalution of Several Methods for Estimating Phylogenetic Trees When Substitution Rates Differ over Nucleotide Sites.

Jones S, Thornton JM: Searching for functional sites in protein structures.
Curr Opin Chem Biol 2004, 8:37. PubMed Abstract  Publisher Full Text

Lichtarge O, Sowa ME: Evolutionary predictions of binding surfaces and interactions.
Curr Opin Struct Biol 2002, 12:2127. PubMed Abstract  Publisher Full Text

Pei J, Grishin NV: AL2CO: calculation of positional conservation in a protein sequence alignment.
Bioinformatics 2001, 17:700712. PubMed Abstract  Publisher Full Text

Lichtarge O, Bourne HR, Cohen FE: An evolutionary trace method defines binding surfaces common to protein families.
J Mol Biol 1996, 257:342358. PubMed Abstract  Publisher Full Text

Tamura K, Nei M: Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees.
Mol Biol Evol 1993, 10:512526. PubMed Abstract  Publisher Full Text

Yang Z: Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods.
J Mol Evol 1994, 39:306314. PubMed Abstract  Publisher Full Text

Bateman A, Birney E, Cerruti L, Durbin R, Etwiller L, Eddy SR, GriffithsJones S, Howe KL, Marshall M, Sonnhammer EL: The Pfam protein families database.
Nucleic Acids Res 2002, 30:276280. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Yang Z: PAML: a program package for phylogenetic analysis by maximum likelihood.
Comput Appl Biosci 1997, 13:555556. PubMed Abstract

Felsenstein J: PHYLIP (Phylogeny Inference Package), version 3.6b.
Department of Genetics, University of Washington, Seattle 2004.

Swofford DL: PAUP*: Phylogenetic Analysis Using Parsimony (*and Other Methods). Version 4.

Berman HM, Westbrook J, Feng Z, Gilliland G, Bhat TN, Weissig H, Shindyalov IN, Bourne PE: The Protein Data Bank.
Nucleic Acids Res 2000, 28:235242. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Murzin AG: OB(oligonucleotide/oligosaccharide binding)fold: common structural and functional solution for nonhomologous sequences.
Embo J 1993, 12:861867. PubMed Abstract

Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ: Gapped BLAST and PSIBLAST: a new generation of protein database search programs.
Nucleic Acids Res 1997, 25:33893402. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Eddy SR: Profile hidden Markov models.
Bioinformatics 1998, 14:755763. PubMed Abstract  Publisher Full Text

Schaffer AA, Aravind L, Madden TL, Shavirin S, Spouge JL, Wolf YI, Koonin EV, Altschul SF: Improving the accuracy of PSIBLAST protein database searches with compositionbased statistics and other refinements.
Nucleic Acids Res 2001, 29:29943005. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

MullerDieckmann HJ, Schulz GE: Substrate specificity and assembly of the catalytic center derived from two structures of ligated uridylate kinase.
J Mol Biol 1995, 246:522530. PubMed Abstract  Publisher Full Text

Whelan S, Goldman N: A general empirical model of protein evolution derived from multiple protein families using a maximumlikelihood approach.
Mol Biol Evol 2001, 18:691699. PubMed Abstract  Publisher Full Text

Press WH, Teukolsky SA, Vetterling WT, Flannery BP: Numerical Recipes in C : The Art of Scientific Computing.
1992, 451455.

Karthikeyan S, Leung T, Birrane G, Webster G, Ladias JA: Crystal structure of the PDZ1 domain of human Na(+)/H(+) exchanger regulatory factor provides insights into the mechanism of carboxylterminal leucine recognition by class I PDZ domains.
J Mol Biol 2001, 308:963973. PubMed Abstract  Publisher Full Text

Read RJ, James MN: Refined crystal structure of Streptomyces griseus trypsin at 1.7 A resolution.
J Mol Biol 1988, 200:523551. PubMed Abstract

Teplyakov A: Highresolution structure of the complex between carboxypeptidase A and Lphenyl lactate.
Acta Crystallogr D Biol Crystallogr 1993, 49:534540. PubMed Abstract  Publisher Full Text

Murzin AG, Brenner SE, Hubbard T, Chothia C: SCOP: a structural classification of proteins database for the investigation of sequences and structures.
J Mol Biol 1995, 247:536540. PubMed Abstract  Publisher Full Text

Hannenhalli SS, Russell RB: Analysis and prediction of functional subtypes from protein sequence alignments.
J Mol Biol 2000, 303:6176. PubMed Abstract  Publisher Full Text

Mirny LA, Gelfand MS: Using orthologous and paralogous proteins to identify specificitydetermining residues in bacterial transcription factors.
J Mol Biol 2002, 321:720. PubMed Abstract  Publisher Full Text

Abele U, Schulz GE: Highresolution structures of adenylate kinase from yeast ligated with inhibitor Ap5A, showing the pathway of phosphoryl transfer.
Protein Sci 1995, 4:12621271. PubMed Abstract

BoriackSjodin PA, Margarit SM, BarSagi D, Kuriyan J: The structural basis of the activation of Ras by Sos.
Nature 1998, 394:337343. PubMed Abstract  Publisher Full Text

Vojtechovsky J, Chu K, Berendzen J, Sweet RM, Schlichting I: Crystal structures of myoglobinligand complexes at nearatomic resolution.
Biophys J 1999, 77:21532174. PubMed Abstract  Publisher Full Text

Doyle DA, Lee A, Lewis J, Kim E, Sheng M, MacKinnon R: Crystal structures of a complexed and peptidefree membrane proteinbinding domain: molecular basis of peptide recognition by PDZ.
Cell 1996, 85:10671076. PubMed Abstract  Publisher Full Text

Ferguson KM, Lemmon MA, Schlessinger J, Sigler PB: Structure of the high affinity complex of inositol trisphosphate with a phospholipase C pleckstrin homology domain.
Cell 1995, 83:10371046. PubMed Abstract  Publisher Full Text

Zhou MM, Ravichandran KS, Olejniczak EF, Petros AM, Meadows RP, Sattler M, Harlan JE, Wade WS, Burakoff SJ, Fesik SW: Structure and ligand recognition of the phosphotyrosine binding domain of Shc.
Nature 1995, 378:584592. PubMed Abstract  Publisher Full Text

Franken SM, Scheidig AJ, Krengel U, Rensland H, Lautwein A, Geyer M, Scheffzek K, Goody RS, Kalbitzer HR, Pai EF, Wittinghofer A: Threedimensional structures and properties of a transforming and a nontransforming glycine12 mutant of p21Hras.
Biochemistry 1993, 32:84118420. PubMed Abstract

Charifson PS, Shewchuk LM, Rocque W, Hummel CW, Jordan SR, Mohr C, Pacofsky GJ, Peel MR, Rodriguez M, Sternbach DD, Consler TG: Peptide ligands of pp60(csrc) SH2 domains: a thermodynamic and structural study.
Biochemistry 1997, 36:62836293. PubMed Abstract  Publisher Full Text

Feng S, Kapoor TM, Shirai F, Combs AP, Schreiber SL: Molecular basis for the binding of SH3 ligands with nonpeptide elements identified by combinatorial synthesis.
Chem Biol 1996, 3:661670. PubMed Abstract  Publisher Full Text

Stoll VS, Eger BT, Hynes RC, Martichonok V, Jones JB, Pai EF: Differences in binding modes of enantiomers of 1acetamido boronic acid based protease inhibitors: crystal structures of gammachymotrypsin and subtilisin Carlsberg complexes.
Biochemistry 1998, 37:451462. PubMed Abstract  Publisher Full Text

Wolf YI, Aravind L, Grishin NV, Koonin EV: Evolution of aminoacyltRNA synthetasesanalysis of unique domain architectures and phylogenetic trees reveals a complex history of horizontal gene transfer events.
Genome Res 1999, 9:689710. PubMed Abstract  Publisher Full Text

Goldman N: Statistical tests of models of DNA substitution.
J Mol Evol 1993, 36:182198. PubMed Abstract