BACKGROUND
in the past decade there has been considerable interest in detecting natural selection in humans and other organisms  <cit>  from dna sequence data.
an often used approach for detecting selection is to use a neutrality test statistic based on allele frequencies, with tajima’s d being the most famous.
such frequency based tests have been used to identify a number of genes, that have undergone selection: lactase  <cit> , abo blood group  <cit> , the hla immune complex  <cit> , and the calcium pathway trpv <dig>  <cit> .
in recent years next-generation sequencing  has revolutionized the field of genetics in general and population genetics in particular.
the underlying technology behind these different high-throughput sequencers are different  <cit> , but common is the unprecedented level of data produced, allowing population geneticists access to full genome data for large population samples at an affordable price.
applying the neutrality tests directly to the genotypes obtained from snp chips can lead to biased results because of the way the snps were selected.
even though work has been done to correct for this  <cit> , corrections are not possible without knowledge of the selection criteria used to ascertain snps for the specific chip.
ngs data do not suffer from such ascertainment biases.
however, with low and medium depth , genotype calling from ngs data can lead to other biases  <cit> .
the biases will depend on the method used for genotype calling.
early genotype calling methods were based on simple criteria, such as calling a heterozygous site, if the proportion of non-reference alleles is between 20% and 80%  <cit> .
most commonly used genotype callers today are probabilistic, i.e.
they are based on genotype likelihoods   <cit> .
a gl is the probability of the observed sequencing data given a genotype.
gls are obtained from the raw ngs read data, typically using information such as position in the read, quality score, and possibly type of mutation  <cit> .
a standard approach for data analysis is to call genotypes, and then estimate neutrality statistics.
most genotype callers are relatively conservative and only call a site to be heterozygous if there is substantial evidence that it actually is heterozygous.
such methods will tend to underestimate the allele frequency of the minor allele.
less conservative methods will tend to produce too many rare variants.
in either case, an unbiased estimate of the allele frequency cannot be obtained for all classes of allele frequencies  <cit> .
in some studies, an attempt to alleviate this problem has made by only calling a genotype when there is substantial statistical evidence supporting the genotype call.
such approaches will generate a considerable amount of missing data which leads to biases if not adequately dealt with  <cit> .
the effect of these errors on downstream analyses can be severe.
for most data produced, the strongest bias has been an excess of singletons .
if just a single individual in a sample is miscalled to be heterozygous in a monomorphic site, the site will appear as polymorphic singleton site.
the effects of genotype calling for downstream population genetic analyses have been extensively illustrated using simulations by  <cit> .
other approaches have been suggested: achaz,  <dig>  <cit>  proposed modified statistics that excludes the singletons to address the problem of high error rate in sequencing data.
liu et al.,  <dig>  <cit>  does not use the quality scores but uses an error rate for each nucleotide for all samples and  <cit>  which uses a fixed error rate.
this is contrasted by the methods of  <cit>  that can use quality scores for snp discovery and incorporating the quality scores into the parameter estimation directly.
we will here show, using both simulated and real human data that, that tajima’s d calculated using genotypes called from ngs data can lead to severely biased results.
the level of bias depends on the sequencing depths and error rates, but disturbingly it also depends on whether or not the data set is neutral or not.
this in effect leads to lower power to detect selection.
we argue, that the solution to this problem is to avoid genotype calling, but instead incorporate uncertainty in the genotypes through direct analyses of the genotype likelihood .
this approach also solves the issue of missing data, and the problem of inferring which sites are variable.
we propose two methods for estimating the neutrality tests statistics from low depth ngs data.
in both methods the estimates are based on gls instead of called genotypes.
the first approach is based on obtaining a maximum likelihood  estimate of the sample allele frequency spectrum.
this approach takes all the uncertainty of the ngs data into account, and provides estimates of the neutrality tests statistics from the estimated sample frequency spectrum.
however, it is too slow for genome-wide window based scans.
the second approach is computationally feasible for entire genomes of any magnitude.
it uses an empirical bayes approach, that also uses an ml estimate of the site frequency spectrum.
this estimate is used as a prior for calculating site specific posterior probabilities for the sample frequency spectrum.
we show, using both simulated and real low depth ngs data, that the test statistics estimated using the empirical bayes approach are at most weakly biased and provide a computationally attractive alternative to genotype calling methods.
methods
estimators of θ and neutrality tests
in a set of n aligned dna sequences, the frequency spectrum  is defined as the vector η = k =  <dig> …,n where ηk is the number of sites with k derived alleles.
a derived allele is a new mutation and the ancestral allele is the nucleotide occurring in the site before mutation.
notice that these population genetic models assume an infinite sites model, and there is therefore an unambiguous definition of ancestral and derived.
in real data, outgroup information is used to estimate which allele is ancestral and which allele is derived.
using this notation, the number of segregating, i.e.
polymorphic, sites is then given by s=∑i=1n−1ηi.
the most commonly applied estimators of θ= <dig> nμ,  are linear functions of the frequency spectrum.
these estimators take the general form: θ^=∑i=0nαiηi, where αi are different vectors of weights used in constructing the estimators.
some of the common estimators include:
• θ^w=a1‒1∑i=1n−1ηi, with a1=∑i=1n−1i.
<cit>
• θ^π=n2−1∑i=1n−1in−iηi.
<cit>
• θ^fl=η <dig>  <cit>
• θ^h=n2−1∑i=1n−1i2ηi.
<cit>
for more background on these estimators, see for example  <cit> .
in brief, notice that the estimators differ in how they weight polymorphisms with different allele frequencies.
the θ^w, and θ^π estimators are symmetric in that they assign equal weight to ancestral and derived alleles of the same frequency.
however, θ^π assigns more weight to alleles segregating at intermediate frequencies while θ^w weights all categories equal.
in contrast, θ^fl only assigns positive weight to derived singletons and θ^h assigns more weight to derived alleles of high frequency.
the first two estimators can be calculated without knowledge of the ancestral by using the folded sfs, η*=ηk*k= <dig> …,n/ <dig> ηk*=ηk+ηn−k, assuming n is even.
the latter two estimators require knowledge regarding ancestral states, i.e.
they use the unfolded sfs .
none of these estimators are maximum likelihood estimators .
however, they are commonly used to construct tests statistics used to test if a predicted sfs fits an observed sfs.
in a standard neutral model e = θ/i <cit> .
selection, as well as other violations of the model assumptions such as population structure or population growth, will cause deviations from this expectation  <cit> .
the sfs based neutrality tests capture this effect by the use of test statistics that are constructed as the standardized difference between two different θ estimators.
they have the general form:
t=θ1‒θ2varθ1‒θ <dig>
where the variance usually is the variance expected under a standard neutral model without recombination.
for example, tajima’s d is given by setting θ^1=θ^π and θ^2=θ^w.
tajima’s d is often used for detecting selective sweeps, i.e.
the effect of an advantageous mutation going to fixation in the population.
after fixation, of an advantageous allele, there is an excess of rare alleles in the population.
this will cause tajima’s d to become negative as the expectation of θ^π will be smaller than the expectation of θ^w, when the sfs contains relatively more rare alleles than expected under a standard neutral model.
genotype likelihoods
the standard method for representing uncertainty in ngs data is in terms of genotype likelihoods .
all methods discussed in this paper are based on such gls.
we, therefore, first introduce the gl calculations before describing the methods used for calculating neutrality statistics.
the gl for an individual in a particular site is defined as the probability of the observed read data  in the site given the genotype type  of the individual in that site:
lg=a <dig> a2|d∝prd|g=a <dig> a <dig> a <dig> a2∈a,c,g,t.
the gl can be calculated by assuming independence of the reads covering a position.
some methods take into account the position within the read and does recalibration which take into account cycle dependencies  <cit> .
the implementation in  <cit>  however assumes that all reads have the same length.
the widely used samtools uses a model derived from the maq program  <cit> .
these gl methods are all single sample based in contrast to the method of  <cit> , which estimates type specific errors for multi samples jointly.
in this paper we base our simulations on a simplified method of  <cit> .
model
full maximum likelihood
several recent papers have suggested methods for estimating the sfs from ngs data  <cit> .
one approach that we will explore for calculating neutrality statistics  for a region, is to estimate the sfs, η = k =  <dig> …,n, using one of these methods, in this case the maximum likelihood method by  <cit>  , and then calculate t using estimated values of η instead of the  true values.
as  has previously been shown to perform well under simple conditions, we would expect this approach for calculating the frequency spectrum to perform well under similar conditions.
we note that the uncertainty regarding genotypes and the effect of sequencing and mapping errors are incorporated in this approach through the calculation of gls as described in the section above, but we will otherwise refer readers to  <cit>  for information regarding algorithmic details.
neutrality tests are normally performed using sliding windows in a genome.
maximum likelihood estimation of the sfs for all windows may be computationally challenging.
for this reason we also explore other approaches below.
empirical bayes
in a second approach we instead calculate the posterior probability of the allele frequencies for each site, and combine them into a joint estimate of the sfs as in  <cit> , i.e.
for each site we calculate:
prx=j|d=prd|x=jprx=j∑i=02nprd|x=iprx=i,
where x is the sample allele frequency in that site, and n being the number of diploid samples.
pr can be calculated fast using a dynamic programming algorithm described in detail in references  <cit> , and we refer readers to these publications for details.
the jth category of the sfs, is then estimated by summing pr over all sites in the region.
the prior, pr can be defined, for example from a ml estimate of the sfs from the entire genome, or from a reasonable subset of the genome.
the computational advantage of this approach is that we avoid optimization of all windows of the genome.
we consider this an empirical bayes approach as the prior is estimated from the data.
another advantage of this approach is that the estimates of θ calculated are actual posterior expectations of θ. this is true because all are linear functions of the ηi ’s, i.e.
eθ^|d=∑i=0nαieηi|d.
calling genotypes
to assess the performance of our method we will compare with a more simple method that estimates the sfs by calling genotypes for each individual using the same genotype likelihood calculations as used to estimate the sfs in the previous approaches.
we evaluate two commonly used gl-based calling procedures, that both assume that the major and minor allele are defined:
gmlike=argmaxg∈ <dig> ,2prd|g=g,ghwe=argmaxg∈ <dig> ,22gf^g1−f^2−gprd|g=g.
notice that we here, for notational convenience, code genotypes as elements of { <dig>   <dig>  2} instead of pairs of nucleotides that each are elements of {a, c, t, g}.
the first method chooses the genotype with the highest gl.
the second procedure  calls genotypes by first estimating the sample allele frequency f from a larger group of individuals, and then uses this estimate of f as a prior on the  <dig> genotypes using the assumption of hardy-weinberg equilibrium .
simply calling all genotypes would result in too many false heterozygotes due to the sequencing errors, resulting in too many variable sites and an inflation of the singleton category of the sfs.
therefore, we also need to select an inclusion criterion for the sites we deem variable, i.e.
sites in which we allow the existence of heterozygotes .
to do so we first perform a likelihood ratio test of the null hypothesis h0: f =  <dig>  combining the data from all individuals  <cit> .
we do so assuming that all sites are di-allelic, and use the approach described in the supplementary of  <cit>  to identify the two most likely alleles for a given site.
throughout we use the maximum likelihood estimator of f by  <cit> .
simulations
to assess the performance of our estimators we simulate genomic data under a standard neutral model and under models that deviate from the neutral model using msms  <cit> .
we set the population size, mutation rate and the recombination rate to realistic values for humans, n =  <dig>   <dig>  μ =  <dig>  · 10−  <dig>  r = 10−  <dig> <cit>  and use an infinite sites model.
as msms only prints out variable sites with binary coding, we insert invariable sites in the sequences and convert from binary coding to nucleotides by sampling randomly with equal probability from all four nucleotides.
we then collapse the simulated haplotypes into genotypes.
for the non-neutral scenarios we simulate strong positive selection under an additive model with a selection coefficient of  <dig> .
to simulate ngs data based on the msms genotypes, we use a model similar to the model assumed for gl calculations in gatk  <cit> .
the simplified gatk model uses only base quality information and calculates the gl for a single site as:
prd|g=a <dig> a2=∏i=1mprbi|g=a <dig> a2=∏i=1m12pbi|a1+12pbi|a <dig> pb|a=e3b≠a1−eb=a.
where m is the sequencing depth, bi the observed base in read i, e is the probability of error and g = {a <dig> a2}.
in our simulations, we assume an equal error rate for all bases and for all sites.
the sequencing depth is sampled from a poisson distribution with mean equal to the specified mean depth.
the simulations then proceed by first simulating g using msms and then simulating d in accordance with the formula given above.
the latter stage is achieved by sampling m nucleotides from g, and then introducing errors independently in each of them with probability e. for each site gls are then calculated according to the gatk model given above.
it should be clarified that we do not actually simulate reads, but sample bases for every site independently.
hence there is a dependency in sequencing depth between adjacent sites in real data that we do not model.
hapmap  <dig> and 1000genomes
for evaluating the performance of the estimators on real data we used  <dig> unrelated ceu samples from hapmap phase  <dig> data  <cit> , sequenced by the  <dig> genomes project  <cit>  using illumina sequencing.
based on the mapped reads we used angsd http://www.popgen.dk/angsd to align the  <dig> mapped samples and calculate the genotype likelihoods using the gatk error model.
ancestral states for all sites were obtained from the multiz46way dataset http://hgdownload.cse.ucsc.edu/goldenpath/hg19/multiz46way/ available from the ucsc browser.
RESULTS
the effect of genotype calling for low or medium coverage data
in order to evaluate the performance of the estimators we simulated multiple genomic regions both without selection and with strong positive selection.
we first attempted to identify an optimal p-value cutoff to use in the likelihood ratio tests applied in the genotype calling methods.
figure  <dig> depicts the distribution of the difference between the estimated and known value of tajima’s d for  <dig> different p-value cutoffs, along with the results using our eb approach.
each box represents the estimate from  <dig> 1 mb regions.
we simulated two neutral scenarios  and one scenario with strong positive selection  and we used a mean depth of  <dig> and  <dig> with an error rate of  <dig> % and  <dig> %.
for the scenario with selection we used an error rate of  <dig> %, but sampled the mean sequencing depth between the different samples from a poisson distribution with mean of  <dig>  for the same three scenarios we also estimated fu and li’s d and f statistics  <cit>  .
we see that it is not possible to choose a single p-value cutoff that is unbiased for all the examined scenarios .
any particular cut-off for genotype calling will result in biases in the estimate of tajima’s d that will depend on the true value of tajima’s d, error rates, and sequencing depth.
these results suggest that, for low or medium coverage data, using called genotypes is in general problematic no matter what cutoff is used and will give rise to biased results that are context dependent.
the effect of snp calling criteria when calling genotypes
to further examine to what extent the bias varies according to the p-value cutoff used in the lrt test for inferring variable sites, we summarized the distribution of estimated and known values in boxplots .
this was done for both genotype calling methods and for three different critical values .
there are some obvious biases due to snp calling.
less conservative snp calling will cause an excess of  low frequency variants, and therefore an underestimation of tajima’s d, and more conservative snp calling will cause a deficiency of rare alleles and, thereby, overestimate tajima’s d. also notice that for both genotype calling methods, the p-value threshold of 10- <dig> has less variance in the selection datasets compared to the neutral datasets, whereas the opposite trend is true with the more relaxed threshold of 5 × 10- <dig> .
in the following sections, we therefore compare our methods to results using several different cutoffs for genotype calling.
estimating tajima’s d using ml estimates of the sfs
we next simulated two neutral scenarios consisting of  <dig> samples and simulated  <dig> simulations of 1 mb regions for each.
for the first neutral scenario we simulated fairly high sequencing depth and assumed a high error rate .
for the second neutral scenario we simulated a lower depth and lower error rates .
we inferred tajima’s d values from the simulated haplotypes, which we will here denote as the true values.
values of tajima’s d were then estimated from the simulated sequencing data using our full ml approach and from genotypes called from the sequencing data .
for called genotypes we only included sites that were likely to be polymorphic with a p-value less than 10- <dig>  in both scenarios our approach gives very similar results to the estimates from the true haplotypes while the approaches based on genotype calling shows large biases as expected from the results presented in the previous section.
not surprisingly, we observe the least accurate results for the low depth scenario, but even at 8x coverage there are substantial biases.
the effect of sample size and sequencing depth
the effect of sequencing depth and error rate is further examined in figure  <dig>  here we show the distribution of the differences between the true and the estimated tajima’s d values.
as above, for each scenario we perform  <dig> simulations of 1 mb regions.
the first  <dig> sets of simulations are conducted with sample sizes of n =  <dig>  and the last set of  <dig> simulations are conducted with sample sizes of n =  <dig>  as can be seen, large error rates and sequencing depth affects the variance of the estimate for the full maximum likelihood method but the estimates remain approximately unbiased.
in contrast, the genotype calling  approaches show large biases that depend on sequencing depths and error rates.
similar results are observed for both sample sizes.
the mean simulated value  for three of the estimators of θ are shown in table  <dig>  the full likelihood method is mostly unbiased whereas the gc based methods are biased to a degree that depends on whether the region is under selection or not.
θ
w
θ
π
θ
fl
θ
w
θ
π
θ
fl
the genotype likelihoods are simulated with an error rate of  <dig> % and an average depth of  <dig> ×.
application to data simulated with a selective sweep
in typical applications to genome-wide data, tajima’s d will usually be calculated separately for multiple smaller regions, often in a sliding window.
we therefore separated the 1 mb regions into  <dig> subregions each of 50 kb.
we used both the full maximum likelihood method for each subregion and applied the empirical bayes  method.
the full maximum likelihood estimate from the whole 1 mb region was used as a prior for the eb method.
all of the methods show a decrease in tajima’s d values around the site under selection .
however, the threshold for inferring polymorphic sites has a strong influence on the performance of the gc methods.
note that because the full likelihood and the eb method do not introduce a bias in the estimate of tajima’s d, they do not have to be standardized using a neutral data set if the underlying demographic model is known.
in figure  <dig> we illustrate the distribution of the difference between the estimated tajima’s d and the true value for every window for the ml estimator and the eb estimator.
we note that the variance is larger for the full ml approach than the eb approach.
this is further examined in additional file 3: figure s <dig> where we have plotted the difference in mean squared error  for the same  <dig> subregions with the ml method and the eb method.
we notice that the mse for the eb method is smaller.
effect of the prior
for the eb method, the prior can have an impact on the estimates of tajima’s d. in figure  <dig> we explore the effect of the prior on the estimated values.
we use three different priors.
1) a prior estimated on a neutral data set of 100 mb, 2) a prior estimated for a 100 mb region under selection 3) a prior based on both types of regions in equal proportion.
as can be seen in figure  <dig> when applying the neutral prior to the neutral datasets or the selection prior on the selection dataset the estimates are approximately unbiased.
however, applying the “wrong prior” tends to either underestimate or overestimate the selection signal.
perhaps a bit surprisingly, using the neutral prior on the selection datasets tends to overestimate the signal, while using the selection prior on the neutral datasets tends to underestimate it.
this is explained by the fact that a region of selection will have less variability than a neutral region.
a prior trained on a neutral dataset will therefore have a higher level of variability than the true level of variability in a region of selection.
applying a neutral prior will therefore allow more variability, and due to sequencing errors these sites will mostly be singletons.
this excess of singletons will have the largest impact on θ^w, and we will therefore underestimate tajima’s d. the opposite argument explains why the selection prior will overestimate tajima’s d when applied to a neutral dataset.
we observe the same trend for other neutrality tests statistics that can be seen in additional file 4: figure s <dig> and additional file 5: figure s <dig>  however, the other investigated neutrality test statistics, fu & li’s d and f, are not estimated with the same accuracy as tajima’s d. for both selection and neutral data sets, we seem to underestimate the statistics with the eb approach and neutral prior .
from figure  we also see that we have problems estimating the true value in the region of the targeted locus using the 50 kb ml approach.
the difference is perhaps caused by the fact that fu & li’s statistics are based on a single category of the frequency spectrum, whereas tajima’s d is based on all categories.
power to detect a selective sweep
often the goal of the investigator is not to estimate the neutrality test statistics in an unbiased manner but instead the goal is to identify regions with extreme values.
to investigate our ability to discriminate between regions with selection and neutral regions, we show receiver operating characteristic  for the different approaches.
these are shown in additional file 6: figure s <dig>  for different depth, error rates and number of individuals.
when the depth is high, all methods perform almost as well as when the genotypes are known without error.
but when the depth decreases or the error rates increases, only the full maximum likelihood approach and the eb approach have similar power as the known genotypes.
as expected, the simulations show that our methods have improved power to discriminate between regions evolving neutrally and under positive selection as more samples are added .
the genotype calling methods perform worse when the error rate is increased and the depth is decreased, especially at low depth with a low p-value cutoff, while the ml method for all scenarios performs almost as well as if the true genotypes where known .
variable depth and variable error rates from ngs data
so far, we have used a simple simulation model assuming sequencing depths follow a poisson distribution and a fixed error rate.
to examine the robustness of our conclusions to these assumptions, we made an additional set of simulations using the observed distribution of quality scores and sequencing depths tabulated from bam files from the  <dig> genomes project .
as before we simulate  <dig> 1 mb regions with and without selection for  <dig> individuals, and apply our two genotype calling methods and the eb method to the simulated data .
for the genotype calling approaches we only included sites that were likely to be polymorphic with a p-value <10- <dig>  we observe results that are highly compatible with the previous results.
the eb method is approximately unbiased in both scenarios and have similarly small mean squared error .
in contrast, the genotype calling methods are biased in both scenarios, and we notice again that the bias depends on whether or not the data is simulated under a neutral or under a selection model.
we observe that the mse for the genotype calling methods are orders of magnitude larger than the eb method, and that the mse for the gc-mlike model is more than twice as high under the selection scenario as under the neutral scenario, whereas we observe the opposite trend for gc-hwe method.
application to real data
we also applied our methods on real ngs data from the  <dig> ceu individuals from the  <dig> genomes project .
bersaglieri et al.,  <dig>  found a strong signal of positive selection surrounding the lct region of chromosome  <dig> .
we estimated tajima’s d values in a sliding window across this chromosome, using a prior estimated from a 50 mb region on chromosome  <dig> .
to compare with earlier published results we used the tajima track <cit>  from the ucsc genome browser  <cit> , which depicts estimated values of tajima’s d from  <dig> individuals of european descent from the perlegen genotyping dataset  <cit> .
for our eb method we performed sliding windows analysis with different window sizes  all using a fixed step size of 10 kb.
we also compared the results for the naïve genotype calling methods using  <dig> different snp inclusion cutoff criteria’s .
notice that the overall estimate of tajima’s d is very positive for the snp data, most likely due to ascertainment biases  <cit> .
also notice that the lowest observed value of tajima’s d is the lct region for the eb approach while there are multiple regions with low tajima’s d values for the gc approaches and the snp chip data.
discussion
we have developed two methods to perform the neutrality tests on ngs data that take the uncertainty of the genotypes into account.
both methods perform better than using called genotypes and in most instances they are approximately unbiased.
the full ml method is computationally slow when applied in sliding windows at a genome-wide scale, which is why we also present a fast empirical bayes method with a prior that is estimated from the data itself, for example the entire genome, or a reasonable subset of the genome.
this makes the method computationally feasible for full genomic data of any magnitude and any windows size.
there is not a single obvious way to identify snp sites and call genotypes from ngs data.
here we have tried different approaches with different cutoffs.
regardless of method and chosen cutoff they all show a large bias in some or all simulated scenarios.
this is evident from the raw theta estimates , and the actual test statistics .
the level of bias varies between the different scenarios, not only for different depths and error rates, but it also depends on whether or not the data set is neutral or affected by selection .
the results from this study suggest/confirms it is not unproblematic to perform neutrality tests on genotypes called from low or medium coverage ngs data.
in contrast, both the ml and the eb approach give fairly accurate estimates for all the examined measures.
when applying a neutral genome-wide prior for our analysis, we observed only small deviations from the true values of tajima’s d even for very low depth data.
when applying the eb approach we did observe a small bias in the regions under selection when the prior was estimated from regions without selection .
however, the bias is always much smaller than the bias of the gc approach.
even though the eb approach can give small biases it can still have an advantage over the full ml approach.
when estimating the neutrality test statistics for small windows, we often obtained a few extreme outliers with positive values of tajima’s d for the ml approach .
since the eb method uses the entire 1 mb region as prior we do not see a similar problem with extreme positive outliers and the variance of the estimates is smaller overall.
if the goal is to identify regions under selection the smaller variance of the eb approach will give fewer regions with extreme values.
this is because regions with little data will increase the variance for the full ml approach while they will give results closer to the prior for the eb approach.
this problem could also be circumvented by using a sliding window approach with window sizes determined by using a fixed number of snps.
we applied the eb method to data from  <dig> individuals from the  <dig> genomes project, and observed trends similar to previous published results for the lct region in europeans.
for snp chip data from the same region, we observed a large over-representation of positive tajima’s d values, presumably due to ascertainment biases introduced by the snp-discovery and selection process from the snp-chip itself.
similar positive values were not observed using the eb method based on the ngs data.
for both the gc and eb methods, we observe negative values around the lct gene, however the estimates are not very extreme for the gc approach.
the eb approach is the only approach for which lct has the most extreme tajima’s d value.
the main advantage of the approaches presented here, is that they, in expectation, have at most a very weak dependence on sequencing depth.
this facilitates the use of genome-wide scans on genomes with varying sequencing depth without introducing biases.
the computational framework suggested here, based on the eb approach, provides a robust and computationally fast method for scanning a genome for regions with outlying or extreme frequency spectrum.
such a method should be of great use in years to come when analyzing population genomic data from a variety of different species.
CONCLUSIONS
in this paper we show through simulations that estimating neutrality test statistics using called genotypes can lead to highly biased result.
the bias is dependent not only on sequencing depth and error rate but interestingly the bias also depends on whether or not the region is under selection or not.
we have developed an empirical bayes method that can calculate the test statistic fast and efficiently.
this method circumvents the problem of snp discovery, genotype calling and missing data, which is a fundamental problem of ngs data.
this is done by working with genotype likelihoods, which contains all relevant information about the uncertainty of the data.
using this approach leads to approximately unbiased estimates in most instances.
availability
all methods discussed in this paper are freely available as part of the analyses of next generation sequencing data  package .
competing interest
the authors declare that they have no competing interests.
authors’ contributions
rn designed the eb model together with tsk and im.
aa helped with the design of the software package and bug checked early versions of the program.
tsk implemented the methods, carried out all analyses and simulations, and drafted the first version of the manuscript with rn.
the manuscript has been thoroughly edited by the remaining authors.
rn supervised the process.
all authors read and approved the final manuscript.
supplementary material
additional file 1: figure s1
the effect of genotype calling for low or medium coverage data using fu & li’s d. the difference between estimated and known fu&li’s d statistic for three different scenarios with  <dig> different p-value cutoffs.
each box is estimated on the basis of  <dig> 1 mb regions.
the top figure is based on genotypes called using the frequency as prior, and the bottom figure is based on genotypes called using a maximum likelihood approach.
notice that no single best cutoff can be chosen across the three different scenarios.
click here for file
additional file 2: figure s2
the effect of genotype calling for low or medium coverage data using fu & li’s f. the difference between estimated and known fu&li’s f statistic for three different scenarios with  <dig> different p-value cutoffs.
each box is estimated on the basis of  <dig> 1 mb regions.
the top figure is based on genotypes called using the frequency as prior, and the bottom figure is based on genotypes called using a maximum likelihood approach.
notice that no single best cutoff can be chosen across the three different scenarios.
click here for file
additional file 3: figure s3
difference in mean squared error  between the full ml and the eb method under a selective sweep.
mse of the estimated tajima’s d  is calculated for every 50 kb sub region of the full 1 mb region.
the figure is based on  <dig> 1 mb regions.
for the eb method we used a prior estimated from the entire 1 mb region on every 50 kb subregion.
click here for file
additional file 4: figure s4
effect of different priors for the eb method using the fu & li’s d. left and center plot are boxplots for the difference between our estimate of fu & li d statistics and the true value, these are based on 100 × 1 mb regions.
right plot is a 50 kb window plot using the 50 kb ml method along with the eb with neutral and mixed prior.
neutral prior is from a genome-wide prior based on a 100 mb region, neu + sel prior is based on a 200 mb prior based on 100 mb selection and 100 mb neutral.
click here for file
additional file 5: figure s5
effect of different priors for the eb method using the fu & li’s f. left and center plot are boxplots for the difference between our estimate of fu & li f statistics and the true value.
right plot is a 50 kb window plot using the 50 kb method along with the neutral and mixed prior.
neutral prior is from a genome-wide prior based on a 100 mb region, neu + sel prior is based on a 200 mb prior based on 100 mb selection and 100 mb neutral.
click here for file
additional file 6: figure s6
power to detect a selective sweep.
roc curve for scenarios, each plot is based on tajima’s d estimate for 100 × 1 mb regions with selection and  <dig> 1 mb regions without selection.
for each scenario we have our eb method along with our two genotype calling methods .
row <dig> is  <dig> individuals row <dig> is  <dig> individuals.
column <dig> is 8× 1% error rate, column2- <dig> is all 4×, but with varying error rates 1%, <dig> % and  <dig> %.
click here for file
additional file 7: figure s7
roc curve for low coverage dataset.
roc curve for a 2× <dig> % error rate.
the lrt criteria is 10- <dig>  this plot is based on  <dig> 1 mb regions with selection, and  <dig> 1 mb neutral regions.
click here for file
additional file 8: figure s8
distribution of quality scores and sequencing depth for a bam file.
the left panel shows the quality score distribution and right panel shows the depth distribution, tabulated for chr <dig> from a bam file from the  <dig> genomes project.
the mean quality score value was approximately  <dig> which corresponds to an average error rate of  <dig> %.
the data covered approximately 87% of the genome, had an average sequencing depth of  <dig> , and contained  <dig>  sites with a sequencing depth above  <dig>  the right panel only shows the first  <dig> observations.
click here for file
additional file 9: figure s9
using observed qscore and depth distributions.
boxplots of the difference between the estimate of tajima’s d and the known value for  <dig> 1 mb regions with our eb method and the two genotype calling methods.
the quality score and depth distributions used for the genotype likelihood calculations are based on the results depicted in figure s <dig>  for the genotype calling methods we used a cut-off for the p-value of the lrt test of 10- <dig>
click here for file
acknowledgements
work for this manuscript was financed through danish national research foundation,  villum foundation and danish council for independent research.

 in the past decade there has been considerable interest in detecting natural selection in humans and other organisms from dna sequence data. an often used approach for detecting selection is to use a neutrality test statistic based on allele frequencies, with tajima’s d being the most famous. such frequency based tests have been used to identify a number of genes, that have undergone selection: lactase, abo blood group, the hla immune complex, and the calcium pathway trpv. in recent years next-generation sequencing has revolutionized the field of genetics in general and population genetics in particular. the underlying technology behind these different high-throughput sequencers are different, but common is the unprecedented level of data produced, allowing population geneticists access to full genome data for large population samples at an affordable price. applying the neutrality tests directly to the genotypes obtained from snp chips can lead to biased results because of the way the snps were selected. even though work has been done to correct for this, corrections are not possible without knowledge of the selection criteria used to ascertain snps for the specific chip. ngs data do not suffer from such ascertainment biases. however, with low and medium depth, genotype calling from ngs data can lead to other biases. the biases will depend on the method used for genotype calling. early genotype calling methods were based on simple criteria, such as calling a heterozygous site, if the proportion of non-reference alleles is between 20% and 80%. most commonly used genotype callers today are probabilistic, i.e. they are based on genotype likelihoods. a gl is the probability of the observed sequencing data given a genotype. gls are obtained from the raw ngs read data, typically using information such as position in the read, quality score, and possibly type of mutation. a standard approach for data analysis is to call genotypes, and then estimate neutrality statistics. most genotype callers are relatively conservative and only call a site to be heterozygous if there is substantial evidence that it actually is heterozygous. such methods will tend to underestimate the allele frequency of the minor allele. less conservative methods will tend to produce too many rare variants. in either case, an unbiased estimate of the allele frequency cannot be obtained for all classes of allele frequencies. in some studies, an attempt to alleviate this problem has made by only calling a genotype when there is substantial statistical evidence supporting the genotype call. such approaches will generate a considerable amount of missing data which leads to biases if not adequately dealt with. the effect of these errors on downstream analyses can be severe. for most data produced, the strongest bias has been an excess of singletons. if just a single individual in a sample is miscalled to be heterozygous in a monomorphic site, the site will appear as polymorphic singleton site. the effects of genotype calling for downstream population genetic analyses have been extensively illustrated using simulations by. other approaches have been suggested: achaz, proposed modified statistics that excludes the singletons to address the problem of high error rate in sequencing data. liu et al., does not use the quality scores but uses an error rate for each nucleotide for all samples and which uses a fixed error rate. this is contrasted by the methods of that can use quality scores for snp discovery and incorporating the quality scores into the parameter estimation directly. we will here show, using both simulated and real human data that, that tajima’s d calculated using genotypes called from ngs data can lead to severely biased results. the level of bias depends on the sequencing depths and error rates, but disturbingly it also depends on whether or not the data set is neutral or not. this in effect leads to lower power to detect selection. we argue, that the solution to this problem is to avoid genotype calling, but instead incorporate uncertainty in the genotypes through direct analyses of the genotype likelihood. this approach also solves the issue of missing data, and the problem of inferring which sites are variable. we propose two methods for estimating the neutrality tests statistics from low depth ngs data. in both methods the estimates are based on gls instead of called genotypes. the first approach is based on obtaining a maximum likelihood estimate of the sample allele frequency spectrum. this approach takes all the uncertainty of the ngs data into account, and provides estimates of the neutrality tests statistics from the estimated sample frequency spectrum. however, it is too slow for genome-wide window based scans. the second approach is computationally feasible for entire genomes of any magnitude. it uses an empirical bayes approach, that also uses an ml estimate of the site frequency spectrum. this estimate is used as a prior for calculating site specific posterior probabilities for the sample frequency spectrum. we show, using both simulated and real low depth ngs data, that the test statistics estimated using the empirical bayes approach are at most weakly biased and provide a computationally attractive alternative to genotype calling methods. methods estimators of θ and neutrality tests in a set of n aligned dna sequences, the frequency spectrum is defined as the vector η = k =  …,n where ηk is the number of sites with k derived alleles. a derived allele is a new mutation and the ancestral allele is the nucleotide occurring in the site before mutation. notice that these population genetic models assume an infinite sites model, and there is therefore an unambiguous definition of ancestral and derived. in real data, outgroup information is used to estimate which allele is ancestral and which allele is derived. using this notation, the number of segregating, i.e. polymorphic, sites is then given by s=∑i=1n−1ηi. the most commonly applied estimators of θ= nμ, are linear functions of the frequency spectrum. these estimators take the general form: θ^=∑i=0nαiηi, where αi are different vectors of weights used in constructing the estimators. some of the common estimators include: • θ^w=a1‒1∑i=1n−1ηi, with a1=∑i=1n−1i. • θ^π=n2−1∑i=1n−1in−iηi. • θ^fl=η • θ^h=n2−1∑i=1n−1i2ηi. for more background on these estimators, see for example. in brief, notice that the estimators differ in how they weight polymorphisms with different allele frequencies. the θ^w, and θ^π estimators are symmetric in that they assign equal weight to ancestral and derived alleles of the same frequency. however, θ^π assigns more weight to alleles segregating at intermediate frequencies while θ^w weights all categories equal. in contrast, θ^fl only assigns positive weight to derived singletons and θ^h assigns more weight to derived alleles of high frequency. the first two estimators can be calculated without knowledge of the ancestral by using the folded sfs, η*=ηk*k= …,n/ ηk*=ηk+ηn−k, assuming n is even. the latter two estimators require knowledge regarding ancestral states, i.e. they use the unfolded sfs. none of these estimators are maximum likelihood estimators. however, they are commonly used to construct tests statistics used to test if a predicted sfs fits an observed sfs. in a standard neutral model e = θ/i. selection, as well as other violations of the model assumptions such as population structure or population growth, will cause deviations from this expectation. the sfs based neutrality tests capture this effect by the use of test statistics that are constructed as the standardized difference between two different θ estimators. they have the general form: t=θ1‒θ2varθ1‒θ where the variance usually is the variance expected under a standard neutral model without recombination. for example, tajima’s d is given by setting θ^1=θ^π and θ^2=θ^w. tajima’s d is often used for detecting selective sweeps, i.e. the effect of an advantageous mutation going to fixation in the population. after fixation, of an advantageous allele, there is an excess of rare alleles in the population. this will cause tajima’s d to become negative as the expectation of θ^π will be smaller than the expectation of θ^w, when the sfs contains relatively more rare alleles than expected under a standard neutral model. genotype likelihoods the standard method for representing uncertainty in ngs data is in terms of genotype likelihoods. all methods discussed in this paper are based on such gls. we, therefore, first introduce the gl calculations before describing the methods used for calculating neutrality statistics. the gl for an individual in a particular site is defined as the probability of the observed read data in the site given the genotype type of the individual in that site: lg=a a2|d∝prd|g=a a a a2∈a,c,g,t. the gl can be calculated by assuming independence of the reads covering a position. some methods take into account the position within the read and does recalibration which take into account cycle dependencies. the implementation in however assumes that all reads have the same length. the widely used samtools uses a model derived from the maq program. these gl methods are all single sample based in contrast to the method of, which estimates type specific errors for multi samples jointly. in this paper we base our simulations on a simplified method of. model full maximum likelihood several recent papers have suggested methods for estimating the sfs from ngs data. one approach that we will explore for calculating neutrality statistics for a region, is to estimate the sfs, η = k =  …,n, using one of these methods, in this case the maximum likelihood method by, and then calculate t using estimated values of η instead of the true values. as has previously been shown to perform well under simple conditions, we would expect this approach for calculating the frequency spectrum to perform well under similar conditions. we note that the uncertainty regarding genotypes and the effect of sequencing and mapping errors are incorporated in this approach through the calculation of gls as described in the section above, but we will otherwise refer readers to for information regarding algorithmic details. neutrality tests are normally performed using sliding windows in a genome. maximum likelihood estimation of the sfs for all windows may be computationally challenging. for this reason we also explore other approaches below. empirical bayes in a second approach we instead calculate the posterior probability of the allele frequencies for each site, and combine them into a joint estimate of the sfs as in, i.e. for each site we calculate: prx=j|d=prd|x=jprx=j∑i=02nprd|x=iprx=i, where x is the sample allele frequency in that site, and n being the number of diploid samples. pr can be calculated fast using a dynamic programming algorithm described in detail in references, and we refer readers to these publications for details. the jth category of the sfs, is then estimated by summing pr over all sites in the region. the prior, pr can be defined, for example from a ml estimate of the sfs from the entire genome, or from a reasonable subset of the genome. the computational advantage of this approach is that we avoid optimization of all windows of the genome. we consider this an empirical bayes approach as the prior is estimated from the data. another advantage of this approach is that the estimates of θ calculated are actual posterior expectations of θ. this is true because all are linear functions of the ηi ’s, i.e. eθ^|d=∑i=0nαieηi|d. calling genotypes to assess the performance of our method we will compare with a more simple method that estimates the sfs by calling genotypes for each individual using the same genotype likelihood calculations as used to estimate the sfs in the previous approaches. we evaluate two commonly used gl-based calling procedures, that both assume that the major and minor allele are defined: gmlike=argmaxg∈,2prd|g=g,ghwe=argmaxg∈,22gf^g1−f^2−gprd|g=g. notice that we here, for notational convenience, code genotypes as elements of { 2} instead of pairs of nucleotides that each are elements of {a, c, t, g}. the first method chooses the genotype with the highest gl. the second procedure calls genotypes by first estimating the sample allele frequency f from a larger group of individuals, and then uses this estimate of f as a prior on the genotypes using the assumption of hardy-weinberg equilibrium. simply calling all genotypes would result in too many false heterozygotes due to the sequencing errors, resulting in too many variable sites and an inflation of the singleton category of the sfs. therefore, we also need to select an inclusion criterion for the sites we deem variable, i.e. sites in which we allow the existence of heterozygotes. to do so we first perform a likelihood ratio test of the null hypothesis h0: f =  combining the data from all individuals. we do so assuming that all sites are di-allelic, and use the approach described in the supplementary of to identify the two most likely alleles for a given site. throughout we use the maximum likelihood estimator of f by. simulations to assess the performance of our estimators we simulate genomic data under a standard neutral model and under models that deviate from the neutral model using msms. we set the population size, mutation rate and the recombination rate to realistic values for humans, n =    μ =   · 10−   r = 10−  and use an infinite sites model. as msms only prints out variable sites with binary coding, we insert invariable sites in the sequences and convert from binary coding to nucleotides by sampling randomly with equal probability from all four nucleotides. we then collapse the simulated haplotypes into genotypes. for the non-neutral scenarios we simulate strong positive selection under an additive model with a selection coefficient of. to simulate ngs data based on the msms genotypes, we use a model similar to the model assumed for gl calculations in gatk. the simplified gatk model uses only base quality information and calculates the gl for a single site as: prd|g=a a2=∏i=1mprbi|g=a a2=∏i=1m12pbi|a1+12pbi|a pb|a=e3b≠a1−eb=a. where m is the sequencing depth, bi the observed base in read i, e is the probability of error and g = {a a2}. in our simulations, we assume an equal error rate for all bases and for all sites. the sequencing depth is sampled from a poisson distribution with mean equal to the specified mean depth. the simulations then proceed by first simulating g using msms and then simulating d in accordance with the formula given above. the latter stage is achieved by sampling m nucleotides from g, and then introducing errors independently in each of them with probability e. for each site gls are then calculated according to the gatk model given above. it should be clarified that we do not actually simulate reads, but sample bases for every site independently. hence there is a dependency in sequencing depth between adjacent sites in real data that we do not model. hapmap and 1000genomes for evaluating the performance of the estimators on real data we used unrelated ceu samples from hapmap phase data, sequenced by the genomes project using illumina sequencing. based on the mapped reads we used angsd http://www.popgen.dk/angsd to align the mapped samples and calculate the genotype likelihoods using the gatk error model. ancestral states for all sites were obtained from the multiz46way dataset http://hgdownload.cse.ucsc.edu/goldenpath/hg19/multiz46way/ available from the ucsc browser.  the effect of genotype calling for low or medium coverage data in order to evaluate the performance of the estimators we simulated multiple genomic regions both without selection and with strong positive selection. we first attempted to identify an optimal p-value cutoff to use in the likelihood ratio tests applied in the genotype calling methods. figure  depicts the distribution of the difference between the estimated and known value of tajima’s d for different p-value cutoffs, along with the results using our eb approach. each box represents the estimate from 1 mb regions. we simulated two neutral scenarios and one scenario with strong positive selection and we used a mean depth of and with an error rate of % and %. for the scenario with selection we used an error rate of %, but sampled the mean sequencing depth between the different samples from a poisson distribution with mean of for the same three scenarios we also estimated fu and li’s d and f statistics. we see that it is not possible to choose a single p-value cutoff that is unbiased for all the examined scenarios. any particular cut-off for genotype calling will result in biases in the estimate of tajima’s d that will depend on the true value of tajima’s d, error rates, and sequencing depth. these results suggest that, for low or medium coverage data, using called genotypes is in general problematic no matter what cutoff is used and will give rise to biased results that are context dependent. the effect of snp calling criteria when calling genotypes to further examine to what extent the bias varies according to the p-value cutoff used in the lrt test for inferring variable sites, we summarized the distribution of estimated and known values in boxplots. this was done for both genotype calling methods and for three different critical values. there are some obvious biases due to snp calling. less conservative snp calling will cause an excess of low frequency variants, and therefore an underestimation of tajima’s d, and more conservative snp calling will cause a deficiency of rare alleles and, thereby, overestimate tajima’s d. also notice that for both genotype calling methods, the p-value threshold of 10- has less variance in the selection datasets compared to the neutral datasets, whereas the opposite trend is true with the more relaxed threshold of 5 × 10-. in the following sections, we therefore compare our methods to results using several different cutoffs for genotype calling. estimating tajima’s d using ml estimates of the sfs we next simulated two neutral scenarios consisting of samples and simulated simulations of 1 mb regions for each. for the first neutral scenario we simulated fairly high sequencing depth and assumed a high error rate. for the second neutral scenario we simulated a lower depth and lower error rates. we inferred tajima’s d values from the simulated haplotypes, which we will here denote as the true values. values of tajima’s d were then estimated from the simulated sequencing data using our full ml approach and from genotypes called from the sequencing data. for called genotypes we only included sites that were likely to be polymorphic with a p-value less than 10- in both scenarios our approach gives very similar results to the estimates from the true haplotypes while the approaches based on genotype calling shows large biases as expected from the results presented in the previous section. not surprisingly, we observe the least accurate results for the low depth scenario, but even at 8x coverage there are substantial biases. the effect of sample size and sequencing depth the effect of sequencing depth and error rate is further examined in figure  here we show the distribution of the differences between the true and the estimated tajima’s d values. as above, for each scenario we perform simulations of 1 mb regions. the first sets of simulations are conducted with sample sizes of n =  and the last set of simulations are conducted with sample sizes of n =  as can be seen, large error rates and sequencing depth affects the variance of the estimate for the full maximum likelihood method but the estimates remain approximately unbiased. in contrast, the genotype calling approaches show large biases that depend on sequencing depths and error rates. similar results are observed for both sample sizes. the mean simulated value for three of the estimators of θ are shown in table  the full likelihood method is mostly unbiased whereas the gc based methods are biased to a degree that depends on whether the region is under selection or not. θ w θ π θ fl θ w θ π θ fl the genotype likelihoods are simulated with an error rate of % and an average depth of ×. application to data simulated with a selective sweep in typical applications to genome-wide data, tajima’s d will usually be calculated separately for multiple smaller regions, often in a sliding window. we therefore separated the 1 mb regions into subregions each of 50 kb. we used both the full maximum likelihood method for each subregion and applied the empirical bayes method. the full maximum likelihood estimate from the whole 1 mb region was used as a prior for the eb method. all of the methods show a decrease in tajima’s d values around the site under selection. however, the threshold for inferring polymorphic sites has a strong influence on the performance of the gc methods. note that because the full likelihood and the eb method do not introduce a bias in the estimate of tajima’s d, they do not have to be standardized using a neutral data set if the underlying demographic model is known. in figure  we illustrate the distribution of the difference between the estimated tajima’s d and the true value for every window for the ml estimator and the eb estimator. we note that the variance is larger for the full ml approach than the eb approach. this is further examined in additional file 3: figure s where we have plotted the difference in mean squared error for the same subregions with the ml method and the eb method. we notice that the mse for the eb method is smaller. effect of the prior for the eb method, the prior can have an impact on the estimates of tajima’s d. in figure  we explore the effect of the prior on the estimated values. we use three different priors. 1) a prior estimated on a neutral data set of 100 mb, 2) a prior estimated for a 100 mb region under selection 3) a prior based on both types of regions in equal proportion. as can be seen in figure  when applying the neutral prior to the neutral datasets or the selection prior on the selection dataset the estimates are approximately unbiased. however, applying the “wrong prior” tends to either underestimate or overestimate the selection signal. perhaps a bit surprisingly, using the neutral prior on the selection datasets tends to overestimate the signal, while using the selection prior on the neutral datasets tends to underestimate it. this is explained by the fact that a region of selection will have less variability than a neutral region. a prior trained on a neutral dataset will therefore have a higher level of variability than the true level of variability in a region of selection. applying a neutral prior will therefore allow more variability, and due to sequencing errors these sites will mostly be singletons. this excess of singletons will have the largest impact on θ^w, and we will therefore underestimate tajima’s d. the opposite argument explains why the selection prior will overestimate tajima’s d when applied to a neutral dataset. we observe the same trend for other neutrality tests statistics that can be seen in additional file 4: figure s and additional file 5: figure s however, the other investigated neutrality test statistics, fu & li’s d and f, are not estimated with the same accuracy as tajima’s d. for both selection and neutral data sets, we seem to underestimate the statistics with the eb approach and neutral prior. from figure we also see that we have problems estimating the true value in the region of the targeted locus using the 50 kb ml approach. the difference is perhaps caused by the fact that fu & li’s statistics are based on a single category of the frequency spectrum, whereas tajima’s d is based on all categories. power to detect a selective sweep often the goal of the investigator is not to estimate the neutrality test statistics in an unbiased manner but instead the goal is to identify regions with extreme values. to investigate our ability to discriminate between regions with selection and neutral regions, we show receiver operating characteristic for the different approaches. these are shown in additional file 6: figure s for different depth, error rates and number of individuals. when the depth is high, all methods perform almost as well as when the genotypes are known without error. but when the depth decreases or the error rates increases, only the full maximum likelihood approach and the eb approach have similar power as the known genotypes. as expected, the simulations show that our methods have improved power to discriminate between regions evolving neutrally and under positive selection as more samples are added. the genotype calling methods perform worse when the error rate is increased and the depth is decreased, especially at low depth with a low p-value cutoff, while the ml method for all scenarios performs almost as well as if the true genotypes where known. variable depth and variable error rates from ngs data so far, we have used a simple simulation model assuming sequencing depths follow a poisson distribution and a fixed error rate. to examine the robustness of our conclusions to these assumptions, we made an additional set of simulations using the observed distribution of quality scores and sequencing depths tabulated from bam files from the genomes project. as before we simulate 1 mb regions with and without selection for individuals, and apply our two genotype calling methods and the eb method to the simulated data. for the genotype calling approaches we only included sites that were likely to be polymorphic with a p-value <10- we observe results that are highly compatible with the previous results. the eb method is approximately unbiased in both scenarios and have similarly small mean squared error. in contrast, the genotype calling methods are biased in both scenarios, and we notice again that the bias depends on whether or not the data is simulated under a neutral or under a selection model. we observe that the mse for the genotype calling methods are orders of magnitude larger than the eb method, and that the mse for the gc-mlike model is more than twice as high under the selection scenario as under the neutral scenario, whereas we observe the opposite trend for gc-hwe method. application to real data we also applied our methods on real ngs data from the ceu individuals from the genomes project. bersaglieri et al., found a strong signal of positive selection surrounding the lct region of chromosome. we estimated tajima’s d values in a sliding window across this chromosome, using a prior estimated from a 50 mb region on chromosome. to compare with earlier published results we used the tajima track from the ucsc genome browser, which depicts estimated values of tajima’s d from individuals of european descent from the perlegen genotyping dataset. for our eb method we performed sliding windows analysis with different window sizes all using a fixed step size of 10 kb. we also compared the results for the naïve genotype calling methods using different snp inclusion cutoff criteria’s. notice that the overall estimate of tajima’s d is very positive for the snp data, most likely due to ascertainment biases. also notice that the lowest observed value of tajima’s d is the lct region for the eb approach while there are multiple regions with low tajima’s d values for the gc approaches and the snp chip data. discussion we have developed two methods to perform the neutrality tests on ngs data that take the uncertainty of the genotypes into account. both methods perform better than using called genotypes and in most instances they are approximately unbiased. the full ml method is computationally slow when applied in sliding windows at a genome-wide scale, which is why we also present a fast empirical bayes method with a prior that is estimated from the data itself, for example the entire genome, or a reasonable subset of the genome. this makes the method computationally feasible for full genomic data of any magnitude and any windows size. there is not a single obvious way to identify snp sites and call genotypes from ngs data. here we have tried different approaches with different cutoffs. regardless of method and chosen cutoff they all show a large bias in some or all simulated scenarios. this is evident from the raw theta estimates, and the actual test statistics. the level of bias varies between the different scenarios, not only for different depths and error rates, but it also depends on whether or not the data set is neutral or affected by selection. the results from this study suggest/confirms it is not unproblematic to perform neutrality tests on genotypes called from low or medium coverage ngs data. in contrast, both the ml and the eb approach give fairly accurate estimates for all the examined measures. when applying a neutral genome-wide prior for our analysis, we observed only small deviations from the true values of tajima’s d even for very low depth data. when applying the eb approach we did observe a small bias in the regions under selection when the prior was estimated from regions without selection. however, the bias is always much smaller than the bias of the gc approach. even though the eb approach can give small biases it can still have an advantage over the full ml approach. when estimating the neutrality test statistics for small windows, we often obtained a few extreme outliers with positive values of tajima’s d for the ml approach. since the eb method uses the entire 1 mb region as prior we do not see a similar problem with extreme positive outliers and the variance of the estimates is smaller overall. if the goal is to identify regions under selection the smaller variance of the eb approach will give fewer regions with extreme values. this is because regions with little data will increase the variance for the full ml approach while they will give results closer to the prior for the eb approach. this problem could also be circumvented by using a sliding window approach with window sizes determined by using a fixed number of snps. we applied the eb method to data from individuals from the genomes project, and observed trends similar to previous published results for the lct region in europeans. for snp chip data from the same region, we observed a large over-representation of positive tajima’s d values, presumably due to ascertainment biases introduced by the snp-discovery and selection process from the snp-chip itself. similar positive values were not observed using the eb method based on the ngs data. for both the gc and eb methods, we observe negative values around the lct gene, however the estimates are not very extreme for the gc approach. the eb approach is the only approach for which lct has the most extreme tajima’s d value. the main advantage of the approaches presented here, is that they, in expectation, have at most a very weak dependence on sequencing depth. this facilitates the use of genome-wide scans on genomes with varying sequencing depth without introducing biases. the computational framework suggested here, based on the eb approach, provides a robust and computationally fast method for scanning a genome for regions with outlying or extreme frequency spectrum. such a method should be of great use in years to come when analyzing population genomic data from a variety of different species.  in this paper we show through simulations that estimating neutrality test statistics using called genotypes can lead to highly biased result. the bias is dependent not only on sequencing depth and error rate but interestingly the bias also depends on whether or not the region is under selection or not. we have developed an empirical bayes method that can calculate the test statistic fast and efficiently. this method circumvents the problem of snp discovery, genotype calling and missing data, which is a fundamental problem of ngs data. this is done by working with genotype likelihoods, which contains all relevant information about the uncertainty of the data. using this approach leads to approximately unbiased estimates in most instances. availability all methods discussed in this paper are freely available as part of the analyses of next generation sequencing data package. competing interest the authors declare that they have no competing interests. authors’ contributions rn designed the eb model together with tsk and im. aa helped with the design of the software package and bug checked early versions of the program. tsk implemented the methods, carried out all analyses and simulations, and drafted the first version of the manuscript with rn. the manuscript has been thoroughly edited by the remaining authors. rn supervised the process. all authors read and approved the final manuscript. supplementary material additional file 1: figure s1 the effect of genotype calling for low or medium coverage data using fu & li’s d. the difference between estimated and known fu&li’s d statistic for three different scenarios with different p-value cutoffs. each box is estimated on the basis of 1 mb regions. the top figure is based on genotypes called using the frequency as prior, and the bottom figure is based on genotypes called using a maximum likelihood approach. notice that no single best cutoff can be chosen across the three different scenarios. click here for file additional file 2: figure s2 the effect of genotype calling for low or medium coverage data using fu & li’s f. the difference between estimated and known fu&li’s f statistic for three different scenarios with different p-value cutoffs. each box is estimated on the basis of 1 mb regions. the top figure is based on genotypes called using the frequency as prior, and the bottom figure is based on genotypes called using a maximum likelihood approach. notice that no single best cutoff can be chosen across the three different scenarios. click here for file additional file 3: figure s3 difference in mean squared error between the full ml and the eb method under a selective sweep. mse of the estimated tajima’s d is calculated for every 50 kb sub region of the full 1 mb region. the figure is based on 1 mb regions. for the eb method we used a prior estimated from the entire 1 mb region on every 50 kb subregion. click here for file additional file 4: figure s4 effect of different priors for the eb method using the fu & li’s d. left and center plot are boxplots for the difference between our estimate of fu & li d statistics and the true value, these are based on 100 × 1 mb regions. right plot is a 50 kb window plot using the 50 kb ml method along with the eb with neutral and mixed prior. neutral prior is from a genome-wide prior based on a 100 mb region, neu + sel prior is based on a 200 mb prior based on 100 mb selection and 100 mb neutral. click here for file additional file 5: figure s5 effect of different priors for the eb method using the fu & li’s f. left and center plot are boxplots for the difference between our estimate of fu & li f statistics and the true value. right plot is a 50 kb window plot using the 50 kb method along with the neutral and mixed prior. neutral prior is from a genome-wide prior based on a 100 mb region, neu + sel prior is based on a 200 mb prior based on 100 mb selection and 100 mb neutral. click here for file additional file 6: figure s6 power to detect a selective sweep. roc curve for scenarios, each plot is based on tajima’s d estimate for 100 × 1 mb regions with selection and 1 mb regions without selection. for each scenario we have our eb method along with our two genotype calling methods. row is individuals row is individuals. column is 8× 1% error rate, column2- is all 4×, but with varying error rates 1%, % and %. click here for file additional file 7: figure s7 roc curve for low coverage dataset. roc curve for a 2× % error rate. the lrt criteria is 10- this plot is based on 1 mb regions with selection, and 1 mb neutral regions. click here for file additional file 8: figure s8 distribution of quality scores and sequencing depth for a bam file. the left panel shows the quality score distribution and right panel shows the depth distribution, tabulated for chr from a bam file from the genomes project. the mean quality score value was approximately which corresponds to an average error rate of %. the data covered approximately 87% of the genome, had an average sequencing depth of, and contained sites with a sequencing depth above the right panel only shows the first observations. click here for file additional file 9: figure s9 using observed qscore and depth distributions. boxplots of the difference between the estimate of tajima’s d and the known value for 1 mb regions with our eb method and the two genotype calling methods. the quality score and depth distributions used for the genotype likelihood calculations are based on the results depicted in figure s for the genotype calling methods we used a cut-off for the p-value of the lrt test of 10- click here for file acknowledgements work for this manuscript was financed through danish national research foundation, villum foundation and danish council for independent research.
we have developed an approach that accommodates the uncertainty of the data when calculating site frequency based neutrality test statistics. using an empirical bayes approach for fast computations, we show that this method produces results for low-coverage ngs data comparable to those achieved when the genotypes are known without uncertainty. the method is implemented in a fast framework which enables researchers to perform these neutrality tests on a genome-wide scale. a number of different statistics are used for detecting natural selection using dna sequencing data, including statistics that are summaries of the frequency spectrum, such as tajima’s d. these statistics are now often being applied in the analysis of next generation sequencing data. a salient feature of this approach is that it implicitly solves the problems of varying sequencing depth, missing data and avoids the need to infer variable sites for the analysis and thereby avoids ascertainment problems introduced by a snp discovery process. however, estimates of frequency spectra from ngs data are strongly affected by low sequencing coverage; the inherent technology dependent variation in sequencing depth causes systematic differences in the value of the statistic among genomic regions.
 a number of different statistics are used for detecting natural selection using dna sequencing data, including statistics that are summaries of the frequency spectrum, such as tajima’s d. these statistics are now often being applied in the analysis of next generation sequencing data. however, estimates of frequency spectra from ngs data are strongly affected by low sequencing coverage; the inherent technology dependent variation in sequencing depth causes systematic differences in the value of the statistic among genomic regions.  we have developed an approach that accommodates the uncertainty of the data when calculating site frequency based neutrality test statistics. a salient feature of this approach is that it implicitly solves the problems of varying sequencing depth, missing data and avoids the need to infer variable sites for the analysis and thereby avoids ascertainment problems introduced by a snp discovery process.  using an empirical bayes approach for fast computations, we show that this method produces results for low-coverage ngs data comparable to those achieved when the genotypes are known without uncertainty. we also validate the method in an analysis of data from the genomes project. the method is implemented in a fast framework which enables researchers to perform these neutrality tests on a genome-wide scale. next-generation sequencingdarwinian selectionneutrality tests 