Abstract
Background
In this article, I propose a modelselectionfree method to map multiple quantitative trait loci (QTL) in variance component model, which is useful in outbred populations. The new method can estimate the variance of zeroeffect QTL infinitely to zero, but nearly unbiased for nonzeroeffect QTL. It is analogous to Xu's Bayesian shrinkage estimation method, but his method is based on allelic substitution model, while the new method is based on the variance component models.
Results
Extensive simulation experiments were conducted to investigate the performance of the proposed method. The results showed that the proposed method was efficient in mapping multiple QTL simultaneously, and moreover it was more competitive than the reversible jump MCMC (RJMCMC) method and may even outperform it.
Conclusions
The newly developed Bayesian shrinkage method is very efficient and powerful for mapping multiple QTL in outbred populations.
Background
There are two kinds of models which can be used to map QTL in outbred populations, the allelic substitution model [13] and the variance component model [47]. In the allelic substitution model, the number of QTL alleles is assumed to be known, and the QTL allelic substitution is estimated by the given linkage phases of parents, which can be inferred from genotypes of family members. The least square [2] and maximum likelihood [1,3,8] of interval mapping are two popular statistical approaches for such models. Compared with the allelic substitute model, the variance component model is more robust because it can handle an arbitrary number of alleles with arbitrary modes of gene actions[9]. Moreover, the linkage phase of parents is unnecessary, which is nice since it is hard to accurately infer, particularly when family size is small, such as with human populations. Therefore, the variance component model is usually used to map QTL in outbred populations [47,912]. In the variance component model, the identitybaseddecent (IBD) matrix may be different for each locus and can provide information to localize the QTL. The least square method [10,11] and the maximum likelihood method [4,13] are also two important statistical methods for handling this model.
Because of the polygenic nature of quantitative traits, multiple QTL mapping is a problem of model selection. The least square method and the maximum likelihood method can nicely handle single QTL model, but is difficult for them to handle multiple QTL model. Recently, the Bayesian reversible jump MCMC (RJMCMC) method has been used to map multiple QTL in the variance component model [9,12]. However, it still has some disadvantages. Because the model dimension is variable, it usually has poor mixing character and is difficult to converge [1416]; moreover, it is also difficult to explore all the model space, especially in genomewide mapping where thousands of possible locus are scanned [16].
Therefore, in this article I proposed a modeldimensionfixed method, in which the estimate of variance is very precise for nonzeroeffect QTL, and gradually converges to zero for zeroeffect QTL. Therefore, special model selection is needless. It is similar to the recent Bayesian shrinkage estimation methods [15,1719], which are based on the allelic substitution model, whereas my method is based on the variance component model. The efficiency of the new method is demonstrated by a series of simulation experiments.
Method
Genetic model
Suppose that one has a sample of n individuals from outbred populations. Assuming that QTL dominant effect and polygenic dominant effect are absent. Then the linear model can be expressed as
where, y is the n × 1 phenotypic vector; β is the k × 1 vector of covariate effects; k is the number of the covariate; X is the n × k design matrix related to the covariate effects; a_{j }~ N(0, Φ_{j}) is the n × 1 vector of random QTL effect, for j = 1,2, ..., q, where Θ_{j }is the IBD matrix and can be inferred by the conditional expectation approach [20]; and is the QTL variance; e ~ N (0, I_{n}) is the vector of random error, where I_{n}, is the n × n identity matrix and is the residual variance; q is the maximum QTL number, which is set beforehand; g ~ N(0, A) is the n × 1 vector of random polygenic effect, here A is the additive relationship matrix and is the polygenic additive variance, the polygenic term g may be excluded from equation (1) in genomewide mapping. The variance component model can be expressed as
Similarly, the term of polygenic variance A should also be excluded from equation (2) in genomewide mapping.
Prior specification and joint posterior distribution
Yi and Xu [9] assigned a uniform prior distribution for QTL variance, , but in my method, the Jeffreys' hyper prior is assumed. The special prior is the key in the new method and will be illustrated in detail later. The prior for polygenic variance and residual variance is assumed to follow scaled inverted chisquare distribution with degree of freedom ω and scaled parameter s^{2 }(see also [21] for detail); and the prior for covariant effect and QTL position λ_{j }are assumed to follow normal distribution, β ~ N (β_{0}, V_{0}), and uniform distribution, respectively. The joint posterior distribution is given in Appendix.
Updating QTL variance by random walk MetropolisHastings algorithm
Because there is no close form for the posterior distribution of QTL variance , the MetropolisHastings algorithm [22,23] is used to simulate it. I firstly propose a new QTL variance and then accept it according to its acceptance probability.
Generating the new proposal QTL variance
I employ the Browne's method [24], a special random walk MetropolisHastings algorithm (RWMH) to update QTL variance. Firstly a new QTL variance is proposed and sampled from a scaled inverted chisquared distribution, conditional on the current value of QTL variance , , with the degree of freedom ν and the scaled parameter that equals the expectation of the current value , i.e. , and then the new QTL variance is accepted according to its accept probability. Since the new generated value closely relies on the old one, this approach is a special case of RWMH, and the degree of freedom ν is equivalent to the tuning parameter [21,24].
Calculating the acceptance probability
The new proposal QTL variance is accepted with probability equal to min (1, r), where,
and represents all elements of θ except . In equation (3), the first term is likelihood, the second is prior and the third is called proposal ratio or Hastings ratio [23]. Because the proposal distribution is not symmetric, , hr must be computed, and
MCMC implementations
The implementations of the MCMC algorithm are summarized as follows:
a. Initialize all parameters from legal values;
b. Update the covariate effect β;
c. Update the QTL variance ;
e. Update the polygenic variance ;
f. Update the residual variance ;
g. Update the QTL positions .
The covariate effect β is updated by the efficient Gibbs sampler; the updating approach of and are similar to that of (see also [21] for detail); the updating of the QTL position is illustrated in Appendix.
PostMCMC analysis
To summarize the posterior probability, I divide the genome into bins with interval of 1 cM and calculate weighted QTL variance for each bin. The weighted QTL variance is defined as the estimate of the QTL variance multiplied by its posterior probability at each bin (see also [21]), which is the modification of the weighted QTL effect [15]. If the profile of the weighted QTL variance generates a notable bump on the genome, the QTL is claimed as detected [15,19,25].
Results
I simulated 500 independent fullsib families with 6 individuals in each one, and therefore 3,000 individuals were investigated in my study. The parents of the fullsib families were randomly sampled from a large outbred populations in HardyWeinberg and linkage equilibrium. One chromosome with the length of 100 cM was simulated, and 11 evenly spaced markers covered the chromosomes with an average marker interval of 10 cM. I assigned 6 alleles for each marker and infinite alleles for each QTL. Three QTL were simulated on the genome positioned at 15 cM, 45 cM and 75 cM, respectively. The additive variances of the three QTL were respectively 0.5, 1.2 and 0.8, and the dominant variances were assumed absent. The residual effect for each offspring was randomly sampled from normal distribution with mean 0 and variance = 1.0. The polygenic variance = 1.0. The simulation method for the polygenic effect has been illustrated in [26]. The population mean equaled to 0. The phenotypic value for each sib was the sum of population mean, QTL effects, polygenic effect and residual effect. Therefore, the heritabilities explained by the three QTL were respectively 11.1%, 26.7% and 17.8%.
Before performing the simulated experiments, I firstly gave a default setup. The excepted QTL number q_{0 }= 2, which may lead to the maximum QTL number, (see [27] for detail); the degree of freedom for generating the proposal variance ν = ν_{A }= ν_{e }= 10; the hyperparameter ω_{A }= ω_{e }= 3 and (are approximately estimated as the phenotypic variance), and hence (equal to the expectation of their variance estimations, see [21]). Because ω_{A }and ω_{e }took a small value, the values of and would have ignorable effect on their estimates. The importance of the hyperparameter also has been illustrated in [21]. The MCMC ran for 21,000 rounds and the data was saved with every 10 rounds after the first 1,000 MCMC was discarded, so that there were 2,000 (20,000/10) posterior samples for posterior analysis.
Performance on simulated data with zeroQTL model
To demonstrate the special character of the proposed method, I firstly analyzed the data from the simulated zeroQTL model. The profiles of QTL intensity and weighted QTL variance are plotted in Figure 1a and Figure 1b, respectively. The profile of weighted QTL variance gives very noisy signals for QTL detection, but the values of weighted QTL variance are very tiny and the profile is much flat, which reflects that the proposed method can effectively shrink the values of the variance of zeroeffect QTL infinitely close to zero. Figure 2 gives the MCMC traces of the polygenic variance and the residual variance, and indicates that the estimates of them are all close to their true values, 1; moreover, the Markov chains of them converge fast and mix well.
Figure 1. Profiles of QTL intensity (a) and weighted QTL variance (b) obtained from the proposed method under simulated zeroQTL model.
Figure 2. Traces of polygenic variance and residual variance obtained from the proposed method.
Investigation into the performance of the special RWMH algorithm
I use a special RWMH algorithm to update the variance components, and the new proposal variance σ^{2 }(QTL variance, polygenic variance or residual variance) is sampled from the scaled inverted chisquared distribution with degree of freedom ν and scaled parameter the variance of the current round. In order to test the influence of ν, I set ν as 3, 15, 30, 50, 100, 150 and 200, respectively. The QTL intensity histogram [28] is plotted in Figure 3a. There are three peaks bumped on the chromosome, but QTL intensity is not used in QTL detection. I also plot the profile of weighted QTL variance, and the general pattern is given in Figure 3b. All other experiments have performed similar pattern, so the figures are not shown. I find the profile of weighted QTL variance is rather flat for the positions that have no QTL, which makes the signals of QTL clearer than QTL intensity. The parameter estimates are listed in Table 1, and there are no clear differences in parameter and standard deviation estimates for different ν. Furthermore, I summarized the acceptance rate of the MH sampler for the variance components. Because it is cumbersome to show them separately, I averaged the acceptance rate over all variance components under different setting of ν. I further plot the profile of the change of the acceptance rate against ν in Figure 4. It shows that the acceptance rate increases by ν, but the rate of change decease by ν. When ν is smaller than 30, the curve is much steeper, but it flatten when ν is larger than 30. The degree of freedom ν may influence the acceptance rate in the special RWMH algorithm, and hence it is equivalent to the tuning parameter in the traditional MetropolisHastings algorithm. Finally, I found that when ν is larger than 200, the shrinkage character is hardly held. The reasons will be addressed in Discussion.
Table 1. The estimates of the QTL parameters and their standard deviations obtained from the proposed method under different levels of degree of freedom.
Figure 3. Typical profiles of QTL intensity (a) and weighted QTL variance (b) from the proposed method.
Figure 4. The change of the acceptance rate against the degree of freedom by the proposed method.
Comparison with the regression method
I also used the software QTL Express [29], an IBD based regression method, to analyze the simulated data. The method is based on singleQTL model. The profile of F statistic is plotted in Figure 5. Only one QTL localized at 44 cM was detected, and the simulated three QTL were combined together. However, the three QTL can be separated successfully by the new method. The results clearly reflect the advantage of the new method that uses multiple QTL model over the regression method that bases on single QTL model.
Figure 5. Profiles of F statistic obtained from QTL Express.
Comparison with the RJMCMC method with repeat experiments
To compare the proposed method with the RJMCMC based method [9], I simulated 30 sets of data. For the RJMCMC method, the maximum QTL number and the expected QTL number were also the same as the default setup; the prior distribution of the QTL variance, polygenic variance and residual variance followed uniform distribution with endpoint being zero and phenotypic variance; the thinning interval was empirically set as 10; the burnin period was 1,000 and the length of the complete chain was 201,000, and hence, there were 10,000 samples saved for posterior analysis. It took ~ 5 hr for the new and RJMCMC method on a Pentium IV PC with a 2.60GHz processor and 1.00 GB RAM.
I list the empirical statistical power and the average estimates of 30 replications for both methods in Table 2, and the results show that: (1) the QTL detecting powers of the proposed method are slightly higher than that of the RJMCMC method; (2) there are no clear differences between the two methods in parameter estimates, and both are very precise.
Table 2. The average estimates of the QTL parameters from 30 replicated experiments in the proposed and RJMCMC method.
Application in genomewide mapping
In the genomewide mapping, I simulated a large genome of 2,000 cM, covered by 201 evenly spaced markers with interval 10 cM. Five QTL were simulated with positions and effects in Table 3. The total heritability was 69.4%, and the heritability explained by these QTL ranged from 8.1% to 22.4%. The residual variance was 1.5 and the polygenic variance was 0. The family structure and other parameters were the same as the previous simulation.
Table 3. QTL parameters and their estimates obtained from the two methods in the genomewide mapping.
Comparison with RJMCMC
The simulated data was analyzed with the proposed method and the RJMCMC method. In the proposed method, the excepted QTL number q_{0 }= 3, which results the maximum QTL number q = 8; the degree of freedom generating the proposal variance ν = ν_{A }= ν_{e }= 10; the hyperparameter ω_{A }= ω_{e }= 3 and for the polygenic variance and the residual variance. The MCMC ran for 51,000 rounds and the data was saved with every 10 rounds after the first 1,000 MCMC was discarded, so that there were 5,000 (50,000/10) samples for posterior analysis. In the RJMCMC method, also q_{0 }= 3 and q_{m }= 8, thinning interval 10 and burnin length 1,000, but the length of the complete chain was 201,000. In the genomewide mapping, for all the QTL that affect the trait are included in the model, the polygenic variance is excluded from the model and thus needn't be estimated.
The QTL intensity profiles of both methods are plotted in Figure 6a. It shows that the profile of the new method is higher than that of the RJMCMC method. But it is not sufficient to prove that the new method is more powerful than the RJMCMC method, because the QTL intensity is not used to detect QTL in shrinkage method. The profile of the weighted QTL variance is given in Figure 6b, and five clear bumps are found around their true simulated positions, which shows that the five simulated QTL are all detected by the new method. However, in the RJMCMC method, the estimated average number of QTL equaled to 3.37. The profile of the posterior QTL intensity is depicted in Figure 7, showing that the trait is mostly affected by three or four QTL with probability 0.565 or 0.359, and the estimated number of QTL is clearly smaller than the true number of QTL. The results suggest that my new method is competitive with the RJMCMC method, and may even outperform it. The computing time of the proposed method and the RJMCMC method were nearly equal and they took ~ 24 hr on a Pentium IV PC with a 2.60GHz processor and 1.00 GB RAM.
Figure 6. Profiles of QTL intensity (a) and weighted QTL variance (b) in the genomewide mapping. The true locations of the simulated QTL are indicated by upward arrows on the horizontal axis.
Figure 7. Posterior distribution of the number of QTL from the RJMCMC method.
Test on the sensitive of the maximum QTL number
The maximum QTL number q is a hyperparameter, which should be ascertained beforehand. In my study I followed the approach of [27] to ascertain it and . For testing the sensitive of the maximum q, I set the expected QTL number q_{0 }= 4 and 5, which led to q = 11 and 14, respectively. The profiles of weighted QTL variance are plotted in Figure 8, and show no clear differences; moreover, they are very similar to the profile in Figure 6b that uses q_{0 }= 3 and q = 8. The results demonstrated that the new method was not very sensitive to the value of the expected QTL number. I also ran the RJMCMC method under q_{0 }= 4 and 5, and the estimated QTL numbers were 3.42 and 3.46, respectively. Clearly, they were also smaller than the true values.
Figure 8. Profiles of weighted QTL variance for q = 11 (a) and q = 14 (b) in the genomewide mapping. The true locations of the simulated QTL are indicated by upward arrows on the horizontal axis.
Discussion
The Jeffreys' hyper prior for QTL variance, , which is much crucial in Bayesian shrinkage analysis for inbred line crosses [15,17,18], is also the key to the new method. Although the two methods handle different statistical models, the behavior is much similar, and they all need not special model selection. However, the vague prior may cause an improper posterior [30,31]. Ter Braak et al. [32] proposed to use the prior with a small value of δ but not the extreme value 0 to avoid generating improper posterior, while this extreme value was just used in [17] and my researches. I also attempted to set δ = 0.0001 and other values and did several times of experiments, and the results at δ = 0.0001 were essentially the same as those at δ = 0 (the results were not shown).
Because there are no close forms for the variance components, the MH algorithm is always used to update them. There are two kinds of MH algorithm, and one is the RWMH, in which the new proposal value is conditionally sampled on the old one; the other is called independent MetropolisHastings algorithm (IMH), and the proposed value is independently sampled on the old one. Generally, the RWMH is more efficient than the IMH, because in the RWMH the proposal values may automatically reach their main support region in the iterations. Another advantage of the RWMH is that the estimate of QTL variance for zeroeffect QTL may be gradually converged to zero. However, the special shrinkage character is hardly held by the IMH algorithm because it is usually low probability that the values of the proposal variance close to zero infinitely. Hence in this article I use a special RWMH to update QTL variance, which is also another key to my method.
The size of the tuning parameter may influence the efficiency of the RWMH algorithm. In this method, the degree of freedom ν is equivalent to the tuning parameter, if ν is smaller, the efficiency of the MH algorithm may decrease due to the low acceptance rate. But if ν is larger, although the acceptance rate increase, it is more difficult for the proposal variance components to explore their posterior distribution. If ν > 200, the chain will be stuck locally and the posterior distribution of variance components will be very difficult to explore, which explains why the shrinkage character is hardly hold by the proposed method when ν >200. In fact, the tuning parameter should be set appropriately, which makes the acceptance rate to be 10~ 40% [33]. Therefore, in my method, the optimal tuning parameter ν should range from 3 to 15 from Figure 4.
I assign a scaled inverted chisquare distribution for polygenic variance and residual variance, which makes the incorporation of the prior information possible, and this has been studied in previous work [21]. Certainly, other priors also can be used and then the formula of acceptance probability should be constructed appropriately.
The polygenic term is excluded in our genomewide mapping because many QTL with relative large effect are investigated in my simulated study. In practice, the trait may be affected by few QTL with substantial effects and many QTL with minor effects, and then it is necessary to include the polygenic effect in the model.
I proposed a basic method for mapping QTL in variance component models. The method is also important in fine mapping [3438] in which both linkage information and linkage disequilibrium (LD) information are utilized. If the markers are densely distributed, fine mapping provides an extremely powerful way for QTL mapping. The new method is also convenient to be modified to the simultaneous fine mapping of multiple QTL as long as the IBD matrix is appropriately constructed. Moreover, the method also can be extended to more complicated situation, such as that involving QTL dominant effect and epistatic effect.
I employed the approach of Yi et al. [27] to ascertain the maximum QTL number and found that the method was not very sensitive to the maximum QTL number. Theoretically, the maximum QTL number may be set as any value as long as it is greater than the actual QTL number. The simplest method is to assume that each marker interval contain one QTL, while it increases the computational burden. The appropriate selection of the maximum QTL number will contribute to saving the computational time.
The new developed method is very simple and easy to implement. The computer program is written in FORTRAN language, and it is also compiled into my software "BayesMapQTL.exe" which can be used to analyze simulated data, as well as field data with variable family size. Both program and software are available for request.
Conclusion
In this research, I proposed a Bayesian shrinkage estimation method for mapping multiple QTL. Different from Xu's (2003) shrinkage method that discriminately estimates QTL substitute effect, my method can shrinkage estimate QTL variance so that it need no special model selection. Simulation experiments show that the proposed method is efficient in simultaneously mapping multiple QTL in outbred populations and may even outperform the RJMCMC method.
Appendix
Joint posterior distribution
Observable parameters include marker information M, additive genetic relationship matrix A, covariate matrix X and phenotypic values ; unobservables include QTL positions λ = , model effects and hyperparameters . Joint posterior probability of unobservables then can be expressed as
where,
and
Update QTL position
The proposal position is moved around the old one, = λ_{j }+ d, where d is a random number sampled from uniform distribution with bound kcM and kcM, where k is a predetermined tuning parameter and equals to 1 for chromosome segment analysis, and 20 for genomewide scan in my study. When the new position is proposed, the IBD matrix is calculated according to marker information, and then the new position is accepted with probability equal to min(1, r),
The proposal ratio or Hastings ratio hr_{p }= 1 due to the symmetric uniform proposal.
Acknowledgements
I am grateful to two anonymous reviewers for their helpful comments which have greatly improved the presentation of the manuscript. Lijun Pu, Weixuan Fu and one of reviewers helped proofread the language. I also thank BMC editorial for the waiver of the publication cost. This research was supported by Chinese Heilongjiang Education Ministry grant 11541254 to M. F.
References

Knott SA, Haley CS: Maximum likelihood mapping of quantitative trait loci using fullsib families.
Genetics 1992, 132(4):12111222. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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

Knott SA, Haley CS: Methods for multiplemarker mapping of quantitative trait loci in halfsib populations.
Theor Appl Genet 1996, 93:7180. Publisher Full Text

Goldgar DE: Multipoint analysis of human quantitative genetic variation.
Am J Hum Genet 1990, 47(6):957967. PubMed Abstract  PubMed Central Full Text

Amos CI: Robust variancecomponents approach for assessing genetic linkage in pedigrees.
Am J Hum Genet 1994, 54(3):535543. PubMed Abstract  PubMed Central Full Text

Xu S, Atchley WR: A random model approach to interval mapping of quantitative trait loci.
Genetics 1995, 141(3):11891197. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Almasy L, Blangero J: Multipoint quantitativetrait linkage analysis in general pedigrees.
Am J Hum Genet 1998, 62(5):11981211. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Liu Y, Jansen GB, Lin CY: Quantitative trait loci mapping for dairy cattle production traits using a maximum likelihood method.
J Dairy Sci 2004, 87(2):491500. PubMed Abstract  Publisher Full Text

Yi N, Xu S: Bayesian mapping of quantitative trait loci under the identitybydescentbased variance component model.
Genetics 2000, 156(1):411422. PubMed Abstract  PubMed Central Full Text

Haseman JK, Elston RC: The investigation of linkage between a quantitative trait and a marker locus.
Behav Genet 1972, 2(1):319. PubMed Abstract  Publisher Full Text

Fulker DW, Cardon LR: A sibpair approach to interval mapping of quantitative trait loci.
Am J Hum Genet 1994, 54(6):10921103. PubMed Abstract  PubMed Central Full Text

Liu J, Liu Y, Liu X, Deng HW: Bayesian mapping of quantitative trait loci for multiple complex traits with the use of variance components.
Am J Hum Genet 2007, 81(2):304320. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Schork NJ: Extended multipoint identitybydescent analysis of human quantitative traits: efficiency, power, and modeling considerations.
Am J Hum Genet 1993, 53(6):13061319. PubMed Abstract  PubMed Central Full Text

Yi N, George V, Allison DB: Stochastic search variable selection for identifying multiple quantitative trait loci.
Genetics 2003, 164(3):11291138. PubMed Abstract  PubMed Central Full Text

Wang H, Zhang YM, Li X, Masinde GL, Mohan S, Baylink DJ, Xu S: Bayesian shrinkage estimation of quantitative trait loci parameters.
Genetics 2005, 170(1):465480. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

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

Xu S: Estimating polygenic effects using markers of the entire genome.
Genetics 2003, 163(2):789801. PubMed Abstract  PubMed Central Full Text

Xu S: Derivation of the shrinkage estimates of quantitative trait locus effects.
Genetics 2007, 177(2):12551258. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Yang R, Xu S: Bayesian shrinkage analysis of quantitative trait Loci for dynamic traits.
Genetics 2007, 176(2):11691185. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Xu S, Gessler DD: Multipoint genetic mapping of quantitative trait loci using a variable number of sibs per family.
Genet Res 1998, 71(1):7383. PubMed Abstract  Publisher Full Text

Fang M, Liu S, Jiang D: Bayesian composite model space approach for mapping quantitative trait Loci in variance component model.
Behav Genet 2009, 39(3):337346. PubMed Abstract  Publisher Full Text

Metropolis N, A RW, Rosenbluth MN, Teller AH, Teller E: Equation of state calculations by fast computing machines.
J Chem Phys 1953, 21:10871091. Publisher Full Text

Hastings WK: Monte Carlo sampling methods using Markov chains and their applications.
Biometrika 1970, 57:97109. Publisher Full Text

Browne WJ: Applying MCMC Methods to Multilevel Models. Bath: University of Bath; 1998.

Xu C, Wang X, Li Z, Xu S: Mapping QTL for multiple traits using Bayesian statistics.
Genet Res 2009, 91(1):2337. PubMed Abstract  Publisher Full Text

Gessler DD, Xu S: Using the expectation or the distribution of the identity by descent for mapping quantitative trait loci under the random model.
Am J Hum Genet 1996, 59(6):13821390. PubMed Abstract  PubMed Central Full Text

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

Sillanpaa MJ, Arjas E: Bayesian mapping of multiple quantitative trait loci from incomplete inbred line cross data.
Genetics 1998, 148(3):13731388. PubMed Abstract  PubMed Central Full Text

Seaton G, Haley CS, Knott SA, Kearsey M, Visscher PM: QTL Express: mapping quantitative trait loci in simple and complex pedigrees.
Bioinformatics 2002, 18(2):339340. PubMed Abstract  Publisher Full Text

Hobert JP, Casella G: The effect of improper priors on Gibbs sampling in hierarchical linear mixed models.

Gelman A, Carlin J, Stern H, Rubin D: Bayesian data analysis. London: Chapman and all/CRC; 2004.

ter Braak CJ, Boer MP, Bink MC: Extending Xu's Bayesian model for estimating polygenic effects using markers of the entire genome.
Genetics 2005, 170(3):14351438. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Roberts GO, Rosenthal JS: Optimal scaling for various MetropolisHastings algorithms.
Statist Sci 2001, 16(4):351367. Publisher Full Text

Meuwissen TH, Goddard ME: Fine mapping of quantitative trait loci using linkage disequilibria with closely linked marker loci.
Genetics 2000, 155(1):421430. PubMed Abstract  PubMed Central Full Text

Meuwissen TH, Karlsen A, Lien S, Olsaker I, Goddard ME: Fine mapping of a quantitative trait locus for twinning rate using combined linkage and linkage disequilibrium mapping.
Genetics 2002, 161(1):373379. PubMed Abstract  PubMed Central Full Text

Meuwissen TH, Goddard ME: Prediction of identity by descent probabilities from markerhaplotypes.
Genet Sel Evol 2001, 33(6):605634. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Meuwissen TH, Goddard ME: Mapping multiple QTL using linkage disequilibrium and linkage analysis information and multitrait data.
Genet Sel Evol 2004, 36(3):261279. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Lee SH, Werf JH: Simultaneous fine mapping of multiple closely linked quantitative trait Loci using combined linkage disequilibrium and linkage with a general pedigree.
Genetics 2006, 173(4):23292337. PubMed Abstract  Publisher Full Text  PubMed Central Full Text