- Methodology article
- Open access
- Published:
Reconstruction of ancestral protein sequences and its applications
BMC Evolutionary Biology volume 4, Article number: 33 (2004)
Abstract
Background
Modern-day 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 distance-based 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 position-specific 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 profile-based 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 non-commercial use. Pre-compiled versions for several platforms can be downloaded from ftp://iole.swmed.edu/pub/ANCESCON/.
Background
Present-day 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 [3–5]. 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], distance-based [9] and parsimony [1]. Distance-based 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 neighbor-joining algorithm [9]. Bruno et al. [10] introduced a distance-based 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 [12–15]. 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 site-by-site 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 branch-and-bound 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 Alignment-Based 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 site-specific 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 sub-trees but exhibit variability among different sub-trees. 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 profile-based 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 Z-score 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 Alignment-Based 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).
a) Correlation between average α AB and average observed α . b) Correlation between average α ML and average observed α . α AB is Alignment-Based 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.
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 site-by-site 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).
Site-specific 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.
The tree used to test ancestral sequence reconstruction. This is an arbitrarily selected evolutionary tree. Evolutionary distances are shown to scale.
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.
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).
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 t-test was used to calculate the one-tail 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 site-specific rate factors estimated by maximum likelihood method (α ML ), Alignment-Based 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 t-test in the third testing set, we show that ANCESCON method with α ML gives significantly better reconstruction than the other 3 methods. Because the site-specific α ML is very close to the true mutation rate at a site (Figure 1b), using the site-specific α 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.
Ancestral sequences used in homology detection
Thirty-eight 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 N-1 sequences (for details see Methods). For each combined alignment (2N-1 sequences), PSI-BLAST [35] searches were performed starting from each sequence and seeded with the combined alignment (-B option in the program BLASTPGP, e-value cutoff 0.01), and all found hits were pooled together.
The benchmark experiment was PSI-BLAST 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 PSI-BLAST 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.
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 non-native sequences to the native alignment results in a change of sequence profile for PSI-BLAST searches. Random sequences can dilute the position-specific 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 PSI-BLAST program uses composition-based statistics with e-value 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 e-values 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 OB-fold structures using reconstructed ancestral sequences.
Prediction of functional sites
Ten well-studied 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 sub-trees and high variation among sub-trees. 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 top-ranking 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.
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).
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 N-terminal 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 well-separated sub-trees. 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 substrate-binding 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.
A partial alignment of the N-terminal 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.
The evolutionary tree for the adenylyl kinase family generated by "Weighbor". The first cutting layer is shown. Evolutionary distances are shown to scale.
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 site-specific evolutionary rates (α), namely Alignment-Based 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(Q t) (2)
Here, Q is the rate matrix. The non-diagonal 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 = (u1, ..., u20) (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:
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.
Alignment-Based Rate Factor α (α AB ) and Rate factor α estimated by Maximum Likelihood (α ML )
Our program supports two methods to estimate a rate factor for each site: Alignment-Based rate factor α (α AB ) and Maximum Likelihood-estimated 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 non-gapped 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 Z-score 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 (Alignment-Based π 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 root-finding 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 site-specific 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 Z-score 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 Z-score 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. OB-fold has a 5-stranded β-barrel structure. In the SCOP (Structure Classification of Proteins) database (version 1.55) [44], there are 7 OB-fold superfamilies. The superfamily of nucleic acid binding proteins is the most populated. Diversity of many OB-fold homologs extends beyond detection by automatic PSI-BLAST searches. Multiple sequence alignments of native sequences were obtained from PSI-BLAST searches starting from the 38 OB-fold 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 N-1 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 N-1 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 N-1 sequences resulted from the shuffling to the native alignment. In the fourth case, N-1 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 sub-tree, but show high variability between different sub-trees. 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; H-Ras) [48], globin (1a6g; HEM) [49], pdz domain (1be9; C-terminal 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) single-linkage clustering and 3) complete-linkage 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 non-gapped 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 least-squares modification of the midpoint rooting procedure to define the root [57].
Tree partitioning
The tree was partitioned into sub-trees at several levels and compared the amino acid usages within each sub-tree and among the sub-trees. For this partitioning, we "cut" the tree into a fixed number of equal-distanced 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 sub-tree under it.
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 L1 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 ∈ {LK+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 sub-trees defined by the first layer. 
to 
reflect the average difference of amino acid compositions within each sub-tree. If the amino acids are highly conserved within each sub-tree but show variability among the sub-trees, 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 top-ranking 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 sub-trees at varying sequence identity thresholds. Sites that are invariant within each individual sub-tree are picked as functional sites. A higher identity threshold gives rise to more sub-trees and, since conserved sites are more frequent in the sub-trees 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 Z-score 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 Z-score 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.
References
Fitch WM: Toward defining the course of evolution: minimum change for a specific tree topology. Syst Zool. 1971, 20: 406-416.
Hartigan JA: Minimum evolution fits to a given tree. Biometrics. 1973, 29: 53-65.
Yang Z, Kumar S, Nei M: A new method of inference of ancestral nucleotide and amino acid sequences. Genetics. 1995, 141: 1641-1650.
Koshi JM, Goldstein RA: Probabilistic reconstruction of ancestral protein sequences. J Mol Evol. 1996, 42: 313-320.
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: 890-896.
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: S139-146.
Yang Z: PAML: a phylogenetic analysis by maximum likelihood. Version 2.0e. Pennsylvania State University, University Park. 1995
Felsenstein J: Evolutionary trees from DNA sequences: a maximum likelihood approach. J Mol Evol. 1981, 17: 368-376.
Saitou N, Nei M: The neighbor-joining method: a new method for reconstructing phylogenetic trees. Mol Biol Evol. 1987, 4: 406-425.
Bruno WJ, Socci ND, Halpern AL: Weighted neighbor joining: a likelihood-based approach to distance-based phylogeny reconstruction. Mol Biol Evol. 2000, 17: 189-197.
Gascuel O: BIONJ: an improved version of the NJ algorithm based on a simple model of sequence data. Mol Biol Evol. 1997, 14: 685-695.
Uzzell T, Corbin KW: Fitting discrete probability distributions to evolutionary events. Science. 1971, 172: 1089-1896.
Yang Z: Maximum-likelihood estimation of phylogeny from DNA sequences when substitution rates differ over sites. Mol Biol Evol. 1993, 10: 1396-1401.
Fitch WM, Margoliash E: Construction of phylogenetic trees. Science. 1967, 155: 279-284.
Fitch WM: The estimate of total nucleotide substitutions from pairwise differences is biased. Philos Trans R Soc Lond B Biol Sci. 1986, 312: 317-324.
Gu X, Zhang J: A simple method for estimating the parameter of substitution rate variation among sites. Mol Biol Evol. 1997, 14: 1106-1113.
Felsenstein J: Taking variation of evolutionary rates between sites into account in inferring phylogenies. J Mol Evol. 2001, 53: 447-455. 10.1007/s002390010234.
Nielsen R: Site-by-site estimation of the rate of substitution and the correlation of rates in mitochondrial DNA. Syst Biol. 1997, 46: 346-353.
Pupko T, Pe'er I, Hasegawa M, Graur D, Friedman N: A branch-and-bound algorithm for the inference of ancestral amino-acid sequences when the replacement rate varies among sites: Application to the evolution of five gene families. Bioinformatics. 2002, 18: 1116-1123. 10.1093/bioinformatics/18.8.1116.
Muller T, Vingron M: Modeling amino acid replacement. J Comput Biol. 2000, 7: 761-776. 10.1089/10665270050514918.
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: 11361-11366. 10.1073/pnas.2034878100.
Yang Z: Evalution of Several Methods for Estimating Phylogenetic Trees When Substitution Rates Differ over Nucleotide Sites. J Mol Evol. 1995, 40: 689-697.
Jones S, Thornton JM: Searching for functional sites in protein structures. Curr Opin Chem Biol. 2004, 8: 3-7. 10.1016/j.cbpa.2003.11.001.
Lichtarge O, Sowa ME: Evolutionary predictions of binding surfaces and interactions. Curr Opin Struct Biol. 2002, 12: 21-27. 10.1016/S0959-440X(02)00284-1.
Pei J, Grishin NV: AL2CO: calculation of positional conservation in a protein sequence alignment. Bioinformatics. 2001, 17: 700-712. 10.1093/bioinformatics/17.8.700.
Lichtarge O, Bourne HR, Cohen FE: An evolutionary trace method defines binding surfaces common to protein families. J Mol Biol. 1996, 257: 342-358. 10.1006/jmbi.1996.0167.
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: 512-526.
Yang Z: Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods. J Mol Evol. 1994, 39: 306-314. 10.1007/BF00178256.
Bateman A, Birney E, Cerruti L, Durbin R, Etwiller L, Eddy SR, Griffiths-Jones S, Howe KL, Marshall M, Sonnhammer EL: The Pfam protein families database. Nucleic Acids Res. 2002, 30: 276-280. 10.1093/nar/30.1.276.
Yang Z: PAML: a program package for phylogenetic analysis by maximum likelihood. Comput Appl Biosci. 1997, 13: 555-556.
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. Sinauer Associates, Sunderland, Massachusetts. 2002
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: 235-242. 10.1093/nar/28.1.235.
Murzin AG: OB(oligonucleotide/oligosaccharide binding)-fold: common structural and functional solution for non-homologous sequences. Embo J. 1993, 12: 861-867.
Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ: Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997, 25: 3389-3402. 10.1093/nar/25.17.3389.
Eddy SR: Profile hidden Markov models. Bioinformatics. 1998, 14: 755-763. 10.1093/bioinformatics/14.9.755.
Schaffer AA, Aravind L, Madden TL, Shavirin S, Spouge JL, Wolf YI, Koonin EV, Altschul SF: Improving the accuracy of PSI-BLAST protein database searches with composition-based statistics and other refinements. Nucleic Acids Res. 2001, 29: 2994-3005. 10.1093/nar/29.14.2994.
Muller-Dieckmann HJ, Schulz GE: Substrate specificity and assembly of the catalytic center derived from two structures of ligated uridylate kinase. J Mol Biol. 1995, 246: 522-530. 10.1006/jmbi.1994.0104.
Whelan S, Goldman N: A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach. Mol Biol Evol. 2001, 18: 691-699.
Press WH, Teukolsky SA, Vetterling WT, Flannery BP: Numerical Recipes in C : The Art of Scientific Computing. 1992, 451-455.
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 carboxyl-terminal leucine recognition by class I PDZ domains. J Mol Biol. 2001, 308: 963-973. 10.1006/jmbi.2001.4634.
Read RJ, James MN: Refined crystal structure of Streptomyces griseus trypsin at 1.7 A resolution. J Mol Biol. 1988, 200: 523-551.
Teplyakov A: High-resolution structure of the complex between carboxypeptidase A and L-phenyl lactate. Acta Crystallogr D Biol Crystallogr. 1993, 49: 534-540. 10.1107/S0907444993007267.
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: 536-540. 10.1006/jmbi.1995.0159.
Hannenhalli SS, Russell RB: Analysis and prediction of functional sub-types from protein sequence alignments. J Mol Biol. 2000, 303: 61-76. 10.1006/jmbi.2000.4036.
Mirny LA, Gelfand MS: Using orthologous and paralogous proteins to identify specificity-determining residues in bacterial transcription factors. J Mol Biol. 2002, 321: 7-20. 10.1016/S0022-2836(02)00587-9.
Abele U, Schulz GE: High-resolution structures of adenylate kinase from yeast ligated with inhibitor Ap5A, showing the pathway of phosphoryl transfer. Protein Sci. 1995, 4: 1262-1271.
Boriack-Sjodin PA, Margarit SM, Bar-Sagi D, Kuriyan J: The structural basis of the activation of Ras by Sos. Nature. 1998, 394: 337-343. 10.1038/28548.
Vojtechovsky J, Chu K, Berendzen J, Sweet RM, Schlichting I: Crystal structures of myoglobin-ligand complexes at near-atomic resolution. Biophys J. 1999, 77: 2153-2174.
Doyle DA, Lee A, Lewis J, Kim E, Sheng M, MacKinnon R: Crystal structures of a complexed and peptide-free membrane protein-binding domain: molecular basis of peptide recognition by PDZ. Cell. 1996, 85: 1067-1076. 10.1016/S0092-8674(00)81307-0.
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: 1037-1046. 10.1016/0092-8674(95)90219-8.
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: 584-592. 10.1038/378584a0.
Franken SM, Scheidig AJ, Krengel U, Rensland H, Lautwein A, Geyer M, Scheffzek K, Goody RS, Kalbitzer HR, Pai EF, Wittinghofer A: Three-dimensional structures and properties of a transforming and a nontransforming glycine-12 mutant of p21H-ras. Biochemistry. 1993, 32: 8411-8420.
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(c-src) SH2 domains: a thermodynamic and structural study. Biochemistry. 1997, 36: 6283-6293. 10.1021/bi970019n.
Feng S, Kapoor TM, Shirai F, Combs AP, Schreiber SL: Molecular basis for the binding of SH3 ligands with non-peptide elements identified by combinatorial synthesis. Chem Biol. 1996, 3: 661-670. 10.1016/S1074-5521(96)90134-9.
Stoll VS, Eger BT, Hynes RC, Martichonok V, Jones JB, Pai EF: Differences in binding modes of enantiomers of 1-acetamido boronic acid based protease inhibitors: crystal structures of gamma-chymotrypsin and subtilisin Carlsberg complexes. Biochemistry. 1998, 37: 451-462. 10.1021/bi971166o.
Wolf YI, Aravind L, Grishin NV, Koonin EV: Evolution of aminoacyl-tRNA synthetases--analysis of unique domain architectures and phylogenetic trees reveals a complex history of horizontal gene transfer events. Genome Res. 1999, 9: 689-710.
Goldman N: Statistical tests of models of DNA substitution. J Mol Evol. 1993, 36: 182-198.
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.
Author information
Authors and Affiliations
Corresponding author
Additional information
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.
Wei Cai, Jimin Pei contributed equally to this work.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Cite this article
Cai, W., Pei, J. & Grishin, N.V. Reconstruction of ancestral protein sequences and its applications. BMC Evol Biol 4, 33 (2004). https://doi.org/10.1186/1471-2148-4-33
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/1471-2148-4-33