A Genetic Algorithm Approach to Detecting Lineage-Specific Variation in Selection Pressure
http://www.100md.com
《分子生物学进展》
Antiviral Research Center, University of California San Diego
Correspondence: E-mail: spond@ucsd.edu.
Abstract
The ratio of nonsynonymous (dN) to synonymous (dS) substitution rates, , provides a measure of selection at the protein level. Models have been developed that allow to vary among lineages. However, these models require the lineages in which differential selection has acted to be specified a priori. We propose a genetic algorithm approach to assign lineages in a phylogeny to a fixed number of different classes of , thus allowing variable selection pressure without a priori specification of particular lineages. This approach can identify models with a better fit than a single-ratio model, and with fits that are better than (in an information theoretic sense) a fully local model, in which all lineages are assumed to evolve under different values of , but with far fewer parameters. By averaging over models which explain the data reasonably well, we can assess the robustness of our conclusions to uncertainty in model estimation. Our approach can also be used to compare results from models in which branch classes are specified a priori with a wide range of credible models. We illustrate our methods on primate lysozyme sequences and compare them with previous methods applied to the same data sets.
Key Words: lineage-specific selection ? model selection ? genetic algorithms ? branch site model ? codon substitution model
Introduction
A common approach to characterizing selection pressure at the protein level is to compare the nonsynonymous substitution rate (dN) with the synonymous substitution rate (dS). Goldman and Yang (1994) and Muse and Gaut (1994) proposed the estimation of the ratio of the nonsynonymous to synonymous substitution rates, , in a maximum likelihood context, using a model of codon substitution. The assumption that is constant throughout the evolutionary history of a sample of sequences may be an unrealistic assumption. Yang (1998) considered a class of models where lineages specified a priori, were assumed to evolve under a different from the other lineages. These models were subsequently extended to consider both site-specific and lineage-specific variation in (Yang and Nielsen 2002).
One disadvantage of these lineage-specific models of substitution is that particular lineages evolving under different models must be specified a priori. Given the large number of possible lineage-to-model assignments, there may exist models that fit the data much better than the a priori hypothesis, and that may even contradict the biological conclusions drawn from the results based on the a priori hypothesis. Hence, it is important to place a pre-specified model in context. Even when there is no underlying biological hypothesis, it is likely that the assumption of a constant selection pressure over the entire phylogenetic tree is unrealistic. However, the large number of possible lineage-specific models make it infeasible to test the fit of every model. In this study, we propose the use of a genetic algorithm approach to assign different classes of to lineages. Genetic algorithms have proven very effective at solving very complex and poorly understood optimization problems by rapid adaptive exploration of the parameter space (Whitley 2001).
Materials and Methods
Branch-Specific Models of Substitution
We use a modification of the MG94 (Muse and Gaut 1994) model to allow for transversion/transition substitution bias and varying nucleotide frequencies between positions in a codon. The general rate matrix element for this model defines the instantaneous rate of substituting a non-stop codon x with a non-stop codon y:
(1)
s and ?s denote instantaneous synonymous and nonsynonymous substitution rates at site s, R is the ratio of transversion/transition substitution rates (we use the inverse of a more common -transition/transversion bias, because parameters with values <1 are easier to estimate numerically) and refers to the frequency of the target nucleotide in codon y in the appropriate codon position (for instance, the target nucleotide in a ACG ACT substitution is T in position 3). While the model does not allow for multiple instantaneous substitutions, because such events are assumed negligibly rare, the transition matrix T(t) = exp(t x MG94) for any t > 0 allows for positive transition probabilities between any two states. The model also provides the equilibrium frequency of a codon composed of the nucleotide triplet ijk. If we denote the frequency of nucleotide n {A, C, G, T} at codon position m = 1, 2, 3 as then the equilibrium frequency of codon ijk is the product of the constituent nucleotide frequencies, scaled to account for the absence of stop codons (TAA, TAG, TGA for the universal genetic code) in the model:
(2)
Consequently, there are nine independent frequency parameters in this model. In practice, they are estimated by the proportions observed in the data, but they could also be inferred by maximum likelihood. This model can be simplified if we assume that nonsynonymous and synonymous substitution rates are uniform across the sequence; i.e., s = = 1 and ?s = ?, such that the ratio of nonsynonymous to synonymous substitution rates, , equals ?. We can extend this model to consider variable nonsynonymous substitution rates in different lineages, by replacing ?s with ?sb, where b refers to a specific lineage. We note that without assuming a molecular clock, we cannot identify variable synonymous substitution rates in different lineages.
This model differs slightly from the GY94 (Goldman and Yang 1994) model used in Yang (1998), in that GY94 uses y in place of in the rate matrix. In our experience, the MG94 x HKY85 yields better likelihood scores than GY94 with the same number of parameters, especially when base composition of the sequences is biased. For comparative purposes, we also included fits for most of the models using the GY94 rate matrix parameterization.
Model Selection Using a Genetic Algorithm
Given n sequences, an unrooted phylogeny with 2n – 3 branches, and the number of desired lineage classes B, we search the state space of lineage-to-class assignments where bi 0 ... B – 1. The ith entry of the vector represents the class assignment for the ith branch in the post-order traversal of the tree. In order not to include equivalent branch partitions more than once (e.g., (0,1,2,1,2) and (0,2,1,2,1) for a 5-branch tree), we require that the first occurrence of each rate assignment in the vector are in increasing order. There are S(2n – 3, k) possible assignments to consider, and if a model with fewer than B classes has a better fitness, it will be included in the search space. Here
denotes the Stirling number of the second kind (Abramowitz and Stegun 1972), which enumerates the number of ways to partition N objects into k non-empty sets. This search space is far too large to allow an exhaustive search, even for a small alignment. For instance, assuming 3 branch classes on a seven sequence alignment yields 29,525 models to consider, and allowing for 4 or fewer branch classes on 19 sequences results in 8.34 x 1015 possibilities.
We employ the CHC—Cross generational elitist selection, Heterogeneous recombination, and Cataclysmic mutation (Eshelman 1991)—an aggressive population-based hill-climber. This algorithm always retains the most fit individual. It uses a radical recombination operator to rapidly explore the parameter space without collapsing it. When the diversity of the population drops below a certain threshold, the algorithm invokes a radical mutation step: a large proportion (15%–35%) of randomly chosen positions in all individuals, except for the most fit individual, are mutated. The mutation step helps to prevent entrapment in local maxima.
The recombination operator of the CHC algorithm generates an offspring from two parent states by randomly choosing one of the two parental values in every position of the state vector which differs between the parents.
The parameter space for our optimization problem is composed of two components: discrete partitioning of tree branches into different classes and a vector of real valued parameters, such as branch lengths, b and R. We used the CHC algorithm to search through the discrete component of the parameter space and conventional numerical optimization techniques to find maximum likelihood estimates of all other model parameters, given a partitioning of branches into rate classes. The fitness of every model was measured by its small sample Akaike's Information Criterion score (AICc; Sugiura 1978; Akaike 1974). The AICc of a model with p parameters fitted to a sample of size s is defined as
For any fixed p, and thus AICc and the standard Akaike Information Criterion coincide for large sample sizes. Clearly, AICc can only be used when s > p + 1. In our case we treat each alignment column as an independent sample, and the method requires at least 2n – 1 + B alignment columns to be applicable. AICc is a second order correction to standard AIC, and its use can be advocated unless the number of samples is much larger (40x or more) than the number of parameters (Burnham and Anderson [2003], p. 66).
The mutation phase is triggered if the difference between the AICc of the best fitting model and the worst fitting model in a generation decreases below a threshold > 0; on average, each 5th position in all model vectors of the current generation (excluding the best fitting one) is randomly replaced with an integer in 0 ... B – 1 (excluding the present state).
Our parallel implementation of the CHC algorithm was run on a 16-node Linux MPI cluster, with each generation consisting of 28 individuals. We used the convergence criterion that requires the AICc of the best fitting individual over all generations to remain constant for 30 generations. We also verified convergence by repeating the runs several times, with the initial population generated randomly every time.
Once we have selected a model using the GA, it is important to assess the confidence we have in the model, both in terms of how much better it fits compared to other candidate models, and in terms of whether we would select the same model again if other independent data sets were available. The goodness-of-fit between the best fitting model identified using the GA, the global (single-ratio) model, and the local (free-ratio) model, can be determined using a likelihood ratio test, as these models are nested. To determine the significance of the difference in fit between models specified a priori, and the best fitting model identified using a GA (which may not necessarily be nested), we employ a Shimodaira-Hasegawa test (Shimodaira and Hasegawa 1999), using the difference in AICc, rather than the difference in log-likelihoods, to compare models with different numbers of parameters.
Rather than just base inference on the single best fitting model, we also consider averaging over a set of models. We calculate the Akaike weights (Akaike 1983) for each model examined during a run of the GA. If we define then the Akaike weight for model mi, where M is the total number of models, is given by the following.
(3)
Model-averaged estimates of a parameter can be obtained using these Akaike weights. In addition, by summing over the Akaike weights from largest to smallest until the sum is just over 0.95, we can rapidly obtain a 95% confidence set of models on which to base inference. In general, it is very difficult to ascertain that any given GA run (or series of GA runs) achieves good coverage of the confidence set. Theoretical results (Whitley 2001) suggest that genetic algorithms sample regions of the search space with probability proportional to their fitness (at least asymptotically), so that highly likely models are going to be sampled with relatively high probability. In addition, model averaged estimates were monitored for apparent convergence, and they were robust to combining models from multiple GA runs. We can also assess the relative likelihoods of model i versus model j simply as the ratio of their Akaike weights, which are called evidence ratios.
Effect of Site-to-Site Rate Variation
Our GA algorithms can be easily adapted to utilize models with site-to-site rate variation while searching for the best branch class assignments. However, using such models will significantly increase run times of individual model fitting and the GA overall. We fitted models with two types of site-to-site rate variation using branch partitioning obtained with the homogeneous site rates model to investigate relative improvements of fit.
First, we consider the scenario where mutational loads vary from site to site. Formally, this scenario can be implemented by using the following rate matrix
(4)
where μs is sampled from a unit mean general discrete distribution (UGDD). Given D rate classes to include in a UGDD, the distribution is completely defined by specifying its values 0 r1 < r2 < ... < rD and weights pi = Pr{μ = ri}, i = 1...D, subject to The unit mean constraint is due to a standard identifiability issue: only products of rates and branch lengths can be estimated.
Second, we test the hypothesis that selective pressure varies across alignment sites. Formally,
(5)
where μs is sampled from a UGDD. The unit mean restriction still applies, because only products of μs and ?sb can be estimated. For the data sets in this study, UGDD with 3 bins and 4 distributional parameters were estimated using maximum likelihood. Even though UGDD has more parameters than a more traditionally used discretized-gamma distribution, it is a much more flexible distribution able, for example, to model invariable sites and bi-modal distributions. Similar distributions have been employed previously, for example in the context of identifying individual amino acid sites under selection (M3 in Yang et al. [2000]).
Note that the MG94 x HKY85 model is nested in both μMG94 x HKY85 and MG94 x HKY85. These models are somewhat similar to the branch-site models (Yang and Nielsen 2002). However, whereas Nielsen and Yang were interested in gaining more power to detect individual sites under selection, we are primarily interested in whether incorporating site-to-site rate variation can significantly alter the improvements in fit obtained by using branch-specific models. Essentially, by asking the question of when selection occurred we cannot ignore the effect that where the selection might have occurred could have on our conclusions, because the two effects may be confounded. Additionally, our branch-specific b serve as rate multipliers at every site in the alignment, where as the branch-site models allow for variable b along pre-specified branches and among sites.
Sequence Data and an A Priori Model of Sequence Evolution
We investigate the primate lysozyme data originally analyzed by Messier and Stewart (1997), and subsequently by Yang (1998) and Yang and Nielsen (2002). Following Yang (1998), we consider a large data set consisting of 19 distinct sequences and a small 7-sequence data set, representing the four major groups (Colobines, Cercopithecines, Hominoids, and New World monkeys) in the phylogeny. As an example of an a priori biological model, we used model E of Yang (1998), which assumed that the branches leading to hominoids and colobines each had an independently estimated value of (c and h), whereas all other branches shared a third value = 0, which was also estimated. Apart from the free-ratio model, model E was the most general of those considered in Yang (1998). We denote this model as Y98.
Implementation
Recent advances in computational power make the use of complex search algorithms eminently feasible. All sequence analyses and model fitting were performed using the HyPhy (Kosakovsky Pond, Frost, and Muse 2000) software package on an 16-node Linux MPI cluster; 14 slave nodes were used to fit various models, and a single master node dispatched the jobs and assembled the results. A single pass of the GA algorithm required from an hour for the small data set to a day for the large data set. Lysozyme alignments used in this article and all HyPhy batch files needed to perform the analyses and result processing can be obtained from http://www.hyphy.org/gabranch/.
Results
Messier and Stewart (1997) suggested that positive selection on lysozyme occurred along the lineage leading to colobine monkeys and the lineage leading to hominoids (see fig. 2), followed by a period of negative selection within each clade. Although Yang (1998) found evidence of an elevated along the lineages leading to colobines and hominoids, as only a limited number of models were tested, it is not clear whether there is a model that fits the data better, and if such a model exists, whether this model supports or refutes the hypothesis of elevated in these lineages. To address this issue, we applied a genetic algorithm approach, which assigns rate classes of to lineages according to their goodness-of-fit (in terms of small sample AICc) to the small and large lysozyme data sets considered by Yang (1998). A representative plot of the AICc of the population, where individuals represent different models over time, is shown in figure 1. This figure illustrates the decrease in AICc over subsequent generations as well as hypermutation events, in which models are randomly mutated when the difference in AICc between the best fitting (lowest AICc) and worst fitting (highest AICc) models in the population falls below a threshold (= 1 in our runs).
FIG. 2.— Branches of the tree estimated using the small and large lysozyme data sets partitioned according to the best fitting GA model under MG94 x HKY85. Branch labels represent model averaged probabilities of dN/dS > 1 for the branch. Percentages for branch classes in the legend reflect the proportion of total tree length (measured in expected substitutions per site per unit time) evolving under the corresponding value of dN/dS. Branches in this tree are unscaled.
FIG. 1.— A run of the branch class (3 classes) assignment genetic algorithm on the large lysozyme data set, using the MG94 x HKY85 model. The solid black line represents the score of the best fitting (in the AICc sense) model, and the gray plot shows the AICc score of every model in the order in which they were tried.
The AICc of the best fitting model assuming different numbers of branch classes are shown in table 1 for the small and large lysozyme data sets. Note that in both cases, moving from 3 to 4 branch classes did not result in a better fit as the models were collapsed to 3 rate classes during the search. The best fitting model identified using the GA can be compared with the single rate model and with the free ratio model by means of a likelihood ratio test. Because of the large number of possible models, stochasticity of observed data may lead to the best fitting model being chosen for spurious reasons. It is currently computationally infeasible to perform a full parametric or nonparametric bootstrap analysis, in which replicates of simulated data are subjected to analysis using the genetic algorithm. However, the robustness of the likelihood ratio test statistic to sampling error can be assessed using the approximate method of Shimodaira and Hasegawa (1999). This test can also be applied to comparison of the models considered by Yang (1998) and the best fitting models.
Table 1 Best AICc Found and Number of Models Examined for a Given Number of Branch Rate Classes under the MG94 x HKY85 Model
For the small data set (table 2), the best fitting model (we used the MG94 x HKY85 substitution model throughout) was significantly better than the single rate model (P = 0.0005), and not significantly different from the free ratio model (P = 0.93). The best fitting model had a lower AICc than the Y98 model—1824.39 versus 1830.72—and this difference was statistically significant as assessed by a Shimodaira-Hasegawa test (P < 0.007 using 10,000 replicates); 21 sites fitted better under the Yang (1998) model, and 109 sites fitted better under the GA model. We also note that these models are not nested.
Table 2 Model Fitting Results for the Small (7 sequences) Lysozyme Data Set
For the large data set (table 3), the best fitting model was also significantly better than the single ratio model (P = 3 x 10–8), and not significantly different from the free ratio model (P = 0.999991). The best fitting model had a significantly lower AICc than the Y98 model, with Shimodaira-Hasegawa P = 0.0156 (based on 10,000 replicates), and the models were not nested.
Table 3 Model Fitting Results for the Large (19 sequences) Lysozyme Data Set
For both the small and the large data sets (tables 2 and 3), the median dN/dS = ?sb for the Y98 model was much higher than than the median ?sb for the GA model. Both model misspecification and overfitting appear to contribute to this disparity. First, the median ?sb for the GA model is close to that of the free ratio model, whereas the median ?sb for the Y98 model is much higher than that of the free ratio model; this suggests that the Y98 model is failing to capture the pattern of variation in ?sb in the data. In addition, both the Y98 and the free ratio model have a class for which ?sb = ; clearly, this is biologically unrealistic and is the consequence of estimating ?sb from branches with no synonymous changes. In contrast, the GA model combines these branches into a single class with a high, but finite, value of ?sb.
There were a large number of models in the 95% confidence set for both the small (211) and large data sets (333). The evidence ratio for the best fitting model and the second best fitting model was 1.07 for the small data set. The two top models differ in class allocation for a single branch—the ancestor of colobines and cercopithecines (labeled as N3 in figure 2). The best fitting model allocated this branch to a class with ?sb = 0.61, while the second best model assigned it to the class with ?sb = 4.39. The third best model had the evidence ratio of 1.92. We can also compute the largest AICc needed for a model to be included in the 95% confidence set. For the smaller data set, this value was 1832.76, and hence the Y98 model is included in the credible set.
For the large data set, there were 19 models with essentially equal support (evidence ratio 1). This is likely due to the relatively few (4) sites per estimated model parameter and concomitant unreliable inference, as well as the the presence of short branches in the phylogeny. The threshold AICc value was found to be 2152.28, and thus the a priori model with AICc of 2167.54 is not included in the confidence set of models.
To assess the robustness of the conclusions drawn from the best fitting model, we calculated the support for an estimate of dN/dS = ?sb > 1 for each branch (fig. 2), averaging over models using their Akaike weights. The support for ?sb greater than 1 in the lineage leading to hominoids and colobines was 94.8% and 88.4%, respectively, for the small data set and 99.9% and 99.6%, respectively, for the large data set.
We next fitted site-to-site rate variation models using the branch assignments derived from running the genetic algorithm using a homogeneous rates model. Adding branch variation to a site-to-site model and site-to-site variation to a branch model yielded substantial AICc improvements for both data sets. Also, AICc scores for μMG94 x HKY85 and MG94 x HKY85 models indicate that variation in site-to-site mutational load fits lysozyme data better than variation in selection pressure over sites. Interestingly, the improvement in goodness-of-fit (in terms of the log likelihood) obtained by allowing both site-to-site variation in mutational load and branch variation was close to the improvement in fit when considering each effect separately. For the large data set, the difference in log-likelihood, logL, between the one ratio model and the GA 3-ratio model with variation in mutational load was 32.15, which is close to the sum of logL when considering branch variation and site-to-site variation separately (17.45 + 14.6 = 32.05). In addition, the estimates for ?sb were similar between models with and without site-to-site variation in mutational load. Thus, we can conclude that branch-specific variation in selection pressure and site-specific variation in mutational load play important (and apparently orthogonal) roles in shaping genetic variation in lysozyme. We also ran the GA using the μMG94 x HKY85 model on the small data set and found that the best fitting branch classes were the same as the ones found by the GA under the MG94 x HKY85 model, consistent with the apparent orthogonality between the models of site-to-site variation in mutational load and branch-specific variation in ?sb.
The model-averaging approach enables us to address issues which cannot be easily resolved using simple hypothesis testing. For example, we could compute probabilities of any two branches evolving under the same selective pressure (table 4), without making a priori assumptions about selective pressure acting on all other branches.
Table 4 Model Averaged Probabilities That Two Branches of the Small Lysozyme Phylogeny Have the Same ?sb = dN/dS Ratio
Discussion
Genetic algorithms have been employed in multiple studies to search for phylogenetic trees (Lewis 1998; Katoh, Kuma, and Miyata 2001; Brauer et al. 2002; Lemmon and Milinkovitch 2002; Kim, Lee, and Moon 2003; Shen and Heckendorn 2004). Our study is the first to our knowledge that uses a genetic algorithm to search for models of substitution in a phylogenetic context. Classical GA algorithms employ large population sizes (hundreds or thousands of individuals), and implicitly rely on the ability to compute the fitness of a single individual (in this case, the AICc of a given model) very quickly. Unfortunately, most molecular evolution problems, such as those considered here, are computationally challenging, and the number of fitness evaluations must be limited. Fortunately, there exist aggressive versions of GA, such as the CHC algorithm used in this study, which work well with smaller population sizes. Alternatively, a reversible jump Markov chain Monte Carlo procedure could be used to simultaneously estimate the number of classes of and the assignments to branches; however, this approach is likely to be far more computationally intensive than the GA approach pursued here. The GA method can be run on a small cluster of computers in a matter of hours and is therefore feasible for immediate use by researchers.
Using Akaike's Information Criterion as a measure of goodness-of-fit has been criticized for being too liberal, and may select models that are overly complex; instead, we have used a small sample AIC, AICc . In principle, information criteria other than AICc could be used, for example Schwarz's Bayesian Information Criterion (BIC) (Schwarz 1978), which is consistent (the probability that it will recover a true low-dimensional model approaches 1 as the sample size tends to infinity). However, the BIC assumes that the true model is contained within the candidate set of models, whereas the AIC does not assume that any of the candidate models is necessarily true; rather, it quantifies the discrepancy between the probability density generated by the model and the data as measured by Kullback-Liebler information (Kullback and Liebler 1951).
Given the large number of possible models and limited data, it is likely that a substantial number—tens or hundreds—of models will be consistent with the data. With the use of Akaike weights, which can be directly interpreted as conditional probabilities of the models, the robustness of conclusions to the choice of model can be assessed by averaging over models. Furthermore, an a priori biological model can be placed in the context of other credible models. If, for example, such a model falls outside the inferred 95% confidence set, it may not be well supported, as there are many more models which fit the data significantly better, even when the a priori model may be significantly better than the single-ratio model, as assessed by a nested comparison using a likelihood ratio test. Additionally we can test whether an a priori biological model is significantly worse than the best fitting model found by the GA using a Shimodaira-Hasegawa test. It is important to consider whether the a priori model is simply a submodel of the best fitting model, and whether the biological conclusions are altered under the best fitting model. Based on the small lysozyme data set, we find that the hypothesis that the lineages leading to Colobines and Hominoids have elevated nonsynonymous to synonymous substitution rates is consistent with our GA-based analysis. However, our GA results based on the large data set suggest that positive selection is ongoing following the radiation of these clades in contrast to the hypothesis of Messier and Stewart (1997).
When an a priori model compares favorably with the set of "good" models, we may be more confident in biological conclusions derived from studying such an a priori model. If, however, there are many alternative models which fit the data much better than our preconceived hypothesis, further exploration of the data, and perhaps a reassessment of our assumptions, is in order. In many practical cases the model itself is a nuisance component of the analysis, for example when one is interested in evolution along a specific lineage, and the rest of the phylogeny forms the "background." By weighing over credible models, quantities of biological interest (such as the probability that along a particular branch dN/dS > 1) may be derived without relying on a specific a priori model and concomitant assumptions, such as uniformity of selective pressure along background branches.
Acknowledgements
We thank Professor A.J. Leigh Brown for many helpful suggestions. This work was supported in part by the National Institutes of Health (grant numbers 5 R01 AI47745 and 5 U01 AI43638 Supp), the University of California Universitywide AIDS Research Program (grant number IS02-SD-701), and by a University of California, San Diego Center for AIDS Research/NIAID Developmental Award to S.D.W.F. (grant no. 2 P30 AI36214).
References
Abramowitz, M., I. Stegun, eds. 1972. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th printing. Dover, New York.
Akaike, H. 1974. A new look at the statistical model identification. IEEE Trans. Automatic Control 119:716–723.
Akaike, H. 1983. Information measures and model selection. Int. Stat. Inst. 44:139–149.
Brauer, M. J., M. T. Holder, L. A. Dries, D. J. Zwickl, P. O. Lewis, and D. M. Hillis. 2002. Genetic algorithms and parallel processing in maximum-likelihood phylogeny inference. Mol. Biol. Evol. 19:1717–1726.
Burnham, K. P., and D. R. Anderson. 2003. Model Selection and Multimodel Inference, 2nd ed. Springer, New York.
Eshelman, L. 1991. Foga-1. Morgan Kaufmann, Los Atlos, Calif.
Goldman, N., and Z. H. Yang. 1994. Codon-based model of nucleotide substitution for protein-coding DNA-sequences. Mol. Biol. Evol. 11:725–736.
Katoh, K., K. Kuma, and T. Miyata. 2001. Genetic algorithm-based maximum-likelihood analysis for molecular phylogeny. J. Mol. Evol. 53:477–484.
Kim, Y. H., S. K. Lee, and B. R. Moon. 2003. Optimizing the order of taxon addition in phylogenetic tree construction using genetic algorithm. Genetic and Evolutionary Computation—Gecco 2003, Pt Ii, Proceedings 2724:2168–2178.
Kosakovsky Pond, S., S. D. W. Frost, and S. V. Muse. 2000. HyPhy: hypothesis testing using phylogenies. Bioinformatics Advance Access published on October 27, 2004, doi:10.1093/bioinformatics/bti079.
Kullback, S., and R. Liebler. 1951. On information and sufficiency. Ann. Math. Stat. 22:79–86.
Lemmon, A. R., and M. C. Milinkovitch. 2002. The metapopulation genetic algorithm: an efficient solution for the problem of phylogeny estimation. Proc. Natl. Acad. Sci. USA 99:10516–10521.
Lewis, P. O. 1998. A genetic algorithm for maximum-likelihood phylogeny inference using nucleotide sequence data. Mol. Biol. Evol. 15:277–283.
Messier, W., and C. B. Stewart. 1997. Episodic adaptive evolution of primate lysozymes. Nature 385:151–154.
Muse, S. V., and B. S. Gaut. 1994. A likelihood approach for comparing synonymous and nonsynonymous nucleotide substitution rates, with application to the chloroplast genome. Mol. Biol. Evol. 11:715–724.
Schwarz, G. 1978. Estimating the dimension of a model. Ann. Stat. 6:461–464.
Shen, J., and R. B. Heckendorn. 2004. Discrete branch length representations for genetic algorithms in phylogenetic search. Appl. Evol. Comput. 3005:94–103.
Shimodaira, H., and M. Hasegawa. 1999. Multiple comparisons of log-likelihoods with applications to phylogenetic inference. Mol. Biol. Evol. 16:1114–1116.
Sugiura, N. 1978. Further analysis of the data by Akaike's information criterion and the finite corrections. Commun. Stat. Theory Methods A7:13–26.
Whitley, D. 2001. An overview of evolutionary algorithms: practical issues and common pitfalls. J. Inf. Software Technol. 43:817–831.
Yang, Z. H. 1998. Likelihood ratio tests for detecting positive selection and application to primate lysozyme evolution. Mol. Biol. Evol. 15:568–573.
Yang, Z. H., and R. Nielsen. 2002. Codon-substitution models for detecting molecular adaptation at individual sites along specific lineages. Mol. Biol. Evol. 19:908–917.
Yang, Z. H., R. Nielsen, N. Goldman, and A. M. K. Pedersen. 2000. Codon-substitution models for heterogeneous selection pressure at amino acid sites. Genetics 155:431–449.(Sergei L. Kosakovsky Pond)
Correspondence: E-mail: spond@ucsd.edu.
Abstract
The ratio of nonsynonymous (dN) to synonymous (dS) substitution rates, , provides a measure of selection at the protein level. Models have been developed that allow to vary among lineages. However, these models require the lineages in which differential selection has acted to be specified a priori. We propose a genetic algorithm approach to assign lineages in a phylogeny to a fixed number of different classes of , thus allowing variable selection pressure without a priori specification of particular lineages. This approach can identify models with a better fit than a single-ratio model, and with fits that are better than (in an information theoretic sense) a fully local model, in which all lineages are assumed to evolve under different values of , but with far fewer parameters. By averaging over models which explain the data reasonably well, we can assess the robustness of our conclusions to uncertainty in model estimation. Our approach can also be used to compare results from models in which branch classes are specified a priori with a wide range of credible models. We illustrate our methods on primate lysozyme sequences and compare them with previous methods applied to the same data sets.
Key Words: lineage-specific selection ? model selection ? genetic algorithms ? branch site model ? codon substitution model
Introduction
A common approach to characterizing selection pressure at the protein level is to compare the nonsynonymous substitution rate (dN) with the synonymous substitution rate (dS). Goldman and Yang (1994) and Muse and Gaut (1994) proposed the estimation of the ratio of the nonsynonymous to synonymous substitution rates, , in a maximum likelihood context, using a model of codon substitution. The assumption that is constant throughout the evolutionary history of a sample of sequences may be an unrealistic assumption. Yang (1998) considered a class of models where lineages specified a priori, were assumed to evolve under a different from the other lineages. These models were subsequently extended to consider both site-specific and lineage-specific variation in (Yang and Nielsen 2002).
One disadvantage of these lineage-specific models of substitution is that particular lineages evolving under different models must be specified a priori. Given the large number of possible lineage-to-model assignments, there may exist models that fit the data much better than the a priori hypothesis, and that may even contradict the biological conclusions drawn from the results based on the a priori hypothesis. Hence, it is important to place a pre-specified model in context. Even when there is no underlying biological hypothesis, it is likely that the assumption of a constant selection pressure over the entire phylogenetic tree is unrealistic. However, the large number of possible lineage-specific models make it infeasible to test the fit of every model. In this study, we propose the use of a genetic algorithm approach to assign different classes of to lineages. Genetic algorithms have proven very effective at solving very complex and poorly understood optimization problems by rapid adaptive exploration of the parameter space (Whitley 2001).
Materials and Methods
Branch-Specific Models of Substitution
We use a modification of the MG94 (Muse and Gaut 1994) model to allow for transversion/transition substitution bias and varying nucleotide frequencies between positions in a codon. The general rate matrix element for this model defines the instantaneous rate of substituting a non-stop codon x with a non-stop codon y:
(1)
s and ?s denote instantaneous synonymous and nonsynonymous substitution rates at site s, R is the ratio of transversion/transition substitution rates (we use the inverse of a more common -transition/transversion bias, because parameters with values <1 are easier to estimate numerically) and refers to the frequency of the target nucleotide in codon y in the appropriate codon position (for instance, the target nucleotide in a ACG ACT substitution is T in position 3). While the model does not allow for multiple instantaneous substitutions, because such events are assumed negligibly rare, the transition matrix T(t) = exp(t x MG94) for any t > 0 allows for positive transition probabilities between any two states. The model also provides the equilibrium frequency of a codon composed of the nucleotide triplet ijk. If we denote the frequency of nucleotide n {A, C, G, T} at codon position m = 1, 2, 3 as then the equilibrium frequency of codon ijk is the product of the constituent nucleotide frequencies, scaled to account for the absence of stop codons (TAA, TAG, TGA for the universal genetic code) in the model:
(2)
Consequently, there are nine independent frequency parameters in this model. In practice, they are estimated by the proportions observed in the data, but they could also be inferred by maximum likelihood. This model can be simplified if we assume that nonsynonymous and synonymous substitution rates are uniform across the sequence; i.e., s = = 1 and ?s = ?, such that the ratio of nonsynonymous to synonymous substitution rates, , equals ?. We can extend this model to consider variable nonsynonymous substitution rates in different lineages, by replacing ?s with ?sb, where b refers to a specific lineage. We note that without assuming a molecular clock, we cannot identify variable synonymous substitution rates in different lineages.
This model differs slightly from the GY94 (Goldman and Yang 1994) model used in Yang (1998), in that GY94 uses y in place of in the rate matrix. In our experience, the MG94 x HKY85 yields better likelihood scores than GY94 with the same number of parameters, especially when base composition of the sequences is biased. For comparative purposes, we also included fits for most of the models using the GY94 rate matrix parameterization.
Model Selection Using a Genetic Algorithm
Given n sequences, an unrooted phylogeny with 2n – 3 branches, and the number of desired lineage classes B, we search the state space of lineage-to-class assignments where bi 0 ... B – 1. The ith entry of the vector represents the class assignment for the ith branch in the post-order traversal of the tree. In order not to include equivalent branch partitions more than once (e.g., (0,1,2,1,2) and (0,2,1,2,1) for a 5-branch tree), we require that the first occurrence of each rate assignment in the vector are in increasing order. There are S(2n – 3, k) possible assignments to consider, and if a model with fewer than B classes has a better fitness, it will be included in the search space. Here
denotes the Stirling number of the second kind (Abramowitz and Stegun 1972), which enumerates the number of ways to partition N objects into k non-empty sets. This search space is far too large to allow an exhaustive search, even for a small alignment. For instance, assuming 3 branch classes on a seven sequence alignment yields 29,525 models to consider, and allowing for 4 or fewer branch classes on 19 sequences results in 8.34 x 1015 possibilities.
We employ the CHC—Cross generational elitist selection, Heterogeneous recombination, and Cataclysmic mutation (Eshelman 1991)—an aggressive population-based hill-climber. This algorithm always retains the most fit individual. It uses a radical recombination operator to rapidly explore the parameter space without collapsing it. When the diversity of the population drops below a certain threshold, the algorithm invokes a radical mutation step: a large proportion (15%–35%) of randomly chosen positions in all individuals, except for the most fit individual, are mutated. The mutation step helps to prevent entrapment in local maxima.
The recombination operator of the CHC algorithm generates an offspring from two parent states by randomly choosing one of the two parental values in every position of the state vector which differs between the parents.
The parameter space for our optimization problem is composed of two components: discrete partitioning of tree branches into different classes and a vector of real valued parameters, such as branch lengths, b and R. We used the CHC algorithm to search through the discrete component of the parameter space and conventional numerical optimization techniques to find maximum likelihood estimates of all other model parameters, given a partitioning of branches into rate classes. The fitness of every model was measured by its small sample Akaike's Information Criterion score (AICc; Sugiura 1978; Akaike 1974). The AICc of a model with p parameters fitted to a sample of size s is defined as
For any fixed p, and thus AICc and the standard Akaike Information Criterion coincide for large sample sizes. Clearly, AICc can only be used when s > p + 1. In our case we treat each alignment column as an independent sample, and the method requires at least 2n – 1 + B alignment columns to be applicable. AICc is a second order correction to standard AIC, and its use can be advocated unless the number of samples is much larger (40x or more) than the number of parameters (Burnham and Anderson [2003], p. 66).
The mutation phase is triggered if the difference between the AICc of the best fitting model and the worst fitting model in a generation decreases below a threshold > 0; on average, each 5th position in all model vectors of the current generation (excluding the best fitting one) is randomly replaced with an integer in 0 ... B – 1 (excluding the present state).
Our parallel implementation of the CHC algorithm was run on a 16-node Linux MPI cluster, with each generation consisting of 28 individuals. We used the convergence criterion that requires the AICc of the best fitting individual over all generations to remain constant for 30 generations. We also verified convergence by repeating the runs several times, with the initial population generated randomly every time.
Once we have selected a model using the GA, it is important to assess the confidence we have in the model, both in terms of how much better it fits compared to other candidate models, and in terms of whether we would select the same model again if other independent data sets were available. The goodness-of-fit between the best fitting model identified using the GA, the global (single-ratio) model, and the local (free-ratio) model, can be determined using a likelihood ratio test, as these models are nested. To determine the significance of the difference in fit between models specified a priori, and the best fitting model identified using a GA (which may not necessarily be nested), we employ a Shimodaira-Hasegawa test (Shimodaira and Hasegawa 1999), using the difference in AICc, rather than the difference in log-likelihoods, to compare models with different numbers of parameters.
Rather than just base inference on the single best fitting model, we also consider averaging over a set of models. We calculate the Akaike weights (Akaike 1983) for each model examined during a run of the GA. If we define then the Akaike weight for model mi, where M is the total number of models, is given by the following.
(3)
Model-averaged estimates of a parameter can be obtained using these Akaike weights. In addition, by summing over the Akaike weights from largest to smallest until the sum is just over 0.95, we can rapidly obtain a 95% confidence set of models on which to base inference. In general, it is very difficult to ascertain that any given GA run (or series of GA runs) achieves good coverage of the confidence set. Theoretical results (Whitley 2001) suggest that genetic algorithms sample regions of the search space with probability proportional to their fitness (at least asymptotically), so that highly likely models are going to be sampled with relatively high probability. In addition, model averaged estimates were monitored for apparent convergence, and they were robust to combining models from multiple GA runs. We can also assess the relative likelihoods of model i versus model j simply as the ratio of their Akaike weights, which are called evidence ratios.
Effect of Site-to-Site Rate Variation
Our GA algorithms can be easily adapted to utilize models with site-to-site rate variation while searching for the best branch class assignments. However, using such models will significantly increase run times of individual model fitting and the GA overall. We fitted models with two types of site-to-site rate variation using branch partitioning obtained with the homogeneous site rates model to investigate relative improvements of fit.
First, we consider the scenario where mutational loads vary from site to site. Formally, this scenario can be implemented by using the following rate matrix
(4)
where μs is sampled from a unit mean general discrete distribution (UGDD). Given D rate classes to include in a UGDD, the distribution is completely defined by specifying its values 0 r1 < r2 < ... < rD and weights pi = Pr{μ = ri}, i = 1...D, subject to The unit mean constraint is due to a standard identifiability issue: only products of rates and branch lengths can be estimated.
Second, we test the hypothesis that selective pressure varies across alignment sites. Formally,
(5)
where μs is sampled from a UGDD. The unit mean restriction still applies, because only products of μs and ?sb can be estimated. For the data sets in this study, UGDD with 3 bins and 4 distributional parameters were estimated using maximum likelihood. Even though UGDD has more parameters than a more traditionally used discretized-gamma distribution, it is a much more flexible distribution able, for example, to model invariable sites and bi-modal distributions. Similar distributions have been employed previously, for example in the context of identifying individual amino acid sites under selection (M3 in Yang et al. [2000]).
Note that the MG94 x HKY85 model is nested in both μMG94 x HKY85 and MG94 x HKY85. These models are somewhat similar to the branch-site models (Yang and Nielsen 2002). However, whereas Nielsen and Yang were interested in gaining more power to detect individual sites under selection, we are primarily interested in whether incorporating site-to-site rate variation can significantly alter the improvements in fit obtained by using branch-specific models. Essentially, by asking the question of when selection occurred we cannot ignore the effect that where the selection might have occurred could have on our conclusions, because the two effects may be confounded. Additionally, our branch-specific b serve as rate multipliers at every site in the alignment, where as the branch-site models allow for variable b along pre-specified branches and among sites.
Sequence Data and an A Priori Model of Sequence Evolution
We investigate the primate lysozyme data originally analyzed by Messier and Stewart (1997), and subsequently by Yang (1998) and Yang and Nielsen (2002). Following Yang (1998), we consider a large data set consisting of 19 distinct sequences and a small 7-sequence data set, representing the four major groups (Colobines, Cercopithecines, Hominoids, and New World monkeys) in the phylogeny. As an example of an a priori biological model, we used model E of Yang (1998), which assumed that the branches leading to hominoids and colobines each had an independently estimated value of (c and h), whereas all other branches shared a third value = 0, which was also estimated. Apart from the free-ratio model, model E was the most general of those considered in Yang (1998). We denote this model as Y98.
Implementation
Recent advances in computational power make the use of complex search algorithms eminently feasible. All sequence analyses and model fitting were performed using the HyPhy (Kosakovsky Pond, Frost, and Muse 2000) software package on an 16-node Linux MPI cluster; 14 slave nodes were used to fit various models, and a single master node dispatched the jobs and assembled the results. A single pass of the GA algorithm required from an hour for the small data set to a day for the large data set. Lysozyme alignments used in this article and all HyPhy batch files needed to perform the analyses and result processing can be obtained from http://www.hyphy.org/gabranch/.
Results
Messier and Stewart (1997) suggested that positive selection on lysozyme occurred along the lineage leading to colobine monkeys and the lineage leading to hominoids (see fig. 2), followed by a period of negative selection within each clade. Although Yang (1998) found evidence of an elevated along the lineages leading to colobines and hominoids, as only a limited number of models were tested, it is not clear whether there is a model that fits the data better, and if such a model exists, whether this model supports or refutes the hypothesis of elevated in these lineages. To address this issue, we applied a genetic algorithm approach, which assigns rate classes of to lineages according to their goodness-of-fit (in terms of small sample AICc) to the small and large lysozyme data sets considered by Yang (1998). A representative plot of the AICc of the population, where individuals represent different models over time, is shown in figure 1. This figure illustrates the decrease in AICc over subsequent generations as well as hypermutation events, in which models are randomly mutated when the difference in AICc between the best fitting (lowest AICc) and worst fitting (highest AICc) models in the population falls below a threshold (= 1 in our runs).
FIG. 2.— Branches of the tree estimated using the small and large lysozyme data sets partitioned according to the best fitting GA model under MG94 x HKY85. Branch labels represent model averaged probabilities of dN/dS > 1 for the branch. Percentages for branch classes in the legend reflect the proportion of total tree length (measured in expected substitutions per site per unit time) evolving under the corresponding value of dN/dS. Branches in this tree are unscaled.
FIG. 1.— A run of the branch class (3 classes) assignment genetic algorithm on the large lysozyme data set, using the MG94 x HKY85 model. The solid black line represents the score of the best fitting (in the AICc sense) model, and the gray plot shows the AICc score of every model in the order in which they were tried.
The AICc of the best fitting model assuming different numbers of branch classes are shown in table 1 for the small and large lysozyme data sets. Note that in both cases, moving from 3 to 4 branch classes did not result in a better fit as the models were collapsed to 3 rate classes during the search. The best fitting model identified using the GA can be compared with the single rate model and with the free ratio model by means of a likelihood ratio test. Because of the large number of possible models, stochasticity of observed data may lead to the best fitting model being chosen for spurious reasons. It is currently computationally infeasible to perform a full parametric or nonparametric bootstrap analysis, in which replicates of simulated data are subjected to analysis using the genetic algorithm. However, the robustness of the likelihood ratio test statistic to sampling error can be assessed using the approximate method of Shimodaira and Hasegawa (1999). This test can also be applied to comparison of the models considered by Yang (1998) and the best fitting models.
Table 1 Best AICc Found and Number of Models Examined for a Given Number of Branch Rate Classes under the MG94 x HKY85 Model
For the small data set (table 2), the best fitting model (we used the MG94 x HKY85 substitution model throughout) was significantly better than the single rate model (P = 0.0005), and not significantly different from the free ratio model (P = 0.93). The best fitting model had a lower AICc than the Y98 model—1824.39 versus 1830.72—and this difference was statistically significant as assessed by a Shimodaira-Hasegawa test (P < 0.007 using 10,000 replicates); 21 sites fitted better under the Yang (1998) model, and 109 sites fitted better under the GA model. We also note that these models are not nested.
Table 2 Model Fitting Results for the Small (7 sequences) Lysozyme Data Set
For the large data set (table 3), the best fitting model was also significantly better than the single ratio model (P = 3 x 10–8), and not significantly different from the free ratio model (P = 0.999991). The best fitting model had a significantly lower AICc than the Y98 model, with Shimodaira-Hasegawa P = 0.0156 (based on 10,000 replicates), and the models were not nested.
Table 3 Model Fitting Results for the Large (19 sequences) Lysozyme Data Set
For both the small and the large data sets (tables 2 and 3), the median dN/dS = ?sb for the Y98 model was much higher than than the median ?sb for the GA model. Both model misspecification and overfitting appear to contribute to this disparity. First, the median ?sb for the GA model is close to that of the free ratio model, whereas the median ?sb for the Y98 model is much higher than that of the free ratio model; this suggests that the Y98 model is failing to capture the pattern of variation in ?sb in the data. In addition, both the Y98 and the free ratio model have a class for which ?sb = ; clearly, this is biologically unrealistic and is the consequence of estimating ?sb from branches with no synonymous changes. In contrast, the GA model combines these branches into a single class with a high, but finite, value of ?sb.
There were a large number of models in the 95% confidence set for both the small (211) and large data sets (333). The evidence ratio for the best fitting model and the second best fitting model was 1.07 for the small data set. The two top models differ in class allocation for a single branch—the ancestor of colobines and cercopithecines (labeled as N3 in figure 2). The best fitting model allocated this branch to a class with ?sb = 0.61, while the second best model assigned it to the class with ?sb = 4.39. The third best model had the evidence ratio of 1.92. We can also compute the largest AICc needed for a model to be included in the 95% confidence set. For the smaller data set, this value was 1832.76, and hence the Y98 model is included in the credible set.
For the large data set, there were 19 models with essentially equal support (evidence ratio 1). This is likely due to the relatively few (4) sites per estimated model parameter and concomitant unreliable inference, as well as the the presence of short branches in the phylogeny. The threshold AICc value was found to be 2152.28, and thus the a priori model with AICc of 2167.54 is not included in the confidence set of models.
To assess the robustness of the conclusions drawn from the best fitting model, we calculated the support for an estimate of dN/dS = ?sb > 1 for each branch (fig. 2), averaging over models using their Akaike weights. The support for ?sb greater than 1 in the lineage leading to hominoids and colobines was 94.8% and 88.4%, respectively, for the small data set and 99.9% and 99.6%, respectively, for the large data set.
We next fitted site-to-site rate variation models using the branch assignments derived from running the genetic algorithm using a homogeneous rates model. Adding branch variation to a site-to-site model and site-to-site variation to a branch model yielded substantial AICc improvements for both data sets. Also, AICc scores for μMG94 x HKY85 and MG94 x HKY85 models indicate that variation in site-to-site mutational load fits lysozyme data better than variation in selection pressure over sites. Interestingly, the improvement in goodness-of-fit (in terms of the log likelihood) obtained by allowing both site-to-site variation in mutational load and branch variation was close to the improvement in fit when considering each effect separately. For the large data set, the difference in log-likelihood, logL, between the one ratio model and the GA 3-ratio model with variation in mutational load was 32.15, which is close to the sum of logL when considering branch variation and site-to-site variation separately (17.45 + 14.6 = 32.05). In addition, the estimates for ?sb were similar between models with and without site-to-site variation in mutational load. Thus, we can conclude that branch-specific variation in selection pressure and site-specific variation in mutational load play important (and apparently orthogonal) roles in shaping genetic variation in lysozyme. We also ran the GA using the μMG94 x HKY85 model on the small data set and found that the best fitting branch classes were the same as the ones found by the GA under the MG94 x HKY85 model, consistent with the apparent orthogonality between the models of site-to-site variation in mutational load and branch-specific variation in ?sb.
The model-averaging approach enables us to address issues which cannot be easily resolved using simple hypothesis testing. For example, we could compute probabilities of any two branches evolving under the same selective pressure (table 4), without making a priori assumptions about selective pressure acting on all other branches.
Table 4 Model Averaged Probabilities That Two Branches of the Small Lysozyme Phylogeny Have the Same ?sb = dN/dS Ratio
Discussion
Genetic algorithms have been employed in multiple studies to search for phylogenetic trees (Lewis 1998; Katoh, Kuma, and Miyata 2001; Brauer et al. 2002; Lemmon and Milinkovitch 2002; Kim, Lee, and Moon 2003; Shen and Heckendorn 2004). Our study is the first to our knowledge that uses a genetic algorithm to search for models of substitution in a phylogenetic context. Classical GA algorithms employ large population sizes (hundreds or thousands of individuals), and implicitly rely on the ability to compute the fitness of a single individual (in this case, the AICc of a given model) very quickly. Unfortunately, most molecular evolution problems, such as those considered here, are computationally challenging, and the number of fitness evaluations must be limited. Fortunately, there exist aggressive versions of GA, such as the CHC algorithm used in this study, which work well with smaller population sizes. Alternatively, a reversible jump Markov chain Monte Carlo procedure could be used to simultaneously estimate the number of classes of and the assignments to branches; however, this approach is likely to be far more computationally intensive than the GA approach pursued here. The GA method can be run on a small cluster of computers in a matter of hours and is therefore feasible for immediate use by researchers.
Using Akaike's Information Criterion as a measure of goodness-of-fit has been criticized for being too liberal, and may select models that are overly complex; instead, we have used a small sample AIC, AICc . In principle, information criteria other than AICc could be used, for example Schwarz's Bayesian Information Criterion (BIC) (Schwarz 1978), which is consistent (the probability that it will recover a true low-dimensional model approaches 1 as the sample size tends to infinity). However, the BIC assumes that the true model is contained within the candidate set of models, whereas the AIC does not assume that any of the candidate models is necessarily true; rather, it quantifies the discrepancy between the probability density generated by the model and the data as measured by Kullback-Liebler information (Kullback and Liebler 1951).
Given the large number of possible models and limited data, it is likely that a substantial number—tens or hundreds—of models will be consistent with the data. With the use of Akaike weights, which can be directly interpreted as conditional probabilities of the models, the robustness of conclusions to the choice of model can be assessed by averaging over models. Furthermore, an a priori biological model can be placed in the context of other credible models. If, for example, such a model falls outside the inferred 95% confidence set, it may not be well supported, as there are many more models which fit the data significantly better, even when the a priori model may be significantly better than the single-ratio model, as assessed by a nested comparison using a likelihood ratio test. Additionally we can test whether an a priori biological model is significantly worse than the best fitting model found by the GA using a Shimodaira-Hasegawa test. It is important to consider whether the a priori model is simply a submodel of the best fitting model, and whether the biological conclusions are altered under the best fitting model. Based on the small lysozyme data set, we find that the hypothesis that the lineages leading to Colobines and Hominoids have elevated nonsynonymous to synonymous substitution rates is consistent with our GA-based analysis. However, our GA results based on the large data set suggest that positive selection is ongoing following the radiation of these clades in contrast to the hypothesis of Messier and Stewart (1997).
When an a priori model compares favorably with the set of "good" models, we may be more confident in biological conclusions derived from studying such an a priori model. If, however, there are many alternative models which fit the data much better than our preconceived hypothesis, further exploration of the data, and perhaps a reassessment of our assumptions, is in order. In many practical cases the model itself is a nuisance component of the analysis, for example when one is interested in evolution along a specific lineage, and the rest of the phylogeny forms the "background." By weighing over credible models, quantities of biological interest (such as the probability that along a particular branch dN/dS > 1) may be derived without relying on a specific a priori model and concomitant assumptions, such as uniformity of selective pressure along background branches.
Acknowledgements
We thank Professor A.J. Leigh Brown for many helpful suggestions. This work was supported in part by the National Institutes of Health (grant numbers 5 R01 AI47745 and 5 U01 AI43638 Supp), the University of California Universitywide AIDS Research Program (grant number IS02-SD-701), and by a University of California, San Diego Center for AIDS Research/NIAID Developmental Award to S.D.W.F. (grant no. 2 P30 AI36214).
References
Abramowitz, M., I. Stegun, eds. 1972. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th printing. Dover, New York.
Akaike, H. 1974. A new look at the statistical model identification. IEEE Trans. Automatic Control 119:716–723.
Akaike, H. 1983. Information measures and model selection. Int. Stat. Inst. 44:139–149.
Brauer, M. J., M. T. Holder, L. A. Dries, D. J. Zwickl, P. O. Lewis, and D. M. Hillis. 2002. Genetic algorithms and parallel processing in maximum-likelihood phylogeny inference. Mol. Biol. Evol. 19:1717–1726.
Burnham, K. P., and D. R. Anderson. 2003. Model Selection and Multimodel Inference, 2nd ed. Springer, New York.
Eshelman, L. 1991. Foga-1. Morgan Kaufmann, Los Atlos, Calif.
Goldman, N., and Z. H. Yang. 1994. Codon-based model of nucleotide substitution for protein-coding DNA-sequences. Mol. Biol. Evol. 11:725–736.
Katoh, K., K. Kuma, and T. Miyata. 2001. Genetic algorithm-based maximum-likelihood analysis for molecular phylogeny. J. Mol. Evol. 53:477–484.
Kim, Y. H., S. K. Lee, and B. R. Moon. 2003. Optimizing the order of taxon addition in phylogenetic tree construction using genetic algorithm. Genetic and Evolutionary Computation—Gecco 2003, Pt Ii, Proceedings 2724:2168–2178.
Kosakovsky Pond, S., S. D. W. Frost, and S. V. Muse. 2000. HyPhy: hypothesis testing using phylogenies. Bioinformatics Advance Access published on October 27, 2004, doi:10.1093/bioinformatics/bti079.
Kullback, S., and R. Liebler. 1951. On information and sufficiency. Ann. Math. Stat. 22:79–86.
Lemmon, A. R., and M. C. Milinkovitch. 2002. The metapopulation genetic algorithm: an efficient solution for the problem of phylogeny estimation. Proc. Natl. Acad. Sci. USA 99:10516–10521.
Lewis, P. O. 1998. A genetic algorithm for maximum-likelihood phylogeny inference using nucleotide sequence data. Mol. Biol. Evol. 15:277–283.
Messier, W., and C. B. Stewart. 1997. Episodic adaptive evolution of primate lysozymes. Nature 385:151–154.
Muse, S. V., and B. S. Gaut. 1994. A likelihood approach for comparing synonymous and nonsynonymous nucleotide substitution rates, with application to the chloroplast genome. Mol. Biol. Evol. 11:715–724.
Schwarz, G. 1978. Estimating the dimension of a model. Ann. Stat. 6:461–464.
Shen, J., and R. B. Heckendorn. 2004. Discrete branch length representations for genetic algorithms in phylogenetic search. Appl. Evol. Comput. 3005:94–103.
Shimodaira, H., and M. Hasegawa. 1999. Multiple comparisons of log-likelihoods with applications to phylogenetic inference. Mol. Biol. Evol. 16:1114–1116.
Sugiura, N. 1978. Further analysis of the data by Akaike's information criterion and the finite corrections. Commun. Stat. Theory Methods A7:13–26.
Whitley, D. 2001. An overview of evolutionary algorithms: practical issues and common pitfalls. J. Inf. Software Technol. 43:817–831.
Yang, Z. H. 1998. Likelihood ratio tests for detecting positive selection and application to primate lysozyme evolution. Mol. Biol. Evol. 15:568–573.
Yang, Z. H., and R. Nielsen. 2002. Codon-substitution models for detecting molecular adaptation at individual sites along specific lineages. Mol. Biol. Evol. 19:908–917.
Yang, Z. H., R. Nielsen, N. Goldman, and A. M. K. Pedersen. 2000. Codon-substitution models for heterogeneous selection pressure at amino acid sites. Genetics 155:431–449.(Sergei L. Kosakovsky Pond)