Estimating Synonymous and Nonsynonymous Substitution Rates Under Realistic Evolutionary Models

Ziheng Yang2,* and Rasmus Nielsen{dagger}

*Department of Biology, University College London, England; and
{dagger}Department of Organismic and Evolutionary Biology, Harvard University


    Abstract
 TOP
 Abstract
 Introduction
 Methods for Estimating...
 Comparison of Methods for...
 Comparison of Human and...
 Discussion
 Program Availability and...
 Acknowledgements
 References
 
Approximate methods for estimating the numbers of synonymous and nonsynonymous substitutions between two DNA sequences involve three steps: counting of synonymous and nonsynonymous sites in the two sequences, counting of synonymous and nonsynonymous differences between the two sequences, and correcting for multiple substitutions at the same site. We examine complexities involved in those steps and propose a new approximate method that takes into account two major features of DNA sequence evolution: transition/transversion rate bias and base/codon frequency bias. We compare the new method with maximum likelihood, as well as several other approximate methods, by examining infinitely long sequences, performing computer simulations, and analyzing a real data set. The results suggest that when there are transition/transversion rate biases and base/codon frequency biases, previously described approximate methods for estimating the nonsynonymous/synonymous rate ratio may involve serious biases, and the bias can be both positive and negative. The new method is, in general, superior to earlier approximate methods and may be useful for analyzing large data sets, although maximum likelihood appears to always be the method of choice.


    Introduction
 TOP
 Abstract
 Introduction
 Methods for Estimating...
 Comparison of Methods for...
 Comparison of Human and...
 Discussion
 Program Availability and...
 Acknowledgements
 References
 
Estimation of synonymous and nonsynonymous substitution rates is important in understanding the dynamics of molecular sequence evolution (Kimura 1983;Citation Gillespie 1991Citation ; Ohta 1995Citation ). As synonymous (silent) mutations are largely invisible to natural selection (but see Akashi 1995Citation ), while nonsynonymous (amino- acid-replacing) mutations may be under strong selective pressure, comparison of the rates of fixation of those two types of mutations provides a powerful tool for understanding the mechanisms of DNA sequence evolution. For example, variable nonsynonymous/synonymous rate ratios among lineages may indicate adaptive evolution (Messier and Stewart 1997Citation ) or relaxed selective constraints along certain lineages (Crandall and Hillis 1997Citation ). Likewise, models of variable nonsynonymous/synonymous rate ratios among sites may provide important insights into functional constraints at different amino acid sites and may be used to detect sites under positive selection (Nielsen and Yang 1998Citation ).

The simplest problem in this regard is estimation of the numbers of synonymous (dS) and nonsynonymous (dN) substitutions per site between two sequences. In the past two decades, a number of intuitive methods have been suggested for this estimation. They involve ad hoc treatments that cannot be justified rigorously, and they will be referred to here as approximate methods. In common, they involve three steps. First, the numbers of synonymous (S) and nonsynonymous (N) sites in the sequences are counted. Second, the numbers of synonymous and nonsynonymous differences between the two sequences are counted. Third, a correction for multiple substitutions at the same site is applied to calculate the numbers of synonymous (dS) and nonsynonymous (dN) substitutions per site between the two sequences. Here, we use the notation of Nei and Gojobori (1986Citation , referred to later as "NG"), which appears to be the most commonly used approximate method; definitions of symbols are given in Table 1 . The reader is referred to Ina (1995, 1996)Citation for a recent discussion of important concepts. While the above strategy appears simple, well-known features of DNA sequence evolution, such as unequal transition and transversion rates and unequal nucleotide or codon frequencies, make it a real challenge to count sites and differences correctly.


View this table:
[in this window]
[in a new window]
 
Table 1 Definitions of Major Symbols

 
Miyata and Yasunaga (1980)Citation and Perler et al. (1980)Citation were the pioneers in this endeavor and developed the basic concepts (see also the simulation study of Gojobori 1983Citation ). The nucleotide substitution (mutation) model underlying the method of Miyata and Yasunaga (1980)Citation and its simplified version (Nei and Gojobori 1986Citation ) is the JC69 model (Jukes and Cantor 1969Citation ). JC69 is also the underlying mutation model in the method of Li, Wu, and Luo (1985)Citation , although the method uses the two-parameter model of Kimura (1980)Citation to correct for multiple hits at the same site. Here, we follow Ina (1995)Citation and use the word "mutation" to refer to DNA-level processes before the operation of natural selection at the protein level, although synonymous mutations may also be under selective constraints (Akashi 1995Citation ). As transitions are more likely to be synonymous at third positions than are transversions, use of the JC69 model to count sites tends to underestimate the number of synonymous sites and overestimate the number of nonsynonymous sites. The method of Li, Wu, and Luo (1985)Citation has been improved by Li (1993)Citation , Pamilo and Bianchi (1993)Citation , and Comeron (1995)Citation to correct for this bias in counting sites. Ina (1995)Citation appears to be the first to fully account for the transition/transversion bias in all steps of the estimation. We note that the underlying mutation (substitution) models used in all of the above methods are no more realistic than that of Kimura (1980)Citation . Moriyama and Powell (1997)Citation made a useful attempt to correct for biased base frequencies. Their modification, however, is limited to multiple-hit correction, and base frequency bias is not considered in the important steps of counting sites and differences.

A maximum-likelihood (ML) method for estimating dS and dN between two sequences was developed by Goldman and Yang (1994)Citation based on an explicit model of codon substitution. The ML method does not involve ad hoc approximations. Furthermore, the ML method is flexible in that knowledge of the substitution process such as transition/transversion bias, codon usage biases, and even chemical differences between amino acids can easily be incorporated into the model.

Relying on insights gained through previous methods, particularly ML estimation, we propose in this paper an approximate method for estimating dS and dN that accounts for two major features of DNA sequence evolution: the transition/transversion bias and the base (codon) frequency bias. We examine the similarities and differences among the ML method, the approximate method of this paper, and two other approximate methods by a consistency analysis of infinitely long sequences and computer simulation of finite data. A real data set is also analyzed using different estimation methods.


    Methods for Estimating Synonymous and Nonsynonymous RatesMaximum-Likelihood Estimation
 TOP
 Abstract
 Introduction
 Methods for Estimating...
 Comparison of Methods for...
 Comparison of Human and...
 Discussion
 Program Availability and...
 Acknowledgements
 References
 
First, we describe the ML method of Goldman and Yang (1994)Citation to introduce the notation and to provide justification for certain steps in our approximate method. An explicit model of codon substitution is required by the ML method. The model we consider in this paper is a simplified version of the model of Goldman and Yang (1994)Citation . The substitution rate from any sense codons i to j (i != j) is given by


Parameter {kappa} is the (mutational) transition/transversion rate ratio, and {omega} = dN/dS is the nonsynonymous/synonymous rate ratio. The equilibrium codon frequencies ({pi}j) are calculated using the nucleotide frequencies at the three codon positions; that is, codon frequencies are proportional to the products of nucleotide frequencies at the three codon positions. This approach was found to produce results similar to those obtained using all codon frequencies as free parameters, although codon frequencies are often quite different from those expected from nucleotide frequencies at codon positions (e.g., Goldman and Yang 1994Citation ; Pedersen, Wiuf, and Christiansen 1998Citation ; Yang and Nielsen 1998Citation ). We note that different assumptions about codon frequencies can be easily incorporated into the ML method and the approximate method of this paper. Equation (1) is similar to the simulation model of Gojobori (1983)Citation and the likelihood model of Muse and Gaut (1994)Citation , although those authors did not consider the transition/transversion rate bias or different base frequencies at the three codon positions.

The diagonal elements of the rate matrix, Q = {qij}, are determined by the mathematical requirement (e.g., Grimmett and Stirzaker 1992Citation , p. 241) that the row sums are zero:


Because time and rate are confounded, we multiply the rate matrix by a scaling factor so that the expected number of nucleotide substitutions per codon is one:


This scaling means that time t (or, equivalently, branch length or sequence distance) is measured by the expected number of (nucleotide) substitutions per codon.

The transition probability matrix is calculated as

where pij(t) is the probability that codon i will become codon j after time t. The probability of observing a codon site with codons i and j in the two sequences is then

The log-likelihood function is given by the multinomial probability


where nij is the number of sites occupied by codons i and j in the two sequences. The codon frequencies {pi}i are estimated using the observed nucleotide frequencies at the three codon positions in the two sequences. Parameters t, {kappa}, and {omega} are estimated by maximizing the likelihood function numerically, and the estimates are used to calculate dS and dN. Specifically, the proportions of synonymous and nonsynonymous substitutions are given as


and {rho}*N = 1 - {rho}*S, respectively. The summation is taken over all codon pairs i and j (i != j) that code for the same amino acid, and aai is the amino acid encoded by codon i. The numbers of synonymous and nonsynonymous substitutions per codon are then t{rho}*S and t{rho}*N, respectively. The proportions of synonymous and nonsynonymous sites are defined as the proportions of synonymous and nonsynonymous "mutations" before the operation of natural selection at the amino acid level (Goldman and Yang 1994Citation ; Ina 1995Citation ). These can be calculated in a manner similar to equation (7) as {rho}1S and {rho}1N (equivalent to {rho}{infty}S and {rho}{infty}N in Goldman and Yang 1994Citation ), using the ML estimate of {kappa} but with {omega} = 1 fixed. The numbers of synonymous and nonsynonymous sites per codon are 3{rho}1S and 3{rho}1N, respectively. The numbers of synonymous and nonsynonymous substitutions per site are then dS = t{rho}*S/(3{rho}1S) and dN = t{rho}*N/(3{rho}1N), respectively (see Table 1 ).

A New Approximate Method
We suggest an approximate method for estimating dS and dN using the strategy adopted by previous authors: counting sites, counting differences, and correcting for multiple hits. In all three steps, we take into account the transition/transversion rate bias and the base (codon) frequency bias. Approximately, our method is based on the HKY85 nucleotide mutation (substitution) model (Hasegawa, Kishino, and Yano 1985Citation ). Although this is not the most general model available, it accounts for two most important features of the mutation process, that is, the transition/transversion bias and unequal base frequencies. Previous results (e.g., Yang 1994aCitation ) suggest that adding further complication is often unnecessary.

Estimating the Transition/Transversion Rate Ratio ({kappa})
We use the fourfold-degenerate sites at the third codon positions and the nondegenerate sites to estimate {kappa}. Mutations at the fourfold-degenerate sites do not change the amino acid, and thus the transition/transversion rate bias at those sites should reflect the mutational bias. Mutations at nondegenerate sites all lead to amino acid changes and can also be used to estimate {kappa} (see eq. 1). Here, we assume that different nonsynonymous substitutions have the same rate irrespective of the pair of amino acids involved, although the assumption is unrealistic (Yang, Nielsen, and Hasegawa 1998Citation ). We calculate an average of {kappa}, weighted by the numbers of nucleotide sites in the two site classes. Since no simple formula is available for estimating {kappa} under the HKY85 model, we use the formula for the F84 model (Yang 1994 bCitation ) instead, relying on the similarity of the two models. We calculate


where T and V are proportions of transitional and transversional differences, respectively, and {pi}Y = {pi}T + {pi}C and {pi}R = {pi}A + {pi}G. Then,

and


(Yang 1994bCitation ). The estimated {kappa}F84 is then transformed to {kappa} (i.e., {kappa}HKY85) using the following formula (see Goldman 1993Citation )


For data of multiple sequences, we suggest estimating a common {kappa} by averaging estimates from all pairwise comparisons and using the combined estimate of {kappa} in the calculation of pairwise dN and dS rates.

Counting Synonymous and Nonsynonymous Sites
Ina’s (1995)Citation Table 1 for counting synonymous and nonsynonymous sites in each codon is correct for mutation models more general than that of Kimura (1980)Citation , although Ina’s Table involves minor errors for codons that can change to stop codons in one step. In general, synonymous and nonsynonymous sites can be counted as in the ML method discussed above for any codon-substitution model (Goldman and Yang 1994Citation ). We count sites using codons in the two compared sequences, rather than the equilibrium codon frequencies expected from the model (see discussion in Yang and Nielsen 1998Citation ). As there should be about 4% loss of sites due to mutations to stop codons, this scaling means that we are slightly underestimating dS and dN, although the {omega} ratio is not affected (Yang and Nielsen 1998Citation ). The numbers of sites (S and N) are scaled so that S + N = 3Lc, where Lc is the number of codons. Nucleotide frequencies at synonymous and nonsynonymous sites are recorded and used later for multiple-hit corrections.

Counting Synonymous and Nonsynonymous Differences
Observed nucleotide differences between the two sequences are classified into four categories: synonymous transitions, synonymous transversions, nonsynonymous transitions, and nonsynonymous transversions. When the two compared codons differ at one position, the classification is obvious. When they differ at two or three positions, there will be two or six parsimonious pathways along which one codon could change into the other, and all of them should be considered. Since different pathways may involve different numbers of synonymous and nonsynonymous changes, they should be weighted differently. Miyata and Yasunaga (1980)Citation and Li, Wu, and Luo (1985)Citation made valuable attempts to weight pathways. They also pointed out that equal weighting of pathways, which is used in later methods such as those of Nei and Gojobori (1986)Citation and Ina (1995)Citation , biases estimates of {omega} toward 1; that is, equal weighting tends to overestimate {omega} when {omega} < 1 and to underestimate {omega} when {omega} > 1.

The most appropriate weights are the relative probabilities of pathways, which we will use in the new approximate method. The probabilities depend on the parameters being estimated: the sequence divergence level (t), the transition/transversion rate ratio ({kappa}), and the dN/dS ratio ({omega}). For given values of t, {kappa}, and {omega}, it is easy to construct the rate matrix Q and calculate the transition probability matrix P(t) (see eqs. 1 and 4). We use the Taylor expansion in this case for its fast speed.


The number of terms used is determined by a preset accuracy level. The weight, that is, the probability for each pathway, is calculated as the product of the probabilities of all changes involved in the pathway. Pathways involving stop codons are given weight 0. If all pathways involve stop codons (for example, between AAG and TGG in the mammalian mitochondrial code), ad hoc decisions have to be made.

An example is given in Figure 1 using a pair of codons in the mitochondrial genes of the human and the orangutan. The concatenated sequence of 12 genes on the H-strand has 3,331 codons (see below). Between the two sequences, 1,198 codons are different at one position, 151 are different at two positions, and 21 are different at all three positions. The 329th codon is TTA in the human and CTC in the orangutan. Transition probabilities for changes involved in each of the two pathways are given in Figure 1 , calculated using the estimates obtained by the new method (t = 0.873, {kappa} = 10.61, and {omega} = 0.057). The probability for the first path is then pTTA,TTC(t) x pTTC,CTC(t) = 0.00541 x 0.04974 = 0.00027, while that for the second path is 0.29219 x 0.08679 = 0.02536. Weights for the two pathways are thus 0.011 and 0.989, and there are 0.022 nonsynonymous and 1.978 synonymous differences between the two codons. Since the nonsynonymous rate is much lower than the synonymous rate ({omega} < 1), the first path is much less likely. Equal weighting of pathways would give one synonymous and one nonsynonymous difference between the two codons. Other codon pairs may not be so extreme, as the different pathways may involve the same numbers of different types of changes.



View larger version (7K):
[in this window]
[in a new window]
 
Fig. 1.—Two parsimonious pathways between codons TTA and CTC. Probabilities are calculated using parameter estimates for the human and orangutan mitochondrial genes.

 
Correcting for Multiple Hits Using Estimated Numbers of Sites and Differences
We use the distance formula (eq. 10) for the F84 model (Tateno, Takezaki, and Nei 1994Citation ; Yang 1994 bCitation ) to correct for multiple substitutions at the same site to calculate dS and dN. For each of the synonymous and nonsynonymous site classes, the proportions of transitional and transversional differences are calculated separately to give T and V for use in equation (8). Different base frequencies are also used for the two site classes. The correction formula used here, as well as those used in all previous approximate methods, is ad hoc. The formula is derived from a Markov process of nucleotide substitution with four states, where each nucleotide can change into one of three other nucleotides. When the formula is used for synonymous (or nonsynonymous) sites only, this basic assumption of the Markov model is violated (Lewontin 1989Citation ). However, a proper treatment of the evolutionary process at synonymous and nonsynonymous sites appears to require the use of a codon substitution model, as in the likelihood method, and an analytical derivation of a correction formula based on such a model seems intracTable.

The Algorithm
Our method for estimating dS and dN can be summarized in the following iterative algorithm.

1. Estimate {kappa} from the fourfold-degenerate sites and the nondegenerate sites under the HKY85-F84 model using base (codon) frequencies from the real data. The estimated {kappa} is used in later steps.

2. Count the numbers of synonymous and nonsynonymous sites (S and N, respectively) using the estimated {kappa} and the observed base (codon) frequencies.

3. Choose starting values for t and {omega} (e.g., using estimates from the NG method).

4. Count the numbers of synonymous and nonsynonymous differences (both transitions and transversions) using {kappa}, the codon frequencies, and the current values of t and {omega}. The transition probability matrix P(t) is calculated by equation (12) and used to weight pathways when the two codons differ at more than one position. This step generates the proportions of transitional (T) and transversional (V) differences for each of the synonymous and nonsynonymous site classes.

5. Correct for multiple hits to calculate dS and dN using counts of sites and differences and base frequencies at synonymous and nonsynonymous sites. This step updates t and {omega}: t = dS x 3S/(S + N) + dN x 3N/(S + N), and {omega} = dN/dS.

6. Repeat steps 4–5 until the algorithm converges.

In general, two or three rounds of iteration are sufficient. Some variations of the above algorithm are possible. For example, one may use an estimate of {kappa} obtained externally. Furthermore, no iteration is needed if pathways are weighted equally when counting differences.


    Comparison of Methods for Estimating dS and dN
 TOP
 Abstract
 Introduction
 Methods for Estimating...
 Comparison of Methods for...
 Comparison of Human and...
 Discussion
 Program Availability and...
 Acknowledgements
 References
 
We examine the performance of the following four methods for estimating dS and dN: ML (Goldman and Yang 1994Citation ), NG (Nei and Gojobori 1986Citation ), the method of Ina (1995Citation ; Method I), and the method of this paper. PAML (Yang 1999Citation ) was used for the ML and NG methods, and Ina’s program was used for Ina’s method. For error-checking, independent programs for the new method were written by the two authors.

Two approaches are taken to evaluate the methods. The first examines infinitely long sequences and may be termed a "consistency analysis." Instead of the observed codons in the two sequences, the data consist of the expected frequencies (fij) of all 61 x 61 codon "site patterns," calculated using equation (5) for given parameters t, {kappa}, {omega}, and codon frequencies {pi}j. ML estimates are known to be statistically consistent when the model is correct (Stuart, Ord, and Arnold 1999Citation , chapter 18). Since the dN and dS rates (and their ratio) are defined as functions of parameters t, {kappa}, {omega}, and {pi}j, ML estimates of dN and dS will also be consistent. Approximate methods, including the method of this paper, involve ad hoc approximations and in general do not give the true values as estimates. They are statistically inconsistent. However, a good approximate method should not deviate too far from the truth with an infinite amount of data. The second approach we take is computer simulation. Finite data sets are generated by simulation and then analyzed by different methods to examine their biases and sampling variances.

We examine effects of the transition/transversion rate ratio ({kappa}), base (codon) frequencies, and the selective pressure on the gene reflected in parameter {omega}. We initially fix t = 1 nucleotide substitution per codon, although the effect of sequence divergence is examined later. For a neutral gene ({omega} = 1), this translates to 1/3 synonymous and 1/3 nonsynonymous substitutions per site. Three values of {omega} are considered: {omega} = 1 (no selection), {omega} = 0.3 (purifying selection) and {omega} = 3 (positive selection). Estimates of {omega} (dN/dS) from real data vary widely from gene to gene, and {omega} = 0.3 appears to represent moderate purifying selection (see, e.g., Ohta 1995Citation ; Li 1997Citation ; Yang and Nielsen 1998Citation ; Eyre-Walker and Keightley 1999Citation ). There are not many genes under positive selection, but estimates at about {omega} = 3 are found in real data (e.g., Lee, Ota, and Vacquier 1995Citation ; Messier and Stewart 1997Citation ). Three sets of base frequencies at codon positions are used. The first set has equal base (codon) frequencies. The second set is from primate mitochondrial protein-coding genes and has very biased base frequencies. The third set is from HIV env genes (Table 2 ). The universal genetic code is used in both the consistency analysis and the computer simulation.


View this table:
[in this window]
[in a new window]
 
Table 2 Base Frequencies at Codon Positions in Two Data Sets

 
Consistency Analysis Using Infinite Data
Consistency is the property that the estimate converges to the true value of the parameter as the amount of data approaches infinity. While consistency is a weak requirement, the approximate methods examined here are all inconsistent. It is nevertheless interesting to examine which steps of the approximate methods (i.e., counting sites, counting differences, and correcting for multiple hits) cause the bias. This is relatively easy since infinite data do not involve sampling errors, and estimates of sites (S and N) and rates (dS and dN) can be directly compared with the correct values.

Estimates of {omega} by different methods are plotted against the transition/transversion rate ratio {kappa} for different values of the dN/dS ratio ({omega}). Results for the three sets of codon frequencies are shown in Figure 2AI.



View larger version (33K):
[in this window]
[in a new window]
 
Fig. 2.—Estimates of {omega} (=dN/dS) as a function of {kappa} for infinite data. The methods are ML (dotted lines, true values), NG (dashed lines), Ina’s (1995) method (mixed lines), and the method of this paper (thick lines). Sequences of 2,000,000 codons were simulated for Ina’s method, while for other methods, the data are the expected frequencies of the 61 x 61 site patterns. The universal genetic code is used. Three selection schemes are used: (A, D, and G) no selection ({omega} = 1), (B, E, and H) purifying selection ({omega} = 0.3), and (C, F, and I) positive selection ({omega} = 3). Three sets of base/codon frequencies are used: (AC) equal frequencies, (DF) frequencies from the primate mitochondrial protein-coding genes, and (GI) frequencies from the HIV-1 env genes.

 
Equal Codon Frequencies
We consider the NG method first. When base (codon) frequencies are equal and transition and transversion rates are equal ({kappa} = 1), assumptions of the NG method are largely satisfied. In this case, NG indeed gives estimates close to the true values. Estimates of {omega} given by NG when {kappa} = 1 are 1.001, 0.318, and 2.523 for the true {omega} = 1, 0.3, and 3, respectively (Fig. 2AC ). The method is biased toward 1 when the true {omega} != 1 due to its use of equal weighting of pathways when counting sites. When there is transition bias ({kappa} > 1), NG underestimates the {omega} ratio, and the bias is more serious when the transition bias is more extreme. This bias is mainly generated in the step of counting sites.

The case of {kappa} = 10 is explored in Table 3 , which shows that NG substantially underestimates {omega} and gives 0.669, 0.216, and 1.812 for {omega} = 1, 0.3, and 3, respectively. In this case, the proportion of synonymous sites (S%) should be 33%, but NG (assuming {kappa} = 1) gives 26% (Yang and Nielsen 1998Citation , Fig. 3 ). Use of equal weighting (assuming {omega} = 1) in counting differences by NG tends to bias the estimate of {omega} toward 1. However, compared with the bias in counting sites, the bias in counting differences is much less important because there may not be many pairs of codons that are different at two or three positions and because different pathways may involve the same numbers of synonymous and nonsynonymous changes. As mentioned above, NG underestimates S considerably (25.5/32.5 = 0.785) when {kappa} = 10, and almost all of this underestimation is translated into the overestimation of dS (0.3333/0.4134 = 0.806 for {omega} = 1). For similar reasons, the underestimation of {omega} by the NG method is more serious for {omega} = 3 than for {omega} = 0.3 (Fig. 2B and C ). In the latter case, equal weighting assuming {omega} = 1 in the NG method counterbalances the effect of ignoring {kappa} in counting sites, while in the former case, the two biases are in the same direction (Table 3 ).


View this table:
[in this window]
[in a new window]
 
Table 3 Estimates of dN, dS, and {omega} in Infinite Data by Different Methods when {kappa} = 10

 


View larger version (14K):
[in this window]
[in a new window]
 
Fig. 3.—Square root of the mean square error, MSE(), of the estimated {omega} ratio, as a function of the sequence length Lc, calculated by computer simulation. Parameter values are (A) t = 1, {kappa} = 20, {omega} = 0.3, with base/codon frequencies from the primate mitochondrial genes (Table 2 ), and (B) t = 1, {kappa} = 2, {omega} = 3, with base/codon frequencies from the HIV-1 env genes (Table 2 ). The universal genetic code is used.

 
Ina’s (1995)Citation method gives quite reliable estimates of {omega} for large values of {kappa} (Fig. 2AC ). For small values of {kappa}, the method tends to overestimate {omega} when {omega} < 1 and to underestimate {omega} when {omega} > 1. This pattern appears to be due to the use of equal weighting of pathways in counting differences. The method of this paper underestimates {omega} slightly when the transition/transversion bias is weak (that is, when {kappa} is close to 1) and when {omega} != 1 (Fig. 2BC ).

Primate Mitochondrial Codon Frequencies
Results obtained using base frequencies at the three codon positions from the primate mitochondrial genes (see Table 2 ) are shown in Figure 2DF. Estimates of S, dS, and dN for {kappa} = 10 are shown in Table 3 . The results are very different from those of Figure 2AC under equal codon frequencies. Except for small values of {kappa} ({kappa} < 2), Ina’s method performs more poorly than NG. The two methods give very different counts of sites (S). While transition bias always leads to more synonymous sites, the effect of base frequency bias is more complicated. Extreme codon usage bias can cause the proportion of synonymous sites to range from 0% (e.g., when only codons TTC and TTA are present in the sequences) to 100% (e.g., when only codons CTT, CTC, CTA, and CTG are present). There tend to be more synonymous sites if the two most frequent nucleotides at third positions are both purines or both pyrimidines.

For the mitochondrial genes, S is much smaller than expected under equal base frequencies, causing NG to overestimate rather than underestimate S. For example, NG gives S% = 26.7%, which is higher than the correct value at either {kappa} = 1 (23.1%) or {kappa} = 10 (23.7%) (Table 3 ). The overestimation of S caused by ignoring the base frequency bias more than compensates for the underestimation caused by ignoring the transition bias. As a result, NG overestimates {omega} when {omega} = 1 (with estimates from 1.1 to 1.3; Fig. 2D ) or {omega} = 0.3 (with estimates from 0.42 to 0.46; Fig. 2E ). When {omega} = 3, equal weighting of pathways (assuming {omega} = 1) in counting differences combined with the assumption of no transition bias ({kappa} = 1) in counting sites cancels the effect of the base frequency bias in counting sites, such that NG produces a quite reliable estimate of {omega} (Fig. 2F ). Ina’s (1995)Citation method, by considering the transition bias alone and ignoring the base/codon frequency bias, substantially overestimates the proportion of synonymous sites and overestimates {omega} (Table 3 and Fig. 2DF ). Nevertheless, it should be noted that the observed pattern depends on the particular set of codon frequencies. For frequencies from other data sets, NG may be considerably worse than Ina’s method.

The method of this paper is slightly better than NG for {omega} = 1 although the two methods have opposite biases. When {omega} = 0.3, the new method has little bias. When {omega} = 3, the new method underestimates {omega}, with estimates from 2.5 to 2.6. Since the new method counts sites correctly (see Table 3 ), the bias must be due to counting of differences and correction for multiple hits. Table 3 suggests that the new method underestimates both dS and dN, but the underestimation of dN is more serious, leading to underestimation of the {omega} ratio. Apart from the case in which {omega} = 3, the new method is better than both NG and Ina’s method.

HIV env Codon Frequencies
Figure 2GI shows estimates of {omega} when base/codon frequencies from the HIV envelope genes (see Table 2 ) are used. Base frequencies in this gene are less biased than are those in the mitochondrial genes, and the effect of ignoring the base frequency bias is minor. For example, the correct proportion of synonymous sites at {kappa} = 1 is 21.9%, while NG gives 24.0%, with very slight overestimation. Patterns in Figure 2GI are quite similar to those for equal codon frequencies (Fig. 2AC ). Exactly at {kappa} = 1, NG gives the estimates 1.105, 0.371, and 2.554 when the true {omega} = 1, 0.3, and 3, respectively. The estimates are biased toward 1, mainly due to the use of equal weighting in counting differences. When {kappa} = 10, NG underestimates the proportion of synonymous sites (24.0% vs. the correct value, 28.5%). The bias is not as extreme as that in the case of equal codon frequencies, as unequal base/codon frequencies appear to counterbalance the effect of transition bias to some extent (Tables 3 and 4 ).

Ina’s (1995)Citation method overestimates the {omega} ratio because it ignores the base frequency bias and thus overestimates the number of synonymous sites. The bias is not as extreme as it is for mitochondrial genes. The new method gives estimates very close to the true values for {omega} = 1 and {omega} = 0.3. When {omega} = 3, the new method slightly underestimates the ratio, as in the case of mitochondrial codon frequencies.

Computer Simulations
The data consist of a pair of codon sequences and are simulated by sampling codon site patterns from the multinomial distribution specified by the site pattern probabilities fij (eq. 5). The sequence has Lc = 100 or 500 codons. Three values of {kappa} are used: 1 (no bias), 2 (small bias), and 20 (large bias). Most estimates of {kappa} from nuclear genes are in the range (1.5, 5), so a value of 2 is typical. Estimates from mitochondrial genes vary considerably among data sets, from 2 or 3 to over 100.

The averages of the {omega} estimates among simulated replicates are listed in Table 4 for the three sets of codon frequencies. Standard errors for the ML estimates are also presented, while those for other methods (not shown) are very small due to the use of many more replicates. Averages of dS and dN are calculated for all methods but not shown. We note that the simulation results are highly consistent with those found for infinite data, discussed above. For example, if a method gives estimates smaller than the true value in infinite data, it tends to have negative biases in finite samples as well.


View this table:
[in this window]
[in a new window]
 
Table 4 Average Estimates of {omega} in Simulated Replicates

 
ML estimates are known to be often biased in small samples. Table 4 shows that MLEs of {omega} are nearly unbiased when {omega} = 1 or 0.3 for all three sets of codon frequencies and for all values of {kappa}. However, it is biased to larger values when {omega} = 3. Although the bias is small in large genes (with 500 codons), it can be quite large for small genes (Lc = 100), especially when {kappa} is small.

The NG method has little bias if codon frequencies are equal and if there is no transition/transversion bias ({kappa} = 1). When {omega} != 1, the method tends to bias toward 1 due to its use of equal weighting in counting sites. The bias is nevertheless small. These results agree well with previous simulations by Ota and Nei (1994)Citation and Muse (1996)Citation , who used similar simple models to examine the performance of NG. However, NG is biased in most other parameter combinations. The biases in general agree with findings of the consistency analysis (Fig. 2 ). In particular, ignoring the transition/transversion bias leads to underestimates of {omega} and ignoring unequal base frequencies leads to overestimates of {omega}. For equal codon frequencies and no selection ({omega} = 1), NG gives severe underestimates of {omega} when {kappa} is large. The effects of the transition/transversion rate bias and base frequency biases tend to cancel each other, such that NG has smaller biases than ML when {omega} = 3 for the mitochondrial codon frequencies. In almost all other cases, ML has smaller biases than NG.

Ina’s (1995)Citation method has small biases when base frequencies are equal. The method tends to overestimate {omega} when {omega} < 1 and to underestimate {omega} when {omega} > 1, probably due to its use of equal weighting in counting sites. This is the same pattern as that found in infinite data (Fig. 2AC ). Ina’s method considerably overestimates {omega} for all values of {omega} under the mitochondrial codon frequencies, probably because it overestimates the number of synonymous sites. For the HIV env codon frequencies, the method overestimates {omega} when {omega} <= 1 and underestimates {omega} when {omega} > 1, as noted for infinite data (Fig. 2GI ).

The new method appears to have little bias over most of the parameter space examined (Table 4 ). When {omega} <= 1, it is less biased than NG or Ina’s (1995)Citation method. When {omega} = 3, it tends to overestimate {omega} in small samples, like the ML method, but the bias seems smaller than that of ML. The new method appears to provide a close approximation of ML over the range of parameter values examined.

Since all methods are biased for at least some parameter combinations, the mean squared error (MSE) may be an appropriate criterion by which to compare methods. The MSE of a parameter estimator is defined as MSE() = E( - {theta})2, where {theta} is the true value. Since MSE() = Var() + [E() - {theta}]2, this measures both bias and variance. The square root of the MSE is plotted in Figure 3 against the sequence length (number of codons). Two parameter combinations are considered. The first is for mitochondrial genes with a strong transition bias ({kappa} = 20) and purifying selection ({omega} = 0.3), and the second is for the HIV env gene with moderate transition bias ({kappa} = 2) and positive selection ({omega} = 3). In the first case (Fig. 3A ), ML and the method of this paper performed very similarly, and both have the smallest MSEs, while Ina’s (1995)Citation method has very large MSEs due to the positive bias of the method (Fig. 2E ). In the second case (Fig. 3B ), ML performed much worse than other methods for short genes (Lc < 300 codons) due to its large positive bias at {omega} = 3, while for large genes (Lc > 500), ML and the new method are better than NG and Ina’s method. In both cases, the new method lies between NG and ML.

We also performed a small-scale simulation to examine the effect of sequence divergence level (t). The results are shown in Figure 4 . We examine two sets of parameter values, with the sequence length fixed at Lc = 500. In the first case (Fig. 4A ), equal codon frequencies are used with {kappa} = 2 and {omega} = 0.3. The new method of this paper overestimates {omega} at small divergences but underestimates {omega} at large divergences. Other methods are insensitive to sequence divergence level. ML and the new method are less biased than NG and Ina’s method. In the second case (Fig. 4B ), mitochondrial codon frequencies are used with {kappa} = 20 and {omega} = 0.3. In this case, ML and the new method have little bias over the whole range of the sequence divergence level. Note that the synonymous rate is quite high, with dS = 0.71 at t = 1, and dS = 1.1 at t = 1.5. NG and Ina’s method involve positive biases, and the biases become more serious when the sequences are more divergent. Although average estimates of both dN and dS by NG increase with t, dN increases at a faster rate, such that the {omega} ratio increases with the increase of t. Muse (1996)Citation discussed the fact that at high sequence divergences, NG does not produce distance estimates linear with time.



View larger version (19K):
[in this window]
[in a new window]
 
Fig. 4.—The averages of estimates of {omega} as a function of sequence divergence level t. The sequence length is Lc = 500 codons. A, Equal codon frequencies are assumed, with {kappa} = 2 and {omega} = 0.3. B, Codon frequencies from primate mitochondrial genes are used, with {kappa} = 20 and {omega} = 0.3. Each average was obtained by simulating 2,000 replicates for NG and YN (the method of this paper), 500 for the method of Ina (1995), and 100 for ML.

 

    Comparison of Human and Orangutan Mitochondrial Genes
 TOP
 Abstract
 Introduction
 Methods for Estimating...
 Comparison of Methods for...
 Comparison of Human and...
 Discussion
 Program Availability and...
 Acknowledgements
 References
 
The concatenated sequences of the 12 protein-coding genes on the H-strand of the mitochondrial genome from the human (Homo sapiens, GenBank accession number D38112) and the orangutan (Pongo pygmaeus p., GenBank accession number D38115) are compared using different methods. The results are shown in Table 5 . We also included the method of Li (1993)Citation in the comparison, implemented in X. Xia’s DAMBE program (available athttp://web.hku.hk/~xxia). ML is applied with different assumptions concerning the transition bias and the codon frequency bias. Note that estimates of the dN/dS ratio vary by up to threefold depending on the method/model used. The pattern is especially revealing for ML estimates under different models. For example, accounting for the transition/transversion rate bias (Fequal, {kappa} estimated) increased the proportion of synonymous sites from 25% to 33%, and increased the dN/dS ratio from 0.093 to 0.130. Accounting further for the codon usage bias (F60, {kappa} estimated) decreased the proportion of synonymous sites to 28%, with an estimate of {omega} = 0.045. The pattern is the same as that in the consistency analysis and the computer simulation discussed above. The results demonstrate that the estimation method is important for estimating dN and dS (and their ratio {omega}) from real data and that methods accounting for both the transition bias and base/codon frequency bias should be used. This conclusion is consistent with Ina (1995)Citation , who found that none of the approximate methods he examined performed well when base/codon frequencies were extreme.


View this table:
[in this window]
[in a new window]
 
Table 5 Proportions of Synonymous Sites (S%) and dS and dN Rates Between the Human and Orangutan Mitochondrial Genes

 
Furthermore, we note that estimates from the NG method are similar to those of ML assuming no transition bias ({kappa} = 1 fixed) and no base frequency bias, and estimates obtained from Li (1993)Citation are similar to ML accounting for the transition bias ({kappa} estimated) and no base frequency bias. Ina’s (1995)Citation program did not run for this data set. Nevertheless, we expect Ina’s method to be close to that of Li (1993)Citation or ML accounting for the transition bias and no base frequency bias. Some minor differences in implementation between ML and the corresponding approximate methods were discussed by Yang and Nielsen (1998)Citation . Those results suggest that ad hoc treatments involved in the approximate methods may not have introduced too much bias and that failure to account for the transition bias and base frequency bias appears to be more important. However, Muse (1996)Citation points out that at high sequence divergence levels, ad hoc treatments (such as those used in multiple-hit correction) in approximate methods may become a more serious problem (see also Fig. 4 ).


    Discussion
 TOP
 Abstract
 Introduction
 Methods for Estimating...
 Comparison of Methods for...
 Comparison of Human and...
 Discussion
 Program Availability and...
 Acknowledgements
 References
 
The approximate method of this paper accounts for two major features of DNA sequence evolution: transition bias and base/codon frequency bias. The consistency analysis of infinite data and the computer simulation of finite data suggest that the new method has smaller biases than either NG or Ina’s (1995)Citation method for almost all parameter combinations examined and produces estimates similar to those obtained with ML. The method may thus be useful for large-scale screening, when ML may be too time-consuming. Our analyses also suggest that NG and Ina’s method may involve large biases when there exist transition/transversion rate bias and base/codon frequency bias, and that it is important to account for those features of DNA sequence evolution.

The ML method for pairwise comparison is less biased and has a lower MSE than the approximate methods for almost all parameter combinations. Only for short sequences and high {omega} ratios does it involve a positive bias and perform more poorly than some of the approximate methods. We suggest that, in general, the ML method, which accounts for both the transition bias and the codon usage bias, should be the preferred method for estimating dS and dN between two sequences. Only in the case of very short sequences may it be advantageous to use simpler models. In the course of this study, we realized that correcting for biases involved in the NG method is extremely complicated, despite the fact that the method is well known for its simplicity. In contrast, ML is conceptually much simpler, mainly because the probability theory employed by the method takes care of the difficult tasks of weighting evolutionary pathways and correcting for multiple hits, with no need for ad hoc approximations. Specifically, the Chapman-Kolmogorov theorem (e.g., Grimmett and Stirzaker 1992Citation , p. 239) states that pij(t) = {Sigma}k pik(s)pkj(t - s) for any 0 <= s <= t; that is, the probability that codon i changes to codon j over time t is a sum over all possible codons (k) at any intermediate time point s. This obvious result ensures that the likelihood calculation (eqs. 4–6) accounts for all possible pathways of changes between the two codons, weighting them appropriately according to their relative probabilities of occurrence.

The major advantage of ML appears to lie in its flexibility in simultaneous comparison of multiple sequences, taking into account their phylogenetic relationship. Hypotheses concerning variable dN/dS ratios among lineages (Yang 1998Citation ; Yang and Nielsen 1998Citation ) or among sites (Nielsen and Yang 1998Citation ) can be tested using the likelihood ratio test. The ML model can easily be extended to include important features of DNA sequence evolution such as the dependence of nonsynonymous rates on the chemical properties of the amino acids (Yang, Nielsen, and Hasegawa 1998Citation ).


    Program Availability and Performance
 TOP
 Abstract
 Introduction
 Methods for Estimating...
 Comparison of Methods for...
 Comparison of Human and...
 Discussion
 Program Availability and...
 Acknowledgements
 References
 
A C program implementing the approximate method of this paper will be included in the PAML package, available at http://abacus.gene.ucl.ac.uk/software/paml.html. On a fast Pentium II, each pairwise comparison takes about 10–20 s by ML and a few seconds by the method of this paper. If pathways are weighted equally in counting differences in the new method, iteration will not be needed, and the method will be about as fast as other approximate methods such as NG, which seem to finish instantaneously.


    Acknowledgements
 TOP
 Abstract
 Introduction
 Methods for Estimating...
 Comparison of Methods for...
 Comparison of Human and...
 Discussion
 Program Availability and...
 Acknowledgements
 References
 
We thank Hinrich Schulenburg, the two referees Keith Crandall and Spencer Muse, and Associate Editor Caro-Beth Stewart for many constructive comments. We thank X. Xia for the analysis using the method of Li (1993)Citation . This study is supported by a BBSRC grant to Z.Y. and NSF grant DEB 9815367.


    Footnotes
 
Caro-Beth Stewart, Reviewing Editor

1 Keywords: synonymous rate, nonsynonymous rate, approximate methods, maximum likelihood, molecular evolution, adaptive evolution, positive selection. Back

2 Address for correspondence and reprints: Ziheng Yang, Department of Biology, 4 Stephenson Way, London NW1 2HE, England. E-mail: z.yang{at}ucl.ac.uk Back


    References
 TOP
 Abstract
 Introduction
 Methods for Estimating...
 Comparison of Methods for...
 Comparison of Human and...
 Discussion
 Program Availability and...
 Acknowledgements
 References
 

    Akashi, H. 1995. Inferring weak selection from patterns of polymorphism and divergence at "silent" sites in Drosophila DNA. Genetics 139:1067–1076.

    Comeron, J. M. 1995. A method for estimating the numbers of synonymous and nonsynonymous substitutions per site. J. Mol. Evol. 41:1152–1159.[ISI][Medline]

    Crandall, K. A., and D. M. Hillis. 1997. Rhodopsin evolution in the dark. Nature 387:667–668.

    Eyre-Walker, A., and P. D. Keightley. 1999. High genomic deleterious mutation rates in hominoids. Nature 397:344–347.

    Gillespie, J. H. 1991. The causes of molecular evolution. Oxford University Press, Oxford, England.

    Gojobori, T. 1983. Codon substitution in evolution and the "saturation" of synonymous changes. Genetics 105:1011–1027.

    Goldman, N. 1993. Statistical tests of models of DNA substitution. J. Mol. Evol. 36:182–198.[ISI][Medline]

    Goldman, N., and Z. Yang. 1994. A codon-based model of nucleotide substitution for protein-coding DNA sequences. Mol. Biol. Evol. 11:725–736.[Abstract/Free Full Text]

    Grimmett, G. R., and D. R. Stirzaker. 1992. Probability and random processes. 2nd edition. Clarendon Press, Oxford, England.

    Hasegawa, M., H. Kishino, and T. Yano. 1985. Dating the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22:160–174.[ISI][Medline]

    Ina, Y. 1995. New methods for estimating the numbers of synonymous and nonsynonymous substitutions. J. Mol. Evol. 40:190–226.[ISI][Medline]

    ———. 1996. Pattern of synonymous and nonsynonymous substitutions: an indicator of mechanisms of molecular evolution. J. Genet. 75:91–115.[ISI]

    Jukes, T. H., and C. R. Cantor. 1969. Evolution of protein molecules. Pp. 21–123 in H. N. Munro, ed. Mammalian protein metabolism. Academic Press, New York.

    Kimura, M. 1980. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J. Mol. Evol. 16:111–120.[ISI][Medline]

    ———. 1983. The neutral theory of molecular evolution. Cambridge University Press, Cambridge, England.

    Lee, Y. H., T. Ota, and V. D. Vacquier. 1995. Positive selection is a general phenomenon in the evolution of abalone sperm lysin. Mol. Biol. Evol. 12:231–238.[Abstract]

    Lewontin, R. 1989. Inferring the number of evolutionary events from DNA coding sequence differences. Mol. Biol. Evol. 6:15–32.[Abstract]

    Li, W.-H. 1993. Unbiased estimation of the rates of synonymous and nonsynonymous substitution. J. Mol. Evol. 36:96–99.[ISI][Medline]

    ———. 1997. Molecular evolution. Sinauer, Sunderland, Mass.

    Li, W.-H., C.-I. Wu, and C.-C. Luo. 1985. A new method for estimating synonymous and non-synonymous rates of nucleotide substitutions considering the relative likelihood of nucleotide and codon changes. Mol. Biol. Evol. 2:150–174.[Abstract]

    Messier, W., and C.-B. Stewart. 1997. Episodic adaptive evolution of primate lysozymes. Nature 385:151–154.

    Miyata, T., and T. Yasunaga. 1980. Molecular evolution of mRNA: a method for estimating evolutionary rates of synonymous and amino acid substitutions from homologous nucleotide sequences and its applications. J. Mol. Evol. 16:23–36.[ISI][Medline]

    Moriyama, E. N., and J. R. Powell. 1997. Synonymous substitution rates in Drosophila: mitochondrial versus nuclear genes. J. Mol. Evol. 45:378–391.[ISI][Medline]

    Muse, S. V. 1996. Estimating synonymous and nonsynonymous substitution rates. Mol. Biol. Evol. 13:105–114.[Abstract]

    Muse, S. V., and B. S. Gaut. 1994. A likelihood approach for comparing synonymous and nonsynonymous nucleotide substitution rates, with application to chloroplast genome. Mol. Biol. Evol. 11:715–724.[Abstract/Free Full Text]

    Nei, M., and T. Gojobori. 1986. Simple methods for estimating the number of synonymous and nonsynonymous nucleotide substitutions. Mol. Biol. Evol. 3:418–426.[Abstract]

    Nielsen, R., and Z. Yang. 1998. Likelihood models for detecting positively selected amino acid sites and applications to the HIV-1 envelope gene. Genetics 148:929–936.

    Ohta, T. 1995. Synonymous and nonsynonymous substitutions in mammalian genes and the nearly neutral theory. J. Mol. Evol. 4056–63.

    Ota, T., and M. Nei. 1994. Variance and covariances of the numbers of synonymous and nonsynonymous substitutions per site. Mol. Biol. Evol. 11:613–619.[Abstract]

    Pamilo, P., and N. O. Bianchi. 1993. Evolution of the Zfx and Zfy genes—rates and interdependence between the genes. Mol. Biol. Evol. 10:271–281.[Abstract]

    Pedersen, A.-M. K., C. Wiuf, and F. B. Christiansen. 1998. A codon-based model designed to describe lentiviral evolution. Mol. Biol. Evol. 15:1069–1081.[Abstract]

    Perler, F., A. Efstratiadis, P. Lomedica, W. Gilbert, R. Kolodner, and J. Dodgson. 1980. The evolution of genes: the chicken preproinsulin gene. Cell 20:555–566.

    Stuart, A., K. Ord, and S. Arnold. 1999. Kendall’s advanced theory of statistics. Vol. 2a, 6th edition. Arnold, London.

    Tateno, Y., N. Takezaki, and M. Nei. 1994. Relative efficiencies of the maximum-likelihood, neighbor-joining, and maximum-parsimony methods when substitution rate varies with site. Mol. Biol. Evol. 11:261–277.[Abstract]

    Yang, Z. 1994a. Estimating the pattern of nucleotide substitution. J. Mol. Evol. 39:105–111.

    ———. 1994b. Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods. J. Mol. Evol. 39:306–314.

    ———. 1998. Likelihood ratio tests for detecting positive selection and application to primate lysozyme evolution. Mol. Biol. Evol. 15:568–573.[Abstract]

    ———. 1999. Phylogenetic analysis by maximum likelihood (PAML). Version 2. University College London, England.

    Yang, Z., and R. Nielsen. 1998. Synonymous and nonsynonymous rate variation in nuclear genes of mammals. J. Mol. Evol. 46:409–418.[ISI][Medline]

    Yang, Z., R. Nielsen, and M. Hasegawa. 1998. Models of amino acid substitution and applications to mitochondrial protein evolution. Mol. Biol. Evol. 15:1600–1611.[Abstract/Free Full Text]

Accepted for publication September 6, 1999.