BACKGROUND
computer simulation is a resource which can be employed to identify optimal breeding strategies to effectively and efficiently achieve specific goals in developing improved cultivars  <cit> .
once identified, optimal strategies can be incorporated in the development ‘process’ to facilitate maximal genetic gains and accelerate the breeding process  <cit> .
through computer simulation, innovative approaches heretofore not feasible without use of current genomic-based technologies can be specifically tailored to meet the need at hand.
in some instances, it is crucial to assess in silico, the options as well as the impact of various crossing schemes and breeding strategies on performance for certain key agronomic traits such as grain yield.
such is the case with introgression  of value-added traits by means of backcross breeding, the goal of which is the recovery of all the performance attributes of the elite variety or hybrid targeted for conversion along with the new genetic elements  associated with the value-added trait.
for this, a means by which the gene effects for the performance attributes of the line to be converted can be modeled in the simulation is particularly advantageous, yet it requires an accurate depiction of the distribution of gene effects.
typically, recovery of the target line  is estimated by the average proportion of genetic material carried through the backcrossing process and this estimate implicitly assumes that the many genes for key quantitative traits like grain yield are dispersed uniformly across the genome, each contributing only small effect.
by including a more realistic representation of gene effects in the genome model to assess backcross breeding strategies, the means to most rapidly and effectively recover not only the germplasm per se, but the important genes contributing to performance of the variety or hybrid targeted for conversion, can be considered in evaluating strategies and approaches.
bennewitz and meuwissen  <cit>  explored the distribution of additive and dominance effects of identified qtl  from three f <dig> populations of pigs evaluated for  <dig> meat quality and carcass traits, recognizing the value of modeling these effects, some of which are large, over use of the infinitesimal model which assumes an infinite number of qtl each with small effect.
capitalizing on the large number of qtl studies, bennewitz and meuwissen  <cit>  conducted a meta-analysis of published qtl mapping data across traits to infer the distribution of additive qtl effects as well as dominance coefficients, fitting a gaussian mixture model .
the idea of utilizing gmm is based on the notion that various qtl and associated genes fall into a number of classes of different-sized effects.
the merit of employing gmm is its flexibility with different combinations of mixtures of normals leading to different shapes of the distribution.
in the finite mixture model, the number of components k must be pre-specified.
the value could be determined based on some specific information or criteria, such as the akaike information criterion or the bayesian information criterion.
this requirement, frequently encountered in parametric statistics, can be sidestepped by introducing a nonparametric dirichlet process which assumes an infinite number of components.
the dirichlet process is defined as a random process by which a sample drawn is a discrete distribution; it can be considered a ‘distribution over distributions’ and has been used widely in the field of population genetics to explain population structure  <cit> .
also desiring to capitalize on the large number of qtl studies, we took the meta-analysis concept a step further.
we devised a method to represent the genomic distribution of additive and dominance gene effects associated with quantitative traits, which utilizes a modified dirichlet process gaussian mixture model   <cit>  to fit the number of mixture components and estimate parameters.
as a departure from traditional dpgmm which only models qtl effects, we modified the model to be able to accommodate both qtl effects and their respective variances.
utilizing previously-identified qtl for a number of quantitative traits in maize, the modified dirichlet process implements a chinese restaurant process  to assign component  membership and uses gibbs sampling to update conditional posterior distributions.
our purpose in devising this method was to facilitate representation of the genetic architecture of grain yield and other key traits for use in computer simulations to optimize breeding strategies for multiple trait introgression .
trait introgression involves backcross breeding and, therefore, bi-parental populations with no more than two potential alleles for a given locus; this was the primary scenario we intended to model.
however, we also explored whether the method would be pertinent to other modeling activities that may involve estimation of breeding value  or prediction of performance based on priors  <cit> , scenarios which could involve multi-allelic populations in various crop or animal species.
besides use in computer simulation, an accurate depiction of the qtl effects could contribute to a better understanding of the overall genetic architecture contributing to variation of expression of a particular trait of interest  <cit> .
methods
description of dpgmm and priors
to begin, we modeled the distribution of additive qtl effects and dominance coefficients using mixtures of normal distributions, namely gmm  <cit>  .
the goal was to assign genetic effects to different mixture components.
two latent variables were introduced: 1) the total number of mixture components  and 2) the assignment of ith qtl effect to components .
the gmm model was modified to accommodate the standard errors of qtl effects:  <dig> pyi|λ <dig> …,λk~∑k=1kπknyi;μk,σk2+τi <dig>  where yi is the ith observed qtl effect, τi is the known standard error of ith effect which is calculated from the qtl analysis, and λk = {πk, μk, σk2} is the kth parameter set, where variables πk, μk and σk <dig> are the mixing proportion, mean and variance of the kth mixture component, respectively.
the dpgmm can be formulated hierarchically as follows  <cit> :  <dig> pyi|ci,Λ~nyi;μci,σci2+τi2ci|π1:k~discreteπ <dig> π <dig> …πkμkσk2~g0π <dig> π <dig> …πk~dirichletα/k,…,α/kΛ=λ <dig> λ <dig> …λk where g <dig> is a joint prior distribution for  and mixing proportions π1 : k are drawn from a symmetric dirichlet distribution with concentration parameter α. conditional on the mixing proportions, the latent indicator variables ci’s were sampled from a discrete distribution, specifically a multinomial distribution, and the prior for the ci in model  <cit>  could be written as a probability conditional on c− i  <cit> : pci=k|c−i,α=n−i,k+α/kn−1+α,
where n− i,k is the number of effects, not including yi that are linked with class k. and as k goes to infinity , the limits of the prior for the ci reach the following:  <dig> pci|c−i,α=n−i,kn−1+αci=k,n−i,k>0αn−1+α∀i≠i',ci≠ci'.
where i’ is the complement of the set i. as k → ∞, the dirichlet distribution becomes a dirichlet process in the limit  <cit> .
thus, infinite limit of model  <cit>  can be written as a dpgmm:  <dig> yi~nyi|θi,σi2μiσi2~gg~dpαg <dig>
where θi ~ n is a nuisance parameter, g is a random discrete distribution drawn from dp, and g <dig> was the base distribution, which specified the joint prior distribution of .
given that the regular choice of priors for the mean and variance of the gaussian are normal and inverse gamma distributions, respectively, conjugate joint priors n* ig were chosen in the model in order to implement the following gibbs sampling.
gibbs sampling
in bayesian framework, unknown variables were sampled and updated from the conditional posterior distribution using markov chain monte carlo   <cit> .
considering the likelihood and priors in formulae  <dig> and  <dig>  the full joint posterior distribution can be written as follows:  <dig> pc,μ,σ2|y∝∏i=1nnyiθiσi2πθiμiσi2ci=∏i=1nnyiθiσi2nθiμiτi2nμiμ0σ02igσi2r1r2pci|c−i,α.
unobservables  were repeatedly sampled and updated from their posteriors, conditional on all other variables.
the gibbs sampler was implemented as follows:initialization: assign initial values for  where k =  <dig> and ci =  <dig>  for i = 1 : n.
update θi: the conditional posterior distribution of θi was
pθi|else∝nyiθiσk2nθiμkτi2∝nθiyiσk2+μkτi21σk2+1τi211σk2+1τi <dig> 3) update cluster indicators ci : the conditional posterior probabilities for ci were:
pci=k|else∝nyiθiσk2nθiμkτi2pci|c−i,α=n−i,k2πσk2τi2exp−yi−θi22σk2−θi−μk22τi <dig>
pci=k+1|else∝α∬nyiθiσk+12nθiμk+1τi2nμk+1μ0σ02igσk+12r1r2dμk+1dσk+12∝α∫nyiθiσk+12igσk+12r1r2dσk+12∫nθiμk+1τi2nμk+1μ0σ02dμk+1=α2πr2r1Γr1Γr1+1212yi−θi2+r2r1+121τi2+σ02exp−θi−μ022τi2+σ <dig> where Γ is the gamma function.
note that constant 1n−1+α was omitted in both probabilities and  were unknown and needed to be integrated out to leave ci as the only variable to be estimated from the markov chain.
the dirichlet process was represented via the crp  <cit> .
effects were assigned to either currently holding cluster or a new cluster based on the above probabilities.
if a new cluster was chosen, then the cluster size was increased, i.e.
k + 1 → k. in case of n− i,k =  <dig>  the kth cluster was eliminated and the cluster indicators were decreased by one, i.e.
k → k −  <dig> ) resample and update  suggested by formulae  <dig> as per neal  <cit>  as follows:
pμk|θi∈kthcluster,else∝∏i=1nknθiμkτi2nμkμ0σ02~nμk∑i=1nkθiτi2+μ0σ02∑i=1nk1τi2+1σ021∑i=1nk1τi2+1σ <dig>  pσk2|yi∈kthcluster,else∝∏i=1nknyiθiσk2igσk2r1r2~igσk2;r1+nk <dig> ∑i=1nkyi−θi2+r <dig>
where nk is the number of effects associated with the kth mixture component.
the derivations of the fully conditional posterior distributions are detailed in the appendix.5) repeat steps  <dig> to  <dig>
gibbs sampler was implemented with  <dig>  iterations of the mcmc to update conditional posterior distributions.
the first  <dig>  samples were discarded as burn-in and the rest of the  <dig>  samples were used to construct joint posterior distribution.
the hyper-parameters in algorithm  <dig> were set to be α =  <dig> , r1 =  <dig>  r2 =  <dig> , μ0 =  <dig>  σ02 =  <dig> .
among hyperparameters, alpha was empirically set to  <dig>  based on the simulation results.
convergence was checked by inspection of negative log-likelihood plots.
after the burn-in period, when the mcmc converges to the stationary distribution, sampled parameters were collected to form the posterior distribution.
we employed posterior means for estimating the mean and variance μ^kσ^k <dig> and posterior modes for estimating ĉi, which was further used to infer π^k.
the bayesian confidence interval , which is the counterpart of the confidence interval in frequentist statistics, was defined as posterior probability that the parameter lies within the interval: ∫−∞apΛ|ydΛ=∫b∞pΛ|ydΛ=α/ <dig>
where α is the significance level.
instead of analytically estimating the confidence interval, the confidence interval for μ^kσ^k <dig> was numerically estimated from quartiles of posterior distribution.
demonstration of method performance with simulated data
to demonstrate the performance of the proposed method, two simulated data sets were processed.
simulation i facilitated evaluation of model performance given complete data.
it was generated from three gmms with respective means - <dig>   <dig> and  <dig> and variances of  <dig> ,  <dig>  and  <dig> , respectively.
a total of  <dig> simulated qtl effects were uniformly distributed  to the three components.
this data set represents the case wherein all true qtl are known.
simulation ii facilitated evaluation of model performance on a truncated distribution.
truncation points were arbitrarily set to ±  <dig> .
the incomplete data set was intended to represent the common situation with qtl mapping data wherein all genetic effects are not detected, especially those with effects of near-zero magnitude.
a truncated gaussian mixture with two mixture components was simulated.
zero mean was assigned to both components.
the first mixture component had mixing proportion π <dig> and variance σ <dig> of  <dig>  and  <dig> , respectively; the second mixture component had π <dig> and σ <dig> of  <dig>  and  <dig> , respectively.
in both simulations, the standard error  τi was generated from a uniform distribution .
implementation with real data
the model was also applied using real data to derive the distributions of additive effects and dominance coefficients.
additive qtl effects were assembled from previous qtl mapping studies performed in corn  .
messmer et al.
<cit>  had evaluated recombinant inbred lines derived from a cross between two subtropical white dent maize lines to map genes controlling yield components and other traits  <cit> , identifying  <dig> qtl in total which are included in data i. briggs et al.
<cit>  had utilized a maize-teosinte backcross  population to explore genes controlling domestication and morphological traits such as plant architecture, primary tassel and lateral inflorescence, identifying  <dig> qtl in total which are included in data ii.
data iii is derived from five maize qtl mapping studies involving segregating populations, all of which share a common parent b <dig>  comprising a total of  <dig> quantitative trait loci including a)  <dig> qtl for kernel oil concentration mapped in an f <dig> population  <cit> ; b)  <dig> qtl for root angle and plant height mapped in an f <dig> population  <cit> ; c)  <dig> qtl for stalk digestibility and kernel composition mapped in a f <dig> population  <cit> ; d)  <dig> qtl for stripe disease resistance mapped in an f <dig> population  <cit> ; and e)  <dig> qtl for grain yield and yield components mapped in an f <dig> population under water-limited conditions  <cit> .
all qtl mapping studies employed either composite interval mapping or multiple interval mapping to detect qtl.
<cit> .
furthermore, these qtl studies all reflect estimates of gene effects in bi-parental populations, which fit with the backcross breeding scenario we intended to simulate.
the histograms of observed additive effects, with values represented in units of phenotypic standard deviations, are shown in fig.
<dig> for all three data sets.table  <dig> qtl associated with various traits across four data sets.
data i, ii, iii, and iv were included in the analysis of qtl additive effects and data iii was used in the analysis of qtl dominance coefficients
to explore application with multi-allelic, large-scale data sets as an alternative to bi-parental qtl data sets, additive qtl effects were assembled from high-resolution genome-wide association studies  with the maize nam  population  <cit> , data provided courtesy of jason g. wallace .
from the  <dig> traits in the data set,  <dig> were chosen to generate data iv: ‘20-kernel weight’  and ‘days to anthesis’  .
only those significant single nucleotide polymorphisms  featured in at least three resample inclusions were utilized to fit the distribution .
the qtl had been detected through composite interval mapping or multiple interval mapping and the association mapping had been performed using the forward-regression genome-wide association method in tassel  <dig> .
<dig>  <cit> .
results were compared to those obtained with qtl mapping data sets from bi-parental populations to contrast multi-allelic vs. bi-allelic, number of qtl, power in detection of the qtl, and single-trait distributions vs. distributions representing multiple traits.
in addition, the distribution of dominance effects was explored.
dominance coefficients, which are defined as the ratio between the observed qtl dominance deviation and absolute value of qtl additive effects, were assembled from data iii.
the absolute value of additive qtl effects was used because the sign of qtl effect only signifies which parent had contributed the favorable allele, not the true direction of specific additive effect.
the se for additive qtl effects and dominance coefficients was measured to take into account the experimental error.
if logarithm  of the odds  scores for qtl were absent, standard errors were generated by taking sample standard deviation of effects from multiple experiments.
the se of data ii and of data from dintinger et al.
<cit>  incorporated in data iii were produced in this way, where only those qtl detected in at least two environments were included in the final dataset.
for the rest of the studies, ses were derived from lod scores as suggested by hayes and goddard  <cit>  or, in the case of the gwas data sets, the ses were determined from the sample errors of the discovered snp effects  <cit> .
standard errors of dominance coefficients were estimated by the delta method suggested by bennewitz and meuwissen  <cit> , assuming no covariance between additive and dominance effects.
specifically, sed/a=d/a*seaa2+sedd <dig>  where sea and sed were standard errors for additive effects a and dominance effects d, respectively.
additive qtl effects from qtl mapping studies were scaled by their corresponding phenotypic standard deviations in order to combine data across traits.
the phenotypic standard deviations were computed using raw data if available from the qtl study.
otherwise, the error variance and heritability of the trait were used to calculate the phenotypic standard deviation or, absent this, the range of phenotype values were used.
phenotypic range was assumed to be  <dig> standard deviations, which covers almost 100 % of the values, considering that most traits follow a normal distribution.
we did not apply phenotypic standardization to the gwas data since an abundance of data points were available.
since for the data sets from which dominance effects were generated, none of the above three conditions was fulfilled to obtain the phenotypic standard deviations for additive qtl effects, the additive effects from those five corn studies were not utilized in analyzing the distribution of additive qtl effects.
note that the scaling process was not necessary for dominance coefficients, because the phenotypic standard deviation canceled out in the d/a ratio.
due to limited statistical power of qtl mapping studies  <cit> , many qtl with near-zero effects were likely not detected in the published studies used in this analysis, which is effectively analogous to a truncation of the additive qtl effect distribution near zero.
faced with this issue, bennewitz and meuwissen  <cit>  suggested a “doubling” process, given the assumption that the true qtl effects occur at the highest frequency around zero.
to compensate for the ‘missing’ qtl effects, a doubling of the data was done to ensure that the mean of each mixture component is estimated at zero, preserving the characteristic greatest density at zero for each cluster.
with doubling, both signs for the same qtl additive effect were created.
for example, for ith effect yi with se τi, − yi was added to the data with the same se.
the above procedure leverages the fact that absolute values of alternative homozygous genotypes at a qtl are the same by definition e.g.
in bi-parental populations  <cit> .
the “doubling” process was not applied to dominance coefficients, because most loci have observed effects around zero.
RESULTS
in this study, we employed a new method, namely dpgmm, to describe the distribution of qtl additive effects and dominance coefficients in the form of mixtures of normals.
although similar to the fitting of a mixture of normals using a modified expectation-maximization  algorithm, this approach differs primarily in the way of dealing with cluster size .
with dpgmm, the number of mixture components does not need to be specified .
in contrast, with use of a finite mixture model, the number of components needs to be preset and later decided under certain circumstances, or determined by some measure, e.g.
akaike information criterion or bayesian information criterion.
the optimum cluster size  will strike a balance between maximum data compression  and maximum accuracy .
in dpgmm, the dirichlet process is represented by the crp, wherein a data point is assigned either to a currently occupied mixture component with probability proportional to the number of data already held in that cluster or to a new cluster with probability proportional to the concentration parameter, g <dig>  by the same token, with each iteration of gibbs sampling, the cluster indicators are also updated along with parameters like the mean and variance.
as such, dpgmm fits the data distribution and explores the potential number of mixture components simultaneously.
to demonstrate the performance of dpgmm in fitting the distribution of qtl effects and to verify accuracy in estimating associated parameters of this distribution, two simulated data sets were processed.
simulation i was structured to represent the case wherein all true qtl are known i.e.
complete data.
simulation i, which featured three components, resulted in a histogram of genetic effects from which it is difficult to infer the number of mixture components visually .
in simulation i, dpgmm clearly fitted the data to three clusters with estimated values close to true values for parameters involving the proportion of mixing among clusters, the mean and the variance of each mixture component .
dpgmm predicted accurately the mean and variance of clusters  <dig> and  <dig>  although missed assignments of cluster membership were observed.
in contrast, the mixing proportion of cluster  <dig> was estimated precisely; however, certain deviations from the true mean and variance were observed.fig.
<dig> histograms of simulated effects from gaussian mixtures: a  three components having mean of - <dig>   <dig> and  <dig>  and variance of  <dig> ,  <dig>  and  <dig> , respectively, and equal mixing proportions for all three components; and   two components having zero means and variance of  <dig>  and  <dig> , respectively, and mixing proportions of  <dig>  and  <dig> , respectively.
distribution in  is truncated at points - <dig>  and  <dig>
k is the mixing proportion in the k
th cluster, and μ
k and σ
k <dig> are the mean and variance of k
th mixture component, respectively.
values expressed in units of phenotypic standard deviation
π
k
π^k
μ
k
μ^k
σ
k2
σ^k2
simulation ii data, based on a truncated mixture of normals featuring two mixture components with zero mean, resemble a scenario common to qtl mapping wherein near-zero genetic effects were not detected.
simulation ii produced a histogram with a pronounced gap around zero as expected in keeping with the data truncation .
in simulation ii, parameters were estimated with accuracy, except for the variance of cluster  <dig>  which was estimated at  <dig>  versus the true value of  <dig>  .table  <dig> true versus estimated  parameters in simulation ii.
π
k is the mixing proportion in the k
th cluster, and μ
k and σ
k <dig> are the mean and variance of k
th mixture component, respectively.
values expressed in units of phenotypic standard deviation
π
k
π^k
σ
k2
σ^k2
given a complete set of data, dpgmm could clearly assign membership to respective clusters with small prediction errors .
in the case of truncated data, dpgmm was still effective in predicting the correct number of mixture components and estimating the variance of components with greater variability; however, dpgmm was less effective in estimating variance of components with small differences among cluster members.
as shown in table  <dig>  the deviation of estimated variance  from true value  was somewhat large and might be attributed to the loss of small-value data in the sample.
this result is in accord with the conclusion of bennewitz and meuwissen that small effects could be missed easily with a mixture model  <cit> .
proceeding to real data, dpgmm was used as a method to fit a mixture of normals for which the number of components is not known, fitting additive effects and dominance coefficients based on previously published qtl mapping data for a number of quantitative traits in maize.
in addition, dpgmm was used to fit additive effects estimated from high-resolution gwas of the maize nam population to compare distributions produced with large-scale, multi-allelic data sets involving a single trait with those obtained with meta-analysis of bi-parental, lower-resolution studies involving multiple traits.
the fitted distributions are the outcome of applying the dpgmm method and these distributions could then become the basis for modelling the genetic architecture of quantitative traits in maize for computer-simulated explorations to identify optimal breeding strategies.
distributions of qtl additive effects
histograms of observed qtl additive effects expressed in units of phenotypic standard deviation were generated from the qtl mapping studies .
we noted that the histogram of data ii  resembled that observed with simulation ii , wherein near-zero effects were not included.
with data i, ii, and iii, it is difficult to infer the number of mixture components for additive effects visually from the histograms.
however, the number of mixture components was inferred by the mode of posterior distribution with regard to cluster indicator ci.
frequency tables of cluster membership clearly suggested fitting all data to one cluster for data i and data ii ; however, for data iii, two components are suggested to fit additive effects.
this likely reflects the types of traits included in each meta-analysis.
data i contains primarily grain yield and yield component traits as well as some developmental traits.
data ii is comprised of mainly domestication and morphological traits.
data iii includes yield and yield component traits under drought stress as well as grain quality traits.
results suggest more similarity among the traits within data i and within data ii, in contrast to data iii that was fitted with two clusters.fig.
<dig> histograms of observed qtl additive effects : a data i; b data ii; and c data iii
fig.
<dig> histograms of the cluster number: a data i; b data ii; and c data iii
for each of the three qtl mapping data sets, the fitted distribution was overlaid on the histogram of “doubled” data .
the fitted distributions of additive effects all have zero mean; variances differed .
the range of the observed qtl effects was tight around the mean for data i, whereas data ii showed some larger effects with absolute value nearing  <dig>  with data iii, the range of observed effects was greater still, with the largest effects nearing  <dig>  in absolute value; nonetheless, the variance of data iii additive effects was similar to that of the distribution of data ii.
the variance was smallest for data i comprised of mainly yield and yield component traits.fig.
<dig> fitted normal distributions to qtl additive effects : a data i; b data ii; c data iii; d data iv featuring traits of 20-kernel weight and days to anthesis
th cluster , the cluster mean and variance  and bayesian confidence interval  for parameters in the distribution of additive effects and dominance coefficients.
values expressed in units of phenotypic standard deviation
σ^12
σ^12
π^1
π^2
σ^12
σ^22
μ^1
σ^12
σ^12
σ^12
comparing the results obtained with the data iv gwas single-trait data sets, the data were fitted to one cluster for both 20-kernel weight and days to anthesis.
each fitted distribution was overlaid on the histogram of “doubled” data .
each distribution had an estimated mean of zero.
the variance associated with 20-kernel weight was estimated as  <dig> , with a narrow confidence interval,  <dig>  to  <dig> ; the variance for days to anthesis was estimated as  <dig> , with a narrow confidence interval,  <dig>  to  <dig>  .
the fitted distribution of 20-kernel weight based on  <dig> qtl shows a high degree of similarity to that for data i which is comprised of  <dig> qtl for yield and yield component traits, maturity, and abiotic stress .
the fitted distribution of days to anthesis based on  <dig> qtl shows a high degree of similarity to that for data ii which is comprised of  <dig> qtl for domestication and morphological traits .
compared to distributions from meta-analyses, the gwas data sets resulted in greater precision as displayed in the relatively smaller bci’s for the estimates of variance; the ratio of the size of the bci’s for data iv 20-kernel weight and data i is  <dig> , and for data iv days to anthesis and data ii is  <dig> .
thus, the accuracy using meta-analysis of multiple-trait qtl and qtl identified through gwas appears to be similar when the type of traits measured are similar.
however, precision is better with the latter, which likely reflects the significantly larger number of qtl in the data iv sets.
despite the significantly greater resolution of the gwas data sets, both distributions produced from data iv displayed a gap around zero: within ± <dig>  for 20-kernel weight and within ± <dig>  for days to anthesis.
bennewitz and meuwissen  <cit>  discussed the potential drawbacks of using meta-analyses of qtl effects detected in bi-parental populations to characterize distributions of gene effects.
a primary concern was for failure to identify all true qtl within the identified qtl sets either because certain alleles were excluded from the data set or because of a lack of statistical power to detect smaller-effect qtl  <cit> .
even in cases where marker density and genome resolution is high, all qtl for a trait of interest may not be identified if mapping methods are utilized that screen detected effects against a significance threshold which basically excludes most of the small effect qtl.
furthermore, in genomic selection, the particular statistical method used to ‘train’ the model can influence the distribution of effects.
despite the higher resolution and multi-allelic nature of the data iv data sets, the inability to detect qtl of near-zero effect was apparent, yet not problematic in fitting a distribution.
clearly, the distributions are centered at zero with highest frequencies observed in the zero vicinity and recognized as single-component.
thus, the results obtained with the gwas data sets demonstrate the robustness of the dpgmm method with use of either meta-analysis of qtl identified across traits or large-scale data sets comprised of qtl for a single trait.
if enough data are available, e.g.
gwas dataset, distributions of qtl effects based on a single trait rather than across traits might be more appropriate and useful.
and with the advent of genotype-by-sequence  and other genotyping technologies that facilitate high-resolution marker sets, very dense marker sets are more available and more widely utilized in qtl identification  <cit> .
for our purpose in exploring the impact of backcrossing strategies on recovery of performance of the hybrid targeted for conversion  <cit> , the meta-analysis involving the data i set based on a bi-parental population met our objectives well in that it focused on yield and yield component traits and additionally included other key traits essential to performance recovery i.e.
maturity and abiotic stress tolerance.
distribution of qtl dominance coefficients
observed dominance coefficients obtained from meta-analysis of five mapping studies varied in magnitude from less than - <dig>  to more than  <dig>  , suggesting that all classes of dominance were represented among the traits measured.
around 50 % of the qtl  displayed an d/a ratio in the range of - <dig>  to  <dig> , indicating partial recessivity, additivity, and partial dominance gene action.
approximately 25 % of the qtl exhibited either partially dominant or dominant gene action  or partially recessive or recessive gene action .
furthermore, 25 % of the qtl exhibited apparent overdominance  or underdominance  gene action.fig.
<dig> histogram of observed dominance coefficients from meta-analysis based on five mapping populations
dominance coefficients were fitted with the normal distribution using dpgmm.
the mode of posterior distribution with regard to k is  <dig>  suggesting that all data could be fitted to a single component ; estimates of the distribution mean and variance were provided by mcmc as well.
figure  <dig> displays the estimated distribution overlaid on the density plot of observed data.
the estimated mean of the distribution was  <dig>  with 95 % bci of  <dig>  to  <dig> , and variance of  <dig>  with 95 % bci of  <dig>  to  <dig>  .
the result that the distribution of dominance coefficients was fitted to a normal distribution with a positive mean conformed to previous studies  <cit> .fig.
<dig> estimation of cluster number, mean, and variance of the fitted distribution of dominance coefficients through mcmc
fig.
<dig> normal distribution fitted to the dominance coefficients, with estimated mean at  <dig>  with 95 % bayesian confidence interval to be  <dig>  and  <dig>  and estimated variance at  <dig>  with 95 % bayesian confidence interval to be  <dig>  and  <dig>
CONCLUSIONS
the dpgmm method offers an alternative to the over-simplified infinitesimal model in computer simulation as a means to better represent the genetic architecture of quantitative traits, which likely involve some large effects in addition to many small effects.
furthermore, it confers an advantage over other methods in that the number of mixture model components need not be known a priori.
the dpgmm method takes advantage of prevalent qtl data to approximate the distributions of additive and dominance gene effects.
the dpgmm method is robust with use of either meta-analysis of smaller-scale qtl analyses involving a number of traits or large-scale, single-trait qtl data sets.
furthermore, qtl data sets from bi-allelic or multi-allelic populations can be utilized.
although the data sets may be missing some near-zero qtl effects that were not resolved in the fundamental analyses, the methodology is able to accommodate this drawback.
r code to facilitate use of the dpgmm method is included in a additional file  <dig> to this paper .
the distributions of qtl additive and dominance effects highlighted through this study were used in modelling the genetic architecture of grain yield and other key performance traits for computer-simulated explorations to identify optimal breeding strategies to facilitate introgression of multiple value-added traits into an elite maize hybrid.
maize grain yield is a complex trait involving dominant and over-dominant gene action.
other traits important to recovery of the essential performance attributes of the hybrid targeted for conversion include maturity, resistance to lodging, and abiotic stress tolerance and these have a bearing on yield performance as well.
readers are directed to the recent work of sun and mumm  <cit>  for an example of how the dpgmm-estimated genetic distribution parameters were deployed in computer simulations to evaluate breeding strategies.
appendix
derivations of the fully conditional posterior distributions are given as follows: pθi|else∝nyiθiσk2nθiμkτi2∝exp−12*yi−θi2σk2*−12*θi−μk2τi2∝exp−12*θi−yi2σk2+θi−μk2τi2∝nθiyiσk2+μkτi21σk2+1τi211σk2+1τi <dig>  pci=k|else∝nyiθiσk2nθiμkτi2pci|c−i,α∝12πσk2exp−yi−θi22σk2*12πτi2exp−θi−μk22τi2*n−i,k=n−i,k2πσk2τi2exp−yi−θi22σk2−θi−μk22τi <dig>   pci=k+1|else∝α∬nyiθiσk+12nθiμk+1τi2nμk+1μ0σ02igσk+12r1r2dμk+1dσk+12∝α∫nyiθiσk+12igσk+12r1r2dσk+12∫nθiμk+1τi2nμk+1μ0σ02dμk+ <dig>
the first integration part: ∫nyiθiσk+12igσk+12r1r2dσk+12=∫12πσk+12−12exp−yi−θi22σk+12*r2r1Γr1σk+12−r1−1exp−r2σk+12dσk+12=∫12πr2r1Γr1σk+12−r1+12−1exp−12yi−θi2+r2σk+12*dσk+12=12πr2r1Γr1Γr1+1212yi−θi2+r2r1+12∫igσk+12;r1+ <dig> yi−θi2+r2*dσk+12=12πr2r1Γr1Γr1+1212yi−θi2+r2r1+ <dig>
the second integration part: ∫nθiμk+1τi2nμk+1μ0σ02dμk+1=∫12πτi2exp−12θi−μk+12τi212πσ02exp−12μk+1−μ02σ02dμk+1=∫12πτi212πσ02exp−12θi−μk+12τi2−12μk+1−μ02σ02dμk+1=12πτi212πσ02∫2π11τi2+1σ02μk+1θiτi2+μ0σ021τi2+1σ0211τi2+1σ02exp−12θi−μ02τi2+σ02dμk+1=exp−12θi−μ02τi2+σ0212πτi212πσ022π11τi2+1σ02=12πτi2+σ02exp−θi−μ022τi2+σ <dig>
together we have: pci=k+1|else=α2πr2r1Γr1Γr1+1212yi−θi2+r2r1+121τi2+σ02exp−θi−μ022τi2+σ <dig>  pμk|θi∈kthcluster,else∝∏i=1nknθiμkτi2nμkμ0σ02∝exp−12∑i=1nkθi−μk2τi2−12μk−μ02σ02~nμk∑i=1nkθiτi2+μ0σ02∑i=1nk1τi2+1σ021∑i=1nk1τi2+1σ <dig>  pσk2|yi∈kthcluster,else=∏i=1nknyiθiσk2igσk2r1r2∝σk2−nk2exp−12∑i=1nkyi−θi2σk2σk2−r1−1exp−r2σk2=σk2−r1+nk2−1exp−12∑i=1nkyi−θi2+r2σk2~igσk2;r1+nk <dig> ∑i=1nkyi−θi2+r <dig>
additional file
additional file 1:
r code for dpgmm.
abbreviations
bcibayesian confidence interval
crpchinese restaurant process
dpgmmdirichlet process gaussian mixture model
emexpectation-maximization
gmmgaussian mixture model
gwasgenome-wide association studies
kcluster number
mcmcmarkov chain monte carlo
namnested association mapping
qtlquantitative trait locus/loci
sestandard error
snpsingle nucleotide polymorphism
competing interests
the authors declare that they have no competing interests.
authors’ contributions
xs and rm conceived and designed the study; xs conducted the meta-analyses of the published marker data, analysis of gwas data sets, and wrote the computer code to embody the method; xs and rm authored the manuscript.
all authors read and approved the final manuscript.
this research was supported in part by a grant from monsanto company, st. louis, mo, usa; x. sun was also supported in his graduate studies as a monsanto fellow in plant breeding through a generous gift of monsanto company to the university of illinois.
furthermore, we thank jianfeng xu  for his helpful suggestions regarding the dpgmm and anonymous reviewers for their constructive feedback which served to improve the manuscript.

 computer simulation is a resource which can be employed to identify optimal breeding strategies to effectively and efficiently achieve specific goals in developing improved cultivars. once identified, optimal strategies can be incorporated in the development ‘process’ to facilitate maximal genetic gains and accelerate the breeding process. through computer simulation, innovative approaches heretofore not feasible without use of current genomic-based technologies can be specifically tailored to meet the need at hand. in some instances, it is crucial to assess in silico, the options as well as the impact of various crossing schemes and breeding strategies on performance for certain key agronomic traits such as grain yield. such is the case with introgression of value-added traits by means of backcross breeding, the goal of which is the recovery of all the performance attributes of the elite variety or hybrid targeted for conversion along with the new genetic elements associated with the value-added trait. for this, a means by which the gene effects for the performance attributes of the line to be converted can be modeled in the simulation is particularly advantageous, yet it requires an accurate depiction of the distribution of gene effects. typically, recovery of the target line is estimated by the average proportion of genetic material carried through the backcrossing process and this estimate implicitly assumes that the many genes for key quantitative traits like grain yield are dispersed uniformly across the genome, each contributing only small effect. by including a more realistic representation of gene effects in the genome model to assess backcross breeding strategies, the means to most rapidly and effectively recover not only the germplasm per se, but the important genes contributing to performance of the variety or hybrid targeted for conversion, can be considered in evaluating strategies and approaches. bennewitz and meuwissen explored the distribution of additive and dominance effects of identified qtl from three f populations of pigs evaluated for meat quality and carcass traits, recognizing the value of modeling these effects, some of which are large, over use of the infinitesimal model which assumes an infinite number of qtl each with small effect. capitalizing on the large number of qtl studies, bennewitz and meuwissen conducted a meta-analysis of published qtl mapping data across traits to infer the distribution of additive qtl effects as well as dominance coefficients, fitting a gaussian mixture model. the idea of utilizing gmm is based on the notion that various qtl and associated genes fall into a number of classes of different-sized effects. the merit of employing gmm is its flexibility with different combinations of mixtures of normals leading to different shapes of the distribution. in the finite mixture model, the number of components k must be pre-specified. the value could be determined based on some specific information or criteria, such as the akaike information criterion or the bayesian information criterion. this requirement, frequently encountered in parametric statistics, can be sidestepped by introducing a nonparametric dirichlet process which assumes an infinite number of components. the dirichlet process is defined as a random process by which a sample drawn is a discrete distribution; it can be considered a ‘distribution over distributions’ and has been used widely in the field of population genetics to explain population structure. also desiring to capitalize on the large number of qtl studies, we took the meta-analysis concept a step further. we devised a method to represent the genomic distribution of additive and dominance gene effects associated with quantitative traits, which utilizes a modified dirichlet process gaussian mixture model to fit the number of mixture components and estimate parameters. as a departure from traditional dpgmm which only models qtl effects, we modified the model to be able to accommodate both qtl effects and their respective variances. utilizing previously-identified qtl for a number of quantitative traits in maize, the modified dirichlet process implements a chinese restaurant process to assign component membership and uses gibbs sampling to update conditional posterior distributions. our purpose in devising this method was to facilitate representation of the genetic architecture of grain yield and other key traits for use in computer simulations to optimize breeding strategies for multiple trait introgression. trait introgression involves backcross breeding and, therefore, bi-parental populations with no more than two potential alleles for a given locus; this was the primary scenario we intended to model. however, we also explored whether the method would be pertinent to other modeling activities that may involve estimation of breeding value or prediction of performance based on priors, scenarios which could involve multi-allelic populations in various crop or animal species. besides use in computer simulation, an accurate depiction of the qtl effects could contribute to a better understanding of the overall genetic architecture contributing to variation of expression of a particular trait of interest. methods description of dpgmm and priors to begin, we modeled the distribution of additive qtl effects and dominance coefficients using mixtures of normal distributions, namely gmm. the goal was to assign genetic effects to different mixture components. two latent variables were introduced: 1) the total number of mixture components and 2) the assignment of ith qtl effect to components. the gmm model was modified to accommodate the standard errors of qtl effects: pyi|λ …,λk~∑k=1kπknyi;μk,σk2+τi where yi is the ith observed qtl effect, τi is the known standard error of ith effect which is calculated from the qtl analysis, and λk = {πk, μk, σk2} is the kth parameter set, where variables πk, μk and σk are the mixing proportion, mean and variance of the kth mixture component, respectively. the dpgmm can be formulated hierarchically as follows: pyi|ci,Λ~nyi;μci,σci2+τi2ci|π1:k~discreteπ π …πkμkσk2~g0π π …πk~dirichletα/k,…,α/kΛ=λ λ …λk where g is a joint prior distribution for and mixing proportions π1: k are drawn from a symmetric dirichlet distribution with concentration parameter α. conditional on the mixing proportions, the latent indicator variables ci’s were sampled from a discrete distribution, specifically a multinomial distribution, and the prior for the ci in model could be written as a probability conditional on c− i: pci=k|c−i,α=n−i,k+α/kn−1+α, where n− i,k is the number of effects, not including yi that are linked with class k. and as k goes to infinity, the limits of the prior for the ci reach the following: pci|c−i,α=n−i,kn−1+αci=k,n−i,k>0αn−1+α∀i≠i',ci≠ci'. where i’ is the complement of the set i. as k → ∞, the dirichlet distribution becomes a dirichlet process in the limit. thus, infinite limit of model can be written as a dpgmm: yi~nyi|θi,σi2μiσi2~gg~dpαg where θi ~ n is a nuisance parameter, g is a random discrete distribution drawn from dp, and g was the base distribution, which specified the joint prior distribution of. given that the regular choice of priors for the mean and variance of the gaussian are normal and inverse gamma distributions, respectively, conjugate joint priors n* ig were chosen in the model in order to implement the following gibbs sampling. gibbs sampling in bayesian framework, unknown variables were sampled and updated from the conditional posterior distribution using markov chain monte carlo. considering the likelihood and priors in formulae and the full joint posterior distribution can be written as follows: pc,μ,σ2|y∝∏i=1nnyiθiσi2πθiμiσi2ci=∏i=1nnyiθiσi2nθiμiτi2nμiμ0σ02igσi2r1r2pci|c−i,α. unobservables were repeatedly sampled and updated from their posteriors, conditional on all other variables. the gibbs sampler was implemented as follows:initialization: assign initial values for where k =  and ci =   for i = 1: n. update θi: the conditional posterior distribution of θi was pθi|else∝nyiθiσk2nθiμkτi2∝nθiyiσk2+μkτi21σk2+1τi211σk2+1τi 3) update cluster indicators ci: the conditional posterior probabilities for ci were: pci=k|else∝nyiθiσk2nθiμkτi2pci|c−i,α=n−i,k2πσk2τi2exp−yi−θi22σk2−θi−μk22τi pci=k+1|else∝α∬nyiθiσk+12nθiμk+1τi2nμk+1μ0σ02igσk+12r1r2dμk+1dσk+12∝α∫nyiθiσk+12igσk+12r1r2dσk+12∫nθiμk+1τi2nμk+1μ0σ02dμk+1=α2πr2r1Γr1Γr1+1212yi−θi2+r2r1+121τi2+σ02exp−θi−μ022τi2+σ where Γ is the gamma function. note that constant 1n−1+α was omitted in both probabilities and were unknown and needed to be integrated out to leave ci as the only variable to be estimated from the markov chain. the dirichlet process was represented via the crp. effects were assigned to either currently holding cluster or a new cluster based on the above probabilities. if a new cluster was chosen, then the cluster size was increased, i.e. k + 1 → k. in case of n− i,k =  the kth cluster was eliminated and the cluster indicators were decreased by one, i.e. k → k −  ) resample and update suggested by formulae as per neal as follows: pμk|θi∈kthcluster,else∝∏i=1nknθiμkτi2nμkμ0σ02~nμk∑i=1nkθiτi2+μ0σ02∑i=1nk1τi2+1σ021∑i=1nk1τi2+1σ pσk2|yi∈kthcluster,else∝∏i=1nknyiθiσk2igσk2r1r2~igσk2;r1+nk ∑i=1nkyi−θi2+r where nk is the number of effects associated with the kth mixture component. the derivations of the fully conditional posterior distributions are detailed in the appendix.5) repeat steps to gibbs sampler was implemented with iterations of the mcmc to update conditional posterior distributions. the first samples were discarded as burn-in and the rest of the samples were used to construct joint posterior distribution. the hyper-parameters in algorithm were set to be α = , r1 =   r2 = , μ0 =   σ02 = . among hyperparameters, alpha was empirically set to based on the simulation results. convergence was checked by inspection of negative log-likelihood plots. after the burn-in period, when the mcmc converges to the stationary distribution, sampled parameters were collected to form the posterior distribution. we employed posterior means for estimating the mean and variance μ^kσ^k and posterior modes for estimating ĉi, which was further used to infer π^k. the bayesian confidence interval, which is the counterpart of the confidence interval in frequentist statistics, was defined as posterior probability that the parameter lies within the interval: ∫−∞apΛ|ydΛ=∫b∞pΛ|ydΛ=α/ where α is the significance level. instead of analytically estimating the confidence interval, the confidence interval for μ^kσ^k was numerically estimated from quartiles of posterior distribution. demonstration of method performance with simulated data to demonstrate the performance of the proposed method, two simulated data sets were processed. simulation i facilitated evaluation of model performance given complete data. it was generated from three gmms with respective means - and and variances of, and, respectively. a total of simulated qtl effects were uniformly distributed to the three components. this data set represents the case wherein all true qtl are known. simulation ii facilitated evaluation of model performance on a truncated distribution. truncation points were arbitrarily set to ± . the incomplete data set was intended to represent the common situation with qtl mapping data wherein all genetic effects are not detected, especially those with effects of near-zero magnitude. a truncated gaussian mixture with two mixture components was simulated. zero mean was assigned to both components. the first mixture component had mixing proportion π and variance σ of and, respectively; the second mixture component had π and σ of and, respectively. in both simulations, the standard error τi was generated from a uniform distribution. implementation with real data the model was also applied using real data to derive the distributions of additive effects and dominance coefficients. additive qtl effects were assembled from previous qtl mapping studies performed in corn. messmer et al. had evaluated recombinant inbred lines derived from a cross between two subtropical white dent maize lines to map genes controlling yield components and other traits, identifying qtl in total which are included in data i. briggs et al. had utilized a maize-teosinte backcross population to explore genes controlling domestication and morphological traits such as plant architecture, primary tassel and lateral inflorescence, identifying qtl in total which are included in data ii. data iii is derived from five maize qtl mapping studies involving segregating populations, all of which share a common parent b comprising a total of quantitative trait loci including a) qtl for kernel oil concentration mapped in an f population; b) qtl for root angle and plant height mapped in an f population; c) qtl for stalk digestibility and kernel composition mapped in a f population; d) qtl for stripe disease resistance mapped in an f population; and e) qtl for grain yield and yield components mapped in an f population under water-limited conditions. all qtl mapping studies employed either composite interval mapping or multiple interval mapping to detect qtl.. furthermore, these qtl studies all reflect estimates of gene effects in bi-parental populations, which fit with the backcross breeding scenario we intended to simulate. the histograms of observed additive effects, with values represented in units of phenotypic standard deviations, are shown in fig.  for all three data sets.table qtl associated with various traits across four data sets. data i, ii, iii, and iv were included in the analysis of qtl additive effects and data iii was used in the analysis of qtl dominance coefficients to explore application with multi-allelic, large-scale data sets as an alternative to bi-parental qtl data sets, additive qtl effects were assembled from high-resolution genome-wide association studies with the maize nam population, data provided courtesy of jason g. wallace. from the traits in the data set, were chosen to generate data iv: ‘20-kernel weight’ and ‘days to anthesis’. only those significant single nucleotide polymorphisms featured in at least three resample inclusions were utilized to fit the distribution. the qtl had been detected through composite interval mapping or multiple interval mapping and the association mapping had been performed using the forward-regression genome-wide association method in tassel.. results were compared to those obtained with qtl mapping data sets from bi-parental populations to contrast multi-allelic vs. bi-allelic, number of qtl, power in detection of the qtl, and single-trait distributions vs. distributions representing multiple traits. in addition, the distribution of dominance effects was explored. dominance coefficients, which are defined as the ratio between the observed qtl dominance deviation and absolute value of qtl additive effects, were assembled from data iii. the absolute value of additive qtl effects was used because the sign of qtl effect only signifies which parent had contributed the favorable allele, not the true direction of specific additive effect. the se for additive qtl effects and dominance coefficients was measured to take into account the experimental error. if logarithm of the odds scores for qtl were absent, standard errors were generated by taking sample standard deviation of effects from multiple experiments. the se of data ii and of data from dintinger et al. incorporated in data iii were produced in this way, where only those qtl detected in at least two environments were included in the final dataset. for the rest of the studies, ses were derived from lod scores as suggested by hayes and goddard or, in the case of the gwas data sets, the ses were determined from the sample errors of the discovered snp effects. standard errors of dominance coefficients were estimated by the delta method suggested by bennewitz and meuwissen, assuming no covariance between additive and dominance effects. specifically, sed/a=d/a*seaa2+sedd where sea and sed were standard errors for additive effects a and dominance effects d, respectively. additive qtl effects from qtl mapping studies were scaled by their corresponding phenotypic standard deviations in order to combine data across traits. the phenotypic standard deviations were computed using raw data if available from the qtl study. otherwise, the error variance and heritability of the trait were used to calculate the phenotypic standard deviation or, absent this, the range of phenotype values were used. phenotypic range was assumed to be standard deviations, which covers almost 100 % of the values, considering that most traits follow a normal distribution. we did not apply phenotypic standardization to the gwas data since an abundance of data points were available. since for the data sets from which dominance effects were generated, none of the above three conditions was fulfilled to obtain the phenotypic standard deviations for additive qtl effects, the additive effects from those five corn studies were not utilized in analyzing the distribution of additive qtl effects. note that the scaling process was not necessary for dominance coefficients, because the phenotypic standard deviation canceled out in the d/a ratio. due to limited statistical power of qtl mapping studies, many qtl with near-zero effects were likely not detected in the published studies used in this analysis, which is effectively analogous to a truncation of the additive qtl effect distribution near zero. faced with this issue, bennewitz and meuwissen suggested a “doubling” process, given the assumption that the true qtl effects occur at the highest frequency around zero. to compensate for the ‘missing’ qtl effects, a doubling of the data was done to ensure that the mean of each mixture component is estimated at zero, preserving the characteristic greatest density at zero for each cluster. with doubling, both signs for the same qtl additive effect were created. for example, for ith effect yi with se τi, − yi was added to the data with the same se. the above procedure leverages the fact that absolute values of alternative homozygous genotypes at a qtl are the same by definition e.g. in bi-parental populations. the “doubling” process was not applied to dominance coefficients, because most loci have observed effects around zero.  in this study, we employed a new method, namely dpgmm, to describe the distribution of qtl additive effects and dominance coefficients in the form of mixtures of normals. although similar to the fitting of a mixture of normals using a modified expectation-maximization algorithm, this approach differs primarily in the way of dealing with cluster size. with dpgmm, the number of mixture components does not need to be specified. in contrast, with use of a finite mixture model, the number of components needs to be preset and later decided under certain circumstances, or determined by some measure, e.g. akaike information criterion or bayesian information criterion. the optimum cluster size will strike a balance between maximum data compression and maximum accuracy. in dpgmm, the dirichlet process is represented by the crp, wherein a data point is assigned either to a currently occupied mixture component with probability proportional to the number of data already held in that cluster or to a new cluster with probability proportional to the concentration parameter, g by the same token, with each iteration of gibbs sampling, the cluster indicators are also updated along with parameters like the mean and variance. as such, dpgmm fits the data distribution and explores the potential number of mixture components simultaneously. to demonstrate the performance of dpgmm in fitting the distribution of qtl effects and to verify accuracy in estimating associated parameters of this distribution, two simulated data sets were processed. simulation i was structured to represent the case wherein all true qtl are known i.e. complete data. simulation i, which featured three components, resulted in a histogram of genetic effects from which it is difficult to infer the number of mixture components visually. in simulation i, dpgmm clearly fitted the data to three clusters with estimated values close to true values for parameters involving the proportion of mixing among clusters, the mean and the variance of each mixture component. dpgmm predicted accurately the mean and variance of clusters and although missed assignments of cluster membership were observed. in contrast, the mixing proportion of cluster was estimated precisely; however, certain deviations from the true mean and variance were observed.fig. histograms of simulated effects from gaussian mixtures: a three components having mean of - and and variance of, and, respectively, and equal mixing proportions for all three components; and two components having zero means and variance of and, respectively, and mixing proportions of and, respectively. distribution in is truncated at points - and k is the mixing proportion in the k th cluster, and μ k and σ k are the mean and variance of k th mixture component, respectively. values expressed in units of phenotypic standard deviation π k π^k μ k μ^k σ k2 σ^k2 simulation ii data, based on a truncated mixture of normals featuring two mixture components with zero mean, resemble a scenario common to qtl mapping wherein near-zero genetic effects were not detected. simulation ii produced a histogram with a pronounced gap around zero as expected in keeping with the data truncation. in simulation ii, parameters were estimated with accuracy, except for the variance of cluster which was estimated at versus the true value of.table true versus estimated parameters in simulation ii. π k is the mixing proportion in the k th cluster, and μ k and σ k are the mean and variance of k th mixture component, respectively. values expressed in units of phenotypic standard deviation π k π^k σ k2 σ^k2 given a complete set of data, dpgmm could clearly assign membership to respective clusters with small prediction errors. in the case of truncated data, dpgmm was still effective in predicting the correct number of mixture components and estimating the variance of components with greater variability; however, dpgmm was less effective in estimating variance of components with small differences among cluster members. as shown in table  the deviation of estimated variance from true value was somewhat large and might be attributed to the loss of small-value data in the sample. this result is in accord with the conclusion of bennewitz and meuwissen that small effects could be missed easily with a mixture model. proceeding to real data, dpgmm was used as a method to fit a mixture of normals for which the number of components is not known, fitting additive effects and dominance coefficients based on previously published qtl mapping data for a number of quantitative traits in maize. in addition, dpgmm was used to fit additive effects estimated from high-resolution gwas of the maize nam population to compare distributions produced with large-scale, multi-allelic data sets involving a single trait with those obtained with meta-analysis of bi-parental, lower-resolution studies involving multiple traits. the fitted distributions are the outcome of applying the dpgmm method and these distributions could then become the basis for modelling the genetic architecture of quantitative traits in maize for computer-simulated explorations to identify optimal breeding strategies. distributions of qtl additive effects histograms of observed qtl additive effects expressed in units of phenotypic standard deviation were generated from the qtl mapping studies. we noted that the histogram of data ii resembled that observed with simulation ii, wherein near-zero effects were not included. with data i, ii, and iii, it is difficult to infer the number of mixture components for additive effects visually from the histograms. however, the number of mixture components was inferred by the mode of posterior distribution with regard to cluster indicator ci. frequency tables of cluster membership clearly suggested fitting all data to one cluster for data i and data ii; however, for data iii, two components are suggested to fit additive effects. this likely reflects the types of traits included in each meta-analysis. data i contains primarily grain yield and yield component traits as well as some developmental traits. data ii is comprised of mainly domestication and morphological traits. data iii includes yield and yield component traits under drought stress as well as grain quality traits. results suggest more similarity among the traits within data i and within data ii, in contrast to data iii that was fitted with two clusters.fig. histograms of observed qtl additive effects: a data i; b data ii; and c data iii fig. histograms of the cluster number: a data i; b data ii; and c data iii for each of the three qtl mapping data sets, the fitted distribution was overlaid on the histogram of “doubled” data. the fitted distributions of additive effects all have zero mean; variances differed. the range of the observed qtl effects was tight around the mean for data i, whereas data ii showed some larger effects with absolute value nearing with data iii, the range of observed effects was greater still, with the largest effects nearing in absolute value; nonetheless, the variance of data iii additive effects was similar to that of the distribution of data ii. the variance was smallest for data i comprised of mainly yield and yield component traits.fig. fitted normal distributions to qtl additive effects: a data i; b data ii; c data iii; d data iv featuring traits of 20-kernel weight and days to anthesis th cluster, the cluster mean and variance and bayesian confidence interval for parameters in the distribution of additive effects and dominance coefficients. values expressed in units of phenotypic standard deviation σ^12 σ^12 π^1 π^2 σ^12 σ^22 μ^1 σ^12 σ^12 σ^12 comparing the results obtained with the data iv gwas single-trait data sets, the data were fitted to one cluster for both 20-kernel weight and days to anthesis. each fitted distribution was overlaid on the histogram of “doubled” data. each distribution had an estimated mean of zero. the variance associated with 20-kernel weight was estimated as, with a narrow confidence interval, to; the variance for days to anthesis was estimated as, with a narrow confidence interval, to. the fitted distribution of 20-kernel weight based on qtl shows a high degree of similarity to that for data i which is comprised of qtl for yield and yield component traits, maturity, and abiotic stress. the fitted distribution of days to anthesis based on qtl shows a high degree of similarity to that for data ii which is comprised of qtl for domestication and morphological traits. compared to distributions from meta-analyses, the gwas data sets resulted in greater precision as displayed in the relatively smaller bci’s for the estimates of variance; the ratio of the size of the bci’s for data iv 20-kernel weight and data i is, and for data iv days to anthesis and data ii is. thus, the accuracy using meta-analysis of multiple-trait qtl and qtl identified through gwas appears to be similar when the type of traits measured are similar. however, precision is better with the latter, which likely reflects the significantly larger number of qtl in the data iv sets. despite the significantly greater resolution of the gwas data sets, both distributions produced from data iv displayed a gap around zero: within ± for 20-kernel weight and within ± for days to anthesis. bennewitz and meuwissen discussed the potential drawbacks of using meta-analyses of qtl effects detected in bi-parental populations to characterize distributions of gene effects. a primary concern was for failure to identify all true qtl within the identified qtl sets either because certain alleles were excluded from the data set or because of a lack of statistical power to detect smaller-effect qtl. even in cases where marker density and genome resolution is high, all qtl for a trait of interest may not be identified if mapping methods are utilized that screen detected effects against a significance threshold which basically excludes most of the small effect qtl. furthermore, in genomic selection, the particular statistical method used to ‘train’ the model can influence the distribution of effects. despite the higher resolution and multi-allelic nature of the data iv data sets, the inability to detect qtl of near-zero effect was apparent, yet not problematic in fitting a distribution. clearly, the distributions are centered at zero with highest frequencies observed in the zero vicinity and recognized as single-component. thus, the results obtained with the gwas data sets demonstrate the robustness of the dpgmm method with use of either meta-analysis of qtl identified across traits or large-scale data sets comprised of qtl for a single trait. if enough data are available, e.g. gwas dataset, distributions of qtl effects based on a single trait rather than across traits might be more appropriate and useful. and with the advent of genotype-by-sequence and other genotyping technologies that facilitate high-resolution marker sets, very dense marker sets are more available and more widely utilized in qtl identification. for our purpose in exploring the impact of backcrossing strategies on recovery of performance of the hybrid targeted for conversion, the meta-analysis involving the data i set based on a bi-parental population met our objectives well in that it focused on yield and yield component traits and additionally included other key traits essential to performance recovery i.e. maturity and abiotic stress tolerance. distribution of qtl dominance coefficients observed dominance coefficients obtained from meta-analysis of five mapping studies varied in magnitude from less than - to more than, suggesting that all classes of dominance were represented among the traits measured. around 50 % of the qtl displayed an d/a ratio in the range of - to, indicating partial recessivity, additivity, and partial dominance gene action. approximately 25 % of the qtl exhibited either partially dominant or dominant gene action or partially recessive or recessive gene action. furthermore, 25 % of the qtl exhibited apparent overdominance or underdominance gene action.fig. histogram of observed dominance coefficients from meta-analysis based on five mapping populations dominance coefficients were fitted with the normal distribution using dpgmm. the mode of posterior distribution with regard to k is suggesting that all data could be fitted to a single component; estimates of the distribution mean and variance were provided by mcmc as well. figure  displays the estimated distribution overlaid on the density plot of observed data. the estimated mean of the distribution was with 95 % bci of to, and variance of with 95 % bci of to. the result that the distribution of dominance coefficients was fitted to a normal distribution with a positive mean conformed to previous studies.fig. estimation of cluster number, mean, and variance of the fitted distribution of dominance coefficients through mcmc fig. normal distribution fitted to the dominance coefficients, with estimated mean at with 95 % bayesian confidence interval to be and and estimated variance at with 95 % bayesian confidence interval to be and  the dpgmm method offers an alternative to the over-simplified infinitesimal model in computer simulation as a means to better represent the genetic architecture of quantitative traits, which likely involve some large effects in addition to many small effects. furthermore, it confers an advantage over other methods in that the number of mixture model components need not be known a priori. the dpgmm method takes advantage of prevalent qtl data to approximate the distributions of additive and dominance gene effects. the dpgmm method is robust with use of either meta-analysis of smaller-scale qtl analyses involving a number of traits or large-scale, single-trait qtl data sets. furthermore, qtl data sets from bi-allelic or multi-allelic populations can be utilized. although the data sets may be missing some near-zero qtl effects that were not resolved in the fundamental analyses, the methodology is able to accommodate this drawback. r code to facilitate use of the dpgmm method is included in a additional file to this paper. the distributions of qtl additive and dominance effects highlighted through this study were used in modelling the genetic architecture of grain yield and other key performance traits for computer-simulated explorations to identify optimal breeding strategies to facilitate introgression of multiple value-added traits into an elite maize hybrid. maize grain yield is a complex trait involving dominant and over-dominant gene action. other traits important to recovery of the essential performance attributes of the hybrid targeted for conversion include maturity, resistance to lodging, and abiotic stress tolerance and these have a bearing on yield performance as well. readers are directed to the recent work of sun and mumm for an example of how the dpgmm-estimated genetic distribution parameters were deployed in computer simulations to evaluate breeding strategies. appendix derivations of the fully conditional posterior distributions are given as follows: pθi|else∝nyiθiσk2nθiμkτi2∝exp−12*yi−θi2σk2*−12*θi−μk2τi2∝exp−12*θi−yi2σk2+θi−μk2τi2∝nθiyiσk2+μkτi21σk2+1τi211σk2+1τi pci=k|else∝nyiθiσk2nθiμkτi2pci|c−i,α∝12πσk2exp−yi−θi22σk2*12πτi2exp−θi−μk22τi2*n−i,k=n−i,k2πσk2τi2exp−yi−θi22σk2−θi−μk22τi pci=k+1|else∝α∬nyiθiσk+12nθiμk+1τi2nμk+1μ0σ02igσk+12r1r2dμk+1dσk+12∝α∫nyiθiσk+12igσk+12r1r2dσk+12∫nθiμk+1τi2nμk+1μ0σ02dμk+ the first integration part: ∫nyiθiσk+12igσk+12r1r2dσk+12=∫12πσk+12−12exp−yi−θi22σk+12*r2r1Γr1σk+12−r1−1exp−r2σk+12dσk+12=∫12πr2r1Γr1σk+12−r1+12−1exp−12yi−θi2+r2σk+12*dσk+12=12πr2r1Γr1Γr1+1212yi−θi2+r2r1+12∫igσk+12;r1+ yi−θi2+r2*dσk+12=12πr2r1Γr1Γr1+1212yi−θi2+r2r1+ the second integration part: ∫nθiμk+1τi2nμk+1μ0σ02dμk+1=∫12πτi2exp−12θi−μk+12τi212πσ02exp−12μk+1−μ02σ02dμk+1=∫12πτi212πσ02exp−12θi−μk+12τi2−12μk+1−μ02σ02dμk+1=12πτi212πσ02∫2π11τi2+1σ02μk+1θiτi2+μ0σ021τi2+1σ0211τi2+1σ02exp−12θi−μ02τi2+σ02dμk+1=exp−12θi−μ02τi2+σ0212πτi212πσ022π11τi2+1σ02=12πτi2+σ02exp−θi−μ022τi2+σ together we have: pci=k+1|else=α2πr2r1Γr1Γr1+1212yi−θi2+r2r1+121τi2+σ02exp−θi−μ022τi2+σ pμk|θi∈kthcluster,else∝∏i=1nknθiμkτi2nμkμ0σ02∝exp−12∑i=1nkθi−μk2τi2−12μk−μ02σ02~nμk∑i=1nkθiτi2+μ0σ02∑i=1nk1τi2+1σ021∑i=1nk1τi2+1σ pσk2|yi∈kthcluster,else=∏i=1nknyiθiσk2igσk2r1r2∝σk2−nk2exp−12∑i=1nkyi−θi2σk2σk2−r1−1exp−r2σk2=σk2−r1+nk2−1exp−12∑i=1nkyi−θi2+r2σk2~igσk2;r1+nk ∑i=1nkyi−θi2+r additional file additional file 1: r code for dpgmm. abbreviations bcibayesian confidence interval crpchinese restaurant process dpgmmdirichlet process gaussian mixture model emexpectation-maximization gmmgaussian mixture model gwasgenome-wide association studies kcluster number mcmcmarkov chain monte carlo namnested association mapping qtlquantitative trait locus/loci sestandard error snpsingle nucleotide polymorphism competing interests the authors declare that they have no competing interests. authors’ contributions xs and rm conceived and designed the study; xs conducted the meta-analyses of the published marker data, analysis of gwas data sets, and wrote the computer code to embody the method; xs and rm authored the manuscript. all authors read and approved the final manuscript. this research was supported in part by a grant from monsanto company, st. louis, mo, usa; x. sun was also supported in his graduate studies as a monsanto fellow in plant breeding through a generous gift of monsanto company to the university of illinois. furthermore, we thank jianfeng xu for his helpful suggestions regarding the dpgmm and anonymous reviewers for their constructive feedback which served to improve the manuscript.
the method, based on meta-analysis of previously-estimated qtl effects following bennewitz and meuwissen, utilizes a modified dirichlet process gaussian mixture model to fit the number of mixture components and estimate parameters of the genomic distribution. to address this need, we devised a method to represent the genomic distribution of additive and dominance gene effects associated with quantitative traits. the method was demonstrated using several maize qtl data sets to provide estimates of additive and dominance effects for grain yield and other quantitative traits for application in maize genome simulations. the dpgmm method offers an alternative to the over-simplified infinitesimal model in computer simulation as a means to better represent the genetic architecture of quantitative traits, which likely involve some large effects in addition to many small effects. thus, the method has high utility in representing the genetic architecture of quantitative traits in computer simulation. in addition, the method is robust with use of large-scale, multi-allelic data sets or with meta-analyses of smaller qtl data sets which may be derived from bi-parental populations in precisely estimating distribution parameters. for this, a means by which gene effects can be represented in the genome model is critical. computer simulation is a resource which can be employed to identify optimal breeding strategies to effectively and efficiently achieve specific goals in developing improved cultivars. furthermore, it confers an advantage over other methods in that the number of mixture model components need not be known a priori. dirichlet processgaussian mixture modelgenetic architecturemeta-analysisqtl effectsmonsanto companyissue-copyright-statement© the author 2016 in some instances, it is crucial to assess in silico the options as well as the impact of various crossing schemes and breeding approaches on performance for traits of interest such as grain yield.
 computer simulation is a resource which can be employed to identify optimal breeding strategies to effectively and efficiently achieve specific goals in developing improved cultivars. in some instances, it is crucial to assess in silico the options as well as the impact of various crossing schemes and breeding approaches on performance for traits of interest such as grain yield. for this, a means by which gene effects can be represented in the genome model is critical.  to address this need, we devised a method to represent the genomic distribution of additive and dominance gene effects associated with quantitative traits. the method, based on meta-analysis of previously-estimated qtl effects following bennewitz and meuwissen, utilizes a modified dirichlet process gaussian mixture model to fit the number of mixture components and estimate parameters of the genomic distribution. the method was demonstrated using several maize qtl data sets to provide estimates of additive and dominance effects for grain yield and other quantitative traits for application in maize genome simulations.  the dpgmm method offers an alternative to the over-simplified infinitesimal model in computer simulation as a means to better represent the genetic architecture of quantitative traits, which likely involve some large effects in addition to many small effects. furthermore, it confers an advantage over other methods in that the number of mixture model components need not be known a priori. in addition, the method is robust with use of large-scale, multi-allelic data sets or with meta-analyses of smaller qtl data sets which may be derived from bi-parental populations in precisely estimating distribution parameters. thus, the method has high utility in representing the genetic architecture of quantitative traits in computer simulation. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords dirichlet processgaussian mixture modelgenetic architecturemeta-analysisqtl effectsmonsanto companyissue-copyright-statement© the author 2016 