Site-Specific Amino Acid Replacement Matrices from Structurally Constrained Protein Evolution Simulations

Evolutionary amino acid replacement rates depend on local structural environment (Overington et al. 1992 ; Koshi and Goldstein 1995 ). Recent models of protein evolution aim to take such site heterogeneity into account by using site-specific amino acid replacement matrices (Lio and Goldman 1998 ; Thorne 2000 ). This is a difficult task, mainly because of a lack-of-data problem: too many sequences would be needed to fit the large number of parameters required to specify the replacement matrix for each site. Two strategies have been used to tackle this problem. A first type of model assumes that protein sites can be classified into a limited number of structural classes. Then, class-specific replacement matrices are obtained by analyzing large databases of sequences (Thorne, Goldman, and Jones 1996 ; Lio et al. 1998 ). These models are promising, provided the assumption of universality of structural classes is a reasonable one. However, such replacement matrices cannot be fitted to the actual sequence family under analysis. A second strategy is to develop models that depend on a reduced number of parameters (Halpern and Bruno 1998 ; Koshi and Goldstein 1998 ). These models have the advantage that they can, in principle, be fitted to the analyzed family. On the other hand, it is possible that the high reduction in the number of parameters oversimplifies such models to be generally applicable. Here we report a third strategy, which allows the determination of site-specific replacement matrices adapted to the sequence family under study, using just one member of the known structure. The lack-of-data problem is overcome by using large amounts of simulated sequence data generated using a structurally constrained protein evolution (SCPE) model developed recently (Parisi and Echave 2001 ). We apply the method to a case taken as an example and compare it with other models of protein evolution. Probabilistic models can be compared using likelihood-based statistical tests (Goldman 1993 ; Posada 2001 ; Whelan, Lio, and Goldman 2001 ). Let Pr(s|T,Q) be the probability of observing a certain sequence data set s = {s1,s2,…,sN}, given a tree T and a model Q. The maximum likelihood of the model, given the data is defined as the result of maximizing such a probability with respect to T (tree topology and branch lengths) and Q (free model parameters): L = maxT maxQPr(s|T,Q). Apart from the maximum likelihood, for model comparison, it is important to consider the number of free parameters. In this work, the relative support given by the data to a model is measured by the Bayesian information criterion (BIC), specified by BIC = − 2 ln L + npln ns, where np is the number of free parameters of the model evaluated and ns is the sequence length (Schwarz 1974 ; Posada 2001 ). The smaller the BIC, the better the fit of the model to the data. As stated previously, our purpose is to obtain evolutionary replacement matrices using the SCPE model. This model is based on introducing random mutations at the gene level and selecting sequences against too much structural variation. Without making a priori assumptions regarding site heterogeneity in substitution patterns, SCPE naturally predicts site-specific amino acid frequency distributions (Parisi and Echave 2001 ). Not being a site-independent model, SCPE cannot be used for maximum likelihood phylogenetic inference. However, it can be employed to generate enough simulated sequence data to determine site-specific substitution matrices. These matrices, in turn, can be used to build an independent-sites model of evolution, that we shall call independent-sites–structurally constrained protein evolution (IS-SCPE). The set of Qc for all site classes constitutes the IS-SCPE model. As an illustrative example, IS-SCPE is applied to the left-handed parallel β helix (LβH) domain of UDP-N-acetylglucosamine acyltransferases (LpxA). From a search of SWISS-PROT (Bairoch and Apweiler 2000 ) 25 LpxA sequences were found: LPXA_AQUAE, LPXA_BURAB, LPXA_CAMJE, LPXA_CHLMU, LPXA_CHLPN, LPXA_CHLTR, LPXA_CHRVI, LPXA_CYACA, LPXA_ECOLI, LPXA_HAEIN, LPXA_HELPJ, LPXA_HELPY, LPXA_NEIMA, LPXA_NEIMB, LPXA_PASMU, LPXA_PROMI, LPXA_PSEAE, LPXA_RICPR, LPXA_RICRI, LPXA_SALTY, LPXA_SYNY3, LPXA_VIBCH, LPXA_XYLFA, LPXA_YEREN, and Q9AQK4. The structure of the LpxA of Escherichia coli is known (Raetz and Roderick 1995 ) (PDB code 1lxa). It shows a left-handed β helix (LβH) domain (sites 1–177) that displays a characteristic hexapeptide motif (Vaara 1992 ; Vuorio et al. 1994 ; Raetz and Roderick 1995 ). Because of this, most sites can be classified into one of the six site classes, according to their position in the hexapeptides. These classes are designated i − 2, i − 1, i, i + 1, i + 2, i + 3. We aligned the 25 LpxA sequences using CLUSTAL W (Thompson, Higgins, and Gibson 1994 ). The columns of the multiple alignment that correspond to the active site (Wyckoff and Raetz 1999 ), loops, and the first and last incomplete coils of the LβH fold were not further considered because they are subject to different constraints than the rest of the hexapeptides. Each of the remaining 116 sites was classified into one of the six hexapeptide classes. With this alignment, we calculated a standard JTT distance matrix (Jones, Taylor, and Thornton 1992 ) and obtained a phylogenetic tree, using module FITCH (Fitch and Margoliash 1967 ) of PHYLIP 3.57c (Felsenstein 1993 ). This tree was used for model comparison, taking advantage of the observation that fixing the tree does not affect the model selection procedure (Posada 2001 ). Such model comparison is shown in table 1 , discussed subsequently. Robustness with respect to tree topology was verified by repeating the calculations on two other trees, obtained using the Neighbor-Joining method and the average linking clustering method, respectively, as implemented in the NEIGHBOR module of PHYLIP 3.57c (data not shown) (Felsenstein 1993 ). All maximum likelihood calculations were performed with PAML (Yang 1997 ). Using LPXA_ECOLI as the starting point, we performed 50 independent SCPE runs of 15,000 mutational steps each. Then we obtained the IS-SCPE replacement matrices for each of the six hexapeptide structural classes, as described previously. These matrices are shown in figure 1 . It is clear from the figure that, as expected, these matrices are strongly site specific. Note that they are particularly sparse for site classes i − 2 and i, which are the most conserved sites that characterize the hexapeptide motif (Raetz and Roderick 1995 ). It is of interest at this point to see whether the IS-SCPE model is able to reproduce reasonably well the evolution of the LpxA family. To do this, we performed a model comparison, shown in table 1 . The second column of this table lists the models compared. JTT, JTT+Γ, JTT+F, and JTT+F+Γ are global models based on using the standard JTT replacement matrix (Jones, Taylor, and Thornton 1992 ) for all site classes. +F indicates that the replacement matrix is modified such that it has equilibrium frequencies in agreement with the analyzed data (19 free parameters: 20 frequencies that must sum up to 1) (Lio and Goldman 1998 ). +Γ points out that a gamma distribution of replacement rates is used to describe site heterogeneity in replacement rates (one extra parameter needed to define the gamma distribution) (Yang 1993 ). The next two models use the site-specific IS-SCPE replacement matrices. We take into account the possibly different overall relative rates of the six classes (designated by +Rc in table 1 . Five free parameters: six rates minus an arbitrary time unit). Also, we consider rate heterogeneity within each class by using a class-dependent gamma distribution (one extra free parameter per class: +Γc). Finally, the last two models of table 1 , JTT+Fc, use site-specific JTT matrices, obtained by adapting JTT to each site class using the class-specific amino acid distributions estimated from the analyzed data. As for IS-SCPE, we consider interclass rate variation (+Rc) and intraclass rate variation (+Γc). The last two columns of table 1 show the logarithm of the maximum likelihood and the BIC values for the different models compared. First of all we note, from comparing the column of BIC values, that IS-SCPE fits the observed sequences better than any of the other models considered. The best fit is obtained with IS-SCPE+Rc+Γc: the model that consists of using site-specific replacement matrices obtained from SCPE simulations as described in this work and taking into account interclass (+Rc) and intraclass (+Γc) rate heterogeneity. IS-SCPE gives better results than the global JTT because it takes into account site-specific replacement patterns. Regarding the site-specific JTT models (JTT+Fc models), the maximum likelihood values seem to show that they are better than the global JTT models (see column 4 of table 1 ). However, the BIC values (column 5 of table 1 ) demonstrate that two global models—JTT+Γ and JTT+Γ+F—provide a better fit than their site-specific counterparts JTT+Fc and JTT+Fc+Γc. IS-SCPE is better than global JTT, whereas JTT+Fc is worse because JTT+Fc has 19 × 6 = 114 more parameters than IS-SCPE: the class-specific equilibrium frequencies (notice the nplnns term of the BIC statistic). In general, JTT+Fc (or any model that treats class frequencies as free parameters) will have 19 × nclasses more parameters than its IS-SCPE counterpart. Therefore, for a general case of a protein that is not as regular as LpxA, for which a larger number of classes may be needed, the gap between IS-SCPE and JTT+Fc models will rapidly increase in favor of IS-SCPE. Moreover, to mimic the site-specific substitution pattern, the JTT+Fc models depend on a reliable estimation of the 19 independent amino acid frequencies for each site class. Therefore, one faces again the lack-of-data problem: such models will be applicable only to large families or a small number of structural classes (or both) to which most sites belong, such as the present case. In contrast, the IS-SCPE replacement matrices are obtained using only one reference protein, being thus applicable independent of the size of the protein family under study or the number of site classes needed. To summarize, we have shown that a whole set of site-dependent evolutionary amino acid replacement matrices can be obtained using just one member of known structure of the protein family under study. The method overcomes the lack-of-data problem by using large sequence data sets simulated using the SCPE model of protein evolution. The replacement matrices obtained constitute an IS-SCPE probabilistic model of the evolution of the family under study. As an illustration, we apply the IS-SCPE model to the LβH domain of the LpxA family, for which we show that it performs better than JTT, or even JTT corrected to take into account site-specific amino acid distributions. William Taylor, Reviewing Editor Keywords: molecular evolution protein evolution simulation model substitution matrix Abbreviations: SCPE, structurally constrained protein evolution; IS-SCPE, independent-sites–structurally constrained protein evolution. Address for correspondence and reprints: Julian Echave, Universidad Nacional de Quilmes, Saenz Peña 180, B1876BXD Bernal, Argentina. je@unq.edu.ar . Table 1 Model Comparison: IS-SCPE Versus JTT Table 1 Model Comparison: IS-SCPE Versus JTT Fig. 1.—Site-specific replacement matrices obtained using SCPE simulations for the six hexapeptide structural classes of LpxA. Classes are identified by their position in the hexapeptide motif relative to site i of each hexapeptide: i − 2, i − 1, i, i + 1, i + 2, i + 3. This designation is indicated at the top, left corner of each panel We wish to acknowledge fruitful discussions with Pietro Lio and Nick Goldman, with whom we are currently working on a rigorous and detailed statistical assessment of IS-SCPE models using likelihood-ratio tests. This work was supported by the Universidad Nacional de Quilmes and the Fundación Antorchas. J.E. is a researcher of CONICET and a Guggenheim Fellow.