Abstract
Background
Although many experiments have measurements on multiple traits, most studies performed the analysis of mapping of quantitative trait loci (QTL) for each trait separately using single trait analysis. Single trait analysis does not take advantage of possible genetic and environmental correlations between traits. In this paper, we propose a novel statistical method for multiple trait multiple interval mapping (MTMIM) of QTL for inbred line crosses. We also develop a novel scorebased method for estimating genomewide significance level of putative QTL effects suitable for the MTMIM model. The MTMIM method is implemented in the freely available and widely used Windows QTL Cartographer software.
Results
Throughout the paper, we provide compelling empirical evidences that: (1) the scorebased threshold maintains proper type I error rate and tends to keep false discovery rate within an acceptable level; (2) the MTMIM method can deliver better parameter estimates and power than single trait multiple interval mapping method; (3) an analysis of Drosophila dataset illustrates how the MTMIM method can better extract information from datasets with measurements in multiple traits.
Conclusions
The MTMIM method represents a convenient statistical framework to test hypotheses of pleiotropic QTL versus closely linked nonpleiotropic QTL, QTL by environment interaction, and to estimate the total genotypic variancecovariance matrix between traits and to decompose it in terms of QTLspecific variancecovariance matrices, therefore, providing more details on the genetic architecture of complex traits.
Keywords:
Genetic architecture; Genotypic variancecovariance; Pleiotropy; Power; QTL by environment interaction; Score statistics; Statistical geneticsBackground
Many traits that are important to agriculture, human health and evolutionary biology are quantitative in nature, influenced by multiple genes. Efficient and robust identification and mapping onto genomic positions of those genes is a very important goal in quantitative genetics. The availability of genomewide molecular markers provides the means for us to locate and map those quantitative trait loci (QTL) in a systematic way. Since the publication of interval mapping method for QTL genomewide scan [1], many statistical methods have been proposed and developed to map multiple QTL with or without epistasis in single trait in a variety of populations [2], e.g. composite interval mapping (CIM) [3,4], least squares [5,6], multiple interval mapping (MIM) [7], and Bayesian interval mapping [8,9].
Although single trait QTL mapping methods have been applied in many studies to estimate the genetic basis and architecture of complex traits, these methods did not utilize the information of genetic and environmental correlations between traits, and are not ideal for data analysis. Multiple trait analysis however can take these into account and also can formally test a number of hypotheses concerning the nature of genetic correlations, such as pleiotropy vs. close linkage and genotype by environment interaction [10]. Moreover, multiple trait analysis can allow the estimation of genetic variancecovariance matrix between traits and its decomposition in terms of individual QTL ( [11,12] pages 109110).
Multiple trait CIM [10], least squares [13] and Bayesian [14,15] methods have been available for multiple trait QTL analysis. However, these methods have not been targeted to multiple QTL for multiple traits, i.e. the whole QTL that contribute to the genetic variances and covariances. Also these methods lack appropriate criteria for assessing genomewide significance level of QTL effects. The multiple trait CIM method uses a genomewide threshold based on either asymptotic approximation of the loglikelihood ratio test (LRT) or permutation [16]. Nevertheless, when applied to multiple QTL models, the permutation test has some limitations in testing some targeted hypotheses. In this study, we have invested efforts in developing: (1) a statistical method for multiple trait multiple interval mapping (MTMIM) of QTL from inbred line crosses, and (2) a scorebased threshold for assessing significance level of QTL that is suitable for MTMIM.
In what follows, we motivate MTMIM modeling from a practical point of view, describe the MTMIM statistical model, build the likelihood function, derive parameter estimators, extend the scorebased threshold method [17] to the MTMIM model, propose a forward selection strategy to build an MTMIM model using the scorebased threshold as a criterion to assess the significance level of QTL effects, and propose a model optimization procedure to fine tune a fitted MTMIM model. We then frame the hypothesis testing of pleiotropic versus closely linked nonpleiotropic QTL, and QTL by environment interaction via the MTMIM model. Next, we implement the MTMIM model and scorebased threshold method and evaluate them with several simulated datasets. More specifically, we evaluate type I error, model fitting, and the efficiency of the test of pleiotropic versus closely linked nonpleiotropic QTL delivered by the MTMIM model. Lastly, we demonstrate the usefulness of the MTMIM model by analyzing data from an experiment with fruit flies Drosophila and draw our final considerations.
We organize this paper in a manner such that a reader less interested in the mathematical aspect of the modeling could skip the analytical derivations while being able to understand the main points regarding multiple trait multiple interval mapping of QTL.
A motivating example
We use data from a cross between fruit flies Drosophila simulans and D. mauritiana to motivate MTMIM modeling. Detailed information about the experiment can be found in [18,19]. Briefly, males from an inbred line of D. mauritiana (Rob A JJ) were crossed to females from an inbred line of D. simulans (13w JJ) to produce F_{1} hybrids. F_{1} females were then crossed to each parental species to produce two backcross populations of males, mauritiana backcross (BM) and simulans backcross (BS). These two crosses were repeated one more time to produce two independent populations from each backcross: BS1 (sample size n=186), BS2 (n=288), BM1 (n=192) and BM2 (n=299). Males from BM1 and BS1 were scored at 45 marker loci for which the two parental lines were homozygous for different alleles. Males from BM2 and BS2 were scored at 42 marker loci out of the same 45 marker loci that BM1 and BS1 were scored. The phenotypic values of each subject are: (1) average over both sides (left and right) of the first principal component of 100 Fourier coefficients of posterior lobe (PC1); (2) area of the posterior lobe (AREA); (3) average over both sides of the first principal component of 100 Fourier coefficients of the rescaled posterior lobe, rescaled so that it has unit area (ADJPC1); and (4) length of the foreleg tibia (TIBIA). While PC1 provides a measure of both size and shape of the posterior lobe, AREA and ADJPC1, on the other hand, provide measures of size and shape somewhat separately. TIBIA provides a measure of overall body size. The genotypic and phenotypic data are freely available at [20].
All variables related to the posterior lobe (PC1, ADJPC1 and AREA) were reported to be highly correlated between themselves in both BM1 and BS1, correlation larger than 0.82 [18]. Therefore, suggesting the presence of pleiotropic and/or closely linked QTL affecting size and shape. However, all variables related to the posterior lobe were weakly correlated with TIBIA. Because posterior lobe shape and size possibly share most of their developmental process components, these two traits could be tightly related mostly due to pleiotropic effects [18]. Results of composite interval mapping analysis of AREA, PC1, and ADJPC1 were very similar to each other, except for the presence of a QTL affecting both AREA and PC1 but not ADJPC1 in the interval between marker loci Ddc and eve. Therefore, this QTL affects size but not shape of the posterior lobe [18]. In this article, we use only PC1 and ADJPC1 traits and the BM1 and BM2 samples. AREA was not analyzed because it is highly correlated (0.99) with PC1 [18], and TIBIA was not analyzed because according to Liu and coauthors [18] it has small correlation with AREA and in general TIBIA is not an important factor governing the variability of posterior lobe shape. Besides, on our single trait analysis no QTL was found for TIBIA. BS1 and BS2 samples were not used for analysis because the main goal of this article is to present a novel method for QTL mapping, rather than to investigate details of the inheritance of posterior lobe shape.
We carried out MIM analysis of PC1 and ADJPC1 in the pooled samples of BM1 and BM2 (n=192+299), hearafter referred as BM data, and we found statistical evidence for seventeen genomic regions harboring QTL (Figure 1), of which twelve genomic regions showed statistical evidence of QTL affecting both traits (perhaps pleiotropic QTL), and five regions showed statistical evidence of QTL affecting either one of the traits (regions 3, 6, 9 , 12 and 15). We want to mention that in all these five regions, expect region 6, even for the trait for which the effect is not statistically significant there is still some evidence of weak putative QTL effect, as shown in the LRT profiles from the MIM analysis of PC1 and ADJPC1. Region 6, which includes marker loci Ddc and eve, was previously reported not to harbor any putative QTL with significant effect on ADJPC1 [18]. Overall, the inferred genomic regions harboring putative QTL in our MIM analysis are in strong agreement with previous inferred QTL in [18,19].
Figure 1. LRT profile of separate MIM analyses of PC1 and ADJPC1, and MTMIM analysis of PC1 and ADJPC1 (Joint) for the BM data. LRT profile of separate MIM analyses of PC1 and ADJPC1, and MTMIM analysis of PC1 and ADJPC1 (Joint) for the BM data with 10% genomewide significance level. Tick marks in the horizontal axis represent positions of genetic markers on chromosomes X, 2 and 3 (from left to right). Bold triagles bellow the horizontal axis indicate positions of mapped QTL in separate and joint analyses. Map distances are expressed in centiMorgans according to Haldane’s mapping function.
Positions of mapped QTL in regions 4, 5, 7, 10, 11, 13, 14, 16 and 17 (Figure 1) did not coincide in the MIM models of PC1 and ADJPC1. Therefore, one could hypothesize the existence of two closely linked nonpleiotropic QTL at each of these regions. In order to test the hypotheses of pleiotropic versus closely linked nonpleiotropic QTL at each one of these regions, a joint analysis of PC1 and ADJPC1 is needed. The joint analysis also allows us to partition the genotypic variancecovariance matrix between traits PC1 and ADJPC1 in terms of QTLspecific variancecovariance matrices. Thus in this motivating example, the main reasons to use the MTMIM model are: (1) to test pleiotropic versus closely linked nonpleiotropic QTL, and (2) to estimate the contribution of each QTL to the total genotypic variancecovariance matrix between traits PC1 and ADJPC1. A third reason for the MTMIM modeling, though not applicable to this specific motivating data, is the possibility to test the hypothesis of QTL by environment interaction [10].
Results and discussion
Type I error
The results show clearly an excellent agreement between estimated type I error and nominal level in the range of 1 to 15% (Figure 2).
Figure 2. Estimated and expected genomewide type I error. Estimated and expected type I error, in percentage, of LRT when using the genomewide scorebased threshold to assess significance level of putative QTL in genomewide scan of 1000 replicates.
Model size (results not shown)
The number of QTL in the MTMIM model of scenario SI was much closer to the simulated parameter (five QTL) when compared to scenario SII, for any genomewide significance level. While a QTL in both scenarios has to exceed very similar thresholds to be declared significant in the forward selection, the number of traits affected by a QTL is rather different between the two scenarios. In scenario SI all QTL have effect on all traits, while in scenario SII a QTL may have effect either on one, two or three traits. Therefore, model overparametrization makes the detection of QTL with effects on one and two traits in scenario SII more difficult. Lastly, our results show that in general the number of mapped QTL is closer to the simulated (five QTL) in the MTMIM than in the MIM model.
FDR
FDR is a very import measure of quality control in statistical analysis [21]. However, FDR is not feasibly estimated in analysis of data from traditional QTL experiments, due to the low discovery rate of putative QTL in such experiments. Nevertheless, in simulation experiments we are able to estimate FDR because we can replicate the experiment many times. We estimate FDR (Table 1) when varying the genomewide significance levels (1, 5, and 10%) and LODd support interval levels (d=1, 1.5 and 2). While FDR is expected to increase with increments in genomewide significance level, our results show that for a fixed LODd level FDR changed little with increments in genomewide significance levels, in both MIM and MTMIM models. Regarding changes in LODd level, our results show that FDR and LODd are negatively correlated, as expected. Higher levels of LODd ultimately translate into wider LODd support intervals, therefore, increasing chances of capturing the true position of QTL. FDR in the MIM and MTMIM models were very similar, except in the MIM model of trait T3 of scenario SII, which was simulated with only one QTL of small effect (heritability of 5%).
Table 1. Estimates of false discovery rate
Power
Results of power for the MIM and MTMIM models of all three scenarios clearly show a remarkable increment in power as genomewide significance levels grow less stringent, for any LODd level (Table 2  results shown for LOD1.5 level only). Based on these results as well as on those that showed almost constance of FDR across genomewide significance levels, we, hereafter, show and discuss results of 10% genomewide significance level only.
Table 2. Power of QTL identification
Results of power (10% genomewide significance level and LOD1.5) to identify QTL in the MTMIM model show that QTL affecting more traits have higher chances of being identified in the forward selection. In scenario SI, which is the most favorable among all three scenarios, all QTL have effects on all traits. Therefore, all QTL were correctly identified very often, power ≥ 97% (Table 2). In scenario SII, Q1 has effect on one trait only, Q2 on two traits, and Q3 on three traits. Power increases from Q1 (78.2%) to Q3 (97.2%) in the MTMIM model. Results also show that the MTMIM model can have lower power to identify QTL that has effects on only a small subset of traits when compared to the MIM model, due to greater genomewide threshold in the MTMIM model. For instance, MTMIM model has less power (78.2%) than MIM model (84.2%) to identify Q1, which affects only T1 (same pattern is seen for Q5). However, as the subset of traits affected by a QTL increases, power of MTMIM model overpasses power of MIM model, even when some traits are not affected by that QTL. For instance, Q2 affects T1 and T2, but not T3, nevertheless, MTMIM model identifies Q2 (95.2%) more frequently than MIM model (81.8%) (same pattern carries over to Q4). The increment in power as the number of traits affected by a QTL increases was also observed in scenario SIII.
In scenarios SII and SIII, we decomposed power of QTL identification (10% genomewide significance level and LOD1.5) into three nonoverlapping subsets (Table 3). In scenario SII, there is a subset of replicates for which a QTL affects T1 only, another subset for which a QTL affects T1 and T2 simultaneously, and finally a subset of replicates for which a QTL affects all traits (T1, T2, and T3) simultaneously. In scenario SIII, there is a subset of replicates for which a QTL affects T1 only, another subset for which a QTL effects T2 only, and finally a subset of replicates for which a QTL affects T1 and T2 simultaneously. These decompositions of power allow us to decompose the total power in the MTMIM model into QTLtrait power, therefore enabling us to measure the frequency in which a nonpleiotropic QTL is mapped as a pleiotropic one. In scenario SII, where all QTL are independent, most of power to identify a QTL is concentrated on the simulated trait affected by that QTL. For instance, in the LOD1.5 level, 66.4 out of 78.2% power (0.85 ratio) to identify Q1 is due to T1 alone, which is the only trait in which Q1 has effect on. In scenario SIII, because linkage between QTL pairs Q1 and Q2, and Q4 and Q5, the contribution of simulated traits affected by these QTL to their overall power is lower than in scenario SII, though the simulated traits still account for a large amount of power. For example, 36.8 out of 70% power (0.53 ratio) to identify Q1 is due to T1 alone, which is the only trait in which Q1 has effect on, and 46 out 68% (0.68 ratio) power to identify Q5 is due to T1 alone, which is the only trait in which Q5 has effect on. Notice that in scenario SIII Q1 was mapped as a pleiotropic QTL (subset (1,1) in Table 3) more often than Q5, i.e. 30.4 out 70% (0.43 ratio) and 20.8 out of 68% (0.31 ratio), respectively. Identification of Q1 as being pleiotropic more often than Q5 is mainly because the distance between Q1 and Q2 is shorter than the distance between Q4 and Q5, 10 and 15 cM, respectively. The smaller the distance between two nonpleiotropic QTL, the harder is to separate them in the MTMIM model. Moreover, separation of nonpleiotropic QTL is also affected by the distance between genetic markers. Linkage maps with markers closely spaced are expected to help in separating nonpleiotropic QTL. On the other hand, separation of nonpleiotropic QTL in linkage maps with sparse markers, such as the linkage map used in our simulations, is a much harder task.
Table 3. Decomposition of total power into QTLtrait power
Mean position of QTL
Our simulations show that mean estimates of QTL position in the MIM and MTMIM models have no qualitative difference and are in close agreement with the simulated parameters (Table 4). There is, though, a trend of smaller variation (measured in terms of standard error of mean) in the MTMIM than in the MIM model. Also, in the MTMIM model there is a trend of smaller variation for those QTL with effects on a larger subset of traits.
Table 4. Means of QTL position, LODd support interval coverage and length
Coverage and length of LODd support interval
In Table 4, we show the results of coverage and length of LODd support interval, and as can be seen, coverage for any LODd level are not remarkably different between the MIM and MTMIM models. However, on average the estimates of LODd support interval length were always larger in the MIM model. Differences in length are only marginal for QTL with effects on only a small subset of traits, but there are considerable differences for those QTL with effects on larger subset of traits. For instance, in scenario SII Q1 affects one trait only and it has LOD1.5 support interval mean length of 29.4 cM in the MIM and 26.4 cM in the MTMIM model. On the other hand, Q2 affects two traits and it has LOD1.5 support interval mean length of 27.7 (T1) and 27.9 (T2) in the MIM models and 21.0 cM in the MTMIM model. An interesting result is that the LOD1.5 support interval produced confidence intervals with approximately 95% coverage in both MIM and MTMIM models.
Mean effect of QTL
The average of effects of QTL in scenario SI (Table 5) shows that estimates of QTL effects in the MTMIM model are overall in close agreement with the simulated parameters, mostly because of high power in this scenario. Results of scenario SII demonstrate the robustness of the MTMIM model in estimating the effects of QTL, whereby QTL without effects on certain traits have estimates near zero, while QTL with nonzero effects have estimates with low bias. However, the robustness of the MTMIM to estimate QTL effect with low bias is less evident in scenario SIII. For instance, notice that while Q2 has zero effect on T1, its effect estimate is not close to zero. In order to understand why this bias is present in Q2 of scenario SIII, we need to understand how we matched a mapped to a simulated QTL. In the forward selection we searched and mapped pleiotropic QTL, then each mapped pleiotropic QTL was tested against the alternative hypothesis of closely linked nonpleiotropic QTL at the neighboring region of the mapped pleiotropic QTL. If the pleiotropic hypothesis was not rejected, we assumed the QTL was pleiotropic. Then, in order to apply our summary statistics, each mapped pleiotropic QTL was matched to its closest (smallest distance) simulated QTL. It could happen that a mapped pleiotropic QTL in the neighboring region of simulated Q1 and Q2 be matched to Q2, even though the major effect of the mapped pleiotropic QTL comes from Q1. Notice that when the previous situation happens, we mistakenly assign the effect of Q1 (which affects only T1) to Q2 (which presumably would not affect T1), therefore, producing biased estimated effect of Q2 on T1. The same explanation of “bias” carries over to Q4 (T1), Q1 (T2) and Q5 (T2) in scenario SIII. We quoted bias to emphasize that the bias observed in scenario SIII is not due to the MTMIM estimation per se, but rather due to our lack of ability to separate closely linked nonpleiotropic QTL or due to our criterion to match mapped to simulated QTL.
Table 5. Mean effect of QTL
The effects of all QTL were overestimated in the MIM model. This phenomena is expected due to estimation conditional on detection, the socalled “Beavis effect” [22]. A qualitative comparison of results show that overall the estimation of QTL effects in the MTMIM model are less biased than in the MIM model.
Pleiotropic versus closely linked nonpleiotropic QTL
In scenario SIII, after selecting an MTMIM model in the forward selection, each mapped pleiotropic QTL was tested against the alternative of closely linked nonpleiotropic QTL. In the bivariate model, we performed a twodimensional search for positions of putative closely linked nonpleiotropic QTL in the neighborhood of the position of each pleiotropic QTL, as suggested in [10]. The model with nonpleiotropic QTL that showed highest likelihood within the twodimensional search region was selected and tested against the model with pleiotropic QTL. We compared two criteria for model selection, the AICc and LRT. The critical value for the LRT at 5% significance level was obtained from a chisquared probability distribution with one degree of freedom.
Because Q3 was simulated as being pleiotropic, rejection of pleiotropic hypothesis for Q3 provides a measure of type I error. On the other hand, Q1 and Q2, and Q4 and Q5 were simulated as pairs of closely linked nonpleiotropic QTL. Therefore, rejection of pleiotropic hypothesis at these QTL provides a measure of power. Under our simulation setting, the LRT performed better than the AICc. The LRT was able to keep the best balance between type I error and power. Estimated frequency of rejecting pleiotropy for Q3 (4%) using the LRT agrees very well with the expected 5% nominal error rate, and estimated frequency of rejecting pleiotropy for Q1 (38%) and Q2 (36%) are satisfactory high, taking into account that Q1 and Q2 are considerably close to each other in a linkage map with markers considerably distant from each other (10 cM from markertomarker). On the other hand, the AICc criterion showed higher power for Q1 (45%) and Q2 (45%), but with a cost of high type I error for Q3 (15%). Moreover, because Q4 and Q5 are 15 cM apart from each other, the frequency of rejecting pleiotropy using LRT for these two QTL (41 and 48%, respectively) is higher than for Q1 (38%) and Q2 (36%), which are 10 cM apart from each other.
Motivating example revisited
Motivated by the fact that the joint analysis of PC1 and ADJPC1 in the Drosophila dataset could provide additional information to distinguish between genetic effects of QTL on size and shape of posterior lobe, we then analyzed these two traits with the MTMIM model. Such additional information are: (1) testing pleiotropic versus closely linked nonpleiotropic QTL, and (2) estimating the contribution of each QTL in the fitted model to the genotypic variancecovariance matrix between PC1 and ADJPC1. In what follows, we show results of the MIM and MTMIM model of the pooled samples from BM1 and BM2 (n=192+299), the BM data. We also take advantage of this dataset to test the GEMNR algorithm for maximizing the likelihood function under the MTMIM model with many QTL. Using data from a genetic experiment would provide more realistic comparisons between the GEMNR and ECM algorithms than a simulated dataset would do.
The LRT profiles of genomewide scan in the BM data (Figure 1) shows that the MTMIM model produced smaller values of LRT than the MIM model for some genomic positions, therefore, seemingly violating the expectation that the MTMIM model would produce greater LRT values than the nested MIM models [10]. Nevertheless, this violation is easily explained because not all positions of putative QTL in the MIM and MTMIM models coincide. Therefore, the MIM models are not nested within the MTMIM model shown here. Seventeen regions in the genome showed statistical evidence of putative QTL in the MTMIM model with 10% genomewide significance level (Figure 1 and Table 6).
Table 6. Estimates of QTL position and main effect on PC1 and ADJPC1 of BM data
MIM models of PC1 and ADJPC1 all together showed statistical evidence of twelve genomic regions with statistical significant QTL affecting both traits, and five regions with statistically significant QTL affecting either one of the traits (regions 3, 6, 9 , 12 and 15 shown in Figure 1 and Table 6). MTMIM model mapped these five regions either exactly or very close to their respective estimated positions in the MIM models. Moreover, the estimated effects of these five regions in the MTMIM model showed small discrepancy from those estimates in the MIM models (Table 6). Nevertheless, empirical results from our simulations suggest that both estimates of positions and effects of QTL in the MTMIM model are more accurate than in the MIM models.
Positions of QTL in regions 4, 5, 7, 10, 11, 13, 14, 16 and 17 (Figure 1 and Table 6) did not coincide with those in the MIM models of PC1 and ADJPC1. Therefore, one could hypothesize the existence of two closely linked nonpleiotropic QTL at each of these regions. We tested the hypothesis of pleiotropic QTL versus closely linked nonpleiotropic QTL at each of these regions, and on the basis of the data available the null hypothesis of pleiotropic QTL could not be rejected for any region. Thus, since PC1 contains attributes of both shape and size of posterior lobe, whereas ADJPC1 contains attributes of size only, the available data provides strong evidence that the genetic mechanisms controlling shape and size of posterior lobe are highly similar.
Partition of the phenotypic variancecovariance matrix between PC1 and ADJPC1 in terms of their environmental and genotypic components, as estimated in the MTMIM model, shows that most of the phenotypic variancecovariance between these traits is due to the genotypic component (Table 6). Moreover, we partitioned the total genotypic variancecovariance matrix in terms of QTLspecific variancecovariance matrices (Table 7) as proposed in [11] and [12] (pages 109110). This decomposition of the genotypic variancecovariance matrix shows how much of the total genotypic variancecovariance is explained by each QTL in the fitted model.
Table 7. Estimated QTLspecific (multiplied by 10^{5}) genotypic variancecovariance matrix between traits PC1 and ADJPC1
The possibility of fitting many traits and many QTL in the MTMIM model imposes severe burden in the estimation of parameters both in terms of reliability of parameter estimates (accuracy) and computation time (speed). The GEMNR and ECM algorithms are two alternative approaches suitable for parameter estimation in such complex models. We evaluate these two algorithms with the BM data by fitting an MTMIM model for PC1 and ADJPC1. The results (Figure 3) show a tremendous gain of GEMNR over ECM in terms of number of iterations, 19 and 52, respectively, as well as in terms of computing time, 8.2 and 30.6 seconds in a desktop PC, respectively. The gain in computation time from GEMNR is even more evident in genomewide scan and model selection because likelihood maximization have to be computed many times. Parameter estimates delivered in the GEMNR and ECM were very similar (Table 6).
Figure 3. Comparison of performances between ECM and GEMNR algorithms. Comparison of performances between ECM and GEMNR algorithms in terms of number of iterations required to the convergence of the likelihood function. Both algorithms were applied to an MTMIM model of traits PC1 and ADJPC1 of the BM data. The algorithms were said to have converged whenever the difference between the natural logarithm of the likelihood function of two consecutive iterations was smaller than or equal to 10^{−4}. (A) shows the values of the natural logarithm of the likelihood function at each iteration [log_{e} (L_{k})] until convergence was reached. The GEMNR algorithm began with 5 iterations of ECM algorithm. Therefore, the first 5 iterations produced identical values in the likelihood function of both algorithms, and because of that we omitted the first 4 iterations. (B) shows the difference between the natural logarithm of the likelihood function of two consecutive iterations until convergence was reached. In (B), the yaxis was rescaled via logarithm of base ten to improve graphical resolution.
Conclusions
A novel statistical method for multiple trait multiple interval mapping (MTMIM) of QTL from inbred line crosses was proposed and developed. We also proposed a novel method for estimating genomewide threshold and assessing the significance level of putative QTL effects in the MTMIM model. The method of genomewide threshold estimation is based on the scorebased resampling framework [17]. The MTMIM model has the advantage of allowing us to map QTL with effects on multiple traits, while taking advantage of information from correlations between traits. The MTMIM model has been implemented in the freely available software Windows QTL Cartographer [23].
The MTMIM model provides a comprehensive framework for QTL inference on multiple traits and the scorebased threshold serves as an essential and elegant tool for computing significance level of effects of putative QTL in the genomewide scan. The MTMIM model and scorebased threshold were evaluated through simulations. Also, we analyzed data from an experiment with Drosophila for the purpose of illustrating the MTMIM model and evaluating the performances of the GEMNR and ECM algorithms.
Results from our simulations showed many interesting features of the MTMIM model and scorebased threshold. First, the scorebased threshold maintained the type I error at a desired nominal level when no QTL effects were present in the simulated datasets. Second, discovery of spurious QTL (false discovery rate) was almost constant across genomewide significance levels of 1, 5 and 10%, while power to identify simulated QTL increased substantially as the significance level grew less stringent. Therefore, a more liberal (10%) genomewide significance level could be used in the genomewide scan, corroborating the results of C. Laurie, S. Wang, L. A. CarliniGarcia and ZB. Zeng as observed in the MIM model (unpublished observations). Third, the MTMIM model could show lower power than the MIM model for QTL with effects on only a small subset of traits. However, as the number of traits affected by a QTL increases, power in the MTMIM model overpasses power in the MIM model even when not all traits under analysis are affected by that QTL. Forth, on average the estimates of QTL position in the MIM and MTMIM models were very similar, but the MTMIM model delivers estimates with smaller sampling variances. Fifth, the LOD1.5 support interval produced confidence intervals for QTL position with approximately 95% coverage in both the MIM and MTMIM models. However, the support interval was much wider in the MIM than in MTMIM model. Overall, a qualitative comparison of results from the MIM and MTMIM models shows that effect estimates in the latter are less biased than in the former. Lastly, the LRT was shown to keep adequate type I error level when testing the null hypothesis of pleiotropic QTL against the alternative of closely linked nonpleiotropic QTL in the bivariate analysis, while it delivered reasonable power when data were generated under the alternative.
Throughout this paper, we provided compelling empirical evidences that the scorebased threshold maintained proper type I error rate and tend to give a false discovery rate within acceptable level, and that the MTMIM model can deliver better parameter estimates and power than the MIM model, and yet the MTMIM model provides a framework to test hypotheses of pleiotropic QTL versus closely linked nonpleiotropic QTL, QTL by environment interaction, and to estimate the total genotypic variancecovariances matrix between traits and to decompose it in terms of QTLspecific variancecovariance matrices. An analysis of phenotypic and genotypic data from an experiment with Drosophila illustrated the new tools present in the MTMIM model. In conclusion, the MTMIM model is a valuable tool to better extract information from experiments with measurements in multiple quantitative traits, therefore, providing more details on the genetic architecture of complex traits.
Methods
In what follows, for any matrix A, its transpose is denoted by , its inverse by A^{−1}, its u^{tℎ} row by A_{[u,·]}, its v^{tℎ} column by A_{[·,v]}, and its element in row u and column v by A_{[u,v]} .
Statistical model
Our statistical model for multiple trait multiple QTL inference for a backcross (BC) population is a linear model, in which the measurement y_{ti} of trait T (T = 1,2,· · ·, T) on each subject i (i = 1,2,· · · n) is regressed on variables x_{ir} (r = 1,2,· · · m). These variables are defined according to Cockerham genetic model [24,25]. For each subject i, x_{ir} takes either value or , depending on whether QTL r has homozygous or heterozygous genotype, respectively. The coefficient β_{tr} is called the main effect of the r^{tℎ} QTL on trait T. The linear model also includes an intercept u_{T} for each trait, it may include a subset p of epistatic effects (w_{trl}) among all pairwise QTL interactions (r and l ∈ {1,2,· · · m}), and it includes a residue e_{ti} . The linear model is:
For each subject i, let be a T × 1 vector of trait measurements, and be a T × 1 random vector assumed to be independent and identically distributed according to a multivariate normal distribution with mean vector zero and positive definite symmetric variancecovariance matrix ∑_{e}, i.e., e_{i} ∼ MVN_{T} (0,∑_{e}). For each r, let be a column vector of main effects. For each pair r and l (r < l, r = 1,2,· · · ,p) of interaction, let be a column vector of epistatic effects (b = 1,2,· · · ,p). Lastly, let be a T × 1 vector of means.
We collect all effect parameters (m main and p epistatic effect vectors) into a T × s (s = m + p) matrix , and all model parameters into a column vector θ = (θ_{1}, θ_{2},· · ·, , where for 1 ≤ b ≤ m and for m <b ≤ s, and vect(∑_{e}) is an operator that stacks the rows of ∑_{e} into a column vector one on the top of the other and then transposes it. Motivated by the fact that a QTL may not have significant effect on all traits under analysis, we allow for the insignificant parameter effects in each vector θ_{b} to be constrained to zero. Therefore, the MTMIM model allows each trait to have its own set of effect parameters, as in the seemingly unrelated regression model [26].
Likelihood function
In order to search the entire genome for significant QTL effects, the genome is partitioned into H points, usually at 1centiMorgan (cM) grid. This partition is denoted by ζ. The set of positions of m putative QTL, λ = {λ_{1}λ_{2},· · · λ_{m} }, is assumed to be a subset of ζ[27]. For any subject i, let M_{i} be the genotypic information of markers flanking the m QTL, and and be the flanking markers on left and right of QTL r, respectively. In a diploid species, a subject from a BC population generated from inbred line crosses has either genotype QQ or Qq for a locus, assuming the recurrent parent has genotype QQ. Therefore, if there are m QTL affecting a trait, there are 2^{m} possible genotypes for any subject i. Genotypes of the form G_{j} = _{Q1}Q_{2} · · ·Q_{m}, where Q_{r} ∈ {QQQq}, r = 1,2,· · · m and j = 1,2,· · · 2^{m}. Then, assuming no crossover interference between marker intervals and no more than one QTL existing within a marker interval, the probability of any genotype G_{j}, conditional on the genotypes of markers flanking the m QTL is , where the probabilities on the right hand side of this equation can be estimated via a Hidden Markov model [28].
We define an s × 2^{m} matrix Z of coded genotypes according to Cockerham genetic model [24,25]. In the matrix Z each row b, Z_{[b,·]}, corresponds to a column of effect parameters in and each column j, Z_{[·,j]}, represents a coded genotype G_{j} . If b ≤ m, Z_{[b,j]} = x_{r}, otherwise Z_{[b,j]} = x_{r} ∗x_{l}, where x_{u} (u ∈ {r, l}) is either or , depending on whether the genotype of QTL Q_{u} in G_{j} is QQ or Qq, respectively.
The individual (L_{i}) and overall likelihood (L) functions of data under the MTMIM model with m QTL are mixtures of 2^{m} multivariate normal distribution functions with different means ( ), assumed same variancecovariance (∑_{e}), and mixing probabilities p_{ij} (j = 1,2,· · · ,2^{m}), i.e., and , where y is a T × n matrix of trait measurements, and is the probability distribution function of a multivariate normal random variable y_{i} with mean and variancecovariance ∑_{e} . In what follows, ℓ_{i} (θy_{i},_{Mi},λ) and ℓ(θy,M,λ) are the natural logarithm of the individual and overall likelihood functions, respectively.
Parameter estimation
Estimation of parameters in the likelihood function is cumbersome due to mixture of distributions. The expectationmaximization (EM) [29] algorithm is very popular for parameter estimation in mixture models. The EM algorithm is very simple to program, given that efficient estimators are available for the “completedata”. Moreover, the EM algorithm guarantees that the likelihood function is nondecreasing in every iteration. However, EM may show slow convergence rate if there are many missing data, and EM does not provide standard errors of parameter estimates.
Many modifications of the EM algorithm and many hybrids of EM and GaussNewton (GN) methods have been proposed [3032]. GN methods are not guaranteed to converge when the logarithm likelihood function is not concave, but if there is convergence its rate is usually quadratic, as opposite to the linear rate of EM. Therefore, speed of convergence of GN may be much faster than EM. We describe two algorithms to obtain the maximum likelihood estimators (MLE) of parameters in the MTMIM model: expectationconditional maximization (ECM) and a hybrid of EM and NewtonRaphson called generalized EMNR (GEMNR).
Expectationconditional maximization algorithm
The EM algorithm [29] solves the incomplete logarithm likelihood function iteratively in terms of the unobserved completedata logarithm likelihood function. If the completedata logarithm likelihood function is messy and the Mstep is complex, then the EM algorithm is no longer attractive. For such cases of complicate Mstep, [33] proposed a class of generalized EM algorithm, called expectationconditional maximization (ECM). The ECM enjoys the convergence properties of the EM while simplifying the estimation of parameters. In the ECM, a complex Mstep is broken down into many simpler CMsteps, each one of them maximizes the expected completedata logarithm likelihood function conditional on some function of the parameters. Besides simplifying the Mstep, the CMstep is often faster and more stable than the Mstep because the conditional maximization are over spaces of smaller dimensions [33].
Estep: The Estep requires computation of the expectation of the completedata logarithm likelihood function, conditional on the observed data y and evaluated at a current value of θ (see Appendix). The Estep at the (v + 1) iteration consists of updating the probabilities Π_{ij} as follows:
It is worth mentioning that in the Estep above, the updating equation at step v + 1 does not use the probabilities from the previous step v, i.e, it uses p_{ij} instead of . This is the case in QTL mapping literature because the a priori probabilities are indeed exellent estimates of the conditional probabilities of QTL given the flanking markers.
The CMstep consists of maximizing the expected complete logarithm likelihood function with respect to the unknown parameters (see Appendix). CMstep without constrained parameters: We split the parameters into the groups , u, and ∑_{e} . Parameters within the same group are estimated simultaneously, while parameters in distinct groups are estimated consecutively. The parameter estimators can be shown to be:
for b ∈ {1,2,· · · ,s}.
CMstep with constrained parameters: The estimator of shown previously is not appropriate if some parameters in are constrained to zero. For instance, when estimating parameters in a model with closely linked nonpleiotropic QTL. If there exist zeroconstrained effect parameters in the MTMIM model, our strategy is to update each element in one at a time. Given the current estimate , the updating equation for the unconstrained effect parameter is:
The E and CMsteps are computed iteratively until convergence of the likelihood function. Our choice of initial values for u and ∑_{e} are the sample mean and the sample variancecovariance, respectively, and all parameters in are set to zero. In the genomewide scan, an alternative efficient choice of initial values is to use converged parameters of a previous position in the search grid. For any small positive real number ε, a stoping rule for the convergence of the likelihood function can be defined as [L(θ^{(v + 1)}Y, M, λ)−L(θ^{(v)}Y, M, λ)] / L(θ^{(v)}Y, M, λ) < ε.
It is worth mentioning that for many combinations of i and j, the probabilities p_{ij} are zero or very close to zero. Therefore, one may choose to ignore unimportant small probabilities in the computations, which may lead to significant improvement on computation time.
Generalized EM algorithm based on NewtonRaphson methods
The generalized EMNewtonRaphson (GEMNR) methods combine the EM algorithm with the NR method for maximizing the completedata logarithm likelihood function [30,31]. The hybrid methods take advantage of the EM algorithm for generating an accurate starting point for the NR algorithm, which usually has faster convergency rate. By introducing a stepsize κ^{(v)} (0 <κ^{(v)} ≤ 1) and by having the incompletedata logarithm likelihood function (ℓ) replaced by the expected completedata logarithm likelihood function (Q_{c}) in the updating NR formula, a modified version of the updating equation [32] (see Appendix) is:
The advantage of using equation (2) is that an appropriate choice of κ^{(v)} guarantees that the logarithm likelihood function increases at each iteration. So long as κ^{(v)} is chosen to make (3) positive definite, the logarithm likelihood function is guaranteed to increase at every iteration (Appendix).
where C is the Cholesky decomposition of the negative of the matrix of second order derivatives of the complete logarithm likelihood function (see Appendix) and I is an identity matrix.
To guarantee that the logarithm likelihood function is nondecreasing, [31] proposed to start the EM algorithm with five iterations to quickly approach the MLE and then to switch to NR until either convergence or decrease of the logarithm likelihood function. If the logarithm likelihood function decreases, they suggested halving the stepsize κ up to five times. If the logarithm likelihood function still decreases, they suggested to return to the EM, run five iterations, and then to switch back to NR. [31] argued that their choice of running the EM algorithm for five iterations is based on previous experiences of [34] that 95% of the change in the initial value of logarithm likelihood function until its maximum value often happens in five EM iterations.
As θ^{(ξ)} lies in the line segment from θ^{(v)} to θ^{(v + 1)}, and θ lives in highdimensional space, the choice of κ^{(v)} to make (3) positive definite may not be easy. We implemented an iterative GEMNR procedure as follows:
1. Run the ECM algorithm a couple of iterations (say five iterations);
2. Let θ^{(v)} be the parameter estimate in the v^{tℎ} EM iteration;
3. Set κ^{(v)} = 1;
4. Estimate θ^{(v + 1)} using equation (2) with the first and second order derivatives of Q_{c} (θy) evaluated at θ^{(v)} ;
5. ● If ℓ(θ^{(v + 1)}Y, M, λ) > ℓ(θ^{(v)}Y, M, λ), then set θ^{(v + 1)} as the updated parameter;
● Otherwise, keep repeating step 4 with smaller and smaller κ^{(v)}, until the likelihood function increases or until κ^{(v)} gets too small, in which case start again in step 1;
In cases in which the completedata logarithm likelihood function does not allow for closed form solution of parameter estimators, [30] have found that the GEMNR can reduce significantly the computation burden when compared to the EM algorithm. In the Appendix, we derived all expressions (first and second order derivatives of the completedata logarithm likelihood function) to implement the GEMNR algorithm.
Genomewide significance level and model selection
Scorebased threshold
We extend the score statistic [17] to assess the genomewide statistical significance level of QTL effect in the MTMIM model. Based on the individual and overall likelihood functions, we derived all required expressions to compute the score statistic to test any effect parameter in the MTMIM model (see Appendix).
Under some regular conditions, the score and LRT statistics are asymptotically equivalent in large sample [35]. But, an interesting characteristic of the score statistic is that it can be approximated by a sum of independent random components. Motivated by this characteristic and based on the decomposition of the score function [17,36] derived the largesample distribution of the score statistic for genomewide QTL mapping.
In multiple trait genomewide scan, a putative pleiotropic QTL is assumed at every position λ ∈ ζ and the significance level of its effects (main or epistatic effects) is tested against the null of no effects. For instance, assume a model with m − 1 QTL with main effects and p epistatic effects between certain QTL pairs. Assume we are scanning for a putative m^{tℎ} QTL. Let l = λ denotes the testing position of the putative QTL coming into the model. Let λ = (λ_{1},λ_{2},· · · ,λ_{(m−1)}, l) be the current positions of all m QTL in the model. Let be a T × 1 vector of effects for the new QTL coming into the model, and let be a column vector of all parameters in the model, where for 1 ≤ b ≤ m and for m <b ≤ s. Let be the column vector of nuisance parameters. Then the hypothesis H_{0} : θ_{m} = 0 versus H_{1} : θ_{m} ≠ 0 is assessed at every position l in the genome by the LRT. The genomic position with the maximum LRT among all l is assessed for the presence of a QTL via the scorebased method.
The score statistic to test H_{0} vs H_{1} can be written as [17,36], where , , and is:
where is the MLE of η under H_{0} (see Appendix for a detailed derivation of first and second order derivatives of the likelihood function).
In order to maintain equal expected variances in the resampled score and score statistic [17], we multiply by random variables z_{i} from the univariate normal distribution with mean zero and unit variance, i.e. z_{i} ∼ N(0, 1). Let be equation (4) evaluated at testing position l. Similarly let and be evaluations of and at testing position l, respectively. Then the steps of the resampling scorebased algorithm are:
1. generate n independent normal variables z_{i} (i = 1,2,· · · ,n) from N(0,1);
2. for each l, compute , . Then, compute ;
3. repeat steps 1 and 2 many times, say N times (resampling), to obtain a sequence ;
4. the scorebased threshold for a given significance αlevel is the 100(1 − α) percentile of the ascending ordered values .
If in and are assumed to be fixed and z_{i} in to be random, then: (I) The conditional distribution of on the observed data is normal with mean zero and limiting covariance as that of ; (II) From I, it follows that the distributions of and are asymptotically equivalent; and, (III) From II, it is possible to approximate the distribution of S(l) by that of s^{∗}(l) under the null hypothesis [17,37].
Model selection
The search for QTL effects on phenotypic traits consists on identifying those subset of genomic regions for which statistical tests are significant. [38] elaborated the problem of finding such a subset of genomic regions as the one of model selection, for which many tools are available in the vast literature of variable selection. However, in QTL studies the identification of a reasonable model, which maximizes the correct number of QTL while controlling the rate of false discovery is predominant over the identification of models with the smallest prediction errors, which is the major criterion for model selection [38].
The scorebased threshold can be used as a criterion to build and refine models with many QTL. Starting with a model with no QTL effect we can select putative QTL and refine the model, by including to or excluding from the MTMIM model any effects, all based on their statistical significance assessed via the scorebased method. We propose an algorithm, analogue to the algorithm described in [11], to build an initial MTMIM model and to refine it upon using the scorebased threshold criterion.
Forward selection
Assuming that model (1) starts with no QTL, one QTL is added at each step of the forward selection. In the m^{tℎ} step of the forward selection, we assume a putative pleiotropic QTL at every position l ∈ ζ (one at the time), but avoiding positions within 5 cM neighboring regions of the m − 1 QTL already in the model and compute the MLE of all parameters. For each position l, we compute the LRT statistic to test the null hypothesis versus . A putative QTL at the position with maximum LRT statistic is added to the model if the LRT statistic is larger than the scorebased threshold. Next, the effect of the selected QTL on each trait is tested individually against the null of no effect using the LRT and critical value from a chisquared probability distribution function with one degree of freedom and prespecified corrected error rate α_{c}, i.e., when T traits are analyzed jointly, the corrected significance level (Bonferroni correction) to test each effect of the m^{tℎ} QTL at an error rate α is α_{c} = α / T. Finally, any nonsignificant effect of the m^{tℎ} QTL is removed from the model, ending the m^{tℎ} step of the forward selection. The forward selection continues until no maximum LRT statistic exceeds the scorebased threshold.
Model optimization
In turns, we update the positions of all QTL in the model. We pick a QTL and hold the other QTL fixed at the positions that they were mapped. The effects of the picked QTL are then removed from the model and a new search is done within the region delimited by its two neighboring QTL, avoiding 5 cM from each neighbor (the search is performed until the end of the chromosome if no neighbor QTL is found on either side of the picked QTL). The new position of the picked QTL is set to the position of the maximum LRT statistic within the searched region and all parameters in the model are updated. This procedure is repeated until the positions of all QTL are updated.
Some suitable hypotheses in the MTMIM model
Testing pleiotropic versus closely linked nonpleiotropic QTL
Although testing for pleiotropic versus closely linked nonpleiotropic QTL is a part of model selection, we preferred to separate it from the model selection because in general this test is performed at the end of the model selection procedure, when the final model is almost fitted.
As previously stated, an advantage of multiple trait analysis is the possibility of testing for a single locus affecting multiple traits versus the alternative of two or more closely linked nonpleiotropic loci. For instance, suppose we have measurements of two traits and a total of three nonepistatic QTL at positions λ_{1}, λ_{2} and λ_{3}. The multiple trait multiple QTL pleiotropic model for a subject i would look like:
The model above assumes that all QTL have the same pattern of pleiotropy, but instead, suppose we want to test whether the last locus in model (5) is indeed two closely linked nonpleiotropic loci. The model with two pleiotropic (positions λ_{1} and λ_{2}) and two closely linked nonpleiotropic QTL (positions λ_{3} and λ_{4}) for a subject i would look like:
Or, suppose we want to test whether the last two QTL in the model (6) are both pleiotropic. The model with four pleiotropic QTL for a subject i would look like:
Many hypotheses can be formulated and tested, for example, the hypotheses of model (5) versus (6) can be stated as H_{0} : λ_{3} = λ_{4} versus H_{1} : λ_{3} ≠ λ_{4}, and the hypotheses of model (6) versus (7) can be stated as H_{0} : β_{14} = β_{23} = 0 versus H_{1} : β_{14} ≠ 0 and β_{23} ≠ 0. In general, testing whether QTL r has pleiotropic main effect or not in a subset S (S ∈ T) of traits in the model, means testing H_{0} : β_{tr} = 0 ∀ T ∈ s versus H_{1} : β_{tr} ≠ 0 for some T ∈ s. And, testing whether QTL r and l has pleiotropic epistatic effect or not in a subset S (S ∈ T) of traits in the model, means testing H_{0} : w_{trl} = 0 ∀ T ∈ s versus H_{1} : w_{trl} ≠ 0 for some T ∈ s. Model (6) illustrates a situation in which parameters are constrained to zero and the parameter estimators derived previously in the CMstep with constrained parameters are suitable.
When models are nested, the critical value to assess the strength of the LRT is straightforward, in the sense that under regular conditions the LRT has asymptotic chisquared distribution with degrees of freedom equal to the difference between the number of parameters in the full and reduced models. However, the pleiotropic and closely linkage models may not be nested (for instance, models (6) and (7)), which then requires some correction for the LRT [39,40]. The parametric bootstrap method [13] is an alternative for computing the empirical distribution of the LRT statistic in QTL mapping when models are not nested. In recognizing the test of pleiotropic versus closely linked nonpleiotropic QTL as one of model selection, we evaluate the performance of Akaike’s Information Criterion corrected (AICc) [41] and LRT, using simulation.
When a QTL has epistasis, testing this QTL for pleiotropy versus close linkage is not trivial because the test not only depends on the QTL being tested but also on any other QTL in the model that might interact with it. In general, we suggest to search for QTL main effects, and upon finishing this search to test for pleiotropy versus close linkage, and finally to search for epistasis and no longer to test pleiotropy or to test solely those QTL without epistasis.
QTL by environment interaction
The possibility of testing for QTL by environment interaction arises as another advantage of the multiple trait analysis. There are two situations in which we are able to study the differential expression of QTL. First, when the same set of genotypes are evaluated phenotypically in different environments (design I), and second when the phenotypic evaluations are done in different sets of genotypes in different environments (design II) [10]. We regard the model for analysis of data in design II as multiple population model, and thus we shall omit further discussion about it while talking about the multiple trait analysis in this paper.
Let us reiterate that in design I we regard the expression of a trait in different environments as different trait states [42]. Therefore, the index T (T = 1,2,· · · T), which was previously defined to index traits, is regarded as the environment index in what follows. With this in mind, testing whether the main effect of QTL r on a trait is statistically different or not in a subset S (S ∈ T) of environments, means testing H_{0} : β_{tr} = β_{r} ∀ T ∈ s versus H_{1} : β_{tr} ≠ β_{r} for some T ∈ s. And, testing whether QTL r and l epistatic effect on a trait is statistically different or not in a subset S (S ∈ T) of environments, means testing H_{0} : w_{trl} = 0 ∀ T ∈ s versus H_{1} : w_{trl} ≠ 0 for some T ∈ s.
The LRT may be used to evaluate the hypotheses above. The cutoff point for the test can be obtained from the chisquared probability distribution function with degrees of freedom being the difference between the number of parameters in the full (H_{1}) and reduced (H_{0}) models.
Evaluation of the MTMIM model by simulation
We implemented the MTMIM model and scorebased threshold method, and evaluated them with several simulated datasets. More specifically, we evaluated type I error, model fitting, and the efficiency of pleiotropic versus closely linked nonpleiotropic QTL testing hypothesis delivered by the MTMIM model.
Genomewide type I error
We use simulation to evaluate the proportion of falsely discovered QTL (type I error) in the analysis of datasets simulated without QTL effects. The LRT statistic is used for hypothesis testing and the scorebased threshold is used as the criterion to assess significance level of QTL effects in a genomewide scan. Each replicate has six chromosomes, each with nine markers evenly spaced 10 cM apart from each other, 300 subjects, and three quantitative traits (see Scenario S0 in Table 8). In the genomewide scan a putative pleiotropic QTL with main effects on all traits, , was assumed at each 1 cM in the genome as the alternative hypothesis. The effects of putative QTL were then tested against the simulated null hypothesis of no effects, (Scenario S0 of Table 8). For each position in the genome, we resampled the score statistic 1000 times to obtain the genomewide scorebased threshold. One thousand replicates were analyzed in this type I error study.
Table 8. Simulated genetic architecture of traits
Model fit evaluations
We use simulation to evaluate the overall performance of the MTMIM model and scorebased threshold as the criterion to assess the significance level of QTL effects in the genomewide scan. We examined the performance of the MTMIM in three different scenarios (SI, SII and SIII shown in Table 8), each evaluated with r = 500 replicates. Each replicate was simulated with six chromosomes, each with nine markers evenly spaced 10 cM apart from each other, and 300 subjects. The genetic architecture of quantitative traits in each scenario is described with details in Table 8. For each replicate we build an MTMIM model using our proposed forward selection and model optimization procedure. The genome was partitioned at 1cM grid for genomewide scan. For the sake of comparison, we also build an MIM model for each trait in each replicate using our proposed forward selection and model optimization procedure. For every position in the genome, the score statistic was resampled 800 times for the purpose of genomewide scorebased threshold estimation.
The general goal of each simulated scenario is: (SI) With a basic and favorable situation, we want to evaluate basic properties of the MTMIM model; (SII) With a mixture of QTL affecting one, two and three traits, we want to evaluate how well the MTMIM model handles the estimation of QTL with effects on only a subset of traits; (SIII) With presence of closely linked nonpleiotropic QTL and a pleiotropic QTL, we want to evaluate the MTMIM model under more complex genetic architecture. In SIII, we build an MTMIM model for each replicate using the forward selection without testing for pleiotropic versus closely linked nonpleiotropic QTL. Each MTMIM model built in the forward selection was then refined with a followup test of pleiotropic versus closely linked nonpleiotropic QTL. The pleiotropic versus closely linked nonpleiotropic test was carried out for every pleiotropic QTL in the MTMIM model.
We evaluated the MTMIM model under three genomewide significance levels: 1, 5 and 10%. For each replicate, all QTL selected in the forward selection are defined as mapped QTL. We summarize the performance of the MTMIM model with measures that are function of the logarithm of odds ratio (LOD) support interval of mapped QTL. The LODd (d = 1, 1.5, and 2) support interval of a mapped QTL is a continuous genomic region that includes the position of the mapped QTL and all positions on its left and right sides with LOD values greater than or equal to the LOD value at the position of the mapped QTL after subtraction of a positive constant d[1]. Let Q_{r}, for r ∈ {1, 2, · · · m = 5}, be a simulated QTL. A simulated QTL is defined as being paired with a mapped QTL if the simulated and mapped QTL are nearby. A mapped QTL is defined as being matched to a paired QTL if the LODd support interval of the mapped QTL includes the paired QTL. A mapped QTL is defined as mismatched if it is not matched. A simulated QTL Q_{r} is defined as identified if it has a matched QTL. For each simulated Q_{r} and for each d, let be the set of replicates for which Q_{r} is identified. We define as the number of elements in . A criterion to match mapped and simulated QTL which uses both LODd support interval and closest distance between mapped and simulated QTL is more appropriate than the usual criterion that uses closest distance alone. Our measures of model fit are: (1) False discovery rate per replicate, FDR_{b}(d), which is the ratio of number of mismatched QTL in replicate b to total number of mapped QTL in replicate b; (2) FDR over all replicates, FDR(d)= ; (3) Power to identify Q_{r}, Power(Q_{r}d)= ; (4) LODd support interval coverage of Q_{r}, c(Q_{r}d), which is the ratio of to the number of replicates for which Q_{r} is paired with a mapped QTL; (5) Mean length of LODd support interval of Q_{r}, which is the average length of LODd support intervals of Q_{r} over replicates in ; (6) Mean effect of Q_{r}, which is the average effects of Q_{r} over replicates in ; (7) Mean position of Q_{r}, which is the average positions of Q_{r} over replicates in ; and (8) Model size, which is the number of mapped QTL. These summary statistics have been proposed by C. Laurie, S. Wang, L. A. CarliniGarcia and ZB. Zeng (unpublished observations).
Appendix
Parameter estimation
Expectationconditional maximization algorithm
Let be a vector with information on “missing” genotypes of m QTL for subject i. Each if i^{tℎ} subject has genotype G_{j} (j = 1,2,· · ·,2^{m} ), otherwise . Let be a matrix containing missing information from all subjects. The joint distribution of observed and missing data for subject i is:
where , and is the probability density distribution of a multivariate normal random vector y_{i} with mean vector and variancecovariance matrix ∑_{e} . The joint distribution of observed and missing data allow us to obtain the completedata logarithm likelihood function (ℓ_{c}):
The Estep requires computation of the expectation of the completedata logarithm likelihood function, conditional on the observed data y and evaluated at current estimated values of θ (denoted here as θ^{(v)} ) [32]:
where
The CMstep consists of maximizing the expected complete logarithm likelihood function with respect to the unknown parameters through derivatives (see Section Derivatives).
NewtonRaphson method
The NR updating formula for parameter estimation [32] is:
The NR method is not very stable for complex functions because it requires accurate initial values of parameters, in certain problems, in order for right convergency. Moreover, the NR method has almost equally chances to move either in the direction of saddle points, local minima or local maxima [32]. Nevertheless, NR method has a major advantage in terms of quadratic convergence rate (when it does converge) and it can provide an estimate of the variancecovariance matrix of parameters at the limiting value of θ, θ^{∗}, through the inverse of the observed Fisher’s information matrix:
Generalized EMNewtonRaphson method
By introducing a stepsize κ^{ (v)} (0 <κ^{ (v)} ≤ 1) and by having the incompletedata logarithm likelihood function (ℓ) replaced by the expected completedata logarithm likelihood function (Q_{c}) in the updating NR formula (8), a modified version of the updating equation [32] is:
The advantage of using the modified version of the updating equation is that an appropriate choice of κ^{(v)} guarantees that the logarithm likelihood function increases at each iteration. The negative of the matrix of second order derivatives is positive definite under usual conditions. Therefore, its inverse has the Cholesky decomposition (10), where c is an upper triangular matrix.
Let θ^{(ξ)} be a point in the line segment from θ^{(v)} to θ^{(v + 1)}, the Taylor’s expansion of the completedata logarithm likelihood function around θ^{(v)} is:
Plugging θ^{(v)} from (2) into (11), and upon making some algebra using (10), we obtain:
where
and I is an identity matrix. From (12), we can see that so long as κ^{(v)} is chosen to make (13) positive definite, the logarithm likelihood function is guaranteed to increase at every iteration.
Derivatives
We provide analytical formulae of the first and second order derivatives of the logarithm of individual and overall likelihood functions of data under the MTMIM model. We borrowed useful ideas from [43,44]. These papers provide many results regarding matrix derivatives as well as their applications in multivariate analysis.
Auxiliary matrices
We assume b = 1, 2, · · · , s, i = 1, 2, · · · ,n and j = 1, 2, · · · ,2^{m} .
J_{uℓ} is a T × T matrix with 1 at positions J_{[u,ℓ]} and J_{[ℓ,u]}, and zero elsewhere I is a T × T identity matrix
First order derivatives of the logarithm of the individual likelihood function
In the following equations we use a shorthand notation ℓ_{i}(θ) = ℓ_{i}(θy_{i}, M_{i}, λ), and assume b = 1, 2, · · · ,s.
Second order derivatives of the logarithm of the overall likelihood function
In the following equations we use a shorthand notation ℓ(θ) = ℓ(θy,M,λ), and assume b = 1,2,· · · ,s, k = 1,2,· · · ,s, u = 1,2,· · · ,T, and ℓ = 1,2,· · · ,T.
First and second order derivatives of the expected completedata logarithm likelihood function
Given current estimated values of , denoted as θ^{(v)}, the first and second order derivatives of the expected completedata logarithm likelihood function are shown bellow. We assume b = 1, 2,· · · ,s, k = 1,2,· · · ,s, u = 1, 2, · · · ,T and ℓ = 1, 2, · · · ,T.
Extension to other crosses
The extension of score statistic to other cross types (for instance, intercross F_{2}, recombinant inbred lines, double haploids) is straightforward, in fact, the auxiliary matrices, expressions of first and second order derivatives of the logarithm of individual and overall likelihood functions can be straightly obtained from the general expressions derived previously. For a specific cross type, the extension consists basically of building an appropriate design matrix Z and matrix of parameters , and substituting 2^{m} in the summations by the appropriate value according to that cross type (for instance, 3^{m} for intercross F_{2}).
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
LDCES derived the analytical equations, wrote computer code, carried out the simulations and data analysis, summarized and interpreted the results, and wrote the first draft manuscript. ZBZ provided some initial results of multiple trait analysis, intellectual support criticizing the derivation and implementation of methods, and helped drafting the manuscript. SW implemented the MTMIM method in the Windows QTL Cartographer software. All authors approved this manuscript.
Acknowledgements
The authors wish to thank the editor Rongling Wu and the two anonymous reviewers for their valuable comments that improved the presentation of this paper. This work was carried out while L.D.C. E Silva was a Ph.D. candidate in Statistics at the North Carolina State University, with a joint fellowship from the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES  Brazil) and Fulbright.
References

Lander ES, Botstein D: Mapping mendelian factors underlying quantitative traits using RFLP linkage maps.
Genetics 1989, 121:185199. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Da Costa E Silva, Zeng ZB: Current progress on statistical methods for mapping quantitative trait loci from inbred line crosses.
J Biopharmaceutical Stat 2010, 20(2):454481. Publisher Full Text

Zeng ZB: Theoretical basis for separation of multiple linked gene effects in mapping quantitative trait loci.
Proc Natl Acad Sci USA 1993, 90(23):1097210976. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zeng ZB: Precision mapping of quantitative trait loci.
Genetics 1994, 136(4):14571468. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Haley CS, Knott SA: A simple regression method for mapping quantitative trait loci in line crosses using flanking markers.
Heredity 1992, 69(4):315324. PubMed Abstract  Publisher Full Text

Haley CS, Knott SA, Elsen JM: Mapping quantitative trait loci in crosses between outbred lines using least squares.
Genetics 1994, 136(3):11951207. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kao CH, Zeng ZB, Teasdale RD: Multiple interval mapping for quantitative trait loci.
Genetics 1999, 152(3):12031216. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Satagopan JM, Yandell BS, Newton MA, Osborn TC: A Bayesian approach to detect quantitative trait loci using Markov chain Monte Carlo.
Genetics 1996, 144(2):805816. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Yi N, Shriner D, Banerjee S, Mehta T, Pomp D, Yandell BS: An efficient Bayesian model selection approach for interacting quantitative trait loci models with many effects.
Genetics 2007, 176(3):18651877. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Jiang C, Zeng ZB: Multiple trait analysis of genetic mapping for quantitative trait loci.
Genetics 1995, 140(3):11111127. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zeng ZB, Kao CH, Basten CJ: Estimating the genetic architecture of quantitative traits.
Genet Res 1999, 74(3):279289. PubMed Abstract  Publisher Full Text

Maia JM: Joint analysis of multiple gene expression traits to map expression quantitative trait loci.
PhD thesis,
North Carolina State University 2007

Knott SA, Haley CS: Multitrait least squares for quantitative trait loci detection.
Genetics 2000, 156(2):899911. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Liu J, Liu Y, Liu X, Deng HW: Bayesian Mapping of Quantitative Trait Loci for Multiple Complex Traits with the Use of Variance Components.
Am J of Human Genet 2007, 81:304320. Publisher Full Text

Banerjee S, Yandell BS, Yi N: Bayesian quantitative trait loci mapping for multiple traits.
Genetics 2008, 179:22752289. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Churchill GA, Doerge RW: Empirical threshold values for quantitative trait mapping.
Genetics 1994, 138(3):963971. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zou F, Fine JP, Hu J, Lin DY: An efficient resampling method for assessing genomewide statistical significance in mapping quantitative trait Loci.
Genetics 2004, 168(4):23072316. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Liu J, Mercer JM, Stam LF, Gibson GC, Zeng ZB, Laurie CC: Genetic analysis of a morphological shape difference in the male genitalia of Drosophila simulans and D. mauritiana.
Genetics 1996, 142:11291145. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zeng ZB, Liu J, Stam LF, Kao CH, Mercer JM, Laurie CC: Genetic architecture of a morphological shape difference between two Drosophila species.
Genetics 2000, 154:299310. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Archive for the genotypic and phenotypic data of the motivating example [ ftp://statgen.ncsu.edu/pub/qtlcart/data/zengetal99/ webcite]

Storey JD, Tibshirani R: Statistical significance for genomewide studies.
Proc Natl Acad Sci USA 2003, 100(16):94409445. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Beavis WD: QTL analyses: power, precision, and accuracy. In Molecular Dissection of Complex Traits. Edited by Paterson AH. New York: CRC Press; 1998:145162.

Wang S, Basten CJ, Zeng ZB: Windows QTL Cartographer 2.51. Department of Statistics, North Carolina State University, Raleigh, NC, 2011 [ http://statgen.ncsu.edu/qtlcart/WQTLCart.htm webcite]

Kao CH, Zeng ZB: Modeling epistasis of quantitative trait loci using Cockerham’s model.
Genetics 2002, 160(3):12431261. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zeng ZB, Wang T, Zou W: Modeling quantitative trait loci and interpretation of models.
Genetics 2005, 169(3):17111725. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zellner A: An efficient method of estimating seemingly unrelated regressions and tests for aggregation bias.
J Am Stat Assoc 1962, 57(298):348368. Publisher Full Text

Yi N, Yandell BS, Churchill GA, Allison DB, Eisen EJ, Pomp D: Bayesian model selection for genomewide epistatic quantitative trait loci analysis.
Genetics 2005, 170(3):13331344. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Jiang C, Zeng ZB: Mapping quantitative trait loci with dominant and missing markers in various crosses from two inbred lines.
Genetica 1997, 101:4758. PubMed Abstract  Publisher Full Text

Dempster AP, Laird N, Rubin D: Maximum likelihood from incomplete data via the EM algorithm.

Rai SN, Matthews DE: Improving the EM algorithm.
Biometrics 1993, 49(2):587591. Publisher Full Text

Aitkin M, Aitkin I: A hybrid EM/GaussNewton algorithm for maximum likelihood in mixture distributions.
Stat Comput 1996, 6:127130. Publisher Full Text

McLachlan GJ, Krishnan T: The EM Algorithm and Extensions. New York: WileyInterscience; 1996.

Meng XL, Rubin DB: Maximum likelihood estimation via the ECM algorithm: a general framework.
Biometrika 1993, 80(2):267278. Publisher Full Text

Redner RA, Walker HR: Mixiture densities, maximum likelihood and the EM algorithm.
SIAM Rev 1984, 26(2):195239. Publisher Full Text

Chang MN, Wu R, Wu SS, Casella G: Score statistics for mapping quantitative trait loci.
Stat Appl Genet Mol Biol 2009, 8:135.
[Article 16]

Cox DR, Hinkley DV: Theoretical Statistics. London: Chapman and Hall; 1974.

Lin D: An efficient Monte Carlo approach to assessing statistical significance in genomic studies.
Bioinformatics 2005, 21(6):781787. PubMed Abstract  Publisher Full Text

Broman KW, Speed TP: A model selection approach for the identification of quantitative trait loci in experimental crosses.
J R Stat Soc Series B (Statl Methodology) 2002, 64(4):641656. Publisher Full Text

Vuong QH: Likelihood ratio tests for model selection and nonnested hypotheses.
Econometrica 1989, 57(2):307333. Publisher Full Text

Kapetanios G, Weeks M: Nonnested models and the likelihood ratio statistic: a comparison of simulation and bootstrap based tests. Technical report: University of London, 2003

Sugiura N: Further analysts of the data by Akaike’s information criterion and the finite corrections.
Commun Stat  Theory and Methods 1978, 7:1326. Publisher Full Text

Falconer DS: The problem of environment and selection.
Am Nat 1952, 86:293298. Publisher Full Text

Dwyer PS, Macphail MS: Symbolic matrix derivatives.
Ann Math Stat 1948, 19(4):517534. Publisher Full Text

Dwyer PS: Some applications of matrix derivatives in multivariate analysis.
J Am Stat Assoc 1967, 62(318):607625. Publisher Full Text