Hostname: page-component-89b8bd64d-rbxfs Total loading time: 0 Render date: 2026-05-07T02:32:15.409Z Has data issue: false hasContentIssue false
  • English
  • Français

Population genetic inference using a fixed number of segregating sites: a reassessment

Published online by Cambridge University Press:  21 January 2008

SEBASTIÁN E. RAMOS-ONSINS ,
SYLVAIN MOUSSET ,
THOMAS MITCHELL-OLDS  and
WOLFGANG STEPHAN
SEBASTIÁN E. RAMOS-ONSINS*
Affiliation:
Max Planck Institute of Chemical Ecology, Hans-Knöll-Strasse 8, 07745 Jena, Germany
SYLVAIN MOUSSET
Affiliation:
Biocenter, Department of Biology II, University of Munich, 82152 Planegg-Martinsried, Germany
THOMAS MITCHELL-OLDS
Affiliation:
Max Planck Institute of Chemical Ecology, Hans-Knöll-Strasse 8, 07745 Jena, Germany
WOLFGANG STEPHAN
Affiliation:
Biocenter, Department of Biology II, University of Munich, 82152 Planegg-Martinsried, Germany
*
*Corresponding author. e-mail: sramosonsins@ub.edu
Rights & Permissions [Opens in a new window]

Summary

Coalescent theory is commonly used to perform population genetic inference at the nucleotide level. Here, we examine the procedure that fixes the number of segregating sites (henceforth the FS procedure). In this approach a fixed number of segregating sites (S) are placed on a coalescent tree (independently of the total and internode lengths of the tree). Thus, although widely used, the FS procedure does not strictly follow the assumptions of coalescent theory and must be considered an approximation of (i) the standard procedure that uses a fixed population mutation parameter θ, and (ii) procedures that condition on the number of segregating sites. We study the differences in the false positive rate for nine statistics by comparing the FS procedure with the procedures (i) and (ii), using several evolutionary models with single-locus and multilocus data. Our results indicate that for single-locus data the FS procedure is accurate for the equilibrium neutral model, but problems arise under the alternative models studied; furthermore, for multilocus data, the FS procedure becomes inaccurate even for the standard neutral model. Therefore, we recommend a procedure that fixes the θ value (or alternatively, procedures that condition on S and take into account the uncertainty of θ) for analysing evolutionary models with multilocus data. With single-locus data, the FS procedure should not be employed for models other than the standard neutral model.

Information

Type
Research Article
Information
Genetics Research , Volume 89 , Issue 4 , August 2007 , pp. 231 - 244
Copyright
Copyright © Cambridge University Press 2007

1. Introduction

Monte Carlo simulation based on the coalescent is widely used in population genetics. This approach enables researchers to generate data for a given sample size under a panmictic neutral model or other evolutionary models. Using tests of neutrality or other summary statistics, the observed and simulated data can then easily be compared.

In order to simulate a sample under a strict neutral panmictic model, it is necessary to know the population mutation parameter θ (where θ=4Nμ, N being the effective population size and μ the mutation rate). The population mutation parameter, however, is generally unknown and must be estimated (Watterson, Reference Watterson1975; Tajima, Reference Tajima1983; Fu, Reference Fu1994; Griffiths & Tavaré, Reference Griffiths and Tavaré1994; Kuhner et al., Reference Kuhner, Yamato and Felsentein1995). To avoid the uncertainty in using an estimate of θ (usual estimation methods are mostly inefficient, e.g. mean pairwise difference is inconsistent, Watterson's estimate converges only asymptotically), several approaches have been proposed (e.g. see Hudson, Reference Hudson, Takahata and Clark1993; Markovtsova et al., Reference Markovtsova, Marjoram and Tavaré2000; Simonsen et al., Reference Simonsen, Churchill and Aquadro1995). The most popular approximate procedure simulates samples by fixing the number of segregating sites (S) instead of using the mutational parameter θ (Hudson, Reference Hudson, Takahata and Clark1993). However, using this approach, population samples with short genealogical trees tend to exhibit high mutation rates, and conversely, long genealogical trees have low mutational rates.

This latter procedure should not be confused with approaches that condition on S taking into account the parameter θ. A procedure that conditions on S (instead of fixing S) considers for a given θ all possible trees weighted by their probabilities of giving rise to the observed number of segregating sites (Tavaré et al., Reference Tavaré, Balding, Griffiths and Donnelly1997; Markovtsova et al., Reference Markovtsova, Marjoram and Tavaré2000, Reference Markovtsova, Marjoram and Tavaré2001; Jakobsson et al., Reference Jakobsson, Hagenblad, Tavaré, Säll, Halldén, Lind-Halldén and Nordborg2006), and thus takes into account the uncertainty of θ (see Kelly, Reference Kelly1997; Depaulis et al., Reference Depaulis, Mousset and Veuille2001, equation 2). In the case of conditioning on S, the internode times are not independent. Therefore, the procedure proposed by Hudson (Reference Hudson, Takahata and Clark1993) might be considered from two different points of view: (i) as an approximation of the standard procedure that is based on a (known) θ value, and (ii) as an approximation of a rigorous procedure that conditions on S and takes into account the uncertainty of θ.

Although fixing the number of segregating sites does not strictly follow the assumptions of coalescent theory (e.g. independence of the genealogical and mutational phases: see Kingman, Reference Kingman1982a, Reference Kingmanb; Hudson, Reference Hudson, Futuyma and Antonovics1990; Donnelly & Tavaré, Reference Donnelly and Tavaré1995; Nordborg, Reference Nordborg, Balding, Bishop and Cannings2001), this approximate procedure is reasonably accurate for obtaining critical values of statistics for the standard neutral model, provided that the true θ value is well supported by the data (Depaulis et al., Reference Depaulis, Mousset and Veuille2001; Markovtsova et al., Reference Markovtsova, Marjoram and Tavaré2001; Wall & Hudson, Reference Wall and Hudson2001). Thus, the probability of obtaining a θ value that is not supported by the data is expected to be low (Depaulis et al., Reference Depaulis, Mousset and Veuille2001), and the critical values of tests obtained by fixing the number of segregating sites appear to be quite accurate. However, for alternative evolutionary models the accuracy of this approximation is unknown, although many publications have used a fixed number of segregating sites for testing purposes (e.g. Braverman et al., Reference Braverman, Hudson, Kaplan, Langley and Stephan1995; Depaulis & Veuille, Reference Depaulis and Veuille1998; Wall, Reference Wall1999; Fay & Wu, Reference Fay and Wu2000; Przeworski, Reference Przeworski2002; Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002; Glinka et al., Reference Glinka, Ometto, Mousset, Stephan and De Lorenzo2003). It is therefore necessary to study the accuracy of this procedure for a wide range of models.

2. Simulation methods

Simulations have been performed with the program mlcoalsim (Ramos-Onsins & Mitchell-Olds, Reference Ramos-Onsins and Mitchell-Olds2007), which is available at http://www.ub.edu/softevol/mlcoalsim.

(i) Evolutionary models

The following models are considered:

  1. (a) The neutral panmictic model with constant population size.

  2. (b) The symmetric finite-island model (Wright, Reference Wright1943) with a constant number of populations (d), equal population size for each population, and a symmetric and constant migration parameter M among islands (with M=4Nm, where N is the effective population size of a deme and m the migration rate).

  3. (c) The logistic growth model (Fu, Reference Fu1997). In this the population size changes with time as follows:

    (1)
    \left. {\matrix{ {N\lpar t\rpar \equals N_{\setnum{0}} } \hfill \tab {{\rm if\ }t\les t_{\setnum{0}} \comma } \hfill \cr {N\lpar t\rpar \equals N_{\setnum{0}} \plus {\displaystyle {{N_{\setnum{1}} \minus N_{\setnum{0}} } \over {1 \plus e^{ \minus \gamma \left( t \minus t_{\setnum{0}} \minus \left( {{t_{\setnum{1}} \minus t_{\setnum{0}} } \over \setnum{2}}\right) \right) }} }}} \hfill \tab {{\rm if\ }t_{\setnum{0}} \lt t \lt t_{\setnum{1} }\comma } \hfill \cr {N\lpar t\rpar \equals N_{\setnum{1}} } \hfill \tab {{\rm if\ }t\ges t_{\setnum{1}} \comma } \hfill \cr} } \right\}

Here N(t) is the population size at time t (expressed in N 0 generations), N 0 is the population size at time t 0, and N 1 is the population size at time t 1. We used γ=10/(t 1−t 0). Coalescence times were calculated by integrating (1) over t. Expansion and reduction of population size was studied using this model. We used 10- and 100-fold differences between N 1 and N 0. We also used recent expansion processes, where t 1−t 0=0·1 N 0 generations, and t 0=0.

  1. (d) A bottleneck model using a constant population size initially for a period of time t d (expressed in 4N generations), then a sharp reduction for a time t b followed by a size increase. Reduction and expansion of population size also follow a logistic model. The conditions used here are bottlenecks of 10- or 100-fold differences in effective population size.

  2. (e) Hitchhiking model: We followed essentially the algorithm described in Braverman et al. (Reference Braverman, Hudson, Kaplan, Langley and Stephan1995) to generate hitchhiking genealogies, and allowed recombination within the locus of interest during the selective and neutral phases (see also Fay & Wu, Reference Fay and Wu2000; Kim & Stephan, Reference Kim and Stephan2002; Przeworski, Reference Przeworski2002). As in Fay & Wu (Reference Fay and Wu2000), we used as parameters the selection coefficient (s), the recombination rate between the selected locus and the studied locus (c), the intragenic recombination rate (r), and the time at which an advantageous mutation is fixed (t f). For the selective phase, we calculated the time to the most recent common ancestor (instead of checking at small time increments) for both the selected and unselected ‘subpopulations’ using the reasoning of Nordborg (Reference Nordborg, Balding, Bishop and Cannings2001, equation 7). The selective phase starts at time t f with a frequency of the selected allele of 1−1/2N, and ends when the frequency of x(t)<1/N. The value x(t) was calculated deterministically using equation 1 in Kim & Stephan (Reference Kim and Stephan2002) (see also equation 3a in Stephan et al., Reference Stephan, Wiehe and Lenz1992). The computer code was tested by comparing the results with those of table 2 in Stephan et al. (Reference Stephan, Wiehe and Lenz1992) and also by comparison with the ssw program (Kim & Stephan, Reference Kim and Stephan2002). The parameter values used in this work were t f=0, 4Ns=2×104, N=106 and ε=1/4N.

(ii) Monte Carlo methods

In all procedures, we used sample data (n lines) obtained from a diploid species. The procedures are as follows:

  1. (a) The Fθ procedure (Fixed θ procedure): This is the original standard coalescent procedure. We used a fixed value of the population mutation parameter θ. We placed a number of mutations on the tree according to the Poisson distribution with the mean value θ times the total length of the tree.

  2. (b) The FS procedure (Fixed S procedure): We placed a fixed number of segregating sites uniformly on each tree generated under the models mentioned above (Hudson, Reference Hudson, Takahata and Clark1993). The procedure for generating trees is identical to the Fθ procedure.

  3. (c) Procedure based on fixing the number of segregating sites and the value of θ (FSθ): This procedure employs the rejection algorithm #2 justified and described by Tavaré et al. (Reference Tavaré, Balding, Griffiths and Donnelly1997).

  4. (d) Procedure based on fixing the number of segregating sites but taking into account the uncertainty of the value of θ (FSθprior procedure): We use the RAU procedure (Rejection Algorithm using a Uniform prior). This procedure employs the rejection algorithm #2 described by Tavaré et al. (Reference Tavaré, Balding, Griffiths and Donnelly1997) but sampling of θ is done from a given prior distribution instead of having a fixed θ value.

(iii) Mutational parameter θ

We assume a uniform distribution over some arbitrarily chosen interval [θmin; θmax] (Depaulis et al., Reference Depaulis, Mousset and Veuille2001). Thus the prior density of θ is

(2)
g_{u} \lpar \theta \rpar \equals {1 \over {\theta _{\max } \minus \theta _{\min } }}.

When we do not have information about the value of θ, the assumption of a uniform density for any value of θ is reasonable. This is a commonly used strategy to estimate θ given observed data (e.g. Watterson, Reference Watterson1975; Wright & Charlesworth, Reference Wright and Charlesworth2004; Wright et al., Reference Wright, Vroh Bi, Schroeder, Yamasaki, Doebley, McMullen and Gaut2005; Haddrill et al., Reference Haddrill, Thornton, Charlesworth and Andolfatto2005). If the researcher has available information about the distribution of θ, then this information should be used instead of assuming a uniform distribution (e.g. Pritchard et al., Reference Pritchard, Seielstad, Perez-Lezaun and Feldman1999; Przeworski, Reference Przeworski2003), but we do not consider other prior distributions in this paper.

In order to avoid biologically unrealistic values of θ, we used as a minimum bound of θ per nucleotide a value of 0·0005, and as a maximum bound a value of 0·05. Numerous publications show that these numbers are realistic. Different bounds might modify some of our results, but do not change the main conclusions.

(iv) Statistical methods

From the simulated trees we calculated Tajima's D (Tajima, Reference Tajima1989), here named TD, Fu & Li's D and F (Fu & Li, Reference Fu and Li1993), here named FD and FF, respectively, and Fay & Wu's H (Fay & Wu, Reference Fay and Wu2000), abbreviated as H, the statistic B (Wall, Reference Wall1999) (considered sensitive to structured populations), F S (Fu, Reference Fu1997) and R2 (Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002), which are sensitive to population size expansion, as well as the number of haplotypes (here divided by the sample size) K w (Strobeck, Reference Strobeck1987; Fu, Reference Fu1996; Depaulis et al., Reference Depaulis, Mousset and Veuille2001; Wall, Reference Wall1999) and the haplotype diversity H w (Depaulis & Veuille, Reference Depaulis and Veuille1998). In total, nine summary statistics were computed. The calculation of these statistics was examined with the software package DnaSP 3.51 (Rozas & Rozas, Reference Rozas and Rozas1999; Rozas et al., Reference Rozas, Sanchez-DelBarrio, Messeguer and Rozas2003).

For multilocus analyses, we treated each locus independently because we assumed that the studied loci are unlinked. When we used the RAU procedure, we chose independent θ values for each locus in order to perform simulations. For a given statistic, we recorded the P value independently for each locus and combined them as in Voight et al. (Reference Voight, Adams, Frisse, Qian, Hudson and Di Rienzo2005) but calculated each tail separately.

To avoid excessively liberal critical values (95% interval of the null sampling distribution) for discrete distributions, we have used the following procedure: for the 2·5% interval of the upper tail, we took the first value that is larger than the observed value at 2·5% of the distribution (see Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002). The same logic was applied to the 2·5% interval of the lower tail and also for comparing discrete distributions. Because of this issue, the realized type I error rates are slightly lower than the expected 5% and have to be assessed for every statistic and method.

(v) Determining the accuracy of the FS procedure in relation to Fθ

We have used the approach described in Wall & Hudson (Reference Wall and Hudson2001) to determine the level of accuracy of the FS procedure compared with the Fθ procedure. In this approach, a number of simulations using a large range of values of S (e.g. from S=1 to S=120) is obtained by the FS procedure, and its critical values and the false positive rate (called the ‘nominal’ size) for each statistic are calculated for a given evolutionary model. Then, the ‘true’ θ value for the same evolutionary model is obtained, which is the value that gives an average of S=20 in 10 000 iterations under the Fθ procedure. Next, a simulation with the Fθ procedure (i.e. the null hypothesis) is performed using the ‘true’ fixed value of θ. For each iteration, the value of S is calculated, and acceptance or rejection of the null hypothesis is determined for every statistic by comparing its value with the critical values of the corresponding FS distribution. The rejection rate (i.e. the false positive rate given the FS critical values, here called the ‘true’ size) for each statistic is stored. The choice of the ‘true’ θ value, although somewhat arbitrary, is suggested by the design of the study, which fixes the number of segregating sites instead of fixing the value of θ (which is unknown to the researcher). Using the following method, we try to be conservative in the sense of minimizing the differences between the two different procedures (see Wall & Hudson, Reference Wall and Hudson2001; Depaulis et al., Reference Depaulis, Mousset and Veuille2003, Reference Depaulis, Mousset, Veuille and Nurminsky2005).

(vi) Determining the accuracy of the FS procedure in relation to FSθU

Here we assume that the θ value is unknown. The FSθU procedure is considered the ‘true’ procedure because it takes into account the uncertainty assumed by the researcher. In this approach, a simulation for a given evolutionary model is performed for the FS procedure (for a fixed S=20) and its critical values and the false positive rate (called the ‘nominal size’) of each statistic are stored. Then, a simulation using the FSθU procedure (for S=20) is performed for the same evolutionary model (the null hypothesis). The ‘true’ false positive rate is calculated in the following way: for each iteration calculated with the FSθU procedure, rejection of the null hypothesis is determined for every statistic by comparing its value with the critical values of the FS distribution. Finally, the rejection rate for each statistic is stored.

(vii) Parameter values

We used n=20 for a locus of 1000 nucleotides. For procedures that fix S, we preferentially used S=20. The reason for choosing n=20 and S=20 is largely historical, as the value of θ (given the neutral equilibrium model) is then approximately 0·005 (which has been used in a number of the theoretical population genetics studies). All θ values under any evolutionary model studied here were in the range we considered biologically realistic. For n=20 and S=50 we found similar results.

3. Analytical approaches

(i) Effects of recombination

Recombination is difficult to take into account in analytical models. However, one can consider the extreme case of free recombination (e.g. a sequence consisting of m freely recombining fragments of equal size such as m nucleotides). This is the limiting case when the recombination parameter R=4Nr becomes very large. Let the parameter T be the vector of the coalescent times T k, then L_{n} \equals \sum\nolimits_{k \equals \setnum{2}}^{n} {kT_{k} } is the total length of the coalescent tree (by summing the length of all branches) measured in units of 4N generations. Assuming m is large, the central limit theorem shows that the prior distribution of the average length \overline{L_{n} } of the m independent coalescent trees will converge towards the normal distribution with mean a 1(n) and variance a 2(n)/m (defined in equations (A3) and (A4) of Appendix 1, respectively). This is,

(3)
\mathop {\lim }\limits_{R \to \infty } f_{\overline{{L_{n}}}}} } \lpar t\rpar \equals \sqrt {{m \over {{\setnum{2}}\pi a_{\setnum{2}} \lpar n\rpar }}} e^{ \minus \displaystyle{ {m \over {\setnum{2}a_{\setnum{2}} \lpar n\rpar }}\lpar t \minus a_{\setnum{1}} \lpar n\rpar \rpar ^{\setnum{2}}} } \comma

where f\hskip2_\overline{L_{_{n}} } \lpar t\rpar is the density of the average length \overline{L_{n} }. This last equation can be used in equations (A7)–(A8) (Appendix 1) to address the effects of the simulation procedures Fθ, FS, FSθ, and FSθprior on the posterior density of the length of the coalescent tree in the limiting case of freely recombining sequences.

(ii) Shape of the coalescent trees

In the present work, we focus on procedures where simulations use a fixed number of segregating sites S. In these procedures, unlike in the Fθ procedure, the total length of a tree does not play a role as long as the shape of the tree remains the same (trees are only scaled). In order to assess the impact of the rejection algorithm on the shape of the trees, we study the ratio of the branch lengths in the upper and lower parts of the tree X/Y n (Fig. 1). The mean value of this ratio is given by

(4)
\eqalign {E\left( {X \over {Y_{n} }}\vert {n\comma S} \right) \equals\tab {1 \over {\eta \lpar n\comma S\rpar }}\cr \tab\hskip-50\times \int\nolimits_{\setnum{0}}^{\infty } {\int\nolimits_{\setnum{0}}^{\infty } {{x \over y}f_{X} \lpar x\rpar f_{Y_{n} } \lpar y\rpar {\bf P}_{{\bf accept}} \lpar x\comma y\rpar } {\rm d}x{\rm d}y}\comma}\hskip-12pt

where f X and f_{Y_{n} } denote the prior densities of X and Y n (see equations (A5) and (A6)) and

\eqalign{{\bf P}_{{\bf accept}\lpar {\bf x}\comma {\bf y}\rpar } \equals\hskip186 \cr\left\{ {\matrix{ 1 \hfill \tab {{\rm \lpar }FS{\rm \ procedure\rpar }} \hfill \cr {\displaystyle{{\theta ^{S} \lpar x \plus y\rpar ^{S} } \over {S\exclam }}e^{ \minus \lpar x \plus y\rpar \theta } } \hfill \tab {{\rm \lpar }FS\theta {\rm \ procedure\rpar }} \hfill \cr {\displaystyle{{\lpar x \plus y\rpar ^{S} } \over {S\exclam }}\int_{\setnum{0}}^{\infty } {e^{ \minus \lpar x \plus y\rpar \theta } g\lpar \theta \rpar {\rm d}\theta } } \hfill \tab {{\rm \lpar }FS\theta _{prior} {\rm \ procedure\rpar }} \hfill \hskip7pt\cr} } \right.}

and \eta \lpar n\comma S\rpar \equals \int_{\setnum{0}}^{\infty } {\int_{\setnum{0}}^{\infty } {f_{X} \lpar x\rpar f_{Y_{n} } \lpar y\rpar {\bf P}_{{\bf accept}} \lpar x\comma y\rpar {\rm d}x{\rm d}y} } is a normalizing constant.

Fig. 1.

Shape index of a coalescent tree. The shape index is the ratio of the branch lengths in the upper and lower parts of the tree, X/Y n (see text).

4. Simulation results

We determined the accuracy of the approximate procedure proposed by Hudson (Reference Hudson, Takahata and Clark1993) (which fixes the number of segregating sites while ignoring the value of the population mutation value θ; the FS procedure) in comparison with the two rigorous procedures: (i) the standard procedure fixing θ (Fθ) and (ii) procedures that while conditioning on S take into account the uncertainty of θ (FSθprior, here using a uniform prior, named FSθU). We consider several different evolutionary models including the neutral panmictic model, population subdivision, population size bottleneck, expansion and genetic hitchhiking. As a measure of accuracy, we calculated the difference in the false positive rate between procedures.

(i) Comparison of the FS procedure with the standard procedure Fθ

We examined the accuracy of the FS procedure for different alternative models. To do this we also re-examined the type I error for the standard neutral model (neutral panmictic population with zero recombination) for the nine statistics studied here, since the type I error for single-locus data was also studied in Depaulis et al. (Reference Depaulis, Mousset and Veuille2001), Markovtsova et al. (Reference Markovtsova, Marjoram and Tavaré2001) and Wall & Hudson (Reference Wall and Hudson2001).

Table 1 shows the difference between the nominal size (the false positive rate of a test given the FS procedure) and the ‘true’ size (the false positive rate of a test given the Fθ procedure) for a 5% critical region. Our analysis confirmed the small difference in the type I error rate observed by Wall & Hudson (Reference Wall and Hudson2001) and Depaulis et al. (Reference Depaulis, Mousset and Veuille2001) under a neutral model for a single locus, although we used some other statistics and slightly different conditions (n=20 and θ=0·0057 for the Fθ procedure; i.e. the estimate of θ for S=20). Fig. 2 A shows the differences in size for each S value among the procedures for the nine studied statistics under the neutral model. Negative values indicate that the critical values of the FS procedure are more liberal than expected (increased type I error). The large fluctuations observed in Fig. 2 are also a consequence of the discrete distribution of values obtained when the number of segregating sites is fixed (see Ramos-Onsins & Rozas, Reference Ramos-Onsins and Rozas2002). Important differences are observed in the proportion of acceptance/rejection when the observed number of mutations is far from the average (S=20), as indicated in Wall & Hudson (Reference Wall and Hudson2001) (see also Depaulis et al., Reference Depaulis, Mousset and Veuille2003, Reference Depaulis, Mousset, Veuille and Nurminsky2005). Nonetheless, the more extreme values contribute very little such that when summed the differences cancel each other out, thus leading to a good accuracy of the FS procedure for the standard neutral model (Table 1). Therefore, we consider the FS procedure as sufficiently accurate under the neutral panmictic model for single-locus data.

Fig. 2.

Differences between the sizes of FS and Fθ procedures for each S separately. (A) Neutral model: In the upper panel the distribution of S values obtained with the Fθ procedure using θ=0·0057 is shown. In the middle and lower panels the size differences for statistics for the upper (97·5%) and the lower (2·5%) tails, respectively, are presented. Note that large S values, although causing large differences, contribute very little to the total. (B) Hitchhiking model: In the upper panel the distribution of S values obtained with the Fθ procedure using θ=0·0127 is shown. In the middle and lower panels the size differences for statistics for the upper (90%) and the lower (10%) tails, respectively, are presented. Abbreviations are explained in Section 2. Critical values are not calculated for B (lower tail), and for K w, and H w (upper tail), because these statistics are not conservative with recombination.

Table 1.

Difference (in percentage) between the FS and Fθ procedures for several critical valuesFootnote a

a Abbreviations are explained in Simulation methods. The critical values studied are 2·5% and 97·5% for the neutral model and 10%, 50% and 90% for alternative models (see Simulation methods). The maximum precision error detected is around ±1%.

b Island model with d=10 subpopulations, 4Nm=0·5 and 4Nr=10. The samples are all obtained from a single population. In the FSθU procedure, the bounds of θ were reduced d times for a better comparison with the other models. In the Fθ procedure, θ=0·00062.

c Logistic expansion model with a 10-fold growth of population size. The expansion process started 0·4N 0 generations ago and finished at present. 4N 0r=10. In the Fθ procedure, θ=0·03384.

d Logistic contraction model: In contrast to the logistic expansion model, the population is reduced 10 fold. In the Fθ procedure, θ=0·00059.

e Bottleneck model: Population size is maintained constant during 0·5N 0 generations since present. It follows a logistic, 100-fold reduction of population size for 0·01N 0 generations, maintained for 0·01N 0 generations, then an instantaneous recovery to the present population size. 4N 0r=10. In the Fθ procedure, θ=0·00860.

f Hitchhiking model: The selective phase is completed at present time (tf=0). N=106, 4Ns=104, the population recombination parameter between the selected and observed locus is 4Nc=400, c/s=0·02, and intragenic recombination is indicated in the table. In the Fθ procedure, θ=0·00127.

Next we consider the accuracy for alternative models. Table 1 shows the difference between nominal (FS procedure) and true (Fθ procedure) distributions of statistical tests for certain critical values. The reason to use different critical values for alternative models (10%, 50% and 90% instead 2·5% and 97·5%) is because we are interested in whether the probability distribution for a given statistic can be accurately represented using the FS procedure. The analysis is performed for five different alternative models (subdivision, expansion, contraction, bottleneck and hitchhiking) and for some arbitrary parameters. The strong differences observed between the two procedures indicate that the FS procedure is inaccurate when alternative models are used. Fig. 2 B shows considerable differences between nominal and true sizes for several of the studied statistics for each S separately, leading to a large difference in total (Table 1).

(ii) Comparison of the FS procedure with a procedure that conditions on S and takes into account the uncertainty of θ

We have examined the accuracy of FS (the nominal procedure) by comparison with the procedure FSθU (the ‘true’ procedure, see Section 2). With regard to the type I error for the standard neutral model, we have studied the 95% interval of the nine neutrality test statistics for the FS procedure and the procedure FSθU by simulation for n=20 and a range of S values from S=2 to 80. Fig. 3 shows the difference between the nominal and ‘true’ sizes. Liberal tests are Fay & Wu's H with regard to the upper tail, and Fu's F S, Ramos-Onsins & Rozas' R2, and Fu & Li's D and F for the lower tail. In the case of Fu & Li's F, the test is, in fact, less liberal than shown, because the critical value for the FS procedure is sometimes much lower than 2·5% (given the discrete distribution of values; not shown). We observed that all statistics at intermediate S values show a difference lower than 1·5%. Therefore, in this comparison the FS procedure is sufficiently accurate for statistical inferences under the neutral panmictic model for single-locus data.

Fig. 3.

Differences between the sizes of the FS and FSθU procedures. Critical values are obtained for each of the nine neutrality tests for the 2·5% upper and lower tails. Abbreviations are indicated in Section 2.

In contrast, for alternative models there are strong differences between the nominal (FS procedure) and true (FSθU procedure) distributions of statistical tests for 10%, 50% or 90% critical values (Table 2). The important differences observed between the two procedures indicate that the FS procedure is also inaccurate when it is compared with the FSθU procedure and when alternative models are used.

Table 2.

Difference (in percentage) between the FS and the FSθU procedures for several critical valuesFootnote a

a As in Table 1.

We performed a multilocus analysis for the nine test statistics using the combined P values for different numbers of loci (see Section 2) for the standard neutral model. The difference between the nominal and true size is examined in Fig. 4. The results indicate that for a larger number of loci the FS procedure becomes inaccurate (through the accumulation of small departures in the single-locus tests). Thus, the small differences observed for one locus between FS and FSθU become very important when multilocus studies are performed. Some statistics, like K w, H w and Fay & Wu's H, are extremely conservative for the lower tail, as are most statistics (except for Fay & Wu's H) for the upper tail (see Fig. 4 A). On the other hand, statistics such as Tajima's D, Fu's F S, Fu & Li's D and F, and Ramos-Onsins & Rozas' R2 are too liberal for the lower tail of the distribution, and Fay & Wu's H is also too liberal for the upper tail of the distribution.

Fig. 4.

Effect of the number of loci on the probability of rejecting the neutral panmictic model for nine neutrality tests. (A) Differences between the sizes of FS and FSθU procedures with no recombination in the upper 2·5% tail, given different numbers of loci. (B) Differences between the sizes of FS and FSθU procedures with no recombination in the lower 2·5% tail, given different numbers of loci. Abbreviations are explained in Section 2. Critical values are not calculated for B (lower tail), and for K w, and H w (upper tail), because these statistics are not conservative with recombination. S was fixed at 20 for each locus. Plots obtained using S values from a distribution compatible with θ=0·0057 gave equivalent results (not shown) although for a small number of loci we observed a large variance.

5. Discussion

We have compared the FS procedure with the standard procedure using a fixed θ value (Fθ) and with a procedure that conditions on S taking into account the uncertainty of θ (FSθU). Our results show that the FS procedure is inaccurate in both comparisons when alternative models are used.

(i) Causes of the discrepancy between the FS and the Fθ or FSθU procedures

In the comparison of FS with Fθ, differences in the length (L n) and topology (X/Y n) of the trees for the FS and FSθ procedures can explain the inaccuracy of FS, given that the critical values of FS are compared with the FSθ distribution for each S. The distribution of the total length of the coalescent trees of a non-recombining sequence is shown in Fig. 5. The shapes of the FS distribution and the distribution obtained with the FSθ procedure (taking S=20) are clearly different. The distribution of L n obtained with the FSθ procedure has lower variance than the FS procedure, which can be explained by the lack of uncertainty in the value of θ (see also Table 3). The mean and standard deviation of the shape index X/Y n for the different simulation procedures are shown in Table 4. Higher values are obtained for trees with long internal branches whereas smaller ratios denote shorter internal branches. We observed that for S=20 the trees obtained with the FSθ procedure have smaller ratios X/Y n than the trees obtained with the FS procedure. These observations can be explained as follows. For a given θ, when S is lower than the average for that θ, this will often have been caused by a tree that is shorter than average. Although the internode times are a priori independent (because the total tree length is dominated by the last few internode times, given a short tree), it will often be that the final internode time is particularly short. Thus X/Y n is smaller than the unconditional expectation. The situation is reversed for S greater than the average for a given θ. Thus, these results show how the simulation procedure may affect the shape of the sampled coalescent trees and may lead to different outcomes of neutrality tests.

Fig. 5.

Posterior densities of the total length L n of a coalescent tree with n=20 are shown for the different procedures and no recombination. Results based on the FS, FSθ and FSθU procedures are displayed. For FSθ the θ parameter was arbitrarily set to 20/a 2(20).

Table 3.

Mean and standard deviation (in parentheses) of the posterior distribution of L 20

a Assuming m=50 recombining fragments (equation 3).

b Assuming θ=20/a1(20).

Table 4.

Mean and standard deviation (in parentheses) of the tree shape

a See Figure 1. Assuming no recombination and n=20.

b Assuming θ=20/a1(20).

To understand the causes of the inaccuracy of the FS procedure in comparison with the procedure FSθU, we examined differences in the trees by considering the total length of the tree (L n) analytically and by simulations. The distribution of the total length of the coalescent tree of a non-recombining sequence is shown in Fig. 5. Some inaccuracies of the FS procedure in comparison with the procedure FSθU are observed under a panmictic neutral model with no recombination. Differences in the distributions of L n are small but apparent, with longer and slightly broader distributions for the procedure FSθU than for FS. It is noteworthy that this procedure samples shorter trees with a lower variance of L n than the FS procedure (Table 3). When we analysed the shape of the coalescent trees (Table 4), the coalescent trees obtained by the FSθU procedure are skewed towards shorter trees with smaller internal branches than trees sampled with the FS procedure.

Figure 6 shows the L n distribution for four of the alternative models analysed with the FS and FSθU procedures. There are strong differences in the distribution of L n for the FS procedure in comparison with FSθU for most of the alternative models. In particular, in the case of subdivided populations, large differences are observed for the parameter values studied. Bottleneck and hitchhiking models show also apparent differences in the trees that are used in FS in relation to FSθU. On the other hand, models of population size expansion exhibit smaller differences among procedures. In summary, differences in the length (and in the topology) of the trees explain the inaccuracy of the FS procedure and these differences are large in most of the alternative models studied.

Fig. 6.

Probability distribution of L n for four alternative models using n=20 and S=20. The FS and FSθU procedures are displayed in each case. Parameter values are the same as in Table 1 but intragenic recombination was set to 4Nr=10 (except for panel D where it was zero). (A) Island model. (B) Logistic expansion model. (C) Bottleneck model. (D) Hitchhiking model.

(ii) The effect of recombination

Analysing the effect of recombination is also important, as the true value of the recombination parameter is generally unknown. Recombination has an important effect on the distribution of segregating sites because it breaks up the correlation among contiguous positions. Differences in L n are observed when recombination is added to the model (Fig. 7). Average values of the parameter L n are constant under the FS procedure, as expected. Indeed, the FS procedure uses all output trees, like the standard coalescent procedure uses a fixed θ value. The expected average value of L n for the FS procedure value is given by equation (A3) in Appendix 1. Therefore, for n=20, E(L n)=3·548. On the other hand, the procedure FSθU has lower average L n values for zero recombination (L n is 3·16 for n=20), and this value increases for larger recombination values. The procedure FSθ leads to a similar pattern as that shown for FSθU, but with fewer differences relative to FS.

Fig. 7.

Effect of recombination on average L n for the FS, FSθ (θ=0·0057) and FSθU procedures in a neutral panmictic population using n=20 and S=20.

The means and standard deviations of L n in the limiting case of free recombination between fragments are given in Table 3. Recombination tends to bring the posterior distributions closer to the prior distribution of L n. Whereas the difference between the FS and FSθU simulation procedures remains clear, the outcomes of the FS and FSθU simulation procedures can hardly be distinguished (independent of the assumed prior for θ), confirming the observation of Fig. 7.

The consequences of having differences in the average L n given no recombination might be important, because it indicates that the zero-recombination neutral panmictic model may have an average of a given statistic that is different from that of a model with recombination, and therefore may not be conservative. Zero recombination in a neutral panmictic model leads to the largest deviation in average L n in relation to the FS procedure. We have observed that for single loci this difference can be tolerated, but becomes too large for multilocus analyses. That is, for multilocus analyses, the FS procedure should not be used unless recombination is quite high for each locus.

The prior distribution of the mutational parameter θ is shown to have an important effect on the posterior distribution of coalescent trees (Fig. 5), although, as previously pointed out by Wall & Hudson (Reference Wall and Hudson2001), the recombination parameter also has a strong effect on this distribution. The improvement obtained from using sophisticated techniques such as the FSθU procedure may, however, be negligible in comparison with the errors caused by the uncertainty in the recombination parameter R (Wall & Hudson, Reference Wall and Hudson2001). A procedure considering also the uncertainty in the recombination parameter would give more appropriate distributions, but is beyond the scope of this paper.

6. Appendix 1

(i) The branch lengths in a coalescent tree

In the standard neutral coalescent process without recombination, the waiting time T n (in units of 4N generations) until two lineages among n have a common ancestor is exponentially distributed with parameter n (n−1) (Kingman, Reference Kingman1982b). The probability density of T n is thus

(A1)
f_{T_{n} } \lpar t\rpar \equals n\lpar n \minus 1\rpar e^{ \minus n\lpar n \minus \setnum{1}\rpar t}.

The density of the sum L n of the lengths of the branches of a coalescent tree with n tips is given by Tavaré (Reference Tavaré1984, unlabelled equation at the top of p. 153) as

(A2)
f_{L_{n} } \lpar t\rpar \equals \lpar n \minus 1\rpar e^{ \minus t} \lpar 1 \minus e^{ \minus t} \rpar ^{n \minus \setnum{2}}.

The mean and variance of L n are

(A3)
E\lpar L_{n} \rpar \equals a_{\setnum{1}} \lpar n\rpar \equals \mathop\sum\limits_{i \equals \setnum{1}}^{n \minus \setnum{1}} {{1 \over i}\comma }
(A4)
V\lpar L_{n} \rpar \equals a_{\setnum{2}} \lpar n\rpar \equals \mathop\sum\limits_{i \equals \setnum{1}}^{n \minus \setnum{1}} {{1 \over {i^{\hskip2\setnum{2}} }}.}

In order to characterize the shape of a coalescent tree, we introduce a shape index X/Y n where X is the length of the two branches from the second-last coalescent to the last (in the upper part of the coalescent tree) and Y n is the length of the rest (the branches in the lower part) of the coalescent tree (see Fig. 1). The density of X follows from equation (A1):

(A5)
f_{X} \lpar t\rpar \equals e^{ \minus t}.

As shown in Appendix 2, the density of Y n is

(A6)
f_{Y_{n} } \lpar t\rpar \equals \lpar n \minus 1\rpar \lpar n \minus 2\rpar e^{ \minus \setnum{2}t} \lpar 1 \minus e^{ \minus t} \rpar ^{n \minus \setnum{3}}.

(ii) Branch lengths in the procedures Fθ and FS

In the standard procedure Fθ, but also in the FS simulation procedures, all simulated coalescent trees are used and thus the posterior density of tree length is the same as the prior density f_{L_{n} } (equation A2).

(iii) Branch lengths in the FSθ and FSθprior procedures

In the FSθ simulation procedure, trees are sampled from the prior distribution according to the probability of observing exactly S mutations given a known mutational parameter θ. The number of mutations in the coalescent tree follows a Poisson distribution with parameter θL n. Thus the probability of observing S mutations in a sample of n lines given θ is (Tavaré, Reference Tavaré1984, unlabelled equation on p. 153)

{\bf P}\lpar S \equals k\vert n\comma \theta \rpar \equals \int_{\setnum{0}}^{\infty } {{{\lpar \theta t\rpar ^{k} e^{ \minus \theta t} } \over {k\exclam }}f_{L_{n} } \lpar t\rpar {\rm d}t} \comma

and the posterior density of L n given S and θ follows from the definition of the posterior density:

(A7)
f_{L_{n} \vert {S\comma \theta }} \lpar t\rpar \equals {{\lpar \theta t\rpar ^{k} e^{ \minus \theta t} } \over {k\exclam {\bf P}\lpar S \equals k\vert n\comma \theta \rpar }}f_{L_{n} } \lpar t\rpar.

In the FSθ procedures that consider the uncertainty of the value of θ (FSθprior), a prior distribution of θ with density g is assumed. Thus the probability of observing S mutations in a sample of n lines is

{\bf P}\lpar S \equals k\vert n\comma \, g\rpar \equals \int_{\setnum{0}}^{\infty } {\int_{\theta _{{\min }} }^{\theta _{{\max }} } {{{\lpar \theta t\rpar ^{k} e^{ \minus \theta t} } \over {k\exclam }}f_{L_{n} } \lpar t\rpar g\lpar \theta \rpar {\rm d}\theta {\rm d}t} \comma }

and the posterior density of L n given S and the prior density g is

(A8)
\eqalign{f_{L_{n} \vert {S\comma g} } \lpar t\rpar \equals \tab {{f_{L_{n} } \lpar t\rpar } \over {{\bf P}\lpar S \equals k \vert {n\comma g} \rpar }}\times \int_{\theta _{{\min }} }^{\theta _{{\max }} } {{{\lpar \theta t\rpar ^{k} e^{ \minus \theta t} } \over {k\exclam }}} g\lpar \theta \rpar {\rm d}\theta.}

If a uniform density g u is assumed for θ, g u is constant (equation 2) and the posterior density in equation (A8) becomes

f_{L_{n} \vert {S\comma g} } \lpar t\rpar \equals {{f_{L_{n} } \lpar t\rpar \int_{\theta _{{\min }} }^{\theta _{{\max }} } {{{\lpar \theta t\rpar ^{k} e^{ \minus \theta t} } \over {k\exclam }}} \, {\rm d}\theta } \over {\int_{\setnum{0}}^{\infty } {\int_{\theta _{{\min }} }^{\theta _{{\max }} } {{{\lpar \theta t\rpar ^{k} e^{ \minus \theta t} } \over {k\exclam }}} f_{L_{n} } \lpar t\rpar {\rm d}\theta {\rm d}t} }}.

Using θmin=0, we see that for n⩾3 a limit of the posterior distribution of L n exists when θmax is very large (this condition is necessary for the integral in the denominator of equation A9 to converge):

(A9)
\lim_{\theta_\max\to\infty}f_{L_{n} \vert S,g_u }\lpar t\rpar\equals {{t^{-1}f_{L_n}\lpar t\rpar}\over{\int_0^\infty t^{-1}f_{L_n}\lpar t \rpar \,{\rm d}t}}\,.

It is noteworthy that this distribution is independent of the observed number of mutations S and may be denoted as f_{L_{n} \vert g_{u} }.

7. Appendix 2. Proof of equation (A6)

(i) Notation and preliminary results

In this section the density function of a random variable X will be denoted as f X. We introduce a new random variable P k=kT k (see equation A1). It is straightforward to show that P k is exponentially distributed with parameter k−1, thus

(A10)
f_{P_{k} } \lpar t\rpar \equals \lpar k \minus 1\rpar e^{ \minus \lpar k \minus \setnum{1}\rpar t} \comma

and Y n is defined by Y_{n} \equals \sum\nolimits_{k \equals \setnum{3}}^{n} {P_{k} } (see Fig. 1). Y n is the sum of n−2 independent exponentially distributed random variables with parameters 2, …, (n−1). As we show below, the density function of Y n can be obtained using simple order statistics.

The following properties of the exponential distribution will be used:

  • • Minimum of independent exponentially distributed random variables: Consider k independent random variables X_{\lambda_{i}} exponentially distributed with parameters λi (i=1, …, k). Then {\rm min}\lpar X_{\lambda_{i}} \rpar is exponentially distributed with parameter \lambda \equals \sum\nolimits_{i \equals \setnum{1}}^{k} {\lambda _{i} }.

  • • Memoryless property of the exponential distribution: Consider a random variable X λ exponentially distributed with parameter λ:

    \hskip9\forall x\ges 0\comma \ {\bf P}\lpar X_{\lambda } \minus x\les t \vert {X_{\lambda } \ges x} \rpar \equals 1 \minus e^{ \minus \lambda t} \equals {\bf P}\lpar X_{\lambda } \les t\rpar.

(ii) Density function of Y n

We consider n−1 independent random variables E i(i=1, …, n−1), exponentially distributed with parameter 1, and denote the smallest one as E (1)⋅ E (1) is exponentially distributed with parameter n−1, and the n−2 remaining random variables are independent and such that E i−E (1) is exponentially distributed with parameter 1. Then the difference between the second and the first smallest random variables E (2)−E (1) is exponentially distributed with parameter n−2 and the n−3 remaining random variables are independent and such that E i−E (2) is exponentially distributed with parameter 1. If we define E (0)=0, it is straightforward to show that the difference between two successive sorted random variables (E k−E (k−1)) is exponentially distributed with parameter n−k. E (n−2) can be rewritten as E_{\lpar n \minus \setnum{2}\rpar } \equals \sum\nolimits_{k \equals \setnum{1}}^{n \minus \setnum{2}} {\lpar E_{\lpar k\rpar } \minus E{}_{\lpar k \minus \setnum{1}\rpar }\rpar}. This is the sum of n−2 independent exponentially distributed random variables with parameters 2, …, (n−1), thus E (n−2) and Y n have the same distribution. Because E (n−2) is the (n−2)th random variable among the n−1 random variables E i sorted in increasing order, the density function of E (n−2) and Y n is (Pitman, Reference Pitman1992, p. 326):

(A11)
f_{Y_{n} } \lpar t\rpar \equals \lpar n \minus 1\rpar \lpar n \minus 2\rpar e^{ \minus \setnum{2}t} \lpar 1 \minus e^{ \minus t} \rpar ^{n \minus \setnum{3}}.

We thank J. M. Comeron, J. Rozas, K. Schmid, J. Zavala and A. Navarro-Quezada for helpful comments on the manuscript. S. R. O. thanks S. Gebauer-Jung and S. Wright for help with the calculation of Hudson's recursive equation, J. Fay for supplying information to check the selective model and Y. Kim for his assistance in the analysis of the selective model. Special thanks go to two anonymous reviewers for their input into this manuscript. This research was supported by the Max Planck Society, a Marie Curie fellowship from the EU (MCFI 2002 01461) to S. M., and grant Ste 325/5 from the DFG to W. S.

References

Braverman, J. M., Hudson, R. R., Kaplan, N. L., Langley, C. H. & Stephan, W. (1995). The hitchhiking effect on the site frequency spectrum of DNA polymorphisms. Genetics 140, 783–796.Google Scholar
Depaulis, F. & Veuille, M. (1998). Neutrality tests based on the distribution of haplotypes under an infinite-site model. Molecular Biology and Evolution 15, 1788–1790.Google Scholar
Depaulis, F., Mousset, S. & Veuille, M. (2001). Haplotype tests using coalescent simulations conditional on the number of segregating sites. Molecular Biology and Evolution 18, 1136–1138.Google Scholar
Depaulis, F., Mousset, S. & Veuille, M. (2003). Power of neutrality tests to detect bottlenecks and hitchhiking. Journal of Molecular Evolution 57(Suppl. 1), S190–200.Google Scholar
Depaulis, F., Mousset, S. & Veuille, M. (2005). Detecting selective sweeps with haplotype tests. In Selective Sweep (ed. Nurminsky, D.), pp. 34–54. Georgetown, TX: Landes Biosciences.Google Scholar
Donnelly, P. & Tavaré, S. (1995). Coalescent and genealogical structure under neutrality. Annual Review of Genetics 29, 401–421.Google Scholar
Fay, J. C. & Wu, C.-I. (2000). Hitchhiking under positive Darwinian selection. Genetics 155, 1405–1413.Google Scholar
Fu, Y.-X. (1994). Estimating effective population size or mutation rate using the frequencies of mutations of various classes in a sample of DNA sequences. Genetics 138, 1375–1386.Google Scholar
Fu, Y.-X. (1996). New statistical tests of neutrality for DNA samples from a population. Genetics 143, 557–570.Google Scholar
Fu, Y.-X. (1997). Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147, 915–925.Google Scholar
Fu, Y.-X. & Li, W.-H. (1993). Statistical tests of neutrality of mutations. Genetics 133, 693–709.Google Scholar
Glinka, S., Ometto, L., Mousset, S., Stephan, W. & De Lorenzo, D. (2003). Demography and natural selection have shaped genetic variation in Drosophila melanogaster: A multi-locus approach. Genetics 165, 1269–1278.Google Scholar
Griffiths, R. & Tavaré, S. (1994). Simulating probability distributions in the coalescent. Theoretical Population Biology 46, 131–159.Google Scholar
Haddrill, P. R., Thornton, K. R., Charlesworth, B. & Andolfatto, P. (2005). Multilocus patterns of nucleotide variability and the demographic and selection history of Drosophila melanogaster populations. Genome Research 15, 790–799.Google Scholar
Hudson, R. R. (1990). Gene genealogies and the coalescent process. In Oxford Surveys in Evolutionary Biology, vol. 7 (ed. Futuyma, D. & Antonovics, J.), pp. 1–45. Oxford: Oxford University Press.Google Scholar
Hudson, R. R. (1993). The how and why of generating gene genealogies. In Mechanisms of Molecular Evolution (ed. Takahata, N. & Clark, A.), pp. 23–36. Sunderland, MA: Sinauer Associates.Google Scholar
Jakobsson, M., Hagenblad, J., Tavaré, S., Säll, T., Halldén, C., Lind-Halldén, C. & Nordborg, M. (2006). A unique recent origin of the allotetraploid species Arabidopsis suecica: evidence from nuclear DNA markers. Molecular Biology and Evolution 23, 1217–1231.Google Scholar
Kelly, J. (1997). A test of neutrality based on interlocus associations. Genetics 146, 1197–1206.Google Scholar
Kim, Y. & Stephan, W. (2002). Detecting a local signature of genetic hitchhiking along a recombining chromosome. Genetics 160, 765–777.Google Scholar
Kingman, J. F. C. (1982 a). The coalescent. Stochastic Processes and their Applications 13, 235–248.Google Scholar
Kingman, J. F. C. (1982 b). On the genealogy of large populations. Journal of Applied Probability 19A, 27–43.Google Scholar
Kuhner, M. K., Yamato, J. & Felsentein, J. (1995). Estimating effective population size and mutation rate from sequence data using Metropolis–Hastings sampling. Genetics 140, 1421–1430.CrossRefGoogle Scholar
Markovtsova, L., Marjoram, P. & Tavaré, S. (2000). The effects of rate variation on ancestral inference in the coalescent. Genetics 156, 1427–1436.Google Scholar
Markovtsova, L., Marjoram, P. & Tavaré, S. (2001). On a test of Depaulis and Veuille. Molecular Biology and Evolution 18, 1132–1133.Google Scholar
Mousset, S., Derome, N. & Veuille, M. (2004). A test of neutrality and constant population size based on the mismatch distribution. Molecular Biology and Evolution 21, 724–731.Google Scholar
Nordborg, M. (2001). Coalescent theory. In Handbook of Statistical Genetics (ed. Balding, D., Bishop, M. & Cannings, C.), pp. 179–212. Chichester: Wiley.Google Scholar
Pitman, J. (1992). Probability. Berlin: Springer.Google Scholar
Pritchard, J. K., Seielstad, M. T., Perez-Lezaun, A. & Feldman, W. (1999). Population growth of human y chromosomes: a study of Y chromosome microsatellites. Molecular Biology and Evolution 16, 1791–1798.Google Scholar
Przeworski, M. (2002). The signature of positive selection at randomly chosen loci. Genetics 160, 1179–1189.Google Scholar
Przeworski, M. (2003). Estimating the time since the fixation of a beneficial allele. Genetics 164, 1667–1676.Google Scholar
Ramos-Onsins, S. & Mitchell-Olds, T. (2007). mlcoalsim: multilocus coalescent simulations. Evolutionary Bioinformatics 2, 41–44.Google Scholar
Ramos-Onsins, S. E. & Rozas, J. (2002). Statistical properties of new neutrality tests against population growth. Molecular Biology and Evolution 19, 2092–2100.Google Scholar
Rozas, J. & Rozas, R. (1999). DnaSP version 3: an integrated program for molecular population genetics and molecular evolution analysis. Bioinformatics 15, 174–175.Google Scholar
Rozas, J., Sanchez-DelBarrio, J. C., Messeguer, X. & Rozas, R. (2003). DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics 19, 2496–2497.Google Scholar
Simonsen, K. L., Churchill, G. A. & Aquadro, C. F. (1995). Properties of statistical tests of neutrality for DNA polymorphism data. Genetics 141, 413–429.Google Scholar
Sokal, R. R. & Rohlf, F. J. (2001). Biometry. New York: W. H. Freeman.Google Scholar
Stephan, W., Wiehe, T. & Lenz, M. W. (1992). The effect of strongly selected substitutions on neutral polymorphism analytical results based on diffusion theory. Theoretical Population Biology 41, 237–254.Google Scholar
Strobeck, C. (1987). Average number of nucleotide differences in a sample from a single subpopulation: a test for population subdivision. Genetics 117, 149–153.Google Scholar
Tajima, F. (1983). Evolutionary relationship of DNA sequences in finite populations. Genetics 105, 437–460.Google Scholar
Tajima, F. (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123, 585–595.Google Scholar
Tavaré, S. (1984). Line-of-descent and genealogical processes, and their applications in population genetics models. Theoretical Population Biology 26, 119–164.Google Scholar
Tavaré, S., Balding, D. J., Griffiths, R. C. & Donnelly, P. (1997). Inferring coalescence times from DNA sequence data. Genetics 145, 505–518.Google Scholar
Voight, B. F., Adams, A. M., Frisse, L. A., Qian, Y., Hudson, R. R., & Di Rienzo, A. (2005). Interrogating multiple aspects of variation in a full resequencing data set to infer human population size changes. Proceedings of the National Academy of Sciences of the USA 102, 18508–18513.Google Scholar
Wall, J. (1999). Recombination and the power of statistical tests of neutrality. Genetical Research 74, 65–79.Google Scholar
Wall, J. & Hudson, R. R. (2001). Coalescent simulations and statistical tests of neutrality. Molecular Biology and Evolution 18, 1134–1135.Google Scholar
Watterson, G. (1975). On the number of segregating sites in genetical models without recombination. Theoretical Population Biology 7, 256–276.Google Scholar
Wright, S. (1943). Isolation by distance. Genetics 28, 114–138.Google Scholar
Wright, S. I. & Charlesworth, B. (2004). The HKA test revisited: a maximum-likelihood-ratio test of the standard neutral model. Genetics 168, 1071–1076.Google Scholar
Wright, S. I., Vroh Bi, I., Schroeder, S. G., Yamasaki, M., Doebley, J. F., McMullen, M. D. & Gaut, B. S. (2005). The effects of artificial selection on the maize genome. Science 308, 1310–1314.Google Scholar
Figure 0

Fig. 1. Shape index of a coalescent tree. The shape index is the ratio of the branch lengths in the upper and lower parts of the tree, X/Yn (see text).

Figure 1

Fig. 2. Differences between the sizes of FS and Fθ procedures for each S separately. (A) Neutral model: In the upper panel the distribution of S values obtained with the Fθ procedure using θ=0·0057 is shown. In the middle and lower panels the size differences for statistics for the upper (97·5%) and the lower (2·5%) tails, respectively, are presented. Note that large S values, although causing large differences, contribute very little to the total. (B) Hitchhiking model: In the upper panel the distribution of S values obtained with the Fθ procedure using θ=0·0127 is shown. In the middle and lower panels the size differences for statistics for the upper (90%) and the lower (10%) tails, respectively, are presented. Abbreviations are explained in Section 2. Critical values are not calculated for B (lower tail), and for Kw, and Hw (upper tail), because these statistics are not conservative with recombination.

Figure 2

Table 1. Difference (in percentage) between the FS and Fθ procedures for several critical valuesa

Figure 3

Fig. 3. Differences between the sizes of the FS and FSθU procedures. Critical values are obtained for each of the nine neutrality tests for the 2·5% upper and lower tails. Abbreviations are indicated in Section 2.

Figure 4

Table 2. Difference (in percentage) between the FS and the FSθU procedures for several critical valuesa

Figure 5

Fig. 4. Effect of the number of loci on the probability of rejecting the neutral panmictic model for nine neutrality tests. (A) Differences between the sizes of FS and FSθU procedures with no recombination in the upper 2·5% tail, given different numbers of loci. (B) Differences between the sizes of FS and FSθU procedures with no recombination in the lower 2·5% tail, given different numbers of loci. Abbreviations are explained in Section 2. Critical values are not calculated for B (lower tail), and for Kw, and Hw (upper tail), because these statistics are not conservative with recombination. S was fixed at 20 for each locus. Plots obtained using S values from a distribution compatible with θ=0·0057 gave equivalent results (not shown) although for a small number of loci we observed a large variance.

Figure 6

Fig. 5. Posterior densities of the total length Ln of a coalescent tree with n=20 are shown for the different procedures and no recombination. Results based on the FS, FSθ and FSθU procedures are displayed. For FSθ the θ parameter was arbitrarily set to 20/a2(20).

Figure 7

Table 3. Mean and standard deviation (in parentheses) of the posterior distribution of L20

Figure 8

Table 4. Mean and standard deviation (in parentheses) of the tree shape

Figure 9

Fig. 6. Probability distribution of Ln for four alternative models using n=20 and S=20. The FS and FSθU procedures are displayed in each case. Parameter values are the same as in Table 1 but intragenic recombination was set to 4Nr=10 (except for panel D where it was zero). (A) Island model. (B) Logistic expansion model. (C) Bottleneck model. (D) Hitchhiking model.

Figure 10

Fig. 7. Effect of recombination on average Ln for the FS, FSθ (θ=0·0057) and FSθU procedures in a neutral panmictic population using n=20 and S=20.