From quantitative trait locus mapping to genomic selection: the roadmap towards a systematic genetics

Traditional QTL mapping methods

Techniques for QTL mapping include single marker mapping, interval mapping, multiple loci mapping as well as composite interval mapping. Principles in these mapping techniques are generally the same and methods used in one population can be extended to other experiment populations. Single marker test examines the segregation of quantitative or qualitative traitswith respect to the examining genotype at a single locus, while multiple loci mapping considers multiple makers, possible high order maker interactions as well as environment factors simultaneously. The interval mapping^[7] and composite interval mapping^[8] are extensions of single marker test and multiple loci test, respectively. In the interval mapping, consecutive two testing markers (or several markers in a testing window in composite interval mapping) are ordered according to their physical location, and those peak testing values in single tests are declared as QTLs^[7]. The various techniques test the genotype and phenotype association with different genetic effects such as additive and dominance effect. In regression models, they can be examined simultaneously by adding dummy variables to encode different effects. A widely used genetic model is the Cockerham model^[9] that defines the values of the additive effect as -1, 0 and 1 for the three genotypes and the values of the dominance effect as -0.5 and 0.5 for homozygotes and heterozygotes, respectively.

Similar to QTL mapping, the goal of association study in natural population is to identify SNPs in individuals that are systematically associated with different disease states. Using the natural occurring DNA variations as markers to trace inheritance in families are similar to QTL mapping, while extra cares are required to handlepopulation structures in genome-wide association study ofcommon human diseases. In the ensuing sections, we will focus on the general computational methods that have been applied to both problems.

Multiple-variant methods

The single variant approach is simple and straightforward for most models such as t-test, analysis of variance and simple linear regression for quantitative traits, the Cochran-Armitage test, Pearson χ² test and generalised regression for qualitative traits. However, single marker methods require multiple test correction such as Bonferroni correction and false discovery rate^[10], to control the overall type I error rate. Given the large number of possible markers, multiple test correction leads to such stringent criterion that most modest effects could not pass the threshold. Moreover, since complex phenotype is controlled by multiple genetic factors and their interactions with environmental effects, marker identified with a single marker method may only explain a small proportion of the phenotype variation.

Providing the promising of high power and reasonable type I error rate, multiple-variant approach needs to take care of several challenges. Firstly, traditional ordinary least square method fails for the case of p >> n. When taking into account of interactions between loci, variable selection is even more demanding. Secondly, it is likely that the increased degree of freedom in a multiple-variant model jeopardises the power gain over single-locus models. Thirdly, in binary-trait or case-control association analysis, complete or semi-complete separation can be a serious problem due to the discrete nature of both marker data and binary outcomes^[11]. Accurate variable selection and sparse model inference methods become critical for the multiple-variant approach, especially when epistatic effects are considered. Variable selection and shrinkage are two techniques in handling the selection problems, and they typically produce sparse models.

Variable selection methods include forward selection, backward elimination and forward stage-wise selection. However, these greedy selection methods may result in suboptimal subset, and are computational expensive even for small number of variables^[12]. On the other hand, variable shrinkage method includes all variables in the model and applies a penalty function or appropriate prior distributions on the variables to automatically shrink most non-effects towards zero. Moreover, including all possible effects into a single model results in whole-genome QTL mapping, which overcomes limitations of the genetic model considered by traditional mapping methods and prevents all problems of model selection^[13].

Whole-genome QTL mapping and genomic selection

Strictly speaking, whole-genome QTL mapping includes additive and dominance main effects and all pair-wise interactions. For a QTL model includes k markers, the total number of effects is $p=2k+4k(k-1)$ . However, in practice,most multiple QTL models cannot handle such large number of effects and epistaticeffects are usually not considered (which is named as genome-wide QTLmapping)^[14]. State-of-the-art whole-genome QTL model can handle up to 10⁸ effects including all main and pairwise interaction effects on a personal computer^[15]. Nevertheless, next generation sequencing technologies can easily genotype up to 10⁷ SNPs for most model organisms and whole-genome QTL mapping at this scale is computationally too demanding and causes memory overflow. In such cases, markers carrying duplicated information content need to be deleted^[16].

Another natural outcome of whole-genome QTL mapping is GS, which is originally proposed in genome-wide QTL mapping that mapping results enable predictions of estimated breeding values (EBV)^[14].GS is rooted from marker-assisted breeding, and can be viewed as a new generation of breeding method different from traditional breeding that relies on phenotypic selection and relative information^[17]. In GS, individual genetic merits are estimated by simultaneously accounting all markers and all types of marker effects. From biological point of view, a marker map must cover all genomic positions such that all QTLs can be assessed for their contributions to EBV, although in reality some adjustment for linkage disequilibriummight be required to avoid colinearity. From computational point of view, the QTL model must simultaneously estimate all genetic effects including main and epistatic effects, genetic and environmental interactions and effects of rare alleles, such that the breeding values are predicted based on all significant effects^[18]. Whole-genome QTL mapping provides more accurate marker effects and phenotypic variance estimation, which results in better performance of predicting genomic EBV.

Variable shrinkage methods

Two well-known methods with ability in handling large number of variables are regularised regression and Bayesian shrinkage method^[19]. Consider the general multiple linear regression problem with n samples and p variables:

$$y=\text{X}\beta +\epsilon ,(1)$$

where y is an n × 1 vector collecting n individual phenotype, X is an n × p designed matrix for p marker effects, ε is the residual error that follows a normal distribution with zero-mean and covariance ${\sigma }_0^2\text{I}:\epsilon ~N(0,{\sigma }_0^2\text{I})$ . The generic penalised regression method infers the model parameters as:

$$\widehat{\beta }={\arg }_{\beta }\min {\Vert y-\text{X}\beta \Vert }_2^2+\lambda {\Vert \eta (\beta )\Vert }_q^q,(2)$$

where η(β)≥0, and ${\Vert \eta \left(\beta \right)\Vert }_q^q={\displaystyle \sum _{j=1}^p\eta {\left({\beta }_j\right)}^q}$ is the l_q penalty term. Here, the intercept is absorbed into matrix X, and it is not penalised. Without loss of generality, variables are assumed to be re-scaled to have variance 1 to obtain invariant penalties. Notice that when $\text{that}\text{when}\eta ({\beta }_j)=\{\begin{array}{l}1,\text{if}{\beta }_j\ne 0\\ 0,\text{otherwise}\end{array},\forall j,$ the model penalises the number of effects,and the solution of (2) is equivalent to the AIC (for λ = 2) or BIC (for λ = ln(n)) in forward selection and backward elimination, and the estimates are the same as the least square solutions. From this point of view, variable shrinkage is a generalisation of variable selection. Let $\eta \left({\beta }_j\right)=\left|{\beta }_j\right|,j=1,2,\cdots ,p$ to penalise the effect size, (2) becomes the bridge regression, which is a generalisation of two most popular penalised regression methods, namely the l₂ penalty of the ridge regression^[20], and the l₁ penalty of the least absolute shrinkage and selection operator (Lasso) that implicitly enforces a sparse solution^[19].

In the Bayesian shrinkage approach, a prior distribution $p(\beta |\lambda )$ is assigned to β. The posterior distribution given the observed data takes the form of

$$p(\beta |y,\lambda )\propto p(y|\beta )p(\beta |\lambda ).(3)$$

If we are looking for β that maximises the posterior distribution, the problem is equivalent to finding

$$\widehat{\beta }={\arg }_{\beta }\max \{\log p(y|\beta )+{\displaystyle \sum _{j=1}^p\log p({\beta }_j|\lambda )}\},(4)$$

which is the penalised ML method. Then $\widehat{\beta }$ is the mode of the posterior distribution and this method is referred as maximum a posteriori (MAP) approach.

The Bayesian interpretation of penalised regression methods has sparked interests in developing Bayesian hierarchical regression models, which can be called as a Bayesian Lasso approach. A direct gain of Bayesian Lassois to incorporate variance information to facilitate significance test, and different techniques are employed for different designs to achieve sparse estimation. Theoretically, prior distributions with spike finite limit at zero and flat tails at two ends can be penalty of the log posterior distribution (Figure 1). In practice, conjugate priors are often chosen. Among them the Normal + inv-χ² (t-family) distribution^[21] and NEG hierarchical priors^[22,23,24] are most oftenly used with both variable shrinkage and computational feasibility considerations. Table 1 listed several popular Bayesian Lasso methods designed.

Figure 1: Prior distributions that penalise posterior distributions.

Table 1 Popular Bayesian Lasso methods

For Bayesian estimation, Markov Chain Monte Carlo (MCMC) can be employed to draw samples from posterior distribution for each parameter. However, for high dimensional data the MCMC method is known to be computationally intensive. The MAP is a modal representation of the posterior distribution that achieves faster computation and easier interpretation. Efficient methods for MAP estimation have been developed, which further integrate out the variance components and employ numerical methods toobtain similar sparse results as Lasso (named as HyperLasso)^[22]. Another method stands in between MCMC and MAP is the expectation-maximization (EM) algorithm that estimates posterior mean and variance simultaneously through iterative expectation (E-step) and maximization (M-step). However, convergence become a serious problem when model dimension increases^[25]. Recently, more efficient algorithm that does not rely on MCMC yet infers the posterior distribution is achieved by empirical Bayesian Lasso (EBlasso) method^[23,24].

Different from the iterative EM algorithm, EBlasso first finds the marginal posterior distribution of the variance components. Due to the shrinkage applied from the prior distribution on the variance components, most variables will have zero variance that maximises the marginal posterior distribution, and only those variables with none-zero variance will stay in the model, result in a sparse presentation. Next, the posterior means of the non-zero variables are estimated with the give variance. Along with other algorithmic techniques, the EBlasso approach is very efficientand is able to handle a model with several million of variables^[26]. Our previous studies demonstrated that EBlasso outperformed several other multiple QTL mapping methods including the empirical Bayes method in^[21], the Bayesian hierarchical generalised linear models(BhGLM)^[25], HyperLasso^[22] and Lasso^[19]. EBlasso has also been applied to whole-genome QTL mapping^[15] and can be easily extended to GS.