# Detecting Selection Using Time-Series Data of Allele Frequencies with Multiple Independent Reference Loci

- Jo Nishino
^{1}

- Center for Information Biology and DNA Data Bank of Japan, National Institute of Genetics, Research Organization of Information and Systems, Mishima, Shizuoka 411-8540, Japan

- 1Corresponding author: Center for Information Biology and DNA Data Bank of Japan, National Institute of Genetics, Research Organization of Information and Systems, 1111 Yata, Mishima, Shizuoka 411-8540, Japan. E-mail: jnishino{at}nig.ac.jp

## Abstract

Recently, in 2013 Feder *et al.* proposed the frequency increment test (FIT), which evaluates natural selection at a single diallelic locus by the use of time-series data of allele frequencies. This test is unbiased under conditions of constant population size and no sampling noise. Here, we expand upon the FIT by introducing a test that explicitly allows for changes in population size by using information from independent reference loci. Various demographic models suggest that our proposed test is unbiased irrespective of fluctuations in population size when sampling noise can be ignored and that it has greater power to detect selection than the FIT if sufficient reference loci are used.

- selection
- neutrality
- time-series data of allele frequencies
- reference loci
- population size fluctuations

In population genetics, most data are obtained from a single point in time. When genetic time-series data are available, the use of such data to detect and estimate natural selection is an attractive concept. Time-series data may have direct information about natural selection because they affect allele frequencies in time. Bollback *et al.* (2008) introduced a statistical framework for estimating and testing natural selection by using time-series data of allele frequencies at a single diallelic locus. These authors applied their framework to an ancient human DNA sequence (Hummel *et al.* 2005) and a sample from an experimental evolution study of the bacteriophage, MS2 (Bollback and Huelsenbeck 2007).

Recent advances in high-throughput sequencing technology, including pooled DNA sequencing, have facilitated the acquisition of time-series data, and Bollback *et al.*’s (2008) method has been extended to more complicated situations. Mathieson and McVean (2013) applied a lattice model of population subdivision that enabled joint estimation of migration rate and spatially varying selection coefficients. The allele age also is an important parameter in selection inferences because it can provide information regarding the origin of a particular phenotype associated with the allele. Malaspinas *et al.* (2012) developed a method to estimate the selection coefficient and the allele age simultaneously. In addition, there has been increasing numbers of studies in which researchers focus on specific settings in each evolutionary experiment (*e.g.*, Illingworth and Mustonen 2011; Gallet *et al.* 2012; Illingworth *et al.* 2012).

Feder *et al.* (2013) reported recently that Bollback *et al.*’s (2008) standard χ^{2}-based test for selection is biased for realistic data with few sampled time points. When the number of sampled time points is sufficiently large, the likelihood ratio statistic (LRS) follows a χ^{2} distribution. However, the actual number of sampled time points rarely exceeds a few dozen. Particularly, when the null hypothesis is composite and the profile likelihood is used, the estimation of nuisance parameters can substantially bias inferences of the parameters of interest (*e.g.*, see Chapter 10 of Pawitan 2001). For the problem described in this report, the nuisance parameter is the population size.

To avoid bias, Feder *et al.* (2013) proposed two methods that both were modeled under conditions of constant population size and no sampling noise. In the empirical likelihood ratio test (ELRT), the population size is preliminarily estimated under neutrality as a first approximation. The estimated population size then is used to generate the empirical distribution of the LRS by computer simulation. Neutrality then can be evaluated by comparing the observed LRS with the empirical distribution. Although the ELRT was shown to be unbiased, this approach can be computationally intensive. To reduce the computational load, Feder *et al.* (2013) proposed the frequency increment test (FIT). The statistic used for FIT is defined as the following: Let be the population frequencies of one allele at a diallelic locus of interest at the sampled time, . The sampling time scales are short compared to the population size. Then the standardized allele frequency increment, (1)is approximately normally distributed with mean 0 under the null model, that is, neutral evolution. The variance of *Y _{i}* is equal to 1/(2

*N*) in the Wright−Fisher model with

*N*diploids. However, the variance is unknown because

*N*is unknown. In such a situation, natural selection can be evaluated by lettingwhere and

*S*are the sample mean and variance of

*Y*, respectively. We then perform a

_{i}*t*-test using the fact that

*t*follows the Student’s

_{FIT}*t*distribution with

*L*− 1 degrees of freedom under the null model. This is the FIT.

The FIT treats the nuisance parameter, *N*, as an unknown parameter instead of estimating it. When *Y _{i}* for any

*i*follows the normal distribution with the same variance under the null model, the FIT is an exact and unbiased test. Feder

*et al.*(2013) verified that actual type I error rates approach the nominal significance level for various parameter settings. These investigators also demonstrated that the power of the FIT is equal to or greater than that of the ELRT. Although the FIT does not account for the sampling process from a population explicitly, the test was shown to work well even when the sampling process exists if the sample size is not small.

The FIT is a simple and bias-controlled method to detect selection. However, it is not clear whether the FIT works well when the population size fluctuates. Theoretically, under the null model with fluctuating population size *Y _{i}* does not follow the same normal distribution for all

*i*, and therefore,

*t*(Data) does not follow the Student’s

_{FIT}*t*distribution.

This study is an extension of Feder *et al.*’s (2013) FIT that allows for fluctuations in population size by using reference loci. First, the FIT’s actual type I error rates in a fluctuating population is investigated. Then, a new test is introduced, the frequency increment test with reference loci (FITR). Given a fluctuating population size, the FITR’s actual type I error rates are almost the same as the nominal significance level. Then, the powers of the FITR and the FIT to detect natural selection were evaluated. Finally, the simulation method used in this study was validated and the properties of the FITR in practical situations were investigated. Model descriptions are presented just below and added before introducing the FITR.

## Materials and Methods

### Model and simulation methods

Let us consider a population evolving according to the Wright–Fisher model with fluctuating population size. The population size fluctuates as a function of generation time, *t*, and is denoted by *N*(*t*). To investigate the actual type I errors and the powers of Feder *et al.*’s (2013) FIT and the FITR introduced in this study, we conducted computer simulations under the five demographic models shown in Figure 1. The two alleles at a diallelic locus of interest are denoted by *A*_{0} and *a*_{0}, respectively. At generation times , the frequencies of *a*_{0} are denoted by . Here, *t*_{0} = 0 and *t _{L}* =

*T*are the first and the last sampling times, respectively, and the number of sampled times is

*L*+ 1. The fitnesses of genotypes

*A*

_{0}

*A*

_{0},

*A*

_{0}

*a*

_{0}, and

*a*

_{0}

*a*

_{0}are assumed to be 1, , and , respectively (

*i.e.*, no dominance is assumed). The population size,

*N*(

*t*), is independent of the frequency of

*a*

_{0}. As described in the next section, the FITR also uses neutral reference loci.

In the Wright–Fisher model, the allele frequency can be obtained exactly every generation as a binomial distribution. However, the generation of these data poses an extreme computational burden that is impractical for large populations (Figure 1). To avoid this burden and simulate changes in allele frequencies, we applied the pseudo-sampling method (Kimura and Takahata 1983), which is an improved version of the methods of Kimura (1980). In this method, to determine the allele frequency every generation, a uniform random number with the same mean and variance as those of the exact binomial distribution is used when the allele frequency is moderate. When the allele frequency is high or low, a Poisson random number with the same mean as that of the exact binomial random number is used. In this study, a frequency of ≤5 minor alleles in the population was used as the criteria for high or low allele frequencies. In addition, the normal distribution was used instead of the uniform distribution because the normal distribution better approximates the binomial distribution, which next-generation allele frequencies follow under the Wright–Fisher model.

## Results and Discussion

### Type I error rate of FIT

Table 1 summarizes the actual type I error rates for Feder *et al.*’s (2013) FIT. In demographic model 1 (constant-size model), as shown by Feder *et al.* (2013), the actual type I error rates approach the nominal level. For model 2 (slow-growth model) and model 3 (moderate-bottleneck model), the test becomes somewhat conservative. For the purposes of controlling type I error, this is a desirable property. However, in model 4 (rapid-growth model) and model 5 (severe-bottleneck model), the test tends to be too conservative, causing it to have less power to detect selection when the population size is fluctuating.

### Frequency increment test with reference loci (FITR)

Here we propose a new test, the FITR. Consider *R* reference loci in addition to the focal locus. It is assumed that these loci are evolving independently and that *R* reference loci are evolving under neutrality. We denote by the population frequencies of one allele at the *h*-th reference diallelic locus at times . Recall that *x*_{0,i} is the allele frequency of the focal locus at *t _{i}*. Suppose that

*N*(

*t*) is a step function and let

*N*be the population size from to such that . Note that although the following FITR discussion assumes is a step function, the same discussion can apply even the case that

_{i}*N*(

*t*) is a continuous function (Figure 1) because

*N*can be interpreted as the variance effective size over the period from . For this reason, the FITR is unbiased irrespective of fluctuations in population size.

_{i}When is not close to 0 or 1 and () is small compared with *N _{i}*, (2)wherefollows the standard normal distribution for

*h*= 0 under the null model and for under the null or alternative models. We then consider whether the allele frequency change from for the focal locus,

*Y*

_{0,i}, is significant. If

*N*is known, we can test for neutrality using the fact that

_{i}*Y*

_{0,i}follows the standard normal distribution. In this case, however,

*N*is unknown. Let us then define a statistic, (3) (4)The statistic is independent of

_{i}*N*, as seen in (4), because

_{i}*N*in (3) is canceled out. In addition, is independent of Δ

_{i}*t*. In (3), the numerator follows the standard normal distribution, and the denominator is equal to the square root of the χ

_{i}^{2}random variable divided by its degrees of freedom,

*R*. Because the numerator and denominator are independent, follows a Student’s

*t*distribution with

*R*degrees of freedom (Fisher 1925). Although we determined the form of the test statistic, , intuitively, can be derived as the exact LRS using data, , in a plausible setting (see Appendix). In other words, the aforementioned

*t*-test is equivalent to the likelihood ratio test under conditions of normality, as observed in several statistical situations (see,

*e.g.*, Lehmann and Romano 2005).

Next, let us define a statistic using all the data from to , , (5) (6) is the standardized sum of overall *i*. The standardization factor allows for a straightforward interpretation of the statistic because asymptotically follows the standard normal distribution as *R* becomes large. The exact distribution of with infinite *R* is difficult to express explicitly, but the distribution of can be obtained empirically by generating *L* Student’s *t* random variables with *R* degrees of freedom and summing them. This approach is valid because each follows a Student’s *t* distribution with *R* degrees of freedom. We obtained using the statistical package R (http://www.R-project.org) with 100,000 simulations of for each combination of *R* and *L*. The test using is the FITR, an exact significance test assuming *Y _{h,i}* follows the standard normal distribution. That is, the actual type I error rate of the test is expected to be close to the nominal significance level regardless of fluctuations in population size. Unlike , the statistic is not the exact LRS, which is very difficult to express explicitly. An ad hoc interpretation of the test statistic,

*t*, is presented in the Appendix.

_{FITR}### Type I error rate of FITR and powers of FITR and FIT

Table 2 shows the actual type I error rates of the FITR. As expected, for all demographic models, the actual type I error rates are close to the nominal level. Figure 2 shows the powers of the FITR and the FIT as a function of the strength of selection. In all demographic models, including the constant-size model, the FITR had more power than the FIT if five or more reference loci were used. For model 2 (slow-growth model) and model 3 (moderate-bottleneck model), the power of the FIT was acceptable. However, for model 4 (rapid-growth model) and model 5 (severe-bottleneck model), the power of the FIT was relatively small, and the FITR demonstrated much larger power than the FIT.

Figure 3A displays the powers of the FITR and the FIT as functions of the number of sampling points, *L*. In many cases, the powers approached certain asymptotic values with increasing *L*. In this case, the values of *L* for which the powers approached their asymptotes were relatively small (*e.g.*, *L* = 10 or 20). The powers of the FITR for *R* = 1 or 2 in Model 1 decreased somewhat as the sampling points increased. This trend also was observed for the FITR with *R* = 1 and for the FIT in Model 5. The powers of the FITR and the FIT as functions of the duration of sampling time, *T*, are given in Figure 3B. For all cases, the powers increased with increasing *T*, as expected. For Figure 3, A and B, the FITR had more power than the FIT if 10 or more reference loci were used even in the constant-size model. This difference in power was more obvious for model 5 (severe-bottleneck model).

### Applying the simulation method and FITR to practical situations

The FITR was developed and evaluated for type I error rate and power under ideal conditions of “selectively neutral” reference loci evolving independently of the focal locus and of each other, allele frequencies at reference loci ≠ 0 or 1, definable FITR statistics, and exactly known population allele frequencies. In practice, these ideal conditions may be violated. Before this section, the computer simulation method used in this study had not been validated. Here, we discuss the applicability of the simulation and describe cases that violate the aforementioned assumptions.

For the determination of allele frequencies in successive generations at *R* + 1 loci, the exact binomial sampling is computationally intensive and impractical for realistic population sizes (Figure 1). Therefore, we used a pseudosampling method to simulate the binomial sampling process (Wright−Fisher model). Even for a small, constant-size (*n* = 20) Wright−Fisher population, the fixation times for a mutant obtained by the pseudosampling were consistent with those obtained by binomial sampling (Kimura 1980). However, we considered only relatively short time scales, and the performance of the pseudo-sampling method was not obvious.

In Table 3, the rows denoted by *r* = “free” correspond with an assumption of free recombination (*i.e.*, evolving independently) among *R* + 1 loci. Rejection rates simulated by the exact binomial sampling and by the pseudosampling are given. Demographic Models 1 and 5 with reduced population sizes were used (see Table 3, legend). We did not observe any differences in the results generated by the binomial sampling *vs.* the pseudosampling under neutral or selective cases. These findings support the applicability of pseudo-sampling to our problem of concern.

Next, we considered the case in which the reference loci and focal loci were not independent. We limited our analysis to the case in which *R* + 1 loci were in linkage equilibrium (LE) at t=0. For closely linked loci, linkage disequilibrium (LD) is a distinct possibility. In addition, selection at the focal locus can drastically promote LD (*e.g.*, Sabeti *et al.* 2002). However, for example, empirical studies of human genomes suggest that LD can be extended, at most, to several megabase pairs from the selective locus (*e.g.*, Saunders *et al.* 2005). Because the genome is large compared with the megabase pairs scale, we can select *R* reference loci such that *R* + 1 loci are in LE. For this reason, our discussion is limited to the case in which *R* + 1 loci are in LE.

In Table 3, the rows indicated by *r* = 0.1, 0.01, and 0 describe results corresponding to a case in which the per-generation recombination fraction between any two adjacent loci are 0.1, 0.01, and 0, respectively. At *t* = 0, *R* + 1 loci are assumed to be in LE. That is, the alleles at *R* + 1 loci are randomly combined to form haplotypes. The simulations were conducted by the exact binomial sampling. For the neutral or selective cases, we observed no obvious differences between free recombination and limited recombination (*r* = 0.1, 0.01, and 0; Table 3). That is, the type I error rates and powers were maintained regardless of recombination fractions.

A case in which the reference loci are under selection is evaluated in Figure 4. The selection model is the same because the focal loci and *R* loci are under the same degree of selection. That is, for all loci, the fitnesses of genotypes *A _{h}A_{h}*,

*A*, and

_{h}a_{h}*a*are assumed to be 1, 1 + 0.5

_{h}a_{h}*s*, and 1 +

_{h}*s*(

_{h}*s*

_{1}=

*s*

_{2}= ··· =

*s*), respectively. The effects of selection at the reference loci are conservative for type I error rates (see the case of

_{R}*s*

_{0}= 0 in Figure 4). The results of Model 1 suggest that if , there is little difference in rejection rates compared to the neutral case. Including Model 5, if the condition

*s*≤ 1/2s

_{h}_{0}is met, the power is not decreased. That is, the power is not highly sensitive to selection at the reference loci. Nevertheless, we recommend using synonymous sites or noncoding regions as references.

We next considered a situation in which allele frequencies at the reference loci were low at *t* = 0 and some allele frequencies could become 0 or 1 by *t* = *T*. When allele frequencies became 0 or 1, the statistic, *t*_{FITR}, in (6) could not be defined. Therefore, it was practical to remove these loci from the calculation of *t*_{FITR}. Table 4 indicates how rejection rates are changed by the data handling. Values in parentheses indicate the average number of reference loci used to calculate the FITR statistics. For *L* = 2 the type I error rate was inflated by a few percent (*e.g.*, 6.69% at most, Model 1) possibly because changes in allele frequencies at reference loci are biased toward smaller values when loci are removed for which the frequencies of alleles become 0 or 1. These apparently reduced changes in allele frequencies could bring about overestimates of change at the selective locus. For *L* = 10 the inflation of the type I error rate becomes small. In general, to prevent inflation of the type I error rate, loci having moderate frequencies of alleles (*e.g.*, ≥10%) should be used in this test.

The effects of sampling error on the type I error rate and power of the FITR are shown in Figure 5. In general, the effects of sampling error on the type I error rate were conservative. As expected, the power decreased as the number of sampled individuals increased. The degree of power reduction differed for different demographic models or values of *L*. This finding reflects that the power is influenced by the relative magnitudes of changes in allele frequencies at *R* + 1 loci and sampling errors. As *L* increased, the relative changes in allele frequencies to the sampling errors decreased. Thus, power was more reduced for larger *L*. Regarding the demographic models, the population size of Model 5 was larger than that of Model 1. Therefore, the relative changes in allele frequencies to the sampling errors were larger in Model 5, and the degree of power is large in Model 5.

In this study, we proposed a neutrality test, the FITR, to accommodate fluctuations in population size using reference loci. Our test is an extension of Feder *et al.*’s (2013) FIT. By computer simulation, the actual type I error rate of the FITR was nearly equal to the nominal significant level regardless of fluctuations in population size when sampling noise could be ignored. The FITR detected selection with remarkable power under conditions of rapid growth (model 4) and severe bottleneck (model 5). Even under a model of constant population size, the FITR using 10 or more reference loci had more power than the FIT.

We also discussed the performance of the FITR in practical situations. The effects of selection at the reference loci were small unless selection was strong. Our findings indicated that when *R* + 1 were in LE, those loci should be considered independent of each other. In addition, loci with moderate frequencies of alleles should be used as references. Our findings may facilitate the development of more sophisticated methods using independent reference loci, including a method that can quantify (estimate) the strength of selection. These methods will enable appropriate inferences about natural selection in real and dynamic populations. Figure 5.

## Acknowledgments

I thank K. Ikeo and S. Mano for useful suggestions on the study. I also thank the associate editor and the two anonymous reviewers for their helpful comments on an earlier version of this manuscript, which improved this work substantially. This work was supported by health labor sciences research grant from The Ministry of Health Labor and Welfare (H23-jituyouka(nanbyou)-006).

## Appendix

Suppose that there are two segregating alleles, denoted by *A _{h}* and

*a*, at the

_{h}*h*-th locus . The fitnesses of genotypes

*A*,

_{h}A_{h}*A*, and

_{h}a_{h}*a*are assumed to be , respectively. As described in the main text, the case of

_{h}a_{h}*h*= 0 corresponds to the focal locus. When is small compared with

*N*

_{i}, the change in the allele frequency of , , from approximately follows a normal distribution,The normalized change, , is approximatelyThen the probability density of , , under is

Considering for all , the joint probability density of , , under .

(A1)i)

**is the exact LRS using data**We will show that given by (3) or (4) is the exact statistic using data to test (for all

*h*) against and in the likelihood ratio test framework.In our model, the likelihood ratio test is a test for which we reject if (A2)Otherwise, we accept

*H*_{0}. Here, is the maximum likelihood estimator (MLE) of*N*under_{i}*H*_{0}, and and are the MLEs of and under*H*_{1}.Under

*H*_{0}, the probability density of is given byfrom (A1). Solvingwe getThenUnder

*H*_{1}, the probability density of is given byfrom (A1). Solvingwe getThenThus, is equivalent toWith some algebra, we let a new constant Finally, we getwhere . Therefore,*t*as given by (3) or (4) is the exact LRS using data ._{FITR(i)}ii)

**Ad hoc interpretation for***t*as a test statistic using the data Δ_{FITR}*x*As

*R*grows, the*t*distribution for*t*approaches the normal distribution with mean and variance 1,Unfortunately, varies with_{FITR(i)}*N*_{i}, , and . When ; 0.50, 0.490, 0.458, 0.400, 0.300 and 0.218, respectively. Nevertheless, the values of are roughly the same if are far from 0 or 1. In addition, the values of and vary with*i*. However, we try to consider as a constant*μ*, such thatConsider a test*H*_{0}:μ = 0 against*H*_{1}:μ ≠ 0 using*L*i.i.d samples from the distribution. The LRS is given by the sum of , . That is, when*R*is large and the variation of the value for is not large,*t*given by (5) or (6) is expected to be close to the LRS._{FITR}

## Footnotes

*Communicating editor: P. Pfaffelhuber*

- Received August 19, 2013.
- Accepted September 26, 2013.

- Copyright © 2013 Nishino

This is an open-access article distributed under the terms of the Creative Commons Attribution Unported License (http://creativecommons.org/licenses/by/3.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.