On a Test of Depaulis and Veuille
Lada Markovtsova, Paul Marjoram, Simon Tavaré · Molecular Biology and Evolution · 2001
In a recent letter to this journal, Depaulis and Veuille (1998) discussed two possible tests of neutrality, the “haplotype number test” and the “haplotype diversity test.” They present in their tables 1 and 2 means and percentage points of the distribution of the number Kn of haplotypes and the sample heterozygosity Hn = 1 − , where are the relative frequencies of those haplotypes in a sample of size n for different values of the number of segregating sites s observed in the data. They assume a neutral infinitely-many-sites model of mutation with no recombination. These percentage points were found by repeatedly simulating a random coalescent tree with n tips, randomly distributing s mutations on the tree, and calculating the observed values of Kn and Hn. See Hudson (1990) for a description of how such simulations can be performed. Depaulis and Veuille (1998) suggested that the percentage points of their statistics could be used as a test of neutrality: for a given sample size n and observed value of s, one compares the observed values of Kn and Hn with the given 95% credible intervals. Values falling outside those intervals lead to rejection of neutrality. Given that the true conditional distributions of Kn and Hn in fact depend on the unknown mutation rate θ, it is likely that the power of this test varies dramatically as a function of θ for given n and s. To assess this hypothesis, we simulated observations from the true joint conditional distribution using equation (1) and then estimated the probability that either statistic would fall outside the limits given by Depaulis and Veuille. This gave an empirical estimate of the probability that their test would reject neutrality for different values of θ. We used a Markov chain Monte Carlo (MCMC) approach to simulate observations from the conditional distribution of (Λ, T) given Sn = s. MCMC methods produce correlated samples, but these samples may be made approximately independent by sampling the output at widely spaced intervals. The results below were generated using the approach in Markovtsova, Marjoram, and Tavaré (2000). An alternative approach is the rejection algorithm of Tavaré et al. (1997) . Table 1 shows the fraction of 10,000 observations that fell inside the DV nominal 95% intervals for three different scenarios. As noted by Fu and Li (1993) , Sn is not a sufficient statistic for θ, so the fraction of observations that fall within the DV nominal 95% intervals varies greatly. It appears that if the true θ value is well supported by the data, the test of neutrality based on the DV intervals will work well. However, if the true θ value is not well supported by the data, the test will be inaccurate, leading to incorrect rejections of neutrality. For further discussion, see Wall and Hudson (2001). As we have seen, the power of the DV test depends on the unknown parameter θ. Even if the mutation rate is known, the compound parameter θ still depends on the underlying effective population size at the time of sampling, and this is not known in practice. Depaulis and Veuille (1998) discuss a data set from the Su(H) locus in Drosophila melanogaster for which n = 20, s = 44. The observed values of K20 and H20 were 7 and 0.76, respectively. The nominal P values were estimated to be 0.011 and 0.017, respectively. We used the simulation approach outlined above to find empirical estimates of these P values for different values of θ, using 10,000 simulated values once more. We used values of θ = 1, 5, 12.4, 50.0, and 100 for illustration; the value 12.4 corresponds to Watterson's (1975) segregating-sites estimator. Results are shown in table 2 . From table 2 , we see that for θ in the range 12.4 or larger, the data are highly unlikely under a neutral scenario. However, for a range of smaller θ values, including, for example, θ = 5, the data become much more likely. Thus, the ability to reject the supposed neutral scenario depends on the true value of θ. It should be noted that one cannot reject neutrality on the basis of this test; rather, one can reject the model upon which the test is based. As Depaulis and Veuille (1998) note, this model does not include recombination (although it is easy to alter it to do so); neither does it include any population demographics such as stratification. If the model is rejected, any of these missing factors could be the cause, not necessarily the assumption of neutrality. Computer programs that implement both the MCMC approach and the rejection method to generate observations from the joint conditional distribution of, for example, (Kn, Hn) given the value of Sn for a given distribution for θ under a variety of demographic scenarios can be obtained from the authors. Yun-Xin Fu, Reviewing Editor Keywords: Markov chain Monte Carlo coalescent Address for correspondence and reprints: Simon Tavaré, Program in Molecular Biology, Department of Biological Sciences, SHS 172, University of Southern California, Los Angeles, California 90089-1340. [email protected]. Table 1 Fraction (%) of 10,000 Observations Falling Within the DV Interval Table 1 Fraction (%) of 10,000 Observations Falling Within the DV Interval Table 2 P Values for Drosophila melanogaster Data Table 2 P Values for Drosophila melanogaster Data We thank Y.-X. Fu and an anonymous reviewer for helpful comments. We were supported in part by NSF grant DBI95-04393 and NIH grant GM 58897.