Mapping Quantitative Trait Loci in the Case of a Spike in the Phenotype Distribution
 Karl W. Broman1
 1 Address for correspondence: Department of Biostatistics, Johns Hopkins University, 615 N. Wolfe St., Baltimore, MD 212052179. Email: kbroman{at}jhsph.edu
Abstract
A common departure from the usual normality assumption in QTL mapping concerns a spike in the phenotype distribution. For example, in measurements of tumor mass, some individuals may exhibit no tumors; in measurements of time to death after a bacterial infection, some individuals may recover from the infection and fail to die. If an appreciable portion of individuals share a common phenotype value (generally either the minimum or the maximum observed phenotype), the standard approach to QTL mapping can behave poorly. We describe several alternative approaches for QTL mapping in the case of such a spike in the phenotype distribution, including the use of a twopart parametric model and a nonparametric approach based on the KruskalWallis test. The performance of the proposed procedures is assessed via computer simulation. The procedures are further illustrated with data from an intercross experiment to identify QTL contributing to variation in survival of mice following infection with Listeria monocytogenes.
THE standard approach for mapping the genetic loci (quantitative trait loci, QTL) contributing to variation in a quantitative trait makes use of the assumption that the residual environmental variation follows a normal distribution (Lander and Botstein 1989). A common departure from this assumption is to observe a spike in the phenotype distribution. For example, in Figure 1, the survival time (in hours) following infection with Listeria monocytogenes is displayed, for 116 female intercross mice (from Boyartchuket al. 2001). Approximately 30% of the mice recovered from the infection and survived to the end of the study (264 hr). Other examples include the density of metastatic tumors, with some individuals exhibiting no metastasis (Hunteret al. 2001), and gallstone weight, with some individuals having no gallstones (Wittenburget al. 2002).
Let us assume, without loss of generality, that the spike in the distribution is at 0 (which we call the null phenotype) and that all other phenotype values are strictly positive. QTL mapping under a normal model can work reasonably well in this situation if the proportion of individuals with the null phenotype is not large and the remainder of the phenotype distribution is not far above 0. However, when this is not the case, maximumlikelihood estimation under a normal mixture model can occasionally produce spurious LOD peaks in regions of low genotype information (e.g., widely spaced markers).
A simple method of analysis is to consider separately the binary trait, defined by whether or not an individual has the null phenotype, and the quantitative trait, for those individuals having a strictly positive phenotype. We develop a parametric, twopart model that allows us to combine these two analyses. In this singleQTL model, an individual with QTL genotype g has probability π_{g} of having a nonzero phenotype; if its phenotype is nonzero, the value is assumed to follow a normal distribution with mean μ_{g} and standard deviation (SD) σ.
We also describe an extension of the KruskalWallis test statistic for nonparametric interval mapping in an intercross (exactly analogous to the extension of the ranksum test described in Kruglyak and Lander 1995, which was suitable for a backcross). While Kruglyak and Lander (1995) had randomized the rank of any tied phenotypes, we assign the average rank to tied phenotypes and apply a standard correction factor. A possible advantage of the nonparametric approach is that the statistical test concerns 2 d.f., while the test derived from use of the twopart model concerns 4 d.f. and so has a larger null threshold.
We illustrate the use of these procedures with data on survival time of mice, following infection with Listeria monocytogenes (Boyartchuket al. 2001). We further study their performance via computer simulations.
METHODS
Consider n F_{2} progeny from an intercross between two inbred strains. Let y_{i} denote the quantitative phenotype for individual i. We assume, without loss of generality, that the spike in the phenotype distribution is at 0. Let z_{i} = 0 if y_{i} = 0 and z_{i} = 1 if y_{i} > 0. Consider data on a set of genetic markers, with a known genetic map. Let m_{i} denote the multipoint marker data for individual i.
Conditional and binary trait analyses: A simple approach for QTL mapping in this situation is to first analyze the quantitative phenotype, y_{i}, using only the individuals for which y_{i} > 0, by standard interval mapping using a normal model (Lander and Botstein 1989), and then separately analyze the binary trait z_{i}.
The analysis of the binary trait deserves further explanation. Xu and Atchley (1996) described maximumlikelihood estimation for a binary trait in the context of composite interval mapping (Zeng 1993, 1994). Visscher et al. (1996) and McIntyre et al. (2001) described approximate methods for analysis of binary traits. We prefer the approach of Xu and Atchley (1996). We briefly describe the special case of no marker covariates.
We consider some fixed position in the genome as the location of a putative QTL and let g_{i} = 1, 2, or 3, according to whether individual i has genotype AA, AB, or BB, respectively, at the QTL. Let us assume that the binary phenotypes, z_{i}, are independent, and let π_{j} = Pr(z_{i} = 1g_{i} = j). Given the marker data, m_{i}, but not knowing the QTL genotypes g_{i}, the z_{i} follow mixtures of Bernoulli distributions (analogous to the mixtures of normals that arise in standard interval mapping).
We assume that we may calculate p_{ij} = Pr(g_{i} = jm_{i}), the QTL genotype probabilities, given the observed multipoint marker data. Under no crossover interference and no genotyping errors, the distribution depends only on the nearest flanking typed markers, but one may also use the approach of Lincoln and Lander (1992) to take account of the presence of genotyping errors.
The likelihood for the parameters π = (π_{j}), given the observed data {(m_{i},z_{i})}, is then
We next calculate a LOD score for the test of H_{0}: π_{j} ∊ π. First note that the MLE, under H_{0}, of the common probability π is the overall proportion,
As with standard interval mapping, the likelihood under H_{0} is calculated once, while the EM algorithm is performed at each position in the genome (in practice, at 1cM steps), producing a LOD curve for each chromosome.
Twopart model: The two separate analyses described above suggest the following twopart, singleQTL model. We again consider n F_{2} progeny and some fixed position in the genome as the location of a putative QTL. Let y_{i}, z_{i}, g_{i}, and m_{i} be defined as above, and again let p_{ij} = Pr(g_{i} = jm_{i}).
We assume that the (m_{i}, y_{i}, z_{i}) are mutually independent, that Pr(z_{i} = 1g_{i} = j) =π_{j}, and that y_{i}(g_{i} = j, z_{i} = 1) ∼ normal(μ_{j}, σ^{2}). In other words, the probability that an individual with QTL genotype j has the null phenotype is 1 π_{j}; if this individual’s phenotype is nonnull, it follows a normal distribution with mean μ_{j}, depending on the QTL genotype, and with SD σ, independent of genotype.
This model contains seven parameters, θ = (π_{1}, π_{2}, π_{3}, μ_{1}, μ_{2}, μ_{3}, σ). The likelihood function is
We may again obtain MLEs with a form of the EM algorithm. Assume at iteration s + 1 we have estimates
We may calculate a LOD score for the test of H_{0}: π_{j} ∊ π, μ_{j} ∊ μ. We first note that, under H_{0}, the MLEs of the three parameters, π, μ, and σ, are
Note that in the case of complete QTL genotype information (i.e., when the putative QTL is at a marker that has been fully typed), the p_{ij} are all either 1 or 0, and the two parts of the model separate fully. As a result, the MLEs under the twopart model are exactly those obtained by the two separate analyses (the analysis of the binary trait and the conditional analysis of the quantitative trait, for those individuals with nonzero phenotype). Further, the LOD score for the twopart model is simply the sum of the LOD scores from the two separate analyses.
Nonparametric analysis: Kruglyak and Lander (1995) described an extension of the Wilcoxon ranksum test for nonparametric interval mapping in a backcross. The ranksum test is a nonparametric version of the twosample ttest. In the case of an intercross, they suggested tests for the additive or dominant effects at a putative QTL. An alternative approach is to extend the KruskalWallis test statistic, a nonparametric version of a oneway analysis of variance, for the comparison of two or more samples (e.g., see Lehmann 1975). We describe such an extension below.
Rank the phenotypes, y_{i}, from 1,..., n, and let R_{i} denote the rank for individual i. In the case of ties, use the average rank within each group of ties. We again consider some fixed position in the genome as the location of a putative QTL and let p_{ij} = Pr(g_{i} = jm_{i}), the QTL genotype probabilities for individual i, given the available multipoint marker data. Whereas, in the KruskalWallis test statistic, one considers the sum of the ranks within each group, here the exact assignment of individuals to QTL genotype groups is not known; rather, individual i has prior probability p_{ij} of belonging to group j. We follow the approach of Kruglyak and Lander (1995) and consider the expected rank sum, S_{j} = Σ_{i}p_{ij}R_{i}. We then consider the statistic
Kruglyak and Lander (1995) had randomized any tied phenotypes, a reasonable approach in the case of very few ties. In our application, however, a large proportion of the individuals share a common phenotype. Thus, rather than randomizing ties, we assign the average rank to each individual within a set of tied phenotypes. A standard correction for the case of ties is to use the statistic H′= H/D, where
EXAMPLE
To illustrate our methods, we consider the data of Boyartchuk et al. (2001), on the time to death following infection with L. monocytogenes in 116 F_{2} mice from an intercross between the BALB/cByJ and C57BL/6ByJ strains. The mice were typed at 133 markers, including 2 on the X chromosome. A histogram of the survival times (in hours) appears in Figure 1. Note that ∼30% of the mice recovered from the infection and survived to 264 hr.
We applied each of the four methods described above to these data: analysis of the binary trait, survived/died (“binary”); standard interval mapping with the log time to death, with only those mice that died (“QT”); use of our twopart model (“twopart”); and the nonparametric intervalmapping method based on the KruskalWallis test statistic (“NP”).
Genomewide LOD thresholds were obtained by permutation tests (Churchill and Doerge 1994), using 11,000 permutation replicates. The estimated 95% genomewide LOD thresholds for the four methods, binary, QT, twopart, and NP, were 3.54, 3.96, 4.91, and 3.27, respectively. The estimated standard errors (SEs) for these thresholds were ∼0.02.
Because of the large differences in the LOD thresholds for the four methods, we converted the LOD curves to a common scale, the estimated experimentwise P values derived from the permutation tests. The results indicated evidence for QTL on chromosomes 1, 5, 13, and 15. In Figure 2, the statistic log_{10}P for each method is displayed for these selected chromosomes.
The locus on chromosome 1 appears to have an effect only on the average time to death among the nonsurvivors. The locus on chromosome 5 appears to have an effect only on the chance of survival. The loci on chromosomes 13 and 15 have an effect on both the chance of survival and the average time to death among nonsurvivors. Note that the locus on chromosome 15 achieved the 5% genomewide significance level only with the nonparametric intervalmapping method.
SIMULATIONS
To better understand the relative performance of these approaches for QTL mapping in the case of a spike in the phenotype distribution, we performed a small simulation study. We first estimated the 95% genomewide LOD threshold for each method, in the case of 250 intercross individuals with 25% having the null phenotype and an autosomal genome modeled after the genetic map for the mouse described in Rowe et al. (1994), consisting of 19 autosomes with total length 1300 cM. Genetic markers were equally spaced on each chromosome, with a marker spacing of 1012 cM. (The intermarker spacing was slightly different for each chromosome, so that the chromosomal lengths could match those in the genetic maps of Roweet al. 1994.) A random 10% of the marker genotype data was missing. We simulated a phenotype that was independent of the marker data. Each individual had probability 25% of having a null phenotype; otherwise their phenotype was drawn from a normal distribution with mean 10 and SD 1.
For each of 10,000 replicates, we simulated such data under the null hypothesis of no QTL, applied each of the four methods, and recorded the maximum LOD score, genomewide, for each method. The 95th percentiles of the maximum LOD score, for the four methods, binary, QT, twopart, and NP, were 3.55, 3.53, 4.64, and 3.41, respectively. Note that the binary, QT, and NP methods have similar LOD thresholds. The LOD threshold for the twopart model is much higher, due to the fact that the corresponding statistical test concerns four free parameters, rather than two.
We also considered a fifth approach, in which one takes the maximum of the LOD scores from the binary and conditional quantitative trait analyses. For this approach, we used a Bonferroni correction and declared significant linkage if the LOD scores for either the binary trait analysis or the conditional quantitative trait analysis exceeded the corresponding 97.5% genomewide LOD thresholds, which were estimated to be 3.88 and 3.86, respectively.
To investigate the power and precision of each of these methods, we simulated data under the twopart model described above, with a single QTL located between two markers near the center of chromosome 1 (of length 103 cM). The QTL was taken to have multiplicative effect ϕ_{π} on the probabilities π_{j} and additive effect Δ_{μ} on the conditional means μ_{j}. The probabilities, π_{j}, were chosen so that π_{2} = ϕ_{π}π_{1} and
We performed 4000 simulations of 250 intercross individuals, for all pairs of effects (ϕ_{π}, Δ_{μ}), except for the case ϕ_{π} = 1, Δ_{μ} = 0. The latter corresponds to the null hypothesis of no QTL; simulations for this case were used to estimate the LOD thresholds (see above). In each case, we applied the four methods to the simulated data on chromosome 1 (containing the QTL), calculated the maximum LOD score on that chromosome, and finally calculated the power of each test, as the proportion of the simulation replicates for which the maximum LOD score exceeded the corresponding 95% genomewide LOD threshold. The power of the fifth procedure, taking the maximum of the binary and conditional quantitative trait LOD scores, was estimated as the proportion of the 4000 replicates in which either the binary or the conditional quantitative trait LOD score exceeded its corresponding 97.5% genomewide LOD threshold.
The estimated power of the procedures appears in Figure 3. In Figure 3, A and D, the QTL had effect only on the probabilities, π_{j}. In these cases, the conditional analysis of the quantitative trait had no power, and the analysis of the binary trait had the greatest power. The twopart model was somewhat inferior to the binary trait analysis, but had greater power than the nonparametric method. Use of the maximum of the binary and conditional quantitative trait LOD scores (with correction for the use of two tests) had somewhat greater power than the twopart model.
In Figure 3, G and H, the QTL had effect only on the conditional means, μ_{j}. In these cases, analysis of the binary trait had no power, and the conditional analysis of the quantitative trait had the greatest power. The results for the other methods were similar to the results in Figure 3, A and D: the twopart model was superior to the nonparametric method, but inferior to either the conditional quantitative trait analysis on its own or the maximum of the binary and conditional quantitative trait analyses.
In Figure 3, B, C, E, and F, the QTL had effect on both the probabilities, π_{j}, and the conditional means, μ_{j}. In these cases, the nonparametric method was best, although the use of the twopart model was competitive; both of these approaches showed considerable gains over either of the two separate analyses and over the maximum of the two separate analyses.
Figure 4 contains the results on the precision of QTL localization for the four basic methods. For each method and for each setting of the parameter values (ϕ_{π}, Δ_{μ}), the rootmeansquare (RMS) of the error in the estimated QTL location, among simulation replicates in which there was significant evidence for the presence of a QTL (i.e., in which the maximum LOD score exceeded the corresponding 95% genomewide LOD threshold), was calculated. Results for the conditional quantitative trait analysis (QT) for Figure 4, A and D, and for the binary trait analysis for Figure 4, G and H, are not shown, since these methods have no power to detect a QTL with the corresponding parameter settings. The results in Figure 4 mirror those in Figure 3. The methods with the highest power have the greatest precision of QTL localization (i.e., the smallest RMS error), while those with the lowest power have the lowest precision.
In summary, if a QTL has an effect only on the probabilities, π_{j}, or the conditional means, μ_{j}, greatest power to detect the QTL is obtained with the separate analysis of that aspect of the data. If a QTL has an effect on both the probabilities, π_{j}, and the conditional means, μ_{j}, the nonparametric method performed best. In all cases, analysis under the twopart model (with which the data were simulated) was second place, in terms of power. Note that further simulations, with 100 rather than 250 intercross individuals and with the proportion of individuals with the null phenotype taken to be 15 or 35% rather than 25%, gave qualitatively similar results (data not shown).
DISCUSSION
We have considered the problem of QTL mapping in the case of a spike in the phenotype distribution, a common departure from the usual normality assumption in standard interval mapping. Standard interval mapping works reasonably well when the spike is not too far from the rest of the phenotype distribution and contains only a small proportion of the individuals. When the spike is well separated and contains an appreciable proportion of the data, maximumlikelihood estimation under a normal mixture model has a tendency to produce spurious LOD score peaks in regions of low genotype information (e.g., widely spaced markers).
We developed a parametric, twopart model for QTL mapping in this situation and have described an extension of the KruskalWallis test statistic for nonparametric interval mapping in the case of an intercross. These approaches serve to combine the analysis of the binary trait with the conditional analysis of the quantitative trait among individuals with positive phenotype.
The interpretation of the results of analysis with the twopart model may deserve further explanation. A QTL identified through the twopart model may influence the probability of having a nonnull phenotype or the average phenotype among individuals with positive phenotypic values or both. Inspection of the estimated QTL effects (the
In our simulation results, most interesting was the comparison among the twopart model, the nonparametric method, and the maximum of the binary and conditional quantitative trait analyses. In the case that QTL have an effect on both the parameters π_{j} (the probability that an individual with QTL genotype j will have a positive phenotype) and μ_{j} (the conditional mean phenotype, among individuals with positive phenotype and QTL genotype j), the nonparametric approach was seen to have greater power than analysis under the twopart model; this is largely due to the fact that the genomewide LOD threshold is considerably larger for the latter method. In the case that QTL have an effect on only the π_{j} or only the μ_{j}, the maximum of the separate analyses will have greatest power, and the nonparametric method will have the least power. Thus, analysis under the twopart model is always second best. On the other hand, the overall average power, across the eight parameter settings considered herein, was greatest for the twopart model. Further, the parametric, twopart model may be more useful in consideration of multipleQTL models.
Thus, while nonparametric interval mapping is a valuable general method, analysis under the twopart model may be preferred for the situation considered here. The extensions of the twopart model for use with multiple QTL (for example, by combining a logistic model for the probabilities with a linear model for the conditional means) deserve exploration.
The methods described in this article have been implemented in the QTL mapping software, R/qtl (http://www.biostat.jhsph.edu/~kbroman/qtl), an addon package for the general statistical software, R (Ihaka and Gentleman 1996).
Acknowledgments
The author thanks Victor Boyartchuk and William Dietrich for providing the Listeria data. This work was supported in part by a Faculty Innovation Fund grant from the Johns Hopkins Bloomberg School of Public Health.
Footnotes

Communicating editor: ZB. Zeng
 Received August 1, 2002.
 Accepted November 26, 2002.
 Copyright © 2003 by the Genetics Society of America