# Statistical Tests for Detecting Positive Selection by Utilizing High-Frequency Variants

^{*}State Key Laboratory of Biocontrol and Key Laboratory of Gene Engineering of the Ministry of Education, Sun Yat-sen University, Guangzhou 510275, China,^{†}Department of Ecology and Evolution, University of Chicago, Chicago, Illinois 60637,^{‡}Human Genetics Center, School of Public Health, University of Texas, Houston, Texas 77030 and^{§}Laboratory for Conservation and Utilization of Bioresources, Yunnan University, Kunming 650091, China

- 1
*Corresponding author:*1101 E. 57th St., Chicago, IL 60637. E-mail: kzeng{at}uchicago.edu

## Abstract

By comparing the low-, intermediate-, and high-frequency parts of the frequency spectrum, we gain information on the evolutionary forces that influence the pattern of polymorphism in population samples. We emphasize the high-frequency variants on which positive selection and negative (background) selection exhibit different effects. We propose a new estimator of θ (the product of effective population size and neutral mutation rate), θ* _{L}*, which is sensitive to the changes in high-frequency variants. The new θ

*allows us to revise Fay and Wu's*

_{L}*H*-test by normalization. To complement the existing statistics (the

*H*-test and Tajima's

*D*-test), we propose a new test,

*E*, which relies on the difference between θ

*and Watterson's θ*

_{L}_{W}. We show that this test is most powerful in detecting the recovery phase after the loss of genetic diversity, which includes the postselective sweep phase. The sensitivities of these tests to (or robustness against) background selection and demographic changes are also considered. Overall,

*D*and

*H*in combination can be most effective in detecting positive selection while being insensitive to other perturbations. We thus propose a joint test, referred to as the

*DH*test. Simulations indicate that

*DH*is indeed sensitive primarily to directional selection and no other driving forces.

DETECTING the footprint of positive selection is an important task in evolutionary genetic studies. Many statistical tests have been proposed for this purpose. Some of them use only divergence data between species (see Nei and Kumar 2000 and Yang 2003 for reviews) while others use both divergence and polymorphism data (for example, Hudson *et al*. 1987; McDonald and Kreitman 1991; Bustamante *et al*. 2002). Those that rely only on polymorphism data can be classified into either “haplotype” tests (for example, Hudson *et al*. 1994; Sabeti *et al*. 2002; Voight *et al*. 2006) or the site-by-site “frequency spectrum” tests (for example, Tajima 1989; Fu and Li 1993).

The frequency spectrum is the distribution of the proportion of sites where the mutant is at frequency *x* and is the focus of this study. We may divide the spectrum broadly into three parts: the low-, intermediate-, and high-frequency variant classes. Comparing the three parts of the spectrum, we shall have a fuller view of the configuration of polymorphisms. Tajima's (1989) *D* and Fu and Li's (1993) *D* ask whether there are too few or many more rare variants than common ones. Fay and Wu's (2000) *H* takes into consideration the abundance of very high-frequency variants relative to the intermediate-frequency ones. Thus far, there is no test that addresses the relative abundance of the very high- and very low-frequency classes, although this comparison might be the most informative about new mutations. After all, new mutations are most likely to be found in the very low class and least abundant in the very high class.

In this article, we study the power to detect selection by applying spectrum tests to a linked neutral locus. In the first section, we briefly summarize some classical results of the frequency spectrum theory. In the second section, we introduce a new estimator of θ (θ = 4*N*μ, where *N* is the effective population size and μ is the neutral mutation rate of the gene). On the basis of this new estimator, we revise Fay and Wu's *H*-test which was not normalized and introduce a new test (*E*) that contrasts the high- and low-frequency variants. We further propose a joint test, *DH*, which is a combination of *D* and *H*. In the rest of the article, we use computer simulation to compare the powers of the tests to detect positive selection. We also consider their sensitivities to balancing selection, background selection, and demographic factors.

## GENERAL BACKGROUND

Let ξ* _{i}* be the number of segregating sites where the mutant type occurs

*i*times in the sample. Following Fu (1995), we refer to the class of mutations with

*i*occurrence as mutations of size

*i*. Given the population mutation rate θ, Fu (1995) showed that(1)for a sample of size

*n*. Since each class of mutations contains information on θ, there can be many linear functions of ξ

*that are unbiased estimators of θ, depending on how each frequency class is weighted. A general form is(2)with*

_{i}*E*(η) = θ, where

*c*'s are weight constants. A few well-known examples are Watterson's (1975) θ

_{i}_{W}, Tajima's (1983) θ

_{π}, Fu and Li's (1993) ξ

*, and Fay and Wu's (2000) θ*

_{e}*:(3)where . Theoretical variances of these estimators are given in Figure 1. Among them, θ*

_{H}_{W}has the smallest variance.

Different estimators have varying sensitivities to changes in different parts of the frequency spectrum. For example, ξ* _{e}* and θ

_{W}are sensitive to changes in low-frequency variants, θ

_{π}to changes in intermediate-frequency variants, and θ

*to high-frequency ones. When the level of variation is influenced by different population genetic forces, different parts of the spectrum are affected to different extents (Slatkin and Hudson 1991; Fu 1996, 1997; Fay and Wu 2000; Griffiths 2003). The difference between two θ-estimators can thus be informative about such forces. For example, rapid population growth tends to affect low-frequency variants more than it affects high-frequency ones. As a result, θ*

_{H}_{W}tends to be larger than θ

_{π}. The first test to take advantage of the differences between estimators is Tajima's

*D*-test (Tajima 1989), as shown below:(4)Others have since been proposed (Fu and Li 1993; Fu 1996, 1997; Fay and Wu 2000). Among the evolutionary forces that may cause the frequency spectrum to deviate from the neutral equilibrium, hitchhiking (Maynard Smith and Haigh 1974) has attracted much attention. A salient feature of positive selection is the excess of high-frequency variants (Fay and Wu 2000).

It should be noted that θ_{W} and θ_{π} can also be written as(5)where *s* is the number of segregating sites and *d _{ij}* is the number of differences between the

*i*th and

*j*th sequences. Equation 5 suggests that θ

_{W}and θ

_{π}can be calculated without using an outgroup sequence to determine the mutant alleles. In contrast, ξ

*and θ*

_{e}*do require knowledge of an outgroup sequence. As consequences, Tajima's*

_{H}*D*can be calculated without using an outgroup sequence, while the

*H*-test needs one (see Fu and Li 1993 for a version of Fu and Li's

*D*-test that uses ξ

*but requires no outgroup sequence in the calculation).*

_{e}## ANALYTICAL RESULTS

#### Measurements of variation based on very high-frequency variants—θ_{H} and θ_{L}:

_{H}

_{L}

To capture the dynamics of high-frequency variants, we need to put most weight on high-frequency variants in the estimation of θ. Fay and Wu (2000) used θ* _{H}* (3). Its variance term, however, is not easy to obtain. Here we propose a new estimator, θ

*. The variance of θ*

_{L}*can be easily obtained and then used to calculate the variance of θ*

_{L}*.*

_{H}Consider a genealogy of *n* genes from a nonrecombining region. The mean number of mutations accumulated in each gene since the most recent common ancestor (MRCA) of the sample can be calculated as(6)Since *E*(ξ* _{i}*) = θ/

*i*, it follows that(7)Although is an asymptotically unbiased estimator of θ, it is more convenient to work with its canonical form(8)which has

*E*(θ

*) = θ. It can be shown that(9)where . Since*

_{L}*i*

^{2}=

*n*·

*i − i*(

*n − i*), it is easy to see that θ

*= 2θ*

_{H}*− θ*

_{L}_{π}. From this property, we obtain(10)(see appendix).

#### The normalized Fay and Wu's *H*-statistic—contrasting high- and intermediate-frequency variants:

Recall that Fay and Wu's (2000) *H* is defined as *H* = θ_{π} − θ* _{H}* = 2(θ

_{π}− θ

*). Now we can write the normalized*

_{L}*H*-statistic as(11)where(12)(see appendix). In practice, θ in (12) can be estimated by θ

_{W}, and θ

^{2}can be estimated by ) (Tajima 1989).

#### A new *E*-test—contrasting high- and low-frequency variants:

As mentioned before, there is no test contrasting the low- and high-frequency parts of the spectrum. Both *D* and *H* use intermediate-frequency variants as a benchmark for comparison. Nevertheless, there are occasions when it is informative to contrast high- and low-frequency variants. For example, when most variants are lost (*e.g*., due to selective sweep), the recovery of the neutral equilibrium is most rapid for the low-frequency ones and slowest for the high-frequency ones (*e.g*., Figure 2B). Taking advantage of the newly derived θ* _{L}*, we propose a new test statistic(13)where(14)(see appendix). In practice, θ and θ

^{2}in (14) can be estimated by methods described after (12). Thus, the three tests,

*D*,

*H*, and

*E*, contrast the three parts of the frequency spectrum, low-, intermediate-, and high-frequency variants, in a pairwise manner. By properly choosing a test statistic, we can have a better chance to detect a particular driving force (

*e.g*., Table 3).

#### The *DH* test—a joint *D* and *H* test:

Simulation results (*e.g*., Table 3) suggest that both *D* and *H* are powerful in detecting selection, but they are sensitive to different demographic factors and are affected by background selection to different degrees. Therefore, one may conjecture that *D* and *H* in combination may be sensitive only to selection and less sensitive to other forces. For example, the sensitivity of *D* to population expansion may be counterbalanced by the insensitivity of *H* to the same factor in the joint *DH* test. To carry out the joint *DH* test, we first define(15)where *X _{s}* is random samples (of size

*n*) under neutrality with

*s*(

*s*≥ 1) segregating sites,

*d*(

*X*) is the value of the

*D*-statistic of sample

*X*and

*h*(·) is that of the normalized

*H*-statistic, and

*d*and

_{p}*h*are critical values such that

_{p}*P*{

*d*(

*X*) ≤

_{s}*d*} =

_{p}*P*{

*h*(

*X*) ≤

_{s}*h*} =

_{p}*p*. In other words, for a given

*s*, when the significance levels of both

*D*and

*H*are

*p*, that of the joint

*DH*test is

*f*(

_{s}*p*). That is, for samples with

*s*segregating sites, the rejection region of the

*DH*test at the significance level of

*f*(

_{s}*p*) is {

*X*|

_{s}*d*(

*X*) ≤

_{s}*d*and

_{p}*h*(

*X*) ≤

_{s}*h*}. In practice, for a given value of

_{p}*p*,

*f*(

_{s}*p*) can be evaluated using neutral simulations with the number of segregating sites

*s*fixed (Hudson 1993). In our implementation, at least 50,000 rounds of simulations were carried out. To determine a value of

*p*(and consequently

*d*and

_{p}*h*) such that the significance level of

_{p}*DH*is, for example, 5%, we solved the equation

*f*(

_{s}*p*) = 0.05 using the bisection algorithm (Press

*et al*. 1992). We denote the solution as

*p**. We also denote

*d*and

_{p}*h*corresponding to

_{p}*p**as

*d**and

*h**, respectively. In Table 1, we give the

*p*-*,

*d*-*, and

*h**-values for three significance levels of

*DH*[

*i.e*.,

*f*(

_{s}*p**)], given various observed values of

*s*. In general,

*p**is close to three times the value of

*f*(

_{s}*p**).

## APPLICATIONS

Here we apply the tests discussed previously to data simulated under a variety of conditions for the purpose of evaluating their statistical powers. Such simulations have been done extensively in the past for studying existing tests (Braverman *et al*. 1995; Simonsen *et al*. 1995; Fu 1997; Przeworski 2002). Here all tests were one-sided and values falling into the lower 5% tail of the null distribution were considered significant. The reason for the one-sided test is that positive selection predictably shifts the pattern to one side of the null distribution. We report only the results of the normalized *H*-test, because it is always more powerful than the original one (varying from ∼1 to 8%). Details of the simulation algorithms and other related issues are described in the appendix.

#### The power to detect positive selection before and after fixation:

For the linked neutral locus, we define , where and are the expected numbers of segregating sites of size *i* under selection and neutrality, respectively (Fu 1997). Figure 2 shows the dynamics of *R*(*i*) before and after the fixation of an advantageous mutation. The reduction in intermediate-frequency variants and the accumulation of high-frequency variants start at the early stage of the sweep (Figure 2A, *f* = 0.3; *f* is the population frequency of the advantageous mutation). When the frequency of the advantageous mutation reaches 60% in the population, the spectrum is already highly deviant from neutrality, losing many of the intermediate-frequency variants (Figure 2A, *f* = 0.6). By the time of fixation, the population has lost virtually all the intermediate-frequency variants, and mutations appear at either very high or very low frequency (Figure 2B, τ = 0; τ is time after fixation of the advantageous mutation and is measured in units of 4*N* generations). In neutral equilibrium, there should be 22.4 mutations in the sample on average for θ = 5 and *n* = 50. However, only 11.4 mutations, on average, are seen in the samples at τ = 0 (Figure 2B), among which 57% have a frequency <10%, and 31% have a frequency >80%. Strikingly, the mean number of mutations of size 49 is 7.2 times that of the neutral expectation. When τ > 0, these high-frequency variants quickly drift to fixation and no longer contribute to polymorphism. For a long period after fixation, fewer-than-expected numbers of intermediate- and high-frequency mutations are observed (Figure 2B).

Figure 3 shows how the tests behave before and after the fixation of the advantageous mutation. Before fixation (Figure 3, A and B, left side), the power of *D*, *H*, and *DH* increases rapidly as the frequency of the advantageous mutation increases. *H* is generally more powerful than *D* in this stage. The behavior of *DH* resembles that of *H*, but note that it is often the most powerful test when the advantageous allele is in high frequency. We should note again that the purpose of the *DH* test is not to increase the power in detecting positive selection over either *D* or *H* but to decrease sensitivity to demographical factors and background selection. This is apparent in later sections. In contrast, the *E*-test has little power. This is expected as the reduction of moderate-frequency variants and the accumulation of high-frequency variants characterize this stage (Figure 2A). Sometimes, *H*, *D*, and *DH* reach their peak of power before fixation. This happens when too much variation is removed by the selective sweep, as is the case in Figure 3A. After fixation (Figure 3, A and B, right side) the *E*-test quickly becomes the most powerful test (τ ≈ 0.1 in Figure 3A; τ ≈ 0.15 in Figure 3B). From Figure 2B, we can see that after fixation the low-frequency part of the spectrum recovers first, and the high-frequency part returns to its equilibrium level last. The result of the *E*-test fits this observation. Furthermore, since the recovery phase is much longer than the selective phase when selection is strong, the *E* test can be useful for detecting sweep. The *D*-test performs reasonably well in the recovery stage. The power of the *H*-test, however, decreases quickly. As a consequence, the power of *DH* also decreases but at a much lower rate than that of *H*.

#### The power to detect balancing selection:

Intuition would suggest that the results of an incomplete sweep (Figure 3) may also apply to the initial stage of balancing selection (*i.e*., a relatively short time period after the birth of the selected allele) under which the selected allele would reach an equilibrium frequency that is high, but <100%. Furthermore, as the selected allele does not reach fixation, the hitchhiking variants also may not have extreme frequencies. As a result, the signature of positive selection and hitchhiking may linger longer under balancing selection.

Simulation results (see supplemental Figure S1 at http://www.genetics.org/supplemental/) indeed suggest that in the phase when the selected allele is approaching the equilibrium frequency (75% in this case), the tests behave very much like the patterns in the left side of Figure 3. After the equilibrium is reached, the powers of *H* and *DH* also remain higher than in the right side of Figure 3 for a longer period of time. Nevertheless, the gain is modest as the power to detect selection generally diminishes when the selected allele stops increasing in frequency.

#### Sensitivity of the tests to other driving forces:

A good test of a particular population genetic force, say, positive selection, should be sensitive to that force and to that force only. Although a test that is sensitive to many different factors can be a useful general tool, it is ultimately uninformative about the true underlying force. Hence, power can in fact be a blessing in disguise. We examine the sensitivity of these tests to forces other than positive selection and balancing selection below.

#### Background selection:

Selection against linked deleterious mutations maintained by recurrent mutation, often referred to as background selection, can have effects on the level of genetic diversity similar to those of a selective sweep (Charlesworth *et al*. 1993). The distinction between selective sweep and background selection is therefore crucial. Fortunately, the two modes of selection often have very different effects on the frequency spectrum (Fu 1997). For example, background selection is not likely to have any effect on high-frequency variants (Fu 1997) and its effect on the low-frequency variants depends on *U* and *N*, where *U* is the deleterious mutation rate per diploid genome and *N* is the effective population size.

Table 2 summarizes the power of the tests to detect background selection. Both *D* and *E* are sensitive, but *D* is slightly more powerful than *E*. In general, the power increases as *U* increases, but this increase also depends on *N*. For a given *U*, the power of the tests decreases as *N* increases. *H* is not affected in all cases and hence is discriminatory between the two modes of selection. Significantly, *DH* is also not sensitive to background selection. It seems that the sensitivity of *D* is well counterbalanced by the insensitivity of *H*.

#### Population growth:

When a population increases in size, it tends to have an excess of low-frequency variants (Slatkin and Hudson 1991; Fu 1997; Griffiths 2003). Both *D* and *E* are sensitive to this type of deviation (Figure 4). *E* is the most sensitive test because high-frequency variants are the last to reach the new equilibrium after expansion. In contrast, *H* is unaffected. Strikingly, *DH* has little power in this case even though *D* exhibits tremendous power.

#### Population shrinkage:

When the population decreases in size, the number of low-frequency variants tends to be smaller than that of the intermediate- and high-frequency ones (Fu 1996). Thus, *H* can be sensitive to population shrinkage, whereas *D*, *E*, and *DH* are largely unaffected (Figure 5).

#### Population subdivision:

When the population is structured, one commonly uses the fixation index *F*_{ST} as a measure of genetic differences among subpopulations. In the symmetric two-deme model, *F*_{ST} = 1/(1 + 16*Nm*), where *m* is the fraction of new migrants in the population and *N* is the population size of a deme (Nordborg 1997). When there is strong population subdivision (*F*_{ST} > 0.2 or 4*Nm* < 1, for example) and samples do not come equally from all subpopulations, the frequency spectrum is deviant from the neutral equilibrium (results not shown). Figure 6A gives an example where all genes are taken from one subpopulation. In this case, *H* is the most sensitive test. The *D*-test is sensitive only when the subdivision is strong (4*Nm* < 0.5 or *F*_{ST} > 0.33). *DH* has some sensitivity. Its power lies in between those of *D* and *H*, but never rises above 15%. Note that, for relatively mild population subdivisions with 4*Nm* > 2, none of these four tests is notably affected. Interestingly, *E* is completely insensitive to population subdivision. Figure 6B shows the effect of sampling on the power of the tests (4*Nm* = 0.1). The *E*-test again shows no sensitivity. Although the power of *H* can be as high as 55%, that of the *DH* test never goes above 14%. This result shows the effectiveness of using *D* and *H* in combination. However, when sampling becomes less and less biased, all tests become progressively less sensitive to population structure.

## CONCLUSION

We summarize our results in Table 3. It is clear that, among the four tests studied, *DH* is sensitive primarily to directional selection. In contrast, other tests are sensitive to two or more other driving forces. *D* is indeed a “general purpose” test as it is sensitive, to various extents, to all the processes we considered (Table 3). Being exclusive is a desirable property because it minimizes false positive results in the search for selectively favored mutations. However, no test is powerful in every stage of the selection process. Therefore, in practice, we should utilize *a priori* information (if available) from other sources to help us choose an appropriate test.

## APPENDIX

#### Analytical results:

We list some useful properties of the θ-estimators (see Y.-X. Fu, unpublished data, for mathematical details):(A1)Since , using (A1), we have(A2) and can be calculated using the above formulas.

#### Simulation algorithms:

##### Positive selection:

We used the coalescent process with a selective phase (Kaplan *et al*. 1989; Stephan *et al*. 1992; Braverman *et al*. 1995) to simulate random samples and estimate the power of various tests. The model assumes that selection at the favored locus is additive and the selective pressure is strong (2*Ns* ≫ 1). In this case, the trajectory of the frequency of the advantageous allele in the interval with endpoints ϵ and 1 − ϵ is approximately deterministic (Stephan *et al*. 1992). We used ϵ = 1/(2*N*) in this study. Note that the choice of ϵ has little effect on the results (Braverman *et al*. 1995). In fact, the results were identical when we ran simulations with ϵ = 100/(2*N*) or ϵ = 1/(4*N*). Our implementation of the algorithm followed the description in Braverman *et al*. (1995), except for the following features: (1) we used a given value of the population mutation rate (θ = 4*N*μ), rather than a fixed number of segregating sites, (2) we allowed intragenic recombination within the neutral locus to happen through the whole process (*i.e*., during both neutral and selective phases), (3) we assumed that only one round of selective sweep had occurred in the history of a sample rather than recurrent sweeps, and (4) the recombination distance between the selective locus and the neutral locus (*c*) was a given value rather than a random variable.

We extended the coalescent simulation algorithm described above to incorporate the case of incomplete hitchhiking. Suppose that the frequency of the favored allele is *f*. By solving the deterministic equation given by Stephan *et al*. (1992; Equation 1 in Braverman *et al*. 1995), one may obtain the time since the birth of this advantageous mutation and the corresponding frequency trajectory. In each replicate of simulation, we assigned the number of genes linked to the advantageous allele by binomial sampling with parameters *n* and *f*. Then the genealogy was constructed following the same procedure as described above.

##### Background selection and demographic scenarios:

We simulated background selection using the coalescent method (Hudson and Kaplan 1994; Charlesworth *et al*. 1995). The software package kindly provided by Hudson (2002) was used to simulate all demographic scenarios.

#### Critical values of the tests:

The tests, *D*, *H*, and *E*, were one-sided; values falling into the lower 5% tail were considered significant. The critical values of these three tests were determined by the Monte Carlo method used by Fu (1997). For a sample of size *n*, we first obtained an estimate of using Watterson's estimator. Then we ran 100,000 rounds of coalescent simulations with no recombination and . The critical values were determined by examining the empirical distribution obtained. The supplemental table (at http://www.genetics.org/supplemental/) lists the critical values of the *E*-test for a number of θ*-*values. Another method to determine critical values was proposed by Hudson (1993). The differences between these two methods have been discussed extensively (Markovtsova *et al*. 2001; Wall and Hudson 2001). The results produced by both methods are quite similar; therefore, we report only those generated by the first method. The critical values of the *DH* test were generated using the second method.

## Acknowledgments

We thank J. Braverman and R. R. Hudson for their helpful discussions about the simulation algorithms. We also thank two anonymous reviewers for their constructive comments. K.Z. is supported by Sun Yat-sen University and the Kaisi Fund. S.S. is supported by grants from the National Natural Science Foundation of China (30230030, 30470119, 30300033, and 30500049). Y.-X. Fu is supported by National Institutes of Health (NIH) grants (GM 60777 and GM50428) and funding from Yunnan University, China. C.-I Wu is supported by NIH grants and an OOCS grant from the Chinese Academy of Sciences.

## Footnotes

Communicating editor: N. Takahata

- Received May 30, 2006.
- Accepted August 16, 2006.

- Copyright © 2006 by the Genetics Society of America