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

Kerstin<EMAIL_ADDRESS>[ Marta
<EMAIL_ADDRESS>[ Andreas Futschik
<EMAIL_ADDRESS>[ Vetmeduni Viennam1, Vienna Graduate School
of Population Geneticsm2 and Johannes Kepler University Linzm3

###### 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 for by 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 from Drosophila, and
investigate how information from intermediate generations can be included when
available. The obtained formulas may also be useful in other situations,
provided that the null variance either is known or can be estimated.

chi-square test,

CMH test,

over-dispersion,

experimental evolution,

evolve and resequence,

genetic drift,

pool sequencing,

###### keywords:

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

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$ is the total number of sequencing reads (coverage), $x_{ij}$ are the reads for allele $j$ in population $i$, $x_{i+}$ is the total number of reads in population $i$, $x_{+j}$ is the total number of allele $j$ in both populations, $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$

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

$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}}.$
(1)

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

$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}},$
(2)

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

$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)}},$
(3)

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

$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{)}}.$
(4)

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

(a) one sampling step

(b) two sampling steps

Figure 1: Sampling schemes for the scenarios with genetic drift and one
sampling step, shown in (a), or genetic drift and two sampling steps, shown in
(b).

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.

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$ is the total number of sequencing reads (coverage) for the considered SNP, $\hat{x}_{ij}$ are the reads for allele $j$ in population $i$, $r_{i}$ is the total number of reads in population $i$, $\hat{x}_{+j}$ is the total number of allele $j$ in both populations, $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$

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

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

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.

Table 3: Estimators $\hat{s}_{1}^{2}$ and $\hat{s}_{2}^{2}$ of $\operatorname{Var}(x_{11})$ and $\operatorname{Var}(x_{21})$ for different scenarios. |  $\hat{s}_{1}^{2}$ | $\hat{s}_{2}^{2}$
---|---|---
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+}\big{(}\hat{p}_{2}\big{(}1-\hat{p}_{2})+\big{(}x_{2+}-1\big{)}\hat{\sigma}_{\text{drift}}^{2}\big{)}$
2 sampling steps |  $\frac{\hat{x}_{11}\,\hat{x}_{12}}{r_{1}}\big{(}1+\frac{r_{1}-1}{x_{1+}}\big{)}$ | $\frac{\hat{x}_{21}\,\hat{x}_{22}}{r_{2}}\big{(}1+\frac{r_{2}-1}{x_{2+}}\big{)}$
2 sampling steps, drift |  $\frac{\hat{x}_{11}\,\hat{x}_{12}}{r_{1}}\big{(}1+\frac{r_{1}-1}{x_{1+}}\big{)}$ | $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{)}$
⋆ This is the situation of the classical chi-square test.
Table 4: Estimators $\hat{s}_{1k}^{2}$ and $\hat{s}_{2k}^{2}$ of $\operatorname{Var}(x_{11k})$ and $\operatorname{Var}(x_{21k})$ for different scenarios. |  $\hat{s}_{1k}^{2}$ | $\hat{s}_{2k}^{2}$
---|---|---
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}\big{(}\hat{p}_{2k}\big{(}1-\hat{p}_{2k})+\big{(}x_{2+k}-1\big{)}\hat{\sigma}_{\text{drift-k}}^{2}\big{)}$
2 sampling steps |  $\frac{\hat{x}_{11k}\,\hat{x}_{12k}}{r_{1k}}\big{(}1+\frac{r_{1k}-1}{x_{1+k}}\big{)}$ | $\frac{\hat{x}_{21k}\,\hat{x}_{22k}}{r_{2k}}\big{(}1+\frac{r_{2k}-1}{x_{2+k}}\big{)}$
2 sampling steps, drift |  $\frac{\hat{x}_{11k}\,\hat{x}_{12k}}{r_{1k}}\big{(}1+\frac{r_{1k}-1}{x_{1+k}}\big{)}$ | $r_{2k}\big{(}\hat{p}_{2k}(1-\hat{p}_{2k})\big{(}1+\frac{r_{2k}-1}{x_{2+k}}\big{)}+(r_{2k}-1)\,\frac{x_{2+k}-1}{x_{2+k}}\hat{\sigma}_{\text{drift-k}}^{2}\big{)}$
⋆ This is the situation of the classical CMH test.

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

$\frac{\frac{x_{11}}{x_{1+}}+\frac{x_{21}}{x_{2+}}}{2}\quad\text{ or
}\quad\frac{\frac{\hat{x}_{11}}{r_{1}}+\frac{\hat{x}_{21}}{r_{2}}}{2}\text{
respectively}.$ (6)

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

$\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}.$ (7)

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

$\frac{\sum_{i=1}^{\gamma}\hat{f}_{i}}{\gamma}.$ (8)

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

$\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{)}.$
(9)

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.

(a) $T_{\text{CMH}}$

(b) $T_{\text{CMH}}^{\text{2s\&d}}$

Figure 2: Distribution of the p-values belonging to the classical CMH test
$T_{\text{CMH}}$ (a) and $T_{\text{CMH}}^{\text{2s\&d}}$ (b). Red line
indicates frequency of $5\%$. Simulation setup: $10^{6}$ neutral loci with
true allele frequencies in base population uniformly distributed on $[0,1]$,
$N_{e}=300$, allele sample size 1000, pool sequencing with coverage Poisson
distributed with mean 80, sequence data for generations 0 and 60, 3 replicate
populations.

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.

(a) $T_{\chi^{2}}^{\text{2s\&d}}$

(b) $T_{\text{CMH}}^{\text{2s\&d}}$

Figure 3: Percentage of loci with p-value smaller than cut-off value against
cut-off value for different minimum values of the allele frequencies of both
alleles in the base population, in (a) for $T_{\chi^{2}}^{\text{2s\&d}}$, and
in (b) for $T_{\text{CMH}}^{\text{2s\&d}}$. The black line is the angle
bisector. Simulation setup: $10^{6}$ neutral loci with true allele frequencies
in base population uniformly distributed on $[0,1]$, $N_{e}=300$, allele
sample size 1000, pool sequencing with coverage Poisson distributed with mean
80, sequence data for generations 0 and 60, 5 replicate populations in (b).

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

(a) $T_{\chi^{2}}^{\text{2s\&d}}$, 60 generations

(b) $T_{\text{CMH}}^{\text{2s\&d}}$, 60 generations

(c) $T_{\chi^{2}}^{\text{2s\&d}}$, 20 generations

(d) $T_{\chi^{2}}^{\text{2s\&d}}$, 200 generations

Figure 4: Power vs selection coefficients for different values of the
effective population size $N_{e}$, in (a), (c), (d) for
$T_{\chi^{2}}^{\text{2s\&d}}$, and in (b) for $T_{\text{CMH}}^{\text{2s\&d}}$.
Simulation setup: $10^{5}$ loci for each selection coefficient with true
allele frequencies in base population uniformly distributed on $[0,1]$, allele
sample size 1000, pool sequencing with coverage 100, 5 replicate populations
in (b); (a) and (b) sequence data generations 0 and 60, (c) sequence data
generations 0 and 20, (d) sequence data generations 0 and 200.

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.

(a) $T_{\chi^{2}}^{\text{2s\&d}}$

(b) $T_{\text{CMH}}^{\text{2s\&d}}$

Figure 5: Power vs sample size of alleles (pool size) for different pool
sequencing coverages, in (a) for $T_{\chi^{2}}^{\text{2s\&d}}$, and in (b) for
$T_{\text{CMH}}^{\text{2s\&d}}$. Simulation setup: $10^{5}$ loci for each
value of the pool size with true allele frequencies in base population
uniformly distributed on $[0,1]$, $N_{e}=300$, selection coefficient 0.05,
sequence data for generations 0 and 60, 5 replicate populations in (b).

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.

Table 5: Comparison of type I error, power, and running times in seconds for
different tests of selection. Simulation setup: $10^{4}$ loci with true allele
frequencies in base population uniformly distributed on $[0,1]$, 10% of the
loci under selection, selection coefficients exponentially distributed with
mean 0.1, $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_{\chi^{2}}^{\text{2s\&d-ig}}$ and $T_{\text{CMH}}^{\text{2s\&d-ig}}$ are
applied, 5 replicate populations in case of $T_{\text{CMH}}^{\text{2s\&d}}$
and $T_{\text{CMH}}^{\text{2s\&d-ig}}$.

Test | Type I Error | Power | Time
---|---|---|---
$T_{\chi^{2}}^{\text{2s\&d}}$ | 0.050 | 0.417 | 0.005
$T_{\chi^{2}}^{\text{2s\&d-ig}}$ | 0.046 | 0.490 | 0.369
$T_{\text{CMH}}^{\text{2s\&d}}$ with 5 replicates | 0.050 | 0.761 | 0.014
$T_{\text{CMH}}^{\text{2s\&d-ig}}$ 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

⋆ Assume diploids, dominance is set to 0.5, the method to estimate selection
is set to ”LLS” and p-values are simulated with N.pval set to 1000, which
means that 1000 simulations are performed to estimate the p-values.

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

(a) $T_{\chi^{2}}^{\text{2s\&d-ig}}$

(b) $T_{\text{CMH}}^{\text{2s\&d-ig}}$

Figure 6: Distribution of the p-values belonging to
$T_{\chi^{2}}^{\text{2s\&d-ig}}$ (a) and $T_{\text{CMH}}^{\text{2s\&d-ig}}$
(b). Red line indicates frequency of $5\%$. Simulation setup: $10^{6}$ neutral
loci with true allele frequencies in base population beta distributed with
parameters 0.2 and 0.2, $N_{e}=300$, allele sample size 1000, pool sequencing
with coverage Poisson distributed with mean 80, sequence data every 10
generations from generation 0 to 60, 3 replicate populations in (b).

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

Figure 7: Manhatten plot for $T_{\text{CMH}}^{\text{1s\&d}}$ applied to real
data from Drosophila simulans from (Barghi et al., 2017).

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

## References

  * Agresti (2002) [author] Agresti, AlanA. (2002). Categorical Data Analysis. John Wiley & Sons.
  * Barghi et al. (2017) [author] Barghi, NedaN., Tobler, RaymondR., Nolte, ViolaV. and Schlötterer, ChristianC. (2017). _Drosophila Simulans_ : A Species with Improved Resolution in Evolve and Resequence Studies. G3: Genes, Genomes, Genetics 7 2337-2343. 10.1534/g3.117.043349
  * Bastide et al. (2013) [author] Bastide, HéloïseH., Betancourt, AndreaA., Nolte, ViolaV., Tobler, RaymondR., Stöbe, PetraP., Futschik, AndreasA. and Schlötterer, ChristianC. (2013). A Genome-Wide, Fine-Scale Map of Natural Pigmentation Variation in _Drosophila_ _Melanogaster_. PLoS Genetics 9 e1003534. 10.1371/journal.pgen.1003534
  * Bollback, York and Nielsen (2008) [author] Bollback, J. P.J. P., York, T. L.T. L. and Nielsen, R.R. (2008). Estimation of 2Nes from Temporal Allele Frequency Data. Genetics 179 497-502. 10.1534/genetics.107.085019
  * Burke et al. (2010a) [author] Burke, Molly K.M. K., Dunham, Joseph P.J. P., Shahrestani, ParvinP., Thornton, Kevin R.K. R., Rose, Michael R.M. R. and Long, Anthony D.A. D. (2010a). Genome-Wide Analysis of a Long-Term Evolution Experiment with _Drosophila_. Nature 467 587-590. 10.1038/nature09352
  * Burke et al. (2010b) [author] Burke, Molly KM. K., Dunham, Joseph PJ. P., Shahrestani, ParvinP., Thornton, Kevin RK. R., Rose, Michael RM. R. and Long, Anthony DA. D. (2010b). Genome-wide analysis of a long-term evolution experiment with Drosophila. Nature 467 587\.
  * Endler et al. (2016) [author] Endler, L.L., Betancourt, A. J.A. J., Nolte, V.V. and Schlötterer, C.C. (2016). Reconciling Differences in Pool-Gwas between Populations: A Case Study of Female Abdominal Pigmentation in _Drosophila_ _Melanogaster_. Genetics 202 843-855. 10.1534/genetics.115.183376
  * Ewens (2004) [author] Ewens, Warren J.W. J. (2004). Mathematical Population Genetics: 1. Theoretical Introduction. Springer Science + Business Media.
  * Falconer (1960) [author] Falconer, Douglas S.D. S. (1960). Introduction to Quantitative Genetics. The Ronald Press Company.
  * Feder, Kryazhimskiy and Plotkin (2014) [author] Feder, A. F.A. F., Kryazhimskiy, S.S. and Plotkin, J. B.J. B. (2014). Identifying Signatures of Selection in Genetic Time Series. Genetics 196 509-522. 10.1534/genetics.113.158220
  * Foll, Shim and Jensen (2015) [author] Foll, MatthieuM., Shim, HyunjinH. and Jensen, Jeffrey D.J. D. (2015). WFABC: A Wright-Fisher ABC-Based Approach for Inferring Effective Population Sizes and Selection Coefficients from Time-Sampled Data. Molecular Ecology Resources 15 87-98. 10.1111/1755-0998.12280
  * Griffin et al. (2017) [author] Griffin, Philippa C.P. C., Hangartner, Sandra B.S. B., Fournier-Level, AlexandreA. and Hoffmann, Ary A.A. A. (2017). Genomic Trajectories to Desiccation Resistance: Convergence and Divergence among Replicate Selected _Drosophila_ Lines. Genetics 205 871-890. 10.1534/genetics.116.187104
  * Illingworth et al. (2012) [author] Illingworth, C. J. R.C. J. R., Parts, L.L., Schiffels, S.S., Liti, G.G. and Mustonen, V.V. (2012). Quantifying Selection Acting on a Complex Trait Using Allele Frequency Time Series Data. Molecular Biology and Evolution 29 1187-1197. 10.1093/molbev/msr289
  * Iranmehr et al. (2017) [author] Iranmehr, AryaA., Akbari, AliA., Schlötterer, ChristianC. and Bafna, VineetV. (2017). CLEAR: Composition of Likelihoods for Evolve and Resequence Experiments. Genetics 1011-1023. 10.1101/080085
  * Jónás et al. (2016) [author] Jónás, A.A., Taus, T.T., Kosiol, C.C., Schlötterer, C.C. and Futschik, A.A. (2016). Estimating the Effective Population Size from Temporal Allele Frequency Changes in Experimental Evolution. Genetics 204 723-735. 10.1534/genetics.116.191197
  * Kofler, Pandey and Schlötterer (2011) [author] Kofler, R.R., Pandey, R. V.R. V. and Schlötterer, C.C. (2011). PoPoolation2: Identifying Differentiation between Populations Using Sequencing of Pooled DNA Samples (Pool-Seq). Bioinformatics 27 3435-3436. 10.1093/bioinformatics/btr589
  * Kofler and Schlötterer (2014) [author] Kofler, RobertR. and Schlötterer, ChristianC. (2014). A Guide for the Design of Evolve and Resequencing Studies. Molecular Biology and Evolution 31 474-483. 10.1093/molbev/mst221
  * Lehmann and Romano (2005) [author] Lehmann, E. L.E. L. and Romano, Joseph P.J. P. (2005). Testing Statistical Hypotheses, 3rd ed ed. Springer Science + Business Media.
  * Levy et al. (2015) [author] Levy, Sasha F.S. F., Blundell, Jamie R.J. R., Venkataram, SandeepS., Petrov, Dmitri A.D. A., Fisher, Daniel S.D. S. and Sherlock, GavinG. (2015). Quantitative Evolutionary Dynamics Using High-Resolution Lineage Tracking. Nature 519 181-186. 10.1038/nature14279
  * Lohr (2009) [author] Lohr, Sharon L.S. L. (2009). Sampling: Design and Analysis. Nelson Education.
  * Malaspinas et al. (2012) [author] Malaspinas, A. S.A. S., Malaspinas, O.O., Evans, S. N.S. N. and Slatkin, M.M. (2012). Estimating Allele Age and Selection Coefficient from Time-Serial Data. Genetics 192 599-607. 10.1534/genetics.112.140939
  * Mathieson and McVean (2013) [author] Mathieson, IainI. and McVean, GilG. (2013). Estimating Selection Coefficients in Spatially Structured Populations from Time Series Data of Allele Frequencies. Genetics 193 973-984. 10.1534/genetics.112.147611
  * McDonald (2014) [author] McDonald, John H.J. H. (2014). Handbook of Biological Statistics, 3rd ed. Sparky House Publishing.
  * Nouhaud et al. (2016) [author] Nouhaud, PierreP., Tobler, RayR., Nolte, ViolaV. and Schlötterer, ChristianC. (2016). Ancestral Population Reconstitution from Isofemale Lines as a Tool for Experimental Evolution. Ecology and Evolution 6 7169-7175. 10.1002/ece3.2402
  * Orozco-terWengel et al. (2012) [author] Orozco-terWengel, PabloP., Kapun, MartinM., Nolte, ViolaV., Kofler, RobertR., Flatt, ThomasT. and Schlötterer, ChristianC. (2012). Adaptation of _Drosophila_ to a Novel Laboratory Environment Reveals Temporally Heterogeneous Trajectories of Selected Alleles. Molecular Ecology 21 4931-4941. 10.1111/j.1365-294X.2012.05673.x
  * R-Core-Team (2018) [author] R-Core-Team (2018). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing.
  * Remolina et al. (2012) [author] Remolina, Silvia C.S. C., Chang, Peter L.P. L., Leips, JeffJ., Nuzhdin, Sergey V.S. V. and Hughes, Kimberly A.K. A. (2012). Genomic Basis of Aging and Life-History Evolution in _Drosophila_ _Melanogaster_. Evolution 66 3390-3403. 10.1111/j.1558-5646.2012.01710.x
  * Schraiber, Evans and Slatkin (2016) [author] Schraiber, Joshua G.J. G., Evans, Steven N.S. N. and Slatkin, MontgomeryM. (2016). Bayesian Inference of Natural Selection from Allele Frequency Time Series. Genetics 203 493-511. 10.1534/genetics.116.187278
  * Steinrücken, Bhaskar and Song (2014) [author] Steinrücken, MatthiasM., Bhaskar, AnandA. and Song, Yun S.Y. S. (2014). A Novel Spectral Method for Inferring General Diploid Selection from Time Series Genetic Data. The Annals of Applied Statistics 8 2203-2222. 10.1214/14-AOAS764
  * Taus, Futschik and Schlötterer (2017) [author] Taus, ThomasT., Futschik, AndreasA. and Schlötterer, ChristianC. (2017). Quantifying Selection with Pool-Seq Time Series Data. Molecular Biology and Evolution 34 3023-3034. 10.1093/molbev/msx225
  * Terhorst, Schlötterer and Song (2015) [author] Terhorst, JonathanJ., Schlötterer, ChristianC. and Song, Yun S.Y. S. (2015). Multi-Locus Analysis of Genomic Time Series Data from Experimental Evolution. PLOS Genetics 11 e1005069. 10.1371/journal.pgen.1005069
  * Tobler et al. (2014) [author] Tobler, RayR., Franssen, Susanne U.S. U., Kofler, RobertR., Orozco-terWengel, PabloP., Nolte, ViolaV., Hermisson, JoachimJ. and Schlötterer, ChristianC. (2014). Massive Habitat-Specific Genomic Response in _D_ _. Melanogaster_ Populations during Experimental Evolution in Hot and Cold Environments. Molecular Biology and Evolution 31 364-375. 10.1093/molbev/mst205
  * Topa et al. (2015) [author] Topa, HandeH., Jónás, ÁgnesA., Kofler, RobertR., Kosiol, CarolinC. and Honkela, AnttiA. (2015). Gaussian Process Test for High-Throughput Sequencing Time Series: Application to Experimental Evolution. Bioinformatics 31 1762-1770.
  * Turner and Miller (2012) [author] Turner, T. L.T. L. and Miller, P. M.P. M. (2012). Investigating Natural Variation in _Drosophila_ Courtship Song by the Evolve and Resequence Approach. Genetics 191 633-642. 10.1534/genetics.112.139337
  * Turner et al. (2011) [author] Turner, Thomas L.T. L., Stewart, Andrew D.A. D., Fields, Andrew T.A. T., Rice, William R.W. R. and Tarone, Aaron M.A. M. (2011). Population-Based Resequencing of Experimentally Evolved Populations Reveals the Genetic Basis of Body Size Variation in _Drosophila_ _Melanogaster_. PLoS Genetics 7 e1001336. 10.1371/journal.pgen.1001336
  * Waples (1989) [author] Waples, Robin S.R. S. (1989). A Generalized Approach for Estimating Effective Population Size from Temporal Changes in Allele Frequency. Genetics 121 379-391.

## 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

$X_{11}\sim\operatorname{Binom}(x_{1+},p_{1})\quad\text{ and }\quad
X_{21}\sim\operatorname{Binom}(x_{2+},p_{2}).$ (A1)

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

$T_{n}:=(X_{11}-\frac{x_{1+}\,X_{+1}}{n}).$ (A2)

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

$\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}$
(A3)

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

$\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}).$
(A4)

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

$\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}$
(A5)

With the central limit theorem we can thus deduce

$\frac{T_{n}}{\sqrt{\operatorname{Var}(T_{n})}}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}\mathcal{N}(0,1),$ (A6)

and with the continuous mapping theorem we conclude

$\frac{T_{n}^{2}}{\operatorname{Var}(T_{n})}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}\chi^{2}_{1},$ (A7)

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

$\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}.$ (A8)

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

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

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

$\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}$ (A10)

###  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),

$\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}),$
(A11)

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

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

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

$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},$ (A13)

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

$X_{11}\sim\operatorname{Binom}(x_{1+},p_{1})\qquad\text{and}\qquad
X_{21}\sim\operatorname{Binom}(x_{2+},P_{2}),$ (A14)

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

$\operatorname{Var}(X_{11})=x_{1+}\,p_{1}\,(1-p_{1})$ (A15)

and

$x_{1+}\,\frac{X_{11}}{x_{1+}}\,\frac{X_{12}}{x_{1+}}=\frac{X_{11}\,X_{12}}{x_{1+}}\stackrel{{\scriptstyle\mathbb{P}}}{{\longrightarrow}}\operatorname{Var}(X_{11}).$
(A16)

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

$\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}$
(A17)

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

$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}).$
(A18)

####  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

$K_{1}\sim\operatorname{Binom}(g_{1},b_{0})\qquad\text{and}\qquad
K_{2}\sim\operatorname{Binom}(g_{2},B_{1})\quad\text{with
}B_{1}=\frac{K_{1}}{g_{1}}.$ (A19)

With the law of total variance we obtain

$\displaystyle\operatorname{Var}(K_{2})$
$\displaystyle=\mathbb{E}[\operatorname{Var}(K_{2}\,|\,K_{1})]+\operatorname{Var}(\mathbb{E}[K_{2}\,|\,K_{1}])$
$\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{)}$
(A20)
$\displaystyle=g_{2}\,b_{0}\,(1-b_{0})\Big{(}1+\frac{g_{2}-1}{g_{1}}\Big{)}.$
(A21)

$\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

$\frac{\hat{X}_{11}\,\hat{X}_{12}}{r_{1}}\Big{(}1+\frac{r_{1}-1}{x_{1+}}\Big{)}\stackrel{{\scriptstyle\mathbb{P}}}{{\longrightarrow}}\operatorname{Var}(\hat{X}_{11}).$
(A22)

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

$\frac{\hat{X}_{21}\,\hat{X}_{22}}{r_{2}}\Big{(}1+\frac{r_{2}-1}{x_{2+}}\Big{)}\stackrel{{\scriptstyle\mathbb{P}}}{{\longrightarrow}}\operatorname{Var}(\hat{X}_{21}).$
(A23)

####  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

$\mathbb{E}[X_{21}]=\mathbb{E}[\mathbb{E}[X_{21}\,|\,P_{2}]]=x_{2+}\,\mathbb{E}[P_{2}].$
(A24)

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

$\begin{split}\operatorname{Var}(\hat{X}_{21})&=\frac{r_{2}}{x_{2+}}\,\Big{(}\mathbb{E}[X_{21}]\Big{(}1-\frac{\mathbb{E}[X_{21}]}{x_{2+}}\Big{)}+\frac{r_{2}-1}{x_{2+}}\operatorname{Var}(X_{21})\Big{)}\\\
&=r_{2}\Big{(}\mathbb{E}[p_{2}](1-\mathbb{E}[p_{2}])\big{(}1+\frac{r_{2}-1}{x_{2+}}\big{)}+(r_{2}-1)\,\frac{x_{2+}-1}{x_{2+}}\operatorname{Var}(P_{2})\Big{)}.\end{split}$
(A25)

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

$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{)}\stackrel{{\scriptstyle\mathbb{P}}}{{\longrightarrow}}\operatorname{Var}(\hat{X}_{21}).$
(A26)

###  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

(a) correction function

(b) correction result

Figure A1: (a) Function $F_{t}(\cdot)$ used to transform small p-values. (b)
Percentage of loci with p-value smaller than cut-off value against cut-off
value for $T^{2s\&d}_{\chi^{2}}$. The black line is the angle bisector.
Simulation setup: $10^{6}$ neutral loci with true allele frequencies in base
population uniformly distributed on $[0,1]$, $Ne=300$, allele sample size
1000, pool sequencing with coverage Poisson distributed with mean 80, sequence
data for generations 0 and 60.

(a) Effective population size $N_{e}$

(b) Pool size

(c) Coverage

(d) Number of generations

(e) Starting allele frequency

(f) Allele frequency range

Figure A2: Percentage of loci with p-value smaller than cut-off value against
cut-off value for $T^{2s\&d}_{\chi^{2}}$ with the correction for the
significant p-values. The black line is the angle bisector. Simulation setup:
For each plot one parameters varies according to the legend and the others
follow the following setup. $10^{6}$ neutral loci with true allele frequencies
in base population uniformly distributed on $[0,1]$, $Ne=300$, allele sample
size 1000, pool sequencing with coverage Poisson distributed with mean 80,
sequence data for generations 0 and 60.

(a) $T_{\chi^{2}}$

(b) $T_{\chi^{2}}^{\text{2s\&d}}$

Figure A3: Distribution of the p-values belonging to the classical chi-square
test $T_{\chi^{2}}$ (a) and $T_{\chi^{2}}^{\text{2s\&d}}$ (b). Red line
indicates frequency of $5\%$. Simulation setup: $10^{6}$ neutral loci with
true allele frequencies in base population uniformly distributed on $[0,1]$,
$N_{e}=300$, allele sample size 1000, pool sequencing with coverage Poisson
distributed with mean 80, sequence data for generations 0 and 60.

(a) $T_{\chi^{2}}^{\text{1s\&d}}$

(b) $T_{\text{CMH}}^{\text{1s\&d}}$

Figure A4: Distribution of the p-values belonging to
$T_{\chi^{2}}^{\text{1s\&d}}$ (a) and $T_{\text{CMH}}^{\text{1s\&d}}$ (b). Red
line indicates frequency of $5\%$. Simulation setup: $10^{6}$ neutral loci
with true allele frequencies in base population uniformly distributed on
$[0,1]$, $N_{e}=300$, pool sequencing of whole population with coverage
Poisson distributed with mean 80, sequence data for generations 0 and 60, 3
replicate populations in (b).

(a) $T_{\chi^{2}}$

(b) chi-square test adapted to 2 sampling steps

(c) $T_{\text{CMH}}$

(d) CMH test adapted to 2 sampling steps

Figure A5: Distribution of the p-values belonging to the classical chi-square
test $T_{\chi^{2}}$ (a), the chi-square test adapted to 2 sampling steps (b),
the classical CMH test $T_{\text{CMH}}$ (c), the CMH test adapted to 2
sampling steps (d). Red line indicates frequency of $5\%$. Simulation setup:
$10^{6}$ neutral loci with true allele frequencies in population 1 uniformly
distributed on $[0,1]$ and the same true frequencies in population 2 (no
drift), allele sample size 1000, pool sequencing with coverage Poisson
distributed with mean 80, 3 replicate populations in (c) and (d).

(a) $T_{\chi^{2}}^{\text{2s\&d-ig}}$

(b) $T_{\text{CMH}}^{\text{2s\&d-ig}}$

Figure A6: Percentage of loci with p-value smaller than cut-off value against
cut-off value in (a) for $T_{\chi^{2}}^{\text{2s\&d-ig}}$, and in (b) for
$T_{\text{CMH}}^{\text{2s\&d-ig}}$. The black line indicates the angle
bisector. Simulation setup: $10^{6}$ neutral loci with true allele frequencies
in base population uniformly distributed on [0,1], Ne = 300, allele sample
size 1000, pool sequencing with coverage Poisson distributed with mean 80,
sequence data for generations 0, 10, 20, 30, 40, 50, 60; 5 replicate
populations in (b).

(a) 20 generations of drift

(b) 200 generations of drift

Figure A7: Power of $T_{\text{CMH}}^{\text{2s\&d}}$ vs selection coefficients
for different values of the effective population size $N_{e}$. Simulation
setup: $10^{5}$ loci for each selection coefficient with true allele
frequencies in base population uniformly distributed on $[0,1]$, allele sample
size 1000, pool sequencing with coverage 100, 5 replicate populations. (a)
sequence data generations 0 and 20, (b) sequence data generations 0 and 200.

(a) $T_{\chi^{2}}^{\text{2s\&d}}$

(b) $T_{\text{CMH}}^{\text{2s\&d}}$

Figure A8: Power of $T_{\chi^{2}}^{\text{2s\&d}}$ (a) and
$T_{\text{CMH}}^{\text{2s\&d}}$ (b) vs selection coefficients for different
values of the effective population size $N_{e}$. Simulation setup: $10^{5}$
loci for each selection coefficient with true allele frequencies in base
population beta distributed with parameters 0.2 and 0.2, allele sample size
1000, pool sequencing with coverage 100, sequence data for generations 0 and
60, 5 replicate populations in (b).