Joint prediction of multiple quantitative traits using a Bayesian multivariate antedependence model.

Predicting organismal phenotypes from genotype data is important for preventive and personalized medicine as well as plant and animal breeding. Although genome-wide association studies (GWAS) for complex traits have discovered a large number of trait- and disease-associated variants, phenotype prediction based on associated variants is usually in low accuracy even for a high-heritability trait because these variants can typically account for a limited fraction of total genetic variance. In comparison with GWAS, the whole-genome prediction (WGP) methods can increase prediction accuracy by making use of a huge number of variants simultaneously. Among various statistical methods for WGP, multiple-trait model and antedependence model show their respective advantages. To take advantage of both strategies within a unified framework, we proposed a novel multivariate antedependence-based method for joint prediction of multiple quantitative traits using a Bayesian algorithm via modeling a linear relationship of effect vector between each pair of adjacent markers. Through both simulation and real-data analyses, our studies demonstrated that the proposed antedependence-based multiple-trait WGP method is more accurate and robust than corresponding traditional counterparts (Bayes A and multi-trait Bayes A) under various scenarios. Our method can be readily extended to deal with missing phenotypes and resequence data with rare variants, offering a feasible way to jointly predict phenotypes for multiple complex traits in human genetic epidemiology as well as plant and livestock breeding.

Introduction

In human genetic epidemiology, accurate prediction of disease risk is vital for disease prevention and personalized medicine. Although thousands of genome-wide association studies (GWAS) have discovered a large number of variants that are significantly associated with complex human traits and diseases (Hindorff ), it is still a challenge to directly implement these findings to predict yet-to-be observed phenotypes and thus to advance preventive and personalized medicine. Partly, this is owing to the so-called ‘missing heritability' that these associated variants can typically account for only a small fraction of the total genetic variance (Manolio ). The issue of ‘missing heritability' probably results from the limitations of statistical techniques commonly used in GWAS, as pointed out by previous studies (Yang ; de los Campos ; Makowsky ); for example, in single-locus association analyses for either SNP markers or gene expression levels, we can merely derive a small number of statistically significant loci even for a complex trait affected by a large number of small-effect variants.

To address the weakness of conventional GWAS on phenotype prediction, some of recent studies have turned to take advantage of the whole-genome prediction (WGP) method initially proposed by Meuwissen . WGP seeks to model genome-wide SNPs simultaneously to predict yet-to-be observed phenotypes, for example, human height (Makowsky ; de Los Campos ), human lifespan (de los Campos ) and skin cancer risk (Vazquez ). These studies have demonstrated the strength of WGP methods for predicting human complex traits of which the underlying genetic architecture likely consists of numerous variants of small effect. In addition, a number of studies suggested that many human diseases, like schizophrenia and bipolar disorder (Purcell ), multiple sclerosis (Bush ) and rheumatoid arthritis (Stahl ), also have such a polygenic architecture, implicating a broad utility for WGP methods to predict human disease risk. Besides phenotype prediction in human genetics, WGP has been receiving much attention and widely employed to genetically improve economically important traits in domestic animals, including milk production performance and fitness traits of dairy cows (VanRaden ).

Since the seminal work of Meuwissen for predicting genomic breeding values in animal and plant breeding, a number of WGP methods have been developed and extensively investigated based on different algorithms (Gianola ; VanRaden, 2008; Gianola ; Habier ; Legarra ). Among these methods, the most popular are the parametric approaches, for example, GBLUP (VanRaden, 2008) and the Bayesian alphabet (Gianola, 2013), which are basically based on the framework of linear regression. Generally, these parametric methods assume that the effect of each marker is independently distributed with a specific prior distribution given by corresponding statistical methods. Clearly, such an assumption of independent distribution for each SNP effect is statistically inappropriate, especially when the adjacent markers are in high linkage disequilibrium (LD) with the same causal gene. This unrealistic assumption potentially sacrifices the prediction accuracy to some extent. To address this issue, Yang and Tempelman (2012) proposed a first-order antedependence model to account for the nonstationary correlations between SNP markers through assuming a linear relationship between the effects of adjacent markers. As expected, the proposed antedependence-based WGP models outperformed their conventional counterparts in the prediction accuracy of genomic merit in the context of single-trait analyses.

Compared with single-trait analyses, multiple-trait joint analyses have been widely confirmed as having obvious advantages in statistical power and parameter estimation accuracy in earlier studies of quantitative trait locus (QTL) mapping and GWAS (Liu ; Liu ; Jiang and Zeng, 1995; Korte ). Following this evidence, it can be speculated that in WGP fields, joint prediction of multiple traits should benefit from genetic correlation between traits and thus have the potential to obtain higher prediction accuracy than single-trait methods, especially for a low-heritability trait that is genetically correlated with a high-heritability trait. Statistical methods for joint prediction of multiple traits have received close attention in recent years. Particularly, Bayes A, Bayes C and Bayes Cπ have been extended to be applicable for multiple trait analyses (Calus and Veerkamp, 2011; Jia and Jannink, 2012), demonstrating clear advantages over single-trait methods.

Logically, a more natural strategy for enhancing WGP is to take into consideration relations between multiple traits, as well as those between SNP effects simultaneously for achieving high prediction accuracy from the viewpoint of statistics. So far, there is still a gap between multiple-trait joint prediction approaches and models considering correlated SNP effects, for example, antedependence-based WGP models. Therefore, aiming at developing an improved prediction methodology, we are in an attempt to construct a multiple-trait antedependence-based WGP model to bridge such a gap aforementioned under the framework of Bayesian algorithm. Specifically, we proposed a novel multivariate method with two different types of prediction models via setting the antedependence parameter as either a matrix or a scalar, respectively. Theoretically, our proposed models can relax the conventional assumption of independence of marker effects while simultaneously taking advantage of the correlation between traits by modeling a linear relationship of effect vector between each pair of adjacent markers. Using simulations as well as the publicly available data sets including the 16th QTL-MAS workshop data and the heterogeneous stock mice real data, we compared our proposed method with the classical approaches including the Bayes A version by Habier and multi-trait Bayes A by Jia and Jannink (2012) to further validate the performance of our proposed method. Our study clearly demonstrated that the proposed multiple-trait antedependence-based WGP method is more accurate and robust than corresponding traditional counterparts for genomic prediction, offering a feasible way to jointly predict complex traits in human genetic epidemiology, as well as plant and livestock breeding.

Materials and methods

Bayesian multivariate antedependence model

We developed a Bayesian version of first-order multivariate antedependence model (Zimmerman and Nunez-Anton, 2010). The performance of proposed method herein was evaluated by comparison with traditional single-trait Bayes A (Meuwissen ; Habier ), as well as the multi-trait Bayes A method initially developed by Jia and Jannink (2012) (see Supplementary Method, Section 1 for details). Considering different types of correlation between SNP effects, we developed two forms of first-order multivariate antedependence models considering the antedependence parameter (Zimmerman and Nunez-Anton, 2010) as a matrix as well as a scalar, named as ‘matrix model' and ‘scalar model', respectively.

Matrix model

Assuming n individuals were genotyped at p SNP markers, the matrix model for joint prediction of m traits is expressed as,

where y is an m-element phenotypic vector for individual i (i=1,…,n); is the vector of overall population mean of m traits; is an m-element vector for the effects of the jth SNP marker on all m traits (j=1,…,p); Z is the SNP genotype code for individual i at marker j; and Σ is an m × m covariance matrix of the residual effects. In this study, all markers are considered to be biallelic, and marker genotypes were coded as 0, 1 or 2 corresponding to the number of copies of an allele at a locus.

As shown in the model (1), each marker effect is assumed to have a linear relationship with that of the preceding adjacent marker based on the physical position or the genetic position along a chromosome, that is, = + . It should be pointed out here that in the traditional univariate model, the prior of antedependence parameter t generally follows a normal distribution (Yang and Tempelman, 2012). For the case of multivariate model, we accordingly adopted a random matrix herein with the prior of normal distribution. It has been demonstrated elsewhere (Minka, 2001) how to use such a matrix normal prior in similar analyses. Following this initial work, we constructed corresponding prior normal distribution for the matrix antedependence parameter in our proposed methodology. Here is a m × m matrix corresponding to the jth maker, following a matrix normal distribution with mean matrix M (m × m), among-row covariance and among-column covariance . δ is an m-element vector for the residual part of the jth marker's effects, following a multivariate normal distribution with mean zero and covariance . Note that is the among-row covariance of the prior distribution of as well as the covariance of the prior distribution of . Introducing to the prior of herein can clearly facilitate the implementation of Gibbs sampling. The reason is that we can readily get a close form of the conditional posterior of under such case (see Page 7 in Supplementary Method). As we know, an explicit expression of posterior probability density function (pdf) is a prerequisite for performing Gibbs sampling. Further, we specified for an inverse Wishart prior distribution with an identity scale matrix and m degrees of freedom. Accordingly, we assigned an inverse Wishart prior distribution with scale matrix Ψ and v degrees of freedom for matrix Σ and an inverse Wishart prior distribution with scale matrix Ψ and v degrees of freedom for (j=1,..,p). Both Ψ and Ψ were further assumed to follow a Wishart distribution with an identity scale matrix and m degrees of freedom.

In our application, we specified the hyperparameters in model (1) as . The matrix model was solved through Gibbs sampling. We give the details on the Gibbs sampling procedures in Supplementary Method, Section 2. In total, 200 000 MCMC iterations were conducted for the matrix model and the first 100 000 iterations were discarded as burn in. Every 10th sample was kept for the follow-up inference.

Scalar model

Basically, the scalar model is a simplified version of the matrix model. The antedependence parameter in the scalar model was specified as a scalar rather than a matrix, that is,

where t is the scalar antedependence parameter of on . Using the scalar means that the antedependence parameter at a marker is identical for all traits. The simplified model can decrease the sampling complexity in MCMC cycles and still reflect the correlation between each pair of adjacent markers. Further, we set t a normal prior as

where was further given a scaled inverse χ2 prior distribution as

The hyperparameter, , is given a Gamma (shape=1, rate=1) prior and need to be estimated in the model. All of other underlying parameters, including (j=1,..,p), (j=1,..,p), Ψ, Σ and Ψ, were separately given a prior as the same as in model (1).

We herein specified the hyperparameters in the scalar model as . The scalar model was solved through Gibbs sampling similar to the matrix model, and the details on the Gibbs sampling procedures were given in Supplementary Method, Section 3.

Prediction of multiple traits

We used a training population to predict the effect of each marker, (j=1,…,p), and then estimated the total genetic values of m traits for an individual as , where g is an m-element vector of additive genetic values of m traits for the ith individual. Prediction accuracy is defined as the correlation between the estimated total genetic values and known true total genetic values in the validation population. If the true total genetic values are unknown in the validation population, we used predictive ability to compare the various methods. Predictive ability is defined as the correlation between the estimated total genetic values and their corresponding phenotypes in the validation population.

To validate the performance of the proposed methods, we conducted extensive simulation and the analyses on the 16th QTL-MAS workshop data set and the heterogeneous stock mice real data to compare it with the single-trait and multi-trait Bayes A models with respect to prediction accuracy or predictive ability.

Simulation analyses

We conducted extensive simulation by considering different levels of genetic and phenotypic correlations between traits, SNP densities and the numbers of underlying QTLs to validate the prediction accuracy and robustness of our proposed models.

Data simulation

A whole-genome simulation program, GPOPSIM (Zhang ), which is based on the mutation-drift equilibrium model, was used to generate biallelic markers and QTLs. The simulation started with a base population of 100 individuals (50 unrelated males and 50 unrelated females), followed by 5000 generations of random mating with the same population size, to generate data with realistic LD patterns caused by mutation and drift. The entire genome was composed of one chromosome of length 1 Morgan. All of 20 001 potential SNP markers were equally spaced on the chromosome with a QTL placed directly in the middle of each interval of adjacent markers. The mutation rate for both SNP markers and QTLs was specified to be 1.0E-4 per locus per generation. In generation 5001, the population size was expanded to 500 by randomly mating 50 males with 50 females from generation 5000. Generation 5002 was generated by randomly mating 50 males with 50 females from generation 5001 and also had a population size of 500. Similar to the study by Yang and Tempelman (2012), we randomly selected 30 QTLs with a minor allele frequency (MAF) >0.05 to generate true breeding values for individuals in generations 5001 and 5002. Following Jia and Jannink (2012), QTL effects on two phenotypic traits were sampled from a standard bivariate normal distribution with correlation 0.5, which assumes a certain level of pleiotropy effects at all loci. The true breeding value for each individual was the genotype-based sum of the QTL effects for each trait. The total genetic covariance matrix for the two traits was determined as , where p and is the MAF and the sampled vector of effect at QTL k, respectively. Normal error deviates were added to achieve heritabilities of 0.5 for trait 1 and 0.1 for trait 2. All individuals have phenotypes on both traits. The covariance of errors between traits was set to zero.

We considered the above simulation scenario as the default validation data set. Accordingly, we perturbed a single-simulation parameter at a time from the default scenario to generate simulated data under other scenarios for comparison. Perturbed parameters for simulations included genetic correlation between traits (0.2, 0.5 and 0.8), error correlation (−0.2, 0 and 0.2), and number of selected QTLs (30 vs 300). Each of these simulation scenarios was repeated 30 times for producing convincing results. For all simulation scenarios, generation 5001 was used as the training population and generation 5002 as the validation population.

Estimation of prediction accuracy in simulated data

For each simulated data set, we filtered the 20 001 SNP markers using a criterion of MAF>0.05, leaving ~4100 markers. The filtered SNPs were sorted merely according to the physical position in base pair and then used to predict genomic breeding values for the validation population. The prediction accuracy for each data set was calculated as the correlation between the estimated genomic breeding values and the known true genomic breeding values in the validation population. In each scenario, we used a paired t-test to test if there was a significant difference of prediction accuracy between our proposed models and the conventional models based on the 30 replicates.

Furthermore, based on the default scenario, we also employed the subsets of SNPs through sampling every 10th and 25th SNPs from the full set of SNPs respectively to predict genomic breeding values for testing the effect of LD between adjacent markers on the prediction accuracy.

Analysis of the 16th QTL-MAS workshop data set

In addition to our simulated data sets, we analyzed the 16th QTL-MAS workshop data set (http://qtl-mas-2012.kassiopeagroup.com/en/index.php) to further validate the prediction accuracy of our proposed methods.

The workshop data were generated to investigate performance of various WGP approaches when dealing with multiple correlated traits. Specifically, a base population of 1020 unrelated individuals (20 males and 1000 females) was generated with a 499.750 Mb long genome consisting of five chromosomes. Each chromosome had a size of 99.95 Mb and carried 2000 equally distributed SNP. A total of 50 QTLs were generated, and each QTL had an effect on at least two traits. The detailed information of the QTLs can be found at http://qtl-mas-2012.kassiopeagroup.com/en/dataset.php. Each of four generations (G1–G4) consisted of 20 males and 1000 females and was generated from the previous generation by randomly mating each male with 51 females. The female prolificacy was set to 1 except for the 20 dams of male that generated two offspring (1 male and 1 female). No Generation overlapped in the process of data simulation. The mimic phenotypes consisted of three genetically correlated milk production traits.

In the analyses, we removed all SNP markers with MAF=0, leaving 9969 markers for further analyses. The SNP markers had been sorted by physical position. The 3000 females from G1–G3 had observations for all three traits, which were used as the training population. The remaining 1020 individuals from G4 with known total genetic values were used as the validation population. The prediction accuracies obtained by the different multi-trait methods were compared with officially reported results (http://qtl-mas-2012.kassiopeagroup.com/en/index.php). Briefly, four different groups independently submitted their results, and all the 14 analytical methods were merely focused on single-trait prediction, including GBLUP, the Bayesian alphabet and various regularized regression methods. The highest prediction accuracies obtained by the methods were 0.794, 0.853 and 0.828 for traits 1, 2 and 3, respectively, which were used as the results of best single-trait methods and as a benchmark for comparison with the results herein.

Real data analyses on heterogeneous stock mice data set

The heterogeneous stock mice data set (http://mus.well.ox.ac.uk/GSCAN) consists of 1940 genotyped mice generated based on intercross mating among eight inbred strains (Valdar ). We randomly selected two pairs of immunological traits, %CD4+/CD3+ and %CD8+, and %CD4+ and %CD8+, to compare prediction accuracy between the methods used in our study. Phenotypes for these traits had been preadjusted for marginally significant covariates.

Following the study by Speed and Balding (2014), SNPs were removed if they had MAF<0.01, HWE P<10-4, or call rate <0.99. Missing SNP genotypes were imputed on the basis of their corresponding allelic frequencies in the data set following Legarra . Individuals were removed unless they had phenotypes of all the selected traits. After quality control, 1404 individuals having both phenotypes and genotypes were kept, and each individual had genotypes of 9159 SNPs mapped to autosomes. In addition, we ordered SNP genotypes on the basis of their genetic positions along the chromosome.

We randomly divided individuals into two nearly equal-sized subsets as training and validation sets, respectively. This was replicated 20 times for performance comparison among different prediction methods.

Results

Results from the simulations

Prediction accuracy under different LD levels of adjacent markers

We investigated the performance of the proposed methods with default simulation parameters while dealing with SNPs with different densities to explore the effect of LD extent on prediction accuracy. Table 1 presented the three average LD levels (r2 = 0.33, 0.22 and 0.14) of adjacent markers corresponding to the full set of markers and two subsets of markers sampled from the full set of markers with intervals of 10 and 25, respectively.

Table 1

Number of used markers and LD level across 30 replicates under varied scenarios

Scenario no.	Genetic correlation	Error correlation	No. of QTLs	Way of using markers	Average no. of used markers (std.)	Average LD level of adjacent markers (std.)a
1 (default)	0.5	0	30	All	4119 (159)	0.333 (0.016)
				Every 10th	411 (16)	0.220 (0.018)
				Every 25th	164 (6)	0.136 (0.015)
2	0.5	0	300	All	4140 (126)	0.333 (0.015)
3	0.2	0	30	All	4044 (162)	0.336 (0.013)
4	0.5	−0.2	30	All	4030 (164)	0.338 (0.015)
5	0.5	0.2	30	All	4071 (141)	0.334 (0.013)

Abbreviation: QTLs, quantitative trait loci.

Average=sum(Ai)/30 and std.=sqrt(sum((Ai-Average)2)/29), where Ai is the average LD level in the ith replicate.

As clearly shown in Figure 1, our proposed methods always exhibited the best performance among all methods across different LD levels. Especially, with the increase of average LD, our methods were becoming more advantageous over other analytical approaches for the low-heritability trait 2.

Figure 1

Prediction accuracies of various models for the high-heritability trait (a) and the low-heritability trait (b) under different LD levels of adjacent markers.

Prediction accuracy under scenarios with a varied number of QTLs

Besides the default 30-QTL scenario, we considered a 300-QTL scenario for further comparison among different analytical methods. Despite the clear superiority of the proposed methods under the 30-QTL scenario, there were no significant differences of prediction accuracy between our proposed models and the traditional multi-trait Bayes A model for the 300-QTL scenario (Figure 2).

Figure 2

Prediction accuracies of various models for the high-heritability trait (a) and the low-heritability trait (b) under scenarios with a varied number of underlying QTLs.

Prediction accuracy under scenarios with varied genetic correlations

It can be clearly shown from Figure 3 that our proposed methods were always among the top rank on the prediction accuracy compared with other analytical models under different levels of genetic correlations. Furthermore, all multi-trait prediction methods significantly outperformed the single-trait analyses on the traits investigated under all scenarios.

Figure 3

Prediction accuracies of various models for the high-heritability trait (a) and the low-heritability trait (b) under scenarios with varied genetic correlations between traits.

It is worth to note that the gain in prediction accuracies of multi-trait over single-trait models were becoming more significant for the low-heritability trait 2 as genetic correlation increased between traits. In contrast, for the high-heritability trait 1, no obvious change in prediction accuracy was observed as the genetic correlation increased from 0.2 to 0.8. The similar trend was also observed in previous studies (Calus and Veerkamp, 2011; Jia and Jannink, 2012). This implied that the superiority of our proposed models over the multi-trait Bayes A model resulted from the feature that the model accounted for the dependence on the effects of adjacent markers.

Prediction accuracy under scenarios with varied error correlations between traits

Besides the default scenario with error correlation of 0.0, we also considered another two scenarios with varied error correlations of −0.2 and 0.2, respectively. From Figure 4, for both low- and high-heritability traits, the prediction accuracies were largely consistent for each individual analytical across the three scenarios. Similar profiles as those aforementioned, our proposed methods still achieved highest prediction accuracies across various scenarios, and all the multi-trait prediction methods rendered obvious advantages over the single-trait approach.

Figure 4

Prediction accuracies of various models for the high-heritability trait (a) and the low-heritability trait (b) under scenarios with varied error correlations between traits.

Results from the analysis of the 16th QTL-MAS workshop data set

We also used the 16th QTL-MAS workshop data set to validate our proposed methods. As shown in Figure 5, multi-trait methods, including the multi-trait Bayes A model, the scalar model and the matrix model, showed considerable advantage over single-trait methods.

Figure 5

Prediction accuracies of various models for the 16th QTL-MAS workshop data set. aMeans the best single-trait method officially reported in the 16th QTL-MAS workshop.

It can be found clearly that our proposed scalar model obtained the best results with regard to the 16th QTL-MAS workshop data set. The prediction accuracies of the scalar model for the three traits increased by as high as 8.6%, 2.8% and 5.5%, respectively, compared with the reported best single-trait methods.

Results from the analysis of real heterogeneous stock mice data set

As shown in Table 2, for both of trait pairs, our proposed models always ranked top two in prediction abilities among all four methods for each of traits, and all three joint prediction methods outperformed the single-trait model significantly (P<0.01).

Table 2

Predictive ability of various methods for the heterogeneous stock mice data across 20 replicates

Trait pair	Trait	Predictive ability of various methods
		Single-trait Bayes A	Multitrait Bayes A	Scalar model	Matrix model
1	%CD4+/CD3+	0.609	0.616	0.620	0.625
	%CD8+	0.663	0.664	0.668	0.671
2	%CD4+	0.358	0.372	0.375	0.375
	%CD8+	0.665	0.667	0.670	0.670

Similar as simulation analyses, our antedependence models obtained predictive abilities significantly higher than the multi-trait Bayes A for both of trait pairs. It is worth noting herein that the proposed matrix model had a better performance than the scalar model for each of trait (P<0.01) in the first trait pair, exhibiting its potential superiority in practice in contrast to its simplified version, the scalar model.

Discussion

In WGP analyses, the inferred adjacent SNP effects should be spatially correlated due to chromosomally proximal effects of potential QTL (Yang and Tempelman, 2012), which is most likely in LD with the SNPs surrounding it, even when no biological mechanism exists among these SNPs. The antedependence model can model such nonstationary correlated effects via incorporating antedependence parameters between adjacent SNPs into the model. This is the theoretical basis that the antedependence model can work better than the traditional independence models given no additional information about markers but genotypes.

In the current study, we firstly developed two types of multi-trait antedepedence models for joint phenotype prediction. We hypothesized that our proposed WGP methods could enhance the prediction performance owing to considering information of correlation among traits and adjacent SNPs simultaneously. This theoretical superiority was further validated herein through extensive simulation and real-data analyses, and the resulting prediction patterns under various scenarios consistently revealed the significant performance gain of our proposed methods compared with the corresponding counterpart methods. Our proposed WGP methods offer a feasible alternative in joint prediction of multiple complex traits in human genetic epidemiology as well as livestock breeding. We have developed C programs to implement the corresponding prediction models, which are freely available upon request.

In our methods, we employed a matrix and a scalar antedependence parameters, respectively, for our proposed multivariate prediction method. The introduced antedepedence parameter matrix considered respective effect correlation between pair of adjacent markers for each of correlated traits, which should outperform the scalar model in theory. However, our simulation data analyses and the analysis of the 16th QTL-MAS data set have not demonstrated obvious advantage of the matrix model over the scalar model. This may be owing to the extreme case that was mimicked in our simulation scenarios as well as the 16th QTL-MAS data set. Specifically, in our simulation, all simulated QTL effects were assumed being sampled from the same distribution and contributing to both traits simultaneously. In the 16th QTL-MAS data set, each simulated QTL contributed to at least two traits out of all three. In such case, the antedependence parameters for each of correlated traits tend to be equal. Accordingly, the scalar model suits the data sets with the such genetic structures better than the matrix model and was prone to obtaining a higher prediction accuracy.

In the simulation scenario when considering all 300 QTLs with very small effects, the antedependence parameters are inefficient to model the dependence between SNPs, lowering its advantage over the traditional methods. Another aspect in our simulation is that a small base population of 100 followed by a large number of generations was simulated. Although this scheme was usually adopted in previous studies (Zhang ; Calus and Veerkamp, 2011; Yang and Tempelman, 2012), it would generate data with extremely small effective population size such that the simulated data could not represent the population structure of real data.

Owing to limited simulation scale in the studies, it is not feasible to simulate various scenarios with genetic structures consistent of all possible built-in mechanisms for performance validation. Hence we turned to further perform real-data analyses on two sets of complex traits from a publicly available heterogeneous stock mice data to explore the feasibility of our methods. The results clearly demonstrated that the proposed multivariate antedependence model outperformed corresponding traditional counterparts, and the matrix model has also shown the advantage over its simplified version, the scalar model, as we expected.

Yang and Tempelman (2012) have reported that the antedependence methods outperformed their corresponding classical counterparts in single-trait prediction, and the increase in prediction accuracy contributed by the single-trait antedependence model was generally <0.05. This increase in prediction accuracy may be not comparable to that for low-heritability traits benefiting from joint prediction of multiple highly correlated traits which can be usually >0.10 (Jia and Jannink, 2012). This has been clearly reflected in our analyses. As shown in Figures 1,2,3,4, the low-heritability trait generally benefited more from joint prediction of multiple traits than the high-heritability trait in various scenarios.

A different aspect in our prediction model from that by Yang and Tempelman (2012) is regarding the expectation of the prior of the antedependence parameter. In the study by Yang and Tempelman (2012), a marker-specific antedependence parameter was assumed a normal distributed prior with unknown expectation which needs to be estimated in the model. In contrast, we directly set the prior expectation as zero. Our consideration for this is that the sign of the additive effect at any marker is merely determined by the way we code its genotypes (for example, genotypes BB can be coded as 0 or 2) and the antedependence parameter should have a sign determined by the signs of SNPs effects therein disregarding the residual term of Equation (2) aforementioned. Thus, as the genotypes of biallelic markers are coded as the number of copies of one arbitrary allele in our study, it is equally probable for the antedependence parameter to be positive or negative, and then it is reasonable to assume its expectation as zero.

It should be pointed out herein that although we focused on the use of common variants in current study, rare variants can also be readily incorporated into our proposed method by drawing idea from the collapsing approach (Morris and Zeggini, 2010). It will be imperative to incorporate rare variants into prediction in future for the situation where resequence data instead of traditional SNP chips are widely used in WGP.

Furthermore, our developed antedependence version of the multi-trait Bayes A models can be modified to other types of Bayes model, for example, Bayes B and Bayes Cπ (Habier ), which have been two popular Bayesian variable selection methods for single-trait WGP. For example, a multi-trait antedependence version of Bayes B can be developed based on model (2) aforementioned with a mixture prior distribution of (j=1,..,p) as . Then the conditional posterior distribution of each unknown parameter can be derived drawing idea from the single-trait antedependence-based Bayes B (Yang and Tempelman, 2012).

Aside from predicting yet-to-be observed phenotypes, our proposed models are also useful in QTL-mapping studies according to previous reports on applying WGP methodology to GWAS (Peters ; Garrick and Fernando, 2013). By adding up genetic variance contributed by SNPs in a SNP window, we can calculate the genetic variance of consecutive SNP windows along the genome and then consider windows contributing relatively large genetic variance to be QTLs (Peters ). Bootstrap analysis can be further used to determine the significance of detected QTLs as described by Peters .

We believed that it is necessary to find a way to integrate multiple advantageous strategies for developing an optimal WGP model, for example, joint prediction of multiple traits, the antedependence model, incorporating dominant and epistatic effects (for example, Wang ; Nishio and Satoh, 2014), modeling genotype × environment interaction (for example, Crossa, 2012). As these strategies are utilizing different principles, it is anticipated that the strength contributed by each of the strategies can be accumulated if these various strategies are integrated in a proper way. Our study can be considered as an example of such integration, which takes advantage of both antedependence and multivariate models for joint prediction.

Data archiving

All simulated data analyzed in the present study have been deposited in Dryad (http://doi.org/10.5061/dryad.dd60v).