## Abstract

Several studies have reported a negative correlation between estimates of the nonsynonymous to synonymous rate ratio (*ω* = *d*_{N}/*d*_{S}) and the sequence distance *d* in pairwise comparisons of the same gene from different species. That is, more divergent sequences produce smaller estimates of *ω*. Explanations for this negative correlation have included segregating nonsynonymous polymorphisms in closely related species and nonlinear dynamics of the ratio of two random variables. Here we study the statistical properties of the maximum-likelihood estimates of *ω* and *d* in pairwise alignments and explore the possibility that the negative correlation can be entirely explained by those properties. We show that the *ω* estimate is positively biased for small *d* and that the bias decreases with the increase of *d*. We also show that the estimates of *ω* and *d* are negatively correlated when *ω* < 1 and positively correlated when *ω* > 1. However, the bias in estimates of *ω* and the correlation between estimates of *ω* and *d* are not enough to explain the much stronger correlation observed in real data sets. We then explore the behavior of the estimates when the model is misspecified and suggest that the observed correlation may be due to protein-level selection that causes very different amino acids to be favored in different domains of the protein. Widely used models fail to account for such among-site heterogeneity and cause underestimation of the nonsynonymous rate and *ω*, with the bias being much stronger for distant sequences. We point out that tests of positive selection based on the *ω* ratio are invariant to the parameterization of the model and thus unaffected by bias in the *ω* estimates or the correlation between estimates of *ω* and *d*.

NATURAL selection acts differently on synonymous and nonsynonymous mutations in protein-coding genes (Kimura 1977). For almost all proteins, the nonsynonymous distance (*d*_{N}, measured by the number of nonsynonymous substitutions per nonsynonymous site) is smaller than the corresponding synonymous distance (*d*_{S}, measured by the number of synonymous substitutions per synonymous site), because the structure and function of a protein impose constraints on the types of nonsynonymous (amino acid) substitutions that can take place, while synonymous substitutions usually accumulate freely with little or no interference from selection. In contrast, a large nonsynonymous distance (*d*_{N} > *d*_{S}) is a strong indicator of positive (diversifying) selection acting on the protein (*e.g.*, Messier and Stewart 1997). The nonsynonymous to synonymous rate ratio, *ω* = *d*_{N}/*d*_{S}, in a protein-coding gene is thus an important measure of the mode and strength of selection acting on the protein, and it is widely used in detecting positive Darwinian selection [indicated by *ω* > 1 (Yang and Bielawski 2000)] and also as a descriptive statistic calculated from pairwise sequence alignments.

A number of recent studies have reported a negative correlation between estimates of *ω* and the evolutionary distance *d* (the total distance separating the sequences, including nonsynonymous and synonymous changes) in pairwise comparisons of the sequences for the same gene from various organisms (Rocha *et al.* 2006; Peterson and Masel 2009; Wolf *et al.* 2009). For example, Rocha *et al.* (2006) found that pairwise estimates of *ω* are larger among closely related bacterial genomes. They proposed that deleterious nonsynonymous polymorphisms explain the larger estimates of *ω* between closely related sequences. As the time of divergence between species increases, nonsynonymous polymorphisms are lost due to purifying selection and the *ω* ratio decays as a function of time. Similarly, Wolf *et al.* (2009) reported a negative correlation between estimates of *ω* and *d* in comparisons among several mammalian and avian genomes. They argued that this correlation is the result of being the ratio of two random variables (*d*_{N} and *d*_{S}) and suggested that current statistical tests based on *ω* (for example, likelihood-ratio tests of positive selection comparing the null hypothesis *ω* = 1 against the alternative hypothesis *ω* > 1) are inappropriate in light of the nonlinear behavior of the ratio of random variables.

The observed pattern of decay of as a function of seems to be common (Rocha *et al.* 2006; Peterson and Masel 2009; Wolf *et al.* 2009). For example, Figure 1A shows pairwise estimates of *ω* and *d* for 12 mitochondrial protein-coding genes from 244 mammal species and Figure 1A′ shows a similar plot for 18 ribosomal protein-coding genes from 49 bacterial species. In both cases, estimates of *ω* are clearly higher for closely related species, and the estimates decay as a function of the estimated distance. Figure 1, B and B′, shows the estimated synonymous distance *d*_{S} *vs.* estimated *d* and Figure 1, C and C′, shows the estimated nonsynonymous distance *d*_{N} *vs.* estimated *d*. Note that estimates of *d*_{S} and *d* are highly correlated in both cases, while estimates of *d*_{N} seem saturated for large *d*.

Pairwise methods ignore the phylogenetic relationships of the species being studied and do not make as efficient use of the information in the sequences as a joint analysis of all sequences on a phylogeny. Nonetheless, pairwise estimates have been widely used, so it is important to understand their limitations. In this article we study the statistical properties of estimates of *ω* in pairwise sequence alignments, using a combination of techniques, such as mathematical analysis, computer simulation, and real data analysis. We derive the biases in estimates of *ω* and *d*, using a fictitious regularized genetic code, and demonstrate similar biases under the universal genetic code, using computer simulation. We then study the impact of heterogeneous selection for amino acids along the protein sequence on estimation of *ω* and *d* by computer simulation. Our analysis suggests that the small-sample biases in estimates of *ω* as well as systematic biases in *ω* estimates caused by model violation may be important factors to explain the correlation between and observed in empirical data sets.

## Asymptotic Bias and Correlation in Estimates of *ω* and *d*

We first consider a simple model of codon evolution that can be treated analytically. Consider a regularized genetic code with 16 amino acids and with every codon fourfold degenerate (Table 1). Substitutions at first and second codon positions are always nonsynonymous and those at the third positions are always synonymous. There are no stop codons. The substitution rate from codon *i* to *j* (*i* ≠ *j*) is (1)The *q _{ij}* are the off-diagonal elements of a 64 × 64 rate matrix,

*Q*, with the diagonal elements,

*q*, set such that each row of

_{ii}*Q*sums to one. Time is measured as the expected number of substitutions per codon site; in other words

*μ*= 1/(6

*ω*+ 3) so that . Because every codon position is always either synonymous or nonsynonymous, evolution of codons can be treated as a mixture of two Jukes–Cantor (Jukes and Cantor 1969) models: the first two positions evolve with rate

*ωμ*and the third position with rate

*μ*. Consider a pairwise sequence alignment with

*n*codons and with

*x*

_{S}and

*x*

_{N}observed synonymous and nonsynonymous differences. Let

*d*

_{S}and

*d*

_{N}be the distances at synonymous and nonsynonymous sites, respectively. The total distance between the two sequences is thus

*d*=

*d*

_{S}+ 2

*d*

_{N}. It follows that

*ω*=

*d*

_{N}/

*d*

_{S}, and

*d*= (1 + 2

*ω*)

*d*

_{S}. The log-likelihood function is given by the fact that evolution at the first two codon positions and at the third codon position is independent, (2)where(3)(4)The maximum-likelihood estimates (MLEs) of the distances at synonymous and nonsynonymous sites are, respectively,(5)(6)where and are the observed proportions of synonymous and nonsynonymous differences, respectively. Note that and , where

*E*denotes expectation. The MLEs of

*d*and

*ω*are (7)(in expected substitutions per codon) and(8)There are only two parameters in the model, and we can work either with

*d*and

*ω*or as a different parameterization with

*p*

_{S}and

*p*

_{N}.

Both and are functions of the binomial proportions and that have variances(9)(10)and . Because and are functions of random variables with known variances, the delta method can be used to calculate the asymptotic expectations(11)(12)and the asymptotic variance/covariance matrix(13)where is the Jacobian matrix of the transform (reparameterization). The Jacobian and the second derivatives of and with respect to and are given in the *Appendix*.

From Equations 11 and 12 the bias of the estimates can be approximated by and . Note that when *n*→∞, and , as expected according to standard theory. However, when *n* is small, the bias in the estimators may be substantial. Figure 2 shows the biases of and for various cases. In particular, there is a positive bias in for small *d* (Figure 2A). On the other hand, there is very little bias in (Figure 2B). Simulations confirm that the approximations (Equation 11 and Equation 12) are very good with ≥500 codons in the alignment (Figure 2, A and B). From Equation 13 the asymptotic correlation can be calculated as . Figure 3 shows the asymptotic correlation for various values of *ω* and *d*. Calculations indicate that when *ω* < 1, the correlation is negative (*ρ* < 0), and when *ω* > 1, it is positive (*ρ* > 1).

The codon substitution model based on the standard genetic code is not tractable. However, we can study the asymptotic properties of the model by simulation. We simulated pairwise sequence alignments of length *n* = 500 codons, with equal codon frequencies (1/61) and with *κ* = 2. We simulated 1000 pairwise alignments for each combination of *d* = 0.01, 0.1, 1, 10 and *ω* = 0.02, 0.2, 1, 2, 10 (1000 × 4 × 5 = 20,000 data sets). Estimates of *ω* and *d* for each data set were then obtained by maximum likelihood (ML). The simulation shows the same pattern of correlation (Figure 4) as in the case of the regularized code.

To assess whether the bias and correlation of the MLEs of *d* and *ω* may be responsible for the strong correlation observed in real data sets, we simulated a 244-species sequence alignment on the mammal phylogeny, using the mitochondrial genetic code. The values of *n*, *ω*, and *κ* and the branch lengths in the simulation were fixed to the MLEs obtained from the joint analysis of the original sequence data. The codon frequencies were fixed to those observed in the real data. Figure 5, A–C, shows the pairwise estimates for this simulation. Although there is a negative correlation between estimates of *ω* and *d*, this correlation is much weaker than that observed for the real data (Figure 1A). Furthermore, when the simulated alignment is much longer (when we set *n* = 10^{6} codons), the correlation almost disappears and the MLE of *ω* converges to the true value (not shown). The simulation model assumes the same *ω* for all sites in the gene, which is not realistic for real genes. We also simulated with a different model where *ω* varies among sites according to a mixture of three site classes [M3 discrete (Yang *et al.* 2000)]. For the real mitochondrial data under this model, the MLEs are , , and , with the estimated frequencies of the site classes being , , and ; and . We then used those parameter values to simulate a 244-species alignment on the mammal phylogeny. Figure 5, A′–C′, shows the pairwise estimates of *ω* and *d* (obtained using a model of a single *ω* for all sites) from this simulated data set. The correlation between and is stronger, but in general much weaker than for the real data (Figure 1A). Therefore, the asymptotic pattern of correlation and bias seems too weak to explain the pattern observed in real data (Figure 1A).

## Variable Selection on Different Domains of the Protein

A major characteristic of protein evolution that has been ignored in the discussion above is among-site heterogeneity. In particular, structural characteristics of real proteins impose constraints on the frequencies of amino acids that vary hugely between protein domains. For example, the core of a protein is composed mainly of hydrophobic amino acids, while the surface of the protein may be composed mainly of hydrophilic amino acids (Baud and Karlin 1999). The active site of an enzyme may tolerate only very few different amino acids that can stabilize a particular substrate and carry out an enzymatic reaction. Halpern and Bruno (1998; see also Tamuri *et al.* 2012) proposed a codon substitution model based on a population genetics model of site-specific amino acid frequencies. To avoid the use of too many parameters we consider a model of site classes. Suppose there are *C* site classes. If the protein has *n* codons, there are *n*/*C* codons for each site class. The substitution rate from codon *i* to *j* in site class *k* is given by (14)where *F _{i}*

_{,}

*is the fitness of the amino acid encoded by codon*

_{k}*i*in site class

*k*, and (15)is the relative fixation probability of an

*i*to

*j*mutation relative to a neutral mutation (Kimura 1983, p. 45; see also Yang and Nielsen 2008), where

*h*(

*F*

_{j}_{,}

*−*

_{k}*F*

_{i}_{,}

*) > 1, = 1, < 1 for*

_{k}*F*

_{j}_{,}

*>*

_{k}*F*

_{i}_{,}

*,*

_{k}*F*

_{j}_{,}

*=*

_{k}*F*

_{i}_{,}

*, and*

_{k}*F*

_{j}_{,}

*<*

_{k}*F*

_{i}_{,}

*, respectively. In other words, selection accelerates or slows down substitutions compared to the rate for a neutral mutation (*

_{k}*i.e.*, when

*F*

_{j}_{,}

*=*

_{k}*F*

_{i}_{,}

*). The equilibrium frequency of a given amino acid*

_{k}*j*in site class

*k*is given by(16)If

*F*

_{j}_{,}

*= −∞, then*

_{k}*π*

_{j}_{,}

*= 0 and the amino acid is not observed in site class*

_{k}*k*. The same value of

*μ*is used to scale all rate matrices for the site classes, and

*μ*is set so that time is measured as the expected number of substitutions per site across the gene (

*i.e.*, so that ).

## The Regularized Genetic Code

We examine whether heterogeneity among sites in amino acid preference may lead to strong correlation between estimates of *ω* and *d* under the model of Equation 1 when the data are generated by the process of Equation 14. First we study a simple case under the regularized genetic code (Table 1). There are four classes of sites in the protein, each class corresponding to a group of four amino acids whose codons differ only at the second position. For example, at sites of class I, we observe only amino acids A–D (Table 1); that is, A–D have the same fitness but the other 12 amino acids have fitness −∞. Similarly, class II is composed of amino acids E–H, class III of I–L, and class IV of M–P. The gene is composed of equal proportions of the four site classes (that is, for a gene with 100 codons, 25 belong to class I, 25 to class II, and so on). Note that the substitution rate from codon *i* to *j* at a site class can be written as (17)where *μ* = 1/(3*ω* + 1).

We note the following about this substitution model:

the

*ω*ratio at the first codon position is*ω*_{1}= 0 and it is*ω*_{2}= 1 at the second. Therefore, the average*ω*across the protein is*ω*= 0.5.Gene-wide, the evolutionary rate is zero at the first codon position and

*μ*at the second and third positions. Therefore, the process can be described by a mixture of two independent Jukes–Cantor processes at the second and third positions.

Let *x*_{N} and *x*_{S} be the numbers of nonsynonymous and synonymous differences between two sequences of length *n* codons. The estimates of the nonsynonymous and synonymous distances are(18)(19)where , instead of as in the model of Equation 1. Then (20)and(21)The factor in the calculation of is necessary to correct for the zero substitution rate at the first codon position. Let us now assume that *ω*, *d*, *d*_{N}, and *d*_{S} are estimated using the model of Equation 1. We note three points:

The synonymous distance

*d*_{S}is correctly estimated.The nonsynonymous distance

*d*_{N}is underestimated and the underestimation is greater for larger sequence divergence (Figure 6A). In fact, the estimate of*d*_{N}approaches , a constant, as*d*→∞. The total distance*d*is also underestimated.As a consequence,

*ω*is underestimated more seriously for longer distances (Figure 6B), and as*d*→∞, . The correct value of*ω*is recovered only if the distance is close to zero:*d*→ 0, .

## The Standard Genetic Code

We perform a simulation study, using the standard genetic code. We use two site classes: a class of buried sites in the hydrophobic core of the protein and a class of exposed sites on the protein’s surface. We perform two simulations. The fitnesses for the 20 amino acids for the two site classes are calculated from the literature (Baud and Karlin 1999, Table 2). The substitution model is given by Equation 14, with *ω* = 0.04 and *n* = 4000 sites with 2000 sites for each site class. We simulated pairwise sequence alignments for different values of *d* (= 0.1, 0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4, 5, 7, and 10). The simulation code was written in R. The simulated alignments of two sequences were then analyzed using CODEML to estimate *ω* and *d*. The results are shown in Figure 7, A–C. It is clear that *ω* is underestimated with progressively long distances, causing a negative correlation between estimates of *ω* and *d*.

The amino acid frequencies of Table 2 are averaged across different proteins (Baud and Karlin 1999) and may not be as extreme as in individual proteins. To mimic extreme amino acid frequencies in specific proteins, we conducted a second simulation in which the fitnesses for both site classes are multiplied by 10. The results are shown in Figure 7, A′–C′. The underestimation of *ω* and the negative correlation between estimates of *ω* and *d* are more pronounced in this case.

## Discussion

Although statistical theory establishes that MLEs are asymptotically unbiased, MLEs often involve substantial biases in small samples. Furthermore, the correlation between the MLEs of two (or more) parameters estimated from the data must be studied for each model of interest. Pairwise estimates of *ω* and *d* exemplify both bias and correlation between MLEs. More problematic, though, is the bias in *ω* when the inference model is misspecified and ignores heterogeneity among sites in amino acid preferences. It appears that the small-sample bias of combined with the heterogeneous amino acid preferences may explain much of the negative correlation between estimates of *ω* and *d*. Models that account for heterogeneity in the substitution process among codons seem necessary to achieve appropriate estimates of *ω*. Bao *et al.* (2007, 2008) have developed fixed- and random-effects models that account for heterogeneity in amino acid frequencies in proteins. However, the models have not been applied to estimation of *ω* in pairwise alignments.

Rocha *et al.* (2006) considered the comparison of two sequences from the same species or from two closely related bacterial genomes and suggested that the decay of *ω* with the time of divergence may be due to deleterious nonsynonymous mutations in the process of being removed by natural selection. Sequences that diverged a long time ago will have more nonsynonymous mutations removed than sequences that diverged a short time ago, thus generating a negative correlation between *ω* and divergence time (or distance). However, their simulation model assumed that initially the population has no mutations, which seems biologically unrealistic. It should be preferable to analyze an equilibrium model, in which at any point in time, there is a distribution of nonsynonymous and synonymous mutations of different ages and frequencies, with the ratio of the two types of polymorphisms being stationary (Sawyer and Hartl 1992: Table 1). In contrast, the biases and correlations studied in this article apply to all levels of sequence divergence and should be considered even in comparisons of closely related sequences. The waiting time for the removal of nonsynonymous mutations discussed by Rocha *et al.* (2006) may be important in comparisons of closely related species, but may not be important in comparisons of distantly related species, such as different mammal species (Figure 1). In any case, at this stage we cannot exclude the possibility that population-level processes may affect the behavior of *ω* at small timescales, and further study is necessary.

It should be noted that the behavior of the MLEs of *ω* and *d* can be understood using standard maximum-likelihood theory, and the likelihood-ratio test of positive selection based on the *ω* ratio is not affected either by the bias in or by the correlation between and . This is because the LRT is invariant to reparameterization whereas bias and correlation are properties that arise and disappear by transformation. For example, the false-positive rates for the test of H_{0}: *ω* = 1 against H_{1}: *ω* ≠ 1 are 6.7%, 4.4%, 4.5%, and 6.1% for *d* = 0.01, 0.1, 1, and 10 in Figure 4, close to the 5% significance level expected according to standard ML theory. The suggestion by Wolf *et al.* (2009) that the LRT cannot control for the dependency between and is incorrect. The pattern of correlation observed in real data (*e.g.*, Figure 1A) may not be fully explained by the statistical properties of and , and some biological variables may play an important role. Nonetheless, we suggest the statistical properties of estimators must be studied carefully for each case of interest, before biological explanations for the observed patterns are offered.

## Acknowledgment

Z.Y. is a Royal Society/Wolfson Merit Award holder.

## Appendix

The Jacobian of Equation 13 is

(A1)The second derivatives of Equations 11 and 12 are(A2)(A3)(A4)(A5)When computing Equations 11–13, the Jacobian and second derivatives are evaluated at the true values and (Equations 3 and 4).

## Footnotes

*Communicating editor: M. A. Beaumont*

- Received April 5, 2013.
- Accepted June 17, 2013.

- Copyright © 2013 by the Genetics Society of America