Modifying the Chi-square and the CMH test for population genetic inference: adapting to over-dispersion

TL;DR

This paper introduces modified chi-square and CMH tests that account for over-dispersion in population genetic data, reducing false positives in genome-wide studies involving allele frequency changes.

Contribution

The authors develop simple adjusted test statistics that incorporate over-dispersion, improving inference accuracy in evolving populations with noise from genetic drift and sequencing.

Findings

Adjusted tests reduce false positive rates.

Application to Drosophila data demonstrates effectiveness.

Formulas applicable to other over-dispersed data scenarios.

Abstract

Evolve and resequence studies provide a popular approach to simulate evolution in the lab and explore its genetic basis. In this context, the chi-square test, Fishers exact test, as well as the Cochran-Mantel-Haenszel test are commonly used to infer genomic positions affected by selection from temporal changes in allele frequency. However, the null model associated with these tests does not match the null hypothesis of actual interest. Indeed due to genetic drift and possibly other additional noise components such as pool sequencing, the null variance in the data can be substantially larger than accounted forby these common test statistics. This leads to p-values that are systematically too small and therefore a huge number of false positive results. Even, if the ranking rather than the actual p-values is of interest, a naive application of the mentioned tests will give misleading results, as the amount of over-dispersion varies from locus to locus. We therefore propose adjusted statistics that take the over-dispersion into account while keeping the formulas simple. This is particularly useful in genome-wide applications, where millions of SNPs can be handled with little computational effort. We then apply the adapted test statistics to real data fromDrosophila, and investigate how in-formation from intermediate generations can be included when avail-able. The obtained formulas may also be useful in other situations, provided that the null variance either is known or can be estimated.

Tables

Table 1. Table 1: Standard contingency table used with chi-square test. Subsequent interpretation in our population genetic application: Entries are allele frequencies for a bi-allelic SNP taken either from two populations or from one population at two time points. n 𝑛 n is the total number of sequencing reads (coverage), x i ​ j subscript 𝑥 𝑖 𝑗 x_{ij} are the reads for allele j 𝑗 j in population i 𝑖 i , x i + subscript 𝑥 limit-from 𝑖 x_{i+} is the total number of reads in population i 𝑖 i , x + j subscript 𝑥 𝑗 x_{+j} is the total number of allele j 𝑗 j in both populations, i , j ∈ { 1 , 2 } 𝑖 𝑗 1 2 i,j\in\{1,2\} . The frequencies are obtained either by individual sequencing of a sample or by pool sequencing applied to the entire population.

	allele 1	allele 2
population 1	$x_{11}$	$x_{12}$	$x_{1+}$
population 2	$x_{21}$	$x_{22}$	$x_{2+}$
	$x_{+1}$	$x_{+2}$	$n$

Table 2. Table 2: Allele frequencies for a bi-allelic SNP taken either from two populations or from one population at two time points, assuming two underlying sampling steps. m 𝑚 m is the total number of sequencing reads (coverage) for the considered SNP, x ^ i ​ j subscript ^ 𝑥 𝑖 𝑗 \hat{x}_{ij} are the reads for allele j 𝑗 j in population i 𝑖 i , r i subscript 𝑟 𝑖 r_{i} is the total number of reads in population i 𝑖 i , x ^ + j subscript ^ 𝑥 𝑗 \hat{x}_{+j} is the total number of allele j 𝑗 j in both populations, i , j ∈ { 1 , 2 } 𝑖 𝑗 1 2 i,j\in\{1,2\} .

	allele 1	allele 2
population 1	${\hat{x}}_{11}$	${\hat{x}}_{12}$	$r_1$
population 2	${\hat{x}}_{21}$	${\hat{x}}_{22}$	$r_2$
	${\hat{x}}_{+1}$	${\hat{x}}_{+2}$	$m$

Table 3. Table 3: Estimators s ^ 1 2 superscript subscript ^ 𝑠 1 2 \hat{s}_{1}^{2} and s ^ 2 2 superscript subscript ^ 𝑠 2 2 \hat{s}_{2}^{2} of Var ⁡ ( x 11 ) Var subscript 𝑥 11 \operatorname{Var}(x_{11}) and Var ⁡ ( x 21 ) Var subscript 𝑥 21 \operatorname{Var}(x_{21}) for different scenarios.

	${\widehat{𝒔}}_{\mathrm{𝟏}}^{\mathrm{𝟐}}$	${\widehat{𝒔}}_{\mathrm{𝟐}}^{\mathrm{𝟐}}$
1 sampling step⋆	$x_{1+}\frac{x_{+1}}{n}\frac{x_{+2}}{n}$	$x_{2+}\frac{x_{+1}}{n}\frac{x_{+2}}{n}$
1 sampling step, drift	$\frac{x_{11}{x}_{12}}{x_{1+}}$	$x_{2+}\left({\hat{p}}_2\left(1-{\hat{p}}_2\right)+\left(x_{2+}-1\right){\hat{\sigma }}_{\text{drift}}^2\right)$
2 sampling steps	$\frac{{\hat{x}}_{11}{\hat{x}}_{12}}{r_1}\left(1+\frac{r_1-1}{x_{1+}}\right)$	$\frac{{\hat{x}}_{21}{\hat{x}}_{22}}{r_2}\left(1+\frac{r_2-1}{x_{2+}}\right)$
2 sampling steps, drift	$\frac{{\hat{x}}_{11}{\hat{x}}_{12}}{r_1}\left(1+\frac{r_1-1}{x_{1+}}\right)$	$r_2\left({\hat{p}}_2\left(1-{\hat{p}}_2\right)\left(1+\frac{r_2-1}{x_{2+}}\right)+\left(r_2-1\right)\frac{x_{2+}-1}{x_{2+}}{\hat{\sigma }}_{\text{drift}}^2\right)$
⋆ This is the situation of the classical chi-square test.

Table 4. Table 4: Estimators s ^ 1 ​ k 2 superscript subscript ^ 𝑠 1 𝑘 2 \hat{s}_{1k}^{2} and s ^ 2 ​ k 2 superscript subscript ^ 𝑠 2 𝑘 2 \hat{s}_{2k}^{2} of Var ⁡ ( x 11 ​ k ) Var subscript 𝑥 11 𝑘 \operatorname{Var}(x_{11k}) and Var ⁡ ( x 21 ​ k ) Var subscript 𝑥 21 𝑘 \operatorname{Var}(x_{21k}) for different scenarios.

	${\widehat{𝒔}}_{\mathrm{𝟏}𝒌}^{\mathrm{𝟐}}$	${\widehat{𝒔}}_{\mathrm{𝟐}𝒌}^{\mathrm{𝟐}}$
1 sampling step⋆	$x_{1+k}\frac{x_{+1k}}{n_k}\frac{x_{+2k}}{n_k-1}$	$x_{2+k}\frac{x_{+1k}}{n_k}\frac{x_{+2k}}{n_k-1}$
1 sampling step, drift	$\frac{x_{11k}{x}_{12k}}{x_{1+k}}$	$x_{2+k}\left({\hat{p}}_{2k}\left(1-{\hat{p}}_{2k}\right)+\left(x_{2+k}-1\right){\hat{\sigma }}_{\text{drift-k}}^2\right)$
2 sampling steps	$\frac{{\hat{x}}_{11k}{\hat{x}}_{12k}}{r_{1k}}\left(1+\frac{r_{1k}-1}{x_{1+k}}\right)$	$\frac{{\hat{x}}_{21k}{\hat{x}}_{22k}}{r_{2k}}\left(1+\frac{r_{2k}-1}{x_{2+k}}\right)$
2 sampling steps, drift	$\frac{{\hat{x}}_{11k}{\hat{x}}_{12k}}{r_{1k}}\left(1+\frac{r_{1k}-1}{x_{1+k}}\right)$	$r_{2k}\left({\hat{p}}_{2k}\left(1-{\hat{p}}_{2k}\right)\left(1+\frac{r_{2k}-1}{x_{2+k}}\right)+\left(r_{2k}-1\right)\frac{x_{2+k}-1}{x_{2+k}}{\hat{\sigma }}_{\text{drift-k}}^2\right)$
⋆ This is the situation of the classical CMH test.

Table 5. Table 5: Comparison of type I error, power, and running times in seconds for different tests of selection. Simulation setup: 10 4 superscript 10 4 10^{4} loci with true allele frequencies in base population uniformly distributed on [ 0 , 1 ] 0 1 [0,1] , 10% of the loci under selection, selection coefficients exponentially distributed with mean 0.1, N e = 300 subscript 𝑁 𝑒 300 N_{e}=300 , allele sample size 1000, pool sequencing with coverage Poisson distributed with mean 80, sequence data for generations 0 and 60, additionally for generations 10, 20, 30, 40, 50 when T χ 2 2s&d-ig superscript subscript 𝑇 superscript 𝜒 2 2s&d-ig T_{\chi^{2}}^{\text{2s\&d-ig}} and T CMH 2s&d-ig superscript subscript 𝑇 CMH 2s&d-ig T_{\text{CMH}}^{\text{2s\&d-ig}} are applied, 5 replicate populations in case of T CMH 2s&d superscript subscript 𝑇 CMH 2s&d T_{\text{CMH}}^{\text{2s\&d}} and T CMH 2s&d-ig superscript subscript 𝑇 CMH 2s&d-ig T_{\text{CMH}}^{\text{2s\&d-ig}} .

Test	Type I Error	Power	Time
	0.050	0.417	0.005
	0.046	0.490	0.369
with 5 replicates	0.050	0.761	0.014
with 5 replicates	0.049	0.756	0.487
$T_{{\chi }^2}$ with empirical FDR by (Orozco-terWengel et al., 2012)	0.374	0.226	0.144
LLS⋆ (Taus, Futschik and Schlötterer, 2017)	0.047	0.456	33568.010

Equations

$$T_{\chi^{2}}:=\sum_{i=1}^{2}\sum_{j=1}^{2}\frac{\big{(}x_{ij}-\frac{x_{i+}\,x_{+j}}{n}\big{)}^{2}}{\frac{x_{i+}\,x_{+j}}{n}}=\frac{n(x_{11}x_{22}-x_{12}x_{21})^{2}}{x_{1+}\;x_{2+}\;x_{+1}\;x_{+2}}.$$

$$T_{\chi^{2}}^{a}(\hat{s}_{1}^{2},\hat{s}_{2}^{2}):=\frac{(x_{11}x_{22}-x_{12}x_{21})^{2}}{x_{2+}^{2}\hat{s}_{1}^{2}+x_{1+}^{2}\hat{s}_{2}^{2}}=\frac{\big{(}x_{11}-\frac{x_{1+}\,x_{+1}}{n}\big{)}^{2}}{\big{(}\frac{x_{2+}}{n}\big{)}^{2}\hat{s}_{1}^{2}+\big{(}\frac{x_{1+}}{n}\big{)}^{2}\hat{s}_{2}^{2}},$$

$$T_{\text{CMH}}:=\frac{\big{(}\sum_{k=1}^{k^{*}}(x_{11k}-\frac{x_{1+k}\,x_{+1k}}{n_{k}})\big{)}^{2}}{\sum_{k=1}^{k^{*}}\frac{x_{1+k}\,x_{+1k}\,x_{2+k}\,x_{+2k}}{n_{k}^{2}(n_{k}-1)}},$$

$$T_{\text{CMH}}^{a}\big{(}\hat{s}_{1k}^{2},\hat{s}_{2k}^{2};k=1,\ldots,{k^{*}}\big{)}:=\frac{\big{(}\sum_{k=1}^{k^{*}}(x_{11k}-\frac{x_{1+k}\,x_{+1k}}{n_{k}})\big{)}^{2}}{\sum_{k=1}^{k^{*}}\big{(}(\frac{x_{2+k}}{n_{k}})^{2}\hat{s}_{1k}^{2}+(\frac{x_{1+k}}{n_{k}})^{2}\hat{s}_{2k}^{2}\big{)}}.$$

$$p_{1}\,(1-p_{1})\,\Big{(}1-\big{(}1-\frac{1}{2\,N_{e}}\big{)}^{t}\Big{)},$$

$$\frac{\frac{x_{11}}{x_{1+}}+\frac{x_{21}}{x_{2+}}}{2}\text{ or }\frac{\frac{{\hat{x}}_{11}}{r_1}+\frac{{\hat{x}}_{21}}{r_2}}{2}\text{ respectively}.$$

$$\frac{x_{11}\,x_{12}}{x_{1+}^{2}}\,\Big{(}1-\big{(}1-\frac{1}{2\,N_{e}}\big{)}^{t}\Big{)}\quad\text{ or }\quad\frac{\hat{x}_{11}\,\hat{x}_{12}}{r_{1}^{2}}\,\Big{(}1-\big{(}1-\frac{1}{2\,N_{e}}\big{)}^{t}\Big{)}\text{ respectively}.$$

$$\frac{{\sum }_{i=1}^{\gamma }{\hat{f}}_i}{\gamma }.$$

$$\sum_{i=1}^{\gamma-1}\hat{f}_{i}\,(1-\hat{f}_{i})\,\Big{(}1-\big{(}1-\frac{1}{2\,N_{e}}\big{)}^{t_{i+1}-t_{i}}\Big{)}.$$

$$X_{11}\sim \mathrm{Binom}(x_{1+},p_1)\text{ and }{X}_{21}\sim \mathrm{Binom}(x_{2+},p_2).$$

$$T_n:=(X_{11}-\frac{x_{1+}{X}_{+1}}{n}).$$

$$\begin{split}\operatorname{Var}(T_{n})&=\operatorname{Var}\Big{(}\frac{n}{n}X_{11}-\frac{x_{1+}}{n}(X_{11}+X_{21})\Big{)}\\ &=\operatorname{Var}\Big{(}\frac{n-x_{1+}}{n}X_{11}-\frac{x_{1+}}{n}X_{21}\Big{)}\\ &=\big{(}\frac{x_{2+}}{n}\big{)}^{2}\operatorname{Var}(X_{11})+\big{(}\frac{x_{1+}}{n}\big{)}^{2}\operatorname{Var}(X_{21}).\end{split}$$

$$\big{(}\frac{x_{2+}}{n}\big{)}^{2}\hat{s}_{1}^{2}+\big{(}\frac{x_{1+}}{n}\big{)}^{2}\hat{s}_{2}^{2}\stackrel{{\scriptstyle\mathbb{P}}}{{\longrightarrow}}\operatorname{Var}(T_{n}).$$

$$\begin{split}\mathbb{E}[T_{n}]&=\mathbb{E}\Big{[}\frac{x_{2+}}{n}X_{11}-\frac{x_{1+}}{n}X_{21}\Big{]}\\ &=\frac{x_{2+}}{n}\mathbb{E}[X_{11}]-\frac{x_{1+}}{n}\mathbb{E}[X_{21}]\\ &=\frac{x_{2+}}{n}\,x_{1+}\,p-\frac{x_{1+}}{n}\,x_{2+}\,p\\ &=0\end{split}$$

$$\frac{T_n}{\sqrt{\mathrm{Var}(T_n)}}\stackrel{d}{⟶}\mathcal{N}(0,1),$$

$$\frac{T_n^2}{\mathrm{Var}(T_n)}\stackrel{d}{⟶}{\chi }_1^2,$$

$$\frac{T_{n}^{2}}{\big{(}\frac{x_{2+}}{n}\big{)}^{2}\hat{s}_{1}^{2}+\big{(}\frac{x_{1+}}{n}\big{)}^{2}\hat{s}_{2}^{2}}=\frac{\big{(}X_{11}-\frac{x_{1+}\,X_{+1}}{n}\big{)}^{2}}{\big{(}\frac{x_{2+}}{n}\big{)}^{2}\hat{s}_{1}^{2}+\big{(}\frac{x_{1+}}{n}\big{)}^{2}\hat{s}_{2}^{2}}=:T_{\chi^{2}}^{a}(\hat{s}_{1}^{2},\hat{s}_{2}^{2})\stackrel{{\scriptstyle d}}{{\longrightarrow}}\chi^{2}_{1}.$$

$$\mathrm{Var}(X_{11})=x_{1+}p(1-p)\text{ and }\mathrm{Var}(X_{21})=x_{2+}p(1-p).$$

$$\begin{split}T_{\chi^{2}}^{a}(x_{1+}\frac{X_{+1}}{n}\frac{X_{+2}}{n},x_{2+}\frac{X_{+1}}{n}\frac{X_{+2}}{n})&=\frac{(X_{11}X_{22}-X_{12}X_{21})^{2}}{x_{2+}^{2}x_{1+}\frac{X_{+1}}{n}\frac{X_{+2}}{n}+x_{1+}^{2}x_{2+}\frac{X_{+1}}{n}\frac{X_{+2}}{n}}\\ &=\frac{(X_{11}X_{22}-X_{12}X_{21})^{2}}{x_{2+}x_{1+}X_{+1}X_{+2}\,\frac{1}{n^{2}}\,(x_{2+}+x_{1+})}\\ &=\frac{\big{(}X_{11}\,X_{22}-X_{12}\,X_{21}\big{)}^{2}}{x_{1+}\;x_{2+}\;X_{+1}\;X_{+2}\,\frac{1}{n}}\\ &=T_{\chi^{2}}.\end{split}$$

$$\big{(}\frac{x_{2+k}}{n_{k}}\big{)}^{2}\hat{s}_{1k}^{2}+\big{(}\frac{x_{1+k}}{n_{k}}\big{)}^{2}\hat{s}_{2k}^{2}\stackrel{{\scriptstyle\mathbb{P}}}{{\longrightarrow}}\operatorname{Var}(T_{nk}),$$

$$\operatorname{Var}\Big{(}\sum_{k=1}^{k^{*}}T_{nk}\Big{)}=\sum_{k=1}^{k^{*}}\operatorname{Var}(T_{nk}).$$

$$T_{\text{CMH}}^{a}\big{(}\hat{s}_{1k}^{2},\hat{s}_{2k}^{2};k=1,\ldots,{k^{*}}\big{)}\stackrel{{\scriptstyle d}}{{\longrightarrow}}\chi^{2}_{1},$$

$$X_{11}\sim \mathrm{Binom}(x_{1+},p_1)\text{and}{X}_{21}\sim \mathrm{Binom}(x_{2+},P_2),$$

$$\mathrm{Var}(X_{11})=x_{1+}{p}_1(1-p_1)$$

$$x_{1+}\frac{X_{11}}{x_{1+}}\frac{X_{12}}{x_{1+}}=\frac{X_{11}{X}_{12}}{x_{1+}}\stackrel{\mathbb{P}}{⟶}\mathrm{Var}(X_{11}).$$

$$\begin{split}\operatorname{Var}(X_{21})&=\mathbb{E}[\operatorname{Var}(X_{21}\,|\,P_{2})]+\operatorname{Var}(\mathbb{E}[X_{21}\,|\,P_{2}])\\ &=\mathbb{E}[x_{2+}\,P_{2}\,(1-P_{2})]+\operatorname{Var}(x_{2+}\,P_{2})\\ &=x_{2+}\,\mathbb{E}[P_{2}]-x_{2+}\,\mathbb{E}[P_{2}^{2}]+x_{2+}^{2}\operatorname{Var}(P_{2})\\ &=x_{2+}\,\mathbb{E}[P_{2}]-x_{2+}\,\big{(}\operatorname{Var}(P_{2})+\mathbb{E}[P_{2}]^{2}\big{)}+x_{2+}^{2}\operatorname{Var}(P_{2})\\ &=x_{2+}\Big{(}\mathbb{E}[P_{2}]\big{(}1-\mathbb{E}[P_{2}])+\big{(}x_{2+}-1\big{)}\operatorname{Var}(P_{2})\Big{)}\end{split}$$

$$x_{2+}\Big{(}\hat{p}_{2}\big{(}1-\hat{p}_{2})+\big{(}x_{2+}-1\big{)}\hat{\sigma}_{\text{drift}}^{2}\Big{)}\stackrel{{\scriptstyle\mathbb{P}}}{{\longrightarrow}}\operatorname{Var}(X_{21}).$$

$$K_1\sim \mathrm{Binom}(g_1,b_0)\text{and}{K}_2\sim \mathrm{Binom}(g_2,B_1)\text{with }{B}_1=\frac{K_1}{g_1}.$$

$$\mathrm{Var}(K_2)$$

$$\displaystyle=\frac{g_{2}}{g_{1}}\,\Big{(}\mathbb{E}[K_{1}]\Big{(}1-\frac{\mathbb{E}[K_{1}]}{g_{1}}\Big{)}+\frac{g_{2}-1}{g_{1}}\operatorname{Var}(K_{1})\Big{)}$$

$$\displaystyle=g_{2}\,b_{0}\,(1-b_{0})\Big{(}1+\frac{g_{2}-1}{g_{1}}\Big{)}.$$

Full text

chi-square test,

CMH test,

over-dispersion,

experimental evolution,

evolve and resequence,

genetic drift,

pool sequencing,

,

and

t1previously: Kerstin Gärtner

1 Introduction

An important question in the field of population genetics is how populations adapt to changes in their environment. Experimental evolution allows to study the adaptation under controlled conditions. If these experiments are combined with high-throughput sequencing, they are called evolve and resequencing (E&R) experiments (Turner et al., 2011). Such experiments are carried out both on microbes and on higher organisms. Due to large population sizes and short generation times, microbes permit to study evolutionary processes based on newly arriving mutations. With higher organisms and sexual reproduction on the other hand, evolution based on standing genetic variation is usually explored. A major goal is to identify genomic positions (a.k.a. loci) that are responsible for the adaptation. For this purpose, the organisms are kept over $t+1$ generations $G_{0},G_{1},\ldots,G_{t}$ under conditions that require adaptation. Allele frequencies are obtained by sequencing the genomes at $G_{0}$ and $G_{t}$ and possibly also at some intermediate time points. Depending on the (time and financial) budget, the organisms are sequenced individually, or as a pool in order to obtain allele frequencies for typically millions of single-nucleotide polymorphisms(SNPs). Individual sequencing can also be implemented using barcoding, where barcode tags that identify the organism are added before sequencing.

Allele frequency changes over time are then tested for signals of selection. For this purpose, usually only bi-allelic SNPs are considered. Indeed, multi-allelic sites are rare in population data and likely caused by sequencing errors (Burke et al., 2010a). Consequently frequencies of the two alleles are used in the base and the evolved population for each tested SNP. Pearson’s chi-square test (for simplicity subsequently called chi-square test) is a very popular method for this purpose (Griffin et al., 2017). Serving the same purpose, Fishers exact test is sometimes used as an alternative (Burke et al., 2010b). Being a generalization of the chi-square test for stratified data, the CMH test is also often used when allele frequency data from several replicate populations are available. See e.g. Orozco-terWengel et al. (2012); Tobler et al. (2014); Nouhaud et al. (2016); Barghi et al. (2017) for applications.

Kofler and Schlötterer (2014) compare several methods for detecting selection by contrasting allele frequencies in different generations. Apart from the CMH test, they consider the pairwise summary statistic “diffStat” (Turner et al., 2011), an association statistic by Turner and Miller (2012), and $F_{ST}$ (Remolina et al., 2012). A comparison of receiver operator curves (ROC) for these tests shows that the CMH test performs best, i.e. has more power than the other tests considered to identify selected SNPs. The aforementioned methods consider two time points, i.e. two generations, only. However, in E&R experiments sometimes organisms are sequenced also at intermediate generations, resulting in repeated measurement data. Further methods are available for detecting selection in this context. The method of Bollback, York and Nielsen (2008), for instance, is based on a hidden Markov model (HMM). Generalizations of this approach are provided by Malaspinas et al. (2012) and Steinrücken, Bhaskar and Song (2014). Also, Mathieson and McVean (2013) adapt HMMs to structured populations. The method CLEAR by Iranmehr et al. (2017) uses Markov chains in a discrete state model and computes the exact likelihood for small populations. Ignoring spatial dependence, linked loci are modeled using composite likelihood statistics. A frequency increment test (FIT) based on an approximation of the allele frequency dynamics by a Gaussian process is proposed by Feder, Kryazhimskiy and Plotkin (2014). Considering all loci separately, Topa et al. (2015) also model the allele frequency trajectories by Gaussian processes, whereas Terhorst, Schlötterer and Song (2015) approximate the joint likelihood for multiple loci. The approach by Taus, Futschik and Schlötterer (2017) is based on linear least square (LLS) regression to fit the allele frequency data to a selection model. Schraiber, Evans and Slatkin (2016) estimate parameters in a Bayesian framework with Markov-chain Monte-Carlo sampling. Another Bayesian approach has been proposed by Levy et al. (2015) for estimating parameters in barcoded lineages. Finally, the Wright-Fisher ABC method proposed by Foll, Shim and Jensen (2015) applies approximate Bayesian computation (ABC).

Besides detecting loci under selection, several of the discussed approaches additionally estimate selection coefficients, often jointly with other parameters like the effective population size (e.g. Bollback, York and Nielsen, 2008), age of alleles (e.g. Malaspinas et al., 2012; Schraiber, Evans and Slatkin, 2016), or allelic dominance (e.g. Taus, Futschik and Schlötterer, 2017). However, such methods are computationally much more demanding than simple methods like the chi-square and the CMH tests. Hence, the latter are still widely used when testing for selection and e.g. implemented in the software tool PoPoolation2 (Kofler, Pandey and Schlötterer, 2011).

When comparing the allele frequencies between pairs of samples, the null model for the classical chi-square, Fishers exact test, and also the CMH test assumes that the probability of sampling a given allele is the same within each given pair. However, the sampling variation is not the only component of variance relevant in E&R experiments. Allele frequency changes between generations happen because of genetic drift, i.e. due to chance. This increases the variance in the data noticeably unless population sizes are large enough to safely ignore drift. Such a situation usually occurs only when working with micro-organisms (e.g. Illingworth et al., 2012). Another potential source of random variation is pool sequencing where the obtained reads can be regarded as a sample from the DNA pool.

When the chi-square or the CMH tests are applied to data, which contain more variance than assumed by the tests (over-dispersion), the resulting values of the test statistic are too large and hence the p-values too small. As a consequence, selection is often inferred for loci where it is not present. In the simulations with drift and pool sequencing, that we carried out, the null hypothesis of neutrality is rejected in up to 80 % of the cases, despite being true. Hence, the additional variance introduced by drift and pool sequencing is by no means negligible. A common procedure to account for it is calculating an empirical false discovery rate (FDR) for false positives due to drift and pool sequencing via computer simulations (Orozco-terWengel et al., 2012). Griffin et al. (2017) chose another approach by applying three different statistical tests and considering only SNPs which are significant with respect to all three methods as candidate loci.

This issue is also well known in the unrelated context of complex surveys, where different strategies have been proposed to obtain appropriate tests of homogeneity. (See e.g. chapter 10 in Lohr (2009).)

Our aim is to adapt the chi-square test and the CMH test in a way that additional sources of variance are directly included in the test statistics, making computer simulations no longer necessary. Thus, we propose a method that is faster than the commonly used ones. When sequence data for intermediate generations is available, the additional information can be included into our test statistics without a considerable increase in computation time. Also in terms of power, our method performs better than other approaches. In particular, our method has considerably more power to detect selected SNPs than the classical CMH test with the empirical FDR correction (Orozco-terWengel et al., 2012). Compared to the approach for detecting selection in Taus, Futschik and Schlötterer (2017), which the authors describe faster than the CLEAR method by Iranmehr et al. (2017) and the Wright-Fisher ABC by Foll, Shim and Jensen (2015), our method has also slightly more power, and is $10^{5}$ orders of magnitude faster. Hence, our method outperforms other common approaches to detect selected loci in speed and power.

In this article, we first present variants of the chi-square and the CMH tests for general underlying variances. The statistics derived in this step can be useful in all situations, where over-dispersion is present.

Further, we provide specific formulas for the test statistics under scenarios with drift and pool sequencing, which are common in E&R experiments, and seen also in other situations. In genome-wide association studies (GWAS) e.g. the CMH test is often used for the inference of an association between a trait and an allele variant, and when the data arise from pool sequencing (e.g. Bastide et al., 2013; Endler et al., 2016), our adapted test could be a good alternative.

The remainder of this article is structured as follows: The test statistics for general underlying variances are presented in section 2 and the scenarios with drift and pool sequencing are considered in section 2.1. In section 3 we apply the adapted tests to simulated data and examine their performance. We apply the tests to real data and present the results in section 4. A discussion and an outlook in section 5 conclude this article.

2 Adapted Tests

In this section, we generalize the chi-square and the CMH tests to work under over-dispersion and derive explicit formulas for scenarios with drift and pool sequencing. Given the application in mind, our focus is on the null model of homogeneity, although our derivations also apply to the test of independence.

We summarize our data in a 2x2 contingency table. Using the notation from table 1, the chi-square test statistic in its standard form is defined as

[TABLE]

As shown in appendix A1 (equation (A10)), we may rewrite $T_{\chi^{2}}$ as follows

[TABLE]

with $\hat{s}_{1}^{2}:=x_{1+}\frac{x_{+1}}{n}\frac{x_{+2}}{n}$, and $\hat{s}_{2}^{2}:=x_{2+}\frac{x_{+1}}{n}\frac{x_{+2}}{n}$.

In order to adapt the test to models that involve different variances, we may now replace $\hat{s}_{1}^{2}$ and $\hat{s}_{2}^{2}$ by consistent estimators of $\operatorname{Var}(x_{11})$ and $\operatorname{Var}(x_{21})$. As for the classical chi-square test (see appendix A1, in particular equation (A8)), $T_{\chi^{2}}^{a}(\hat{s}_{1}^{2},\hat{s}_{2}^{2})$ converges in distribution to a $\chi^{2}$-distribution with 1 degree of freedom under the null hypothesis of homogeneity.

The CMH test is based on a 2x2x${k^{*}}$ contingency table, where the ${k^{*}}$ partial 2x2 tables are assumed to be independent. We use the same notation as for the chi-square test and indicate a variable belonging to the $k^{th}$ partial table with the additional index $k$. The null hypothesis is that both true proportions within each partial table are the same, i.e. the odds ratio in each partial table equals 1 (McDonald, 2014). The classical CMH test statistic is

[TABLE]

see Agresti (2002) 111 This is the test statistic as proposed by Mantel and Haenszel, which is commonly considered for the CMH test. The statistic proposed by Cochran differs by the factor $\frac{1}{n_{k}}$ instead of $\frac{1}{n_{k}-1}$ in each term of the denominator. Asymptotically this difference is negligible. .

Analogous to the chi-square test, one can adapt the CMH test to general underlying variances. As a first step, we write the test statistic of the CMH test as

[TABLE]

and insert $x_{1+k}\frac{x_{+1k}}{n_{k}}\frac{x_{+2k}}{n_{k}-1}$ for $\hat{s}_{1k}^{2}$ and $x_{2+k}\frac{x_{+1k}}{n_{k}}\frac{x_{+2k}}{n_{k}-1}$ for $\hat{s}_{2k}^{2}$, $k=1,\ldots,{k^{*}}$. (See appendix A2.) As with the chi-square test, the formula assumes one sampling step only, which is not appropriate for more complex models. Again, however, $\hat{s}_{1k}^{2}$ and $\hat{s}_{2k}^{2}$ can be replaced by consistent estimators of $\operatorname{Var}(x_{11k})$ and $\operatorname{Var}(x_{21k})$, $k=1,\ldots,{k^{*}}$. In the next section we present suitable variance estimators for situations with drift and pool sequencing.

2.1 Adaptation of the tests to drift and pool sequencing

We first focus on the chi-square test and assume that allele frequency data is available for a single population at two time points. Under the null hypothesis, the population allele frequency at the later time point $p_{2}$ for a given bi-allelic SNP arises from $p_{1}$ according to the Wright-Fisher model of genetic drift (see e.g. chap. 3 in Ewens, 2004). Therefore $p_{2}$ is modeled as a random variable. Usually one cannot observe $p_{1}$ and $p_{2}$ directly. Indeed both quantities are frequently estimated from population samples. If experimenters use pool sequencing as a further sampling step, we model this by binomial sampling as e.g. in (Waples, 1989) or (Jónás et al., 2016). The random sample size is known as coverage, and the success probability is taken as the frequency of allele 1 in the underlying DNA material. If only a sample of the population is sequenced, we assume again binomial sampling for simplicity: An extension to hypergeometric sampling is straightforward.

Typically genomic selection is taken as alternative hypothesis that leads to differences between our estimates for $p_{1}$ and $p_{2}$ that cannot be explained by sampling and drift.

Figure 1(a) summarizes the scenario with drift and one underlying sampling step, which is either taking a sample from the population for (individual) sequencing, or applying pool sequencing (to the whole population). The scenario with drift and two sampling steps (sampling from the population and pool sequencing) is outlined in figure 1(b).

Table 1 summarizes the notation for a contingency table based on one sampling step (sampling from the population or pool sequencing), while table 2 is for two sampling steps (sampling plus pool sequencing). Note that in the scenario with two sampling steps, only the sample sizes $x_{1+}$ and $x_{2+}$ but not the allele frequencies are known for the first step. In the situations where one population descends from another, population 1 is the base population and population 2 is the evolved population.

Conditional on $p_{1}$, the variance in allele frequency due to drift after $t$ generations can be calculated as

[TABLE]

where $p_{1}$ is the frequency of allele 1 in the base population and $2\,N_{e}$ is the effective population size considering gametes (Falconer, 1960).

If we have data from ${k^{*}}$ replicate populations, the above considerations hold analogously for every replicate $k\in\{1,\ldots,{k^{*}}\}$ and we indicate the $k^{\text{th}}$ replicate by an additional index $k$ such as in $p_{ik}$, $i\in\{1,2\}$.

We derive variance estimators for the described scenarios in Appendix A3. In table 3, we present our estimators for $\operatorname{Var}(x_{11})$ and $\operatorname{Var}(x_{21})$, and in table 4 for $\operatorname{Var}(x_{11k})$ and $\operatorname{Var}(x_{21k})$. To obtain the adapted test statistics for the different scenarios, they can be inserted for $\hat{s}_{1}^{2}$ and $\hat{s}_{2}^{2}$ in (2), or $\hat{s}_{1k}^{2}$ and $\hat{s}_{2k}^{2}$ in (4) respectively. The proposed formulas use estimators $\hat{p}_{2}$ and $\hat{p}_{2k}$ of $\mathbb{E}[p_{2}|p_{1}]$ and $\mathbb{E}[p_{2k}|p_{1k}]$. Also, $\hat{\sigma}_{\text{drift}}^{2}$ and $\hat{\sigma}_{\text{drift-k}}^{2}$ are estimators of $\operatorname{Var}(p_{2}|p_{1})$ and $\operatorname{Var}(p_{2k}|p_{1k})$. Choices for these quantities are discussed below.

Notice that different models may apply at different time points. If an experiment involves for instance individual sequencing of a sample from the base population and pool sequencing of a sample of the evolved population, the variance estimators should be chosen accordingly: One would take (2) as test statistic with $\hat{s}_{1}^{2}$ replaced by $\frac{x_{11}\,x_{12}}{x_{1+}}$ and $\hat{s}_{2}^{2}$ replaced by r_{2}\big{(}\hat{p}_{2}(1-\hat{p}_{2})\big{(}1+\frac{r_{2}-1}{x_{2+}}\big{)}+(r_{2}-1)\,\frac{x_{2+}-1}{x_{2+}}\hat{\sigma}_{\text{drift}}^{2}\big{)}.

A simple estimator for $\mathbb{E}[p_{2}|p_{1}]$ is $\frac{x_{11}}{x_{1+}}$ or $\frac{\hat{x}_{11}}{r_{1}}$, depending on the number of underlying sampling steps. However, our simulations show that often the distribution of the corresponding p-values is closer to a uniform distribution on $[0,1]$ when $\mathbb{E}[p_{2}|p_{1}]$ is estimated as

[TABLE]

For a consistent variance estimator of $p_{2}$ after $t$ generations of drift, $\hat{\sigma}_{\text{drift}}^{2}$, we can approximate (5) by

[TABLE]

When sequence data for intermediate generations between 0 and $t$ are available, we can use this additional information for the estimation of $\mathbb{E}[p_{2}|p_{1}]$ and $\operatorname{Var}(p_{2}|p_{1})$.

Let $t_{1}=0,t_{2},\ldots,t_{\gamma}\,=\,t$ be the generations for which sequence data is available, and let $p_{1}=f_{1},f_{2},\ldots,f_{\gamma}=p_{2}$ be the corresponding population frequencies of allele 1. Estimating these frequencies in each generation by the corresponding relative sample frequencies $\hat{f}_{1},\ldots,\hat{f}_{\gamma}$, we may proceed as in (6) and use this additional information to estimate $\mathbb{E}[p_{2}|p_{1}]$ by

[TABLE]

Extending (7), the drift variance may also be estimated as

[TABLE]

Analogous estimators may be used for the estimation of $\mathbb{E}[p_{2k}|p_{1k}]$ and $\operatorname{Var}(p_{2k}|p_{1k})$ in the situation with replicate populations.

We explore the behavior of the adapted test statistics by computer simulations and present the results in the following section. Hereby we focus on the scenario with two sampling steps and genetic drift, since there the additional variance is the largest.

We introduce the following notation: If we don’t have data for intermediate generations, we estimate $\mathbb{E}[p_{2}|p_{1}]$ and $\operatorname{Var}(p_{2}|p_{1})$, as well as $\mathbb{E}[p_{2k}|p_{1k}]$ and $\operatorname{Var}(p_{2k}|p_{1k})$ by the estimators given in (6) and (7) or their analogs for the $k^{\text{th}}$ of ${k^{*}}$ replicates. With drift and one sampling step we denote the adapted tests then $T_{\chi^{2}}^{\text{1s\&d}}$ and $T_{\text{CMH}}^{\text{1s\&d}}$, with drift and two sampling steps we name them $T_{\chi^{2}}^{\text{2s\&d}}$ and $T_{\text{CMH}}^{\text{2s\&d}}$. If data for intermediate generations is available, we apply the estimators (8) and (9) or their analogs for the $k^{\text{th}}$ replicate. In the case of drift, we denote the adapted tests by $T_{\chi^{2}}^{\text{1s\&d-ig}}$ resp. $T_{\text{CMH}}^{\text{1s\&d-ig}}$ (one sampling step), and $T_{\chi^{2}}^{\text{2s\&d-ig}}$ resp. $T_{\text{CMH}}^{\text{2s\&d-ig}}$ (two sampling steps).

3 Simulation Results

We carried out extensive simulations in R (R-Core-Team, 2018) in order to explore the behavior of the adapted tests described in the previous section. We simulated genetic drift using the package poolSeq (Taus, Futschik and Schlötterer, 2017). When we encountered loci with frequency 0 for one allele in the base population but a positive frequency in a later generation, we changed the allele count from 0 to 1 in order to always obtain a well-defined test statistic. When these situations occur with real data, there are different possible explanations: Either a mutation arose in a later generation and the allele really was not present before, or the frequency of the respective allele is low, but not 0, in the base population and the allele was just by chance not sampled or amplified in the sequencing process. Since mutation rates are usually low over such a time span (Burke et al., 2010a), the latter scenario is the more likely one. Finally, if the frequency in the later generation is very low, the nonzero frequency may also be due to a sequencing error. Overall, our method to deal with this phenomenon seems to be a pragmatic compromise.

We first provide results under the null hypothesis, then we examine the power of our adapted test statistics. After that we compare our adapted tests to other state-of-the-art methods. We first simulate allele frequencies in generation 0 uniformly distributed on $[0,1]$ to give all possible true allele frequencies the same weight. At the end of the section we also consider an allele frequency distribution in generation 0 that resembles the one encountered in our experimental data.

Null distribution

In our simulations we set $N_{e}=300$, and used 1000 as sample size of alleles that were sequenced at generations 0 and 60; the sequencing coverage was chosen Poisson distributed with mean 80. These parameter choices were motivated by the real data for Drosophila taken from (Barghi et al., 2017) and discussed in section 4.

To control the type I error it is desirable that the p-values of a test are uniformly distributed on $[0,1]$ or at least stochastically larger than uniform if the null hypothesis is true. Indeed, under the neutral Wright–Fisher model we observe that the distribution of the p-values belonging to the adapted tests is close to a uniform distribution and that the tests control the 5 % significance level. In contrast, the non-adapted tests show a huge excess of small p-values. Based on $10^{6}$ simulated loci, figure 2 displays the distribution of the p-values for the classical CMH test with 3 replicates and the test adapted to drift, sampling, and pool sequencing $T_{\text{CMH}}^{\text{2s\&d}}$. For the adapted chi-square test $T_{\chi^{2}}^{\text{2s\&d}}$, the distribution of the p-values is slightly further from a uniform distribution but controls the 5% significance threshold (appendix figure A3). We obtain similar results in situations with drift and only one sampling step (appendix, figure A4). With two underlying sampling steps but without drift, a situation occurring e.g. in GWAS using pool sequencing data (e.g. Bastide et al., 2013; Endler et al., 2016), the adapted tests show again an improved performance (appendix, figure A5).

In genome-wide applications and corresponding corrections for multiple testing, also significance levels much below 0.05 become important. To check the type I error control in such a situation, we did a simulation analysis for $10^{6}$ loci based on sequence data with the parameter choices as above (but 5 replicate populations in case of the CMH test). As shown in figure 3, especially $T_{\chi^{2}}^{\text{2s\&d}}$ turns out to be anti-conservative for very small significance levels. The red lines correspond to the full, unfiltered data sets.

To understand this issue, it should be noted that the distributional approximations involved in these tests are less reliable for loci with a very small or very high allele frequency.

An obvious remedy would be to introduce a threshold value $\gamma$ and only consider loci with frequency of allele 1 as well as allele 2 larger than $\gamma$ in the base population. According to figure 3, filtering has the desired effect, if $\gamma$ is chosen large enough. Indeed, the asymptotic $\chi^{2}$-distribution of the test statistic under the null model for the chi-square test as well as for the CMH test is based on the normal approximation of the binomial distribution, which is more accurate for intermediate frequencies for finite samples.

One disadvantage of this approach is that by filtering out SNPs with small and large allele frequencies, we exclude a lot of potentially selected loci. We therefore explore also other approaches to resolve the issue: If sequence data is available not only from two time points but also from intermediate generations, we can modify the adapted tests by taking the additional information into account, resulting in the test statistics $T_{\chi^{2}}^{\text{1s\&d-ig}}$, $T_{\chi^{2}}^{\text{2s\&d-ig}}$, $T_{\text{CMH}}^{\text{1s\&d-ig}}$, and $T_{\text{CMH}}^{\text{2s\&d-ig}}$, see section 2.1. Simulations with the same parameters as in figure 3, but additionally with sequence data every 10 generations, show that $T_{\chi^{2}}^{\text{2s\&d-ig}}$ and $T_{\text{CMH}}^{\text{2s\&d-ig}}$ hold the 5% level or are conservative for all significance levels without filtering out small and large allele frequencies (appendix, figure A6).

If time series data is not available, a p-value correction may be applied by fitting the distribution $F_{p}$ of p-values simulated under the null hypothesis and transforming the p-values to uniform using $F_{p}(\cdot)$. In appendix A4 we propose a parametric choice of $F_{p}$ for correcting small (potentially significant) p-values, and show that it leads to a substantial improvement of the null distribution of p-values obtained with $T^{2s\&d}_{\chi^{2}}.$

Power

We additionally carried out simulations involving $10^{5}$ selected loci in order to examine the power of the adapted tests at a significance threshold of $\alpha=0.05$. We first consider a realistic set of standard parameter choices: $N_{e}=300$, sample size: 1000, coverage: 100, and sequence data available for generations 0 and 60. To also investigate the influence of these model parameters, we considered a set of alternative values for each of them.

Not surprisingly, the power of $T_{\text{CMH}}^{\text{2s\&d}}$ is higher than the power of $T_{\chi^{2}}^{\text{2s\&d}}$ because the amount of information increases with more replicates. As also expected, the power increases with the selection coefficient $s$ and the effective population size $N_{e}$. Figure 4(a) shows the power of $T_{\chi^{2}}^{\text{2s\&d}}$ for different values

of $s$ and $N_{e}$. Figure 4(b) displays the same for $T_{\text{CMH}}^{\text{2s\&d}}$ with 5 replicates. Comparing figures 4(c) and 4(d), which are for the same scenario as in figure 4(a) but with 20 and 200 generations of evolution, we see that for small selection coefficients the power increases with the number of generations. This is since more generations of evolution lead to larger frequency differences between base and evolved populations unless the selection coefficients are large. Under strong selection many alleles soon reach frequency close to 1 (fixation), and hence the signal of selection does not become stronger anymore with more generations. On the other hand, the drift variance that we calculate for the denominator of the test statistic increases with every generation, which reduces the power. The gain in power due to more generations of evolution and the loss in power due to fixation may cancel each other out, leading to the plateauing effect at around $0.8$ we observe in figure 4(d). The impact of the number of generations on the power is qualitatively the same for $T_{\text{CMH}}^{\text{2s\&d}}$. (appendix, figure A7).

The influence of sample size and coverage is shown in figure 5, in (a) for $T_{\chi^{2}}^{\text{2s\&d}}$, and in (b) for $T_{\text{CMH}}^{\text{2s\&d}}$. The power increases with the sample size and with the coverage. Since the coverage values are an order of magnitude smaller than the values for the sample size, the effect is much more pronounced for the coverage.

The power of $T_{\chi^{2}}^{\text{2s\&d-ig}}$ and $T_{\text{CMH}}^{\text{2s\&d-ig}}$ is similar to the power of $T_{\chi^{2}}^{\text{2s\&d}}$ and $T_{\text{CMH}}^{\text{2s\&d}}$, see table 5. In general, $T_{\chi^{2}}^{\text{2s\&d-ig}}$ and $T_{\text{CMH}}^{\text{2s\&d-ig}}$ are affected in the same way by the different parameters as $T_{\chi^{2}}^{\text{2s\&d}}$ and $T_{\text{CMH}}^{\text{2s\&d}}$ (not shown).

Method comparison

We compared the performance of our adapted test statistics to other state-of-the-art methods. As summaries we considered the type I error, the power, and the run time. We looked at the classical chi-square test with an empirical FDR correction (Orozco-terWengel et al., 2012) and the LLS approach for detecting selection of Taus, Futschik and Schlötterer (2017). As stated in Taus, Futschik and Schlötterer (2017), the LLS method is faster than the methods CLEAR by Iranmehr et al. (2017) and Wright-Fisher ABC by Foll, Shim and Jensen (2015). Still the computation times with LLS method are high and we restricted our simulations therefore here to $10^{4}$ loci of which 10 % are under selection. The other parameters were chosen again $N_{e}=300$, allele sample size 1000, coverage Poisson distributed with mean 80, 60 generations of drift.

Compared with the LLS method by Taus, Futschik and Schlötterer (2017), $T_{\chi^{2}}^{\text{2s\&d}}$ and $T_{\chi^{2}}^{\text{2s\&d-ig}}$ have a similar power to detect selection while being much faster, see table 5. The chi-square test with the empirical FDR correction (Orozco-terWengel et al., 2012) is not much slower than our tests (although the empirical FDR correction requires additional computations). As the test statistic does not use the correct variance terms, it performs worse, however, than our corrected tests both in terms of the power and the type I error.

Alternative starting allele frequency distribution

In the previous simulations the true allele frequencies in the base populations were chosen uniformly distributed on $[0,1]$. In general, uniformly distributed allele frequencies are not common in natural populations. Therefore we also looked for a more realistic distribution of allele frequencies. As an example we consider the u-shaped beta distribution proposed in (Jónás et al., 2016) for the fruit fly Drosophila. The u-shape depicts the observed excess of high and low allele frequencies in a folded site frequency spectrum where it is not known which allele is ancestral and which one is derived.

We therefore chose the true allele frequencies in generation 0 from a beta-distribution with parameters 0.2 and 0.2. The other parameters remained unchanged compared to the previous simulations. As the chosen distribution produces a higher proportion of allele frequencies close to the boundaries, the distribution of the p-values becomes less uniform (see figure 6). The spikes are caused by discreteness phenomena occurring for very low starting allele frequencies. However, $T_{\chi^{2}}^{\text{2s\&d-ig}}$ and $T_{\text{CMH}}^{\text{2s\&d-ig}}$ turned out to be quite conservative (even for significance levels smaller than 0.05) when sequence data every 10 generations is available. At the same time power values of the adapted tests were lower than when under uniformly distributed starting allele frequencies (appendix, figure A8). An explanation would be that even SNPs with positively selected alleles have a large probability of being lost due to drift in early generations, if the initial allele frequency is very low. For such SNPs we do not have any power of detection.

In particular, we calculated a power of approximately 0.59 for $T_{\text{CMH}}^{\text{2s\&d}}$ based on $10^{5}$ simulations with selection coefficient 0.1, allele sample size 1000, pool sequencing in generations 0 and 60 with coverage Poisson distributed with mean 80, and 5 replicates. When additional sequence data for generations 10, 20, 30, 40, and 50 was simulated, we obtained a power of approximately 0.60 using $T_{\text{CMH}}^{\text{2s\&d-ig}}$ as test statistic.

4 Application to experimental data from Drosophila

Here we consider data from an evolve and resequence experiment on Drosophila simulans as described in (Barghi et al., 2017). In this publication the classical CMH test has been used to infer candidates of selection. Allele frequency measurements were taken from three replicate populations at generations 0 and 60. All flies were maintained under a cycling routine of 12 hours at 18∘C in a dark environment to mimic night and 12 hours at 28∘C with light for the day.

In the original paper neutral simulations have been used to define a cut-off that leads to 2% false positive SNPs under the simulated global null model. This cut-off has then been taken as a threshold for the p-values obtained with the CMH test applied to the real data, see (Orozco-terWengel et al., 2012) for a more detailed description. Notice, however, that with approximately 4 million SNPs (Barghi et al., 2017), this approach will lead to approximately 80000 false positive SNPs and an unclear false discovery rate.

An advantage of our approach is that the resulting proper p-values can be combined with a standard procedure such as Bonferroni or Benjamini-Hochberg that controls for multiple testing. Another advantage of our method is that it will lead to a more proper ranking of the p-values, as the amount of error incurred with the classical test statistics depends on the relative magnitudes of the variance components (drift variance, sequencing coverage and sample size) and will therefore vary between SNPs.

Since the DNA from the whole population was used in the pool sequencing step, we applied $T_{\chi^{2}}^{\text{1s\&d-ig}}$ and $T_{\text{CMH}}^{\text{1s\&d-ig}}$ as test statistics using the model parameters specified in (Barghi et al., 2017). We only present the results for the modified CMH test on the entire data set. A Manhattan plot of the SNP positions versus the logarithm (base 10) of the p-values corrected for multiple testing with the Benjamini-Hochberg method is shown in figure 7. The computation time needed for this analysis was only about 20 seconds on a standard laptop. We infer more significant SNPs than (Barghi et al., 2017) ($0.0049\%$ and $0.0002\%$ of the total number of tested SNPs, respectively). This result is concordant with our simulations showing that our method has higher power compared to the empirical FDR correction (Orozco-terWengel et al., 2012) which is also applied in (Barghi et al., 2017). The population genetics group at the University of Veterinary Medicine in Vienna plans to apply our tests on new data sets that will become available soon.

5 Discussion

With population genetic applications in mind, we propose modified test statistics for the chi-square and the CMH tests in scenarios with over-dispersion, i.e. more variance in the data than considered by the original tests. Compared with the classical versions of these tests that are still commonly applied, our approach does not require simulations to find a cut-off for the test statistic. Our proposed approach can also be used instead of Fishers exact test which faces the same problems in this context. Our inclusion of proper variance terms leads to a better performance both with respect to power and type I error. Using the classical tests does not even lead to a proper ranking of the SNPs as the amount of over-dispersion varies between the SNPs. While more sophisticated testing procedures have also been proposed, they are usually considerably more time consuming, especially when applied on a genome wide scale.

While this issue is also known in other applications (such as complex surveys), our underlying model requires a different adaptation of the test statistics. We therefore expressed the test statistics in dependence of the variances of entries in the contingency tables as a first step. Using this general form, suitable test statistics can be calculated in any situation with over-dispersion provided that the required variances can be properly estimated.

We then derived explicit formulas for the adapted test statistics for use in different types of E&R experiments under scenarios with one or two sampling steps, resulting from actual sampling or pool sequencing, and genetic drift.

Our test statistics do not only provide a more appropriate error control, but have also a much larger power than the the classical tests combined with computer simulations such as proposed in Orozco-terWengel et al. (2012).

Compared with more sophisticated methods, our adapted tests have approximately the same power but need much less computation time, which is an important factor when the whole genome is scanned for traces of selection. Our tests are e.g. $10^{5}$ times faster than the LLS method by Taus, Futschik and Schlötterer (2017).

When sequence data is available for two time points only, our tests are not always conservative for very small allele frequencies in the base population. This problem is related with inaccuracies of the normal approximation of proportions close to zero and one. In this case, we propose a correction of the p-values based on simulations under the null model. Notice, however, that this additional correction is usually not needed when data is available also at intermediate time points, as supported by our simulated scenarios.

We implemented the adapted tests in an R package called ACER, which can be downloaded on https://github.com/MartaPelizzola/ACER. If many loci or even the whole genome is tested for selection, we recommend to control for the FDR by the Benjamini-Hochberg or a similar procedure.

Our results suggest that our adapted test statistics provide fast, reliable and powerful methods to detect selection. Hence, they have the potential to considerably facilitate the inference of selected loci in population genetics, in particular in the context of E&R.

Acknowledgments

We thank Neda Barghi and Christian Schlötterer for providing the data and helpful comments, as well as Thomas Taus for his input on the acceleration of the R-code. This work has been supported by the Austrian Science Fund (FWF Doctoral Program “Vienna Graduate School of Population Genetics”, DK W1225-B20).

Appendix

Supplemental material

In our derivations we either assume $p_{1}$ to be not random, or do our computations conditional on $p_{1}$.

A1 Derivation of $T_{\chi^{2}}^{a}(\hat{s}_{1}^{2},\hat{s}_{2}^{2})$

We restrict our analysis here on 2x2 tables since we consider two samples222For time series data and hence more sequence samples, the samples from the intermediate time points are considered only for a variance estimator within the here described framework of two sequence samples. and bi-allelic SNPs. However, our results could be generalized to contingency tables with larger dimensions by expressing the test statistic in dependence of the respective variances and replacing these variances by appropriate estimators, which we will do now for 2x2 tables.

We follow the rationale of chapter 14.3 in (Lehmann and Romano, 2005) for the goodness-of-fit test, but adopt it to the tests of independence and homogeneity, and derive a test statistic $T_{\chi^{2}}^{a}(\hat{s}_{1}^{2},\hat{s}_{2}^{2})$ for the chi-square test depending on the consistent estimators $\hat{s}_{1}^{2}$ and $\hat{s}_{2}^{2}$ of $\operatorname{Var}(X_{11})$ and $\operatorname{Var}(X_{21})$.

We assume $X_{11}$ and $X_{21}$ to be independent. As it is usually done when testing for homogeneity, we consider the total number of alleles in each generation, $x_{1+}$ and $x_{2+}$, as known and assume

[TABLE]

Note that $p_{2}$ will often be a random variable, such as under genetic drift. Under the null hypothesis we have that $E(p_{2}|p_{1})=p_{1}=:p$. We assume that the null hypothesis is true.

We define an auxiliary variable

[TABLE]

With $X_{+1}=X_{11}+X_{21}$ we obtain

[TABLE]

Hence, with the consistent estimators $\hat{s}_{1}^{2}$ and $\hat{s}_{2}^{2}$ of $\operatorname{Var}(X_{11})$ and $\operatorname{Var}(X_{21})$, we obtain

[TABLE]

Above we considered $X_{+1}$ (and hence also $X_{+2}=n-X_{+1}$) as random variable. If, however, we condition on $x_{+1}$ (and $x_{+2}$), which is e.g. done for the test for independence, we are in a hypergeometric framework and can estimate the variances from there. The estimator for $\operatorname{Var}(T_{n})$ from the hypergeometric framework yields a factor $\frac{1}{n-1}$ where in the binomial framework with $X_{+1}$ random there is a factor $\frac{1}{n}$. Since we are interested in the asymptotic behavior of the test statistic, this difference is negligible.

With (A3) and (A1) it is further

[TABLE]

With the central limit theorem we can thus deduce

[TABLE]

and with the continuous mapping theorem we conclude

[TABLE]

where $\chi^{2}_{1}$ stands for the $\chi^{2}$-distribution with 1 d.f. Finally, we obtain with (A4), (A7), and again the continuous mapping theorem

[TABLE]

Now we show that setting $\hat{s}_{1}^{2}=x_{1+}\frac{X_{+1}}{n}\frac{X_{+2}}{n}$ and $\hat{s}_{2}^{2}=x_{2+}\frac{X_{+1}}{n}\frac{X_{+2}}{n}$ yields the original test statistic $T_{\chi^{2}}$.

From (A1) we know that

[TABLE]

Further, $\frac{X_{+1}}{n}$ is a consistent estimator for $p$, and hence $1-\frac{X_{+1}}{n}=\frac{X_{+2}}{n}$ is a consistent estimator for $(1-p)$. This means that $x_{1+}\frac{X_{+1}}{n}\frac{X_{+2}}{n}$ and $x_{2+}\frac{X_{+1}}{n}\frac{X_{+2}}{n}$ indeed are consistent estimators for $\operatorname{Var}(X_{11})$ and $\operatorname{Var}(X_{21})$.

[TABLE]

A2 Derivation of $T_{\text{CMH}}^{a}(\hat{s}_{1k}^{2},\hat{s}_{2k}^{2})$

In Agresti (2002) the test statistic is already presented in a form similar to (3).

The derivation can be carried out in a similar way as for the chi-square test. For each of the ${k^{*}}$ partial 2x2 tables we obtain analogously to the chi-square test, compare (A4),

[TABLE]

with $T_{nk}:=(X_{11k}-\frac{x_{1+k}\,X_{+1k}}{n_{k}})$. Due to the independence of the partial 2x2 tables it is

[TABLE]

Hence with the analogous considerations as carried out in section A1 we obtain

[TABLE]

with T_{\text{CMH}}^{a}\big{(}\hat{s}_{1k}^{2},\hat{s}_{2k}^{2};k=1,\ldots,{k^{*}}\big{)} as defined in (4).

A3 Derivation of the variance estimators for the different scenarios

For the different scenarios we derive in the following the estimators of $\operatorname{Var}(X_{11})$ and $\operatorname{Var}(X_{21})$ for the chi-square test. The respective estimators for the CMH test can be deduced analogously.

A3.1 Genetic drift and one sampling step

We consider a bi-allelic SNP and assume to have allele counts from individual sequencing for allele 1 and allele 2 in the base population and a population that has evolved for $t$ generations. Table 1 is the corresponding contingency table and the sampling scheme is summarized in figure 1(a).

We assume that the null hypothesis of homogeneity is true. Let be

[TABLE]

corresponding to (A1). Due to drift, the frequency of allele 1 after $t$ generations, $P_{2}$, is a random variable with $\mathbb{E}[P_{2}]=p_{1}$. This is in contrast to the original test, where $P_{2}=p_{1}$ under the null hypothesis. Note that still holds $E[T_{n}]=0$ as in equation (A5).

Recall that we assume $X_{11}$ and $X_{21}$ to be independent. $p_{1}$ can consistently be estimated as $\frac{X_{11}}{x_{1+}}$. Drift affects only the evolved population, so we can simply infer from the Binomial distribution

[TABLE]

and

[TABLE]

For the calculation of $\operatorname{Var}(X_{21})$ we use the law of total variance:

[TABLE]

With the consistent estimators $\hat{p}_{2}$ of $\mathbb{E}[P_{2}]$ and $\hat{\sigma}_{\text{drift}}^{2}$ of $\operatorname{Var}(P_{2})$ we obtain

[TABLE]

A3.2 No drift and two sampling steps

We now assume that all sequence data is from the same generation, and hence there is no drift, and we obtain the data from pool sequencing of a sample of the population. We model this as two binomial sampling steps. For a given SNP we have $r_{1}$ reads for the first and $r_{2}$ reads for the second population, i.e. $r_{1}$ and $r_{2}$ are the sequencing coverage values. Our model is to sample $r_{1}$ reads from a binomial distribution with probability $\frac{X_{11}}{x_{1+}}$ for the first population and $r_{2}$ reads with probability $\frac{X_{21}}{x_{2+}}$ for the second population. Table 2 is the corresponding contingency table.

Note that we only observe the data in table 2, while we don’t know the underlying data from table 1, except for the numbers of sequenced individuals $x_{1+}$ and $x_{2+}$.

Since we do two-step binomial sampling, we describe first a general two-step binomial sampling scenario. In (Jónás et al., 2016) similar calculations are carried out. Let be

[TABLE]

With the law of total variance we obtain

[TABLE]

$\hat{X}_{11}$ corresponds to $K_{2}$. In the first step we sample $x_{1+}$ alleles, and the probability for allele 1 is $p_{1}$. The second step is the sampling of the read counts as described above. With (A21) and estimating $p_{1}$ as $\frac{\hat{X}_{11}}{r_{1}}$ we can set

[TABLE]

Note that without drift, we have $X_{21}\sim\operatorname{Binom}(x_{2+},p_{2})$, where $p_{2}$ is not a random variable. For population 2 we have hence analogously

[TABLE]

A3.3 Genetic drift and two sampling steps

We now consider a scenario where additional variance is induced by drift as well as pool sequencing. We model drift and pool sequencing as previously described, the sampling scheme is summarized in figure 1(b). As drift does not affect the base population, we can estimate $\operatorname{Var}(\hat{X}_{11})$ again as in (A22).

For the estimator of $\operatorname{Var}(\hat{X}_{21})$ we use equation (A20). Since $P_{2}$ is a random variable like in (A14), the variance of $X_{21}$ is different from the binomial sampling variance and (A21) does not hold here. Similar calculations can be found in (Jónás et al., 2016).

Further, with the law of total expectation it is

[TABLE]

With (A20), (A17), and (A24) we obtain

[TABLE]

With the consistent estimators $\hat{p}_{2}$ of $\mathbb{E}[P_{2}]$ and $\hat{\sigma}_{\text{drift}}^{2}$ of $\operatorname{Var}(P_{2})$ we obtain then

[TABLE]

A4 P-value Correction

To obtain a more uniform tail distribution for the null p-values obtained with the chi-square test, we propose a correction that uses the following functions:

• •

$F_{A}(x)$ and $F_{A2}(x)$ are beta cumulative distribution functions with suitably chosen parameters

• •

$F_{B}(x):=\begin{cases}sF_{A2}(x/s)&0\leq x<s\\ \delta+\frac{x-\delta}{1-\delta}&s\leq x\leq 1\end{cases}$

where $\delta=sF_{A2}(1)$, and suitably chosen $s\in[0,z]$, $z$ in $[0,1]$.

For a p-value $\pi_{z}$, smaller than the chosen threshold $z$, we transform $\pi^{\ast}_{z}:=F_{t}(\pi_{z}):=F_{B}(z\cdot F_{A}(\pi_{z}/z))$. The parameters of the beta distributions have been obtained using method of moments estimates for the rescaled data, i.e. $\pi_{z}/z$ for $F_{A}(\cdot)$, and $\pi_{s}/s$ with $\pi_{s}=z\cdot F_{A}(\pi_{z}/z))$ for $F_{A2}(\cdot).$ Figure 1(a) displays the transformation $F_{t}(\cdot)$, and in figure 1(b) the improvement may be seen achieved by the transformation. Figure A2 suggests that some models and circumstances require different transformations.

A5 Additional figures