Skip to main content

Multitrait analysis of quantitative trait loci using Bayesian composite space approach



Multitrait analysis of quantitative trait loci can capture the maximum information of experiment. The maximum-likelihood approach and the least-square approach have been developed to jointly analyze multiple traits, but it is difficult for them to include multiple QTL simultaneously into one model.


In this article, we have successfully extended Bayesian composite space approach, which is an efficient model selection method that can easily handle multiple QTL, to multitrait mapping of QTL. There are many statistical innovations of the proposed method compared with Bayesian single trait analysis. The first is that the parameters for all traits are updated jointly by vector or matrix; secondly, for QTL in the same interval that control different traits, the correlation between QTL genotypes is taken into account; thirdly, the information about the relationship of residual error between the traits is also made good use of. The superiority of the new method over separate analysis was demonstrated by both simulated and real data. The computing program was written in FORTRAN and it can be available for request.


The results suggest that the developed new method is more powerful than separate analysis.


Multitrait analysis is defined as a method that includes all traits simultaneously in a single model [1], and can take into account the correlation among all traits. Many methods have been developed for mapping QTL by combining information of multiple traits. Jiang and Zeng [2] proposed a maximum likelihood approach, and concluded that joint analysis could improve the precision of parameter estimates and had higher QTL detecting power than separate analysis. A multitrait least-square approach was proposed by Knott and Haley [3] to detect QTL. It is a method that programs easily and computes fast, and compared with separate analysis of each trait, can increase the power to detect a pleiotropic QTL and improve the precision of the location estimate. Xu et al. [1] developed a maximum likelihood approach for jointly mapping multiple binary traits, which is implemented via EM algorithm. They found that the QTL detecting power of joint analysis was higher than the sum of those of separate analysis. But after the QTL detecting power for separate analysis was redefined more reasonably by a combined power (see also [1]), the power of joint analysis was almost equal to the combined power, that is, joint analysis had almost the same power as separate analysis. For QTL parameter estimation, joint analysis can improve the precision of the QTL position estimates, but the QTL effects and their standard deviations have no obvious difference. Another class of approaches for multitrait analysis that use a dimension reduction technique was proposed by Korol et al. [4]. Mangin et al. [5] used this technique to analyze independent PCA (principal components analysis) trait, and used the PCA test values to detect QTL, which was proved to be asymptotically equivalent to the multivariate maximum-likelihood ratio test. However, the parameters of this kind of methods are often too difficult to interpret biologically. A maximum-likelihood method for multitrait mapping of QTL under outbred population was developed by Eaves et al. [6], which based on identity-by-descent (IBD) variance components model approach, and QTL effects were treated as random.

All the joint mapping approaches mentioned above were based on one-QTL model. Recently, Bayesian methodology has been used for mapping QTL [717], and the main advantage is that it can easily handle multiple QTL simultaneously. Currently, Bayesian reversible jump MCMC (RJMCMC) has become a usual method for mapping multiple QTL. Liu et al. [7] applied the method to multitrait mapping of QTL in outbred population under random effect model. However, because the dimension of RJMCMC is variable, it is always subject to poor mixing and hard to converge. Godsill [18] developed an effective Bayesian composite space method for model selection which keeps the model dimension fixed in each round of updating, and therefore it converges faster and is much easier to program. Yi et al. [1517] successfully applied the novel approach to map QTL. In this article, we extend Bayesian composite space approach to multitrait analysis under inbred line crosses, and use both simulated data and real data to demonstrate the advantages and disadvantages of the proposed method.


Simulation Study

We simulated 200 backcross individuals, and each has marker information and phenotypic records for three traits. One chromosome with length of 600 cM was investigated. Twenty-one markers were put on the genome with an average distance of 20 cM. Marker genotypes were observed for all the individuals. Thirteen QTL were added onto the genome, of which locus 96, 423, 487 and 584 had pleiotropic effects, and locus 250, 253 and 256, and locus 535 and 537 were closely linked and controlled different traits respectively. The positions and the effects of QTL for each trait are listed in Table 1. The population means for all traits were set to zero. The residual (co)variances are listed in Table 2. The heritability of each trait can be calculated as 0.728 for trait 1, 0.691 for trait 2 and 0.598 for trait 3.

Table 1 QTL Parameters and their estimates obtained from the simulated data
Table 2 The true values and their estimates of residual error (co)variance obtained from the simulated data

In order to investigate the performance of our approach, two methods were used to analyze the simulated data. The first method was the proposed multitrait analysis; the second is single-trait analysis. In single-trait analysis, we use the method 1 of [16], for the proposed method was a direct extension from it. In both multitrait analysis and single-trait analysis, the prior variance and degree of freedom of the residual error was set to zero, because no prior information was available. The prior expected number of QTL lk was 3 and the maximum number of QTL Lk equaled to the number of marker intervals (30). Therefore, the prior inclusion probability of the model indicator variable equaled to 0.1. For both methods, the MCMC ran for 1000 cycles as burn-in period (deleted) and then for additional 20,000 cycles after the burn-in. The chain was then thinned to reduce serial correlation by one observation saved every 10 cycles. The posterior sample contained 2000 (20, 000/10 = 2000) observations for the post-MCMC analysis.

The estimates of the QTL parameters for multitrait analysis and separate analysis are listed in Table 1 and Table 2. The results showed that there were no clear differences of the two methods in the estimates of the QTL positions, QTL effects and the corresponding standard deviation. Both methods can estimate QTL positions and effects, all closed to the true values.

Figure 1 and 2 respectively show the profiles of the posterior probability of the QTL positions and the 2logeBF statistic for multitrait analysis, and Figure 3 and 4 for separate analysis. From these figures, we found that both profiles of the posterior probability of QTL positions and the 2logeBF statistic for multitrait analysis are generally higher than those for separate analysis. Moreover, two additional QTL located at 483 and 245 were detected by multitrait analysis. These suggested that multitrait analysis may be more powerful than separate analysis.

Figure 1
figure 1

The profiles of the posterior probability for multitrait analysis using the simulated data. The profiles of the posterior probability obtained from multitrait analysis using the simulated data: (a) for trait 1; (b) for trait 2; (c) for trait 3. The true locations of the simulated QTL are indicated with an arrow (↑).

Figure 2
figure 2

The profiles of Bayes factors for multitrait analysis using the simulated data. The profiles of the Bayes factors (rescaled as 2logeBF and negative values are truncated as zero) obtained from multitrait analysis using the simulated data: (a) for trait 1; (b) for trait 2; (c) for trait 3. The true locations of the simulated QTL are indicated with an arrow (↑). The horizontal line indicates the critical value.

Figure 3
figure 3

The profiles of the posterior probability for single trait analysis using the simulated data. The profiles of the posterior probability obtained from separate analysis using the simulated data: (a) for trait 1; (b) for trait 2; (c) for trait 3. The true locations of the simulated QTL are indicated with an arrow (↑).

Figure 4
figure 4

The profiles of Bayes factors for single trait analysis using the simulated data. The profiles of Bayes factors (rescaled as 2 logeBF and negative values are truncated as zero) obtained from separate analysis using the simulated data: (a) for trait 1; (b) for trait 2; (c) for trait 3. The true locations of the simulated QTL are indicated with an arrow (↑). The critical value is given as horizontal line.

Real data analysis

We applied the new method to analyze the data from the North American Barley Genome Mapping Project [22]. The DH population included 150 lines (n = 150), each of which was genotyped for 223 codominant markers. These markers covered ~1500 cM of the genome along seven linkage groups with an average marker interval of ~7 cM. Eight traits, grain yield, lodging, height, heading data, grain protein, alpha amylase, diastatic power, and malt extract, were investigated in this project. Agronomic traits were measured in 16 areas, and malting quality traits in 9 areas. In our research, only three traits were studied, grain yield, height, and alpha amylase, and only the records in Crookston and Minnesota were used.

In the analysis, the prior expected number of QTL was taken as 3 for each trait, then the maximum number of QTL was calculated as L k ≈ 3 + 3· l k MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaOaaaeaacqWGSbaBaSqabaGcdaWgaaWcbaGaem4AaSgabeaaaaa@2EE6@ or L k = 8. Therefore, the prior inclusion probability of the model indicator variable equals to 0.375. To reduce the model space, we assumed each chromosome contain at most one QTL, except that the 7th was divided into two parts at the middle point and each part contains one QTL, for the results of other analysis (IM, CIM) always show signals of two QTL on 7th chromosome for some traits. Also two methods, multitrait analysis and Bayesian single-trait analysis (method 1 in [16]), were used to analyze the real data. The MCMC ran for 5 × 104 cycles after the first 2000 was discarded. The chain was thinned by every 10 cycles one observation being saved, which yielded 5000 samples for posterior Bayesian analysis.

Figure 5 and Figure 6 show the profiles of 2logeBF statistic with real data by multitrait analysis and separate analysis. The profiles of Figure 5 are generally higher than that of Figure 6. For trait 1 (grain yield), no QTL was detected by separate analysis (Figure 6a), while eight QTL were detected by multitrait analysis (Figure 5a); for trait 2 (height), three QTL located on chromosomes 1, 2, and 7 were detected by separate analysis, however by multitrait analysis, not only much stronger signals of these three QTL, but also four additional QTL on chromosome 3, 4, 5 and 6 were detected; for trait 3 (alpha amylase), two additional QTL located on chromosome 1, 3 were detected by multitrait analysis. The results of real data analysis also supported the conclusion that multitrait analysis was more powerful than separate analysis.

Figure 5
figure 5

The profiles of Bayes factors for multitrait trait analysis using real data. The profiles of Bayes factors (rescaled as 2 logeBF and negative values are truncated as zero) obtained from multitrait analysis using the real data: (a) for trait 1; (b) for trait 2; (c) for trait 3. The dotted vertical lines on the horizontal axis separate the chromosomes. The critical value is given as horizontal line. On the x-axis, inner tick marks represent markers.

Figure 6
figure 6

The profiles of Bayes factors for single trait analysis using real data. The profiles of Bayes factors (rescaled as 2 logeBF and negative values are truncated as zero) obtained from separate analysis using the real data: (a) for trait 1; (b) for trait 2; (c) for trait 3. The dotted vertical lines on the horizontal axis separate the chromosomes. The critical value is given as horizontal line. On the x-axis, inner tick marks represent markers.


The selection of hyper-parameter of the QTL effect is important in Bayesian analysis, which can influence the efficiency of the model selection. For example, with Bayesian shrinkage method [14], the hyper-parameter is a variable and assigned a special distribution so that no model selection is need. In Bayesian composite space approach, the updating of model indicator variables is closely dependent on QTL effects, but the selection of hyper-parameter is not much strict as Bayesian shrinkage analysis. Many approaches have been proposed for selection of hyper-parameter, and our method is only an extension of the approach of Yi et al. [15]. Moreover, we followed the approaches developed by Yi et al. [15] to obtain the prior probability for model indicator variables. However, we didn't investigate the influence of different prior probability on the results, because the proposed method is very computationally intensive. In addition, we suggested to use CIM-based multitrait analysis [2] to obtain the prior of variance-covariance of residual, but if prior information is not indeed known, we may take the noninformative prior [19], p ( Σ e ) Σ e 1 MathType@MTEF@5@5@+=feaagaart1ev2aqatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiCaaNaeiikaGccceGae83Odm1aaSbaaSqaaiabdwgaLbqabaGccqGGPaqkcqGHDisTcqWFJoWudaqhaaWcbaGaemyzaugabaGaeyOeI0IaeGymaedaaaaa@385E@ . In this simulation study, the noninformative prior is used and proved to be able to bring a precise estimate for variance-covariance of residual error.

The proposed multitrait analysis is based on Bayesian composite space approach, while other popular model selection approaches such as Bayesian shrinkage method [14] and Bayesian SSVS method [23] are also very easily extended, and the details will be demonstrated in another paper. We used BC and DH population as examples to demonstrate the efficiency of the method. The new method can be modified to be applied to other experiment designs, such as RIL, F2 design, etc. In addition, we only take the main effect into account, while the epistatic effect also can be included into the model. In that case, the model should be written as: y i = b 0 + q = 1 p Φ q X i q b q + q 1 < q 2 p Φ q 1 q 2 X i q 1 X i q 2 w q 1 q 2 + e i MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCyEaK3aaSbaaSqaaiabdMgaPbqabaGccqGH9aqpcqWHIbGydaWgaaWcbaGaeGimaadabeaakiabgUcaRmaaqahabaacceGae8NPdy0aaSbaaSqaaiabdghaXbqabaGccqWHybawdaWgaaWcbaGaemyAaKMaemyCaehabeaakiabhkgaInaaBaaaleaacqWGXbqCaeqaaaqaaiabdghaXjabg2da9iabigdaXaqaaiabdchaWbqdcqGHris5aOGaey4kaSYaaabCaeaacqWFMoGrdaWgaaWcbaGaemyCae3aaSbaaWqaaiabigdaXaqabaWccqWGXbqCdaWgaaadbaGaeGOmaidabeaaaSqabaGccqWHybawdaWgaaWcbaGaemyAaKMaemyCae3aaSbaaWqaaiabigdaXaqabaaaleqaaOGaeCiwaG1aaSbaaSqaaiabdMgaPjabdghaXnaaBaaameaacqaIYaGmaeqaaaWcbeaakiabhEha3naaBaaaleaacqWGXbqCdaWgaaadbaGaeGymaedabeaaliabdghaXnaaBaaameaacqaIYaGmaeqaaaWcbeaaaeaacqWGXbqCdaWgaaadbaGaeGymaedabeaaliabgYda8iabdghaXnaaBaaameaacqaIYaGmaeqaaaWcbaGaemiCaahaniabggHiLdGccqGHRaWkcqWHLbqzdaWgaaWcbaGaemyAaKgabeaaaaa@6B6B@ , where q is main effect, q1 and q2 is two interacting QTL, and w q 1 q 2 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeC4DaC3aaSbaaSqaaiabdghaXnaaBaaameaacqaIXaqmaeqaaSGaemyCae3aaSbaaWqaaiabikdaYaqabaaaleqaaaaa@32A4@ is (1 × m) column vectors of epistatic effect between QTL q1 and q2. Certainly, the implementation will be complicated and quite time-consuming, but nevertheless, the extension is feasible and expected to be very efficient for mapping interacting QTL.

In this paper, we have not given a test procedure to distinguish closely linked and pleiotropic QTL which cause the genetic correlations between each trait. There have been some of literatures about it, and generally, the likelihood ratio (LR) statistic [1, 2] and Bayesian factor (BF) statistic [7] always have been used to solve the problem [7]. In our multitrait analysis, although the LR testing procedure in [2] is completely applicable, it is not optimal, because it is based on single-QTL model. Also Bayesian approach can be used for such testing, but the computing time is a big factor of concern. Hopefully, an efficient and fast approach will be developed that could solve the problem nicely.


Bayesian composite space approach [18] is an effective method for model selection. Yi [16] firstly used it for QTL mapping and proved it to be effective for mapping multiple QTL. In this article, we extended this novel statistical method to multitrait mapping of QTL. Compared with separate analysis, joint analysis is optimal, because the parameters are updated by vector or matrix and the correlation information between multiple traits can be made good use of. The powerful of the proposed multitrait method also be proved by both simulation experiments and real data analysis, and they all showed that the multitrait analysis tends to give higher statistical power than the single trait analysis.


Multivariate linear model

Consider n individuals derived from a backcross population crossed from two inbred lines with observations on some densely distributed codominant markers and on m quantitative traits. Supposed that the maximum number of QTL is p, the phenotypic value y ki of individual i for k th trait can be described by the following multivariate linear model:

y k i = b k 0 + j = 1 p γ k j x k i j b k j + e k i , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemyEaK3aaSbaaSqaaiabdUgaRjabdMgaPbqabaGccqGH9aqpcqWGIbGydaWgaaWcbaGaem4AaSMaeGimaadabeaakiabgUcaRmaaqahabaGaeq4SdC2aaSbaaSqaaiabdUgaRjabdQgaQbqabaGccqWG4baEdaWgaaWcbaGaem4AaSMaemyAaKMaemOAaOgabeaakiabdkgaInaaBaaaleaacqWGRbWAcqWGQbGAaeqaaaqaaiabdQgaQjabg2da9iabigdaXaqaaiabdchaWbqdcqGHris5aOGaey4kaSIaemyzau2aaSbaaSqaaiabdUgaRjabdMgaPbqabaGccqGGSaalaaa@51DC@

for i = 1, 2, ..., n and k = 1, 2, ..., m, where γ kj is model indicator variable, indicating the j th QTL of k th trait included (1) or excluded (0) from the model; bk 0is population mean; b kj is QTL effect; x kij is QTL genotype, if QTL genotype is homozygote x kij = 1, otherwise -1; e ki is residual error and assumed to follow multivariate normal distribution. If we denote equation (1) by matrix, it can be expressed as:

y i = b 0 + j = 1 p Φ j X i j b j + e i , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCyEaK3aaSbaaSqaaiabdMgaPbqabaGccqGH9aqpcqWHIbGydaWgaaWcbaGaeGimaadabeaakiabgUcaRmaaqahabaacceGae8NPdy0aaSbaaSqaaiabdQgaQbqabaGccqWHybawdaWgaaWcbaGaemyAaKMaemOAaOgabeaakiabhkgaInaaBaaaleaacqWGQbGAaeqaaaqaaiabdQgaQjabg2da9iabigdaXaqaaiabdchaWbqdcqGHris5aOGaey4kaSIaeCyzau2aaSbaaSqaaiabdMgaPbqabaGccqGGSaalaaa@494E@

for i = 1, 2, ..., n, where y i = [y1i, y2i, ..., y mi ]T, b0 = [b10, b20, ..., bm 0]T, b j = [b1j, b2j, ..., b mj ]T, e i = [e1i, e2i, ..., e mi ]T. They are all (1 × m) column vectors. Equation (3) is QTL genotype matrix and Equation (4) is model indicator matrix, they are all (m × m) diagonal matrix.

X i j = [ x 1 i j 0 0 0 x 2 i j 0 0 0 x m i j ] MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCiwaG1aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqGH9aqpdaWadaqaauaabeqaeqaaaaaabaGaemiEaG3aaSbaaSqaaiabigdaXiabdMgaPjabdQgaQbqabaaakeaacqaIWaamaeaacqWIVlctaeaacqaIWaamaeaacqaIWaamaeaacqWG4baEdaWgaaWcbaGaeGOmaiJaemyAaKMaemOAaOgabeaaaOqaaiabl+UimbqaaiabicdaWaqaaiabl6Uinbqaaiabl6UinbqaaiablgVipbqaaiabl6UinbqaaiabicdaWaqaaiabicdaWaqaaiabl+UimbqaaiabdIha4naaBaaaleaacqWGTbqBcqWGPbqAcqWGQbGAaeqaaaaaaOGaay5waiaaw2faaaaa@56FB@
Φ j = [ γ 1 j 0 0 0 γ 2 j 0 0 0 γ m j ] MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae8NPdy0aaSbaaSqaaiabdQgaQbqabaGccqGH9aqpdaWadaqaauaabeqaeqaaaaaabaGaeq4SdC2aaSbaaSqaaiabigdaXiabdQgaQbqabaaakeaacqaIWaamaeaacqWIVlctaeaacqaIWaamaeaacqaIWaamaeaacqaHZoWzdaWgaaWcbaGaeGOmaiJaemOAaOgabeaaaOqaaiabl+UimbqaaiabicdaWaqaaiabl6Uinbqaaiabl6UinbqaaiablgVipbqaaiabl6UinbqaaiabicdaWaqaaiabicdaWaqaaiabl+Uimbqaaiabeo7aNnaaBaaaleaacqWGTbqBcqWGQbGAaeqaaaaaaOGaay5waiaaw2faaaaa@525B@

Prior specification

The prior distribution of each QTL effect vector b j is multivariate normal distribution, p(b j ) ~ N(0, Σ B j MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae83Odm1aaSbaaSqaaiabdkeacnaaBaaameaajugWaiabdQgaQbadbeaaaSqabaaaaa@3166@ ), where Σ B j MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae83Odm1aaSbaaSqaaiabdkeacnaaBaaameaajugWaiabdQgaQbadbeaaaSqabaaaaa@3166@ is the hyper-parameter, and We take Σ B j = [ X . j T Σ e 1 X . j ] 1 n MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae83Odm1aaSbaaSqaaiabdkeacnaaBaaameaajugWaiabdQgaQbadbeaaaSqabaGccqGH9aqpdaWadaqaaiabhIfayjabc6caUmaaDaaaleaacqWGQbGAaeaacqWGubavaaGccqWFJoWudaqhaaWcbaGaemyzaugabaGaeyOeI0IaeGymaedaaOGaeCiwaGLaeiOla4YaaSbaaSqaaiabdQgaQbqabaaakiaawUfacaGLDbaadaahaaWcbeqaaiabgkHiTiabigdaXaaakiabgwSixlabd6gaUbaa@47AA@ , which is simply an extension from Bayesian single trait analysis [15]. The importance of the choice of the hyper-parameter will be discussed later. In a large backcross population and under the definition of x mij (-1 or 1), Σ B j MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae83Odm1aaSbaaSqaaiabdkeacnaaBaaameaajugWaiabdQgaQbadbeaaaSqabaaaaa@3166@ can be simplified as Σ B j MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae83Odm1aaSbaaSqaaiabdkeacnaaBaaameaajugWaiabdQgaQbadbeaaaSqabaaaaa@3166@ = Σ e . The prior of the covariance matrix of residual error follows Inverse Wishart distribution, Σ e ~ Wishart-1(v e , S e 2 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeC4uam1aa0baaSqaaiabdwgaLbqaaiabikdaYaaaaaa@2F7A@ ), where, v e and S e 2 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeC4uam1aa0baaSqaaiabdwgaLbqaaiabikdaYaaaaaa@2F7A@ are prior degree of freedom and covariance matrix of residual error, respectively, and can be obtained from other method, such as CIM based multitrait analysis [2], etc. The prior distribution of population mean b0 is normal distribution with mean and variance equal to those calculated by phenotypic values. The prior probability distribution of QTL position λ kj is uniform distribution with bounds of two flanking markers, p(λ kj ) = 1/d j , where d j is length of the interval where j th QTL is confined. Assuming that epistatic effect is absent, the prior inclusion probability for j th effect can be expressed as p(γ kj = 1) = 1 - l k /L k ]1/N(see also [15]), where l k is the prior expected number of main-effect QTL, and could be roughly estimated with the use of standard genome scans; N is the number of possible main effects for each QTL and equal to 1 in BC family [15]; L k is the upper bound of QTL number, and equals to the number of marker interval in our simulation study, while in another approach suggested by Yi [15]L k is taken as 3 + 3· l k MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaWaaOaaaeaacqWGSbaBaSqabaGcdaWgaaWcbaGaem4AaSgabeaaaaa@2EE6@ , which causes the model space to reduce dramatically [15].

Joint posterior density

The observable variables include phenotypic values, y = { y i } i = 1 n MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCyEaKNaeyypa0Jaei4EaSNaeCyEaK3aaSbaaSqaaiabdMgaPbqabaGccqGG9bqFdaqhaaWcbaGaemyAaKMaeyypa0JaeGymaedabaGaemOBa4gaaaaa@394D@ and marker information, m = { m i j } i = 1 , j = 1 n , p MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCyBa0Maeyypa0Jaei4EaSNaemyBa02aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqGG9bqFdaqhaaWcbaGaemyAaKMaeyypa0JaeGymaeJaeiilaWIaemOAaOMaeyypa0JaeGymaedabaGaemOBa4MaeiilaWIaemiCaahaaaaa@40F2@ . The unobservable variables include population mean, b 0 = { b k 0 } k = 1 m MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCOyai2aaSbaaSqaaiabicdaWaqabaGccqGH9aqpcqGG7bWEcqWHIbGydaWgaaWcbaGaem4AaSMaeGimaadabeaakiabc2ha9naaDaaaleaacqWGRbWAcqGH9aqpcqaIXaqmaeaacqWGTbqBaaaaaa@3B09@ ; QTL effects, b = { b j } j = 1 p MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCOyaiMaeyypa0Jaei4EaSNaeCOyai2aaSbaaSqaaiabdQgaQbqabaGccqGG9bqFdaqhaaWcbaGaemOAaOMaeyypa0JaeGymaedabaGaemiCaahaaaaa@38F9@ ; QTL genotypes, X = { X i j } i = 1 , j = 1 n , p MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCiwaGLaeyypa0Jaei4EaSNaeCiwaG1aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqGG9bqFdaqhaaWcbaGaemyAaKMaeyypa0JaeGymaeJaeiilaWIaemOAaOMaeyypa0JaeGymaedabaGaemOBa4MaeiilaWIaemiCaahaaaaa@40A2@ ; model indicator variables, Φ = { Φ j } j = 1 p MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae8NPdyKaeyypa0Jaei4EaSNae8NPdy0aaSbaaSqaaiabdQgaQbqabaGccqGG9bqFdaqhaaWcbaGaemOAaOMaeyypa0JaeGymaedabaGaemiCaahaaaaa@394A@ ; (co)variance of residual error, Σ e , and QTL positions, λ = { λ k j } k = 1 , j = 1 m , p MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaaccmGae83UdWMaeyypa0Jaei4EaSNaeq4UdW2aaSbaaSqaaiabdUgaRjabdQgaQbqabaGccqGG9bqFdaqhaaWcbaGaem4AaSMaeyypa0JaeGymaeJaeiilaWIaemOAaOMaeyypa0JaeGymaedabaGaemyBa0MaeiilaWIaemiCaahaaaaa@419E@ . Let θ be the vector of hyper-parameters, Θ = {b0, b, Σ e , λ, X, Φ}, then the joint prior density of the unobservable variables is denoted by p(Θ|θ). The joint posterior probability of Θ, given the observable variables y and m, can be expressed as:

p(Θ|y, m) p(Θ|θp(y, m|Θ), (2)

where, p(y, m|Θ) is the likelihood and can be written as:

p(y, m|Θ) = p(y|Θp(m|Θ), (6)

where p(y|Θ) is multivariate normal density, and p(m|Θ) can be derived from a Markov model [14].

MCMC sampling

MCMC algorithm generates samples from Markov chains which converge to the posterior distribution of parameters, without the constant of proportionality being calculated. From these posterior samples, summary statistic of the posterior distribution can be calculated. MCMC algorithm proceeds as follows:

a. Initialize all parameters with values in their legal domain.

b. Update the population mean b0.

c. Update the QTL effects vectors { b j } j = 1 p MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaei4EaSNaeCOyai2aaSbaaSqaaiabdQgaQbqabaGccqGG9bqFdaqhaaWcbaGaemOAaOMaeyypa0JaeGymaedabaGaemiCaahaaaaa@36A2@ .

d. Update the variance-covariance matrix Σ e of the residual error.

e. Update the QTL genotype indicator matrices { X i j } i = 1 n MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaei4EaSNaeCiwaG1aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqGG9bqFdaqhaaWcbaGaemyAaKMaeyypa0JaeGymaedabaGaemOBa4gaaaaa@37E3@ and the QTL location vectors { λ k j } k = 1 m MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaei4EaSNaeq4UdW2aaSbaaSqaaiabdUgaRjabdQgaQbqabaGccqGG9bqFdaqhaaWcbaGaem4AaSMaeyypa0JaeGymaedabaGaemyBa0gaaaaa@3860@ jointly, for j = 1, 2,..., p.

f. Update the model indicator variable matrices { Φ j } j = 1 p MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaei4EaShcceGae8NPdy0aaSbaaSqaaiabdQgaQbqabaGccqGG9bqFdaqhaaWcbaGaemOAaOMaeyypa0JaeGymaedabaGaemiCaahaaaaa@36D0@ .

The conditional posterior distribution of the population mean b0 is multivariate normal with mean

b ¯ 0 = [ i = 1 n ( Σ e 1 ) ] 1 i = 1 n Σ e 1 ( y i j = 1 p Φ j X i j b j ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafCOyaiMbaebadaWgaaWcbaGaeGimaadabeaakiabg2da9maadmaabaWaaabCaeaacqGGOaakiiqacqWFJoWudaqhaaWcbaGaemyzaugabaGaeyOeI0IaeGymaedaaOGaeiykaKcaleaacqWGPbqAcqGH9aqpcqaIXaqmaeaacqWGUbGBa0GaeyyeIuoaaOGaay5waiaaw2faamaaCaaaleqabaGaeyOeI0IaeGymaedaaOWaaabCaeaacqWFJoWudaqhaaWcbaqcLbqacqWGLbqzaSqaaiabgkHiTiabigdaXaaakiabcIcaOiabhMha5naaBaaaleaacqWGPbqAaeqaaaqaaiabdMgaPjabg2da9iabigdaXaqaaiabd6gaUbqdcqGHris5aOGaeyOeI0YaaabCaeaacqWFMoGrdaWgaaWcbaGaemOAaOgabeaakiabhIfaynaaBaaaleaacqWGPbqAcqWGQbGAaeqaaOGaeCOyai2aaSbaaSqaaiabdQgaQbqabaaabaGaemOAaOMaeyypa0JaeGymaedabaGaemiCaahaniabggHiLdGccqGGPaqkcqGGSaalaaa@6544@

and variance-covariance matrix

Σ b 0 = [ i = 1 n ( Σ e 1 ) ] 1 . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae83Odmvcfa4aaSbaaSqaaKqzaeGaemOyai2cdaWgaaadbaqcLbmacqaIWaamaWqabaaaleqaaOGaeyypa0ZaamWaaeaadaaeWbqaaiabcIcaOiab=n6atnaaDaaaleaacqWGLbqzaeaacqGHsislcqaIXaqmaaGccqGGPaqkaSqaaiabdMgaPjabg2da9iabigdaXaqaaiabd6gaUbqdcqGHris5aaGccaGLBbGaayzxaaWaaWbaaSqabeaacqGHsislcqaIXaqmaaGccqGGUaGlaaa@4620@

The conditional posterior distribution of the QTL effect b j is sampled from multivariate normal distribution with mean

b ¯ j = [ Σ B 1 + i = 1 n ( X i j T Φ j T Σ e 1 Φ j X ) i j ] 1 i = 1 n X i j T Φ j T Σ e 1 ( y i j 1 p Φ j X i j b j b 0 ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGafCOyaiMbaebadaWgaaWcbaGaemOAaOgabeaakiabg2da9maadmaabaacceGae83Odm1aa0baaSqaaiabdkeacbqaaiabgkHiTiabigdaXaaakiabgUcaRmaaqahabaGaeiikaGIaeCiwaG1aa0baaSqaaiabdMgaPjabdQgaQbqaaiabdsfaubaakiab=z6agnaaDaaaleaacqWGQbGAaeaacqWGubavaaGccqWFJoWudaqhaaWcbaGaemyzaugabaGaeyOeI0IaeGymaedaaOGae8NPdy0aaSbaaSqaaiabdQgaQbqabaGccqWHybawdaWgbaWcbaGaemyAaKMaemOAaOgabeaakiabcMcaPaWcbaGaemyAaKMaeyypa0JaeGymaedabaGaemOBa4ganiabggHiLdaakiaawUfacaGLDbaadaahaaWcbeqaaiabgkHiTiabigdaXaaakmaaqahabaGaeCiwaG1aa0baaSqaaiabdMgaPjabdQgaQbqaaiabdsfaubaakiab=z6agnaaDaaaleaacqWGQbGAaeaacqWGubavaaGccqWFJoWudaqhaaWcbaGaemyzaugabaGaeyOeI0IaeGymaedaaOGaeiikaGIaeCyEaK3aaSbaaSqaaiabdMgaPbqabaGccqGHsislaSqaaiabdMgaPjabg2da9iabigdaXaqaaiabd6gaUbqdcqGHris5aOWaaabCaeaacqWFMoGrdaWgaaWcbaGaemOAaOgabeaakiabhIfaynaaBaaaleaacqWGPbqAcqWGQbGAaeqaaOGaeCOyai2aaSbaaSqaaiabdQgaQbqabaaabaGaemOAaOMaeyiyIKRaeGymaedabaGaemiCaahaniabggHiLdGccqGHsislcqWHIbGydaWgaaWcbaGaeGimaadabeaakiabcMcaPiabcYcaSaaa@8951@

and variance-covariance matrix

Σ b j = [ Σ B 1 + i = 1 n ( X i j T Φ j T Σ e 1 Φ j X ) i j ] 1 . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae83Odmvcfa4aaSbaaSqaaKqzaeGaemOyaiwcfa4aaSbaaWqaaKqzaeGaemOAaOgameqaaaWcbeaakiabg2da9maadmaabaGae83Odmvcfa4aa0baaSqaaiabdkeacbqaaKqzaeGaeyOeI0IaeGymaedaaOGaey4kaSYaaabCaeaacqGGOaakcqWHybawdaqhaaWcbaGaemyAaKMaemOAaOgabaGaemivaqfaaOGae8NPdy0aa0baaSqaaiabdQgaQbqaaiabdsfaubaakiab=n6atnaaDaaaleaacqWGLbqzaeaacqGHsislcqaIXaqmaaGccqWFMoGrdaWgaaWcbaGaemOAaOgabeaakiabhIfaynaaBeaaleaacqWGPbqAcqWGQbGAaeqaaOGaeiykaKcaleaacqWGPbqAcqGH9aqpcqaIXaqmaeaacqWGUbGBa0GaeyyeIuoaaOGaay5waiaaw2faamaaCaaaleqabaGaeyOeI0IaeGymaedaaOGaeiOla4caaa@5D9A@

The posterior distribution of the residual error follows inverted Wishart distribution,

Σ e ~ W i s h a r t 1 ( d f e + ν e , Ω T Ω + S e 2 ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae83Odm1aaSbaaSqaaiabdwgaLbqabaGccqGG+bGFcqWGxbWvcqWGPbqAcqWGZbWCcqWGObaAcqWGHbqycqWGYbGCcqWG0baDdaahaaWcbeqaaiabgkHiTiabigdaXaaakiabcIcaOiabdsgaKjabdAgaMnaaBaaaleaacqWGLbqzaeqaaOGaey4kaSIaeqyVd42aaSbaaSqaaiabdwgaLbqabaGccqGGSaalcqWFPoWvdaahaaWcbeqaaiabdsfaubaakiab=L6axjabgUcaRiabhofatnaaDaaaleaacqWGLbqzaeaacqaIYaGmaaGccqGGPaqkcqGGSaalaaa@511D@

where Ω = y i j = 1 p Φ j X i j b j b 0 MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae8xQdCLaeyypa0JaeCyEaK3aaSbaaSqaaiabdMgaPbqabaGccqGHsisldaaeWbqaaiab=z6agnaaBaaaleaacqWGQbGAaeqaaOGaeCiwaG1aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqWHIbGydaWgaaWcbaGaemOAaOgabeaaaeaacqWGQbGAcqGH9aqpcqaIXaqmaeaacqWGWbaCa0GaeyyeIuoakiabgkHiTiabhkgaInaaBaaaleaacqaIWaamaeqaaaaa@46CC@ and df e = n.

In step e, the QTL locations and QTL genotype matrices are updated jointly. For locus j, we can firstly sample a new QTL position for each trait from their prior distribution (described later), then sample the QTL genotype matrices { X i j } i = 1 n MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaei4EaSNaeCiwaG1aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqGG9bqFdaqhaaWcbaGaemyAaKMaeyypa0JaeGymaedabaGaemOBa4gaaaaa@37E3@ on the new position using equation (15), and finally, they are updated by the efficient Metropolis-Hastings algorithm [20, 21]. Because the sampling of X ij is too complicate and we are going to firstly describe it. Due to the QTL genotype x kij has two possible values (-1 or 1) in BC line, if m traits are investigated jointly, X ij has 2mkinds of possible formations, and the general pattern of X ij can be written as:

H i j , z 1 z 2 z m = [ x 1 i j = z 1 0 0 0 x 2 i j = z 2 0 0 0 x m i j = z m ] , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCisaG0aaSbaaSqaaKqzaeGaemyAaKMaemOAaOMaeiilaWIaemOEaOxcfa4aaSbaaWqaaKqzadGaeGymaedameqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGOmaidameqaaKqzaeGaeS47IWKaemOEaOxcfa4aaSbaaWqaaKqzadGaemyBa0gameqaaaWcbeaakiabg2da9maadmaabaqbaeqabqabaaaaaeaacqWG4baEdaWgaaWcbaGaeGymaeJaemyAaKMaemOAaOgabeaakiabg2da9iabdQha6naaBaaaleaacqaIXaqmaeqaaaGcbaGaeGimaadabaGaeS47IWeabaGaeGimaadabaGaeGimaadabaGaemiEaG3aaSbaaSqaaiabikdaYiabdMgaPjabdQgaQbqabaGccqGH9aqpcqWG6bGEdaWgaaWcbaGaeGOmaidabeaaaOqaaiabl+UimbqaaiabicdaWaqaaiabl6Uinbqaaiabl6UinbqaaiablgVipbqaaiabl6UinbqaaiabicdaWaqaaiabicdaWaqaaiabl+UimbqaaiabdIha4naaBaaaleaacqWGTbqBcqWGPbqAcqWGQbGAaeqaaOGaeyypa0JaemOEaO3aaSbaaSqaaiabd2gaTbqabaaaaaGccaGLBbGaayzxaaGaeiilaWcaaa@74EC@

where, z1, z2, ..., z m {-1, 1}. For clarity, we omit the subscript ij from H i j , z 1 z 2 z m MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCisaG0aaSbaaSqaaKqzaeGaemyAaKMaemOAaOMaeiilaWIaemOEaOxcfa4aaSbaaWqaaKqzadGaeGymaedameqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGOmaidameqaaKqzaeGaeS47IWKaemOEaOxcfa4aaSbaaWqaaKqzadGaemyBa0gameqaaaWcbeaaaaa@4197@ and present formulas H z 1 z 2 z m MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCisaG0aaSbaaSqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGymaedameqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGOmaidameqaaKqzaeGaeS47IWKaemOEaOxcfa4aaSbaaWqaaKqzadGaemyBa0gameqaaaWcbeaaaaa@3DFF@ to denote the genotype matrix of i th individual and j th loci. Because the QTL genotypes x kij of i th individual in the j th interval for all traits may be correlated, the joint prior probability of the genotype matrix X ij can't be simply expressed by the following equation:

p ( X i j = H z 1 z 2 z m | λ j , m i , j , m i , j + 1 ) = p ( x 1 i j = z 1 , x 2 i j = z 2 , , x m i j = z m | λ j , m i , j , m i , j + 1 ) = k = 1 m p ( x k i j = z k | m i , j , m i , j + 1 ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbiqaaGHefaqabeGabaaabaGaemiCaaNaeiikaGIaeCiwaG1aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqGH9aqpcqWHibasdaWgaaWcbaqcLbqacqWG6bGEjuaGdaWgaaadbaqcLbmacqaIXaqmaWqabaqcLbqacqWG6bGEjuaGdaWgaaadbaqcLbmacqaIYaGmaWqabaqcLbqacqWIVlctcqWG6bGEjuaGdaWgaaadbaqcLbmacqWGTbqBaWqabaaaleqaaOWaaqqaaeaaiiqacqWF7oaBdaWgaaWcbaGaemOAaOgabeaakiabcYcaSiabd2gaTnaaBaaaleaacqWGPbqAcqGGSaalcqWGQbGAaeqaaOGaeiilaWIaemyBa02aaSbaaSqaaiabdMgaPjabcYcaSiabdQgaQjabgUcaRiabigdaXaqabaaakiaawEa7aiabcMcaPiabg2da9iabdchaWjabcIcaOiabdIha4naaBaaaleaacqaIXaqmcqWGPbqAcqWGQbGAaeqaaOGaeyypa0JaemOEaO3aaSbaaSqaaiabigdaXaqabaGccqGGSaalcqWG4baEdaWgaaWcbaGaeGOmaiJaemyAaKMaemOAaOgabeaakiabg2da9iabdQha6naaBaaaleaacqaIYaGmaeqaaOGaeiilaWIaeS47IWKaeiilaWIaemiEaG3aaSbaaSqaaiabd2gaTjabdMgaPjabdQgaQbqabaGccqGH9aqpcqWG6bGEdaWgaaWcbaGaemyBa0gabeaakmaaeeaabaGae83UdW2aaSbaaSqaaiabdQgaQbqabaGccqGGSaalcqWGTbqBdaWgaaWcbaGaemyAaKMaeiilaWIaemOAaOgabeaakiabcYcaSiabd2gaTnaaBaaaleaacqWGPbqAcqGGSaalcqWGQbGAcqGHRaWkcqaIXaqmaeqaaaGccaGLhWoacqGGPaqkaeGabaqegiaaxMaacqGH9aqpdaqeWbqaaiabdchaWjabcIcaOiabdIha4naaBaaaleaacqWGRbWAcqWGPbqAcqWGQbGAaeqaaOGaeyypa0JaemOEaO3aaSbaaSqaaiabdUgaRbqabaGcdaabbaqaaiabd2gaTnaaBaaaleaacqWGPbqAcqGGSaalcqWGQbGAaeqaaOGaeiilaWIaemyBa02aaSbaaSqaaiabdMgaPjabcYcaSiabdQgaQjabgUcaRiabigdaXaqabaaakiaawEa7aiabcMcaPaWcbaGaem4AaSMaeyypa0JaeGymaedabaGaemyBa0ganiabg+GivdaaaOGaeiOla4caaa@B64F@

Instead, it can be derived from the Markov model (see Equation 14), assuming that the order of markers and QTL is M j Q1Q2 ... Q m Mj+1(see Figure 7), where, Q1, Q2, ..., and Q m denote the QTL respectively affecting trait 1, trait 2, ..., and trait m in j th marker interval. Indicator variables x1ij, x2ij, ..., and x mij denote the genotypes of these QTL.

Figure 7
figure 7

The positions of markers and QTL and their sequence ranged on a certain marker interval.

p ( X i j = H z 1 z 2 z m | m i , j , λ j , m i , j + 1 ) = p ( x 1 i j = z 1 , x 2 i j = z 2 , , x m i j = z m | m i , j , λ j , m i , j + 1 ) = p ( x 1 i j = z 1 | m i , j , λ 1 j , m i , j + 1 ) p ( x 2 i j = z 2 | m i , j , λ 2 j , x 1 i j , m i , j + 1 ) × p ( x m i j = z m | m i , j , x 1 i j , x 2 i j , , x ( m 1 ) i j , λ m j , m i , j + 1 ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGceaqabeaacqWGWbaCcqGGOaakcqWHybawdaWgaaWcbaGaemyAaKMaemOAaOgabeaakiabg2da9iabhIeainaaBaaaleaajugabiabdQha6LqbaoaaBaaameaajugWaiabigdaXaadbeaajugabiabdQha6LqbaoaaBaaameaajugWaiabikdaYaadbeaajugabiabl+UimjabdQha6LqbaoaaBaaameaajugWaiabd2gaTbadbeaaaSqabaGcdaabbaqaaiabd2gaTnaaBaaaleaacqWGPbqAcqGGSaalcqWGQbGAaeqaaOGaeiilaWccceGae83UdW2aaSbaaSqaaiabdQgaQbqabaGccqGGSaalcqWGTbqBdaWgaaWcbaGaemyAaKMaeiilaWIaemOAaOMaey4kaSIaeGymaedabeaaaOGaay5bSdGaeiykaKIaeyypa0JaemiCaaNaeiikaGIaemiEaG3aaSbaaSqaaiabigdaXiabdMgaPjabdQgaQbqabaGccqGH9aqpcqWG6bGEdaWgaaWcbaGaeGymaedabeaakiabcYcaSiabdIha4naaBaaaleaacqaIYaGmcqWGPbqAcqWGQbGAaeqaaOGaeyypa0JaemOEaO3aaSbaaSqaaiabikdaYaqabaGccqGGSaalcqWIVlctcqGGSaalcqWG4baEdaWgaaWcbaGaemyBa0MaemyAaKMaemOAaOgabeaakiabg2da9iabdQha6naaBaaaleaacqWGTbqBaeqaaOWaaqqaaeaacqWGTbqBdaWgaaWcbaGaemyAaKMaeiilaWIaemOAaOgabeaakiabcYcaSiab=T7aSnaaBaaaleaacqWGQbGAaeqaaOGaeiilaWIaemyBa02aaSbaaSqaaiabdMgaPjabcYcaSiabdQgaQjabgUcaRiabigdaXaqabaaakiaawEa7aiabcMcaPaqaauaabeqaceaaaeGabaalqiaaxMaacqGH9aqpcqWGWbaCcqGGOaakcqWG4baEdaWgaaWcbaGaeGymaetcLbqacqWGPbqAcqWGQbGAaSqabaGccqGH9aqpcqWG6bGEdaWgaaWcbaGaeGymaedabeaakmaaeeaabaGaemyBa02aaSbaaSqaaiabdMgaPjabcYcaSiabdQgaQbqabaGccqGGSaalcaa57oWaaSbaaSqaaiabigdaXiabdQgaQbqabaGccqGGSaalcqWGTbqBdaWgaaWcbaGaemyAaKMaeiilaWIaemOAaOMaey4kaSIaeGymaedabeaaaOGaay5bSdGaeiykaKIaeyyXICTaemiCaaNaeiikaGIaemiEaG3aaSbaaSqaaiabikdaYKqzaeGaemyAaKMaemOAaOgaleqaaOGaeyypa0JaemOEaO3aaSbaaSqaaiabikdaYaqabaGcdaabbaqaaiabd2gaTnaaBaaaleaacqWGPbqAcqGGSaalcqWGQbGAaeqaaOGaeiilaWIaaq+UdmaaBaaaleaacqaIYaGmcqWGQbGAaeqaaOGaeiilaWIaemiEaG3aaSbaaSqaaiabigdaXiabdMgaPjabdQgaQbqabaaakiaawEa7aiabcYcaSiabd2gaTnaaBaaaleaacqWGPbqAcqGGSaalcqWGQbGAcqGHRaWkcqaIXaqmaeqaaOGaeiykaKcabiWaai+baycdaORdcaWLjaGaaCzcaiaaxMaacqGHxdaTcqGHflY1cqWIVlctcqGHflY1cqWGWbaCcqGGOaakcqWG4baEdaWgaaWcbaGaemyBa0wcLbqacqWGPbqAcqWGQbGAaSqabaGccqGH9aqpcqWG6bGEdaWgaaWcbaGaemyBa0gabeaakmaaeeaabaGaemyBa02aaSbaaSqaaiabdMgaPjabcYcaSiabdQgaQbqabaGccqGGSaalcqWG4baEdaWgaaWcbaGaeGymaetcLbqacqWGPbqAcqWGQbGAaSqabaGccqGGSaalcqWG4baEdaWgaaWcbaGaeGOmaitcLbqacqWGPbqAcqWGQbGAaSqabaGccqGGSaalcqWIVlctcqGGSaalcqWG4baEdaWgaaWcbaqcLboacqGGOaakcqWGTbqBcqGHsislcqaIXaqmcqGGPaqkcqWGPbqAcqWGQbGAaSqabaaakiaawEa7aiabcYcaSiaaKV7adaWgaaWcbaGaemyBa0MaemOAaOgabeaakiabcYcaSiabd2gaTnaaBaaaleaacqWGPbqAcqGGSaalcqWGQbGAcqGHRaWkcqaIXaqmaeqaaOGaeiykaKIaeiilaWcaaaaaaa@268D@

If no segregation interference is considered, the joint prior probability can be factorized into equation (14), and each term in equation (14) can be derived from Haldane map function. Only the first term in equation (14) is conditional on two flanking markers; others are not only conditional on two flanking markers but also on the genotypes of all the QTL prior to the interested one. If double recombination is ignored [2], each term in equation (14) can be inferred only by the genotype of the left nearest loci (marker or QTL) and the right marker, then equation (14) can be simplified as:

p ( X i j = H z 1 z 2 z m | m i , j , λ j , m i , j + 1 ) = p ( x 1 i j = z 1 , x 2 i j = z 2 , , x m i j = z m | m i , j , λ j , m i , j + 1 ) = p ( x 1 i j = z 1 | m i , j , λ 1 j , m i , j + 1 ) p ( x 2 i j = z 2 | x 1 i j , λ 2 j , m i , j + 1 ) × p ( x m i j = z m | x ( m 1 ) i j , λ m j , m i , j + 1 ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGceiqabiaayrbauvraaiabdchaWjabcIcaOiabhIfaynaaBaaaleaacqWGPbqAcqWGQbGAaeqaaOGaeyypa0JaeCisaG0aaSbaaSqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGymaedameqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGOmaidameqaaKqzaeGaeS47IWKaemOEaOxcfa4aaSbaaWqaaKqzadGaemyBa0gameqaaaWcbeaakmaaeeaabaGaemyBa02aaSbaaSqaaiabdMgaPjabcYcaSiabdQgaQbqabaGccqGGSaaliiqacqWF7oaBdaWgaaWcbaGaemOAaOgabeaakiabcYcaSiabd2gaTnaaBaaaleaacqWGPbqAcqGGSaalcqWGQbGAcqGHRaWkcqaIXaqmaeqaaaGccaGLhWoacqGGPaqkcqGH9aqpcqWGWbaCcqGGOaakcqWG4baEdaWgaaWcbaGaeGymaeJaemyAaKMaemOAaOgabeaakiabg2da9iabdQha6naaBaaaleaacqaIXaqmaeqaaOGaeiilaWIaemiEaG3aaSbaaSqaaiabikdaYiabdMgaPjabdQgaQbqabaGccqGH9aqpcqWG6bGEdaWgaaWcbaGaeGOmaidabeaakiabcYcaSiabl+UimjabcYcaSiabdIha4naaBaaaleaacqWGTbqBcqWGPbqAcqWGQbGAaeqaaOGaeyypa0JaemOEaO3aaSbaaSqaaiabd2gaTbqabaGcdaabbaqaaiabd2gaTnaaBaaaleaacqWGPbqAcqGGSaalcqWGQbGAaeqaaOGaeiilaWIae83UdW2aaSbaaSqaaiabdQgaQbqabaGccqGGSaalcqWGTbqBdaWgaaWcbaGaemyAaKMaeiilaWIaemOAaOMaey4kaSIaeGymaedabeaaaOGaay5bSdGaeiykaKcabaGaaCzcaiabg2da9iabdchaWjabcIcaOiabdIha4naaBaaaleaacqaIXaqmjugabiabdMgaPjabdQgaQbWcbeaakiabg2da9iabdQha6naaBaaaleaacqaIXaqmaeqaaOWaaqqaaeaacqWGTbqBdaWgaaWcbaGaemyAaKMaeiilaWIaemOAaOgabeaakiabcYcaSiaaKV7adaWgaaWcbaGaeGymaeJaemOAaOgabeaakiabcYcaSiabd2gaTnaaBaaaleaacqWGPbqAcqGGSaalcqWGQbGAcqGHRaWkcqaIXaqmaeqaaaGccaGLhWoacqGGPaqkcqGHflY1cqWGWbaCcqGGOaakcqWG4baEdaWgaaWcbaGaeGOmaitcLbqacqWGPbqAcqWGQbGAaSqabaGccqGH9aqpcqWG6bGEdaWgaaWcbaGaeGOmaidabeaakmaaeeaabaGaemiEaG3aaSbaaSqaaiabigdaXiabdMgaPjabdQgaQbqabaaakiaawEa7aiabcYcaSiaaKV7adaWgaaWcbaGaeGOmaiJaemOAaOgabeaakiabcYcaSiabd2gaTnaaBaaaleaacqWGPbqAcqGGSaalcqWGQbGAcqGHRaWkcqaIXaqmaeqaaOGaeiykaKcabaGaaCzcaiaaxMaacqGHxdaTcqGHflY1cqWIVlctcqGHflY1cqWGWbaCcqGGOaakcqWG4baEdaWgaaWcbaGaemyBa0wcLbqacqWGPbqAcqWGQbGAaSqabaGccqGH9aqpcqWG6bGEdaWgaaWcbaGaemyBa0gabeaakmaaeeaabaGaemiEaG3aaSbaaSqaaKqzGdGaeiikaGIaemyBa0MaeyOeI0IaeGymaeJaeiykaKIaemyAaKMaemOAaOgaleqaaaGccaGLhWoacqGGSaalcaa57oWaaSbaaSqaaiabd2gaTjabdQgaQbqabaGccqGGSaalcqWGTbqBdaWgaaWcbaGaemyAaKMaeiilaWIaemOAaOMaey4kaSIaeGymaedabeaakiabcMcaPiabcYcaSaaaaa@083D@

Each term in equation (15) can be easily inferred.

It is worth mentioning that we assume the sequence of markers and QTL is M j Q1Q2 ... Q m Mj+1, and in fact, the sequence of QTL may be variable in each round of updating. Therefore, we should firstly ascertain the sequence in each round, and then construct the appropriate formula to calculate the joint prior probability of the QTL genotype p(X ij = H z 1 z 2 z m MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCisaG0aaSbaaSqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGymaedameqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGOmaidameqaaKqzaeGaeS47IWKaemOEaOxcfa4aaSbaaWqaaKqzadGaemyBa0gameqaaaWcbeaaaaa@3DFF@ |mi,j,λ j,mi,j+1) according above rules. For clarity, we take an example to demonstrate it. Consider 3 QTL Q1, Q2, and Q3 that affect 3 traits respectively in an interval. Assuming that in a certain round the sequence of markers and QTL is M j Q3Q1Q2Mj+1, then the formula for calculating the joint prior probability of the QTL genotype can be written as:

p ( X i j = H z 1 z 2 z 3 | m i , j , λ j , m i , j + 1 ) = p ( x 1 i j = z 1 , x 2 i j = z 2 , x 3 i j = z 3 | m i , j , λ j , m i , j + 1 ) = p ( x 3 i j = z 3 | m i , j , λ 3 j , m i , j + 1 ) p ( x 1 i j = z 1 | x 3 i j , λ 1 j , m i , j + 1 ) × p ( x 2 i j = z 2 | x 1 i j , λ 2 j , m i , j + 1 ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGceiqabeaaumraaiabdchaWjabcIcaOiabhIfaynaaBaaaleaacqWGPbqAcqWGQbGAaeqaaOGaeyypa0JaeCisaG0aaSbaaSqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGymaedameqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGOmaidameqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeG4mamdameqaaaWcbeaakmaaeeaabaGaemyBa02aaSbaaSqaaiabdMgaPjabcYcaSiabdQgaQbqabaGccqGGSaaliiqacqWF7oaBdaWgaaWcbaGaemOAaOgabeaakiabcYcaSiabd2gaTnaaBaaaleaacqWGPbqAcqGGSaalcqWGQbGAcqGHRaWkcqaIXaqmaeqaaaGccaGLhWoacqGGPaqkcqGH9aqpcqWGWbaCcqGGOaakcqWG4baEdaWgaaWcbaGaeGymaeJaemyAaKMaemOAaOgabeaakiabg2da9iabdQha6naaBaaaleaacqaIXaqmaeqaaOGaeiilaWIaemiEaG3aaSbaaSqaaiabikdaYiabdMgaPjabdQgaQbqabaGccqGH9aqpcqWG6bGEdaWgaaWcbaGaeGOmaidabeaakiabcYcaSiabdIha4naaBaaaleaacqaIZaWmcqWGPbqAcqWGQbGAaeqaaOGaeyypa0JaemOEaO3aaSbaaSqaaiabiodaZaqabaGcdaabbaqaaiabd2gaTnaaBaaaleaacqWGPbqAcqGGSaalcqWGQbGAaeqaaOGaeiilaWIae83UdW2aaSbaaSqaaiabdQgaQbqabaGccqGGSaalcqWGTbqBdaWgaaWcbaGaemyAaKMaeiilaWIaemOAaOMaey4kaSIaeGymaedabeaaaOGaay5bSdGaeiykaKcabaqbaeqabiqaaaqaceaaumHaaCzcaiabg2da9iabdchaWjabcIcaOiabdIha4naaBaaaleaajugabiabiodaZiabdMgaPjabdQgaQbWcbeaakiabg2da9iabdQha6naaBaaaleaacqaIZaWmaeqaaOWaaqqaaeaacqWGTbqBdaWgaaWcbaGaemyAaKMaeiilaWIaemOAaOgabeaakiabcYcaSGGaaiab+T7aSnaaBaaaleaacqaIZaWmcqWGQbGAaeqaaOGaeiilaWIaemyBa02aaSbaaSqaaiabdMgaPjabcYcaSiabdQgaQjabgUcaRiabigdaXaqabaaakiaawEa7aiabcMcaPiabgwSixlabdchaWjabcIcaOiabdIha4naaBaaaleaajugabiabigdaXiabdMgaPjabdQgaQbWcbeaakiabg2da9iabdQha6naaBaaaleaacqaIXaqmaeqaaOWaaqqaaeaacqWG4baEdaWgaaWcbaGaeG4mamJaemyAaKMaemOAaOgabeaaaOGaay5bSdGaeiilaWIae43UdW2aaSbaaSqaaiabigdaXiabdQgaQbqabaGccqGGSaalcqWGTbqBdaWgaaWcbaGaemyAaKMaeiilaWIaemOAaOMaey4kaSIaeGymaedabeaakiabcMcaPaqaciaaRTpaozGaaCzcaiabgEna0kabdchaWjabcIcaOiabdIha4naaBaaaleaajugabiabikdaYiabdMgaPjabdQgaQbWcbeaakiabg2da9iabdQha6naaBaaaleaacqaIYaGmaeqaaOWaaqqaaeaacqWG4baEdaWgaaWcbaqcLboacqaIXaqmcqWGPbqAcqWGQbGAaSqabaaakiaawEa7aiabcYcaSiab+T7aSnaaBaaaleaacqaIYaGmcqWGQbGAaeqaaOGaeiilaWIaemyBa02aaSbaaSqaaiabdMgaPjabcYcaSiabdQgaQjabgUcaRiabigdaXaqabaGccqGGPaqkcqGGUaGlaaaaaaa@F736@

Once we obtain the joint prior probability of the QTL genotype, the joint conditional posterior probability of X ij can be expressed as:

p ( X i j = H z 1 z 2 z m | y i , ) = f ( y i | X i j = H z 1 z 2 z m , ) p ( X i j = H z 1 z 2 z m | λ j , m i j , m i , j + 1 ) h 1 = 1 1 h 2 = 1 1 h m = 1 1 f ( y i | X i j = H h 1 h 2 h m , ) p ( X i j = H h 1 h 2 h m | λ j , m i j , m i , j + 1 ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiCaaNaeiikaGIaeCiwaG1aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqGH9aqpcqWHibasdaWgaaWcbaqcLbqacqWG6bGEjuaGdaWgaaadbaqcLbmacqaIXaqmaWqabaqcLbqacqWG6bGEjuaGdaWgaaadbaqcLbmacqaIYaGmaWqabaqcLbqacqWIVlctcqWG6bGEjuaGdaWgaaadbaqcLbmacqWGTbqBaWqabaaaleqaaOWaaqqaaeaacqWH5bqEdaWgaaWcbaGaemyAaKgabeaakiabcYcaSiabl+UimbGaay5bSdGaeiykaKIaeyypa0tcfa4aaSaaaeaacqWGMbGzcqGGOaakcqWH5bqEdaWgaaqaaiabdMgaPbqabaWaaqqaaeaacqWHybawdaWgaaqaaiabdMgaPjabdQgaQbqabaGaeyypa0JaeCisaG0aaSbaaeaacqWG6bGEdaWgaaqaaiabigdaXaqabaGaemOEaO3aaSbaaeaacqaIYaGmaeqaaiabl+UimjabdQha6naaBaaabaGaemyBa0gabeaaaeqaaiabcYcaSiabl+UimbGaay5bSdGaeiykaKIaemiCaaNaeiikaGIaeCiwaG1aaSbaaeaacqWGPbqAcqWGQbGAaeqaaiabg2da9iabhIeainaaBaaabaGaemOEaO3aaSbaaeaacqaIXaqmaeqaaiabdQha6naaBaaabaGaeGOmaidabeaacqWIVlctcqWG6bGEdaWgaaqaaiabd2gaTbqabaaabeaadaabbaqaaGGabiab=T7aSnaaBaaabaGaemOAaOgabeaacqGGSaalcqWGTbqBdaWgaaqaaiabdMgaPjabdQgaQbqabaGaeiilaWIaemyBa02aaSbaaeaacqWGPbqAcqGGSaalcqWGQbGAcqGHRaWkcqaIXaqmaeqaaaGaay5bSdGaeiykaKcabaWaaabmaeaadaaeWaqaaiabl+UimnaaqadabaGaemOzayMaeiikaGIaeCyEaK3aaSbaaeaacqWGPbqAaeqaamaaeeaabaGaeCiwaG1aaSbaaeaacqWGPbqAcqWGQbGAaeqaaiabg2da9iabhIeainaaBaaabaGaemiAaG2aaSbaaeaacqaIXaqmaeqaaiabdIgaOnaaBaaabaGaeGOmaidabeaacqWIVlctcqWGObaAdaWgaaqaaiabd2gaTbqabaaabeaacqGGSaalcqWIVlctaiaawEa7aiabcMcaPiabdchaWjabcIcaOiabhIfaynaaBaaabaGaemyAaKMaemOAaOgabeaacqGH9aqpcqWHibasdaWgaaqaaiabdIgaOnaaBaaabaGaeGymaedabeaacqWGObaAdaWgaaqaaiabikdaYaqabaGaeS47IWKaemiAaG2aaSbaaeaacqWGTbqBaeqaaaqabaWaaqqaaeaacqWF7oaBdaWgaaqaaiabdQgaQbqabaGaeiilaWIaemyBa02aaSbaaeaacqWGPbqAcqWGQbGAaeqaaiabcYcaSiabd2gaTnaaBaaabaGaemyAaKMaeiilaWIaemOAaOMaey4kaSIaeGymaedabeaaaiaawEa7aiabcMcaPaqaaiabdIgaOnaaBaaabaGaemyBa0gabeaacqGH9aqpcqGHsislcqaIXaqmaeaacqaIXaqmaiabggHiLdaabaGaemiAaG2aaSbaaeaacqaIYaGmaeqaaiabg2da9iabgkHiTiabigdaXaqaaiabigdaXaGaeyyeIuoaaeaacqWGObaAdaWgaaqaaiabigdaXaqabaGaeyypa0JaeyOeI0IaeGymaedabaGaeGymaedacqGHris5aaaaaaa@E9B4@

where f ( y i | X i j = H z 1 z 2 z m , ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOzayMaeiikaGIaeCyEaK3aaSbaaSqaaiabdMgaPbqabaGcdaabbaqaaiabhIfaynaaBaaaleaacqWGPbqAcqWGQbGAaeqaaOGaeyypa0JaeCisaG0aaSbaaSqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGymaedameqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGOmaidameqaaKqzaeGaeS47IWKaemOEaOxcfa4aaSbaaWqaaKqzadGaemyBa0gameqaaaWcbeaakiabcYcaSiabl+UimbGaay5bSdGaeiykaKcaaa@4DB3@ is likelihood, and follows multivariable normal distribution,

f ( y i | X i j = H z 1 z 2 z m , ) = 1 ( 2 π ) m / 2 | Σ e | 1 / 2 exp { 1 2 ( y i j = 1 p Φ j X i j b j b 0 ) T Σ e 1 ( y i j = 1 p Φ j X i j b j b 0 ) } MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemOzayMaeiikaGIaeCyEaK3aaSbaaSqaaiabdMgaPbqabaGcdaabbaqaaiabhIfaynaaBaaaleaacqWGPbqAcqWGQbGAaeqaaOGaeyypa0JaeCisaG0aaSbaaSqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGymaedameqaaKqzaeGaemOEaOxcfa4aaSbaaWqaaKqzadGaeGOmaidameqaaKqzaeGaeS47IWKaemOEaOxcfa4aaSbaaWqaaKqzadGaemyBa0gameqaaaWcbeaakiabcYcaSiabl+UimbGaay5bSdGaeiykaKIaeyypa0tcfa4aaSaaaeaacqaIXaqmaeaacqGGOaakcqaIYaGmcqaHapaCcqGGPaqkdaahaaqabeaadaWcgaqaaiabd2gaTbqaaiabikdaYaaaaaWaaqWaaeaaiiqacqWFJoWudaWgaaqaaiabdwgaLbqabaaacaGLhWUaayjcSdWaaWbaaeqabaGaeGymaeJaei4la8IaeGOmaidaaaaakiGbcwgaLjabcIha4jabcchaWnaacmaabaGaeyOeI0scfa4aaSaaaeaacqaIXaqmaeaacqaIYaGmaaGccqGGOaakcqWH5bqEdaWgaaWcbaGaemyAaKgabeaakiabgkHiTmaaqahabaGae8NPdy0aaSbaaSqaaiabdQgaQbqabaGccqWHybawdaWgaaWcbaGaemyAaKMaemOAaOgabeaakiabhkgaInaaBaaaleaacqWGQbGAaeqaaaqaaiabdQgaQjabg2da9iabigdaXaqaaiabdchaWbqdcqGHris5aOGaeyOeI0IaeCOyai2aaSbaaSqaaiabicdaWaqabaGccqGGPaqkdaahaaWcbeqaaiabdsfaubaakiab=n6atnaaDaaaleaacqWGLbqzaeaacqGHsislcqaIXaqmaaGccqGGOaakcqWH5bqEdaWgaaWcbaGaemyAaKgabeaakiabgkHiTmaaqahabaGae8NPdy0aaSbaaSqaaiabdQgaQbqabaGccqWHybawdaWgaaWcbaGaemyAaKMaemOAaOgabeaakiabhkgaInaaBaaaleaacqWGQbGAaeqaaOGaeyOeI0IaeCOyai2aaSbaaSqaaiabicdaWaqabaaabaGaemOAaOMaeyypa0JaeGymaedabaGaemiCaahaniabggHiLdGccqGGPaqkaiaawUhacaGL9baaaaa@A4F3@

Once we have calculated 2mpossible posterior probabilities for the corresponding QTL genotype matrices, we are going to sample one genotype matrix according to their posterior probabilities. We firstly constructed the cumulative probability function F(d) by accumulating the 2mprobabilities in an arbitrary sequence for d = 1, 2, ..., 2mand F(0) = 0, which is a discrete distribution; then sampled a random number from uniform distribution, u ~ U[0,1]; and compared u with F(d), if F(d - 1) <uF(d), then the d th genotype matrix is accepted.

The new sampled QTL genotype matrices { X i j } i = 1 n MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaei4EaSNaeCiwaG1aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqGG9bqFdaqhaaWcbaGaemyAaKMaeyypa0JaeGymaedabaGaemOBa4gaaaaa@37E3@ are only the proposal value, which should be updated along with the proposal QTL position vector λ j = [λ1j, λ2j, ..., λ mj ] by the Metropolis-Hastings algorithm [20, 21]. For each trait, the new proposal position is sampled around the existing one from uniform distributions, λ k j MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeq4UdW2aa0baaSqaaiabdUgaRjabdQgaQbqaaiabgEHiQaaaaaa@3161@ ~ [λ kj - δ, λ kj + δ), where δ is tuning parameter, usually taking a value of 1 or 2 cM. The new position vector is denoted by λ j = [ λ 1 j , λ 2 j , , λ m j ] MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae83UdW2aa0baaSqaaiabdQgaQbqaaiabgEHiQaaakiabg2da9maadmaabaGaeq4UdW2aa0baaSqaaiabigdaXiabdQgaQbqaaiabgEHiQaaakiabcYcaSiabeU7aSnaaDaaaleaacqaIYaGmcqWGQbGAaeaacqGHxiIkaaGccqGGSaalcqWIVlctcqGGSaalcqaH7oaBdaqhaaWcbaGaemyBa0MaemOAaOgabaGaey4fIOcaaaGccaGLBbGaayzxaaaaaa@4782@ ; then the new QTL genotype matrix X i j MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeCiwaG1aa0baaSqaaiabdMgaPjabdQgaQbqaaiabgEHiQaaaaaa@30E6@ is sampled conditionally on the new position using equation (16); finally, the position vector λ j MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae83UdW2aa0baaSqaaiabdQgaQbqaaiabgEHiQaaaaaa@3008@ and genotype matrices { X i j } i = 1 n MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaei4EaSNaeCiwaG1aaSbaaSqaaiabdMgaPjabdQgaQbqabaGccqGG9bqFdaqhaaWcbaGaemyAaKMaeyypa0JaeGymaedabaGaemOBa4gaaaaa@37E3@ are accepted jointly with probability equal to min(1,α), where

α = i = 1 n p ( y i | X i j , λ j , ) p ( X i j | λ j , ) p ( λ j ) i = 1 n p ( y i | X i j , λ j , ) p ( X i j | λ j , ) p ( λ j ) q ( X i j | y i , ) q ( λ j ) q ( X i j | y i , ) q ( λ j ) , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeqySdeMaeyypa0tcfa4aaSaaaeaadaqeWaqaaiabdchaWjabcIcaOiabhMha5naaBaaabaGaemyAaKgabeaadaabbaqaaiabhIfaynaaDaaabaGaemyAaKMaemOAaOgabaGaey4fIOcaaiabcYcaSGGabiab=T7aSnaaDaaabaGaemOAaOgabaGaey4fIOcaaiabcYcaSaGaay5bSdGaeS47IWKaeiykaKIaemiCaaNaeiikaGIaeCiwaG1aa0baaeaacqWGPbqAcqWGQbGAaeaacqGHxiIkaaWaaqqaaeaacqWF7oaBdaqhaaqaaiabdQgaQbqaaiabgEHiQaaaaiaawEa7aiabcYcaSiabl+UimjabcMcaPiabdchaWjabcIcaOiabZV7aSnaaDaaabaGaemOAaOgabaGaey4fIOcaaiabcMcaPaqaaiabdMgaPjabg2da9iabigdaXaqaaiabd6gaUbGaey4dIunaaeaadaqeWaqaaiabdchaWjabcIcaOiabhMha5naaBaaabaGaemyAaKgabeaadaabbaqaaiabhIfaynaaBaaabaGaemyAaKMaemOAaOgabeaacqGGSaalaiaawEa7aiab=T7aSnaaBaaabaGaemOAaOgabeaacqGGSaalcqWIVlctcqGGPaqkcqWGWbaCcqGGOaakcqWHybawdaWgaaqaaiabdMgaPjabdQgaQbqabaWaaqqaaeaacqWF7oaBdaWgaaqaaiabdQgaQbqabaGaeiilaWIaeS47IWeacaGLhWoacqGGPaqkcqWGWbaCcqGGOaakcqWF7oaBdaWgaaqaaiabdQgaQbqabaGaeiykaKcabaGaemyAaKMaeyypa0JaeGymaedabaGaemOBa4gacqGHpis1aaaakiabgwSixNqbaoaalaaabaGaemyCaeNaeiikaGIaeCiwaG1aaSbaaeaacqWGPbqAcqWGQbGAaeqaamaaeeaabaGaeCyEaK3aaSbaaeaacqWGPbqAaeqaaiabcYcaSiabl+UimbGaay5bSdGaeiykaKIaemyCaeNaeiikaGIae83UdW2aaSbaaeaacqWGQbGAaeqaaiabcMcaPaqaaiabdghaXjabcIcaOiabhIfaynaaDaaabaGaemyAaKMaemOAaOgabaGaey4fIOcaamaaeeaabaGaeCyEaK3aaSbaaeaacqWGPbqAaeqaaiabcYcaSiabl+UimbGaay5bSdGaeiykaKIaemyCaeNaeiikaGIae83UdW2aa0baaeaacqWGQbGAaeaacqGHxiIkaaGaeiykaKcaaiabcYcaSaaa@C05F@

p( λ j MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaacceGae83UdW2aa0baaSqaaiabdQgaQbqaaiabgEHiQaaaaaa@3008@ ) and p(λ j ) is the prior probability of new and old position respectively, and they are cancelled out under uniform prior distribution; p ( X i j | λ j , ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiCaaNaeiikaGIaeCiwaG1aa0baaSqaaiabdMgaPjabdQgaQbqaaiabgEHiQaaakmaaeeaabaacceGae83UdW2aa0baaSqaaiabdQgaQbqaaiabgEHiQaaaaOGaay5bSdGaeiilaWIaeS47IWKaeiykaKcaaa@3CAA@ and p(X ij |λ j , ...) is the prior probability of QTL genotype conditional on new and old position, which has been described detailed previously; q ( X i j | y i , ) q ( X i j | y i , ) = p ( X i j | y i , ) p ( X i j | y i , ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqWGXbqCcqGGOaakcqWHybawdaWgaaqaaiabdMgaPjabdQgaQbqabaWaaqqaaeaacqWH5bqEdaWgaaqaaiabdMgaPbqabaGaeiilaWIaeS47IWeacaGLhWoacqGGPaqkaeaacqWGXbqCcqGGOaakcqWHybawdaqhaaqaaiabdMgaPjabdQgaQbqaaiabgEHiQaaadaabbaqaaiabhMha5naaBaaabaGaemyAaKgabeaacqGGSaalcqWIVlctaiaawEa7aiabcMcaPaaakiabg2da9KqbaoaalaaabaGaemiCaaNaeiikaGIaeCiwaG1aaSbaaeaacqWGPbqAcqWGQbGAaeqaamaaeeaabaGaeCyEaK3aaSbaaeaacqWGPbqAaeqaaiabcYcaSiabl+UimbGaay5bSdGaeiykaKcabaGaemiCaaNaeiikaGIaeCiwaG1aa0baaeaacqWGPbqAcqWGQbGAaeaacqGHxiIkaaWaaqqaaeaacqWH5bqEdaWgaaqaaiabdMgaPbqabaGaeiilaWIaeS47IWeacaGLhWoacqGGPaqkaaaaaa@6A3D@ and q ( λ j ) q ( λ j ) = k = 1 m p ( λ k j ) k = 1 m p ( λ k j ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaqcfa4aaSaaaeaacqWGXbqCcqGGOaakiiqacqWF7oaBdaWgaaqaaiabdQgaQbqabaGaeiykaKcabaGaemyCaeNaeiikaGIae83UdW2aa0baaeaacqWGQbGAaeaacqGHxiIkaaGaeiykaKcaaOGaeyypa0tcfa4aaSaaaeaadaqeWaqaaiabdchaWjabcIcaOiabeU7aSnaaBaaabaGaem4AaSMaemOAaOgabeaacqGGPaqkaeaacqWGRbWAcqGH9aqpcqaIXaqmaeaacqWGTbqBaiabg+GivdaabaWaaebmaeaacqWGWbaCcqGGOaakcqaH7oaBdaqhaaqaaiabdUgaRjabdQgaQbqaaiabgEHiQaaacqGGPaqkaeaacqWGRbWAcqGH9aqpcqaIXaqmaeaacqWGTbqBaiabg+Givdaaaaaa@591A@ , are all proposal ratio.

In step f, block sampling of the indicator variable matrix Φ j is expected to have a better performance than separately updating each γ kj in Φ j . Due to there are two possible values (0 or 1) for each model indicator γ kj , if m traits are investigated jointly, each model indicator matrix Φ j has 2mkinds of formations. The general formula of it can be written as:

W j , w 1 w 2 w m = [ γ 1 j = w 1 0 0 0 γ 2 j = w 2 0 0 0 0 γ m j = w m ] , MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaem4vaC1aaSbaaSqaaKqzaeGaemOAaOMaeiilaWIaem4DaCxcfa4aaSbaaWqaaKqzGdGaeGymaedameqaaKqzaeGaem4DaC3cdaWgaaadbaqcLboacqaIYaGmaWqabaqcLbqacqWIVlctcqWG3bWDjuaGdaWgaaadbaqcLboacqWGTbqBaWqabaaaleqaaOGaeyypa0ZaamWaaeaafaqabeabeaaaaaqaaiabeo7aNnaaBaaaleaacqaIXaqmcqWGQbGAaeqaaOGaeyypa0Jaem4DaC3aaSbaaSqaaiabigdaXaqabaaakeaacqaIWaamaeaacqWIVlctaeaacqaIWaamaeaacqaIWaamaeaacqaHZoWzdaWgaaWcbaGaeGOmaiJaemOAaOgabeaakiabg2da9iabdEha3naaBaaaleaacqaIYaGmaeqaaaGcbaGaeGimaadabaGaeGimaadabaGaeSO7I0eabaGaeSO7I0eabaGaeSy8I8eabaGaeSO7I0eabaGaeGimaadabaGaeGimaadabaGaeS47IWeabaGaeq4SdC2aaSbaaSqaaiabd2gaTjabdQgaQbqabaGccqGH9aqpcqWG3bWDdaWgaaWcbaGaemyBa0gabeaaaaaakiaawUfacaGLDbaacqGGSaalaaa@6EDD@

where, w k {0,1}, for k = 1, 2, ..., m. Because the prior probability of each γ kj is independent, the joint prior probability for all possible formations can be written as p ( Φ j = W l ) = k = 1 m p ( γ k j = w k ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiCaaNaeiikaGccceGae8NPdy0aaSbaaSqaaiabdQgaQbqabaGccqGH9aqpcqWGxbWvdaWgaaWcbaGaemiBaWgabeaakiabcMcaPiabg2da9maarahabaGaemiCaaNaeiikaGIaeq4SdC2aaSbaaSqaaiabdUgaRjabdQgaQbqabaGccqGH9aqpcqWG3bWDdaWgaaWcbaGaem4AaSgabeaakiabcMcaPaWcbaGaem4AaSMaeyypa0JaeGymaedabaGaemyBa0ganiabg+Givdaaaa@498B@ . Then the conditional posterior probability of Φ j can be written as

p ( Φ j = W j , w 1 w 2 w m | ) = p ( Φ j = W j , w 1 w 2 w m ) i = 1 n f ( y i | Φ j = W j , w 1 w 2 w m , ) g 1 { 0 , 1 } g 2 { 0 , 1 } g m { 0 , 1 } ( p ( Φ j = W j , g 1 g 2 g m ) i = 1 n f ( y i | Φ j = W j , g 1 g 2 g m , ) ) . MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xI8qiVKYPFjYdHaVhbbf9v8qqaqFr0xc9vqFj0dXdbba91qpepeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaemiCaaNaeiikaGccceGae8NPdy0aaSbaaSqaaiabdQgaQbqabaGccqGH9aqpcqWGxbWvdaWgaaWcbaqcLbqacqWGQbGAcqGGSaalcqWG3bWDjuaGdaWgaaadbaqcLboacqaIXaqmaWqabaqcLbqacqWG3bWDlmaaBaaameaajug4aiabikdaYaadbeaajugabiabl+UimjabdEha3LqbaoaaBaaameaajug4aiabd2gaTbadbeaaaSqabaGcdaabbaqaaiabl+UimbGaay5bSdGaeiykaKIaeyypa0tcfa4aaSaaaeaacqWGWbaCcqGGOaakcqWFMoGrdaWgaaqaaiabdQgaQbqabaGaeyypa0Jaem4vaC1aaSbaaeaacqWGQbGAcqGGSaalcqWG3bWDdaWgaaqaaiabigdaXaqabaGaem4DaC3aaSbaaeaacqaIYaGmaeqaaiabl+UimjabdEha3naaBaaabaGaemyBa0gabeaaaeqaaiabcMcaPmaaradabaGaemOzayMaeiikaGIaeCyEaK3aaSbaaeaacqWGPbqAaeqaamaaeeaabaGae8NPdy0aaSbaaeaacqWGQbGAaeqaaiabg2da9iabdEfaxnaaBaaabaGaemOAaOMaeiilaWIaem4DaC3aaSbaaeaacqaIXaqmaeqaaiabdEha3naaBaaabaGaeGOmaidabeaacqWIVlctcqWG3bWDdaWgaaqaaiabd2gaTbqabaaabeaacqGGSaalcqWIVlctaiaawEa7aiabcMcaPaqaaiabdMgaPjabg2da9iabigdaXaqaaiabd6gaUbGaey4dIunaaeaadaaeqaqaamaaqababaGaeS47IW0aaabeaeaadaqadaqaaiabdchaWjabcIcaOiab=z6agnaaBaaabaGaemOAaOgabeaacqGH9aqpcqWGxbWvdaWgaaqaaiabdQgaQjabcYcaSiabdEgaNnaaBaaabaGaeGymaedabeaacqWGNbWzdaWgaaqaaiabikdaYaqabaGaeS47IWKaem4zaC2aaSbaaeaacqWGTbqBaeqaaaqabaGaeiykaKYaaebmaeaacqWGMbGzcqGGOaakcqWH5bqEdaWgaaqaaiabdMgaPbqabaWaaqqaaeaacqWFMoGrdaWgaaqaaiabdQgaQbqabaGaeyypa0Jaem4vaC1aaSbaaeaacqWGQbGAcqGGSaalcqWGNbWzdaWgaaqaaiabigdaXaqabaGaem4zaC2aaSbaaeaacqaIYaGmaeqaaiabl+UimjabdEgaNnaaBaaabaGaemyBa0gabeaaaeqaaiabcYcaSiabl+UimbGaay5bSdGaeiykaKcabaGaemyAaKMaeyypa0JaeGymaedabaGaemOBa4gacqGHpis1aaGaayjkaiaawMcaaaqaaiabdEgaNnaaBaaabaGaemyBa0gabeaacqGHiiIZcqGG7bWEcqaIWaamcqGGSaalcqaIXaqmcqGG9bqFaeqacqGHris5aaqaaiabdEgaNnaaBaaabaGaeGOmaidabeaacqGHiiIZcqGG7bWEcqaIWaamcqGGSaalcqaIXaqmcqGG9bqFaeqacqGHris5aaqaaiabdEgaNnaaBaaabaGaeGymaedabeaacqGHiiIZcqGG7bWEcqaIWaamcqGGSaalcqaIXaqmcqGG9bqFaeqacqGHris5aaaacqGGUaGlaaa@E17D@

The approach to sample Φ j is similar to QTL genotypes sampling previously mentioned.

Post-MCMC analysis

For summarizing the posterior sample, we use the mean of the posterior sample to estimate the QTL effect and the residual (co)variance, and the mode of the posterior probability or the peak of the 2logeBF statistic to localize QTL. 2logeBF statistic was introduced by Yi et al.[17] into QTL mapping, and BF statistic is defined as the ratio of the posterior odds to the prior odds for inclusion against exclusion of the locus [24]. The critical value of BF is 3 or 2logeBF = 2.1 for declaring the existence of a QTL.

In single-trait analysis, we can pick the QTL by plotting the profile of the posterior probability or 2logeBF statistic against the genome. In multitrait analysis, if only two traits are considered jointly, we can use a three-dimension graph to summarize the statistic for all traits jointly (e.g., Figure 2 in [19]). However, if the number of trait is greater than 2, we can't plot them in one graph. Instead, we can solve the problem by plotting the marginal posterior probability distribution. If we divide the genome into H bins, and denote each bin of k th trait with ζ kg , for g = 1,2, ..., H, then the marginal posterior probability distribution of ζ kg is defined as p(ζ kg |y) = p[(ζ kg = λ kq ) ∩ (γ kq = 1)], where, q indicates the q th interval that locus ζ kg resides in. Then BF ( ζ k g ) = p ( ζ k g | y ) 1 p ( ζ k g | y ) 1 p ( ζ k g ) p ( ζ k g ) MathType@MTEF@5@5@+=feaagaart1ev2aaatCvAUfKttLearuWrP9MDH5MBPbIqV92AaeXatLxBI9gBaebbnrfifHhDYfgasaacPC6xNi=xH8viVGI8Gi=hEeeu0xXdbba9frFj0xb9qqpG0dXdb9aspeI8k8fiI+fsY=rqGqVepae9pg0db9vqaiVgFr0xfr=xfr=xc9adbaqaaeGaciGaaiaabeqaaeqabiWaaaGcbaGaeeOqaiKaeeOrayKaeiikaGIaeqOTdO3aaSbaaSqaaiabdUgaRjabdEgaNbqabaGccqGGPaqkcqGH9aqpjuaGdaWcaaqaaiabdchaWjabcIcaOiabeA7a6naaBaaabaGaem4AaSMaem4zaCgabeaadaabbaqaaiabhMha5bGaay5bSdGaeiykaKcabaGaeGymaeJaeyOeI0IaemiCaaNaeiikaGIaeqOTdO3aaSbaaeaacqWGRbWAcqWGNbWzaeqaamaaeeaabaGaeCyEaKhacaGLhWoacqGGPaqkaaGccqGHflY1juaGdaWcaaqaaiabigdaXiabgkHiTiabdchaWjabcIcaOiabeA7a6naaBaaabaGaem4AaSMaem4zaCgabeaacqGGPaqkaeaacqWGWbaCcqGGOaakcqaH2oGEdaWgaaqaaiabdUgaRjabdEgaNbqabaGaeiykaKcaaaaa@6180@ , which can be calculated at each possible locus for each trait, respectively.


  1. Xu CW, Li ZK, Xu S: Joint mapping of quantitative trait loci for multiple binary characters. Genetics. 2005, 169: 1045-1059. 10.1534/genetics.103.019406.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  2. Jiang C, Zeng ZB: Multiple trait analysis of genetic mapping for quantitative trait loci. Genetics. 1995, 140: 1111-1127.

    PubMed Central  CAS  PubMed  Google Scholar 

  3. Knott SA, Haley CS: Multitrait least squares for quantitative trait loci detection. Genetics. 2000, 156: 899-911.

    PubMed Central  CAS  PubMed  Google Scholar 

  4. Korol AB, Ronin YT, Itskovich AM, Peng J, Nevo E: Enhanced efficiency of quantitative trait loci mapping analysis based on multivariate complexs of quantitative traits. Genetics. 2001, 157: 1789-1803.

    PubMed Central  CAS  PubMed  Google Scholar 

  5. Mangin B, Thoquet P, Grimslev N: Pleiotropic QTL analysis. Biometrics. 1998, 54: 88-99. 10.2307/2533998.

    Article  Google Scholar 

  6. Eaves LJ, Neale MC, Maes H: Multivariate multipoint linkage analysis of quantitative trait loci. Behav Genet. 1996, 26: 519-525. 10.1007/BF02359757.

    Article  CAS  PubMed  Google Scholar 

  7. Liu JF, Liu YJ, Liu XG, Deng H-W: Bayesian mapping of quantitative trait loci for multiple complex traits with the use of variance components. Am J Hum Genet. 2007, 81: 304-320. 10.1086/519495.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  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: 805-816.

    PubMed Central  CAS  PubMed  Google Scholar 

  9. Yi N, Xu S: Bayesian mapping of quantitative trait loci for complex binary traits. Genetics. 2000, 155: 1391-1403.

    PubMed Central  CAS  PubMed  Google Scholar 

  10. Yi N, George V, Allison DB: Stochastic search variable selection for identifying multiple quantitative trait loci. Genetics. 2003, 164: 1129-1138.

    PubMed Central  CAS  PubMed  Google Scholar 

  11. Yi N, Xu S, Allison DB: Bayesian model choice and search strategies for mapping multiple epistatic quantitative trait loci. Genetics. 2003, 165: 867-883.

    PubMed Central  CAS  PubMed  Google Scholar 

  12. Yi N, Xu S, Allison DB: Bayesian model choice and search strategies for mapping interacting quantitative trait loci. Genetics. 2003, 165: 867-883.

    PubMed Central  CAS  PubMed  Google Scholar 

  13. Xu S: Derivation of the shrinkage estimates of quantitative trait locus effects. Genetics. 2007, 177: 1255-1258. 10.1534/genetics.107.077487.

    Article  PubMed Central  PubMed  Google Scholar 

  14. Wang H, Zhang YM, Li XM, Masinde GL, Mohan S, Baylink DJ, Xu S: Bayesian shrinkage estimation of quantitative trait loci parameters. Genetics. 2005, 170: 465-480. 10.1534/genetics.104.039354.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  15. 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: 1333-1344. 10.1534/genetics.104.040386.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  16. Yi N: A unified Markov chain Monte Carlo framework for mapping multiple quantitative trait loci. Genetics. 2004, 167: 967-975. 10.1534/genetics.104.026286.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  17. 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: 1865-1877. 10.1534/genetics.107.071365.

    Article  PubMed Central  PubMed  Google Scholar 

  18. Godsill SJ: On the relationship between MCMC model uncertainty methods. J Comput Graph Stat. 2001, 10: 230-248. 10.1198/10618600152627924.

    Article  Google Scholar 

  19. Gelman A, Carlin J, Stern H, Rubin D: Bayesian Data Analysis. 2004, London, Chapman & Hall

    Google Scholar 

  20. Hastings WK: Monte Carlo sampling methods using markov chains and their applications. Biometrika. 1970, 57: 97-109. 10.1093/biomet/57.1.97.

    Article  Google Scholar 

  21. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E: Equations of state calculations by fast computing machines. J Chem Phys. 1953, 21: 1087-1091. 10.1063/1.1699114.

    Article  CAS  Google Scholar 

  22. Tinker NA, Mather DE, Rossnagel BG, Kasha KJ, Kleinhofs A, Hayes PM, Falk DE, Ferguson T, Shugar LP, Legge WG, Irvine RB, Choo TM, Briggs KG, Ullrich SE, Franckowiak JD, Blake TK, Graf RJ, Dofing SM, Saghai Maroof MA, Scoles GJ, Hoffman D, Dahleen LS, Kilian A, Chen F, Biyashev RM, Kudrna DA, Steffenson BJ: Regions of the genome that affect agronomic performance in two-row barley. Crop Sci. 1996, 36: 1053-1062.

    Article  Google Scholar 

  23. Yi N, George V, Allison DB: Stochastic search variable selection for identifying multiple quantitative trait loci. Genetics. 2003, 164: 1129-1138.

    PubMed Central  CAS  PubMed  Google Scholar 

  24. Kass RE, Raftery AE: Bayes factors. J Am Stat Assoc. 1995, 90: 773-795. 10.2307/2291091.

    Article  Google Scholar 

Download references


We deeply thank four anonymous reviewers for their criticisms and comments which have greatly improved the presentation of the manuscript. This work was partly supported by Heilongjiang August First Land Reclamation University.

Author information

Authors and Affiliations


Corresponding author

Correspondence to Ming Fang.

Additional information

Authors' contributions

MF coordinated the study, developed the foundational principle of the method and wrote the computing program and the paper. Others were responsible for the simulation experiment, carried out the analysis of results and helped to consummate the whole paper.

Authors’ original submitted files for images

Rights and permissions

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Fang, M., Jiang, D., Pu, L.J. et al. Multitrait analysis of quantitative trait loci using Bayesian composite space approach. BMC Genet 9, 48 (2008).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: