Abstract
The nonparametric bootstrap approach is known to be suitable for calculating central confidence intervals for the locations of quantitative trait loci (QTL). However, the distribution of the bootstrap QTL position estimates along the chromosome is peaked at the positions of the markers and is not tailed equally. This results in conservativeness and large width of the confidence intervals. In this study three modified methods are proposed to calculate nonparametric bootstrap confidence intervals for QTL locations, which compute noncentral confidence intervals (uncorrected method I), correct for the impact of the markers (weighted method I), or both (weighted method II). Noncentral confidence intervals were computed with an analog of the highest posterior density method. The correction for the markers is based on the distribution of QTL estimates along the chromosome when the QTL is not linked with any marker, and it can be obtained with a permutation approach. In a simulation study the three methods were compared with the original bootstrap method. The results showed that it is useful, first, to compute noncentral confidence intervals and, second, to correct the bootstrap distribution of the QTL estimates for the impact of the markers. The weighted method II, combining these two properties, produced the shortest and less biased confidence intervals in a large number of simulated configurations.
SEVERAL statistical methods were proposed for detecting quantitative trait loci (QTL) using marker information. Most of them are based either on maximum-likelihood procedures (Lander and Botstein 1989) or on regression methods (Haley and Knott 1992; Martinez and Curnow 1992). However, none of these methods led to a straightforward calculation of a confidence interval. These are intervals on the chromosomes that contain the QTL with a given coverage probability P (usually 90 or 95%). Calculating small confidence intervals that hold the stated coverage probability is of fundamental importance for practical breeding purposes as well as for molecular biologists. For example, breeders who wish to use the marker information in marker-assisted selection (MAS) breeding schemes need to know the optimum length of the chromosome segment of interest, and biologists need this knowledge for any further decisions regarding experimental fine-mapping strategy.
Different methods were described to calculate confidence intervals in QTL mapping. Conneally et al. (1985) and Lander and Botstein (1989) proposed the LOD drop-off method. This is based on the conversion of the distribution of the likelihood-ratio test statistic of linkage vs. no linkage between a marker and a putative QTL into a chi-square distribution. Using this method, the endpoints of the confidence interval are calculated by finding the locations at either side of the estimated QTL position having the maximum LOD score that have a LOD score of 1 unit less than the estimated QTL position. The confidence interval is then defined by all those points on the chromosome where the LOD score has fallen within 1 unit from the maximum (i.e., the positions between the two interval endpoints). However, this method often produces intervals that are too short; i.e., the intervals contain the QTL with a lower coverage probability than stated. Van Ooijen (1992) and Mangin et al. (1994) showed this to be the case for small and medium-sized populations and for QTL of small effects because the regular conditions for the conversion of the likelihood-ratio test toward a chi-square distribution do not hold in these genetic configurations. To attain a higher probability of containing the QTL, Van Ooijen (1992) suggested that the LOD drop-off confidence intervals should be based on two LOD differences. In a simulation study it was found that the sizes of these confidence intervals were large and variable (Van Ooijen 1992) and, in population sizes used in real QTL mapping procedures, unsatisfactory for QTL with smaller effects.
Mangin et al. (1994) and Mangin and Goffinet (1997) developed complex formulas for calculating confidence intervals. These are based on a maximum-likelihood-ratio test using statistics in which asymptotic distribution does not depend on the effect of the QTL. The major difficulty of this method is the computation of the correct threshold for the likelihood-ratio test for which Mangin and Goffinet (1997) proposed an approximate analytical solution.
Visscher et al. (1996) introduced a nonparametric bootstrap resampling approach (Efron 1979; Efron and Tibshirani 1993) to construct a confidence interval for a QTL position. The evaluation of generated bootstrap samples provided a distribution of QTL position estimates along the chromosome that forms the base for the confidence interval calculation. Within a simulation study they compared the nonparametric bootstrap method with the LOD drop-off method. The nonparametric bootstrap procedure seemed a good alternative to the LOD drop-off for defining confidence intervals for QTL positions, but tended to be unsuitably large when the QTL effect was small.
To reduce the width of the estimated bootstrap confidence interval a selective nonparametric bootstrap method was proposed by Lebreton and Visscher (1998). Their selection was based on criteria related to the assumed genetic model. By comparing both methods in a simulation study they found that the selective bootstrap method produced smaller confidence intervals with a reduced conservativeness. Moreover, they showed that there is scope for improvement of the classical bootstrap approach for calculating confidence intervals as proposed by Visscher et al. (1996).
Walling et al. (1998) provided a parametric bootstrap resimulation method (Efron 1982) as an alternative to the nonparametric bootstrap approach and tested this method in a simulation study. As the authors themselves pointed out, the parametric bootstrap method produced poor results: The derived confidence intervals varied from conservative intervals when the QTL position was at the location of the marker to anticonservative intervals when the QTL position was between two markers.
A main result from the study of Walling et al. (1998) was the discovery of a bias in calculating confidence intervals from the distribution of QTL position estimates from nonparametric bootstrap samples. This bias is defined as a significant difference between the observed and the stated coverage probability and is deemed to have been caused by the locations of the markers. Eliminating this bias would probably result in smaller confidence intervals. As mentioned by Lebreton and Visscher (1998), an additional reduction of the interval width could be possible by computing noncentral confidence intervals. Taking this into consideration, we present three modified nonparametric bootstrap schemes to calculate confidence intervals: (i) noncentral, (ii) marker corrected, and (iii) noncentral and marker corrected. Within a simulation study these schemes were compared with the nonparametric bootstrap method presented by Visscher et al. (1996).
MATERIALS AND METHODS
Simulation: A first generation backcross population consisting of one half-sib family and derived from two inbred lines was generated and contained 300 or 500 (N = 300, 500) progeny. For each individual a single 100-cM chromosome (L = 100 cM) evenly spaced with markers every 20 cM (Δ = 20 cM) was simulated. For calculating the crossover probabilities the Kosambi mapping function was used. Each chromosome contained a single segregating diallelic QTL with a position midway between two flanking markers (d = 30, 50 cM), closer to one of the flanking markers (d = 25, 55 cM), or directly at a marker position (d = 20, 60 cM). The effect of the QTL was expressed in units of the polygenic standard deviation, and to cover a wide range of effects it was varied between 0.1, 0.5, and 1.0 (a = 0.1, 0.5, 1.0). Additional simulations were carried out to investigate the impact of a longer chromosome (L = 140 and 180 cM, Δ = 20 cM, N = 300, 500 with a = 0.5 and d in the middle of the chromosome) and a higher marker density [Δ = 10 cM and 10/5 cM (additional markers at 35, 45, 55, and 65 cM), L = 100 cM, a = 0.5 with N = 300, 500, respectively]. To avoid any confounding effects due to different detection power (Darvasiet al. 1993) the QTL position was located exactly midway between two flanking markers when marker spacing was reduced (d = 55 cM when Δ = 10 cM and d = 52.5 cM when Δ = 10/5 cM).
The total genetic variance was the sum of the polygenic variance and the variance due to the QTL. The polygenic genetic component was normally distributed. The heritability of the trait, defined as the polygenic variance divided by the sum of polygenic variance and the nongenetic variance, was assumed to be 0.25. Random nongenetic factors were normally distributed. For each simulated genetic configuration 1000 replicates were generated. The simulation program was written in FORTRAN 77, supplemented with routines from the NAG library (Numerical Algorithms Group 1990). QTL analyses were done by using the programs BIGMAP and ADRQLT (briefly described in Reinsch 1999) and were based on the multimarker regression described by Knott et al. (1994). Confidence intervals were calculated with FORTRAN 77 programs and SAS procedures (SAS Institute 1992).
Generating the bootstrap samples: A total of 250 nonparametric bootstrap samples for each genetic configuration were generated as described by Visscher et al. (1996). Briefly, for each single bootstrap sample N observations out of the pool of N original observations were drawn with replacement, so that each observation consisted of an individual's marker genotype and the corresponding phenotype. Some observations may appear more than one time in a bootstrap sample while others may be not included at all. Each single bootstrap sample was analyzed and the best estimate of the position of the QTL was recorded. After 250 bootstrap samples the distribution of the 250 estimates of the QTL position was derived by ordering them along the chromosome. Because the QTL is linked to two flanking markers, this distribution is termed the linkage distribution.
Modeling marker impact on bootstrap distributions: As reported by Walling et al. (1998) the reasons for the marker impact on the distribution of the bootstrap estimates for the QTL position are an accumulation of the type I error rate at the marker loci and the general trend of the mapping procedure to put the estimate of a QTL position toward a marker locus. In the following section a closer look at this marker impact is presented by a simple model, and the theoretical background for a marker correction is given. At first we focus on the expected distribution of the best estimates for a QTL position out of a bootstrap sample when a QTL may be present but not linked with any marker. This distribution is termed the null distribution, because every estimate indicates a type I error. In this null distribution the probability of each position achieving the best estimate would be
Weighted method I: In the weighted method the marker correction approach was followed. First, the probability in the linkage distribution to attain the best estimate P(i)ld and the selection coefficient S(i) had to be determined for each position on the chromosome. The values for P(i)ld were the corresponding observed frequencies of best QTL estimates in the linkage distribution [f(i)ld]:
The confidence intervals with a coverage probability P of 95 and 90% were calculated by using the top and the bottom 2.5th and 5th percentiles as the upper and lower interval endpoints, respectively. The obtained confidence intervals were central.
Weighted method II: For the weighted method II the corrected distribution was computed in the same manner as for the weighted method I. Whereas for the weighted method I one-half of the given error rate of 1 − P was subtracted from each side of the corrected distribution during the computation of the confidence intervals, in the weighted method II an analog highest posterior density (HPD) method was used to obtain noncentral intervals when the corrected distribution is not tailed equally. This method was adapted from the standard HPD approach as described in Box and Tiao (1992) and is presented in detail in the appendix.
Uncorrected method I: This method was the original one as proposed by Visscher et al. (1996) and did not correct for the marker impact. The central confidence intervals were calculated directly from the linkage distribution by taking the top and bottom 2.5th and 5th percentiles of the linkage distribution as the upper and lower interval endpoints, respectively.
Uncorrected method II: In this method confidence intervals were computed from the linkage distribution with the analog HPD method described in the appendix. The calculated intervals were noncentral and not marker corrected.
The confidence intervals with the stated coverage probability P of 95 and 90% were calculated for each replicate by the four methods described above, and it was assessed whether the QTL was located within these intervals or not. The rate of the confidence intervals for each simulated genetic configuration that did not include the simulated position of the QTL was termed the noninclusion rate (NI), and this should be not greater than 1 − P. For each configuration the average width of the confidence interval out of 1000 replicates was also calculated, and this is ideally small.
RESULTS
The initial results for all four bootstrap schemes are presented in Tables 1, 2, 3 and 4. Summarized over all genetic configurations analyzed, the uncorrected method I produced the largest and most conservative confidence intervals. A reduction of the interval width and conservativeness was possible using the uncorrected method II. Further reduction of the width was achieved using the weighted method I and, even more, by the weighted method II without lifting the noninclusion rate seriously above 1 − P.
QTL effect: As expected, an increase of the QTL effect resulted in a smaller confidence interval for all bootstrap schemes (Table 1 and 2). It is very unlikely that a QTL with an effect of a = 0.1 (this QTL would explain only 0.13% of the total variance in our simulated populations) would be detected in real mapping experiments and, hence, the uncorrected method I produced very conservative confidence intervals, which cover almost the whole chromosome. For this low QTL effect each confidence interval produced by this method was >93 cM. However, the uncorrected method II (weighted method I, weighted method II) produced 90% intervals for this small QTL effect that were significantly <90 cM (85 cM, 80 cM) for the 90% confidence intervals and were less conservative. All confidence intervals for situations of a = 0.1 were only marginally influenced by the remaining genetic parameters.
Elevating the QTL effect from a = 0.1 to 0.5 (this QTL would explain 3.2% of the total variance) reduced the average 90% interval width by 15–25 cM. The respective change by stepping 0.5 to 1.0 (11.8% of the total variance) was 40–50 cM, depending on the population size and QTL position. This reduction was remarkably constant for the four methods.
Effect of QTL position within marker interval, population size, and QTL effect on the confidence intervals for the different bootstrap methods—QTL position proximal
Marker spacing: Changing marker distances from 20 to 10 cM (top of Table 2 and Table 3) and adding four additional markers around the QTL position (Table 3) resulted in a slightly reduced interval width for all bootstrap schemes. The reduction for the 90% confidence intervals varied between 2 and 4 cM. The corresponding changes were between 1 and 5 cM when marker distance was reduced to 10 cM and when additional markers were added around the QTL position, depending on the population size and the bootstrap scheme. No general dependency of the noninclusion rate on the marker density was observed for any method.
Population size: Tables 1, 2, 3 and 4 show that the elevation of the population size from N = 300 to N = 500 led to significantly smaller confidence intervals for all four bootstrap schemes. This reduction was only marginal for a = 0.1 but for a = 0.5 it was between 12 and 14 cM and for a = 1.0 it was between 10 and 15 cM. The uncorrected method I was most sensitive for increased population size as the reduction of the interval width was greatest. This might be due to the fact that intervals calculated with this method were largest and consequently there was the most space for reduction. It was not possible to assign a change of noninclusion rates to the effect of an increased population size for the uncorrected method II and for the weighted methods I and II; it appears that the intervals computed with the uncorrected method I were in general less conservative when the population size was increased.
Effect of QTL position within marker interval, population size, and QTL effect on the confidence intervals for the different bootstrap methods—QTL position central
In general, a gain of information, provided by a greater population size, resulted in smaller intervals. This is consistent with other studies concerning QTL bootstrap confidence intervals (Visscheret al. 1996; Lebreton and Visscher 1998; Wallinget al. 1998). The effect of the larger population size was stronger than the effect of the higher marker density.
Effect of an increased marker density on the confidence intervals for the different bootstrap methods
Length of chromosome: The length of the chromosome had a large impact on the size of the confidence intervals and a small impact on the noninclusion rates for all bootstrap schemes; i.e., the increase of the chromosome length resulted in larger intervals and, in most cases, in slightly lower noninclusion rates (Table 4). The absolute difference between the length of the intervals produced by the four methods increased slightly with a longer chromosome without changing the ranking order.
QTL position: The effect of the QTL position can be attributed to two factors: the QTL position within a marker interval (midway between two markers, closer to one marker, or at the same position as the marker) and the chromosomal QTL position, whether proximal or central. To avoid any confounding effects these two factors were analyzed separately.
Effect of the length of the chromosome on the confidence intervals for the different boostrap methods
It seems that when the QTL is located centrally (position 60 instead of 20, 55 instead of 25, and 50 instead of 30; Tables 1 and 2) the computed intervals are smaller, regardless of the method. This is in close agreement with Lebreton and Visscher (1998) and Walling et al. (1998) and may be caused by the fact that the mapping procedure tends to estimate the QTL position closer to the middle of the chromosome, as reported by Hyne et al. (1995). Additionally, the noninclusion rates were in general lower for a = 0.1 and a = 0.5 but higher for a = 1.0 for all methods.
Changing the QTL position within the marker interval from midway between two flanking markers toward a marker position (top, middle, and bottom of Tables 1 and 2) resulted in significantly smaller confidence intervals for all four methods. For example, the reduction for the 90% interval for a = 0.5 and N = 500 was between 1 and 2 cM when d was changed from 50 to 55 cM and ~4 cM when d was changed from 55 to 60 cM, regardless of the method. The uncorrected methods I and II produced intervals that were more conservative when the QTL and the marker positions are identical compared to situations with QTL positions midway between two flanking markers. This is consistent with the findings of Walling et al. (1998). The noninclusions were in general increased for the weighted method I when d was moved toward a marker position whereas no general dependency could be observed between the QTL position and the noninclusion rates for the weighted method II.
Another remarkable topic might be the differences between the interval widths from the uncorrected methods I and II on the one hand and the differences between the weighted methods I and II on the other hand when moving the QTL closer to the start of the chromosome. To investigate the impact of the QTL position on these differences in detail, additional simulations with N = 500, a = 0.5, Δ = 20 cM, L = 100 cM, and a QTL position at d = 10 or 15 cM were carried out. The average 90% interval widths computed with the uncorrected method II (weighted method II) were subtracted from the corresponding intervals computed with the uncorrected method I (weighted method I) for the genetic configurations described above and the same configurations but with d = 20 and 30 cM (Table 1) and d = 50 cM (Table 2). Additionally, the ratio of the subtracted upper tail and the subtracted lower tail of the linkage distribution and the corrected distribution during the computation of the 90% confidence interval with the analog HPD method were calculated. The results show that both the differences between the central and noncentral confidence intervals and the ratio of the subtracted upper and lower tail from the distributions increased when the QTL was located closer to the start of the chromosome (Table 5). The increase of the ratio was greater for the intervals computed from the corrected distribution than for those computed from the linkage distribution. When the QTL was midway on the chromosome (d = 50 cM) the ratio became, as expected, exactly one or near to one. It is important to note that the computed intervals for d = 10 cM were slightly anticonservative (10.4, 12.5, 12.9, and 13.9% noninclusions for the 90% intervals computed with uncorrected methods I and II and the weighted methods I and II, respectively; results not shown).
Effect of the QTL position on the differences between the widths of central and noncentral confidence intervals
Selection step: In practical QTL mapping experiments confidence intervals would be computed only when the QTL is declared significant. To investigate whether the rejection of nonsignificant replicates would change the results a selection step for two genetic configurations similar to Visscher et al. (1996) was introduced. First, according to Churchill and Doerge (1994), threshold levels that gave a nominal type I error rate of 10 and 5% were calculated by taking the 90th and 95th quantiles of the test statistic of 10,000 analyzed simulated populations with no segregating QTL. Confidence intervals were computed only for replicates that showed QTL test statistics above the threshold. The results show that all computed confidence intervals were smaller and less conservative when introducing a selection step (top of Table 2 and Table 6) and when reducing the nominal type I error rate (Table 6). This finding is in agreement with Visscher et al. (1996). Generally, the change of the nominal error rate from 10 to 5% reduced the interval width between 1 and 3 cM and increased the noninclusion rate by ~1–2%. No differences in sensitivity for the selection step between the four bootstrap methods could be observed.
Effect of rejection of nonsignificant replicates on the confidence intervals for the different bootstrap methods
Number of bootstrap samples: In computing confidence intervals from the linkage distribution Visscher et al. (1996) showed that there is only a small impact on the number of bootstrap samples. To examine whether this holds also for the computation of intervals from the corrected distribution, the number of bootstrap samples was varied for two genetic configurations (Δ = 20 cM, a = 0.5, L = 100, d = 50 cM, and N = 300, 500, respectively). It seems that there is a lower bound for the number of bootstrap samples to generate the corrected distribution for interval computation that is between 100 and 200 and dropping it below 100 results in smaller intervals with a higher noninclusion rate. Increasing the number of bootstrap samples above 200 did not change the results (Table 7). Following this, the number of bootstrap samples (250) used in this study is appropriate.
DISCUSSION
The classical bootstrap method (in this study the uncorrected method I) for calculating confidence intervals for a QTL position as proposed by Visscher et al. (1996) seems a suitable and practical approach, and in some recent real QTL mapping projects this method was used to obtain confidence intervals for the QTL position (e.g., De Koninget al. 1998; Zhanget al. 1998; Wallinget al. 2000). However, this method produces intervals that tend to be rather large and conservative, and therefore we proposed three methods (uncorrected method II and weighted methods I and II) to improve this classical approach to calculate smaller and less conservative confidence intervals and compared these methods within a simulation study with the classical approach.
Bias of confidence intervals: The bias of confidence intervals for QTL position estimates (biasci) can be defined as where NI is the noninclusion rate. The biasci for the intervals computed with the uncorrected method I was in general greatest and was upward except when the QTL was at the start of the chromosome and significantly increased for this method when the QTL position was located at a marker position. This is in agreement with the findings of Visscher et al. (1996) and Walling et al. (1998). When comparing the uncorrected method I with the weighted method I and the uncorrected method II with the weighted method II, respectively, it follows that the marker correction of the linkage distribution leads in general to a reduced bias of the confidence intervals as the noninclusion rate becomes closer to 1 − P. Note that this bias was for the weighted method I, and even more, for the weighted method II, slightly downward, but not substantial, when the QTL effect was small and the QTL was located at the start of the chromosome but upward for the remaining configurations. Together with the fact that the weighted methods I and II produced significantly smaller confidence intervals than the uncorrected methods I and II it follows that the described model of characterizing the marker impact (Equations 1, 2, 3 and 4) and the marker correction in the weighted methods I and II is appropriate to obtain smaller and less biased bootstrap confidence intervals for the location of QTL that appear to be significant. Note that the probability of making a type I error is a function of the chosen significance level, which is not affected by any of the presented methods to compute confidence intervals.
Effect of number of bootstrap samples on the confidence intervals for the weighted methods I and II
Marker correction: Although it may be possible to calculate the influence of the markers on the linkage distribution, in terms of the selection coefficient for each possible QTL position (each centimorgan in our study), in an analytical way, in this study it was done empirically by computing the null distribution for each genetic configuration, i.e., the distribution of the best QTL position estimates when there is no cosegregation between the QTL and any marker, and, hence, every estimate indicates a type I error per definition. This was possible because in the model it is assumed that the selection coefficient is the same for the linkage distribution as for the null distribution. The selection coefficient for each position was the corresponding frequency of the best estimate in the null distribution (Equation 6). In Figure 1A a null distribution for population size of N = 500, marker space of 20 cM, chromosome length of L = 100 cM, and QTL effect of a = 0.5 is presented. When the QTL is not linked with any marker, the probability for each position on the chromosome to achieve the best estimate is significantly higher at the marker loci. Over 45% of the estimates were placed there (under an equal distribution it would be only 5.94%). With an increased marker density (10 cM) 64% of the estimates would be placed at the marker loci. Except for the marker density and chromosome length it seems that the null distribution is of little variability, since varying the QTL effect and the population size resulted in almost the same null distribution (not shown).
The extent of the influence of the markers on the linkage distribution depends on the level of the QTL effect as it is shown in (3). Figure 1B shows an average linkage distribution from all bootstrap samples and all replicates for N = 500, L = 100 cM, Δ = 20 cM, a = 0.5, and for QTL position of d = 30 cM from the start of the chromosome. When the QTL effect was a = 0.1 (a = 0.5, a = 1.0), 53% (32%, 9%) of the best estimates were located at the marker loci. The same genetic configurations but with a marker density of 10 cM gave roughly 10% higher values for the hits per marker locus, regardless of the QTL effect (results not shown), and this confirmed the presence of the marker influence on the linkage distribution as first reported by Walling et al. (1998). The two higher peaks at 20 and 40 cM are a result of the high value of the product of PQTL and selection coefficients at these positions (Equation 3). Both are relatively high because the distance to the position with the simulated QTL is only 10 cM and a marker is located at these positions, respectively. It follows from this and from the empirical results of the simulation (Tables 1 and 2) that the smaller the QTL effect (which is unknown a priori in real mapping experiments) and the lower the population size the more advisable it is to correct for the marker impact expressed in selection coefficients.
The correction for the marker was done by dividing the frequency of the best estimate in the linkage distribution by the value of the corresponding selection coefficient (see Equation 4) whereby the number of bootstrap samples should be >200 to compute a corrected distribution suitable for calculating confidence intervals (Table 7). Figure 1C shows an average corrected distribution again from all bootstrap samples and all replicates for N = 500, Δ = 20 cM, L = 100 cM, d = 30 cM, and a = 0.5. The peaks at the marker loci observed in Figure 1, A and B, disappeared and the highest frequency of the best estimate was in the interval of the markers flanking the QTL. Right from position 60 cM the frequency became nearly equal for all positions (including marker positions).
Distributions of bootstrap QTL estimates along the simulated chromosome. (A) Empirical distribution of the QTL estimates when the recombination rate between the QTL and every marker was set to 0.5 (null distribution). (B) Empirical distribution of the QTL estimates when the QTL position was put at d = 30 cM (linkage distribution). (C) Empirical distribution of the marker-corrected QTL estimates when the QTL position was put at d = 30 cM (corrected distribution). (D) Empirical distribution of the QTL estimates when the QTL position was put at d = 20 cM (linkage distribution). (E) Empirical distribution of the marker-corrected QTL estimates when the QTL position was put at d = 20 cM (corrected distribution). A half-sib backcross family consisting of N = 500 individuals was simulated. The QTL effect was a = 0.5 and the markers were put at 0, 20, 40, 60, 80, and 100 cM, respectively. Data are from all bootstrap samples and all replicates (250 × 1000 estimates).
An average linkage distribution for a situation QTL located at a marker position (d = 20 cM) is presented in Figure 1D (N = 500, Δ = 20 cM, L = 100 cM, a = 0.5). The high peak at the QTL position can be attributed to the relatively high values of PQTL and to the selection coefficient at this point. In contrast to Figure 1B, the accumulation of the effect of a high PQTL value and a high selection coefficient is focused only on one position (20 cM) and this might be the reason for the shorter intervals with an increased conservativeness when computed with the uncorrected methods I and II. Correcting this average distribution results in Figure 1E. Right from position 60 cM the frequency of best estimates in this distribution became again nearly equal for all positions. Most of the probability mass is put between the markers at the start of the chromosome and position 40, but at the marker coinciding with the QTL position (20 cM) the frequency is significantly reduced (~1%). Similar patterns were observed at the positions of the markers flanking the QTL in the average corrected distribution in Figure 1C (QTL midway between flanking markers). Note that for this genetic configuration the null distribution follows the distribution presented in Figure 1A rather well and is therefore not shown.
Computing the null distribution: In the simulation the null distributions were derived in a simplified manner by setting the recombination rate between the QTL and any marker to 0.5. This had the advantage that the required CPU time was minimal, but is not possible in real QTL mapping experiments. Alternatively, the permutation approach as proposed by Churchill and Doerge (1994) can be used to compute the null distribution. The authors introduced this method to estimate empirical threshold values in QTL mapping, which are tailored to the experiment. By repeatedly randomly shuffling the trait values while keeping the marker data constant, they uncoupled the phenotype-genotype association, and after applying the QTL-mapping procedure, every estimate would indicate a type I error per definition. There is no reason why this distribution of the QTL estimates should not be influenced by the markers expressed in the selection coefficient for each position on the chromosome as described in our model, and so this distribution should be a very similar null distribution as obtained in this simulation. To show that the permutation method is suitable to calculate the null distribution we performed a further simulation with 1000 permutations of the phenotypic data for every replicate. The null distribution was the distribution from the QTL position estimates out of the evaluation of the permuted data and was calculated for each replicate individually. To avoid a frequency of zero in the null distribution it was assumed that each position on the simulated chromosome would receive at least one QTL position estimate from the evaluation of the permuted data. The population size was N = 500, marker spacing was 20 cM, the QTL was put at d = 30 cM, and the QTL effects were a = 0.1, 0.5, and 1.0, respectively. To save CPU time the total number of replicates for each genetic configuration was limited to 200. The confidence intervals were calculated with the weighted methods I and II as described above and the obtained results were compared with the corresponding genetic configuration in the initial simulation, where the null distribution was derived in a simplified manner. The results are presented in Table 8 and show that there were no significant differences between the two simulated series although the null distribution was built out of 1000 estimates when computed with the permutation method, compared to 250,000 estimates in the initial simulation, and although different numbers of replicates were chosen. This confirmed the finding of the low variability of the null distribution as mentioned above.
Effect of the method for calculating the null distribution on the confidence intervals for the different corrected bootstrap methods
In real QTL mapping experiments the permutation method has become a standard to calculate the error probabilities. Therefore, the null distribution can be derived as a by-product without extra computational work. It is only necessary to retain the position of the highest test statistic for each permutation and chromosome and to ensure that this distribution covers the whole chromosome. This is necessary to avoid a division by zero (Equation 4). This can be done by choosing either a sufficiently high number of permutations or, as in the present study, by adding one estimate to those positions not present in the initial distribution of the estimates from the evaluated permuted data. In the simulation, 1000 permutations were shown to be enough to compute marker-corrected confidence intervals for a 100-cM chromosome with six informative markers (Table 8) when adding additional hits. The number of permutations carried out to compute threshold values is in practice usually 10-fold higher.
Analog HPD method: The use of the analog HPD method to compute noncentral confidence intervals led to a reduced width and a reduced conservativeness of the intervals not only when the QTL was located closer to one end but also when located in the middle of the chromosome (Table 5). Intuitively one would not expect significant differences between central and noncentral intervals when the QTL is in the middle of the chromosome and the markers are located symmetrically on the chromosome. The likely explanation is that the bootstrap distributions, whether corrected for the markers or not, are tailed equally only in exceptional cases, regardless of the position of the QTL (not shown). A prerequisite for the analog HPD method to compute shorter confidence intervals is that the bootstrap distribution is not tailed equally (see appendix and Box and Tiao 1992).
Conclusions: In summary, it was shown that in the computation of bootstrap confidence intervals for estimated QTL locations it is useful, first, to correct the bootstrap linkage distribution for the marker impact and, second, to calculate noncentral intervals. Therefore, the use of the weighted method II is recommended, which combines the properties of the marker correction by the weighted method and the properties of the computation of noncentral intervals by the analog HPD method. The numerical results of the simulation showed that this method produced the shortest confidence intervals while maintaining approximately the noninclusion rate at 1 − P for a wide range of simulated genetic configurations. Provided permutation testing is used in an experiment, this method does not require extra computational work in comparison to the original bootstrap method. The combination of the permutation approach and the bootstrap approach can be termed permutation bootstrapping.
Acknowledgments
This study has benefited enormously from the critical comments of two anonymous referees. It was supported by the German Cattle Breeders Federation (ADR) and the German Ministry of Education, Science, Research and Technology (BMBF).
APPENDIX: THE ANALOG HPD METHOD
Box and Tiao (1992) define a HPD region as a region in which the probability density of every point inside is at least as large as that of any point outside. The analog HPD method was adapted from the standard HPD to compute noncentral confidence intervals from the distribution of the best QTL estimate (from the linkage distribution in the uncorrected method II and from the corrected distribution in the weighted method II, respectively). When the distribution was not tailed equally these intervals were noncentral because unequal parts were subtracted from each tail. In contrast to the standard HPD method it does not make any assumptions regarding the form of the distribution, whether it is unimodal or even multimodal. If the distribution has the latter form the initial definition of a HPD region cited above does not hold.
The 95 and 90% confidence intervals were calculated iteratively as follows:
1. Draw a line parallel to the x-axis (the simulated chromosome) with an initial low value of ε (the distance between the line and the x-axis), label the position of the x-axis from the lower outstanding intercept point between the distribution and the line as to and the corresponding position from the upper outstanding intercept point as tu.
2. Calculate the sum of the frequencies of the distribution left from to and right from tu,
3. Decide three possible cases:
Case 1. The sum is smaller than the given error rate of 1 − P. In this case repeat steps 1 and 2 with a marginally increased value of ε, and label the new positions of to and tu subsequently as
Case 2. The sum is equal to the given error rate of 1 − P. The lower interval endpoint is the position of to and the upper interval endpoint is the position of tu. The frequencies of both interval endpoints are equal. See Figure A1A.
Case 3. The sum is greater than the given error rate of 1 − P and too much was subtracted from the distribution. If it is the first round of iteration start step 1 again with a lower initial value of ε. If it is not the first round of iteration this is not caused by a too low value for ε (a marginally lower value for ε would have resulted in case 2 and not in case 1 in the previous round of iteration); it is caused by the form of the distribution, which is in this case not unimodal. To determine the interval endpoints, two differences (Δ1 and Δ2) must be calculated:
Calculating confidence intervals from bootstrap distributions with the analog HPD method. (A) The hypothetical distribution is unimodal. The lower and upper interval endpoints are built by positions on the simulated chromosome of 
and 
, respectively. (B) The hypothetical distribution is not unimodal. The upper interval endpoint is the position of 
and the position of the lower interval endpoint is calculated iteratively using Equation A1.
Footnotes
-
Communicating editor: C. Haley
- Received July 30, 2001.
- Accepted January 9, 2002.
- Copyright © 2002 by the Genetics Society of America
