A fast algorithm for BayesB type of prediction of genome-wide estimates of genetic value

Genomic selection uses genome-wide dense SNP marker genotyping for the prediction of genetic values, and consists of two steps: (1) estimation of SNP effects, and (2) prediction of genetic value based on SNP genotypes and estimates of their effects. For the former step, BayesB type of estimators have been proposed, which assume a priori that many markers have no effects, and some have an effect coming from a gamma or exponential distribution, i.e. a fat-tailed distribution. Whilst such estimators have been developed using Monte Carlo Markov chain (MCMC), here we derive a much faster non-MCMC based estimator by analytically performing the required integrations. The accuracy of the genome-wide breeding value estimates was 0.011 (s.e. 0.005) lower than that of the MCMC based BayesB predictor, which may be because the integrations were performed one-by-one instead of for all SNPs simultaneously. The bias of the new method was opposite to that of the MCMC based BayesB, in that the new method underestimates the breeding values of the best selection candidates, whereas MCMC-BayesB overestimated their breeding values. The new method was computationally several orders of magnitude faster than MCMC based BayesB, which will mainly be advantageous in computer simulations of entire breeding schemes, in cross-validation testing, and practical schemes with frequent re-estimation of breeding values.


Introduction
The recent detection of thousands to millions of SNP markers and the dramatic improvements in high-throughput, cost effective genotyping technology have made it possible to apply marker assisted selection at a genome wide scale, which is termed genomic selection [1]. These authors suggested three methods for the estimation of genetic value from dense SNP marker data, namely GS-BLUP, BayesA, and BayesB. GS-BLUP applies the BLUP approach to the estimation of the effects of the marker alleles, which assumes a normal prior distribution for the marker effects, where the variance of the prior distribution was assumed equal for all the markers. Since an equal variance for each of the marker effects seems unrealistic, the BayesA method extended the GS-BLUP method by estimating the variance of every marker separately, and an inverse chi-square prior was used for the estimation of these variances. In the BayesB method it was assumed that many of the markers will actually have no effect, and the prior distribution of the variances was a mixture of a distribution with zero variance and an inverse chi-squared distribution. In a simulation study where the genetic model included a finite number of loci with exponentially distributed effects, BayesB provided more accurate prediction of genetic value than BayesA, which in turn was more accurate than GS-BLUP.
Although BayesB has the potential for the development of more faithful genetic models, and so seems the method of choice for estimating genome wide breeding values (GW-EBV), its calculation requires the use of computer intensive MCMC techniques. For practical applications and for computer simulations of genomic selection breeding schemes, where many selection rounds and replications are required, it would be advantageous if a much faster algorithm for the calculation of BayesB GW-EBV would be available. Thus, our aim here is to present a fast non-MCMC based algorithm for the calculation of BayesB type estimates of GW-EBV. By using a mixture of a distribution with zero effects and an exponential distribution as a prior for the marker effects, the integration involved in calculating the expectation of the breeding values given the data can be solved analytically, which makes non-MCMC estimation of GW-EBV possible.

The model
We will first develop a model for the estimation of the effect of one SNP. The model will be extended to m SNPs in the fourth section of Methods, where m will be assumed to be much larger than the number of records n. When estimating one SNP effect, we generally do not need to use prior information, since the prior will usually be overwhelmed by the information from the data. However, since we will be extending the model to the estimation of many SNP effects, we will use prior information here. Also, in order to facilitate expansion to multi-SNP models we assume that the SNP considered is a randomly chosen SNP with a random effect, instead of a particular SNP with some particular effect. We will assume that a SNP marker has two alleles, 0 and 1, with allele 1 the reference allele having frequency p. For the purpose of estimating gene effects, the SNP genotypes were standardised to form an n × 1 vector of covariates, b, defined for animal i as: (1) b i = -2p/SD, if i has homozygous genotype '0_0'; (2) Thus, the b i are standardised such that their mean is 0 and their variance is 1 in a random mating population. If g is the effect of the SNP in the population, and we assume that it is drawn from a distribution of marker effects with mean 0, then the variance due to the marker is Var(b i g) = Var(b i )Var(g) = Var(g). Therefore by model-ling the variance of the effect of the SNP, we also model the variance directly due to the SNP, irrespective of its frequency. This way of modelling the variance due to the SNP will be discussed in detail in Discussion. It may be noted however that the standardisation of the variance of the b i is not essential to our derivations below, i.e. when considered more appropriate we could choose the b i as -2p, (1-2p), and 2(1-p), respectively. The model for the records denoted by the n × 1 vector y is: where g is the effect of the SNP; and e is an n × 1 vector of residuals which are assumed normally distributed with variance  e 2 . We will ignore the estimation of an overall mean for now since, even if it were included in the model, its estimate would not affect the equation for the estimation of g since 1 ' b = 0, following from the standardisation of b.

Prior distribution of g
The prior distribution of the SNP effect g, (g), is assumed to be a mixture of a distribution with a discrete probability mass of zero and an exponential distribution reflected about g = 0 (see [2]): where  is the fraction of the SNPs that are in linkage disequilibrium (LD) with a quantitative trait locus (QTL) and may consequently have a non-zero effect; and  is the parameter of the exponential distribution. As described in the Appendix 1 this may be re-expressed as where (g) is the Dirac delta function. Therefore since the variance of the zero-reflected exponential distribution is 2 -2 , for m markers the total genetic variance is = 2m -2 . Covariances between marker effects are expected to be zero, because any non-zero covariance will depend on the coding of the marker alleles which is arbitrary, i.e. a positive covariance can be changed in a negative one by recoding the marker alleles. The hyper-parameters of the prior distribution ( and ) are assumed known here. If one is willing to assume that a fraction  of the markers are in LD with QTL, the variance per marker in LD with a QTL is  a 2 / (m), and an estimate of the  hyper-parameter is (2 m/  a 2 ).

The posterior distribution and the expectation of g | y
For the likelihood of the data y given the genetic effect g a multi-variate normal distribution is assumed with mean bg and variance I e 2 from (2), where I is the n × n identity matrix. Thus, the likelihood function is (y; bg, I ), where (y; , V) is the multi-variate normal density function with mean  and (co)variance matrix V. In what follows the dimensionality of the multi-variate (.;.,.) will be left implicit from the parameters and will include the univariate density. The summary statistic for g is Y = (b ' b) -1 b ' y = b ' y/n and, consequently, all information on g contained in y is also contained in Y. The variance of Y is  2 = (b ' b) -1  e 2 =  e 2 /n since b'b = n. Therefore, the likelihood function (y; bg, I ) is proportional to the univariate density (Y; g,  2 ), as shown in Appendix 2, and so we will use the univariate version of the likelihood function for simplicity. The posterior distribution of g|y now becomes: The posterior distribution is not affected by the use of either (y; bg, I ) or (Y; g,  2 ) since they are proportional to each other, i.e. (y; bg, I )/(Y; g,  2 ) = constant, and this constant is a factor of both the numerator and the denominator. The expectation of g given y is (e.g. [2]) as E [g|y] =  gp(g|y)dg and thus: Both the integral in the numerator and that in the denominator are analytically derived in full in Appendix 1. This results in a closed form for E [g|y], which is presented in the Appendix 1 in Equation (B3).

Extension to m SNPs
The extension to the estimation of the effects of m SNPs is obtained by the use of a modification of the Iterative Conditional Mode (ICM) algorithm [3], which we will call the Iterative Conditional Expectation (ICE) algorithm. The ICE algorithm uses the expectation/mean instead of the mode of the posterior, because (1) the mean of g maximises the correlation between g and [2]; and (2) due to the spike at zero the posterior may be bi-modal (see Results), in which case the mode may be quite far away from the mean. The ICE algorithm iteratively calculates E(g|Y) for each SNP in turn, using the current solutions of all the other SNPs as if they were true effects, which is known to only approximately converge to the correct solution. We will first describe the algorithm and next its approximate nature. The effects of the SNPs are denoted g i (i = 1, .., m), and in vector notation by g. When estimating the effect of the i th SNP, the current estimates of all other SNPs are used to calculate Y i and  2 . Iterating from a starting solution, e.g. = 0, the algorithm performs within each iteration the following steps: For all SNPs i = 1, .., m, Step 1: calculate 'adjusted' records, y -i , which are corrected for all the other SNPs so . Estimate the sufficient statistics , and  2 =  e 2 /n.
Step 2: Equation B3 of Appendix 1 is used to calculate , which is used to update the solution for marker i.
The iteration is stopped when the changes in the solutions become small, i.e.
, where subscript q denotes the q th iteration. In the Step 1, it is computationally more efficient to calculate Y i directly as . The advantage being that the elements of the normal equations matrix may be stored, which speeds up the calculation of Y i .
As mentioned before, if the only fixed effect is the overall mean, no correction for fixed effects is needed. If it is required to estimate fixed effects  with design matrix X, the calculation of y -i in the Step 1 becomes and each iteration also updates the solutions for the fixed effects e.g. by .
The approximate nature of the ICE algorithm is due to y -i and thus Y i not being known, but being estimated. For ease of notation, this was not denoted by a hat elsewhere in this paper, but in this paragraph it is necessary to make the distinction between the true values of where Var( ) is due to prediction error variances and covariances of the g j . It follows that The latter is probably only the case if there are few markers j  i. The first section of Results investigates the non-linearity, i.e. second and higher order derivatives, of the E(g i |Y i ) function.

Comparison of non-and MCMC based BayesB
The BayesB algorithm developed here will be denoted by fBayesB (for 'fast BayesB') and the BayesB algorithm using MCMC (as described in [1]) will be denoted by MBayesB. fBayesB and MBayesB were compared in twenty simulated populations. The simulation of the population and analysis by MBayesB was conducted by Solberg et al. [4] and their paper gives a detailed description of the simulation procedures. In brief, the populations were simulated for 1000 generations at an effective size of 100 in order to create mutation drift balance and LD between the markers and the genes. The genome consisted of 10 chromosomes of 100 cM with a total of 8010 equidistant marker loci including markers at each end of the chromosome. The mutation rate per marker locus per meiosis was high (2.5*10 -3 ) to ensure that virtually all the loci were segregating. If more than one mutation occurred at a marker locus, the mutation that resulted in the highest Minor Allele Frequency (MAF) was considered as 'visible', whereas the others were considered 'unvisible' and were thus ignored. With these parameters, the visibility procedure turned 99% of markers into biallelic markers, even if several mutations had occurred, with the remaining 1% being monoallelic. There were 1000 equidistant putative QTL positions, which were chosen such that a QTL position was always in the middle between two marker loci. Whether or not a putative QTL position had an effect on the trait depended on the existence of a mutation during the simulated generations, which occurred with a mutation rate of 2.5*10 -5 . The effects of QTL alleles were sam-pled from a gamma distribution with scale parameter 1.4 and shape parameter 4.2, and were considered with equal probabilities of being positive or negative.
The scale parameter of the gamma distribution was chosen such that the expected genetic variance was 1 (as in [1]). In generation 1001 and 1002 the population size was increased to 1000, and the animals were genotyped for the 8010 markers. In generation 1001, the animals also had phenotypic records, i.e. the phenotype of animal i was: where a ij(p) is the effect of the paternal (m, maternal) allele that animal i inherited at QTL locus j; e i~N (0, ), where was set equal to the realised genetic variance in the replicate, such that the heritability was 0.5 in each replicate. Marker effects were estimated using the phenotypes and genotypes obtained from generation 1001. The total genetic values for the animals of generation 1002 were predicted (GW-EBV) from their marker genotypes and the estimates of the marker effects. The accuracy of these estimates was calculated as the correlation between the GW-EBV and the true breeding value, which is known in this simulation study. The coefficient of the regression of the true genetic value on the GW-EBV was used as a measure of the bias of the EBV, and a regression coefficient of 1 implies no bias.

The non-linear regression curve
In Figure 1, E [g|Y] is plotted against the value of Y with  2 = 1; since Y is the sufficient statistic for g given the data y, E [g|y] = E [g|Y]. Figure 1 shows also the regression curve when the integrals in Equation (4) were numerically evaluated, as was done by Goddard [2]. The empirical curve of Goddard has similar characteristics, which is relatively flat at Y = 0, but approaches a derivative of 1 for extreme values of Y. However as a result of the closed expression in Appendix 1 (B3) it is possible to explore the full solution space and Figure 2 The expectation of the genetic value given the summary statistic of the data Y, E(g|Y), as a function of Y Figure 1 The expectation of the genetic value given the summary statistic of the data Y, E(g|Y), as a function of Y. The parameter of the exponential distribution is  = 1,  2 = 1, and the probability of a marker having a true effect is  = 0.05; E(g|Y) calculated by numerical integration is represented by black dots and the analytical solution is shown as white dots.    (5) suggests that the estimates of g will be conservative in the sense that they will show a bias towards zero. This bias will be increased for smaller , i.e.
for a smaller number of QTL over number of markers ratio. Figure 3 depicts the prior, likelihood and posterior distribution of g, in a situation where the posterior is bimodal.
For smaller values of Y, the likelihood curve moves to the left and the two peaks of the posterior merge into one. For larger values of Y, the likelihood moves to the right, hence its value at 0 reduces, and the peak of the posterior at zero disappears. In the latter case, the posterior becomes approximately a symmetric distribution, and the use of The probability density function (PDF) of the prior (...), likelihood (_ _ _), and posterior (___) of g

The accuracy and bias of GW-EBV using BLUP, fBayesB and MBayesB
MBayesB is slightly more accurate than fBayesB (Table 1). The difference in accuracy is 0.011 with s.e. of difference of 0.005, which is statistically significant. For comparison, the accuracy of BLUP GW-EBV is also shown, calculated with a prior that each marker contributes an equal amount of variance, namely 1/m, i.e. the infinitesimal model is assumed to hold, at least approximately. The accuracy of the BLUP GW-EBV is considerably lower than that of MBayesB and fBayesB, which is probably because the genetic model underlying the simulations is quite far from the infinitesimal model. Both MBayesB and fBayesB yielded biased EBV, but, interestingly, the biases are opposite in direction. MBayesB yields EBV that are too variable, so the EBV require to be shrunk in order to predict the TBV without bias, hence the regression of the TBV on the EBV is <1. In contrast, fBayesB yields EBV with too little variance, so the EBV require to be scaled up in order to predict the TBV without bias, i.e. the regression is >1. This conservative behaviour of , i.e. differences between the estimates are smaller than the real life differences, is expected based on the nonlinearity observed in Figure 2 and Equation (5) (see first section of Results). Especially, the bias of the fBayesB is considerable, and should be corrected for by rescaling the EBV when used in a breeding scheme where EBVs based on different amounts of information are to be compared.
In order to investigate any systematic deviations of fBayesB-EBV from the MBayesB-EBV, which is considered the golden standard, Figure 4 plots both types of EBV against each other. The regression of the EBV from both types of non-linear regression methods against each other seems pretty linear, which indicates that both methods seem to non-linearly regress the phenotypic data in a very similar way.

Computer time
The difference in computer time is large: whereas fBayesB takes 2 to 5 minutes MBayesB takes 47 h to compute. No attempt was made to implement parallel computing witĥ g   MBayesB-EBV

fBayesB-EBV
MBayesB, although it is readily amenable to parallel procedures through running multiple shorter MCMC chains (each chain should be longer than the burn-in period), which reduces the required wall-time for the calculations. The impact of parallel computing procedures of fBayesB is less clear, but there also appears to be no need for this.

Discussion
A fast, non-MCMC based algorithm, called Iterative Conditional Expectation (ICE), was derived for the calculation of GW-EBV using dense SNP genotype and phenotypic data. The speed of improvement is due to the analytical integration of the integrals involved in the calculation of E [g|Y]. The Bayesian estimation model used here has very similar characteristics to BayesB as described by [1] and denoted MBayesB here: the prior distribution is a mixture of a heavy-tailed distribution (fBayesB: exponential distribution; MBayesB: a normal distribution whose variance is sampled from an inverse chi-squared) and a distribution with zero effects. The latter mixture distribution is also called a 'spike and slab' mixture [5]. The QTL effects had a Gamma distribution, and thus differed from both that of fBayesB and MBayesB. It may be noted that the prior distribution of fBayesB and MBayesB does not apply to the QTL effects but to the marker effects, for which the true prior distribution will be hard to derive. Although we cannot rule out that the slightly better accuracy of MBayesB compared to fBayesB (Table 1) is due to the prior distribution of MBayesB being better than that of fBayesB, the non-linear regressions of MBayesB and fBayesB seemed very similar since they resulted in a linear regression of the fBayesB-EBV on the MBayesB-EBV ( Figure 4). The difference in accuracy seems to be due to that the ICE algorithm ignores higher order derivatives of E(g|Y) function to Y (Eqn. 5). Relaxing this assumption requires (1) taking second order derivatives of the E(g|Y) function, (2) calculation of prediction error variances of , and more research is needed to perform these calculations computationally efficiently.
The justification for the 'spike and slab' prior distribution is that many of the SNPs will not be in LD with a QTL and thus have no effect, whereas the SNPs that are in LD with a QTL have a distribution of effects that is similar to that of the QTL, albeit smaller in magnitude due to the need for several markers to predict the effect of the true QTL genotypes. The true distribution of QTL is often reported to be exponential or gamma [6]. Hayes and Goddard [6] found a shape parameter for the gamma distribution of 0.4, i.e. a leptokurtic shape similar to that of the exponential distribution. Where the marker is not in perfect LD with the gene, it will pick up only a fraction of the gene effect and the impact of this on the distribution of marker effects is included within the assumptions concerning their prior distribution. The variance due to a marker is Var("Marker") = Var(b i )Var(g). Here we standardised the variance of the genotypes to Var(b i ) = 1, i.e. the prior variance assumed for the marker effects, g, applies directly for the variance due to the marker, and thus does not depend on marker frequency. We prefer this parameterisation, assuming that the variance of the marker effects is frequency dependent, because (1) QTL with large effects are expected to be at rare frequencies, which implies that the variance of the QTL is roughly constant (at least considerably more constant than when QTL effects were not frequency dependent); (2) if we assume that the QTL variance needs to be above a certain threshold before the markers pick up its effect, these QTL will have a much more constant variance than randomly picked QTL. The algorithm is equally capable of handling the coding of b i as 0, 1, and 2, (or after subtracting the mean -2p, (1-2p), and 2(1-p)) for the genotypes mm, Mm, and MM, respectively. The accuracies of the GW-EBV were virtually the same as those in Table 1 when the latter parameterisation of the b i was used (result not shown).
The computational advantage of our fast algorithm for the BayesB approach to GW-EBV will not outweigh the reduced accuracy observed, if confirmed for typical trait architecture, when used in practical breeding schemes and the computational time and effort can be afforded. If, in practice, breeding schemes wish to select upon GW-EBV that require frequent updating, then a more appropriate comparison is between frequently updated fBayesB estimates of marker effects and the use of 'old' MCMC based Σb g j ĵ estimates of marker effects, where the GW-EBV of animals are calculated without updating the estimates of marker effects, because of computational constraints. In the latter case there will be a loss of accuracy. It will depend on the amount of new information coming into the breeding value evaluation system, which of these alternatives should be favoured. In simulations of breeding schemes and in cross-validation testing of GW-EBV, the large number of EBV evaluations required may make our fast algorithm the only means to implement BayesB type genome-wide breeding value estimation.