Bayesian Mapping of Multiple Quantitative Trait Loci From Incomplete Inbred Line Cross Data

Corresponding author: Mikko J. Sillanpää, Rolf Nevanlinna Institute, Research Institute of Mathematics, Statistics, and Computer Science, P.O. Box 4, FIN-00014 University of Helsinki, Finland, mjs{at}rolf.helsinki.fi (E-mail).

Communicating editor: Z-B. ZENG

	ABSTRACT

TOP ABSTRACT MODEL SIMULATION ANALYSIS DISCUSSION APPENDIX A: ESTIMATION OF... APPENDIX B LITERATURE CITED

A novel fine structure mapping method for quantitative traits is presented. It is based on Bayesian modeling and inference, treating the number of quantitative trait loci (QTLs) as an unobserved random variable and using ideas similar to composite interval mapping to account for the effects of QTLs in other chromosomes. The method is introduced for inbred lines and it can be applied also in situations involving frequent missing genotypes. We propose that two new probabilistic measures be used to summarize the results from the statistical analysis: (1) the (posterior) QTL-intensity, for estimating the number of QTLs in a chromosome and for localizing them into some particular chromosomal regions, and (2) the locationwise (posterior) distributions of the phenotypic effects of the QTLs. Both these measures will be viewed as functions of the putative QTL locus, over the marker range in the linkage group. The method is tested and compared with standard interval and composite interval mapping techniques by using simulated backcross progeny data. It is implemented as a software package. Its initial version is freely available for research purposes under the name Multimapper at URL http://www.rni.helsinki.fi/~mjs.

WHEN two purely homozygous, inbred, very genetically divergent lines are crossed, all offspring (F₁ generation) are genetically identical, being heterozygous at each locus. In the haplotypes of the F₁ individuals, the locus next to each quantitative trait locus (QTL) has the same allele as it had in the parental haplotype. This is because, in this ideal case, parents are homozygous at each locus and recombination events cannot change haplotypic arrangements. Therefore, linkage disequilibrium (nonrandom allelic association) in this group is maximal.

When an F₂ or backcross generation is produced, linkage disequilibrium is reduced slightly but still remains at a high level. The degree of reduction depends on the distance and on the recombination fraction between the considered QTL and the nearby marker locus. If mating is continued till F₃ and the succeeding generations, disequilibrium area surrounding a QTL is reduced further in each generation. This is why the backcross or F₂ intercross data from inbred lines is particularly suitable for QTL mapping.

The commonly used QTL mapping methods for plants and animals introduced recently use offspring data from divergent inbred lines in backcross or F₂ intercross design. The reason for using such a design is to maximize linkage disequilibrium and the amount of heterozygosity (information content) in meioses. Examples of such methods are interval mapping (LANDER and BOTSTEIN 1989 ), least squares method (HALEY and KNOTT 1992 ), and composite interval mapping (JANSEN 1993 ; JANSEN and STAM 1994 ; ZENG 1993 , ZENG 1994 ). For some species, it is very impractical, time consuming and also expensive to produce inbred lines. In such cases, methods have been developed for considering crosses between outbred lines that are genetically divergent and show two very separate groups of phenotypic values, for example, due to different selection histories. One such procedure was presented in HALEY et al. 1994 , where the analysis was done in terms of line origins, with the assumption that crossed lines are "fixed" for different genes (or alleles) and would then show homozygosity in most of the QTL loci. This method (design) requires genotypic data from the parental and grandparental generations in addition to genotypic and phenotypic offspring data. Recently JANSEN 1996 introduced a general method for line crosses by applying the EM-algorithm (DEMPSTER et al. 1977 ), where the evaluation of the expectation step was conducted by a Markov chain Monte Carlo (MCMC) technique analogous to that of GUO and THOMPSON 1992 . All these methods are based on the assumption that the distribution of the phenotypes is Gaussian. A robust method in this respect was developed recently by KRUGLYAK and LANDER 1995 . Modifications of (composite) interval mapping to binary traits were presented by VISSCHER et al. 1996A and XU and ATCHLEY 1996 , and a QTL mapping method for ordinal categorical traits by HACKETT and WELLER 1995 .

In interval mapping (LANDER and BOTSTEIN 1989 ), which is currently used routinely for QTL mapping (Mapmaker/QTL by LINCOLN et al. 1992 ), it is possible to calculate likelihood scores for a putative QTL placed in any position between two adjacent flanking markers. By changing the flanking markers one at a time, it is possible to determine the likelihood curve over the whole genome. The procedure is based on regression of phenotypes on QTL genotypes, and because QTL genotypes are unknown, results are obtained by using an iterative EM algorithm in which convergence to a local maximum is guaranteed. For each fixed location separately, the algorithm searches the parameter vector value giving the highest likelihood score. These (profile likelihood) scores are then used to draw the LOD-score curve corresponding to different QTL positions.

The composite interval mapping procedure introduced by ZENG 1993 , ZENG 1994 , and the multiple-QTL mapping of JANSEN 1993 and JANSEN and STAM 1994 , are in principle similar to interval mapping except that also some markers outside the tested interval (also in other chromosomes) are fitted to the model as covariates in order to reduce background noise caused by other QTLs and/or polygenic variation. The most significant markers for such reduction can be chosen in a preliminary analysis by using stepwise regression. These background control markers lie near the QTLs and they are used instead of the true QTLs (whose locations are unknown), yielding a better resolution than would be possible if those QTL effects were not considered at all. Theoretical support is provided by the fact that when a QTL has an effect on the trait one can see the same effect indirectly through the closest marker which is in linkage disequilibrium with the QTL (TANKSLEY 1993 ). The degree by which the corresponding phenotypic effect is reduced is determined by the recombination fraction between the marker and the QTL. However, there should be at least some effect, because linkage disequilibrium is almost at its maximum value in the backcross and F₂ intercross population. The improvement in resolution of the composite interval mapping over standard interval mapping is sometimes huge (e.g., KUITTINEN et al. 1997 ), thus yielding more putative QTL findings. However, the method does not yet seem to be widely used in practice, even though it is implemented in some software packages (e.g., QTL Cartografer by BASTEN et al. 1996, and MapQTL by VAN OOIJEN and MALIEPAARD 1996, and PLABQTL BY UTZ and MELCHINGER 1996 ).

The purpose of this paper is to present an approach for high resolution QTL mapping in inbred line crosses from incomplete data when the phenotypic distribution is assumed to be Normal. Our modeling approach is based on regression and the assumption that individual QTL effects are additive. The method belongs to the general framework of variable dimensional Bayesian models (e.g., GREEN 1995 ), applying MCMC algorithms (see APPENDIX A) to obtain the posterior distribution of the number of influential QTLs as well as estimating their locations in the (analyzed) chromosome and the corresponding phenotypic effects. Such a possibility was hinted at in a conference presentation by A. F. M. SMITH 1996 , and it has been explicitly implemented by SATAGOPAN and YANDELL 1996 , STEPHENS and FISCH 1996 , and UIMARI and HOESCHELE 1997 .

MCMC methods (METROPOLIS et al. 1953 ; HASTINGS 1970) are not new techniques, but their widespread application in statistics did not start before the introduction of Gibbs sampling (GEMAN and GEMAN 1984 ). For a review of applications in gene mapping see THOMPSON 1994 , THOMAS and GAUDERMAN 1995 and references therein. An excellent introduction to the Gibbs sampling and to the Metropolis-Hastings (M-H) algorithm can be found from CASELLA and GEORGE 1992 , and CHIB and GREENBERG 1995 . More advanced papers are GEYER 1992 , and BESAG et al. 1995 . Variable dimensional parameterizations are considered in ARJAS and GASBARRA 1994 , and GREEN 1995 .

Bayesian inference for gene mapping has been considered by TAI 1989 , THOMAS and CORTESSIS 1992 , SMITH and ROBERTS 1993 , and STEPHENS and SMITH 1993 . A Bayesian QTL mapping method was proposed recently by SATAGOPAN et al. 1996 , where a prespecified number of QTLs was assumed to be present in the considered chromosome. The method did not take into account effects of QTLs in other chromosomes. Judgement concerning the actual numbers of QTLs (in the model) was proposed to be made by using Bayes factors from separate MCMC runs with different numbers of QTLs in each. SATAGOPAN and YANDELL 1996 , and STEPHENS and FISCH 1996 considered all chromosomes simultaneously, treating the number of QTLs as a random variable. [The example in SATAGOPAN and YANDELL 1996 included only one chromosome, however.] SATAGOPAN and YANDELL 1996 used an M-H within Gibbs scheme in estimation, in contrast to an M-H scheme applied by STEPHENS and FISCH 1996 , as well as here. STEPHENS and FISCH 1996 simulated 10 chromosomes in 6 different datasets, with a different heritability in each. They also tested several priors. However, they did not consider missing data.

Bayesian QTL mapping in outbred livestock population, using granddaughter design, has been studied by UIMARI et al. 1996A , UIMARI et al. 1996B and UIMARI and HOESCHELE 1997 . A Gibbs sampler was used in single and bi-QTL models which included both major gene and polygenic effects. UIMARI and HOESCHELE 1997 also tested the idea of a variable dimensional model, considering the cases of either zero, one or two QTLs and using the convention that the second QTL, if present, was always "to the right" of the first one. They concluded that their method was sensitive to the way in which the parameters of the second QTL were chosen. Effects of major QTLs in the other chromosomes were not taken into account in these studies, because composite interval mapping related techniques in multi-generational pedigrees are difficult to apply. Apart from UIMARI et al. 1996A , techniques for handling missing genotypic data were not considered.

The contents of this paper are as follows. Next, we describe our statistical model. The results from simulation experiments are described thereafter. The final section contains a discussion of the method and of the experiences we had. In APPENDIX A we outline the MCMC algorithm used in the estimation.

	MODEL

TOP ABSTRACT MODEL SIMULATION ANALYSIS DISCUSSION APPENDIX A: ESTIMATION OF... APPENDIX B LITERATURE CITED

QTL search is usually concentrated on a given chromosome (or chromosomal segment). Therefore everything in the following is with respect to such a chosen fixed linkage group (chromosome) unless the contrary is stated. Let y = (y₁,y₂, ..., y_Nind) denote the vector of known phenotypes (missing phenotypes are not considered here), where N_ind is the number of individuals in the experiment and y_i is the phenotype of the ith individual. Suppose that the observations y_i are distributed according to a normal law resulting from some design in inbred linecross data. We shall view the unknown number of QTLs, denoted by N_qtl, as an unobserved random variable. Denote the QTL locations by l = (l₁,l₂, ..., l_Nqtl) and let be an N_qtl x N_ind matrix, where the qth column x_q = (x_q1,x_q2, ...,x_qNind)' is the QTL genotype vector in location l_q, with element x_qi referring to the ith individual. Let G* and G be the corresponding complete and incomplete (observed) marker information respectively; G* and G are taken to be N_ind x N matrices, where N is the number of markers. Let I be the chromosomal interval with the first and the last markers of the chromosome as endpoints. A fixed marker map (i.e., known recombination fractions between markers whose order is known), denoted by m, is assumed known before the analysis.

We denote complete and incomplete (observed) marker information in other chromosomes respectively by G^*_o and G_o. Let X^*_o be a subset of the complete marker information in other chromosomes, consisting of selected columns (marker genotype vectors) of G^*_o . Using an obvious set notation, X^*_o G^*_o . Similarly, let X_o be a subset of incomplete marker information, X_o G_o. We assume that X^*_o is chosen to correspond to known background control information (a selected set of markers that are hoped to be close to influential QTLs outside the interval I). Again, we arrange X^*_o into the form of an N_ind x N_bc matrix and denote its (i,k)th element by X^*_ik , corresponding to the genotype of the ith individual's kth background control. Here, N_bc is the number of background controls.

In the chosen design (i.e., experimental cross), let N_gen be the number of possible genotypes and let = (₁, ..., _Ngen) be the vector including all possible genotypes at any (marker or QTL) locus, ordered so that all homozygote genotypes come before heterozygotes. In this setting, genotypes AB and BA are considered to be the same. The regression parameters are the following: a is the regression intercept (mean value), b_q = (b_q1, ..., b_qNgen) is a vector of regression coefficients (classification variables) where b_qj is the regression coefficient for the qth QTL genotype _j at location l_q, ² = Var(e_i) is the residual variance and C = (c_kj) is an N_bc x N_gen matrix of regression coefficients c_kj for background controls. We consider the following statistical model for y:

(1)

Here 1_{xqi=j} is the indicator variable (dummy), which takes value one if the ith individual's qth QTL genotype x_qi at location l_q is _j, and otherwise its value is zero. Similarly, 1_{X^*ik=j} is the indicator variable taking value one if the ith individual's marker genotype in the kth background control is _j and otherwise its value is zero. Here we assume that the residuals e_i are independent and normally distributed according to N(0, ²). In order to maintain the traditional way of considering gene effects one can make the following convention. In the case of three possible genotypes, say = (AA,BB,AB), the constraint b_q2 = -b_q1 will produce an additive effect for homozygotes and make b_q3 correspond to a dominance effect of a heterozygote for each q. For the background control parameters, the constraint c_k2 = -c_k1 for k = 1, ..., N_bc will have a similar interpretation. In case of backcross, the corresponding constraints are b_q1 = 0 and c_k1 = 0, respectively.

We use the shorthand notation = (a,b₁, ..., b_Nqtl, ²,C) and = (,,l,G*,X^*_o,N_qtl) . Notation A* ~ A means that A* is consistent with A in cases where A is incomplete (observed) and A* is complete information. From Bayes' theorem, we get p(|y,G,X_o,m) = p(y,,G,X_o|m) , where the joint density p(y,, G,X_o|m) can be factored into a likelihood and a (joint) prior density as p(y,,G,X_o|m) = p(y,G,X_o|,m)p(|m). Here the likelihood can be further written into the form of the product p(y,G,X_o|, m) = p(G|,m)p(X_o|G,,m)p(y|,m,G,X_o) = 1_{G*~G}1_{X^*o~Xo}p(y|,m) . In this, 1_{G*~G} and 1_{X^*o~Xo} are indicators taking values one and zero depending on whether the complete genotypes G* in the chromosome and in the background control sites X^*_o (in the other chromosomes) are consistent with their observed incomplete counterparts or not.

As for the joint prior distribution p(|m), given the marker map m, we make the following (conditional) independence assumptions:

1. Given N_qtl, the vector consisting of the parameters of the linear model (1) is independent of the other coordinates of , i.e., of (,l,G*,X^*_o) , as well as of m;

2. X^*_o , the vector consisting of the true background control genotypes in other chromosomes, is independent of (,l,G*,N_qtl), the true marker and QTL-information in the considered chromosome as well as of the marker map m;

3. the true genotypes at marker locations, G*, are conditionally independent of the number of QTLs (N_qtl) and their locations (l) given the marker map m.

Then the joint prior density function can be factored further and be presented in the product form

(2)

The density p(X^*_o) is not conditioned on the fixed marker map, because the prior for genotypes can be thought not to be dependent on the marker order or distances. The posterior density of is then proportional to the joint density

(3)

(4)

In specifying prior densities it is important to note that the number of possible alleles and genotypes (N_gen) depends on the (experimental) design in question. Therefore, also the prior densities p(|G*,l,m,N_qtl), p(G*|m) and p(X^*_o) should reflect the design. Crosses between inbred lines have two alleles, forming two (three) different genotypes in backcross (F₂ intercross) designs.

Let M^*_s be the complete genotype vector in the sth marker position, i.e., the sth column in G*, and let M^*_s,i be its ith component. In case there are some unobserved genotypic data, i.e., G*\\G 0, we consider the following conditional independence structure for the prior p(G*|m): we assume that p(M^*_s|G^*_-s,m) = p(M^*_s|M^*_s-1,M^*_s+1,m) = ^Nind_i=1 p(M^*_s,i|M^*_s-1,i, M^*_s+1,i,m) , where G^*_-s includes all the other columns in G* except the sth. In a backcross design, the prior for an individual i with complete genotype information G* can be computed as the product

(5)

(6)

In the case of F₂ intercross, for each individual the transition probabilities p(M^*_s+1,i|M^*_s,i,m) from position s to s + 1 can be arranged into the 3 x 3 matrix containing all possible transitions between states AA, BB, and AB

(7)

The prior distribution of N_qtl, the number of QTLs, is here assumed to be truncated Poisson, where the Poisson mean and the maximum number N_qtlmax of QTLs are fixed control parameters such that 0 N_qtlmax. The upper bound N_qtlmax is introduced for computational reasons.

In the following, we shall use the generic term "object" for any marker or QTL in the chromosome. The prior distribution of QTL genotypes at locations l is assumed to have the following product form:

(8)

Here, G^*q_iL and G^*q_iR are the genotypes of the left and the right flanking object (marker or QTL) for the qth QTL in individual i, chosen among the complete set of the markers in the chromosome and the QTLs at positions l₁,.,l_q-1. When the location of a QTL and the corresponding flanking object genotypes are known, the QTL genotype is independent of the genotypes of other objects (markers or QTLs) in this list. We denote by r_q = (r_q1, r_q2) the resulting recombination fractions between the QTL at l_q and the corresponding flanking objects. As is often done in QTL mapping applications, the same recombination rates for male and female meioses are assumed also here. The algorithms for constructing the probabilities of different QTL genotypes (last term in Equation 8) under backcross and F₂ designs can be found in APPENDICES B1 and B2, respectively. In this construction, Haldane's map function is used for converting the distance between l_q and the left flanking object to a corresponding recombination fraction for each q. Haldane's formula r_q2 = then gives the recombination fraction between l_q and the right flanking object, with r^q_fm being the recombination fraction between the two flanking objects of the qth QTL. Note that Equation 8 allows more than a single QTL within the same marker interval. Apart from being more general than models in which at most one QTL is allowed, this feature is thought to improve the mixing properties of the MCMC sampling algorithm.

An obvious choice for the priors of all the QTL locations is the uniform distribution on the considered chromosome, corresponding to the assumption that no prior knowledge concerning the QTL loci is available. However, if some `non-data dependent' knowledge has been obtained, for example, using cytogenetical methods (e.g., physical exclusion mapping), one can specify a prior which has most of its mass on some narrower chromosomal area. (Note that the uniform prior densities here do not cause any integrability problems, because all chromosomes are of finite length.) We shall assume equal prior probabilities for background control genotypes and use uniform prior distributions for all regression parameters. The natural range for the residual variance, for example, is between zero and the phenotypic variance.

The main features of the proposed model are summarized graphically in Figure 1. Our main interest is in the number of QTLs and in their positions in the considered chromosome. In order to arrive at a meaningful description of the results from the estimation we consider the QTLs as forming a nonhomogenous Poisson process over the chromosome. The results of the statistical analysis can then be expressed in an intuitive and coincise manner in terms of the corresponding estimated intensity. In practice, we divide the chromosome into intervals (bins) ₁, ₂, ...,_Nbins of equal length (according to the Haldane distance). The interval length _j chosen by the analyst reflects the resulting mapping resolution. Denote the number of MCMC cycles (sampling iterations) by N_cycs, and let

(9)

{l^(k)_qj}

{sj}

View larger version (23K): In this window In a new window Download PPT slide   Figure 1. —Hierarchical structure of the model. Boxes refer to fixed values and ellipses to random variables. Layer one is observed, layer five given, and the others are unknown (sampled). Solid arrows indicate the direction of hierarchical dependency. Dotted arrows describe direct functional relationship.

For assessing phenotypic effects in the backcross design, let D_j(d) be the cumulative distribution function (c.d.f.) associated with the phenotypic effect of a putative QTL in bin _j. There will then be one such c.d.f. for each bin. In the backcross design, let

(10)

{l^(k)_qj, b^(k)_q1-µ^(k)_qd}

{l^(k)_qj, b^(k)_q3-µ^(k)_qd}

{l^(k)_qj, b^(k)_q1d}

{l^(k)_qj, b^(k)_q3d}

View larger version (44K): In this window In a new window Download PPT slide   Figure 2. —The complete data analysis. On the left, the results from interval mapping (IM, solid line) and composite interval mapping with five background controls (CIM/05, broken line) are shown. Simulated true QTL locations are indicated with an asterix (*). The histogram corresponds to the (approximate) posterior QTL intensity over the chromosome, with binlength 1 cM. The left (right) y-axis corresponds to the likelihood ratio statistic (posterior QTL intensity). On the right, corresponding phenotypic effect estimates are shown. The solid line is the posterior median and the grey lines are the 2.5- and 97.5 percent quantiles of the posterior distribution of the phenotypic effect of a putative QTL. Shaded regions are suggested credible intervals for QTL localization (see Table 5). The phenotypic effect estimates (posterior median and quantiles) are reliable only in these regions.

	SIMULATION ANALYSIS

TOP ABSTRACT MODEL SIMULATION ANALYSIS DISCUSSION APPENDIX A: ESTIMATION OF... APPENDIX B LITERATURE CITED

In order to test the performance of this method, a data set was simulated by the QTL Cartografer software (BASTEN et al. 1996 ). In particular, we wanted to compare its performance in the case of a simple backcross design to interval mapping (IM), and to composite interval mapping with five background controls (CIM/05). A backcross population (N_ind = 250) was generated of individuals having three 100-cM length chromosomes and in each 11 equally spaced markers 10 cM apart from each other. The trait was assumed to have heritability (the proportion of phenotypic variance explained by the simulated QTLs) 0.7 and phenotypic variance 1.0. The effects and the locations of the six simulated QTLs can be found in Table 1. A second data set was generated from this complete simulated set by randomly deleting 30 percent of the marker genotypes. The second set was used to test how well our method is capable of handling situations where a large proportion of genotypes are unknown.

First, the two data sets were analyzed by using the IM and CIM/05 methods. The background controls for the analyses were chosen by standard stepwise regression (QTL Cartografer software). The background control markers for the complete data were marker two in chromosome 1, markers zero, one, and six in chromosome 2, and marker two in chromosome 3. In the analyses where 30 percent of genotypes were missing, the background controls were marker two and three in chromosome 1, markers one and six in chromosome 2 and marker three in chromosome 3. The same background controls were also used in the corresponding Bayesian analyses. In the CIM, "window width" was 10 cM (i.e., the background controls less than 10 cM from the analyzed interval were not included in the model).

Several test runs were made in order to carefully adjust the range of the proposal distributions (i.e., parameters that control the maximal step size) of the Metropolis-Hastings algorithm. Such a range is specified for each of the following: (1) QTL locations, (2) regression mean, (3) residual variance, and regression coefficients of both (4) the QTLs and (5) background control genotypes. These control parameters influence directly the rejection rates; if they are chosen carelessly the chain may not convergence to the correct limiting distribution within a reasonable time. The values used in the final analyses are given in Table 2.

The (C-program implementing a Metropolis-Hastings-Green) chain was run 500,000 rounds for all analyses in chromosomes 1 and 2 and 1,000,000 and 1,500,000 rounds respectively for the complete and incomplete data in chromosome 3. Computations were made on an UltraSparc Model 170 workstation, with running times varying between 1 hr and 30 min, and 6 hr and 20 min, depending on the chromosome and on other work load on the computer. The initial value for the number of QTLs was three, and the corresponding locations were 20.0 cM, 50.0 cM, and 80.0 cM in all MCMC runs. The Poisson mean (hyperparameter) was set to = 2 and the maximum number of QTLs (in the analyzed chromosome) to N_qtlmax = 3. The prior of the residual variance was chosen to be uniform over the range [0.0, 0.89], the right endpoint being equal to the phenotypic variance estimate from the data. The prior of the intercept was taken to be uniform over [-100, 100], and that for QTL and background control genotypic regression coefficients uniform over [-2, 2]. In all cases the chosen ranges were certain to cover all realistic parameter values. Finally, the prior for QTL locations was uniform over [0, 100].

Results:
The likelihood ratio statistic (LRS) curves (in base 10 logarithmic scale) from IM and CIM/05 runs, and the Bayesian posterior QTL intensities in the considered three chromosomes, are shown on the left side of Figure 2 (complete data) and 3 (incomplete data). The phenotypic effect estimates given by the IM and CIM/05 methods, as well as the curves consisting of the pointwise (i.e., in different locations of the putative QTLs) medians and the 2.5- and 97.5-percent quantiles of the posterior distributions of phenotypic effects are shown in the same figures on the right. For an obvious reason, the phenotypic effect estimates deserve serious consideration only in those chromosomal regions in which the statistical analysis suggests that there actually might be a QTL. Judging by the level of the QTL-intensity, we have shaded such areas in the figures.

Approximate posterior probability distributions of the number of QTLs in different chromosomes, for both the complete and the incomplete data, are given in Table 3 and Table 4. The posterior expectation of the number of QTLs in each chromosome (shown in Table 3 and Table 4) coincides with the area under the corresponding QTL-intensity curve (in Figure 2 and Figure 3).

View larger version (47K): In this window In a new window Download PPT slide   Figure 3. —The incomplete (30% of genotypes missing) data analysis. On the left, the results from interval mapping (IM, solid line) and composite interval mapping with five background controls (CIM/05, broken line) are shown. Simulated true QTL locations are indicated with an asterix (*). The histogram corresponds to the (approximate) posterior QTL intensity over the chromosome, with binlength 1 cM. The left (right) y-axis corresponds to the likelihood ratio statistic (posterior QTL intensity). On the right, corresponding phenotypic effect estimates are shown: the solid line is the posterior median and the grey lines are the 2.5- and 97.5-percent quantiles of the posterior distribution of the phenotypic effect of a putative QTL. Shaded regions are suggested credible intervals for QTL localization (see Table 5). The phenotypic effect estimates (posterior median and quantiles) are reliable only in these regions.

Perhaps the most natural question to ask in this context is the following: "What is the (posterior) probability that some particular chromosomal area, say , contains at least one QTL, given the evidence in the data?" Denoting the number of QTLs in by N_qtl , one can calculate various MCMC approximations of that probability as shown in Equation 11 below, where (s) is the QTL-intensity at point s, and ^(s) is its estimator:

(11)

The last approximation, where the integral on the right is actually an expression for the (posterior) expected number of QTLs in , is reasonable only if the integral is small. Table 5 gives a few such approximations for some chromosomal areas, based on either the complete or the incomplete data. These regions represent a moderate to high posterior QTL-intensity region surrounding the mode. Making the interval longer will obviously always increase both the probability P(N_qtl 1 | data) and the expectation E(N_qtl | data) , but this will be at the cost of less accurate localization of the genes. In other words, given the evidence in the data, there is always a trade-off between accuracy and probability, just as in confidence intervals in classical statistical inference. In this sense, Table 5 represents only one possible summary of our findings. We have not made an attempt to establish a standard for forming such intervals or corresponding cut-off points here. Indeed, we think that the intensity curve in itself is the best summary of the information concerning the number of QTLs and their locations in the chromosome as obtained from the statistical analysis.

View this table: In this window In a new window   Table 5. Approximate (posterior) probability (1 - exp{ -I^(s)ds}) that a given chromosomal area contains at least one QTL calculated for different areas

We do not display empirical threshold values (corresponding to permutation tests) for IM and CIM in Figure 2 and Figure 3 because they were not the main issues here. However, the LOD score 3.0 would correspond to the likelihood ratio statistic 3.0 x (e) 13.82 , which can be used as a rule-of-thumb threshold when examining the figures (even though it was originally derived for monogenic traits in human). In practice the CIM method needs a somewhat higher threshold than IM (ZENG 1994 ). Note that the LRS curves and posterior QTL-intensity decrease sharply close to some marker positions. This is because genotype information concerning putative QTLs tends to be more accurate close to markers, and unless a QTL actually coincides with a marker locus, there will be strong evidence against placing a QTL in exactly, or very close to, these positions (see KUITTINEN et al. 1997 ).

In Figure 2, the IM and CIM curves, apart from CIM in chromosome 3, contain enough evidence (in the sense that the LRS is greater than the threshold value 13.82) of QTL activity. From Figure 2 it can be seen that our method performs well in all chromosomes, by giving well localized high intensities for QTLs close to their true locations. Note that the highest QTL-intensities are often found near the modes of CIM or IM curves. The weakest QTLs in chromosomes 1 and 2 remained undetected by any of these methods, apparently because of their small phenotypic effects. All three methods estimated the phenotypic effect of the most influential QTLs very accurately, but IM and CIM/05 do not give confidence intervals for the estimates [unless they are estimated separately with a permutation test (CHURCHILL and DOERGE 1994 ) or by bootstrapping]. Our method gave posterior credible intervals for the phenotypic effects in all cases.

In the case of incomplete data, the mode of the intensity given by our method differs by 12 cM from the simulated true location in chromosome 3 (Figure 3). This is apparently a consequence of reduced genotype information in the incomplete data set, as well as of the weak phenotypic effect attributable to this particular QTL. Note in this case the low maximal level of the posterior QTL-intensity, compared to the levels reached in the other two chromosomes.

For a more explicit comparison of the three methods, one-lod-support intervals (see OTT 1991 , pp. 66–67) were determined around the modes of the IM and CIM/05 curves estimated from the complete data. The threshold values defining the one-lod-support intervals are the maximal value (the mode) minus 1.0 LOD (i.e., (e) units in the LRS scale). An alternative would have been to estimate confidence intervals for the QTL locations, by applying a permutation test (see CHURCHILL and DOERGE 1994 ) or by bootstrapping, as in VISSCHER et al. 1996B . The results from the comparison are summarized in Table 6. (The numerical accuracy of the estimated support interval depends on the chosen step size which specifies how frequently the LRS is evaluated along the chromosome.)

View this table: In this window In a new window   Table 6. True locations (QTL), estimated locations and one-lod-support intervals for IM and CIM, and Bayesian point estimates (modes of the QTL intensity) together with the support intervals () from Table 5

Considering first the complete data set and the first chromosome with the LRS evaluated at every 0.1 cM, we obtained the one-lod-support interval [23.0 cM, 29.2 cM] using the IM, and [21.6 cM, 27.3 cM] using the CIM method. The latter almost coincides with the interval = [21 cM, 28 cM] (see Table 5 and Table 6), obtained from Figure 2 so that it contains practically all moderate to high posterior QTL-intensity values. All these intervals covered the true QTL at 24.15 cM. The posterior probability that the region contains at least one QTL is 0.63. The true phenotypic effect of the second QTL (at 72.21 cM) in chromosome 1 was so small that this QTL was not detected by any of the three methods discussed here.

At the "left end" of chromosome 2, we obtained the one-lod-support interval [2.9 cM, 9.0 cM] when using the IM, and [10.0 cM, 11.9 cM] when using the CIM method. The latter interval is actually so narrow that the true QTL at 7.77 cM falls just outside it. The posterior QTL-intensity in this region is concentrated almost completely on the interval = [4 cM, 10 cM] which again contains the true QTL. At the "right end" of chromosome 2, the LRS curve arising from CIM is bimodal, resulting in two overlapping one-lod-support intervals, [53.6 cM, 59.6 cM] and [52.9 cM, 67.3 cM]. The shorter CIM interval is so narrow that it does not contain the true QTL at 61.67 cM, but the wider one does. In this case, a natural choice for a Bayesian credible region would be the interval = [53 cM, 68 cM], which is very close to that obtained by CIM and for which the posterior probability of there being at least one QTL in the region is 0.64. The IM method failed to reveal any statistically significant QTL activity in this part of chromosome 2.

In chromosome 3, we obtained the one-lod-support interval [0.0 cM, 39.8 cM] with the IM method. With the CIM method, we could not determine a corresponding support interval because the maximum LRS value was less than the critical value LOD 3.0. Because of the small phenotypic effect of the simulated QTL, the Bayesian method did not localize this QTL as well as in chromosomes 1 and 2. Perhaps a natural interval to suggest, if at all, would be [0 cM, 37 cM], which would be nearly identical to the interval suggested by the IM-analysis and which would cover 90.2 percent of the total area under the QTL-intensity curve in that chromosome.

Whenever possible, all three methods estimated the locations of the QTLs fairly well (see Table 6). In all these analyses the Bayesian estimates (posterior medians) of individual phenotypic effects were close to the true values, and they were practically the same as what was obtained when applying the CIM method. In this respect, the IM method turned out to be much inferior.

The results from analysing the incomplete data were largely similar, see Table 5 and Table 6, and Figure 3 for details. The main difference was that the one-lod-support intervals become somewhat wider, as did the corresponding Bayesian credible regions. The differences were not very large, however.

	DISCUSSION

TOP ABSTRACT MODEL SIMULATION ANALYSIS DISCUSSION APPENDIX A: ESTIMATION OF... APPENDIX B LITERATURE CITED

We have presented here a novel method for high resolution mapping of multiple QTLs and for the estimation of their phenotypic effects, using the general framework of Bayesian variable dimensional models. For the estimation of the parameters (see APPENDIX A) we use the Metropolis-Hastings algorithm with "reversible jumps" (GREEN 1995 ) between models with different numbers of QTLs. In the case of more than one QTL, this construction improved the mixing properties of the sampler when compared to a single-QTL model (results not shown). The effects on the QTLs in other chromosomes are taken into account indirectly, through nearby markers. When the genotype of a marker is missing, information from surrounding markers is utilized by applying the conditional probabilities as specified in Equation 5.

The main advantage in using Bayesian, instead of the more traditional frequentist inferential methods, is that they enable the analyst to quantify probabilistically the uncertainty involved in each claim made about QTLs, without needing to use problematic mental constructs such as "the relative frequency of incorrect decisions made in a long sequence of trials repeated under similar conditions." The transformation of the prior distribution into the posterior through Bayes' formula, corresponds directly to the natural intuition of "learning from the data." Depending on the goals of the study, the specification of the prior can be "neutral" if the goal is to present a statistical summary of the information in the data, or, if available, it can also reflect an expert's prior knowledge about unobserved quantities, in which case the posterior will be a synthesis of such expert knowledge and empirical evidence coming from the data. Another major advantage of the Bayesian approach is the relative ease by which missing data (such as missing genotypes) problems can be handled, together with the estimation of all other unobservables. Finally, the application of MCMC methods gives considerable freedom in building large hierarchical statistical models corresponding to the analyst's perception of the underlying genetic structures and dependencies.

The results of the statistical analysis are summarized by two new measures: the posterior QTL-intensity, considered as a function of its location in the chromosome, and the posterior distribution of the phenotypic effect of the corresponding putative QTL. Such probabilistic summary measures seem to correspond directly to the immediate objectives of QTL mapping, that is, localizing the important QTLs in different chromosomes and estimating their effects on the phenotype(s). In particular, in situations where one can expect that there are several QTLs in the considered chromosomal area, the posterior QTL-intensity captures the essential information about their number and positions in an easily interpretable probabilistic form. In this way we can avoid completely the difficult inferential problems concerning the "correct" threshold values of a LOD score (or a LRS) which arise in testing multiple QTL hypotheses. Moreover, our method does not appear to produce false positives easily, as one can expect to get a low QTL-intensity in regions where there is no, or is only little, QTL activity.

The performance of our method was compared to interval mapping (IM) and composite interval mapping (CIM) by using a simulated backcross population of 250 offspring. A second data set was obtained by randomly deleting 30 percent of the marker genotypes in the complete set.

In the execution of the MCMC sampling we used an overparametrized regression model which has one extra coefficient for each QTL and for each background control locus. Therefore the model intercept and the genotypic coefficients are not identifiable as such, but their contrasts (phenotypic effects) are. We also tested the alternative updating scheme used by SATAGOPAN and YANDELL 1996 where, when the number of QTLs was proposed to be changed, its effect was first balanced against a corresponding change in the overall mean. Our overparametrized model seemed to have better mixing properties, however.

The MCMC methods require some amount of computational effort and will, in practice, need the capacity of a workstation. To ensure sufficient mixing, we performed some relatively long test runs before a final QTL analysis. Another possibility is to apply some diagnostic tools, such as CODA (BEST et al. 1995 ). Because we ran long simulation trials, no sampled values were rejected because of burn-in. The mixing properties of the sampling algorithm do not seem to be very sensitive to the prespecified proposal probabilities. In the case of adding or deleting a QTLs, the only restriction appeared to be that the proposal probabilities of changing the dimension should not be too small. As an illustration of the degree of mixing which was typically encountered, we show (Figure 4) how N_qtl was varying in a simulation run of length 1,000,000. In the final runs, the observed rejection rates for both adding and deleting steps varied between 0.998 and 0.999. By comparison, the rejection rates for updating QTL locations were rather low in general because the proposals were usually limited to a fairly narrow interval around the current value. (Larger jumps were allowed when a QTL was either deleted or added.)

View larger version (70K): In this window In a new window Download PPT slide   Figure 4. —The sample path k N(k)qtl from a simulation trial of 1,000,000 iterations with complete data set (Chromosome 3).

An important issue which has so far come up only implicitly in our analysis is that the number of QTLs and their phenotypic effects cannot be considered in isolation from each other. The idea that any gene with a non-zero effect on phenotypes is a QTL which in principle could be detected from data seems like an idealization far from reality, and as a consequence, "the correct number of QTLs" will exist as an objectively defined quantity only in simulated data. Here we have controlled this problem by bounding the number of QTLs in each chromosome by a given constant. An alternative would be to modify our method by paying attention only to influential QTLs, in the sense that their phenotypic effect exceeds some given threshold T. The only change this would make into our formulas is that, in Equation 9 and Equation 10, we would have to add into indicators the restriction |b^(k)_q2 - b^(k)_q1| T .

Another question we have not discussed is how to choose the markers which are to be used as covariates in the analysis. Instead of using the results from a preliminary stepwise regression analysis as we have done here, one could "learn" from Bayesian QTL analyses of other chromosomes and choose certain markers from high posterior QTL-intensity areas as QTL representatives (background controls). Alternatively, one could choose covariates from amongst a set of candidate genes in case such genes are available. A final possibility would be to consider the entire genome in a single variable dimensional QTL analysis, always using the "current" QTLs as controls (SATAGOPAN and YANDELL 1996 ; STEPHENS and FISCH 1996 ). For this our method would need some adjustments, however.

An initial version of the program source code (written in C language) is freely available for research purposes from Rolf Nevanlinna Institute's web page (http://www.rni.helsinki.fi/~mjs). A more user-friendly documented version is currently being developed. The present framework will be extended later to cover outbred linecross and (human) pedigree data. Epistatic effects (interactions between QTLs) and multiple trait analysis would also be worth considering in the future.

	ACKNOWLEDGMENTS

M.S. thanks KARI AURANEN, JUKKA RANTA and DARIO GASBARRA for many useful discussions about Metropolis-Hastings algorithms, and MATTI TASKINEN for helpful hints in the programming and computer work. We are also grateful to OUTI SAVOLAINEN, LEENA PELTONEN, CLAUS VOGL, JANNE PITK¨ANIEMI, PEKKA UIMARI and four anonymous referees for their constructive comments on an earlier version of the paper. This work was supported by a research grant from the Academy of Finland.

Manuscript received April 7, 1997; Accepted for publication November 17, 1997.

	APPENDIX A: ESTIMATION OF MODEL PARAMETERS VIA MARKOV CHAIN MONTE CARLO

TOP ABSTRACT MODEL SIMULATION ANALYSIS DISCUSSION APPENDIX A: ESTIMATION OF... APPENDIX B LITERATURE CITED

Markov chain Monte Carlo (MCMC) methods enable one to calculate expectations with respect to the posterior in an approximate manner in situations where computation by traditional methods, because of the high dimension of the parameter space, would be complicated or impossible. In MCMC, each expectation is approximated by a sample average where the sample is drawn by simulation from an ergodic Markov chain constructed in such a way that its limiting distribution coincides with the posterior. The fact that the number of QTLs has not been fixed in advance, and is actually estimated in conjunction with the other model parameters, leads us to consider this problem within the general framework of variable dimensional parameter estimation.

We have preferred to use the Metropolis-Hastings algorithm (M-H) instead of its special case, the Gibbs sampler. The main motivation behind this choice has been that the M-H algorithm is so easy to implement, without needing to work with analytically involved full conditional distributions. An additional bonus of M-H is its greater flexibility and the relative ease by which it can be extended to more complex designs and pedigree structures.

For algorithmic reasons (mixing properties), we do not use any constraints for genotypic coefficients, even if this kind of overparametrization represents an uncommon modeling practice. Genotypic coefficients are here not identifiable as such, but their contrasts can be estimated.

Following GREEN'S 1995 construction of a Markov chain with reversible jumps, we apply acceptance probabilities of the form min{ 1 , (posterior ratio) x (proposal ratio) x (Jacobian)}. Throughout our consideration here, the Jacobian is one, being the determinant of an identity matrix. This is because the location(s) of existing QTL(s) do not determine the location of a proposed new QTL, nor does the deletion of one QTL influence the position(s) of the remaining QTL(s). In the following, we describe the steps of our Metropolis-Hastings-Green cycle. Initial values of the number of QTLs, of their locations, as well as the Poisson mean , the maximal number of QTLs allowed and the ranges of the uniform priors, are all given by the analyst (see simulation analysis). Initial values of the regression parameters are generated to be close to the center of their prior range. Initial values for the QTL genotypes and missing marker genotypes are generated from their priors. Figure 1, showing the hierarchical structure of the model, may help to understand the relationships between the parameters in the following updating scheme.

Step 1:
We consider here three different move types: (1.1) modify the location(s) and configuration(s) of existing QTL(s), (1.2) add one QTL to the model, and (1.3) delete one QTL from the model. These move types have proposal probabilities p_m, p_a, and p_d, respectively, such that p_a = 1_{{Nqtl<Nqtlmax}}c, p_d = 1_{Nqtl>0}c , and p_m = 1 - p_a - p_d . Here c is a given positive constant in [0, ]. In step 1.1 we do not fix the order of QTLs unlike in RICHARDSON and GREEN 1997 .

Step 1.1:
N^(t)_qtl: = N^(t-1)_qtl . The following cycle is repeated for each QTL, q = 1, ..., N^(t)_qtl : A new proposal for a QTL location, l^new_q , is sampled from a symmetric uniform density around the previous value. The proposal is accepted with probability

When nothing else is changed, the likelihood remains unchanged even if the proposal is accepted and therefore these two likelihoods cancel from the acceptance probability expression. If a proposal is accepted then l^(t)_q = l^new_q , and otherwise l^(t)_q = l^(t-1)_q . New QTL genotypes are proposed separately for each individual so that x^new_qi is sampled from p(x_qi | G^*(t-1)_qi,l, G^*(t-1)_qi,r, r_q) (constructed as show in APPENDIX B). Individual proposals are accepted with probabilities {1, }, where the likelihoods L_i1 = p(y_i | ^(t-1), x^new_qi, ^(t-1)_-qi, l^(t), X^*(t-1)_i,o, N^(t)_qtl, m) and L_i2 = p(y_i | ^(t-1), x^(t-1)_qi, ^(t-1)_-qi, l^(t), X^*(t-1)_i,o, N^(t)_qtl, m) for individual i are evaluated at the new and the old QTL genotypes, respectively. If the proposal for individual i is accepted then x^(t)_qi = x^(new)_qi , and otherwise x^(t)_qi = x^(t-1)_qi . Here ^(t-1)_-qi includes all except the qth QTL genotypes in individual i in round t for QTLs having indices lower than q and in round t - 1 for higher indices.

Step 1.2:
Add one new QTL. A proposal (N^(t)_qtl = N^(t-1)_qtl + 1) for the number of QTLs is made. The new QTL location, l_{(N^(t-1)qtl+1)} , is proposed from the uniform density on I. The QTL genotypes at location l_{(N^(t-1)qtl+1)} are proposed from p(x_{(N^(t-1)qtl+1)}|G^*(t-1),x^(t-1)₁,., x^(t-1)_Nqtl, l^(t-1), l_{(N^(t-1)qtl+1)},m) . The regression coefficients of the new QTL genotypes are drawn from their priors. The proposal is accepted with probability

The term in the denominator is squared only because we are not fixing the order of the QTLs (p. 23 in GEYER 1996 , where the order of quantities is fixed).

Step 1.3:
Delete one QTL. The proposal is made that the number of QTLs is decreased from N^(t-1)_qtl by one, each of the deletions being equally likely. The deletion is accepted with probability

Step 2:
The marker genotype proposals G^*_new , denoted by G^*_i,new for individual i, in the chromosome are sampled for all individuals from the distribution, where all the consistent genotypes are considered as equally likely. p(G^*_new | m) is evaluated for all the individual according to Equation 5. The marker genotype proposals are accepted for each individual i separately with probability

[Note that p(M^*_1,i) for each individual i in Equation 5 cancels from the Hastings ratio]. If the proposals for individual i are accepted, then G^*(t)_i = G^*_i,new , and otherwise G^*(t)_i = G^*(t-1)_i .

Step 3:
New proposals for the regression parameters are sampled from the symmetric uniform densities around their previous values (random walk). Denoting the likelihoods L₁ = p(y|^new,^(t),l^(t),X^*(t-1)_o,m,N^(t)_qtl) and L₂ = p(y|^(t-1),^(t),l^(t), X^*(t-1)_o, m,N^(t)_qtl) , the proposals (in the prior range) are accepted simultaneously with probability min{1,L₁/L₂}. If new regression parameter values are accepted, then ^(t) = ^new, and otherwise ^(t) = ^(t-1).

Step 4:
The genotype proposals for the background control markers in other chromosomes are sampled directly from the uniform prior density p(X^*_o) for all the individuals. All proposals which are consistent with the incomplete observations are considered as equally likely. The proposals are accepted, separately for each individual i, with probability min{1,L₁/L₂}, where the corresponding likelihoods L₁ = p(y_i | ^(t),^(t)_i,l^(t),X^*new_i,o,m,N^(t)_qtl) and L₂ = p(y_i|^(t),^(t)_i,l^(t),X^*(t-1)_i,o, m,N^(t)_qtl) are evaluated at the new and the old parameter values, respectively. If the proposals for individual i are accepted, then X^*(t)_i,o = X^*new_i,o , and otherwise X^*(t)_i,o = X^(t-1)_i,o .

Step 5:
Go back to the start, Step 1, until a prespecified number of cycles has been reached.

	APPENDIX B

TOP ABSTRACT MODEL SIMULATION ANALYSIS DISCUSSION APPENDIX A: ESTIMATION OF... APPENDIX B LITERATURE CITED

This appendix contains pseudo code algorithms for calculating the conditional probabilities of different QTL genotypes given the flanking objects (Equation 8) for each individual i in backcross (APPENDIX B1) and F₂ intercross (APPENDIX B2) designs. The main idea of the algorithms is described in detail in the backcross case.

B₁:

• for i = 1, ..., N_ind

• total = 0

• for j = 1, ..., N_gen