Directional Selection and the Site-Frequency Spectrum (original) (raw)

Abstract

In this article we explore statistical properties of the maximum-likelihood estimates (MLEs) of the selection and mutation parameters in a Poisson random field population genetics model of directional selection at DNA sites. We derive the asymptotic variances and covariance of the MLEs and explore the power of the likelihood ratio tests (LRT) of neutrality for varying levels of mutation and selection as well as the robustness of the LRT to deviations from the assumption of free recombination among sites. We also discuss the coverage of confidence intervals on the basis of two standard-likelihood methods. We find that the LRT has high power to detect deviations from neutrality and that the maximum-likelihood estimation performs very well when the ancestral states of all mutations in the sample are known. When the ancestral states are not known, the test has high power to detect deviations from neutrality for negative selection but not for positive selection. We also find that the LRT is not robust to deviations from the assumption of independence among sites.

THE Poisson random field (PRF) models of Sawyer and Hartl (1992) and Hartl et al. (1994) afford a general likelihood framework for estimating mutation and selection parameters in various population genetic settings. This is currently the only method available for estimating selection directly from DNA sequence data in population genetics. The PRF model has proven quite useful in studying selection on “silent” site polymorphism in Drosophila (Akashi 1995; Akashi and Schaeffer 1997) and Escherichia coli (Hartl et al. 1994) as well as amino acid variation in E. coli and Salmonella enterica (Hartl et al. 2000). It has also been used to study the power of various tests of neutrality on the basis of polymorphism and divergence data (Akashi 1999).

Since the PRF model assumes independence among sites, it is a limiting case for population genetic models that incorporate mutation, selection, and recombination. For typical DNA sequence data, in which many sites segregate simultaneously, the assumption of independence among sites is equivalent to the assumption that there is free recombination between segregating nucleotides. Studying statistical inference in the PRF framework can therefore illuminate efforts to build models that relax the assumption of independence among sites. The reason for this is that any model that hopes to estimate mutation, selection, and recombination parameters can do as well as or worse than the PRF if the data conform to the PRF model.

It would be of great interest to know, therefore, how well maximum-likelihood methods for point estimation and hypothesis testing in the PRF framework perform under varying levels of selection and mutation. One property of interest is the probability of rejection neutrality if the locus is truly under selection, which addresses the power of the likelihood-ratio tests (LRTs) of neutrality. Another statistical issue of interest concerns the robustness of these models to deviations from the assumption of independently evolving sites, since this assumption may be unrealistic for certain types of data.

This article treats in detail the PRF model for site-frequency data proposed by Hartl et al. (1994) for both folded and unfolded site configurations. In particular, we (1) derive the asymptotic variances and covariance of the maximum-likelihood estimates (MLEs) of the selection and mutation parameters, (2) explore the power of the likelihood-ratio test of no selection under varying levels of mutation and selection, (3) explore the coverage of two types of confidence intervals for the selection parameter, and (4) explore the robustness of the model to deviations from the assumption of independence among sites.

THEORY

Poisson random field model: First let us consider a simple example with which to illustrate the type of data being modeled. The data in Table 1 represent a sample of six individuals sequenced at five variable DNA sites. Each column is an aligned DNA site and each row is an individual sampled from the population. At each site where the individual carries the derived nucleotide, the resulting entry in the matrix is 1; otherwise the entry is 0. For each site, then, the sum of the column gives

TABLE 1

A simple example of the type of data modeled in this article

The site-frequency spectrum for this sample is (3, 1, 0, 1, 0) as explained in the text.

TABLE 1

A simple example of the type of data modeled in this article

The site-frequency spectrum for this sample is (3, 1, 0, 1, 0) as explained in the text.

the number of individuals in the sample that carry the derived form of the mutation.

Define the unfolded site-frequency spectrum as the random vector X = (_X_1, _X_2, … , _Xn_−1) of sample configurations Xi, where Xi represent the number of sites that have n − i ancestral and i derived nucleotides for 1 ≤ i ≤ n − 1 among the n aligned nucleotide sequences. For the example in Table 1 the unfolded site-frequency spectrum is (3, 1, 0, 1, 0) since there are three sites where one individual carries the derived mutation, one site where two individuals carry the derived mutation, no sites where three individuals or five individuals carry the derived mutations, and one site where four individuals carry the derived mutation.

So far we have assumed that for each site we know which state is derived and which is ancestral. If this information is unknown, we would model the folded site-frequency spectrum, _X_′, defined as the number of sites where 1 individual or (n − 1) individuals carry a mutation, 2 individuals or (n − 2) individuals carry a mutation, etc., so that Xi′=Xi+Xn−i for i < n − i and Xn∕2′=Xn∕2 for i = n/2. For the example above, the folded site-frequency spectrum is (3, 2, 0).

We wish to model the site-frequency spectrum for DNA sites in an infinite-sites (i.e., irreversible mutation) independent (i.e., freely recombining) model with (weak) directional haploid selection. We assume that all mutations that enter the population have the same selective effect and that there are no epistatic interactions among mutations (i.e., independent fitness effects of all mutations). We also assume that there is no population structure to the data and the population has achieved statistical stationarity.

Fisher (1930) and Wright (1938, 1969) derived a formula for what Wright called the “transient distribution” for the frequency of a single newly arisen mutation under selection. Letting f(q, γ) be the density of the frequency of a single mutation, q, in the range dq, their formula is

f(q,γ)=1−e−2γ(1−q)1−e−2γ2q(1−q),

(1)

where γ is the product of the selection coefficient (a in Fisher's notation, usually s in Wright's notation) and the effective haploid population size. Here s is the pergeneration increase in log(p/q) so that Δ_q_ = s p(1 − p) if s is small.

Sawyer and Hartl (1992) showed that for recurrent mutation at a locus with independently evolving sites, the expected density of a random field composed of the individual frequencies of the derived mutations is (θ/2)f(q, γ), where θ/2 is the per-locus mutation rate scaled by the effective population size. They also showed that this random field is a PRF, which means that the number of sites whose derived mutation frequency q satisfies 0 < a < q < b < 1 has a Poisson distribution for given a, b and also that the number of sites for nonoverlapping intervals (a, b) is independent (Kingman 1993).

For a particular site, given the frequency q of the mutation in the population, the probability of sampling i sequences of one type and n − i of another is binomially distributed with parameters n and q. If mutations enter the population as a Poisson process with rate θ/2, the Xi's are independent Poisson-distributed random variables with mean, θ_F_(i, γ), where

F(i,γ)=∫01f(q,γ)2⋅Pr(i∣q)dq=∫011−e−2γ(1−q)1−e−2γ1q(1−q)(ni)qi(1−q)n−idq.

(2)

Therefore, the probability of observing xi sites that have i derived and (n − i) ancestral mutations is

p(Xi=xi∣θ,γ)=e−θF(i,γ)(θF(i,γ))xixi!.

(3)

The PRF model outline above leads directly to a likelihood-ratio test of neutrality, which compares the null hypothesis γ = 0 with the alternative hypothesis that γ ≠ 0. Since the Xi's are independent, the likelihood function is the product of the individual p(Xi = xi|θ, γ). Letting Lf(θ, γ | x) be the likelihood function for the folded model and Lu be the likelihood function for the unfolded model,

Lu(θ,γ∣x)=Πi=1n−1e−θF(i,y)(θF(i,γ))xixi!

(4)

and

Lf(θ,γ∣x′)=Πi=1⌊n∕2⌋e−θG(i,γ)(θG(i,γ))xi′xi′!,

(5)

where ⌊i_⌋ is the greatest integer ≤_i and

G(i,γ)={F(i,γ+F(n−i,γ))ifi<n−i,F(n2,γ)ifi=n2,0otherwise.}

(6)

To perform the likelihood-ratio test, we need to maximize (4) or (5) for θ and γ under (A) the unconstrained state space and (B) the constraint that γ = 0. Letting γ^ and θ^ be the unconstrained MLEs of γ and θ, and θ^W the MLE of θ under neutrality, the likelihood-ratio test statistic, Λ, is

Appealing to large sample results (Kendall 1987), 2 ln(Λ) is ~χ12 distributed since the neutral and selected models differ by 1 d.f.

Maximum-likelihood estimation and asymptotic results: Finding the maximum-likelihood estimates of θ and γ for both the folded and unfolded site frequency models is relatively straightforward. Let l be the log-likelihood function for the unfolded model

l(θ,γ∣x)∝∑i=1n−1xilog(θF(i,γ))−θ∑i=1n−1F(i,γ).

(8)

The maximum-likelihood estimates of γ and θ (i.e., γ^ and θ^ are usually found by setting the partial derivatives of the log-likelihood function equal to zero. Letting S represent the total number of segregating sites (⁠S=∑i=1n−1xi⁠),

∂l∂γ=−θ∑i=1n−1F′(i,γ)+∑i=1n−1xiF′(i,γ)F(i,γ),

(10)

where

F′(i,γ)=∫01f′(q,γ)2(ni)qi(1−q)n−idq

(11)

and

f′(q,γ)=4e−2γ(1−q)(1−e−2γ)q−4e−2γ(1−e−2γ(1−q))(1−e−2γ)2q(1−q).

(12)

An efficient approach to solving for the MLEs is to maximize the profile likelihood that has the same global maximum but requires only maximization in one dimension. The log-profile-likelihood function, l*(γ), is l(θ~,γ)⁠, where θ~ is the maximum-likelihood estimate of θ given γ. By setting (9) to zero and solving for θ~ one can easily show that

It follows from (2) that

which proves that under neutrality

Here θ^W is Watterson's (1975) estimator of θ. This shows that the total number of segregating sites is a sufficient statistic for θ given γ. The two preceding results hold, regardless of the folding of the site configurations since ∑i=1n−1F(i,γ)=∑i=1[n∕2]G(i,γ)⁠.

We can now maximize the log-profile-likelihood function l*(γ) numerically, using a standard optimization technique such as Newton-Raphson iteration (Lange 1999). To implement Newton-Raphson we need to find the first and second derivatives of l*(γ) with respect to γ. To do this we first substitute (13) into (8) and take the necessary derivatives. The results are

dl∗dγ=∑i=1n−1xiF′(i,γ)F(i,γ)−S∑i=1n−1F′(i,γ)∑i=1n−1F(i,γ)

(16)

and

d2l∗dγ2=∑i=1n−1xiF′′(i,γ)F(i,γ)−F′(i,γ)2F(i,γ)2−ST0T2−T12T02,

(17)

where

and

F′′(i,γ)=∫01f′′(i,γ)2(ni)qi(1−q)n−idq,

(21)

where

f′′(q,γ)=(1−e−2γ(1−q)q(1−q))(16e−4γ(1−e−2γ)3+8e−2γ(1−e−2γ)2)−16e−2γ−2γ(1−q)(1−e−2γ)2q−8(1−q)e−2γ(1−q)(1−e−2γ)q.

(22)

To find the equivalent results for the folded model, substitute n by [n/2] and _F, F_′, and _F_″ by _G, G_′, and _G_″, respectively, in the equations above, where

G′(i,γ)={F′(i,γ)+F′(n−i,γ)ifi<n−i,F′(n2,γ)ifi=n2,0otherwise}

(23)

G′′(i,γ)={F′′(i,γ)+F′′(n−i,γ)ifi<n−i,F′′(n2,γ)ifi=n2,0otherwise}

(24)

To find the asymptotic covariance matrix of the maximum-likelihood estimates, we compute the Fisher information matrix I(θ, γ), which is minus the expected value of the matrix of second derivations of the log-likelihood function. Since the components of I are linear in xi and S, this is equivalent to replacing xi by E(Xi) = θ_F_(i, γ) and S by θ_T_0(γ). The asymptotic covariance matrix of θ^ and γ^ is the inverse of the empirical Fisher information matrix I (⁠θ^,γ^⁠).

The components of the Fisher information matrix I (θ, γ) are easily shown to be

I12=I21=−E[∂2l∂γ∂θ]=T1(γ)

(26)

I22=−E[∂2l∂γ]=θ(T2(γ)−t(γ)),

(27)

where

t(γ)=∑i=1n−1F′′(i,γ)F(i,γ)−F′(i,γ)2F(i,γ).

(28)

The asymptotic covariance of the MLEs is the inverse matrix of I (⁠θ^,γ^⁠), which by Cramer's rule equals

M=1I11I22−I122(I22−I12−I12I11)

(29)

evaluated at θ^,γ^⁠. It then follows that

VAR(θ^)=M11=θ^T0(γ^)−[T1(γ^)2∕(T2(γ^)−t(γ^))]

(30)

COV(γ^,θ^)=M12=1[T0(γ^)(T2(γ^)−t(γ^))]∕T1(γ^)−T1(γ^)

(31)

VAR(γ^)=M22=1θ^(T2(γ^)−t(γ^)−T1(γ^)2∕T0(γ^)).

(32)

Note that the negative inverse of the second derivative of the log-profile-likelihood function (17), with xi replaced by E(Xi) = θ_F_(i, γ) and θ by S/_T_0, is formally the same as (32). this is consistent with the fact that both expressions give the asymptotic variance of γ^⁠.

From general likelihood theory, we expect (⁠θ^,γ^⁠) in large samples to be distributed as a multivariate normal with mean (θ, γ) and variance-covariance matrix M. Using this result and the fact that 2 ln(Λ) is approximately χ12 distributed, we can define two types of confidence sets for γ at given confidence level 1 −α.

The first confidence set consists of γ in the interval γ^±Z(1−α∕2)VAR(γ^)⁠, where Z(1 −α/2) is the standard normal quantile corresponding to level 1 − α/2 (e.g., Z(0.975) = 1.96). The second confidence set corresponds to the set of values of γ for which we would not reject the likelihood-ratio test: γ such that l∗(γ)≥l∗(>^)−0.5χ1,1−α2 (Kendall 1987).

SIMULATIONS

The two properties of the LRT that we are interested in exploring are the size and power of the test.

The size, confidence level, or probability of type I error of a test is defined as the probability of rejecting a true null hypothesis and is denoted α. The power of a test is the probability of rejecting a false null hypothesis and is denoted 1 −β, where β is the probability of not rejecting a false null hypothesis (i.e., the probability of a type II error).

In particular, we explore whether the LRT maintains the proper size even if the assumption of independence among segregating sites is not met. If the test is robust to the assumption of free recombination and the null hypothesis is true, then in repeated sampling we should reject the null hypothesis (100 · α)% of the time. If the test is not robust to the assumption of free recombination, the realized size of the test will vary with the rate of recombination. We might also expect the realized size of the test to vary with the mutation rate if the test is not robust to recombination, since a larger mutation rate will produce a test with a higher probability of rejecting the null hypothesis. It would also be of interest to know how the power of the LRT changes as the strength of selection changes for data that conform to the independence assumption.

To explore these two issues, as well as the coverage of the two different types of confidence intervals for γ^ we simulated data under two types of models. The first type of data was generated using coalescent simulations for neutral data under varying rates of recombination and mutation (R. Hudson, personal communication). These data sets were used to explore the robustness of the test to deviations from independence among sites. The second type of data was generated using the PRF model outlined above by sampling independent Poisson random variables with rates given by mutation and selection parameters of interest. These data were used to explore the power of the test under varying levels of mutation and selection, as well as the sampling distributions of θ^W⁠, θ^⁠, and γ^⁠. All of the data sets generated had

TABLE 2

Proportion of LRTs that reject neutrality at the 5% level for neutral loci evolving under varying levels of recombination

R refers to the population recombination rate, 4Ne_r_, set in the coalescence simulation.

TABLE 2

Proportion of LRTs that reject neutrality at the 5% level for neutral loci evolving under varying levels of recombination

R refers to the population recombination rate, 4Ne_r_, set in the coalescence simulation.

Figure 1.

Sampling distribution for unfolded site configurations under varying levels of recombination. The _x_-axis for each graph is the number of sites that were found at frequencies 1/10, 2/10, … , 9/10. Each line is a replicate data set of the 500 simulated data sets for each level of recombination.

n = 10 sequences and consisted of 500 replicates. The same data were used for all folded and unfolded analysis, with the folded site configurations generated by summing 1 and (n − 1), 2 and (n − 2), etc.

All numerical integration was done using Romberg integration, maximization by Newton-Raphson iteration with bisection, and random number generation using standard numerical algorithms (Press et al. 1988). Software for calculating the PRF likelihood-ratio test is available from the authors upon request.

RESULTS AND DISCUSSION

Size and power of the LRT: In Table 2 we summarize the results for the analysis of the size of the PRF test for both folded and unfolded site configurations under two different mutation rates and 10 levels of recombination. From this analysis we see that the Poisson random field LRT is extremely sensitive to the assumption of independence among sites. If the test were not sensitive to this assumption, we would expect the entries in the table to fluctuate ~5% but we see rejection levels that reach upward of 50% for tight linkage. While the strength of the effect declines quickly as the rate of recombination increases, the tests do not attain level α until the rate of recombination is more than an order of magnitude greater than the mutation rate.

The reason for this becomes clear when one examines the effect of recombination rate on the distribution of the site-frequency spectrum, X. In Figure 1 we present representative replicates of X from the coalescent simulations (γ = 0) under varying levels of recombination for θ = 10. The expected value of the number of sites at a given frequency i/n is the same for each of the six graphs presented [E(Xi) = θ/_i_]. The variance of Xi, though, increases as the recombination rate tends to

Figure 2.

Power of the LRT under varying levels of selection and mutation for folded (F) and unfolded (U) data configurations. On the _x_-axis we plot γ, the value of the selection parameter in the PRF model under which the data were simulated, and the lines connect the proportion of data sets (n = 500 for each point) that reject γ = 0.

Figure 3.

Distributions for the maximum-likelihood estimates of γ (⁠γ^⁠) as a function of the simulated value of γ for folded (F) and the unfolded (U) frequency data at two levels of mutation. The open regions correspond to the space containing 95% of the observed γ^'s and the dashed lines connect the median estimates across values of γ.

zero. The reason, therefore, that we reject neutrality more than expected in Table 2 with increasing linkage is that we underestimated the variance for the site frequencies. Likewise, if neutral sites are linked to selected sites there will also be an increase in the variance in the distribution of the X. These results suggest that the analysis of single genes using this method is invalid if the rate of recombination is not extremely large, since one rejects the null hypothesis (which of course includes the assumption of free recombination) if the data are neutral but not independent.

Performance of maximum likelihood for point estimation in the PRF framework: While the method may not be applicable for linked regions, the analysis of genome-wide single-nucleotide polymorphism may have the quasi-independence required for the test to maintain level α. It would be desirable, therefore, to know how well the method performs in moderate sample sizes when the data conform to the PRF model. For the remainder of this article, all of the analysis is performed on data generated under the PRF model.

In Figure 2 we summarize the results for the analysis of the power of the PRF test for both folded and unfolded site configurations under the same mutation rates above for γ ε [−10, 10]. We see from that that the LRT has very good power to detect deviations from

Figure 4.

The log-profile-likelihood functions of γ for the folded and unfolded site configurations of the data set in the text with peaks at γ = ∞ for the unfolded model and at γ = −0.4734 and γ = ∞ for the folded model (true γ = 10 as explained in the text).

Figure 5.

Data and fitted values for the example discussed in the text. The bar graph is the data and the line graphs are the fitted values under the neutral model (dashed line), infinitely strong positive selection model (solid line), and under the MLEs of the PRF model (circles).

neutrality even in moderate sample sizes when the ancestral states of variable sites are known (i.e., for unfolded site configurations). A surprising result is that for folded site configurations the test has strong power to detect negative selection but very weak power to detect positive selection. The reason for this is related to the issue of point estimation in general in the PRF framework.

From general-likelihood theory, we expect the maximum-likelihood estimate for a parameter under certain regularity conditions to be consistent, unbiased, and efficient as the sample size used to estimate the parameter tends to infinity. By “consistent” we mean that the parameter estimate will converge in probability to the true value of the parameter (i.e., the variance in the parameter estimate will tend toward 0), by “unbiased” we mean that the expected value of the parameter estimate is centered on the true value of the parameter, and by “efficient” we mean that the estimator has the lowest possible variance of all possible estimators. In practice,

Figure 6.

Distribution for the maximum-likelihood estimates of θ as a function of the simulated value of γ. The open region corresponds to the space containing 95% of the observed θ^'s.

Figure 7.

Coverage of 95% confidence intervals for γ as a function of the simulated value of γ for two mutation levels for both folded (F) and unfolded (U) site-frequency spectra. Dotted lines connect the percentage of confidence intervals that contain the true γ for confidence intervals based on l* and dashed lines connect the percentage of confidence intervals that contain the true value for confidence intervals based on the normal approximation and asymptotic variance.

though, maximum-likelihood estimation can often behave very poorly if the sample size is not large enough or if the regularity conditions are not met.

For the levels of mutation examined in our generated data, we find that the distributions for γ^ are reasonably centered around the true value of γ for the unfolded model (Figure 3) with increasing variance as one moves away from neutrality. We note that in regions of extreme selection (|γ| > 5) the distribution for γ^ becomes quite skewed with γ^=∞ or γ^=−∞ for certain replicate data sets. The reason for this is that the log-likelihood function becomes very flat when the magnitude of γ is very large. In the case of negative selection, samples consisting of only singletons are consistent with extreme negative selection and very high mutation, so that γ^ becomes indistinguishable from −∞ and we show shortly why the likelihood function asymptotes to ∞ for any given data set. These results indicate that these extrema represent either a local maximum or minimum that needs to be checked once an MLE for γ has been estimated. If it is found that the likelihood at one of these two extremes is higher than the estimate γ^⁠, this is problematic for point estimation but not for hypothesis testing since we can still test whether the observed sample configuration is consistent with a neutral model (for most of the cases we examined where γ^=∞ or γ^=−∞ for unfolded data, the LRT is significant).

While maximum-likelihood estimation works reasonably well for point estimation using unfolded data, the method performs very poorly for folded data in the region of positive selection as alluded to above in the section on power. As we see in Figure 3, for the folded data the method performs well in terms of estimating γ for the region of negative selection, performs reasonably well for a small region of positive selection, and then asymptotes to neutrality for arbitrarily large γ.

For a better understanding of this phenomenon, let us focus on one of the replicate data sets we generated under θ = 10 and γ = 10, x = (11, 6, 8, 1, 2, 6, 2, 3, 10). For these data, θ^=8.66⁠, γ^=∞⁠, and LRT = 18.08 (P < 0.0001) for the unfolded model. We note that the profile-likelihood function for the unfolded data is monotonic with a peak at ∞ (Figure 4). As expected from the LRT, we see in Figure 5 that the predictions from the infinite selection model are quite different from the predictions of the neutral model. For the folded model, θ^=20.03⁠, γ^=−0.473⁠, and LRT = 0.18 (P = 0.328). As expected from the LRT, we see in Figure 5 that the folded site-frequency spectrum is fitted very well by a neutral model. What may be surprising is that the data are equally well fit by a model with infinitely strong positive selection! Another interesting property of these data is that the likelihood function for folded data configurations has a bimodal nature with a region in between that can be significantly different from either peak, implying that the data are consistent with either neutrality or strong selection but not weak selection (Figure 4).

The reason for this phenomenon is that when the site configurations are folded, the limiting distribution as γ → ∞ is indistinguishable from neutral evolution with a mutation rate that is twice as large as the true mutation rate. This can be shown by noting that the log-likelihood function for a given data set will reach finite asymptotic value for γ → ∞ since

Recalling that G(i, γ) = F(i, γ) + F(n − i, γ) and the result derived earlier that F(i, 0) = 1/i, we can easily see that

Figure 8.

Log-profile-likelihood function (solid line) and its normal approximation (dashed line) under varying sample size. (S represents the total number of segregating sites in the sample). The dotted line corresponds to 1.92 log-likelihood units from the maximum. The intersection of this line and the log-profile likelihood or its approximation corresponds to the 95% confidence interval constructed from the method.

leading to a bimodal log-likelihood function in the extreme case with maxima near 0 and ∞ as seen in the example above.

It is important to note that this is a general problem for methods that attempt to detect selection using folded site configurations [_e.g._, Tajima's (1989) _D_-statistic]. This also explains the odd result reported in Akashi (1999) that Tajima's _D_-statistic has little to no power to detect positive selection. Since this is a general problem for estimating selection from population genetics data, the inclusion of a parameter for recombination would make the problems of inference even worse, since in this case we know that the data conform to the assumptions of the model.

One interesting point related to this issue is the effect of the estimation of γ on the estimation of θ. We see from Figure 6 that estimation of θ is excellent for the unfolded model for regions of positive selection. It may seem paradoxical that while the variance in the estimates of γ increases with increasing level of selection, the variance in the estimates of θ decreases. The reason for this is that once selection becomes large the predictions of the PRF model for the site-frequency spectrum become indistinguishable (this is exactly the reason for the asymptote at ∞ for positive selection), but only a small range of θ is consistent with the observed number of segregating sites. In the case of negative selection, the results are reversed since as selection increases in the negative selection the values for F(i, γ) become very small, leading to large variance in the estimate of both selection and mutation.

For the folded model, in Figure 6 we begin to see the effect of the maximum near 0 for γ^ on the estimation of mutation. Watterson's estimator of θ will always be larger than the estimated value of θ under positive selection, since F(i, γ) increases with increasing γ and reaches its maximum as explained above at twice the value of the partial harmonic sum in Watterson's estimator. This means that since we tend to underestimate selection in the folded model for positive selection, we will tend to overestimate θ.

Confidence sets for γ**:** The last issue we explore is the construction of confidence intervals for the selection parameter in the PRF framework. As we outlined above our two methods for constructing confidence intervals are (1) the region that contains (1 −α) · 100% of the area in the normal approximation to likelihood function and (2) the region in the profile likelihood space that is α0.5χ1,1−α2 likelihood units from the maximum-likelihood point. For the first confidence set, at α = 0.05 this corresponds to the area that is 1.96 standard deviations away from γ^ where the standard deviation is given by VAR(γ^)⁠. Likewise for the second set, at α = 0.05 this corresponds to the region where l∗(γ^)−l∗(γ)>3.84∕2 likelihood units.

We see from Figure 7 that confidence intervals based on the normal approximation to the likelihood tend to undercover (i.e., they tend to be too small), particularly for the folded model and for strong positive selection. The confidence intervals based on the profile likelihood, however, have excellent coverage regardless of the folding of the data or the strength of selection.

The reason for this is seen in Figure 8, where we plot the log-profile-likelihood function and its normal approximation for a data set with the same γ^ but varying sample sizes (in terms of the number of segregating sites for the sample). We see that the number of segregating sites needs to be quite large before the normal approximation begins to have the same coverage as the profile likelihood in a region of 1.92 log-likelihood units. This occurs because a confidence interval based on the normal approximation will, by definition, be symmetric, and the log-likelihood function is quite right skewed (as we explained above). Since we need an immensely large sample for the region around the maximum-likelihood point to be normally distributed (upward of 1000 segregating sites), confidence intervals based on the normal approximation will tend to undercover.

CONCLUSIONS

The Poisson random field model for site-frequency data holds great potential for untangling the effects of mutation and selection on standing genetic variation. Our analysis suggests, however, that when the data do not conform to the assumptions of very loose linkage, the estimates from this model may be misleading. Likewise, if the data set being analyzed is small, one should be cautious about the asymptotic assumptions implicit in maximum likelihood as a method of point estimation. This is a general problem for methods that attempt to estimate selection using site-frequency data, and methods that attempt to incorporate recombination into point estimation will only make some of these problems worse. The solution may be to focus on genetic variation that is independent (such as the growing number of single nucleotide polymorphisms for many organisms) rather than focusing on methods that detect selection at single loci. We also find that folding the site-frequency spectrum obscures positive selection and suggest that every effort should be made to incorporate information on the ancestry of mutations being studies, if one wishes to estimate selection.

We also point out that simply adding information on the number of fixed differences in the region between the study population and a close species will not aid in point estimation for selection. Such a method would add exactly one data point and one parameter (time since divergence), which can be scaled to fit that point exactly. A method that uses two classes of sites (e.g., silent and replacement) and fixed differences circumvents this problem and will add useful information for the estimation of selection. We have developed a method that uses divergence data and the site-frequency spectrum to estimate selection. Likewise, information on variation in the site-frequency spectrum among different genomic regions also holds great potential for refining point estimation for single regions of interest. We have also developed a method that uses Markov chain Monte Carlo to perform Gibbs and Metropolis sampling for a hierarchical model of this form. We will expand upon these two methods in subsequent publications.

Acknowledgement

We thank Richard Lewontin, Rasmus Nielsen, and Steve Palumbi for many helpful comments. This work was supported by a Howard Hughes Medical Institute Graduate Fellowship to C.D.B. and National Science Foundation grant DMS9707045 to S.A.S.

Footnotes

Communicating editor: N. Takahata

LITERATURE CITED

Akashi

H

,

1995

Inferring weak selection from patterns of polymorphism and divergence at “silent” sites in Drosophila DNA

.

Genetics

139

:

1067

–

1076

.

Akashi

H

,

1999

Inferring the fitness effects of DNA mutations from polymorphism and divergence data: statistical power to detect directional selection under stationarity and free recombination

.

Genetics

151

:

221

–

238

.

Akashi

H

,

Schaeffer

S W

,

1997

Natural selection and the frequency distributions of “silent” DNA polymorphism in Drosophila

.

Genetics

146

:

295

–

307

.

Fisher

R A

,

1930

The Genetical Theory of Natural Selection

.

Clarendon Press

,

Oxford

.

Hartl

D L

,

Moriyama

E N

,

Sawyer

S A

,

1994

Selection intensity for codon bias

.

Genetics

138

:

227

–

234

.

Hartl

D

,

Boyd

E F

,

Bustamante

C D

,

Sawyer

S

,

2000

The glean machine: what can we learn from DNA sequence polymorphism

, pp.

37

–

49

in

Genomics and Proteomics: Functional and Computational Aspects

, edited by

Suhai

S

.

Plenum Press

,

New York

.

Kendall

M

,

1987

Kendall's Advanced Theory of Statistics

, Vol.

2

.

Oxford University Press

,

Oxford

.

Kingman

J

,

1993

Poisson Processes

.

Oxford University Press

,

Oxford

.

Lange

K

,

1999

Numerical Analysis for Statisticians

.

Springer-Verlag

,

New York

.

Press

W H

,

Flannery

B P

,

Vetterling

S A T W T

,

1988

Numerical Recipes in C

.

Cambridge University Press

,

Cambridge, UK

.

Sawyer

S

,

Hartl

D L

,

1992

Population genetics of polymorphism and divergence

.

Genetics

132

:

1161

–

1176

.

Tajima

F

,

1989

Statistical method for testing the neutral mutation hypothesis by DNA polymorphism

.

Genetics

123

:

585

–

595

.

Watterson

G A

,

1975

On the number of segregating sites in genetic models without recombination

.

Theor. Popul. Biol.

7

:

256

–

276

.

Wright

S

,

1938

The distribution of gene frequencies under irreversible mutation

.

Proc. Natl. Acad. Sci. USA

24

:

253

–

259

.

Wright

S

,

1969

Evolution and the Genetics of Populations, Vol. 2: The Theory of Gene Frequencies

.

University of Chicago Press

,

Chicago

.

© Genetics 2001