A Bayesian Algorithm for Functional Mapping of Dynamic Complex Traits

Functional mapping of dynamic traits measured in a longitudinal study was originally derived within the maximum likelihood (ML) context and implemented with the EM algorithm. Although ML-based functional mapping possesses many favorable statistical properties in parameter estimation, it may be computationally intractable for analyzing longitudinal data with high dimensions and high measurement errors. In this article, we derive a general functional mapping framework for quantitative trait locus mapping of dynamic traits within the Bayesian paradigm. Markov chain Monte Carlo techniques were implemented for functional mapping to estimate biologically and statistically sensible parameters that model the structures of time-dependent genetic effects and covariance matrix. The Bayesian approach is useful to handle difficulties in constructing confidence intervals as well as the identifiability problem, enhancing the statistical inference of functional mapping. We have undertaken simulation studies to investigate the statistical behavior of Bayesian-based functional mapping and used a real example with F2 mice to validate the utilization and usefulness of the model.


Introduction
Because of polygenic control and environmental modification, many traits of agricultural, biological and biomedical importance vary in a quantitative way [1].For this reason, the genetic analysis of these so-called quantitative traits has been one of the most difficult issues in genetic research for many decades.With the advent of powerful molecular biotechnologies and statistical methodologies, it has now been possible to precisely dissect the phenotypic variation of a quantitative trait into individual loci at the molecular level (known as quantitative trait loci or QTLs) and understand the effects of interactions between QTLs and environments on the trait phenotype [2][3][4].Thanks to tremendous efforts by geneticists worldwide, hundreds of thousands of QTLs have now been identified for a variety of complex quantitative traits in plants, animals, and humans.Meanwhile, a vast body of literature has been available about the development of statistical methods for QTL mapping [5][6][7], although most of these methods ignore the developmental or dynamic features of a trait in time course.
More recently, a library of statistical models, called functional mapping, has been developed to map QTLs that control dynamic traits [8,9].Functional mapping incorporates fundamental biological principles behind trait growth and development into a mapping framework.It capitalizes on the mathematical aspect of a ubiquitous growth law that every biological trait experiences a growth and developmental change with time through altering the ratio of the anabolic or metabolic rate and the rate of catabolism [10].The advantages of this approach lie in its increased biological relevance by embedding biological principles into the estimation process and its flexibility to generate a number of quantitative testable hypotheses about the developmental and genetic regulation of growth.From a statistical perspective, functional mapping estimates parameters that determine the shape of a genotype-specific growth curve and/or the covariance structure, thus strikingly increasing its statistical power for QTL detection.
Functional mapping was originally constructed within the maximum likelihood (ML) context, incorporated by a finite mixture model and implemented with the EM algorithm [11,12].Although ML-based approaches have many favorable statistical properties for parameter estimation, they may often face some significant drawbacks when applied to functional mapping.First, functional mapping usually uses nonlinear mathematical equations to model the mean and covariance structures, it is extremely difficult or eventually impossible to derive log-likelihood equations for nonlinear parameters in the maximization step.Second, functional mapping concerns genetic analyses of longitudinal data whose intrinsic highdimensional complexity makes computation increasingly prohibitive, especially when permutation tests are used to determine critical thresholds.Third, but not specific to functional mapping, ML-based approaches do not automatically provide confidence interval estimates of the estimators, and thus affect the inference of parameter estimation.
Many of the problems described above for ML can be overcome by using a Bayesian method.In ML, the unknown parameters are treated as unknown variables (unobservables) and the likelihood function is maximized in these variables.In the Bayesian paradigm, each unobservable parameter is given a prior distribution, and we then infer the posterior distribution of each unobservable conditional on the data (the observables).The summary statistics of the posterior distribution, e.g., the mean, the mode or the median, can be regarded as Bayesian estimates of unobservables [13].The interval estimate can be obtained simply by examining the posterior distribution.The mean of this marginal posterior distribution is a candidate Bayesian estimator of an unknown parameter.Although this marginal distribution rarely has an explicit form, and numerical integration is often prohibited because of high dimensionality of parameters, a Markov chain Monte Carlo (MCMC) algorithm can be used to simulate the sample from the joint posterior distribution.The potential of the Bayesian approach implemented with the Gibbs sampler or Metropolis-Hastings algorithm for QTL mapping has been explored for various genetic designs [14][15][16][17][18]. Yang and Xu [19] incorporated a Bayesian approach to map QTLs involved in a dynamic trait based on nonparametric Legendre polynomials, although they did not take a full advantage of functional mapping in quantifying biologically meaningful hypothesis tests about the genetic control of development and modeling the structure of the covariance matrix.For such biological processes as growth curve, HIV dynamics and pharmacodynamics, in which explicit mathematical functions exist to specify their dynamic changes, functional mapping based on parametric modeling has proven to be powerful for asking and addressing the biological questions at the interplay between genetic actions and developmental patterns.
In longitudinal data analysis, parametric covariance modeling has several advantages compared to conventional multivariate approaches ignoring the covariance structure.The most significant advantage is that parametric modeling generally results in more efficient estimation of the covariance matrix (and therefore the mean structure).Also, it can deal more effectively with data when the number of measurement times is relatively large.Lastly, it can handle the data more effectively in the cases that the measurement times are not common across subjects.For those reasons, the development of explicit parametric models for the variance-covariance structure has been an important topic over the last two decades [20,21].
In this article, we will develop a general Bayesian framework for functional mapping of complex dynamic traits based on parametric modeling of the mean-covariance structures.This framework is constructed by a mixture model in which multiple mixture components corresponding to the genotypes of the underlying QTLs are involved.We will implement the MCMC algorithm to estimate the posterior distribution of each parameter contained within the mixture model including those that define curve shapes and the covariance structure.A real example for the F 2 mouse progeny is used to demonstrate the utilization of the model and validate its usefulness in a practical genomic project of dynamic QTLs.We perform simulation studies to investigate the statistical properties of the Bayesian functional mapping model in terms of its convergence rate, estimation precision and power for QTL detection.

Linear Model
Suppose there is an F 2 population of size n, initiated with two inbred lines.In this F 2 population, many markers are genotyped to construct a genetic linkage map, aimed to identify the QTLs that affect growth curves.For each progeny, a particular growth trait, such as body weight, tail length or cell number, is measured at a series of time points (1, • • • , T ).Consider a putative QTL with genotypes qq (denoted as 0), Qq (denoted as 1) and QQ (denoted as 2) that affects the shape of growth curves.At a specific time point t, the phenotypic value of the growth trait for progeny i due to the QTL may be expressed by a linear model, i.e., where ξ ij is an indicator variable for progeny i carrying a QTL genotype and defined as 1 if a particular QTL genotype j is indicated and 0 otherwise, u j (t) is the expected phenotypic value for QTL genotype j at time t, and e i (t) is independently and identically distributed as N (0, σ 2 (t)).Note that measurements within an individual are likely to be correlated across times, with covariance σ(t 1 , t 2 ) between times t 1 and t 2 (t 1 , t 2 = 1, . . ., T ).These variances and covariances form a (T × T ) matrix Σ.

Modeling the Mean-Covariance Structures
Functional mapping models u j (t) by a biologically meaningful mathematical equation and the timedependent covariance matrix composed of σ 2 (t) and σ(t 1 , t 2 ) by statistically robust approaches.For example, the growth of a living entity can be defined as the irreversible increase of size with time.A series of mathematical models have been proposed to describe growth curves.Among these models, the sigmoidal (or logistic) growth function is regarded as being nearly universal in living systems to capture age-specific changes [10].Specifically, this S-shaped curve for a QTL genotype j can be mathematically expressed as: where Θ j = (α j , β j , γ j ) is set of parameters that determine curve shape of genotype j.Parameter α is the limiting value of growth as time t goes to infinity, α/(1 + β) is the initial value of growth, and γ is the relative growth rate.Many biologically important features of growth curves, such as the time of maximal growth rate and the duration of maximal growth, can be described by these parameters or their combinations.
A number of statistical methods have been derived to model the structure of covariance between longitudinal measurements [21].The first-order autoregressive (AR(1)) model has been applied to model the structure of the within-subject covariance matrix for functional mapping.This model uses two simplified assumptions, i.e., variance stationarity -the residual variance (σ 2 ) is constant over time, and covariance stationarity -the correlation between different measurements decreases proportionally (in ρ) with increased time interval.These two assumptions facilitate the computation of functional mapping.Wu et al. [22] embedded Carrol and Rupert's [23] transform-both-sides (TBS) model into the growthincorporated finite mixture model, in order to reduce the heteroscedasticity of the residual variance and preserve original biological means of curve parameters after the data are transformed.Núñez-Antón and colleagues proposed a series of so-called structured antedependence (SAD) models to approximate the age-specific change of correlation in the analysis of longitudinal variables [20].The SAD model has been employed in several studies and displays many favorable properties for genetic mapping of dynamic traits [24].

Likelihood
Although only phenotypic values and marker genotypes are observable, the probability distribution of the QTL genotypes in a mapping population can be expressed in terms of the location of the putative QTL, the marker genotypes, and the distance between the markers.Let y i be the phenotypic value of progeny i at different time points, M i = {M ik } m k=1 be the m-marker genotype of progeny i, λ be the QTL position measured by the distance of the QTL from the first marker of an ordered linkage group, and D = {D k } m k=1 be the distances between markers 1 and k.Suppose this QTL is located between markers k and k + 1.Then, we use ω j|i to denote the distribution of the QTL genotype of progeny i, i.e., which is the conditional probability of QTL genotype j given the marker genotype of progeny i.
The likelihood of parameters λ, Θ and Σ given observations can be written as: where y = {y i } n i=1 , Θ = {Θ j } 2 j=0 and π is assumed to follow a multivariate normal distribution with the genotype-specific mean vector expressed as: (5) and covariance matrix Σ.

Estimation theory
We will derive a Bayesian approach to estimate the unknown parameters.This approach needs specifying prior distributions for the unknowns.Given the data and the priors specified, the posterior distribution over all unknown parameters can be obtained by using Bayes' theorem.Let Q = {Q i } n i=1 denote the QTL genotypes for all n progeny.Then, the posterior density of λ, Θ, Σ, and Q is given by: where π(y|Q, Θ, Σ) = π(y i |Q i = j, Θ j , Σ) denotes the probability mass of the observation y given the QTL genotypes, π(Q|λ) = π(Q i |λ) is the probability mass of the QTL genotypes of all n progeny given their marker genotypes (M) and the QTL position (λ), and π(λ, Θ, Σ) is the prior imposed on the genetic parameters.It is reasonable to assume the priors are independent for the parameters.Thus, we have: The priors can be chosen from related studies.In principle, if there is reliable information for parameters, such as Θ, priors with a small variance may be used.For a parameter without enough information used to determine the prior, such as Σ and λ, noninformative priors or priors with a large variance may be utilized.Here, we choose a uniform distribution on [0, D m ] as a prior of λ.Multivariate normal priors with moderate variances are used for Θ j can be obtained.The standard prior distribution for the inverse of the covariance matrix Σ −1 is the Wishart (R, ρ) [25,26], where the so-called scale matrix C = R −1 represents prior structural information about Σ and ρ is the degree of freedom, greater than T − 1.A small value of ρ gives a relative flat distribution.The Wishart prior with low degrees of freedom and a specified R is regarded as a reference (or noninformative) proper prior.Despite less flexibility, this distribution offers the advantage of being a conjugate prior, leading to a relative simple form for the posterior.
Theoretically, the marginal posteriors of the unknown parameters can be obtained from the joint posterior ( 6) by integrating over the other unknowns.Unfortunately, in practice, the evaluation of such high-dimensional integrals in a closed form is not possible.However, it is straightforward to derive either full conditional posterior distributions for some parameters, or the explicit expressions that are proportional to the corresponding full conditional posterior distributions for other parameters, i.e., and Wi(D −1 , n + ρ) (10) where Θ −j = {Θ j : j = 0, 1, 2, j = j}, y ij contains the observations from those individuals with genotype j, and

Algorithm implementation
Within the Bayesian framework, with the explicit expressions described in the preceding section, a Markov chain Monte Carlo (MCMC) technique can be used to draw samples from the joint posterior distributions of unknown parameters (6).We will use a hybrid scheme of Gibbs sampler and Metroplolis-Hastings (M-H) algorithm [27,28].In particular, Gibbs sampling steps update Σ −1 , while the M-H algorithm updates Θ j and the QTL position.After generating a random sequence of states , the final MCMC samples are collected, from which we are able to make the inference about the unknown parameters.The construction of the Markov chain process was described in Appendix A.

Estimation issues
For parameters λ and Θ, the empirical means of their marginal posteriors are the Bayesian estimators.As justified by Tierney [28], for any unknown parametric family f , which is a square integrable with respect to the stationary distribution π, we have: under the assumption that (λ ) are the samples from the Markov chain.In other words, the empirical averages of the corresponding MCMC samples may be regarded as the consistent estimators for the unknown parameters.
The marginal posterior densities of these parameters are determined from the kernel density estimator [29], since the closed form of their full conditional posteriors are not available.For example, the histogram kernel density estimator for λ is given by: A second important estimation issue is to obtain the confidence intervals for the unknowns.Box and Tao [30] suggested that highest posterior density (HPD) regions can be constructed to give the confidence intervals for the parameters of interest and the detailed method for developing an approximated HPD via MCMC samples can be seen in Ritter and Tanner [31].Alternatively, we can obtain the approximate HPD for the parameters directly by from their corresponding smooth density estimators.
Finally, in order to estimate the Monte Carlo error via the central limit theorem, Geyer [32] suggested three types of consistent estimators, the window estimators, the method of standardized time series and the specialized Markov Chain estimators.Among them, the windows estimators probably provides the best estimates, although it requires stronger regularity conditions for consistency.

Structuring the Covariance Matrix
In longitudinal data analysis, statistical modeling of the covariance structure is generally desirable, given that time-dependent variances and correlations follow a certain pattern.Below, we incorporate two commonly used approaches for the covariance structure in functional mapping into the Bayesian framework.
The closed forms of the full conditional posterior distributions of σ 2 and ρ are not available, but the explicit expressions that are proportional to these posteriors can be derived, which are: for σ 2 , and: for ρ.Based on these expressions, the corresponding Metroplolis-Hastings steps can be developed to update σ 2 and ρ within the MCMC estimation scheme (Appendix B).

Structured Antedependence Model
According to Gabriel [33], the error at the current time depends on those at the previous times and the innovative error generated at the current time.Such a structured antedependence (SAD) model can be described by antedependence coefficients and innovation variances.The first-order SAD or SAD(1) model only contains one antedependence coefficient (φ).Assuming a constant innovation variance (ν 2 ), Jaffrézic et al. [34] derived a SAD(1)-structured covariance as: It can be seen that neither the variance nor the correlation function is stationary for the SAD(1) model, i.e. the variance can change with time, and the correlation does not depend only on lag time We pose an inverse-gamma prior on innovation variance ν 2 and a normal prior on antedependence coefficient φ.In a real data analysis, the priors are selected as π(ν 2 ) = IG(α, β) = IG(1, 1) and π(φ) = N (µ φ , η φ ) = N (0, 10).The explicit expressions that are proportional to the full conditional posteriors of these two parameters can be derived as: for ν 2 , and: for φ.
The detailed M-H steps for updating ν 2 and φ are given in Appendix C.

Bayes Factor
The estimate of the number of QTLs is a very important but difficult issue for QTL mapping.In ML, the QTL number can be estimated on the basis of an assumed number by changing the dimensionality of the mode.The best QTL number is chosen using some model-selection criterion such as Akaike's AIC or by hypothesis testing.However, the theory underlying AIC could not apply in the mixture mode context because of the absence of a single dominating local maximum in the likelihood.The unsuitability of the theory used to construct the null distribution of the likelihood ratio statistic for mixture model means that hypothesis testing is not straightforward in this context.Moreover, because the number of QTL (κ) cannot be estimated as a parameter, it is not possible to estimate the sampling error and confidence interval of the estimate for QTL number.But in the Bayesian paradigm, inference is based on the posterior distribution of the parameters.Inference for κ is then based on the marginal posterior distribution, P r(κ = l|y), l = 1, 2, ..., .
Models κ = l 1 and κ = l 2 may be compared via the ratio of their posterior probabilities: where the ratio of marginal probabilities of y, B l 1 l 2 = P (y|κ=l 1 ) P (y|κ=l 2 ) , is known as the Bayes factor for comparing κ = l 1 with κ = l 2 .The Bayes factor does not depend on the prior distribution of κ.Recently, the Bayesian analysis of mixtures with an unknown number of components has received great attention [13,35,46].In practice, a Bayes factor larger than 100 can often be regarded as an evidence for the preference of l 1 QTLs over l 2 QTLs.

Mapping Population
Cheverud et al. [36] constructed a linkage map with 76 microsatellite markers for 535 F 2 mice derived from two strains, the Large (LG/J) and Small (SM/J).The total length of this map is ∼1,500 cM (in Haldane's units) and an average marker interval length is ∼27.5 cM.The same experiment was repeated by Vaughn et al. [37] in which 502 F 2 mice were generated and a linkage map of 1,780 cM long was constructed.In both experiments, each F 2 progeny was measured for their body mass at 10 weekly intervals starting at age seven days.The raw weights were corrected for the effects of each covariate due to dam, litter size at birth, parity and sex.We combine these two F 2 populations for QTL mapping of growth trajectories.Overall, about 10% of the marker genotypes were randomly missing.The mice with missing data were excluded from the analyses.[37] data set by using functional mapping with maximum-likelihood based methods.They showed that body masses in the F 2 mice follow a logistic curve which can be described by Equation 2, but display substantial variation in the shape of curves.Bayesian-based functional mapping was used to genome-wide search for the QTLs that control mouse growth curves by estimating the QTL location (λ), genotype-specific curve parameters (Θ), and the structure or unstructured covariance matrix (Σ).The prior for the locations of s QTLs (λ 1 , ..., λ s ) along each chromosome is simply assumed to be uniform over [0, D m ].Maximum likelihood estimates of growth curve parameters by Zhao et al. [24,37] provide the information about priors of these parameters.The priors for Θ j 1 j 2 ...j s (j 1 , ..., j s = 0, 1, 2) are a multivariate normal, centered at (30, 10, 0.6) T with a large dispersion diag(9, 4, 1).The prior for the covariance matrix is set to be inverse-Wishart(R −1 , ρ), where R is given by the sample covariance matrix.

Zhao et al. [23 38] first analyzed Vaughn et al.'s
We used the distance of 40 cM from the first marker on each chromosome as the initial value of λ.The initial values of the other parameters were generated from their priors.For each analysis, the Markov chain was run for 60,000 cycles after discarding the first 2,000 cycles for the burn-in period.In order to reduce the serial correlation among samples, the chain was thinned by saving one iteration in every 60 cycles so that the total number of samples stored for calculating the targeted posterior distributions of the unknowns was 1,000 [26].
Bayesian functional mapping was used to genome-wide scan for the existence of QTLs for body mass growth trajectories in the F 2 mouse population.Estimated marginal posteriors of the QTL locations over all 19 mouse chromosomes are illustrated in Figure 1. from which the posterior peaks were observed on chromosome 6, 7, 10, 11, and 15.By calculating the Bayes factors of the model with equation 28 by assuming one QTL over the model assuming no QTL, chromosomes 6, 7 and 10 were each detected to harbor a QTL given that the logarithmic scaled Bayes factors are 12.91, 13.47 and 7.99 for the three chromosomes, respectively.The estimated locations of these QTLs are between between markers D6nds and D6M it58 on chromosome 6, marker D7nds1 and D7M it17 on chromosome 7 and marker D10M it133 and D10M it14 on chromosome 10.  1 tabulates Bayesian estimates for the locations of the three QTLs detected and their genotypespecific curve parameters from the marginal posterior means of these parameters.From the estimated 95% equal-tail confidence intervals as HPD regions, Bayesian estimates are considered to be reasonably precise (Table 1).
The estimated curve parameters were used to draw the growth curves of three different genotypes at each of the three detected QTLs (Figure 2).The three QTLs exert an effect on body mass growth curves in a similar pattern.The homozygote for the LG/J allele has the best growth during the time course of measurement, followed by the heterozygote for the LG/J and SM/J alleles and the homozygote for the SM/J allele.Based on quantitative genetic theory, we can partition time-dependent genotypic values into time-dependent additive (a(t)) and time-dependent dominance effect components (d(t)) for each QTL, illustrated in Figure 3.The additive effect due to the substitution of the SM/J allele with the LG/J allele is positive and increases with age for all the three detected QTLs.The dominant effect due to the interaction between the LG/J and LG/J alleles is consistently small in time course.

Monte Carlo Simulation
Simulation studies were performed to study the statistical behavior of Bayesian-based functional mapping and demonstrated its applicability and utilization.The simulation design is based on an F 2 population containing 450 individuals.A linkage group of length 100 cM was simulated with 11 equallyspaced, ordered codominant markers.A QTL that affects growth curves was assumed at 34 cM from the first marker.The true values of curve parameters for three QTL genotypes are given as: Θ 0 = (33.4,11.2, 0.65), Θ 2 = (36.7,11.9, 0.65), Θ 1 = (35.6,11.2, 0.64), which are close to the estimates for body mass growth from the real mouse example used above.A total of 10 evenly spaced time points were assumed to measure growth curves.The residual covariance matrix among different time points was set to be the same as the sample covariance matrix estimated from the mouse example.Given this hypothesized covariance matrix and curves parameters for the three QTL genotypes, the heritability of the simulated growth trait at a middle time point is about 0.1.The time-dependent phenotypic values of the growth trait are assumed to follow a multivariate normal distribution.
The priors for curve parameters Θ j (j = 0, 1, 2) are assumed to be a multivariate normal distribution centered at (30, 10, 0.7) with covariance matrix Λ = diag{10, 5, 4}.We used a uniform distribution as a prior for the QTL location over the simulated linkage group.Three different approaches were used to estimate the residual covariance matrix Σ: (1) estimating the unstructured covariance matrix based on an inverse Whishart prior with T = 10 degrees of freedom, (2) estimating the structured covariance matrix based on the SAD(1) model by imposing IG(1, 1) as a prior for innovation variance ν 2 and N ormal(0, 10) as a prior for antedependence parameter φ, and (3) estimating the structured covariance matrix based on the AR(1) model by imposing IG(10, 1) as a prior for variance σ 2 and U nif orm(−1, 1) as a prior for correlation ρ.For each MCMC-implemented run, 10,000 initial "burn-in" iterations were discarded and, after then, samples are collected once for every 60 cycles so that a working set of 1,000 states were used to estimate the posteriors of parameters.The MCMC experiment was repeated several times with the same simulated data set to ensure the chain to converge to a stationary distribution.Table 2 gives the Bayesian estimates of the QTL location and growth curves for different QTL genotypes, along with the 95% empirical highest posterior density (HPD) confidence regions, under three different matrix-estimating approaches described above, respectively.It can be seen that our Bayesian-based functional mapping provides reasonably good estimates of the parameters in terms of accuracy and precision.The estimates of parameters, especially the QTL location (Figure 4), are affected by matrix-estimating approaches.Overall, the SAD(1)-structured approach is more precise in parameter estimation than the unstructured approach, whereas the unstructured approach is more precise than the AR(1)-structured approach.The SAD(1)-structured approach gives a narrowest confidence interval [33.01, 36.30] for the localization of QTL, followed by the unstructured [28.02, 38.16] and AR(1)-structured approach [25.54, 41.39].Also, the AR(1)-structured approach generates a small peak for the marginal posterior distribution at a wrong location (Figure 4), thus increasing a risk to detect a spurious QTL.The simulated data set was further analyzed by ML-based functional mapping in which the likelihood ratio (LR) test statistic was computed and plotted across the linkage group (Figure 5).We found two drawbacks for the ML-based approach.First, the LR profile has two small peaks, giving a chance to claim a false positive QTL.On the other hand, none of these peaks is sharp over a flatted profile, suggesting less power and precision for QTL detection by ML-based functional mapping.Second, the significance test for the ML-based approach is based on computationally expensive permutation tests.ML-based functional mapping costs about 20 times in computing times (including 100 permutation tests) as much as does Bayesian-based functional mapping.

Discussion
The genetic architecture of complex traits can well be understood by incorporating their underlying developmental features described by mathematical functions.Functional mapping that integrates genetics, statistics and developmental biology as a whole can be useful for deciphering the ontogenetic development of the genetic control of a complex quantitative trait [8,9].Original models for functional mapping were derived within the maximum likelihood (ML) context and implemented with the EM algorithm.A similar approach for mapping logistic QTLs was used by Malosetti [8].Although ML-based approaches possess many favorable statistical properties in parameter estimation, they may not be powerful enough to handle the complexity of high-dimensional QTL mapping models, as often seen in functional mapping.As an increasingly popular approach, Bayesian methods display remarkable capacity to estimate genetic parameters in QTL mapping [17-19, 39, 40].The profile of the log-likelihood ratio (LR) test statistics between the full (there is a QTL) and reduced (there is no QTL) models across a simulated linkage group.The covariance matrix was structured by the SAD(1) model.
In this article, we derived a general Bayesian framework for functional QTL mapping of dynamic traits and implemented Markov chain Monte Carlo (MCMC) algorithms to locate genomic positions of QTLs, and estimate the mathematical parameters that define a biological process and the statistical parameters that model the covariance structure.The Bayesian-based model allows the estimation of these parameters and their confidence intervals based on posterior distributions, and has great power to handle complex estimation issues related to functional mapping in an effective way.Like original parametric functional mapping [9], the new model allows to approximate the ontogenetic changes of the genetic effects triggered by a QTL.Because many biological processes, such as growth, follow a particular pattern of development, the ontogenetic control of a QTL can be mathematically described and, thereby, tested by estimating the parameters that define a biological process.The new model also take an advantage of functional mapping to model the structure of covariance matrix by a stationary or nonstationary approach.
The application of Bayesian approaches to functional mapping was also considered by other authors.Yang and Xu [19] integrated Bayesian shrinkage approaches to map dynamic QTL, but their model was based on a nonparametric Legendre polynomial fitting.Although this treatment may be statistically flexible, its biological relevance may be limited because no biologically sensible functions are deployed.Also, Yang and Xu's [19] approach did not utilize the high-efficiency characteristic of functional mapping through modeling the structure of an autocorrelated covariance matrix.Given a simulated time-dependent covariance matrix, Bayesian-based functional mapping obtains different results about the precision of parameter estimation, depending on whether and how such a covariance matrix is structured.This suggests the importance of estimating the residual covariance matrix effectively and efficiently.The unstructured approach captures the full information of time-dependent variances and covariances, which is effective but not efficient.The SAD(1)-structured model approximates the changes of variance and correlation over time [20], which is not only efficient (through estimating fewer parameters) but also effective (by reflecting the reality of this simulated data set to great extent).The AR(1)-structure model assuming the stationarity in both variance and correlation is efficient but not effective as much as the SAD(1) model.Overall, the SAD(1)-structured approach performs best for this particular simulated longitudinal data, followed by the unstructured and AR(1)-structured approach in an order.If the covariance does not have a structure, an unstructured approach, such as one based on the Wishart prior, should be used.In Appendix D, we describe a different approach for estimating the full matrix based on a reference prior.
The model was used to reanalyze a published data set on the growth of body mass in mice, confirming the discovery of a few QTL detected by conventional ML-based functional mapping.Yet, the new model provides estimates of confidence intervals of curve parameters, thus allowing better statistical inferences about the genetic control of dynamic QTLs.Simulation studies show that Bayesian-based functional mapping is more robust than ML-based functional mapping, in that the former provides more reasonable estimates of QTL positions and dynamic QTL effects when the heritability of growth curves is modest (0.1) compared to the latter.However, a unique issue related to the Bayesian approach is about its stability or sensitivity in parameter estimation to different priors.One approach for testing the model's stability is to initiate the chain with several different starting points and/or different priors under the same MCMC sampling scheme [41], and further examine how the estimates depend on the choices of initial values and priors.
The model proposed in this article can be modified by considering a network of genetic control.As a basic Bayesian framework, our single-QTL interval mapping model is not adequate to explore the effects of interactions on developmental variation for a growth trait between different QTLs from the same genome [43] or different genomes [44], and between QTLs and environments [24].One of the significant advantages of Bayesian approaches lies in the estimation of the optimal number of QTLs involved.Variable selection via stepwise regression is used in ML mapping, but it is highly computationally expensive.Corresponding to this variable selection procedure, reversible jump MCMC is proposed in Bayesian analysis [45,46], although it is subject to poor mixing and a slow convergence to the stationary distribution [47,48].More efficient methods based on Bayesian shrinkage analysis [16] and stochastic search variable selection [17] have now been proposed.These methods do not rely upon any explicit form of variable selection; rather they proceed implicitly by shrinking the effects of excessive QTLs to zero.Our Bayesian-based model modified by considering complicated features of growth and development will certainly prove its value in elucidating the genetic architecture of dynamic traits and will probably be the beginning of detecting the driving forces behind developmental genetics and its relationship to the organism as a whole.
Appendix A: Estimating (λ, Q, Θ, Σ) The Markov chain for Bayesian functional mapping described in the Parameter Estimation and Algorithm section is constructed as follows: Step 1. Initialize the iteration at an arbitrary point (λ 0 , Q 0 , Θ 0 , Σ 0 ), which has a positive posterior density; Step 2. Modify four blocks of the unknowns parameters and move to a new state from the previous step and Σ k .More specifically, given the values of the unknowns (λ, Q, Θ, Σ) from the current state, we proceed as follows: U pdating λ : In each step, following the idea of Satagopan et al. [14], λ is updated by using the Metropolis algorithm.A new value of λ * is generated from Uniform(max(0, λ − δ), min(λ + δ, D m )) (where δ > 0 is the tuning parameter), and this proposed distribution is denoted by q(λ, λ * ).This proposed λ is accepted with probability min(α λ , 1), and the state keeps the current value if the proposal is rejected, where α λ is given as: Note that, π(λ * |y, Q, Θ, Σ) = π(λ * |y, Q, Θ, Σ, M) Similarly, we have: Hence, the acceptance probability (A 1) can be simplified as: U pdating Q : Because of independence among n progeny, Q is updated by separately updating each Q i .For each progeny i and QTL genotype j, the full conditional density is in the form of a multinormial with cell probabilities: Hence, at each cycle, we can sample the QTL genotype Q i directly from this full conditional density.U pdating Θ : We update each Θ j successively by a Metropolis-Hastings algorithm.For each QTL genotype, a new value Θ * j is generated from a proposed density q(Θ j , Θ * j ), given the current Θ.Evaluate the acceptance probability of the move is min(1, α Θ j ).In general, α Θ j can be expressed as: Note that the choice of the Metropolis kernel q is essentially arbitrary, and a symmetric q in the sense that q(Θ j , Θ * j ) = q(Θ j * , Θ j ) is usually preferred.And in that case, the ratio q(Θ * j , Θ j )/q(Θ j , Θ * j ) is canceled in the above expression (A 6).Here, for the proposed density, we use a multivariate normal distribution centered at the current Θ, with variance-covariance matrix given by an information-type matrix [49] whose inverse has (u, v)th element: 7) This expression combines the expected information (the first term) with the prior information (the second term) and offers an advantage of avoiding singular information matrices.Unfortunately, a tedious initial analysis has to be conducted to obtain estimates of Θ j and Σ from which to evaluate (A 7).So, the Metropolis algorithm described above can be carried out by using an arbitrary variance-covariance matrix.The posterior means of Θ j and Σ are then plugged into (A 7) from the subsequent analysis.U pdating Σ : We generate a new value of Σ −1 directly from its full conditional posterior distribution.This is straightforward since it has an explicit expression for its full conditional posterior distribution.Comparing the reference prior with the Jeffrey's prior, it is noted that the posterior given in equation (D3) is always proper.Also note that, since this reference prior put more mass near the region for equal eigenvalues, it can produce an estimator with better eigenstructure.Again, it is very difficult to analytically evaluate the posterior (D3).Yang and Berger [53] suggested using a Metroplisized hit-and-run sampler algorithm to obtain the integration.The detail sampling procedure at the kth iteration is given as follows: (1) Given the current positive-definitive matrix Σ k , we set W k = logΣ k , in the sense that (2) Randomly generate a symmetric p by p matrix T, with elements t ij = z ij /( l m z 2 lm ) 1/2 , where z ij i.i.d.N (0, 1), f or i j; (3) Set W * = W k + νT where ν is generated from N (0, 1); (4) Update W k with an acceptance probability min(1, π * (W * |y, Θ, Q, λ)/π * (Σ k |y, Θ, Q, λ)).

Figure 1 .
Figure 1.A profile of Estimated marginal posterior distribution of the QTL location by assuming that exactly one QTL is located on one of the chromosome respectively.

Figure 2 .
Figure 2. Fitted growth curves for the three QTL genotypes (qq, red; Qq, blue; QQ, black) assuming a single QTL is located on mouse chromosome 6, 7 and 10.

Figure 3 .
Figure 3. Dynamic changes of the additive and dominant effect due to the QTL located on mouse chromosome 6,7,and 10 respectively.

Figure 5 .
Figure 5.The profile of the log-likelihood ratio (LR) test statistics between the full (there is a QTL) and reduced (there is no QTL) models across a simulated linkage group.The covariance matrix was structured by the SAD(1) model.

Table 1 .
Bayesian estimates of QTL locations and genotype-specific growth curves for the QTLs detected on mouse chromosome 6, 7 and 10.Numbers in parentheses are the 95% equal-tail confidence intervals.

Table 2 .
Bayesian estimates of QTL locations and genotype-specific growth curves for an assumed QTL from the simulated data set for an F 2 population of 450 individuals based on different covariance-structuring approaches.Numbers in parentheses are the 95% equal-tail confidence intervals.