Abstract
Quantitative trait loci (QTL) hotspots (genomic locations affecting many traits) are a common feature in genetical genomics studies and are biologically interesting since they may harbor critical regulators. Therefore, statistical procedures to assess the significance of hotspots are of key importance. One approach, randomly allocating observed QTL across the genomic locations separately by trait, implicitly assumes all traits are uncorrelated. Recently, an empirical test for QTL hotspots was proposed on the basis of the number of traits that exceed a predetermined LOD value, such as the standard permutation LOD threshold. The permutation null distribution of the maximum number of traits across all genomic locations preserves the correlation structure among the phenotypes, avoiding the detection of spurious hotspots due to nongenetic correlation induced by uncontrolled environmental factors and unmeasured variables. However, by considering only the number of traits above a threshold, without accounting for the magnitude of the LOD scores, relevant information is lost. In particular, biologically interesting hotspots composed of a moderate to small number of traits with strong LOD scores may be neglected as nonsignificant. In this article we propose a quantilebased permutation approach that simultaneously accounts for the number and the LOD scores of traits within the hotspots. By considering a sliding scale of mapping thresholds, our method can assess the statistical significance of both small and large hotspots. Although the proposed approach can be applied to any type of heritable highvolume “omic” data set, we restrict our attention to expression (e)QTL analysis. We assess and compare the performances of these three methods in simulations and we illustrate how our approach can effectively assess the significance of moderate and small hotspots with strong LOD scores in a yeast expression data set.
QUANTITATIVE trait loci (QTL) hotspots, groups of traits comapping to the same genomic location, are a common feature of genetical genomics studies. Genomic locations associated with many traits are biologically interesting since they may harbor influential regulators. Therefore, statistical procedures aiming to assess the significance of such hotspots are of key importance.
Brem et al. (2002) and Schadt et al. (2003) detected hotspots by dividing the genome of an organism into equally spaced bins, counting the number of expression traits with a QTL in each bin. A hotspot was considered significant if it had more traits with QTL than expected if the expression QTL were randomly distributed across the genome. Darvasi (2003) and PerezEnciso (2004) pointed out that these hotspots may arise as an artifact of high correlation among expression traits, rather than by the action of a common master regulator. Nongenetic mechanisms, uncontrolled environmental factors, and unmeasured variables are capable of inducing strong correlations among clusters of traits. Hence, whenever a trait shows a spurious linkage, many correlated traits will likely map to the same locus, creating a spurious expression (e)QTL hotspot. Furthermore, multipletesting and relaxedmapping thresholds may inflate the hotspots (Darvasi 2003).
West et al. (2007) and Wu et al. (2008) proposed a permutation test where the positions of the eQTL detected in the original data set are permuted across the genome separately by trait, using the distribution of the maximal number of expression traits across the genome to assess hotspot significance. This Qmethod permutes QTL positions, not phenotype or genotype data, and improves upon the permutation approaches of Brem et al. (2002) and Schadt et al. (2003) since it accounts for multiple testing across the genome. However, the Qmethod implicitly assumes traits are uncorrelated and hence underestimates the clustered pattern of spurious eQTL for correlated traits.
Breitling et al. (2008) proposed a permutation test that randomized rows in the marker data relative to rows in the trait data, preserving the correlation structure among phenotypes. The null distribution for this Nmethod of hotspot sizes depends on the number N of traits with LOD score exceeding a predetermined LOD threshold at each locus. The choice of LOD threshold is important: higher LOD thresholds yield smallersized spurious hotspots by chance under the null hypothesis of no hotspots. Two natural LOD threshold choices are the Churchill–Doerge (Churchill and Doerge 1994) singletrait LOD threshold, controlling genomewide error rate (GWER) for one trait, and a conservative permutation threshold for the maximum LOD score across all traits and all genomic locations. The former allows large hotspots by chance under the null distribution (Breitling et al. 2008). The latter favors small hotspots composed of traits with high LOD scores under the null. Which threshold is more appropriate?
We propose a quantilebased permutation approach, the NLmethod, with a sliding scale of thresholds ranging from the conservative to the singletrait threshold, jointly considering hotspot size and the distribution of LOD scores among correlated traits. Hence, even a small hotspot, with a modest number of correlated traits all having high LOD scores at a location, can be significant. The NLmethod controls the genomewide error rate across a range of possible hotspot sizes. Explicitly, we examine spurious hotspot size n, ranging from 1 to N, with N the hotspot size threshold delivered by the Nmethod using the singletrait LOD threshold. While the Nmethod yields the minimum significant hotspot size for a fixed LOD threshold, we turn the problem around and determine an empirical LOD threshold given a spurious hotspot size.
We assessed and compared the performances of the NL, N, and Qmethods using (i) simulated examples, where we generated hotspots with varying LODscore distributions for data sets with correlated and uncorrelated traits, and (ii) simulation studies, where we generated null data sets, i.e., data sets where none of the phenotypes had any QTL, and assessed the error rates of the three procedures under different levels of correlation among the traits. Application of the NLmethod to a yeast data set detected additional moderate and small hotspots considered nonsignificant by the Nmethod and avoided spurious hotspots detected by the Qmethod. This ability to assess the statistical significance of hotspots with varying sizes and LODscore distributions has the potential to yield important additional biological discoveries.
Methods
The Q and Nmethods
The now standard permutation threshold method for QTL mapping (Churchill and Doerge 1994) estimates the null distribution of the genomewide maximum LOD score by shuffling the phenotypes relative to the genotype data, breaking the association between the phenotype and the genotypes. Our interest, though, is to assess the significance of QTL hotspots. This section presents two different permutation schemes that have been used in hotspot analysis. Supporting Information, Figure S1 shows a schematic of the genotype data, the phenotype data, and the output of hotspot analysis in genetical genomics experiments.
The first permutation scheme for the Qmethod (Figure S2) derives the null distribution of hotspot sizes by permuting the cells of the observed QTL matrix along its rows, independently for each column; that is, the QTL are permuted across genomic locations separately by trait. The Qmethod does not account for the correlation structure among the phenotypes and, contrary to the Churchill–Doerge and Breitling’s permutation tests, does not break the connection between phenotypes and genotypes. The Qmethod permutation null distribution is generated under the assumption that phenotypes are uncorrelated. Violation of this assumption leads to a severe underestimation of the null distribution of hotspot sizes and to detection of many spurious hotspots, as shown by Breitling et al. (2008) in the reanalysis of the Wu et al. (2008) data and illustrated in the simulation study and examples below.
The second permutation scheme (Breitling et al. 2008) (Figure S3) used for the N and NLmethods breaks the connection between genotypes and phenotypes while preserving the correlation structure separately within each of the types of data by permuting the rows of the phenotype data matrix relative to the rows of the genotype data. Mapping analysis is redone for all traits with the permuted data. This scheme, a direct extension of Churchill and Doerge (1994), preserves the correlation structure among the phenotypes, accounting for spurious hotspots due to nongenetic correlation.
The NLmethod
In linkage analysis of a single phenotype we are usually interested in controlling the GWER of falsely detecting a QTL. For a given error rate α, we determine a singletrait mapping LOD threshold λ such that
Which threshold is more adequate depends on the underlying situation. We discard small hotspots with strong LOD scores when we adopt λ, but we miss hotspots composed of many traits with linkages barely reaching the singletrait mapping threshold using λ_{c}. Therefore, we propose a sliding scale of thresholds λ_{n}_{,}_{α} for n varying from 1 to N, where N represents the hotspot size (i.e., the number of traits with significant LOD at a given genomic location) expected by chance under the Nmethod’s permutation scheme, computed using the LODscore threshold λ.
Given a fixed mapping threshold, the Nmethod determines the hotspot size expected by chance. We turn the problem around: given a fixed spurious hotspot size n, the NLmethod determines the associated mapping threshold λ_{n}_{,}_{α}. We adopt max_{k}{qLOD_{k}(n)} as a test statistic where qLOD_{k}(n) corresponds to the nth LOD value of an ordered sample of T LOD scores, ordered from highest to lowest. Note that by taking the maximum of qLOD_{k}(n) across the genome we are able to control the genomewide error rate associated with the qLOD_{k}(n) statistic. Explicitly, we control
Observe that when n = N, the LOD threshold λ (that controls the detection of a false QTL at a GWER α) matches λ_{N}_{,}_{α} (that controls the detection of a false hotspot of size N or higher), and qLOD_{k}(N) corresponds to LOD_{k}. Therefore, when n = N, the quantity in (3) reduces to (1). Similarly, when n = 1, we have that λ_{c} = λ_{1,}_{α} and qLOD_{k}(1) = max_{t}{LOD_{t}_{,}_{k}}, so that max_{k}{qLOD_{k}(1)} = max_{t}_{,}_{k}{LOD_{t}_{,}_{k}} and (3) reduces to (2).
Finally, the quantity in (3) is the probability of detecting at least one spurious hotspot of size exactly n somewhere in the genome given that none of the traits have a QTL anywhere in the genome. However, under the null hypothesis of no QTL, detecting a hotspot of size n* > n is less likely than detecting a hotspot of size n; therefore, if a threshold λ_{n}_{,}_{α} controls the fullnull GWER for a hotspot of size n, it will also control the fullnull GWER for a hotspot of size larger than n. Below, we detail the permutation algorithm for LOD quantiles.
NLmethod algorithm:
For a fixed hotspot size, n = 1, … , N, we obtain the permutation LOD threshold that controls the genomewide error rate of detecting at least one hotspot of size n or higher, at a fixed α level as follows:
Permute the data according to the Nmethod to break the associations among genotypes and phenotypes, while keeping the correlation structure among phenotypes intact. Compute the LOD scores for all phenotypes across all genomic locations.
Process the LOD profile of each trait as follows: (i) determine the LOD peak for each chromosome, (ii) compute the LOD support interval around the peak (Lander and Botstein 1989; Dupuis and Siegmund 1999; Manichaikul et al. 2006), and (iii) set to zero the LOD scores outside the LOD support interval (and below the singletrait mapping threshold).
For a fixed hotspot size n, compute qLOD_{k}(n) for genomic positions k = 1, … , K and store its maximum.
Repeat steps 1–3, B times. The histogram of the Bpermutation samples of max_{k}{qLOD_{k}(n)} is an estimate of the null distribution of the test statistic for at least one spurious hotspot of size n or higher anywhere in the genome, given that none of the traits have a QTL anywhere in the genome.
The upper (1 − α)quantile of the permutation sample generated in step 4 is our threshold (denoted by λ_{n}_{,}_{α}).
This algorithm is analogous to the traditional permutation test, replacing LOD scores by LOD quantiles. Chen and Storey (2006) perform permutation tests for distinct quantilebased statistics, in a different context, where they consider a set of relaxed significance thresholds to detect multiple QTL for a single trait.
The LODscore processing step described in step 2 of the NLmethod algorithm confines the hotspot location on the chromosome. LOD support intervals are the most commonly used interval estimates for the location of a QTL. Following Manichaikul et al. (2006) we adopt 1.5LOD support intervals for a backcross and 1.8LOD support intervals for an intercross, decreasing the spread of the hotspot as illustrated in Figure S4.
Finally, note that instead of running a genuine, but infeasible, multipletrait joint analysis to account for the correlation structure among the traits, our strategy is to perform multiple singletrait mapping analyses with an appropriate multipletrait permutation threshold. Justification for permutation tests in the context of QTL mapping of a single phenotype is given by Churchill and Doerge (1994). A sufficient condition for a permutation test to have type I error rate held at the nominal level is that the observations are exchangeable (Good 1994, p. 203; Lehmann 1986, p. 231). Violation of the exchangeability assumption can lead to an inflation of type I error rates (Churchill and Doerge 2008). Observations are exchangeable if the joint probability of any outcome is the same irrespective of the order in which the observations are considered (Lehmann 1986, p. 231). Permutation tests remain valid for a multivariate response that can be reduced to a singlevalued test statistic (Good 1994, Chap. 5). Exchangeability of subjects under the null distribution follows by the construction of an experimental cross. At a fixed genomic location our test statistic corresponds to qLOD_{k}(n). Across the whole genome, we adopt max_{k}{qLOD_{k}(n)} as our genomewide and singlevalued test statistic.
Results
Simulated examples
In this section we illustrate the application of the Q, N, and NLmethods to two simulated data sets: one with highly correlated traits and the other with uncorrelated traits. We generated data from backcrosses composed of 112 individuals with 16 chromosomes of length 400 cM containing 185 equally spaced markers each, and phenotype data on 6000 traits. The phenotype data were generated according to the models
In both examples we simulated three hotspots: (i) a small hotspot located at 200 cM on chromosome 5 showing high LOD scores (see Figure S5, A and D), (ii) a big hotspot located at 200 cM on chromosome 7 showing LOD scores ranging from small to high (see Figure S5, B and E), and (iii) a big hotspot located at 200 cM on chromosome 15 showing LOD scores ranging from small to moderate (see Figure S5, C and F).
In both simulations we set σ^{2} = 1 and γ = 2. QTL analysis was performed using Haley–Knott regression (Haley and Knott 1992) with the R/qtl software (Broman et al. 2003). We adopted Haldane’s map function and genotype error rate of 0.0001. Because we adopted a dense genetic map (our markers are ∼2.16 cM apart), we did not consider putative QTL positions between markers.
In the first example, denoted simulated example 1, we adopted a latent effect of 1.5. In the second example, denoted simulated example 2, we adopted a latent effect of 0 and simulated uncorrelated traits. Figure S6, A and B, shows the distribution of all pairwise correlations among the 6000 traits for both simulated examples. These extreme examples illustrate the effect of phenotype correlation on QTL hotspot sizes. The correlations of the real data are actually intermediate (see Figure S6C).
Figure 1 shows the results for the Q and Nmethods for simulated example 1, using α = 0.05. Figure 1A shows the hotspot architecture computed using a singletrait LOD threshold of 3.65; i.e., at each genomic location the plot shows the number of traits with LOD score >3.65. In addition to the simulated hotspots on chromosomes 5, 7, and 15, Figure 1A shows a few spurious hotspots, including a big hotspot on chromosome 8. The blue and red lines show the N and Qmethods’ thresholds, 560 and 7, respectively. In this example the Nmethod was unable to detect any hotspots, whereas the Qmethod detected false hotspots on chromosomes 3, 6, 8, 9, 12, and 16. Figure 1, B and C, shows the hotspot size null distributions and the 5% significance thresholds for the N and Qmethods, respectively.
Figures 2 and 3 show the NLmethod analysis results for simulated example 1, using α = 0.05. Figure 2, A–D, presents the hotspot architecture inferred using four different quantilebased permutation thresholds. Figure 2A presents the hotspot architecture inferred using a LOD threshold of 7.07. Only the true hotspots (on chromosomes 5, 7, and 15) were significant by this conservative threshold. Figure 2B presents the hotspot architecture computed using a LOD threshold of 4.93 that aims to control GWER ≤ 0.05 for spurious hotspots of size 50. The hotspots on chromosomes 5, 7, and 15 were detected by this threshold. Figure 2, C and D, shows the hotspot architectures using LOD thresholds of 4.21 and 3.72, respectively. Only the hotspot on chromosome 7 was detected as significant for these thresholds. Note that neither the big spurious hotspot on chromosome 8 nor any of the other spurious hotspots shown in Figure 1A were picked up by the quantilebased thresholds.
Figure 3 connects hotspot size to quantilebased threshold. This hotspot size significance profile depicts a sliding window of hotspot size thresholds ranging from n = 1, … , N, where N = 560 corresponds to the hotspot size threshold derived from the N method. For each genomic location, the hotspot size (left axis) is significant for the LOD threshold (right axis). For example, the chromosome 5 hotspot was significant up to size 49, meaning that >1 trait mapped to the hotspot locus with LOD > 7.07, >2 traits mapped to the hotspot locus with LOD > 6.46, and so on up to hotspot size 49 where >49 traits mapped to the hotspot locus with LOD > 4.93. The hotspot on chromosome 7 was significant up to size 499, and the hotspot on chromosome 15 (higher peak) was significant for hotspot sizes 2–129 and 132–143.
The NLmethod detected only the real hotspots on chromosomes 5, 7, and 15, whereas the Nmethod did not detect any hotspots and the Qmethod detected 6 spurious hotspots, in addition to the real hotspots. The sliding window of quantilebased thresholds detected the small hotspot composed of traits with high LOD scores on chromosome 5 as well the big hotspots on chromosomes 7 and 15. Equally important, the NLmethod dismissed spurious hotspots, such as chromosome 8, composed of numerous traits with LOD scores <5.57.
Figure 4 shows the results for the Q and Nmethods for simulated example 2, using α = 0.05. Figure 4A shows the hotspot architecture. The blue and red lines show the N and Qmethod’s thresholds, 19 and 8, respectively. In this example, both the N and Qmethods were able to correctly pick up the hotspots on chromosomes 5, 7, and 15.
Comparison of Figures 1A and 4A shows that the spurious hotspots tend to be much smaller when the traits are uncorrelated (compare chromosome 8 on both plots), leading to much smaller Nmethod thresholds (compare the blue lines). The Qmethod thresholds, on the other hand, are quite close. This is expected since the Qmethod threshold depends on the number of significant QTL (we observed 3162 significant linkages in simulated example 1, against 3586 significant linkages in example 2) and not on the correlation among the traits.
Figure S7 displays the hotspot size significance profile for simulated example 2. The NLmethod also detected the hotspots on chromosomes 5, 7, and 15.
Simulation study
In this simulation study we assess and compare the error rates of the Q, N, and NLmethods under three different levels of correlation among the traits. To determine whether the methods are capable of controlling the GWER at the target levels, we conduct separate simulation experiments as follows:
We generate a “null genetical genomics data set” from a backcross composed of (i) 6000 traits, none of which is affected by a QTL, but that are nevertheless affected by a common latent variable to generate a correlation structure among the traits, and (ii) genotype data on 2960 equally spaced markers across 16 chromosomes of length 400 cM (185 markers per chromosome). Any detected QTL hotspot is spurious, arising from correlation among the traits.
We perform QTL mapping analysis, and 1.5LOD support interval processing, of the 6000 traits. For each one of the following singletrait QTL mapping permutation thresholds (that control GWER at the α = 0.01, 0.02, … , 0.10 levels, respectively), we do the following:

a. We compute the observed QTL matrix and generate the Qmethod hotspot size threshold on the basis of 1000 permutations of the observed QTL matrix. We record whether or not we see at least one spurious hotspot of size greater than the Qmethod threshold anywhere in the genome.

b. For each genomic location we count the number of traits above the singletrait LOD threshold. We compute the Nmethod hotspot size threshold on the basis of 1000 permutations of the null data set. We record whether at least one spurious hotspot of size greater than the Nmethod threshold is anywhere in the genome.

c. We compute the NLmethod LOD thresholds for spurious hotspot size thresholds ranging from 1 to the Nmethod threshold. For each NLmethod LOD threshold, λ_{n}_{,}_{α}, where n = 1, … , N, we count, at each genomic location, how many traits mapped to that genomic location with a LOD > λ_{n}_{,}_{α} and record whether there is at least one spurious hotspot of size greater than n anywhere in the genome.

We repeat the first two steps 1000 times. For each one of the three methods, the proportion of times we recorded spurious hotspots, out of the 1000 simulations, gives us an estimate of the empirical GWER associated with the method.
QTL analysis was performed as described above. Figure 5 shows the simulation results for null data sets generated using latent variable effects of 0.0, 0.25, and 1.0. The Q and Nmethods, with observed GWER (red), and target error rate (black), have two αlevels, α_{1} for QTL mapping and α_{2} for the tail area of the hotspot size permutation null distribution. Figure 5 displays the results when α_{1} = α_{2} = 0.01, 0.02, … , 0.10. The NLmethod has a single αlevel; the red curves are the observed GWERs for spurious hotspot sizes n = 1, … , N, where N represents the N method’s permutation threshold.
Figure 5, A–C, shows that for uncorrelated traits the Q and Nmethods were conservative, below target levels, whereas the NLmethod shows error rates about the right target levels for most hotspot sizes. Figure 5, D and G, shows that error rates for the Qmethod are higher than target levels when the traits are correlated and increase with correlation strength among the phenotypes. These results are expected since the Qmethod’s thresholds depend on the number of QTL detected in the unpermuted data and tend to increase with the number of phenotypes. Because we generated the same number of phenotypes in the three simulation studies, the Qmethod’s thresholds were similar. Therefore, the number and the size of the spurious QTL tend to be proportional to the correlation strength of the phenotypes. The N and NLmethods, on the other hand, are designed to cope with the correlation structure among the phenotypes and show error rates close to the target levels as shown in Figure 5, E, F, H, and I.
Yeast data set example
In this section we illustrate and compare the Q, N, and NLmethods using data generated from a cross between two parent strains of yeast: a laboratory strain and a wild isolate from a California vineyard (Brem and Kruglyak 2005). The data consist of expression measurements on 5740 transcripts measured on 112 segregant strains, with dense genotype data on 2956 markers. Processing of the expression measurements of raw data was done as described in Brem and Kruglyak (2005), with an additional step of converting the processed measurements to normal quantiles by the transformation
Hotspot analysis of the yeast data, based on the Nmethod (Figure 6A), detected significant eQTL hotpots on chromosomes 2 (second peak), 3, 12 (first peak), 14, and 15 (first peak), at a GWER of 5% according to null distribution of hotspot sizes shown in Figure 6B. The blue line represents the N method’s significance threshold of N = 96. The maximum hotspot size on chromosome 8 was 95 and almost reached significance. Nonetheless, Figure 6A also shows suggestive (although substantially smaller) peaks on chromosomes 1, 4, 5, 7, 9, 12 (second peak), 13, 15 (second peak), and 16 that did not reach significance according to the Nmethod’s significance threshold.
The red line in Figure 6A represents the Qmethod’s significance threshold of 28, derived from the null distribution of hotspot sizes shown in Figure 6C. The Qmethod detected significant hotspots on chromosomes 2 (both peaks), 3, 4, 5 (both peaks), 7, 8, 12 (both peaks), 13, 14, and 15 (both peaks).
Figure 7 shows the hotspot significance profile for the NLmethod. The major hotspots on chromosomes 2, 3, 12 (first peak), 14, and 15 (first peak) were significant across all thresholds tested up, and the hotspot on chromosome 8 was significant up to size 93. Furthermore, the NLmethod showed that the small hotspots detected by the Qmethod on chromosomes 5, 12 (second peak), 13, and 15 (second peak) might indeed be real. Nonetheless, the small hotspots on chromosomes 4 and 7, detected by the Qmethod, are less interesting than the small hotspot on chromosome 1 that was actually missed by the Qmethod.
Discussion
A common feature in genetical genomics studies of expression traits is the presence of eQTL hotspots where a single polymorphism leads to widespread downstream changes in the expression of distant genes. These genomic loci associated with many distant genes are biologically interesting since they may harbor important regulators. Statistical procedures aiming to assess the significance of such hotspots are of key importance.
Breitling et al. (2008) were the first to propose a permutation test (the Nmethod) for eQTL hotspots that accounts for the correlation structure among phenotypes due to the effect of confounders. However, the authors restricted their attention to the singletrait empirical threshold only and may have overlooked interesting hotspots composed of moderate to small numbers of traits with strong LOD scores.
In this article, we adopt the Breitling et al. permutation scheme and propose a method to determine a range of quantilebased permutation thresholds (the NLmethod) that allows us to assess the significance of hotspots on the basis of the number and the linkage strength of the traits composing those hotspots. For a fixed error rate α, our approach investigates the significance of a hotspot, using a range of N distinct mapping thresholds, where N is the smallest hotspot size that is significant by the Nmethod. For each n = 1, … , N we determine the LOD threshold that controls the genomewide error rate of detecting at least one spurious hotspot of size n or higher somewhere in the genome, at an error rate ≤α.
Our simulated examples and simulation studies show that the Qmethod performs well when the traits are uncorrelated, but detects spurious hotspots at high rates when the traits are correlated. This result is not surprising since the Qmethod implicitly assumes that the traits are uncorrelated. Molecular traits such as mRNA expression levels, metabolite concentrations, and protein levels are often highly correlated and the Qmethod is not adequate in these situations. On the other hand, our simulations suggest that the N and N and NLmethods perform adequately for correlated or uncorrelated traits, showing genomewide error rates close the target levels.
The advantage of the NLmethod over the Nmethod is that it can assess the significance of hotspots with any type of LODscore distribution. For instance, (i) a hotspot composed of many traits with moderate LOD scores will be found with thresholds close to the singletrait threshold, (ii) a hotspot consisting of a few traits with strong LOD scores will be detected with thresholds close to the conservative threshold, and (iii) a large hotspot with a range of moderate to large LOD scores will be significant at all thresholds in our sliding scale. The ability to assess the significance of these different types of hotspots can lead to important additional biological findings that might be overlooked by previous approaches, while still avoiding the detection of spurious hotspots. In the analysis of the yeast data, the hotspots on chromosomes 5, 8, 12 (second peak), 13, and 15 (second peak) have a LOD distribution of type ii. The hotspots on chromosomes 2, 3, 12 (first peak), 14, and 15 (first peak) have LOD distributions of type iii. No hotspot with LOD distribution of type i is present in the yeast data set. Note that hotspots composed of moderate to small numbers of traits with moderate LOD scores will be missed by all thresholds in our sliding scale and will be discarded as nonsignificant by our analysis. Application of the Nmethod detected only the five big hotspots on chromosomes 2, 3, 12, 14, and 15. Additionally, the simulated example 1 shows an example were the NLmethod was able to pick up the three simulated hotspots missed by the Nmethod.
The NLmethod is in a certain sense analogous to the approach proposed by Chen and Storey (2006). In the same way that Chen and Storey relax the singletrait mapping threshold by controlling the probability that a trait falsely maps to k or more genomic locations, we relax the conservative threshold by controlling the probability that ≥n traits falsely map to a common genomic location.
Even though the sliding window of thresholds delivered by the NLmethod method is more informative than the singlehotspot size threshold of the N method, these approaches have the same computational complexity. They use exactly the same permutations but summarize the results differently. Both methods are computationally intensive: reliable results require 1000 or more permutations and for each permuted data set we perform mapping analysis of several thousand traits. Thus, in general, parallel computations on a cluster are required. To reduce the computational burden, we adopted Haley–Knott regression and mapped traits with common missing phenotype data patterns as blocks. An R package called qtlhot is being submitted to Comprehensive R Archive Network (CRAN).
The approach in this article relied on singleQTL mapping methods. To examine whether an apparent hotspot could be an artifact, such as a ghost QTL (Haley and Knott 1992), we used multipleQTL methods (Manichaikul et al. 2009) for some smaller hotspots (data not shown). Most traits from these hotspots continued to map to the same location detected by singletrait analysis when we allowed for other possible QTL on the same or other chromosomes. It would be possible to extend our quantilebased permutation approach to multipleQTL mapping (Jansen 1993; Jansen and Stam 1994; Manichaikul et al. 2009) by considering the LOD profile for each QTL adjusted for all other QTL [e.g., using the addqtl or the multipleQTL mapping functions in R/qtl (Broman and Sen 2009; Arends et al. 2010)]. However, this would require considerably more computation and is left for future research.
The analysis of data sets containing groups of repetitive traits (that is, distinct traits representing slightly different measurements of a same “baseline” phenotypic trait) must be conducted with care. Repetitive traits are artifacts of the experimental design rather than indications of underlying biological processes. For instance, traits derived from oligos of same gene are often highly correlated simply because they arise from the same gene and might be picked up as a hotspot. Thus, repetitive traits can introduce artifactual hotspots that are indistinguishable statistically from biologically driven hotspots, unless this is addressed by attention to the design. Other examples of repetitive traits include (i) protein traits where one protein can exist in many variants due to posttranslational modifications and the abundance of each variant is measured and used as a separate trait and (ii) classical phenotypic traits such as flowering in Arabidopsis, where a major QTL has been investigated in a number of independent studies, under different environmental conditions, leading to a group of repetitive traits strongly mapping to the same QTL (see supplement in Fu et al. 2009). If repetitive traits are known ahead of time, they should be removed or otherwise accounted for in the analysis. For example, Fu et al. (2009) proposed organizing repetitive classical traits into disjunct phenotypic groups on the basis of trait annotations and performed hotspot analysis on the average trait per category.
Fu et al. (2009) point out that large eQTL hotspots may or may not persist when examining proteomic QTL, metabolic QTL, and phenotypic QTL gene mapping. Now that we can infer smaller hotspots composed of any of these QTL types, it may be possible to find more connections. A small hotspot could in fact be quite important to reveal genetic effects on wholebody phenotypes.
Acknowledgments
We thank Rachel Brem for sharing the yeast data set and Bill Taylor from the Center for High Throughput Computing of University of Wisconsin (Madison) for his assistance with cluster computation. This work was supported by National Council for Scientific and Technological Development (CNPq) Brazil (E.C.N.); National Cancer Institute Integrative Cancer Biology Program grant U54CA149237 and National Institutes of Health (NIH) grant R01MH090948 (to E.C.N.); National Institute of Diabetes and Digestive and Kidney Diseases grants DK66369, DK58037, and DK06639 (to A.D.A., M.P.K., A.F.B., B.S.Y., and E.C.N.); National Institute of General Medical Sciences grants PA02110 and GM06943001A2 (to B.S.Y.); The Netherlands Bioinformatics Centre, Distinguished Scientist Traveling Stipend (to R.C.J.); and by NIH grant R01GM074244 (to K.W.B.).
Footnotes
Communicating editor: J. Borevitz
 Received February 8, 2012.
 Accepted May 18, 2012.
 Copyright © 2012 by the Genetics Society of America
Available freely online through the authorsupported open access option.