Integrated Framework for Selection of Additive and Nonadditive Genetic Markers for Genomic Selection

1. Introduction

Genomic selection (GS) was proposed by Meuwissen et al. (2001) as an improved variant of marker-assisted selection for efficient breeding of animals and plants. Marker information covering whole genome is used to estimate the breeding value in GS, even without information of exact gene location. With the advancement of next-generation sequencing technology, it is now easier to genotype single-nucleotide polymorphisms (SNPs) within a short period of time, which is being used for GS of superior parents based on genomic-estimated breeding value (GEBV). To estimate GEBV, phenotypic and genotypic information of the reference population is collected to build a prediction model. Initially, all the individuals in the reference population are genotyped to get SNPs from the entire genome. The precision of prediction model increases by increasing the sample size of reference population. The collected information on genotype and phenotype from reference population is used to train a predictive model for estimation of GEBV. This developed model is being used to obtain GEBV of an individual from the breeding population based on marker information.

Many models are available in literature for selection of genetic markers. Most of them are developed to capture additive genetic effects. However, some models are also available, which can be used for modeling nonadditive genetic effects, that is, epistasis. In general, it has been observed that these models are able to efficiently capture either the additive or nonadditive genetic effects. Selection of appropriate genetic markers is an important step in GS for estimation of GEBV. The selection of most relevant markers depends on adopting the most efficient feature selection method and the estimation of GEBV for trait of interest. There is hardly any model in literature that is capable to capture both additive and nonadditive genetic effects efficiently. Therefore, there is a strong need to appropriately develop an integrated framework for selection of the most appropriate additive and nonadditive genetic markers for estimation of GEBV. Therefore, first, we have evaluated statistical performance of important models for selection of additive effect markers, namely, linear least-squared regression (Meuwissen et al., 2001), least absolute shrinkage and selection operator (LASSO; Tibshirani, 1996), ridge regression (Hoerl and Kennard, 1970a,b), BLUP (best linear unbiased prediction) (Henderson, 1975), sparse additive models (SpAM; Ravikumar et al., 2009) along with two important procedures for selection of nonadditive genetic markers, namely, mRMR (minimum redundancy maximum relevance) (Peng et al., 2005) and Hilbert–Schmidt Independence Criterion (HSIC) LASSO (Yamada et al., 2014) through a simulation study.

In this empirical study, it was observed that SpAM performs very well for additive genetic architecture and HSIC LASSO is comparatively efficient in case of an epistatic effect in data for selection of genetic markers. It may be noted that, in the presence of an epistatic effect, the performance of SpAM is very poor, whereas for additive genetic effect, the performance of HSIC LASSO is encouraging. Keeping in view these facts, it is pertinent to develop a robust framework for selection of genetic markers for estimation of GEBV for GS, which can account for both additive and epistatic genetic effects in selection of parents from the breeding population. In this article, an attempt has been made to propose an integrated framework for selection of genetic markers to estimate GEBV. The comparative performance of the proposed integrated framework has been examined through a simulation study.

2. Methods

An integrated framework for selection of genetic markers to estimate GEBV has been developed for GS by combining selected genetic markers from an additive model, that is, SpAM and a nonadditive model, that is, HSIC LASSO as these models are found to be superior for selection of additive and nonadditive genetic markers, respectively. In the first instance, a brief description of SpAM and HSIC LASSO models is given below.

2.1. Sparse additive models

The SpAM is useful for high-dimensional feature selection (Liu et al., 2009; Ravikumar et al., 2009; Raskutti et al., 2012; Suzuki and Sugiyama, 2013). The SpAM optimization problem can be expressed as follows:

where n is the total number of individuals, p is the total number of variables (i.e., markers), y is the vector of phenotypic value for all the individuals, are regression coefficient vectors, is gram matrix, and , are the elements of the incidence matrix, is kernel function, is a coefficient for , and λ > 0 is a regularization parameter. SpAM is similar to the hierarchical multiple kernel learning (Bach, 2009), which implements an SpAM with an alternative sparsity-inducing regularization. SpAM is a convex method and can be efficiently optimized by the back-fitting algorithm.

A disadvantage of SpAM is that it can only deal with additive effects. In case of epistatic effects in the data, SpAM may fail to select significant markers. Another disadvantage is that the output y should be a real number in SpAM, which indicates that SpAM cannot deal with structured outputs such as multilabel and graph data. Moreover, SpAM is a computationally expensive procedure.

2.2. HSIC least absolute shrinkage and selection operator

A feature-wise nonlinear LASSO of the following form was proposed by Yamada et al. (2013), which is called as HSIC (Gretton et al., 2005) LASSO.

where n is the total number of individuals, p is the total number of variables (markers), ‖ · ‖_Frob is the Frobenius norm, and are centered gram matrices, and are gram matrices, and are the elements of the incidence matrix, y_i and y_j are the elements of the vector of phenotypic value for all the individuals, and are kernel functions, is the centering matrix, is the n-dimensional identity matrix, and is the n-dimensional vector with all ones. Here, a non-negativity constraint is used so that meaningful markers are selected. As output, gram matrix L is used to select markers in HSIC LASSO, and it is possible to naturally incorporate structured outputs via kernels. Moreover, we can perform feature selection even if the training data set consists of input x and its affinity information L link structures between inputs. Differences from the original formulation of LASSO are that, in this case, kernel functions K and L are different and non-negativity constraint is imposed. The first term in this equation means that we are regressing the output kernel matrix by a linear combination of feature-wise input kernel matrices.

2.3. Proposed integrated framework for marker selection

In this study, Linear least-squared regression (Meuwissen et al., 2001), LASSO (Tibshirani, 1996), ridge regression (Hoerl and Kennard, 1970), BLUP (Henderson, 1975), SpAM (Ravikumar et al., 2009), mRMR (Peng et al., 2005), and HSIC LASSO (Yamada et al., 2013) have been evaluated through a simulation study for identification of additive and nonadditive genetic markers based on prediction accuracy (PA), fraction of correctly selected features, and redundancy rate (RED). It was found that SpAM and HSIC LASSO are reasonably superior for identification of additive and nonadditive genetic markers (Guha Majumdar et al., in press). Therefore, it is proposed to estimate GEBV for GS based on additive and nonadditive genetic markers identified through these procedures. To estimate GEBV following an integrated model is proposed:

(3)

where is the predicted phenotype (GEBV) of the integrated model, w is , where and are the error variances of models HSIC LASSO and SpAM, respectively, is the predicted GEBV from SpAM and is the GEBV from HSIC LASSO. Let the error variance of be denoted by . The can be written after optimizing w, as follows:

(4)

The estimation of the above error variance is highly complex in nature. Therefore, it is proposed to estimate this by following the approach of Fan et al. (2012). The detail of the same is described in the following section.

2.4. Estimation of error variance of proposed integrated framework

There are several methods to estimate error variances and . In this study, we have considered the refitted cross-validation (RCV) method of error variance estimation (Fan et al., 2012). The algorithm of the RCV method is described below:

i.	Consider a data set . We denote it as .
ii.	is divided into two groups randomly, each subsample having size of n/2. where and
iii.	Variable selection has been performed on using SpAM or HSIC LASSO and let be a set of selected variables from subsample .
iv.	Using only the selected variables, error variance is estimated through with ordinary least-squares estimation. where
v.	In the next stage, variable selection has been performed on using SpAM or HSIC LASSO and is a set of selected variables from subsample .
vi.	Using only the selected variables, variance is estimated in with ordinary least-squares estimation. where
vii.	The final estimator of variance is (5)

Here, represents both and . This means both and will be obtained by following the above algorithm of RCV. Let error variance estimator of SpAM or HSIC LASSO be represented by and , respectively.

Then, estimator of error variance of integrated model, that is, , is given by substituting and with estimated values of and , that is, and in expression (iv) of .

2.5. Simulation study

A genome with 10 chromosomes with specified length was simulated using R package QTL (quantitative trait locus) Bayesian interval mapping (“qtlbim”; Yandell et al., 2007; Yandell et al., 2012). The package qtlbim follows Cockerham's model for simulating quantitative trait loci with epistasis (Kao and Zeng, 2002). The Cockerham's model for quantitative trait was controlled by two epistatic genes A and B, from a sample of size n of an F2 population, and the trait value of the k^th individual with genotype ij can be described as follows:

(6)

where denotes the genotypic value of the genotype , is the mean, a₁ is the additive effect of locus A, d₁ is the dominance effect of locus A, a₂ is the additive effect of locus B, d₂ is the dominance effect of locus B, is additive × additive effect of loci A and B, is additive × dominance effect of loci A and B, is dominance × additive effect of loci A and B, and is dominance × dominance effect of loci A and B. The coded variables are defined as follows:

and is residual of the k^th individual with genotype ij.

The genome contains 1000 markers, which were equally spaced over the chromosomes and each chromosome has 100 markers. The phenotypic values are normally distributed. For the additive model, there is one QTL in each chromosome with either a positive or negative additive effect and no epistatic interaction. For the epistatic model, we assumed two QTLs on each of the five chromosomes and the remaining five chromosomes have no QTL resulting in 5 two-way epistatic interactions. In this way, 24 different data sets are generated with genotypic and phenotypic information for the F2 population. Each data set contains 200 individuals and 1000 biallelic markers. Each data set has a different combination of heritability (narrow sense) and epistatic effect. Details of simulated data sets are given in Table 1. For each data set with a specific genetic architecture, 50 replicates were generated and from each replicate, 80% individuals were selected randomly for training and the remaining 20% individuals are kept for testing, that is, to predict phenotype (GEBV). Again for each of these 50 replicates, this random selection of training/testing data sets was repeated 50 times; this results in 2500 training/testing data sets in total for a specific genetic architecture.

From these simulated data sets, GEBV was estimated using above models, that is, SpAM, HSIC LASSO, and proposed integrated model framework using R (R Core Team, 2017). For SpAM, we have used SAM package (Zhao et al., 2014) and applied samQL function with default parameter values. Ten highly significant markers have been selected and used for the prediction of GEBV of the testing data set with the help of predict function using SpAM fitted model on the training data set. For HSIC LASSO, a function in R has been developed to kernelize the input variable matrix and output vector. Then, the penalized function of penalized package (Goeman et al., 2010) has been used to fit this kernelized LASSO model and 10 highly significant markers were selected to predict phenotype (GEBV) from testing data set using predict function. Furthermore, to predict the GEBV from the proposed integrated model framework, R code has been written to combine the selected markers from these two models.

Then the performance of SpAM, HSIC LASSO, and the proposed integrated model framework has been evaluated on the basis of PA, fraction of correctly selected features (F), and the RED. PA can be described as the correlation between the actual phenotypic () values and the predicted phenotypic () values (Howard et al., 2014).

The RED score (Zhao et al., 2010) is obtained by the following:

where is the correlation coefficient between the k-th and l-th markers. A high RED score implies that selected markers are strongly correlated to each other, which means many redundant markers are selected. So, lower RED is desirable for feature selection.

Fraction of correctly selected features (Yamada et al., 2013) is expressed by the following:

where F is fraction of correctly selected features, N_c is the number of correctly selected features, and N_t is the number of total selected features.

3. Results

The error variance estimation method RCV has been applied to the simulated data sets using R platform. The values of the error variances of SpAM, HSIC LASSO, and integrated model framework are given in Table 2.

Table 2 shows the estimated error variances of GEBV estimates obtained from different models. The following points may be noted from the results of Table 2:

Also, the performances of each model for selection of genetic markers and estimation of GEBV on the basis of PA, fraction of correctly selected features (F), and RED, along with their standard error (SE) of mean, are given in Supplementary Table S1.

The relationship among PA, F, and RED is presented in Figures 1 to 3 for different levels of heritability and epistatic effects. The following points may be noted from these figures.

4. Discussion

Statistical models are extensively used for determination of superior parents for breeding on the basis of their phenotypic performance. With the advent of high-throughput genotyping technology, it becomes possible to estimate GEBV of individuals from the breeding population by combining genotypic and phenotypic information. Statistical models are available for estimation of GEBV for GS. The accuracy of GEBV estimate depends on selection of relevant genetic markers associated with trait of interest. Most of the procedures for identification of these genetic markers are able to capture additive genetic effects; recently, attempts have been made to develop techniques to identify markers that can capture nonadditive genetic effects. It is pertinent to note that the estimation of GEBV of an individual using any of these set of markers, that is, additive or nonadditive may lead to underestimate the genetic potential of an individual. Therefore, it is pertinent to use both categories of markers to estimate GEBV reliably. Considering this, several additive and nonadditive models are evaluated on the basis of their performance for marker identification. These procedures include linear least-squared regression, LASSO, ridge regression, BLUP, mRMR, SpAM, and HSIC LASSO.

In this article, an attempt has been made to integrate selected markers from most efficient additive and nonadditive marker selection procedures available in the literature for GS. From the above results, it can be observed that the performance of SpAM, that is, additive effect model, is superior in terms of SE of estimates in the absence of epistatic effect especially at a low heritability level, which is the most practical situation in the real world. It may be noted that in the presence of nonadditive effects, the performance of HSIC LASSO, that is, model for nonadditive effect, is far superior to SpAM due to obvious reasons. In the proposed integrated framework that captures both additive and nonadditive genetic effects, the SE is quite low in comparison with both of these estimates of GEBV as this is derived by optimally combining the estimates from these models. In terms of PA, the integrated framework performs very well in case of dominance of additive effects in the genetic architecture, which is the normally occurring phenomenon in the real world. The proposed integrated model also performs well in terms of fraction of correctly selected features in the presence of moderate and high epistatic effects, which shows that it is robust against nonlinear dominant effects. The performance of all three models in terms of RED score is comparable. This shows that all three models are efficient in selection of independent features, that is, markers. These results are based on empirical comparison through a simulation study, and therefore, there is a need to validate the same on a real data set.

5. Conclusion

In this article, attempt has been made to identify superior techniques for selection of additive and nonadditive genetic markers among competitive methods. Also, a robust integrated approach is proposed based on selected genetic markers to estimate GEBV, which performs well for both additive and nonadditive genetic architectures. The proposed framework has been developed by optimally combining selected genetic markers obtained from the additive model, that is, SpAM and the nonadditive model, that is, HSIC LASSO, which are found to be better among all other methods considered in this study. The performance of the proposed framework is found to be satisfactory in terms of important evaluation measures, that is, SE, PA, F, and RED in most practical situations of the real world.

Authors' Contributions

A.R. conceptualized, guided, and wrote the article. S.G. and D.M. performed the computational and analysis part. All authors read and approved the final article.

Author Disclosure Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. No competing financial interests exist.

Funding Information

The authors received no specific funding for this work.

Supplementary Material

Supplementary Table S1