Email updates

Keep up to date with the latest news and content from BMC Genetics and BioMed Central.

Open Access Methodology article

Multiple trait multiple interval mapping of quantitative trait loci from inbred line crosses

Luciano Da Costa E Silva1, Shengchu Wang1 and Zhao-Bang Zeng2*

Author Affiliations

1 Department of Statistics & Bioinformatics Research Center, North Carolina State University, Raleigh, NC 27695-7566, USA

2 Department of Genetics, North Carolina State University, Raleigh, NC 27695-7566, USA

For all author emails, please log on.

BMC Genetics 2012, 13:67  doi:10.1186/1471-2156-13-67

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2156/13/67


Received:2 March 2012
Accepted:28 June 2012
Published:1 August 2012

© 2012 Da Costa E Silva et al.; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 score-based method for estimating genome-wide 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 score-based 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 variance-covariance matrix between traits and to decompose it in terms of QTL-specific variance-covariance matrices, therefore, providing more details on the genetic architecture of complex traits.

Keywords:
Genetic architecture; Genotypic variance-covariance; Pleiotropy; Power; QTL by environment interaction; Score statistics; Statistical genetics

Background

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 genome-wide 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 genome-wide 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 variance-covariance matrix between traits and its decomposition in terms of individual QTL ( [11,12] pages 109-110).

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 genome-wide significance level of QTL effects. The multiple trait CIM method uses a genome-wide threshold based on either asymptotic approximation of the log-likelihood 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 score-based 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 score-based threshold method [17] to the MTMIM model, propose a forward selection strategy to build an MTMIM model using the score-based 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 non-pleiotropic QTL, and QTL by environment interaction via the MTMIM model. Next, we implement the MTMIM model and score-based 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 F1 hybrids. F1 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].

thumbnailFigure 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% genome-wide 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 variance-covariance matrix between traits PC1 and ADJPC1 in terms of QTL-specific variance-covariance 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 variance-covariance 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).

thumbnailFigure 2. Estimated and expected genome-wide type I error. Estimated and expected type I error, in percentage, of LRT when using the genome-wide score-based threshold to assess significance level of putative QTL in genome-wide 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 genome-wide 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 genome-wide significance levels (1, 5, and 10%) and LOD-d support interval levels (d=1, 1.5 and 2). While FDR is expected to increase with increments in genome-wide significance level, our results show that for a fixed LOD-d level FDR changed little with increments in genome-wide significance levels, in both MIM and MTMIM models. Regarding changes in LOD-d level, our results show that FDR and LOD-d are negatively correlated, as expected. Higher levels of LOD-d ultimately translate into wider LOD-d 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 genome-wide significance levels grow less stringent, for any LOD-d level (Table 2 - results shown for LOD-1.5 level only). Based on these results as well as on those that showed almost constance of FDR across genome-wide significance levels, we, hereafter, show and discuss results of 10% genome-wide significance level only.

Table 2. Power of QTL identification

Results of power (10% genome-wide significance level and LOD-1.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 genome-wide 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% genome-wide significance level and LOD-1.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 QTL-trait 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 LOD-1.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 QTL-trait 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, LOD-d support interval coverage and length

Coverage and length of LOD-d support interval

In Table 4, we show the results of coverage and length of LOD-d support interval, and as can be seen, coverage for any LOD-d level are not remarkably different between the MIM and MTMIM models. However, on average the estimates of LOD-d 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 LOD-1.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 LOD-1.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 LOD-1.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 so-called “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 two-dimensional 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 two-dimensional 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 chi-squared 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 marker-to-marker). 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 variance-covariance 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 GEM-NR 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 GEM-NR and ECM algorithms than a simulated dataset would do.

The LRT profiles of genome-wide 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% genome-wide 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 variance-covariance 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 variance-covariance between these traits is due to the genotypic component (Table 6). Moreover, we partitioned the total genotypic variance-covariance matrix in terms of QTL-specific variance-covariance matrices (Table 7) as proposed in [11] and [12] (pages 109-110). This decomposition of the genotypic variance-covariance matrix shows how much of the total genotypic variance-covariance is explained by each QTL in the fitted model.

Table 7. Estimated QTL-specific (multiplied by 105) genotypic variance-covariance 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 GEM-NR 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 GEM-NR 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 GEM-NR is even more evident in genome-wide scan and model selection because likelihood maximization have to be computed many times. Parameter estimates delivered in the GEM-NR and ECM were very similar (Table 6).

thumbnailFigure 3. Comparison of performances between ECM and GEM-NR algorithms. Comparison of performances between ECM and GEM-NR 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 [loge (Lk)] until convergence was reached. The GEM-NR 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 y-axis 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 genome-wide threshold and assessing the significance level of putative QTL effects in the MTMIM model. The method of genome-wide threshold estimation is based on the score-based 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 score-based threshold serves as an essential and elegant tool for computing significance level of effects of putative QTL in the genome-wide scan. The MTMIM model and score-based 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 GEM-NR and ECM algorithms.

Results from our simulations showed many interesting features of the MTMIM model and score-based threshold. First, the score-based 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 genome-wide 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%) genome-wide significance level could be used in the genome-wide scan, corroborating the results of C. Laurie, S. Wang, L. A. Carlini-Garcia and Z-B. 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 LOD-1.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 score-based 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 variance-covariances matrix between traits and to decompose it in terms of QTL-specific variance-covariance 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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M17">View MathML</a>, its inverse by A−1, its ut row by A[u,·], its vt 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 yti of trait T (T = 1,2,· · ·, T) on each subject i (i = 1,2,· · · n) is regressed on variables xir (r = 1,2,· · · m). These variables are defined according to Cockerham genetic model [24,25]. For each subject i, xir takes either value <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M18">View MathML</a> or <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M19">View MathML</a>, depending on whether QTL r has homozygous or heterozygous genotype, respectively. The coefficient βtr is called the main effect of the rt QTL on trait T. The linear model also includes an intercept uT for each trait, it may include a subset p of epistatic effects (wtrl) among all pairwise QTL interactions (r and l ∈ {1,2,· · · m}), and it includes a residue eti . The linear model is:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M20">View MathML</a>

(1)

For each subject i, let <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M21">View MathML</a> be a T × 1 vector of trait measurements, and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M22">View MathML</a> 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 variance-covariance matrix e, i.e., eiMVNT (0,∑e). For each r, let <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M23">View MathML</a> be a column vector of main effects. For each pair r and l (r < l, r = 1,2,· · · ,p) of interaction, let <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M24">View MathML</a> be a column vector of epistatic effects (b = 1,2,· · · ,p). Lastly, let <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M25">View MathML</a> 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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M26">View MathML</a>, and all model parameters into a column vector θ = (θ1, θ2,· · ·, <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M27">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M28">View MathML</a> for 1 ≤ bm and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M29">View MathML</a> for m <bs, 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 1-centiMorgan (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 Mi be the genotypic information of markers flanking the m QTL, and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M30">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M31">View MathML</a> 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 2m possible genotypes for any subject i. Genotypes of the form Gj = Q1Q2 · · ·Qm, where Qr ∈ {QQQq}, r = 1,2,· · · m and j = 1,2,· · · 2m. Then, assuming no crossover interference between marker intervals and no more than one QTL existing within a marker interval, the probability of any genotype Gj, conditional on the genotypes of markers flanking the m QTL is <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M32">View MathML</a>, where the probabilities on the right hand side of this equation can be estimated via a Hidden Markov model [28].

We define an s × 2m 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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M33">View MathML</a> and each column j, Z[·,j], represents a coded genotype Gj . If bm, Z[b,j] = xr, otherwise Z[b,j] = xrxl, where xu (u ∈ {r, l}) is either <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M34">View MathML</a> or <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M35">View MathML</a>, depending on whether the genotype of QTL Qu in Gj is QQ or Qq, respectively.

The individual (Li) and overall likelihood (L) functions of data under the MTMIM model with m QTL are mixtures of 2m multivariate normal distribution functions with different means ( <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M36">View MathML</a>), assumed same variance-covariance (e), and mixing probabilities pij (j = 1,2,· · · ,2m), i.e., <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M37">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M38">View MathML</a>, where y is a T × n matrix of trait measurements, and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M39">View MathML</a> is the probability distribution function of a multivariate normal random variable yi with mean <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M40">View MathML</a> and variance-covariance e . In what follows, i (θ|yi,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 expectation-maximization (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 “complete-data”. 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 Gauss-Newton (GN) methods have been proposed [30-32]. 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: expectation-conditional maximization (ECM) and a hybrid of EM and Newton-Raphson called generalized EM-NR (GEM-NR).

Expectation-conditional maximization algorithm

The EM algorithm [29] solves the incomplete logarithm likelihood function iteratively in terms of the unobserved complete-data logarithm likelihood function. If the complete-data logarithm likelihood function is messy and the M-step is complex, then the EM algorithm is no longer attractive. For such cases of complicate M-step, [33] proposed a class of generalized EM algorithm, called expectation-conditional maximization (ECM). The ECM enjoys the convergence properties of the EM while simplifying the estimation of parameters. In the ECM, a complex M-step is broken down into many simpler CM-steps, each one of them maximizes the expected complete-data logarithm likelihood function conditional on some function of the parameters. Besides simplifying the M-step, the CM-step is often faster and more stable than the M-step because the conditional maximization are over spaces of smaller dimensions [33].

E-step: The E-step requires computation of the expectation of the complete-data logarithm likelihood function, conditional on the observed data y and evaluated at a current value of θ (see Appendix). The E-step at the (v + 1) iteration consists of updating the probabilities Πij as follows:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M41">View MathML</a>

It is worth mentioning that in the E-step above, the updating equation at step v + 1 does not use the probabilities from the previous step v, i.e, it uses pij instead of <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M42">View MathML</a>. 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 CM-step consists of maximizing the expected complete logarithm likelihood function with respect to the unknown parameters (see Appendix). CM-step without constrained parameters: We split the parameters into the groups <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M43">View MathML</a>, 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:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M44">View MathML</a>

for b ∈ {1,2,· · · ,s}.

CM-step with constrained parameters: The estimator of <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M45">View MathML</a> shown previously is not appropriate if some parameters in <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M46">View MathML</a> are constrained to zero. For instance, when estimating parameters in a model with closely linked nonpleiotropic QTL. If there exist zero-constrained effect parameters in the MTMIM model, our strategy is to update each element in <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M47">View MathML</a> one at a time. Given the current estimate <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M48">View MathML</a>, the updating equation for the unconstrained effect parameter <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M49">View MathML</a> is:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M50">View MathML</a>

The E- and CM-steps 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 variance-covariance, respectively, and all parameters in <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M51">View MathML</a> are set to zero. In the genome-wide 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 pij 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 Newton-Raphson methods

The generalized EM-Newton-Raphson (GEM-NR) methods combine the EM algorithm with the NR method for maximizing the complete-data 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 step-size κ(v) (0 <κ(v) ≤ 1) and by having the incomplete-data logarithm likelihood function () replaced by the expected complete-data logarithm likelihood function (Qc) in the updating NR formula, a modified version of the updating equation [32] (see Appendix) is:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M52">View MathML</a>

(2)

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).

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M53">View MathML</a>

(3)

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 step-size κ 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 high-dimensional space, the choice of κ(v) to make (3) positive definite may not be easy. We implemented an iterative GEM-NR procedure as follows:

1. Run the ECM algorithm a couple of iterations (say five iterations);

2. Let θ(v) be the parameter estimate in the vt EM iteration;

3. Set κ(v) = 1;

4. Estimate θ(v + 1) using equation (2) with the first and second order derivatives of Qc (θ|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 complete-data logarithm likelihood function does not allow for closed form solution of parameter estimators, [30] have found that the GEM-NR 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 complete-data logarithm likelihood function) to implement the GEM-NR algorithm.

Genome-wide significance level and model selection

Score-based threshold

We extend the score statistic [17] to assess the genome-wide 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 large-sample distribution of the score statistic for genome-wide QTL mapping.

In multiple trait genome-wide 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 mt 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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M54">View MathML</a> be a T × 1 vector of effects for the new QTL coming into the model, and let <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M55">View MathML</a> be a column vector of all parameters in the model, where <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M56">View MathML</a> for 1 ≤ bm and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M57">View MathML</a> for m <bs. Let <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M58">View MathML</a> be the column vector of nuisance parameters. Then the hypothesis H0 : θm = 0 versus H1 : θm0 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 score-based method.

The score statistic to test H0 vs H1 can be written as <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M59">View MathML</a>[17,36], where <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M60">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M61">View MathML</a>, and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M62">View MathML</a> is:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M63">View MathML</a>

(4)

where <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M64">View MathML</a> is the MLE of η under H0 (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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M65">View MathML</a> by random variables zi from the univariate normal distribution with mean zero and unit variance, i.e. zi ∼ N(0, 1). Let <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M66">View MathML</a> be equation (4) evaluated at testing position l. Similarly let <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M67">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M68">View MathML</a> be evaluations of <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M69">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M70">View MathML</a> at testing position l, respectively. Then the steps of the resampling score-based algorithm are:

1. generate n independent normal variables zi (i = 1,2,· · · ,n) from N(0,1);

2. for each l, compute <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M71">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M72">View MathML</a>. Then, compute <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M73">View MathML</a>;

3. repeat steps 1 and 2 many times, say N times (resampling), to obtain a sequence <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M74">View MathML</a>;

4. the score-based threshold for a given significance α-level is the 100(1 − α) percentile of the ascending ordered values <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M75">View MathML</a>.

If <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M76">View MathML</a> in <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M77">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M78">View MathML</a> are assumed to be fixed and zi in <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M79','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M79">View MathML</a> to be random, then: (I) The conditional distribution of <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M80">View MathML</a> on the observed data is normal with mean zero and limiting covariance as that of <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M81','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M81">View MathML</a>; (II) From I, it follows that the distributions of <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M82','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M82">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M83','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M83">View MathML</a> 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 score-based 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 score-based 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 score-based 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 mt 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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M84','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M84">View MathML</a> versus <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M85','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M85">View MathML</a>. A putative QTL at the position with maximum LRT statistic is added to the model if the LRT statistic is larger than the score-based 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 chi-squared probability distribution function with one degree of freedom and pre-specified corrected error rate αc, i.e., when T traits are analyzed jointly, the corrected significance level (Bonferroni correction) to test each effect of the mt QTL at an error rate α is αc = α / T. Finally, any nonsignificant effect of the mt QTL is removed from the model, ending the mt step of the forward selection. The forward selection continues until no maximum LRT statistic exceeds the score-based 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:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M86">View MathML</a>

(5)

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:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M87">View MathML</a>

(6)

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:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M88">View MathML</a>

(7)

Many hypotheses can be formulated and tested, for example, the hypotheses of model (5) versus (6) can be stated as H0 : λ3 = λ4 versus H1 : λ3 ≠ λ4, and the hypotheses of model (6) versus (7) can be stated as H0 : β14 = β23 = 0 versus H1 : β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 H0 : βtr = 0 ∀ T ∈ s versus H1 : βtr ≠ 0 for some T ∈ s. And, testing whether QTL r and l has pleiotropic epistatic effect or not in a subset S (ST) of traits in the model, means testing H0 : wtrl = 0 ∀ T ∈ s versus H1 : wtrl ≠ 0 for some Ts. Model (6) illustrates a situation in which parameters are constrained to zero and the parameter estimators derived previously in the CM-step 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 chi-squared 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 H0 : βtr = βr ∀ T ∈ s versus H1 : β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 (ST) of environments, means testing H0 : wtrl = 0 ∀ Ts versus H1 : wtrl ≠ 0 for some Ts.

The LRT may be used to evaluate the hypotheses above. The cut-off point for the test can be obtained from the chi-squared probability distribution function with degrees of freedom being the difference between the number of parameters in the full (H1) and reduced (H0) models.

Evaluation of the MTMIM model by simulation

We implemented the MTMIM model and score-based 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.

Genome-wide 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 score-based threshold is used as the criterion to assess significance level of QTL effects in a genome-wide 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 genome-wide scan a putative pleiotropic QTL with main effects on all traits, <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M89','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M89">View MathML</a>, 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, <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M90','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M90">View MathML</a> (Scenario S0 of Table 8). For each position in the genome, we resampled the score statistic 1000 times to obtain the genome-wide score-based 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 score-based threshold as the criterion to assess the significance level of QTL effects in the genome-wide 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 1-cM grid for genome-wide 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 genome-wide score-based 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 follow-up 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 genome-wide 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 LOD-d (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 Qr, 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 LOD-d 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 Qr is defined as identified if it has a matched QTL. For each simulated Qr and for each d, let <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M91">View MathML</a> be the set of replicates for which Qr is identified. We define <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M92','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M92">View MathML</a> as the number of elements in <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M93">View MathML</a>. A criterion to match mapped and simulated QTL which uses both LOD-d 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, FDRb(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)= <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M94','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M94">View MathML</a>; (3) Power to identify Qr, Power(Qrd)= <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M95','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M95">View MathML</a>; (4) LOD-d support interval coverage of Qr, c(Qrd), which is the ratio of <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M96">View MathML</a> to the number of replicates for which Qr is paired with a mapped QTL; (5) Mean length of LOD-d support interval of Qr, which is the average length of LOD-d support intervals of Qr over replicates in <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M97','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M97">View MathML</a>; (6) Mean effect of Qr, which is the average effects of Qr over replicates in <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M98','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M98">View MathML</a>; (7) Mean position of Qr, which is the average positions of Qr over replicates in <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M99','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M99">View MathML</a>; and (8) Model size, which is the number of mapped QTL. These summary statistics have been proposed by C. Laurie, S. Wang, L. A. Carlini-Garcia and Z-B. Zeng (unpublished observations).

Appendix

Parameter estimation

Expectation-conditional maximization algorithm

Let <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M100','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M100">View MathML</a> be a vector with information on “missing” genotypes of m QTL for subject i. Each <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M101','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M101">View MathML</a> if it subject has genotype Gj (j = 1,2,· · ·,2m ), otherwise <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M102','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M102">View MathML</a>. Let <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M103','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M103">View MathML</a> be a matrix containing missing information from all subjects. The joint distribution of observed and missing data <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M104','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M104">View MathML</a> for subject i is:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M105','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M105">View MathML</a>

where <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M106','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M106">View MathML</a>, and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M107','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M107">View MathML</a> is the probability density distribution of a multivariate normal random vector yi with mean vector <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M108','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M108">View MathML</a> and variance-covariance matrix e . The joint distribution of observed and missing data allow us to obtain the complete-data logarithm likelihood function (c):

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M109','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M109">View MathML</a>

The E-step requires computation of the expectation of the complete-data logarithm likelihood function, conditional on the observed data y and evaluated at current estimated values of θ (denoted here as θ(v) ) [32]:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M110','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M110">View MathML</a>

where

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M111','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M111">View MathML</a>

The CM-step consists of maximizing the expected complete logarithm likelihood function with respect to the unknown parameters through derivatives (see Section Derivatives).

Newton-Raphson method

The NR updating formula for parameter estimation [32] is:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M112','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M112">View MathML</a>

(8)

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 variance-covariance matrix of parameters at the limiting value of θ, θ, through the inverse of the observed Fisher’s information matrix:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M113','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M113">View MathML</a>

Generalized EM-Newton-Raphson method

By introducing a step-size κ (v) (0 <κ (v) ≤ 1) and by having the incomplete-data logarithm likelihood function () replaced by the expected complete-data logarithm likelihood function (Qc) in the updating NR formula (8), a modified version of the updating equation [32] is:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M114','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M114">View MathML</a>

(9)

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.

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M115','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M115">View MathML</a>

(10)

Let θ(ξ) be a point in the line segment from θ(v) to θ(v + 1), the Taylor’s expansion of the complete-data logarithm likelihood function around θ(v) is:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M116','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M116">View MathML</a>

(11)

Plugging θ(v) from (2) into (11), and upon making some algebra using (10), we obtain:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M117','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M117">View MathML</a>

(12)

where

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M118','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M118">View MathML</a>

(13)

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, · · ·  ,2m .

Juℓ is a T × T  matrix with 1 at positions J[u,] and J[,u], and zero elsewhere I is a T × T identity matrix

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M119','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M119">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M120','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M120">View MathML</a>

First order derivatives of the logarithm of the individual likelihood function

In the following equations we use a short-hand notation i(θ) = i(θ|yi, Mi, λ), and assume b = 1, 2, · · · ,s.

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M121','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M121">View MathML</a>

Second order derivatives of the logarithm of the overall likelihood function

In the following equations we use a short-hand notation (θ) = (θ|y,M,λ), and assume b = 1,2,· · · ,s, k = 1,2,· · · ,s, u = 1,2,· · · ,T, and  = 1,2,· · · ,T.

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M122','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M122">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M123','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M123">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M124','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M124">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M125','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M125">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M126','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M126">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M127','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M127">View MathML</a>

First and second order derivatives of the expected complete-data logarithm likelihood function

Given current estimated values of <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M128','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M128">View MathML</a>, denoted as θ(v), the first and second order derivatives of the expected complete-data logarithm likelihood function are shown bellow. We assume b = 1, 2,· · · ,s, k = 1,2,· · · ,s, u = 1, 2, · · · ,T and  = 1, 2, · · · ,T.

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M129','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M129">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M130','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M130">View MathML</a>

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M131','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M131">View MathML</a>

Extension to other crosses

The extension of score statistic to other cross types (for instance, intercross F2, 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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/67/mathml/M132','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/67/mathml/M132">View MathML</a>, and substituting 2m in the summations by the appropriate value according to that cross type (for instance, 3m for intercross F2).

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

  1. Lander ES, Botstein D: Mapping mendelian factors underlying quantitative traits using RFLP linkage maps.

    Genetics 1989, 121:185-199. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  2. 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):454-481. Publisher Full Text OpenURL

  3. Zeng ZB: Theoretical basis for separation of multiple linked gene effects in mapping quantitative trait loci.

    Proc Natl Acad Sci USA 1993, 90(23):10972-10976. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  4. Zeng ZB: Precision mapping of quantitative trait loci.

    Genetics 1994, 136(4):1457-1468. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  5. Haley CS, Knott SA: A simple regression method for mapping quantitative trait loci in line crosses using flanking markers.

    Heredity 1992, 69(4):315-324. PubMed Abstract | Publisher Full Text OpenURL

  6. Haley CS, Knott SA, Elsen JM: Mapping quantitative trait loci in crosses between outbred lines using least squares.

    Genetics 1994, 136(3):1195-1207. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  7. Kao CH, Zeng ZB, Teasdale RD: Multiple interval mapping for quantitative trait loci.

    Genetics 1999, 152(3):1203-1216. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  8. 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):805-816. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  9. 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):1865-1877. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  10. Jiang C, Zeng ZB: Multiple trait analysis of genetic mapping for quantitative trait loci.

    Genetics 1995, 140(3):1111-1127. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  11. Zeng ZB, Kao CH, Basten CJ: Estimating the genetic architecture of quantitative traits.

    Genet Res 1999, 74(3):279-289. PubMed Abstract | Publisher Full Text OpenURL

  12. Maia JM: Joint analysis of multiple gene expression traits to map expression quantitative trait loci.

    PhD thesis,

    North Carolina State University 2007

    OpenURL

  13. Knott SA, Haley CS: Multitrait least squares for quantitative trait loci detection.

    Genetics 2000, 156(2):899-911. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  14. 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:304-320. Publisher Full Text OpenURL

  15. Banerjee S, Yandell BS, Yi N: Bayesian quantitative trait loci mapping for multiple traits.

    Genetics 2008, 179:2275-2289. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  16. Churchill GA, Doerge RW: Empirical threshold values for quantitative trait mapping.

    Genetics 1994, 138(3):963-971. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  17. Zou F, Fine JP, Hu J, Lin DY: An efficient resampling method for assessing genome-wide statistical significance in mapping quantitative trait Loci.

    Genetics 2004, 168(4):2307-2316. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  18. 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:1129-1145. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  19. 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:299-310. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

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

  21. Storey JD, Tibshirani R: Statistical significance for genomewide studies.

    Proc Natl Acad Sci USA 2003, 100(16):9440-9445. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

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

  23. 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]

  24. Kao CH, Zeng ZB: Modeling epistasis of quantitative trait loci using Cockerham’s model.

    Genetics 2002, 160(3):1243-1261. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  25. Zeng ZB, Wang T, Zou W: Modeling quantitative trait loci and interpretation of models.

    Genetics 2005, 169(3):1711-1725. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  26. Zellner A: An efficient method of estimating seemingly unrelated regressions and tests for aggregation bias.

    J Am Stat Assoc 1962, 57(298):348-368. Publisher Full Text OpenURL

  27. Yi N, Yandell BS, Churchill GA, Allison DB, Eisen EJ, Pomp D: Bayesian model selection for genome-wide epistatic quantitative trait loci analysis.

    Genetics 2005, 170(3):1333-1344. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  28. Jiang C, Zeng ZB: Mapping quantitative trait loci with dominant and missing markers in various crosses from two inbred lines.

    Genetica 1997, 101:47-58. PubMed Abstract | Publisher Full Text OpenURL

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

    J R Stat Soc Ser B (Methodological) 1977, 39:1-38. OpenURL

  30. Rai SN, Matthews DE: Improving the EM algorithm.

    Biometrics 1993, 49(2):587-591. Publisher Full Text OpenURL

  31. Aitkin M, Aitkin I: A hybrid EM/Gauss-Newton algorithm for maximum likelihood in mixture distributions.

    Stat Comput 1996, 6:127-130. Publisher Full Text OpenURL

  32. McLachlan GJ, Krishnan T: The EM Algorithm and Extensions. New York: Wiley-Interscience; 1996. OpenURL

  33. Meng XL, Rubin DB: Maximum likelihood estimation via the ECM algorithm: a general framework.

    Biometrika 1993, 80(2):267-278. Publisher Full Text OpenURL

  34. Redner RA, Walker HR: Mixiture densities, maximum likelihood and the EM algorithm.

    SIAM Rev 1984, 26(2):195-239. Publisher Full Text OpenURL

  35. Chang MN, Wu R, Wu SS, Casella G: Score statistics for mapping quantitative trait loci.

    Stat Appl Genet Mol Biol 2009, 8:1-35.

    [Article 16]

    OpenURL

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

  37. Lin D: An efficient Monte Carlo approach to assessing statistical significance in genomic studies.

    Bioinformatics 2005, 21(6):781-787. PubMed Abstract | Publisher Full Text OpenURL

  38. 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):641-656. Publisher Full Text OpenURL

  39. Vuong QH: Likelihood ratio tests for model selection and non-nested hypotheses.

    Econometrica 1989, 57(2):307-333. Publisher Full Text OpenURL

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

  41. Sugiura N: Further analysts of the data by Akaike’s information criterion and the finite corrections.

    Commun Stat - Theory and Methods 1978, 7:13-26. Publisher Full Text OpenURL

  42. Falconer DS: The problem of environment and selection.

    Am Nat 1952, 86:293-298. Publisher Full Text OpenURL

  43. Dwyer PS, Macphail MS: Symbolic matrix derivatives.

    Ann Math Stat 1948, 19(4):517-534. Publisher Full Text OpenURL

  44. Dwyer PS: Some applications of matrix derivatives in multivariate analysis.

    J Am Stat Assoc 1967, 62(318):607-625. Publisher Full Text OpenURL