CodonTest: modeling amino acid substitution preferences in coding sequences

Abstract

Codon models of evolution have facilitated the interpretation of selective forces operating on genomes. These models, however, assume a single rate of non-synonymous substitution irrespective of the nature of amino acids being exchanged. Recent developments have shown that models which allow for amino acid pairs to have independent rates of substitution offer improved fit over single rate models. However, these approaches have been limited by the necessity for large alignments in their estimation. An alternative approach is to assume that substitution rates between amino acid pairs can be subdivided into rate classes, dependent on the information content of the alignment. However, given the combinatorially large number of such models, an efficient model search strategy is needed. Here we develop a Genetic Algorithm (GA) method for the estimation of such models. A GA is used to assign amino acid substitution pairs to a series of rate classes, where is estimated from the alignment. Other parameters of the phylogenetic Markov model, including substitution rates, character frequencies and branch lengths are estimated using standard maximum likelihood optimization procedures. We apply the GA to empirical alignments and show improved model fit over existing models of codon evolution. Our results suggest that current models are poor approximations of protein evolution and thus gene and organism specific multi-rate models that incorporate amino acid substitution biases are preferred. We further anticipate that the clustering of amino acid substitution rates into classes will be biologically informative, such that genes with similar functions exhibit similar clustering, and hence this clustering will be useful for the evolutionary fingerprinting of genes.

Figure 1. Simulation studies used to derive the appropriate penalty term for .

Each panel plots the difference in log likelihood () normalized by the logarithm of the sample size (number of characters), between best fitting GA models with and rates (), against the number of sites in the alignment. For simulations with a single rate class we plotted , top right. Figures for multiple rate simulations (2–5 rates) show as black dots (left column); and as blue dots (right column). Values to the right of row report simulated rates for each class. The left column is a reflection of power, whereas the right column – of the degree of over-fitting. For the case where a single rate was simulated, the degree of over-fitting is the rate of false positives. The desired behavior for is achieved when the model with rate classes is preferred to models with , and rate classes. For a modified BIC criterion with , the former happens if (more definitively with increasing sample size), and the latter if (regardless of sample size).

Figure 2. Evolutionary rate estimation as “curve fitting.”

An example from HIV-1 polymerase gene alignment for which the inferred 7 non-synonymous rate classes. The idealized biological rate distribution (unobservable) is depicted by the dashed line. The goodness of fit, the complexity of the models, and the range of maximum likelihood parameter estimates are listed in the table.

Figure 3. Neighbor-joining trees built from matrices of pairwise substitution spectrum distances (Eq. 2) computed between different models fitted to the HIV-1 group M pol alignment, and between models inferred from different alignments.

Figure 4. Evolutionary rate clusters in structured GA models () inferred from the HIV-1 group M pol alignment.

Each cluster is labeled with the maximum likelihood estimate of its rate inferred under . The residues (nodes) are annotated by their biochemical properties and Stanfel class, and the rates (edges) are labeled with model-averaged () rate estimates. The style of an edge is determined by its cluster affinity, where high cluster affinities indicate that a large proportion of models in the credible set were consistent with the structured model.

Figure 5. Evolutionary rate clusters in structured GA models () inferred from the vertebrate rhodopsin protein alignment.

Each cluster is labeled with the maximum likelihood estimate of its rate inferred under . The residues (nodes) are annotated by their biochemical properties and Stanfel class, and the rates (edges) are labeled with model-averaged () rate estimates. The style of an edge is determined by its cluster affinity, where high cluster affinities indicate that a large proportion of models in the credible set were consistent with the structured model.

Figure 6. Correlations of lower substitution rates and property preservation in the HIV-1 group M pol alignment.

Model-averaged rates were stratified by whether or not they involved a change in polarity, charge or Stanfel class, the medians of two rate distributions were compared using a one sided Wilcoxon rank-sum test. We further correlated the magnitude of substitution rates with one of five property-based distances between the corresponding residues (defined in [18]) using a one-sided (negative correlation) Pearson product-moment correlation test.

Figure 7. Correlations of lower substitution rates and property preservation in the vertebrate rhodopsin alignment.

Model-averaged rates were stratified by whether or not they involved a change in polarity, charge or Stanfel class, the medians of two rate distributions were compared using a one sided Wilcoxon rank-sum test. We further correlated the magnitude of substitution rates with one of five property-based distances between the corresponding residues (defined in [18]) using a one-sided (negative correlation) Pearson product-moment correlation test.