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

Abstract

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 χ²-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 χ² 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.

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_FIT (Data) does not follow the Student’s 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₀ and a₀, respectively. At generation times $t0=0,t1,…,tL=T$⁠, the frequencies of a₀ are denoted by $x0,0,x0,1, …,x0,L$⁠. Here, t₀ = 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₀A₀, A₀a₀, and a₀a₀ are assumed to be 1, $1+0.5s0$⁠, and $1+s0$⁠, respectively (i.e., no dominance is assumed). The population size, N(t), is independent of the frequency of a₀. As described in the next section, the FITR also uses neutral reference loci.

Figure 1

Open in new tabDownload slide

Demographic models used in this study. Model 1: constant-size model (N(t) = 10⁴); Model 2: slow-growth (grows exponentially from N(0) = 10⁴ to N(T) = 10⁵); Model 3: moderate-bottleneck model (⁠$N(0≤t<0.5T)=N(0.75≤t≤T)=5×104 and N(0.5≤t<0.75T)=104$⁠); Model 4: rapid-growth model (⁠$grows exponentially fromN(0)=104 to N(T)=108$⁠); Model 5: severe-bottleneck model (⁠$N(0≤t<0.5T)=N(0.75≤t≤T)=106 and N(0.5≤t<0.75T)=104$⁠).

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.

Actual type I error rates (%) of FIT

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 $R+1$ loci are evolving independently and that R reference loci are evolving under neutrality. We denote by $xh,0,xh,1,…,xh,L (h=1,2,…,R)$ the population frequencies of one allele at the h-th reference diallelic locus at times $t0=0,t1,…,tL=T$⁠. 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_i be the population size from $ti−1$ to $ti$ such that $N(ti−1<t≤ti)=Ni$⁠. Note that although the following FITR discussion assumes $N(t)$ is a step function, the same discussion can apply even the case that N(t) is a continuous function (Figure 1) because N_i can be interpreted as the variance effective size over the period from $ti−1 to ti$⁠. For this reason, the FITR is unbiased irrespective of fluctuations in population size.

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.

Actual type I error rates (%) of FITR

Figure 2

Open in new tabDownload slide

Powers of the FITR and the FIT in various demographic models. Powers of the FIT (black line) and the FITR with R reference loci (colored lines) are shown as functions of the selection coefficients in the five demographic models. Each point corresponds to the power obtained by 100,000 simulations at the 5% significance level. The duration of sampling time and the number of sampled points were T = 1000 and (L + 1) = 11, respectively. The intervals between any two adjacent sampled points were the same at $Δt=ti−ti−1=100 (i=1,2,…,L)$⁠. The initial frequency for all R+ 1 loci, $xh,0$⁠, was assumed to be 0.5.

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).

Figure 3

Open in new tabDownload slide

Powers of the FITR and the FIT as functions of (A) the number of sampling points and (B) the duration of sampling. Powers of the FIT (black line) and the FITR with R reference loci (colored lines) are shown for the demographic models 1 and 5. Each point corresponds to the power obtained by 100,000 simulations at the 5% significance level. The selection coefficients were s = 0.002 for model 1 and s = 0.0005 for model 5. The intervals between any two adjacent sampled points were the same, $Δt=ti−ti−1=T/L (i=1,2,…,L)$⁠. The initial frequency for all R+1 loci, $xh,0$⁠, was assumed to be 0.5. (A) The duration of sampling time was fixed at T = 1000. (B) The number of sampled points was fixed at L = 10. (Note: The FIT curve nearly overlapped with the FITR curve for L = 5 in Model 1.)

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.

Results obtained by binomial sampling with various recombination fractions

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 $h(≠0)$ loci, the fitnesses of genotypes A_hA_h, A_ha_h, and a_ha_h are assumed to be 1, 1 + 0.5s_h, and 1 + s_h (s₁ = s₂ = ··· = s_R), respectively. The effects of selection at the reference loci are conservative for type I error rates (see the case of s₀ = 0 in Figure 4). The results of Model 1 suggest that if $Nsh<5$⁠, there is little difference in rejection rates compared to the neutral case. Including Model 5, if the condition s_h ≤ 1/2s₀ 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.

Figure 4

Open in new tabDownload slide

The effects of selection at reference loci on the power of the FITR. The powers of the FITR under various selection strengths, s₀, at focal loci are shown as functions of the selection coefficient, $sh(h≠0)$⁠, at reference loci for demographic models 1 and 5. Each point corresponds to the power obtained by 100,000 simulations at the 5% significance level. The number of reference loci, the duration of sampling time, and the number of sampled points were R = 10, T = 1000, and (L + 1) = 11, respectively. The intervals between any two adjacent sampled points were the same at $Δt=ti−ti−1=100 (i=1,2,…,L)$⁠. The initial frequency for all $R+1$ loci, $xh,0$⁠, was assumed to be 0.5.

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.

Effects of allele frequencies at the reference loci

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.

Figure 5

Open in new tabDownload slide

The effects of sampling error on the power of the FITR. The powers of the FITR under various sampling regimes are shown as functions of the selection coefficients for demographic models 3 and 5. All, or 5000, 1000, 500, or 100 individuals in a population were assumed to be sampled. Sampling was assumed to be conducted by binomial sampling at each (R + 1) locus and at each (L + 1) time point. Each point corresponds to the power obtained by 100,000 simulations at the 5% significance level. The number of reference loci, the duration of sampling time, and the number of sampled points were R = 10, T = 1000, and (L + 1) = 3 (upper graphs) or 11 (lower graphs), respectively. The intervals between any two adjacent sampled points were the same at Δ_t = t_i–t_i-1 = 500 (top graphs) or 100 (bottom graphs) (i = 1,2,…,L). The initial frequency for all $R+1$ loci, $xh,0$⁠, were assumed to be 0.5.

Appendix

Considering $Δxh,i$ for all $h∈(0,1,2,…,R)$⁠, the joint probability density of $Δxi=(Δx0,i, Δx1,i⋯,ΔxR,i)$⁠, $fi(Δxi,s,Ni)$⁠, under $s=(s0,s1,⋯,sR) and Ni$⁠.

• i)

  $tFITR(i)$  is the exact LRS using data  $Δxi$

  We will show that $tFITR(i)$ given by (3) or (4) is the exact statistic using data $Δxi$ to test $H0:sh=0$ (for all h) against $H1:s0≠0$ and $sh=0$  $(h≠0)$ in the likelihood ratio test framework.

  In our model, the likelihood ratio test is a test for which we reject $H0$ if

  $λ=fi(Δxi,0,Nˇi)fi(Δxi,(s^0,0),N^i)<c (constant).$

  (A2)

  Otherwise, we accept H₀. Here, $Nˇi$ is the maximum likelihood estimator (MLE) of N_i under H₀, and $N^i$ and $s^0$ are the MLEs of $Ni$ and $s0$ under H₁.
  Under H₀, the probability density of $Δxh,i$ is given by

  $fi(Δxi,0,Ni)=(NiπΔti)R+12exp(−NiΔti∑h=0RΔxh,i2)$

  from (A1). Solving
  $∂lnfi∂Ni=R+12Ni−1Δti∑h=0RΔxh,i2=0,$

  we get

  $Nˇi=(R+1)Δti2∑h=0RΔxh,i2.$

  Then

  $fi(Δxi,0,Nˇi)=(NˇiπΔti)R+12exp(−R+12)$

  Under H₁, the probability density of $Δxh,i$ is given by

  $fi(Δxi,(s0,0),Ni)=(NiπΔti)R+12exp(−NiΔti(Δx0,i−Δti2s0x0,i−1(1−x0,i−1))2−NiΔti∑h=1RΔxh,i2)$

  from (A1). Solving

  $∂lnfi∂Ni=R+12Ni−1Δti(Δx0,i−Δti2s0x0,i−1(1−x0,i−1))2−1Δti∑h=1RΔxh,i2=0 and$

  $∂lnfi∂s0=Ni(Δx0,ix0,i−1(1−x0,i−1)−Δti2s0x0,i−1(1−x0,i−1))=0,$

  we get

  $N^i=(R+1)Δti2∑h=1RΔxh,i2 and$

  $s^0=2Δx0,iΔtix0,i−1(1−x0,i−1).$

  Then

  $fi(Δxi,(s^0,0),N^i)=(N^iπΔti)R+12exp(−R+12).$

  Thus, $(A2)$ is equivalent to

  $λ=fi(Δxi,0,Nˇi)fi(Δxi,(s^0,0),N^i)=(NˇiN^i)R+12=(1∑h=0RΔxh,i21∑h=1RΔxh,i2)R+12=(∑h=0RΔxh,i2∑h=1RΔxh,i2)−R+12=(1+Δx0,i2∑h=1RΔxh,i2)−R+12<c.$

  With some algebra, we let a new constant $c′=Rc−2R+1−R$

  $Δx0,i21R∑h=1RΔxh,i2>c′.$

  Finally, we get

  $Δx0,i1R∑h=1RΔxh,i2=tFITR(i)<−c″, c″<Δx0,i1R∑h=1RΔxh,i2=tFITR(i)$

  where $c″=c′$⁠. Therefore, t_FITR(i) as given by (3) or (4) is the exact LRS using data $Δxi$⁠.
• ii)

  Ad hoc interpretation for t_FITR as a test statistic using the data Δx

  As R grows, the t distribution for t_FITR(i) approaches the normal distribution with mean $s0NiΔtix0,i−1(1−x0,i−1)/2$ and variance 1,

  $tFITR(i)=Y0,i1R∑h=1RYh,i2∼N(s0NiΔtix0,i−1(1−x0,i−1)/2, 1).$

  Unfortunately, $s0NiΔtix0,i−1(1−x0,i−1)/2$ varies with N_i, $Δti$⁠, and $x0,i−1$⁠. When $x0,i−1=0.5, 0.4, 0.3,0.2,0.1, and 0.05$⁠; $x0,i−1(1−x0,i−1)=$ 0.50, 0.490, 0.458, 0.400, 0.300 and 0.218, respectively. Nevertheless, the values of $x0,i−1(1−x0,i−1)$ are roughly the same if $x0,i−1$ are far from 0 or 1. In addition, the values of $Ni$ and $Δti$ vary with i. However, we try to consider $s0NiΔtix0,i−1(1−x0,i−1)/2$ as a constant μ, such that

  $tFITR(i)∼N(μ,1).$

  Consider a test H₀:μ = 0 against H₁:μ ≠ 0 using L i.i.d samples from the distribution. The LRS is given by the sum of $tFITR(i)$⁠, $∑i=1LtFITR(i)$⁠. That is, when R is large and the variation of the value for $NiΔtix0,i−1(1−x0,i−1)$ is not large, t_FITR given by (5) or (6) is expected to be close to the LRS.

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).

Footnotes

Literature Cited