COALESCENT SIMULATION

This study requires a coalescent simulation in which both intragenic recombination and directional selection take place during the ancestry of a DNA sample. The ancestral recombination graph (ARG) described by Griffiths and Marjoram (1997) allows the realization of intragenic recombination for any length of a DNA sequence. We modified the ARG to incorporate the effect of genetic hitchhiking. Figure 1 illustrates an ARG with hitchhiking for a sample of three sequences. The ancestral history of the sample is divided into neutral phases and a selective phase.

During a neutral phase the ARG is constructed as described by Griffiths and Marjoram (1997), with the following modifications. Edges of the ARG represent ancestral sequences at a given time (T) in the past measured in units of 2N generations, where N is the effective number of individuals in a diploid population. If there are k edges at a given time (k = n at T = 0), either coalescent or recombination events can occur with rates k(k − 1)/2 and kR/2, respectively. R is 4N times the recombination rate between both ends of the sequence. Therefore, if the sequence is L nucleotides long, R = 4NLρ, where ρ is the per-nucleotide recombination rate. Each of the k edges is labeled by a pair of integers (I_i, J_i) (i = 1, … , k). This pair of integers delimits the region within which sequences ancestral to sample sequences are found. Therefore, recombination outside this region can be ignored. At T = 0, (I_i, J_i) = (1, L) for all n edges. At a coalescence event, two randomly chosen edges (for example, the lth and mth edges) join to a new edge, which is then labeled by (Min(I_l, I_m), Max(J_l, J_m)). At a recombination event, an edge is chosen randomly and a random uniform integer, U, is drawn between 1 and L − 1. If an edge labeled by (I, J) was chosen, it joins to two parental edges only if I ≤ U < J. Then, the two parental edges are labeled as (I, U) and (U + 1, J). If U < I or > J, no change is made at the edge. This procedure is necessary to minimize k in the simulation.

The selective phase is the period when a substitution of a beneficial mutation that causes a hitchhiking effect takes place. The beneficial allele B has a genic selective advantage s over the parent allele b. This substitution occurs at a site M nucleotides away from the left end of the sequence and the fixation of B is completed at T = τ. The allele frequency of B, x, is assumed to change deterministically from 1 − ζ to ζ. Therefore, x at T = τ + t is given by x(t)=ξξ+(1−ξ)eα(t−ts)(0≤t≤ts) (1) (Stephanet al. 1992), where α = 2Ns and t_S = −(2/α)log(ζ), which is the length of the selective phase. The choice of ζ does not change the resulting genealogy significantly (Bravermanet al. 1995). We use ζ = 100/(2N) for the simulations. During the selective phase, B and b edges exist, indicating whether an ancestral sequence includes the beneficial allele or not. Therefore, all edges are B edges at the beginning of the selective phase (T = τ). The system of labeling edges is also changed: (I, J) at the end of the neutral phase at T = τ changes to (Min(I, M), Max(J, M)). This change means that recombination between the site of directional selection and the ancestral sequence should be followed during the selective phase. There are four possible events during the selective phase: (1) coalescence between B edges; (2) coalescence between b edges; (3) recombination in a B edge; (4) recombination in a b edge. The probability of these four events during the time interval [t, t + Δt] is given by λ1(t)Δt=kB(kB−1)2x(t)Δt, (2a) λ2(t)Δt=kb(kb−1)2(1−x(t))Δt, (2b) λ3(t)Δt=kBR2Δt, (2c) λ4(t)Δt=kbR2Δt, (2d) respectively, where k_B and k_b are the numbers of B and b edges at T = τ + t, respectively. The waiting time, Δt, between events is randomly drawn from an exponential distribution with parameter λ_T(t) = λ₁(t) + λ₂(t) + λ₃(t) + λ₄(t). Then, one of the four events is allowed to occur according to its probability. This method should be used only when waiting time, Δt, is short («1/α) such that the change of x(t) between events is negligible. In this study, due to large values of R, values of Δt are sufficiently small. With lower values of R, a rejection method such as the one by Braverman et al. (1995) should be used. When a recombination event occurs in a B edge with (I, J), it joins to two parental edges only if a random integer U is between I and J. Then, one parental edge is labeled by (Min(I, M), Max(U, M)) and the other one by (Min(U + 1, M), Max(J, M)). If M ≤ U, the former parental edge must become a B edge, since the beneficial allele has descended from the ancestral sequence in this edge. The other parental edge, however, becomes either a B edge with probability x(t) or a b edge with probability 1 − x(t). Likewise, if M > U, the parental edge with (U + 1, Max(J, M)) becomes a B edge, with the other parental edge becoming either B or b. The same principle is applied to a recombination event in a b edge. The selective phase, which ends when x(t) ≤ ζ or the combined number of B and b edges becomes 1, is followed by another neutral phase where the distinction between B and b edges is erased.

The coalescent for each nucleotide site (or the “marginal tree”) is embedded in the ARG. The marginal tree is extracted as described in Griffiths and Marjoram (1997). The number of edges, k, at a given time changes stochastically during the construction of the ARG. Theoretically, k eventually hits 1, and the ARG is completed. However, if R »10, k fluctuates at high numbers, and waiting until k = 1 in the simulation is then practically impossible. We therefore stop the construction of the ARG at an arbitrary point, T_limit, which is chosen to be a large number relative to the mean time to the most recent common ancestor (MRCA) for a marginal tree. If a marginal tree cannot find the MRCA until T_limit, the remaining branches of the tree are forced to coalesce at T_limit. This procedure has a negligible effect on our results since most nucleotide sites find their MRCAs before T_limit. T_limit = 7.0 was used in this study. Source codes written in C for the simulation of ARGs and other procedures described in this article are available upon request.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Average nucleotide diversity along a recombining chromosome under the model of genetic hitchhiking, with n = 5, ρ = 10⁻⁸, N = 2 × 10⁵, φ = 0.01, and T_limit = 7.0. Squares represent average heterozygosity at single nucleotide sites averaged over 50,000 replicates of the simulations. The expected π value (continuous curve) was calculated using Equation 13 of Kim and Stephan (2000), with r = ρ |i − 20,000| as the recombination rate between a nucleotide site i and the site of selection. Directional selection occurs at position 20 kb with s = 0.001 and τ = 0.005.

PATTERNS OF GENETIC VARIATION ALONG A CHROMOSOME WITH HITCHHIKING

The simulated patterns of sequence polymorphism are obtained by introducing mutations into the marginal tree for each nucleotide site. To verify that the simulation procedure generates the correct ancestral genealogy expected under the model of hitchhiking, nucleotide diversities at many fixed sites along the sequence were summarized over 50,000 replicates of the ARG for a set of parameters (Figure 2). The simulation results agreed well with the expectation on the basis of the analytic solutions by Stephan et al. (1992) and Kim and Stephan (2000). Without the selective phase, the mean and variance of coalescent times of marginal trees agreed well with the expectation under the standard neutral model (data not shown). The number of shifts from one MRCA to another along the sequence was also observed and compared to the prediction (Equation 5 of Wiuf and Hein 1999). With R = 20 (2N = 10⁵, L = 10⁴, and ρ = 10⁻⁸), the observed numbers of shifts for sample sizes 2 and 10, each averaged over 200 replicates, were 13.72 and 19.74, respectively. The corresponding expectations are 13.33 and 19.64, respectively.

We assume that the derived allele can be distinguished from the ancestral allele, which is defined to be the allele at the root of the marginal tree. If more than one mutant is segregating at one site, all mutant alleles are classified as the derived allele and not distinguished from each other. To examine the pattern of variation, three different estimators [θ^π (Tajima 1983), θ^w (Watterson 1975), and θ^H (Fay and Wu 2000)] of φ = 4Nμ were calculated for the simulated sequences. Differences among the three estimators reveal deviations from neutrality (Tajima 1989; Fay and Wu 2000). Figure 3 shows the values of the three estimators for sliding windows along the sequence. Four replicates are shown for each set of parameter values, where each replicate was produced from an ARG of a 40-kb sequence. The ARGs were generated using the neutral model and the hitchhiking model with s = 0.001 (α = 100 or 1000) and τ = 0.001 or 0.2. As only four examples randomly chosen from the simulations are shown for each model, one may not be allowed to draw a general conclusion from these figures. However, some features of hitchhiking effects on sequence variation could be consistently identified from these examples. We use these examples mainly to illustrate these features.

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Patterns of genetic variation along a recombining chromosome. Simulated data were generated with (a) a neutral model with N = 5 × 10⁵ (R = 800); hitchhiking models with (b) N = 5 × 10⁵, s = 0.001, R = 800, α = 1000, and τ = 0.001; (c) N = 5 × 10⁵, s = 0.001, R = 800, α = 1000, and τ = 0.2; and (d) N = 5 × 10⁴, s = 0.001, R = 80, α = 100, and τ = 0.001. The values of the other parameters are n = 10, ρ = 10⁻⁸, φ = 0.005, and T_limit = 7.0. Selection occurs at position 20 kb. For each model, four replicates are shown. (black), (blue), and (red) were calculated along the sequence with window size 1 kb and step size 0.25 kb. Positions on the sequence are shown in units of 1 kb. *, ↓, and + indicate the positions of specific features described in the text.

A local reduction or valley of heterozygosity (θ^π ) along the sequence is the most important pattern expected under the model of genetic hitchhiking. Under neutral evolution (Figure 3a), the stochastic change of θ^π along the sequence occasionally generates deep valleys of variation (for example, regions indicated by *). However, in this case valleys are usually narrow compared to those under hitchhiking (Figure 3, b–d). The stochastic spatial pattern of variation is influenced by R. Using the same parameter values as in Figure 3a but smaller N, deeper and wider valleys were frequently observed (data not shown). With hitchhiking (Figure 3, b–d), a deep valley always appears at or around the site of directional selection. However, the “shape” of the valleys varies considerably among realizations for a given value of s. Valleys are rather asymmetrical around the site of the beneficial mutation, which implies that the shape of the valley may provide imprecise information about the location of the target of selection (see below). The asymmetry gets larger as N decreases. Figure 3d uses the same values of s, τ, and ρ as Figure 3b but a 10 times smaller N. As a result, the stochastic noise in the spatial pattern has been dramatically increased.

Compared to neutrality (Figure 3a), the relative level of θ^H versus θ^π increased immediately after the hitchhiking event (Figure 3b), as expected by Fay and Wu (2000). However, at τ = 0.2 (Figure 3c), i.e., 0.4N generations after the hitchhiking event, a higher relative level of θ^H as shown in Figure 3b is not observed. Especially, θ^H is distinctively lower than θ^π around the site of selection where the level of nucleotide diversity has only partially recovered since the selective sweep (regions labeled by ↓). This is consistent with the observation that the excess of high frequency variants appears suddenly after a hitchhiking event but soon disappears through the fixation of these alleles, reversing the excess of high frequency mutants (Kim and Stephan 2000). θ^w is expected to become larger than θ^π due to hitchhiking (Bravermanet al. 1995). In Figure 3, the increase of θ^w is not as obvious as in the case of θ^H . However, we could identify regions where θ^w became distinctively larger than θ^π in Figure 3b (labeled by +). The generality of these observations drawn from the examples of Figure 3 is further investigated by a statistical test applied to larger simulated datasets (see below).

The stronger the hitchhiking effect, the larger is the region that is expected to be affected. To find the relationship between the mean length of the region of reduced variation and the parameter values of the hitchhiking model, we generated 200 simulated datasets for a fixed combination of N, ρ, s, τ, and φ. A 1-kb-long window moves from the left end of the 40-kb-long sequence with an increment of one nucleotide and calculates θ^π at each position. Regions of reduced variation are defined by the centers of windows for which θ^π<θ∕2 . Therefore, a segment of the affected area is delimited by the centers of the two windows that mark the beginning and the end of the stretch of nucleotides with θ^π<θ∕2 . The length of the longest of such segments found on the 40-kb sequence is defined as W_φ/2. The mean and standard deviation of W_φ/2 over 200 replicates are shown in Table 1. The proportion, p_within, of the simulated datasets in which the site of the beneficial mutation is included in the largest segment that defines W_φ/2 is also recorded (Table 1). Examples 1–3, 7, and 8 show that an increase of the mutation rate per nucleotide, given by φ, does not lead to a proportional decrease in W_φ/2. Therefore, the mutation rates used in Table 1 are high enough to “saturate” and reveal the underlying stochastic patterns of coalescent times along the sequence. As expected, W_φ/2 with hitchhiking is significantly larger than that without. The mean of W_φ/2 is roughly proportional to s/ρ, as expected from the solutions of the hitchhiking effect (Maynard Smith and Haigh 1974; Stephanet al. 1992). As predicted by these solutions, population size N has also an effect on W_φ/2; W_φ/2 decreases with increasing N (examples 6 and 12). However, this effect is not as large as that determined by the parameter s/ρ, in particular for large values of N (Equation 19 in Stephanet al. 1992). Barton (2000) offered the following explanation of this effect of N. A beneficial mutation in a large population takes longer to reach a high frequency, by which time the lineages originally associated with it have become separated by recombination. This is equivalent to the statement that the “effective” length of the selective phase, −(2/s)log(ε) generations, is longer in large populations, where ε (≈1/(2Ns)) is the frequency of the beneficial mutation when it starts increasing deterministically. Looking backward in time, the rate of coalescence for two gene lineages during the selective phase gets sufficiently high only when the product of population size and beneficial allele frequency becomes low. Therefore, the waiting time (in generations) until the coalescence event is longer in large populations. However, the recombination rate per generation is independent of the population size. Therefore, the probability of the recombination event being the first event is higher in a larger population. This explains the smaller effect of a single hitchhiking event in a larger population as shown in Table 1.

Examples 6 and 7 show that almost identical W_φ/2's are obtained with the same values of Nρ and α, but with different ρ and s values. Therefore, for a given τ, Nρ and α are the two principal parameters governing the pattern of variation caused by a hitchhiking event. We compared the mean W_φ/2 to the theoretical prediction, E[W_φ/2] (Table 1), which is based on the expectation of θ^π along the sequence from Equation 13 of Kim and Stephan (2000). Mean W_φ/2 is consistently smaller than E[W_φ/2]. This discrepancy occurs partly because the calculation of E[W_φ/2] assumes that the duration of the selective phase is negligible on the timescale of 2N generations. However, the lengths of the selective phase, t_S = −(2/α)log(ζ), are 0.076 and 0.017 for examples 6 and 12, respectively, which are not much smaller than τ = 0.05. In the early part of the selective phase (when the frequency of the beneficial mutation is high), the behavior of the genealogy is similar to that in the neutral phase. Therefore, the length of the first neutral phase, i.e., the time since the last hitchhiking event, is effectively longer than τ. It should be noted that the standard deviation of W_φ/2 is considerably large, as suggested by Figure 3. Furthermore, in >20 of 200 realizations, the site of the beneficial mutation is not included in the largest segment of reduced variation (see p_within in Table 1) even with strong hitchhiking (examples 11 and 12). These results again indicate a large amount of stochasticity in the pattern of variation shaped by hitchhiking effects.

After a selective sweep, the level of genetic variation is slowly restored due to new neutral mutations. Therefore, with given values of Nρ and α, W_φ/2 should become smaller with increasing τ, as is indeed observed in examples 11–13, for which τ = 0.001, 0.05, and 0.2, respectively. The level of variation around the site of selection, which is zero immediately after the sweep, should be characterized by τ (Wiehe and Stephan 1993). As the site of selection is expected to be found at the center of the largest segment of reduced variation (defined above), we examined the average nucleotide diversity, θ^π∗ , for the middle one-third of this segment. The average values (±standard deviation) of θ^π∗ for examples 11–13 are 3.88 × 10⁻⁴ (±4.09 × 10⁻⁴), 5.48 × 10⁻⁴ (± 3.85 × 10⁻⁴), and 9.85 × 10⁻⁴ (±3.86 × 10⁻⁴), respectively. The corresponding theoretical values obtained by numerical integration of Equation 13 of Kim and Stephan (2000) are 5.45 × 10⁻⁴, 7.24 × 10⁻⁴, and 1.23 × 10⁻³, respectively. Therefore, the “depth” of the valley of reduced variation contains information about the time of the last hitchhiking event. However, the large standard deviations of θ^π∗ suggest that a correct estimation of τ from polymorphism data will be very difficult.

So far the patterns of genetic variation based on the segregation at single sites were examined. Another important aspect of sequence variation is the association of polymorphisms between neighboring loci. We calculated r² (Hill and Robertson 1968) for all pairs of segregating sites in each simulated set and obtained the average, r2¯ , over the entire 40-kb region. Hitchhiking caused an increase in r2¯ (Table 1). It might be possible that this increase in r2¯ was caused by the excess of rare alleles (i.e., singletons) generated by the hitchhiking effect, because r² frequently becomes large by chance when the allele frequencies at both loci are extreme. Unfortunately, we could not exclude singletons from the analysis since not many segregating sites are left if singletons are removed from the data generated under the hitchhiking model. However, a visual inspection of the raw hitchhiking data reveals that there are several extensive haplotype structures that are not likely to be created by chance alone. Large stretches of polymorphic sites share an identical pattern of segregation; i.e., there are only two haplotypes observed in such a stretch. We recorded the maximum number, S_max, of such consecutive sites found in each dataset. Table 1 shows that S_max increases with hitchhiking. With strong selection, the increase of S_max is very large (for instance, compare examples 5 and 11). The increase of r2¯ and S_max by hitchhiking can be explained by a coalescent argument. Ancestral histories of two neutral loci become identical if no recombination event occurs before the MRCA for both loci is found (Griffiths and Marjoram 1997; Wiuf and Hein 1999). During the selective phase the rate of coalescence is increased while that of recombination remains the same. Therefore, hitchhiking removes opportunities for recombination between two linked loci. This leads to an increase in the correlation of gene genealogies between segregating sites around the site of directional selection. However, the buildup of linkage disequilibrium due to hitchhiking quickly decays as τ increases (examples 11–13), because recombination events during that neutral phase break up associations created in the selective phase.

STATISTICAL TEST OF A LOCAL SIGNATURE CAUSED BY GENETIC HITCHHIKING

In the following, a maximum-likelihood method is developed to examine the significance of a local reduction of genetic variation and to estimate the strength of directional selection. The probability of observing a certain frequency of derived alleles at a site after a recent hitchhiking event can be obtained by previously used analytic approximations. Under neutrality, the expected number of sites where the derived variant is in the frequency interval [p, p + dp] in the population is given by ϕ0(p)dp=θpdp (3) (Kimura 1971). Immediately after a hitchhiking event, this distribution is transformed approximately to ϕ1(p)={θp−θCfor0<p<CθCfor1−C<p<1} (4) (Fay and Wu 2000), where C is given approximately by 1 − ε^r/s (appendix). Here, r is the recombination fraction between the neutral locus and the selected locus and ε is the frequency of the beneficial allele when it begins to increase deterministically. It should be noted that Equation 4 is obtained by assuming deterministic changes of allele frequencies during the selective phase. The probability of observing a site where k derived alleles are found in a sample of size n is given by Pn,k=∫01(nk)pk(1−p)n−kϕ(p)dp(k=1,…,n−1) (5) and Pn,0=1−(Pn,1+…+Pn,n−1), where φ(·) = φ₀(·) under the neutral model, and φ(·) = φ₁(·) under the hitchhiking model. P_n,k was found to be sensitive to the choice of ε. We used ε = 1/α, which gave the best fit to the simulation results (Figure 4). The likelihood of all data under the model of genetic hitchhiking is obtained by multiplying the probabilities for all nucleotide sites under consideration. This is a composite likelihood because there is a correlation of P_n,k between sites due to shared ancestral histories. Therefore, it should be distinguished from the conventional likelihood-ratio test that is based on exact likelihoods. A statistical test in this analysis thus depends on an empirical distribution of the test statistic obtained by simulation. Composite likelihood is frequently used when the derivation of exact likelihoods is difficult (e.g., Rannala and Slatkin 2000). It should also be noted that higher-order structures in the polymorphism data, such as linkage disequilibria, are neglected in this analysis.

As it is currently unrealistic to have polymorphic data from a reasonably large sample of long continuous sequences, we apply the test to a region for which only short segments are sequenced, interspaced with larger nucleotide stretches for which no data are available. That is, we consider a survey in which 11 1-kb-long segments distributed over a 40-kb region are sequenced. The distances between segments are uniformly 2.9 kb. The sample size is 10 for all segments. Simulated data used for Table 1 (examples 9 and 11–13) were reused, but only sites from the 11 segments were included. The maximum composite likelihood under the neutral model (L₀) and that under the hitchhiking model (L₁) were obtained for each simulated dataset. Then, the likelihood ratio is given by L₁/L₀. L₀ is a function of φ and L₁ is a function of N, μ, ρ, s, and the location of the selected locus, X. X is allowed to vary in the middle 10-kb region of the sequence (15 kb < X < 25 kb); i.e., we consider a situation where a candidate region for the site of selection has already been inferred. It is difficult practically to allow all these parameters to vary freely until a unique combination that maximizes L₁ is found. Therefore, we chose only s and X as free variables and assumed that separate estimates of N, μ, and ρ are available. Thus, in one test (test A), the same values of N, μ, and ρ specified in the simulation are used in the calculation of L₀ and L₁. In the other test (test B), to be conservative, we let the mutation rate be inferred from the data by using the average heterozygosity (θ^π ) over all 11 segments of the sequence as the fixed prior estimate of φ. Therefore, the standing level of variation is simply the level observed in the data. But we still used the true value of N for the calculation of α = 2Ns. We also assumed either that the derived neutral allele is distinguished from the ancestral allele at each site (option 1) or that they are not distinguished (option 2). For the latter, there are only five ratios of segregating variants in the sample of 10 sequences. Let Q_k,n (k = 1, … , n/2) be the probability of observing a [k(n − k) segregation ratio. Then the likelihood ratios are calculated simply by using Q_n,k = P_n,k + P_n,n−k.

The null distribution of likelihood ratios was obtained by applying tests to datasets generated under the neutral model (200 replicates corresponding to example 9 of Table 1). To reduce the problem of local optima, eight different initial guesses of X between 15 and 25 kb, with s = 0.01, were used to start the maximization procedure using Powell's method (Presset al. 1992; program written by D. L. Swofford and kindly provided by J. P. Huelsenbeck). The null distribution of log(L₁/L₀) for test A/option 1 is shown in Figure 5. All test methods (tests A and B and options 1 and 2) produced almost identical null distributions (data not shown). About 30% of the log(L₁/L₀) values are negative. This may suggest a failure of obtaining the global maximum of L₁ for some null datasets, since L₁, with two more free variables, should be larger than L₀. However, increasing the number of initial guesses did not improve the likelihood ratios (data not shown). The negative values are more likely due to the restricted range of s and the assumption of τ = 0. As α should be large enough to generate a hitchhiking effect (also ε = 1/α should remain small), s was not allowed to be <6 × 10⁻⁵ (α = 60). Only when s is very small, the hitchhiking model with τ = 0 can fit neutral data in which a local reduction of variation is not found. Therefore, there was a limit in maximizing L₁.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Frequency of derived alleles at individual nucleotide sites in a sample of five sequences under genetic hitchhiking. Three nucleotide sites were chosen such that r/s = 0.01, 0.1, and 0.2. For each site, the numbers of occurrences of one to four derived alleles were recorded out of 50,000 replicates. The values of the other parameters are φ = 0.01, N = 2 × 10⁵, and s = 0.001. Theoretical expectations were obtained using Equation 5, with ε = 1/α.

Table 2 summarizes the power of the test and the point estimates of s and X for each hitchhiking model. Power is the proportion of replicates that produce log(L₁/L₀) values greater than the 95th percentile of the corresponding null distribution. Test A yielded very high power of rejecting neutral evolution, even for larger values of τ. The main reason for obtaining large likelihood ratios from test A is that it uses the “true” standing level of variation (φ) for the calculation of L₁ and L₀. As the average heterozygosity has been reduced below φ due to selective sweep, the neutral model based on the true value of φ cannot fit the data. The negligible differences in the power between options 1 and 2 (except at τ = 0.2) indicate that the additional information obtained by distinguishing between ancestral and derived alleles contributed little in test A, whereas a significant reduction of heterozygosity played a major role in increasing the likelihood ratio.

On the other hand, a reduction of heterozygosity is not a major factor for increasing the likelihood ratio in test B. To obtain a higher likelihood ratio in this test for a given number of segregating sites, the spatial distribution of those sites along the sequence and the allele frequency spectrum should be close to the expectation under the hitchhiking model. As expected, the power of test B is smaller than that of test A (Table 2). However, it is still high (84–97%) for small values of τ (0.001 and 0.05). Power declines as τ increases, since the spatial pattern and the frequency spectrum of segregating sites approach those under neutrality as time passes after the selective sweep. Tests using option 1 yield higher power than those using option 2 for τ = 0.001 and 0.05, which means that the skew toward high-frequency-derived alleles at segregating sites is observed as described by P_n,k (Figure 4). For τ = 0.2, however, both tests A and B had higher power with option 2. This is obvious from the fact that, at τ = 0.2, the proportion of high-frequency-derived alleles is lowered below its level under neutrality (Figure 3c). Therefore, distinguishing the derived allele from the ancestral allele in these tests has an advantage for detecting very recent hitchhiking events only. However, it should be noted that our analytic prediction of the frequency spectrum is based on the assumption of τ = 0. Complete solutions for P_n,k for any value of τ may make option 1 still useful for detecting more distant hitchhiking events.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Distribution of log(L₁/L₀) obtained by applying test A/option 1 to neutral datasets (R = 800).

Maximum (composite)-likelihood estimates of s and X were also obtained. Test A with option 1 produced the most unbiased estimates of s, although the accuracy is quite low for all combinations of the test methods (Table 2). Joint estimates of s and X using test A with option 1 for datasets generated under neutrality and under hitchhiking with s = τ = 0.001 are shown in Figure 6. From the neutral data, joint estimates were clustered in the parameter space of small s (close to the lower limit) and X between the “sequenced” segments. This can be expected since the hypothesized valley due to hitchhiking should be sufficiently narrow to fit between the sequenced segments where the level of polymorphism is high. On the other hand, the joint estimates from the hitchhiking datasets were centered around the true value (s = 0.001, X = 20 kb). It is also shown that estimates of X tend to cluster on the sequenced segments in this case.

To further investigate the performance of the composite likelihood-ratio test, we produced additional but shorter (10-kb) sequences by simulation. Unlike in the previous analysis, polymorphism data are assumed to be obtained from the entire continuous region and X is also allowed to vary over the entire region. Only option 1 is used. To make results comparable to the previous cases, R was adjusted to be 800 by setting ρ = 4 × 10⁻⁸. In the simulation of selective sweeps, selection occurs at position 5 kb with τ = 0.001, but with various s values (Table 3). With s = 0.001 (α = 1000), the powers of tests A and B were 0.97 and 0.915, respectively. These can be compared with 0.995 and 0.97 from the previous analysis (Table 2, s = 0.001 and τ = 0.001). Considering a slight reduction of the surveyed region (11–10 kb) where informative segregating sites were observed, the power of the tests with this new scheme appears to remain as high as in the previous one using discontinuous regions. With decreasing s, the power of detecting the hitchhiking event and the accuracy of the parameter estimates decrease, as expected (Table 3). However, power declined also when s increased from 0.002 to 0.005. Examination of simulated data showed that, with s = 0.005, the number of segregating sites is highly reduced and those sites are frequently found clustered on one side of the site of selection. This produced very low likelihood ratios in both tests A and B. This effect should disappear if a wider region of the chromosome is surveyed.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

Joint estimates of s and X using test A and option 1 for simulated neutral and hitchhiking (s = 0.001, i.e., α = 1000; τ = 0.001) datasets. The substitution of the beneficial allele occurs at position 20 kb. Solid squares represent estimates obtained from hitchhiking datasets and open circles represent those from neutral datasets. A total of 200 replicates were used for each model. Shaded segments over the y-axis represent regions from which polymorphism data were obtained.

We also conducted a few tests to assess the effect of uncertainty in prior estimates of N and ρ. Tests A and B were performed for a dataset described above (Table 3, s = 0.001) but with a prior estimate of N = 10⁵, fivefold lower than the true value. New null distributions of the likelihood ratio were obtained accordingly. The power of tests A and B decreased to 0.935 and 0.865, respectively, from 0.97 and 0.915 (Table 3) due to the incorrect assumption of N. Average ŝ decreased slightly (8.5 × 10⁻⁴ from 1.15 × 10⁻³ for test A and 5.1 × 10⁻⁴ from 7.1 × 10⁻⁴ for test B). Next, we used the correct value of N but a fivefold lower prior estimate of ρ (8 × 10⁻⁹). Average ŝ decreased about fivefold (1.78 × 10⁻⁴ and 1.19 × 10⁻⁴ for tests A and B, respectively) as expected. Power decreased to 0.855 and 0.78 for tests A and B, respectively.