Bayes factors for detection of Quantitative Trait Loci

A fundamental issue in quantitative trait locus (QTL) mapping is to determine the plausibility of the presence of a QTL at a given genome location. Bayesian analysis offers an attractive way of testing alternative models (here, QTL vs. no-QTL) via the Bayes factor. There have been several numerical approaches to computing the Bayes factor, mostly based on Markov Chain Monte Carlo (MCMC), but these strategies are subject to numerical or stability problems. We propose a simple and stable approach to calculating the Bayes factor between nested models. The procedure is based on a reparameterization of a variance component model in terms of intra-class correlation. The Bayes factor can then be easily calculated from the output of a MCMC scheme by averaging conditional densities at the null intra-class correlation. We studied the performance of the method using simulation. We applied this approach to QTL analysis in an outbred population. We also compared it with the Likelihood Ratio Test and we analyzed its stability. Simulation results were very similar to the simulated parameters. The posterior probability of the QTL model increases as the QTL effect does. The location of the QTL was also correctly obtained. The use of meta-analysis is suggested from the properties of the Bayes factor.


INTRODUCTION
Mapping of quantitative trait loci (QTLs) is a rapidly evolving topic in Statistical Genomics. Several procedures have been described for mapping QTLs in experimental crosses [10,20,21] and in outbred populations [1,14,33]. In all these settings, hypothesis testing is one of the most delicate and controversial issues.
From a Bayesian perspective, a procedure was described by Hoeschele and van Raden [16,17]. It allows the estimation of QTL effects, and it has been implemented using Monte Carlo methods in crosses [27,29] and in outbred populations [18,28]. In a Bayesian setting, QTL detection involves the calculation of the Bayes factor (BF) or the posterior probability of the models [19,22]. The Bayes factor provides a rigorous framework for model testing in terms of probability, and it does not require assuming any asymptotic property as it does for the Likelihood Ratio Test (LRT). Unfortunately, the exact calculation of general BF is not feasible for relatively complex models [19]. For this reason, Monte Carlo methods, such as the Harmonic Mean Estimation [24] or the Monte Carlo marginal likelihood [3], have been developed, as reviewed by Gelman and Meng [7] and Han and Carlin [11]. Moreover, some other alternatives for providing posterior probabilities have been suggested [4,8]. Among these methods, the Reversible Jump Markov Chain Monte Carlo [8] has been used in the scope of QTL detection [13,18,28,30,32]. This method provides a useful tool for calculating the posterior probability of each model, although it becomes more difficult as the complexity of the models increases (multiple markers or multiple alleles at the QTL).
Following the point null Bayes factor approach [2], García-Cortés et al. [6] described a procedure to compare nested variance component models from the perspective of a Dirac Delta approach. The objective of the present paper is to describe a point null approach to calculate the Bayes factor using a Markov Chain Monte Carlo method. The method was compared with LRT and its performance and stability in QTL mapping.

Theory
We compare models that only differ by the presence of a QTL. These are considered as nested models because the parameters of the simple model (ω) are a subset of the parameters of the complex model (θ, ω). Following the procedure described in the Appendix, if we compare two nested models, one complete (A), and one reduced (B), BF can be calculated from the following simple expression: where p A (θ = 0) and p A (θ = 0|y) are the prior and posterior densities of θ. First, we will apply this procedure to a simple QTL model, and, later on, we will analyze a mixed QTL model which also includes polygenic effects.

Calculation of Bayes factor
Now, we present the Bayes factor for a model containing a QTL effect over a no-QTL model. Consider the following model (model 1): where y contains the phenotypic records, µ is the overall mean, Z is the incidence matrix relating observations to QTL effects (q) and e is the vector of residuals, q and e are assumed to be normally distributed: with σ 2 q being the variance explained by the QTL, σ 2 e , the residual variance, and Q, the relationship matrix between QTL effects. Model 1 can be reparameterized as: Consequently, where h 2 q = σ 2 q /σ 2 p is the proportion of phenotypic variation explained by the QTL, and σ 2 p = σ 2 q + σ 2 e is the phenotypic variance. The joint distribution of all variables in model 1 is: where: and 0 otherwise, where k 1 and k 2 are two small enough values to ensure a flat distribution over the parametric space. The null hypothesis model is the no-QTL model (model 2): e ∼ N(0, Iσ 2 p ). Then, the joint distribution of records and parameters is: where we can assume that prior distributions p 2 (µ) and p 2 (σ 2 p ) are identical to equations (2) and (3), respectively, and (1): because p 1 (h 2 q = 0) = 1.

Mixed QTL model
Let us now consider a mixed inheritance model (model 3) that includes polygenic effects (u): where u ∼ N(0, Aσ 2 u ), A being the polygenic relationship matrix and σ 2 u the polygenic genetic variance, Z 1 and Z 2 are incidence matrices. Notation and distribution of random QTL effects (q) and residuals (e) are assumed to be the same as in model 1.
This model can again be reparameterized as: consequently, where h 2 u = σ 2 u /σ 2 p is the proportion of phenotypic variation explained by polygenes and σ 2 p is the phenotypic variance σ 2 u + σ 2 q + σ 2 e . Records and parameters are jointly distributed as: and 0 otherwise, 1 k 2 and otherwise. (6) Note that, assuming prior independence, marginal priors of h 2 q and h 2 u are: (1,2). Model 3 will be compared to the following null hypothesis model (model 4): where: e * = Z 1 u + e, consequently e * ∼ N(0, V) where priors for µ and σ 2 p are the same as in model 3, equations (5) and (6), respectively. Prior distribution for h 2 u is p 4 h 2 u = U (0, 1) = p 3 h 2 u |h 2 q = 0 . U denotes a uniform distribution. As before, model 4 is a particular case of model 3 when h 2 q = 0. The BF of model 3 versus model 4:

Simple QTL model a) Simulation
A two-generation pedigree was simulated, 15 sires were mated to 5 dams each, with 5 offspring per dam. Four different cases were simulated as described in Table I, with different heritabilities and locations of the QTL. A single chromosome of 60 cM in length was simulated with four completely informative markers located at 0, 20, 40 and 60 cM. Phenotypes and marker genotypes were assumed to be known in all animals. Simulation of phenotypic records was performed by an overall mean (µ), a random QTL effect (q) and a residual (e). Twenty replicates were run per case, except in case II, where 1 000 replicates were run to compare BF with the Likelihood Ratio Test (LRT).

b) Calculation of the Marker Relationship Matrix (Q)
The (co)variance matrix (Q) at the candidate QTL position was obtained as the probabilities for individuals of sharing alleles identical by descent [23]. The genetic origin of marker alleles was unambiguously known. In this case, the probability of identity by descent was easy to calculate by comparing the haplotypes of the flanking markers between both half-and full-sibs. In these cases, the relationship matrix between sibs (i and j) at position x can be calculated from: is the probability for chromosomes H i and H j of sharing a replicate of the allele at position x. Several cases can be considered in relation to the structure of markers between parents and offspring, where λ is the genetic distance between markers. Probabilities of identity by descent at position x are: 1. Both haplotypes present the same alleles at the flanking markers and in the same phase as their parents where r x , r λ−x , r λ are the recombination fraction between the right marker and position x, between the x and the left marker and between both markers, respectively. 2. Both haplotypes share both markers but in a different phase to their parents Both haplotypes do not share any markers and the haplotypes are in the same phase as their parents 4. Both haplotypes do not share any markers but they are in a different phase to their parents Both haplotypes only share the right marker Both haplotypes only share the left marker The coefficient of relationship between parents and progeny is always 0.5. Relationship matrices in cases involving more complicated pedigrees or noninformative markers can be calculated after an explicit analysis [15,31] or numerically by using MCMC [9,25].

c) Calculation of the Bayes factor
Density p 1 ( h 2 q = 0 y) suffices to obtain BF (equation (4)). This value can be obtained from the Gibbs sampler output by averaging the full conditional densities of each cycle at h 2 q = 0 using the Rao-Blackwell argument. The Gibbs sampler algorithm involves updating samples from the full conditional distributions, which are: where n is the number of records. Note that h 2 q is involved in the structure of V, and this is not a standard probability distribution. Thus, a Metropolis-Hastings step [12] within each Gibbs sampling cycle was performed. The length of the Gibbs sampler was 10 000 cycles after discarding the first 1 000 iterations. A genomic scan was performed, in which, BF was computed every cM.

d) Meta-analysis
From the definition of BF where PO is the Posterior odds between models and PrO is the Prior odds. Let us consider the successive simulated replicates (n different data sets) as a sequential number of experiments. Then, the joint posterior odds is where BF i is the Bayes factor calculated from the ith replicate.

e) Likelihood Ratio Test
In case II of simulation (10% of phenotypic variation explained by a QTL), 1 000 replicates were simulated. In every replicate, BF and LRT were calculated. LRT was computed according to the following expression: where L 1 μ,ĥ 2 q , σ 2 p is the likelihood under the model 1 at maximum likelihood estimates μ,ĥ 2 q , σ 2 p and L 2 μ,σ 2 p is the likelihood under the model 2 at maximum likelihood estimates under this model. Maximum likelihood estimates were obtained through a simplex algorithm [26].
Twice the logarithm of the Likelihood Ratio Test (LLRT) was calculated to compare with limits of significance with a chi square distribution of 1 and 2 degrees of freedom as suggested by Grignola et al. (1996). Later on, LLRT was compared to the logarithm of the Bayes factor (LBF).

Mixed QTL model
The population structure was as in the previous model with the simulation parameters given in Table I. The simulation model included a random polygenic effect, and in all cases σ 2 q + σ 2 u = 0.5σ 2 p . Bayes factors were calculated at positions of 10, 30 and 50 cM. The Bayes factor was computed from the output of a Gibbs Sampler using the argument of Rao-Blackwell, as before. The calculation of Q matrix was performed as in the previous chapter. The numerator relationship matrix (A) between polygenic effects was calculated from the pedigree information [23].
Conditional distributions involved are the same as in model 1, except that here , and the conditional sampling for h 2 u requires an extra Metropolis-Hastings step at every iteration. Twenty replicates were performed for each of the four different cases of simulation.

Stability Analysis
Two replicates of case II (10% of variation was located on the QTL) were analyzed 1 000 times with Monte Carlo chains of 20, 100, 500, 2 500 and 10 000 iterations. Means and variances of BF and posterior probability were calculated for every case.

Simple QTL model
The results of the single QTL model are presented in Table II for the four different cases of simulation. Following Kass and Raftery [19], values of the Bayes factors were classified into five categories according to posterior probability: a) smaller than 0.5 (BF < 1), b) between 0.5 and 0.762 (1 < BF < 3.2), c) between 0.762 and 0.909 (3.2 < BF < 10), d) between 0.909 and 0.990 (10 < BF < 100), and e) greater than 0.990 (BF > 100). The posterior probability of the presence of a QTL depended on its effect rather than on its relative position on the chromosome, because the simulation assumed equallyinformative and spaced markers. In case I (h 2 q = 0), the no-QTL model had a higher probability than the QTL model in all replicates, and the percentage of replicates, when the QTL model was more likely, increased with the effect of the QTL (cases II, III and IV).
In the context of the simulation study, the properties of posterior estimates by repeated sampling are also presented in Table II. It is interesting to note that both the average of posterior mean estimates of h 2 q and the position were close to the simulated values, especially as the QTL effect increased. The posterior mean estimates of h 2 q were biased upwards when the QTL effects were small, because of the effect of the lower bound of the parametric space. The average position at the maximum Bayes factor was close to the simulated value, and the average posterior probability of the QTL model increased to 0.96 in cases III and IV (h 2 q = 0.20). Meta-analysis results from the joint analysis of the 20 replicates are presented in Figures 1 to 4. Conclusive evidence for a QTL together with an accurate estimation of its location were observed in cases II, III and IV. In case I, when the no-QTL effect was simulated, the maximum PO was 2 × 10 −25 , and the no-QTL model was far more likely than the QTL model.
Finally, we compared the log-likelihood criteria (LLRT) with the logarithm of BF (LBF) in 1 000 replicates of case II (h 2 q = 0.10). As can be observed in Figure 5, both criteria were strongly related. In replicates, the correlation coefficient between these two criteria was higher than 0.99. An LLRT greater than 5.99 is exhibited by 62.1% of replicates which represented the 5% of the first type error, when chi-square with 2 degrees of freedom was assumed. Moreover, 78.4% of replicates presented an LLRT greater than 3.84, corresponding to the   Table III shows the results obtained for cases without QTLs. In these cases, the most probable model was the "no-QTL" model in almost all replicates. Nevertheless, in 3 out of 60 replicates, the model including QTL effects had larger posterior probabilities than the "no-QTL" model. The presence of polygenic genetic variance may lead to wrong estimates of the QTL effect, because of similarity of relationship matrices. It can also be observed that a spurious estimate of σ 2 q appeared when the mixed inheritance model was used. As in likelihood procedures, variances in Bayesian methods were constrained within the positive values, but σ 2 q + σ 2 u was accurately estimated. The most sensible action is to test whether the probability of presence of a QTL is small enough to justify the use of the simple infinitesimal model. The Meta-analysis shows that evidence against a QTL is conclusive along the chromosome.

Mixed QTL model
Consider the second case of simulation (Tab. IV). It can be observed that the probability of the presence of a QTL was smaller than in the equivalent simulation case when σ 2 u = 0. The power of the analysis decreased because of the complexity of the model, with the presence of two genetics effects (QTL and polygenic). However, when all replicates were analyzed jointly via the meta-analysis, the posterior probability of QTL presence is almost 1. As in Table III, the posterior mean estimates of σ 2 q were still greater than the simulated values, but the difference between simulated and estimated values was smaller.
Results of the third and fourth cases of simulation are presented in Tables V and VI, respectively. In these cases, the probability of the presence of a QTL was greater than 0.5 at the true position of the QTL, and the probability decreased as the distance between the true position of the QTL and the position being analyzed increased. If the replicate estimates were pooled in a meta-analysis, the position of the QTL was estimated accurately, although the posterior mean estimates were still greater than the corresponding simulated values. If the QTL was located in a centromeric position, then any scanned position along the chromosome suggested its presence (Tab. V). In contrast, if the QTL was located in a telomeric position, then distant positions did not support the existence of a QTL (Tab. VI). Finally, a stability analysis was performed with two replicates of case II with the mixed QTL model. As can be observed in Table VII, the Monte Carlo approach described here is stable and accurate to estimate the Bayes factor or posterior probability, when the number of iterations is moderately large (2 500 or greater). Estimates of Bayes factor are unbiased, even when a small number of iterations are considered. Posterior probabilities are slightly biased with a small number of iterations, because of the range limits between 0 and 1. In the present study, all replicates were analyzed with 10 000 iterations after discarding the first 1 000. Thus we can conclude that the Bayes factor or posterior probabilities are accurately estimated.

DISCUSSION
We have developed a stable procedure to calculate the Bayes factor in a QTL analysis framework. The percentage of replicates that assigns strong evidence of QTL presence increases with the QTL effect. BF also allows to determine the position of the QTL.
Equation (1) avoids the instability of other MCMC approaches to obtaining the BF. The BF estimate from (1) is stable and can be computed with a relatively short chain, as shown in Table VII. The results are consistent with the rapid mixing of the variables observed by García-Cortés et al. [6], after integrating out the random effects. This fact represents a great advantage over other numerical approximations to the Bayes factor or posterior probabilities such as the harmonic mean [24] or the Reversible Jump Markov Chain Monte Carlo [8]. However, the procedure is not general in the sense that it can be used only in the context of nested models. This is not a serious disadvantage in QTL studies, where the interest is usually centered in ascertaining whether a QTL is segregating at a given position. Comparison of nested models (QTL vs. no-QTL models) is required. This approach cannot be generally applied to other situations, i.e., testing a non-linear model vs. a linear model. In relation to other procedures, such as the Likelihood Ratio Test [9,34], the Bayes factor provides a rigorous and clear framework to compare competing models. Its results can be expressed in terms of probability. It means that the calculation of significance levels either with simulation [5] or with theoretical approximations are unnecessary. In the scope of the simulation study, the correlation between both criteria was very high, and the power of the test was similar to LRT, when a 5% type I error was considered. However, the Bayes factor does not depend on the asymptotic properties and it can be used safely even with small samples. The classical hypothesis tests try to discard the null hypothesis in favour of an alternative hypothesis, while the Bayes factor provides a ratio of probabilities between models, without any requirement to define the null or the alternative model.
Another important property when using meta-analysis with different sources of information is to calculate the overall posterior odds by multiplication of BFs from different experiments, in contrast with alternative procedures, in which the only way to combine information is to jointly analyze all data. A strong concordance between simulation and results from meta-analysis was observed. It must be noticed that each meta-analysis was carried out using 300 sire families and a total of 9 300 records. However, it must be taken into account that meta-analysis can only be carried out when the conditions for trait measurements in all the experiments are similar.
Certain aspects need to be highlighted in relation to the use of the Bayes factor. First, the Bayes factor strongly depends on the prior distributions assumed for all the parameters in the model. For that reason, some caution must be practiced and a sensitivity analysis is fully recommended. In this study, we considered flat priors for h 2 q for the simple QTL model and a flat prior for the pair (h 2 q , h 2 u ) in the mixed QTL model. However, the procedure can be applied to any other prior distribution on intraclass correlations. It must also be highlighted that for simplicity purposes it is necessary to assume independent prior distributions for heritabilities and phenotypic variances in calculating the Bayes factor.
In this study, we compared the model with and without a QTL at a given location. If we are interested in testing the QTL at any position along the chromosome, the following approach can be considered. Let BF = p( QTL| l) p (no-QTL) be the BF of presence of a QTL at a given location l, then the BF c of the presence of a QTL at any position of the chromosome over the non-presence of a QTL is obtained by computing the following integral along its parametric space (Ω l ): over any predetermined prior distribution of location of the QTL (p(l)), such as uniform distribution along the chromosome or any other distribution defined by the density of candidate genes or other criteria. An alternative approach is to include the location of the QTL in the model as an extra variable, and the marginal posterior distribution of the location will also be obtained. This approach is equivalent to calculating the above integral over marginal distribution. In practice, more or less dense genotyping along the genome is available, and the question arise whether a given chromosome contains QTLs above a prefixed effect. In this case a series of BFs can be formulated, i.e., a model in which a set of chromosomes contains QTLs vs. a model in which only a subset of these chromosomes contains QTLs. This does not require any novel theoretical developments.
In conclusion, the proposed method is able to split σ 2 p into σ 2 q and σ 2 e and correctly identifies whether a particular location substantially contributes to covariance between individuals. The ability to detect QTLs in individual experiments is relatively low, thus meta-analysis will be necessary for practical purposes. The proposed procedure allows us to easily combine information from different experiments.