Introduction

Advancements in sequencing and genotyping techniques have enabled the direct use of DNA information for the selection of genetically superior individuals. Within this context, Genome-Wide Selection (GWS), as initially proposed by Meuwissen, Hayes, and Goddard (2001), has demonstrated exceptional efficiency in predicting genomic breeding values (GEBVs). This methodology facilitates early individual selection, ultimately enhancing genetic gain per unit of time (Crossa et al., 2017; Voss-Fels, Cooper, & Hayes, 2019). GWS methods can analyze and estimate the effects of thousands of Single Nucleotide Polymorphism (SNP) markers on phenotypes. However, only markers in linkage disequilibrium (LD) with quantitative trait loci (QTLs) are relevant for predicting GEBVs and explaining genetic variations in traits of interest.

The primary statistical methods employed in GWS typically assume distributions for marker effects in the models or rely on implicit regressions or dimensionality reduction techniques (Resende, Silva, & Azevedo, 2014). Nonetheless, there is a growing interest in nonparametric methods, including Artificial Neural Networks (ANNs), for handling GWS data. Specifically, the Multilayer Perceptron Networks (MLP) class of ANNs has shown promising results in genetic breeding studies conducted by researchers (Crossa et al., 2017; Cruz & Nascimento, 2018; Rosado et al., 2020; Sousa et al., 2022). ANNs operate akin to biological neurons, acquiring knowledge through experiences and capable of capturing all available information, including complex feature architectures, to generate decision-making criteria (Long et al., 2010; Howard, Carriquiry, & Beavis, 2014; Cruz & Nascimento, 2018; Montesinos López, Montesinos López, & Crossa, 2022; Xu et al., 2022).

Despite the potential of ANN procedures in genomic analysis, numerous markers can limit their applicability. High-dimensional scenarios often result in substantial computational costs and challenges in the learning process since many resources (markers) represent unnecessary aspects of the ANN search space (Long, Gianola, Rosa, & Weigel, 2011; Ehret, Hochstuhl, Gianola, & Thaller, 2015). In other words, molecular markers not in LD with the QTLs related to the trait of interest can introduce unnecessary complexity when applying ANNs. However, selecting a subset of markers potentially associated with the traits of interest can mitigate the high-dimensionality problem, leading to improved learning and enhanced predictive power of ANN models (Long et al., 2010; Crossa et al., 2017; Sant’Anna et al., 2020a; Silva et al., 2022; Aono et al., 2022).

Among the proposed methods for dimensionality reduction, those based on principal components and independent components merit mention (Azevedo, Resende, Silva, Lopes, & Guimarães, 2013; James, Witten, Hatie, & Tibshirani, 2013; Long et al., 2010; Resende et al., 2014; Costa et al., 2020; Paixão et al., 2022; Fialho et al., 2023).

Additionally, methods for direct selection of variables, such as Bayes B, Bayes Cπ, and Multivariate Adaptive Regression Splines (MARS) (Friedman, 1991; Huang et al., 2020; Costa et al., 2022), and indirect selection through variable importance, such as Boosting and Random Forest (Ho, Schierding, Wake, Saffery, & O’Sullivan, 2019; Sousa et al., 2022), have been proposed. Marker selection has already demonstrated its effectiveness in genomic prediction, as demonstrated by Sant’Anna et al. (2020a), who employed stepwise regression, and Sousa et al. (2022), who used indirect selection via Bagging.

Based on the considerations outlined above, this study aims to assess the accuracy, mean squared error (MSE), and computational effort of a two-step procedure for genomic prediction. Initially, this approach involves either the direct selection of a group of SNPs using Multivariate Adaptive Regression Splines (MARS) or the indirect selection of different sets of the most relevant SNPs identified by Boosting. Subsequently, these selected markers are used for genomic prediction using ANNs. Predictions are made using simulated data, considering two distinct scenarios (SC1 and SC2). In SC1, ten loci control the trait with a heritability of 0.4, while in SC2, 100 loci control the trait with a heritability of 0.2.

Material and methods

Data simulation

The data set was simulated using the GENES software, as described by Cruz (2013). The genome comprised ten linkage groups, each spanning 20 cM and containing 200 SPNs. Thus, approximately 2000 SNPs were distributed equidistantly across the genome to ensure comprehensive coverage. For genetic linkage analysis, an F1 generation was simulated, consisting of one homozygous dominant parent and one homozygous recessive parent. The genotypes of the F1 population were used to generate an F2 mapping population comprising 1,000 individuals. Two traits were simulated, both with a zero degree of dominance (d/a = 0), where "a" and "d" represent the genotypic values of homozygotes and heterozygotes, respectively, and a mean value of 100. However, different narrow-sense heritabilities (h ² _a ) and genetic architectures were combined to create scenarios SC1 and SC2 (Table 1). In SC1, the trait was simulated with a moderate heritability (h ² _a = 0.40) and controlled by ten loci, with each chromosome containing a single QTL. In SC2, the trait had a lower heritability (h ² _a = 0.20), with each chromosome harboring 10 QTLs. Thus, two genetic architectures, explaining equal parts of the genetic variance, were created, with the QTLs distributed within the regions covered by SPNs. Additionally, the proportion of genetic variation associated with the QTLs explained by the markers (r² _mg ) was calculated according to Goddard, Hayes, and Meuwissen (2011), as follows:

Where n is the number of SNPs and n _QTL is the number of QTL.

Marker selection using multivariate adaptive regression splines

Multivariate Adaptive Regression Splines (MARS) is a nonparametric regression method used for modeling the relationship between predictive and dependent variables via basis functions, primarily in high-dimensional problems (Friedman, 1991; Huang et al., 2020).

The model proposed by Friedman (1991) is:

Where c ₀ is the intercept, B _i (X) is a basis function, c _i is the coefficient of B _i and M is the number of basis functions automatically defined by the MARS algorithm (Abdulelah Al-Sudani et al., 2019), and εis the random error. The coefficients c _i (i = 0, 1, ..., M) are estimated via minimization of the residual sum of squares (RSS) (Hastie, Tibshirani, & Friedman, 2009).

MARS involves two phases: forward and backward. In the forward phase, the model incorporates pairs of basic functions that minimize the Residual Sum of Squares (RSS) (Hastie et al., 2009; Park & Kim, 2018). In the backward phase, certain basis functions are excluded to create a more parsimonious model that avoids overfitting. The MARS algorithm can select variables during model building, and this approach was employed to select SPNs, which were then used as inputs for ANN training.

Marker selection using boosting

As described by James et al. (2013), the Boosting method is a statistical learning approach that employs regression trees to adjust the residuals of an initial model. The residuals are updated in each tree, grown sequentially from the previous tree, and the response variable involves a combination of many trees, such that $\widehat{f}\left(x\right)={\sum }_{b=1}^B\lambda {\widehat{f}}^b\left(x\right)$. The function $\widehat{f}(.)$ refers to the final tree resulting from the sequential combination of the previously adjusted trees. The parameter λ is the shrinkage that controls the learning rate of the method. This method also needs to be adjusted with several splits in each tree. This parameter controls the complexity of the boosting and is known as depth.

In this method, the importance of variables, as described by Friedman (2001), can be used to indirectly select sets of the most relevant SPNs. The sets of markers comprised 5, 20, and 50% of the markers, equating to 100, 400, and 1,000 markers used, respectively.

Artificial neural network (ANN)

In this study, the Multilayer Perceptron Network (MLP), initially proposed by Rosenblatt (1958), was employed as the ANN class. The error backpropagation learning algorithm was chosen for adjusting weights during ANN training. The network topology varied, featuring a single hidden layer with one to 30 neurons. The logistic function served as the activation function, and the identity function was used as the output function. The training phase encompassed 1,000 iterations.

The output of the hidden layer is computed by taking a linear combination of the I input variables and then applying the logistic activation function, φ(.). Specifically, the output of the k ^th neuron in the layer is determined as follows:

where K is the number of neurons in the hidden layer, w _ik is the synaptic weight between the ith marker and the kth neuron, and w _0k is the synaptic weight of the bias b ₁ for the kth neuron.

The ANN output is given as follows:

Where w _i ’’ is the synaptic weight between the ith neuron of the hidden layer and the network output layer, and w ₀ ’’ is the synaptic weight of the bias b ₂ .

Cross-validation

The performance of ANN models was evaluated through a k-fold cross-validation process, dividing the population (1000 observations) into five groups, each comprising 200 observations. In each iteration, one group served as the validation population, while the remaining four constituted the training population. This process was repeated five times, ensuring all individuals participated in the validation population exactly once.

Two-step genomic prediction and comparison of marker selection approaches

Initially, the training sets (800 observations) underwent direct marker selection via MARS and indirect marker selection via Boosting in each fold. Subsequently, the selected markers were used as input variables for ANN. Accuracy and mean square error (MSE) were measured using the validation set (200 observations), except for time, which was recorded during the training phase. Accuracy was determined using the Pearson correlation between observed GEBVs and predicted GEBVs, while MSE was calculated based on the mean squared difference between observed and predicted GEBVs (Figure 1).

The analyses were conducted using the R software (R Core Team, 2022).

For MARS adjustment, the "earth" function from the "earth" package was employed. For Boosting, the "gbm" function from the "gbm" package was used. Finally, the ANN analysis was performed using the "mlp" function from the "RSNNS."

Comparison of GWS methodologies

Following the k-fold process and obtaining performance metrics, the network topology with the highest accuracy value (i.e., a network model with the number of neurons ranging from 1 to 30 in the hidden layer) was selected to represent the network with the best predictive capacity within each k-fold iteration. Accuracy, MSE, and time values were summarized by calculating their mean and standard error (SE) to compare the performance of different two-step genomic prediction methods to ANNs used without prior marker selection.

Results

Figures 2, 3, and 4 illustrate the accuracy, MSE, and time required for accessing the results for both scenarios, SC1 and SC2. During the MARS fitting, a fixed number of ten markers was consistently selected, irrespective of the scenario or the training set used in the k-fold process. Consequently, only ten markers out of the total 2,000 were used for the subsequent training and validation of the ANN models.

In SC1, a reduction in accuracy was observed as the number of selected markers increased. In this scenario, the highest accuracy (0.96 ± 0.00) was achieved with markers directly selected via MARS (10 markers), representing a substantial 30% increase compared to the accuracy (0.73 ± 0.01) obtained with no marker selection (Figure 2).

Similarly, in SC2, a pattern analogous to SC1 was observed. However, the highest accuracy (0.98 ± 0.00) was attained when the top 5% most important markers (100 markers) were selected by Boosting, leading to a notable 15% increase in accuracy relative to the accuracy (0.85 ± 0.01) obtained with no marker selection (Figure 2).

In both scenarios, marker selection resulted in a sparse set of markers covering a significant portion of the genome. Supplementary Figures S1 and S2 present the positions of markers selected by the best methodology for SC1 and SC2, along with the corresponding QTLs. The supplementary material is available for download at https://www.licae.ufv.br/pesquisa/.

As expected, predictions obtained with no marker selection exhibited the highest MSE values in both scenarios (Figure 3). However, in SC1, the ANN predictions achieved the lowest MSE when MARS selected the markers (5.33 ± 0.34), representing a substantial 90% reduction in MSE compared to no marker selection. In SC2, the top 5% most important markers identified by Boosting yielded the lowest MSE (11.65 ± 1.29), resulting in an approximately 83% decrease in MSE compared to no marker selection (Figure 3).

In SC1, the computational time for training the ANN varied from 2.38s (with MARS-selected markers) to 27.54s (with no marker selection). In SC2, the range was 3.55s (with the top 5% most important markers selected via Boosting) to 38.59s (with no marker selection). As expected, computational time was longer without marker selection and was generally slightly higher in SC2 compared to SC1 (Figure 4).

Discussion

This study proposed a two-step genomic prediction procedure that combines either MARS or Boosting with an ANN to enhance computational efficiency and predictive capacity. Initially, MARS and Boosting were employed for marker selection in two scenarios (SC1 and SC2). Scenario SC1 involved a trait with h ² _a = 0.4 and ten trait-controlling loci, while scenario SC2 involved a trait with h ² _a = 0.2 and 100 trait-controlling loci. Subsequently, the selected markers were used as input for the ANN to establish the genomic prediction model.

While ANNs are favored in genetic breeding for their flexibility and lack of a priori assumptions about the genotype-phenotype relationship (Glória et al., 2016), they are disadvantaged by their high computational cost compared to traditional genomic prediction methods. For instance, ANNs can be up to nine times slower than GBLASSO (Sousa et al., 2020).

A potential solution to overcome this limitation is to employ the proposed two-step genomic prediction approach, as ANNs alone do not perform marker selection. Marker selection is crucial to reduce the impact of non-informative markers and improve the accuracy of genomic prediction. Additionally, marker selection supports the assumption that not all markers in the genome are in linkage disequilibrium (LD) with the QTLs.

MARS was selected for marker selection in this study due to its nonparametric regression capabilities, which have not been extensively utilized in genetic breeding studies. Nevertheless, MARS offers several advantages, including flexibility in model adjustment and the ability to detect important variables in high-dimensional scenarios (Zabihi, Pourghasemi, Motevalli, & Zakeri, 2019; Nayana, Kumar, & Chesneau, 2022). This potential has also made MARS highly promising for predicting and selecting genomic markers, outperforming conventional statistical models such as GBLUP (Costa et al., 2022). Furthermore, Boosting, an ensemble learning method, has demonstrated success in genetic breeding studies (Silveira, Lima, Nascimento, Nascimento, & Silva, 2020; Sousa et al., 2020; Westhues et al., 2021). It is particularly indicated for regression problems due to its high predictive capabilities, surpassing methods like Random Forest and Support Vector Machine. Additionally, it incorporates automatic indirect marker selection (Ogutu, Piepho, & Schulz-Streeck, 2011).

This study demonstrates the efficiency of marker selection for genomic prediction in both evaluated scenarios (SC1 and SC2). Compared to using ANN without marker selection, a significant improvement of over 30% in predictive capacity was observed when using markers selected with MARS in SC1. Similarly, in SC2, approximately a 15% improvement in predictive capacity was achieved when the top 5% most important markers were indirectly selected via Boosting. It is noteworthy that MARS and Boosting (5%) selected ten markers in SC1 and 100 markers in SC2, which align with the QTL numbers in these scenarios. These results suggest that MARS is an intriguing method for selecting markers in oligogenic genetic architectures. Conversely, the Boosting approach performs better in the SC2 scenario, which comprises 100 QTLs.

The markers selected by MARS and Boosting, as displayed in Supplementary Figures S1 and S2, underscore the effectiveness of using only a few markers to predict GEBVs. Markers with high LD among themselves represent redundant information and may be discarded. Therefore, incorporating information into the model requires only one marker from those with high LD. This implies that high-density panels may not necessarily improve predictive capacity (Song & Hu, 2022).

Sousa et al. (2022) found consistent results with this study when assessing the effect of marker selection on genomic prediction in Coffea canephora. Li et al. (2018) employed machine learning methods to select 3000 markers, achieving a predictive capacity of 0.41 using Random Forest for marker selection and 0.46 using Gradient Boosting Machine (GBM), compared to 0.43 using the entire available panel of SPN markers (38,082). These researchers observed that the predictive capacity increased as the number of markers decreased, reaching optimal values with the use of 3.4% to 6.9% of the total number of available markers.

Sant’Anna et al. (2020a) demonstrated that by using simulated data and selecting markers via stepwise regression, the radial basis function network had improved reliability in predicting and reduced the ANN fitting time by up to 20 times. Additionally, the efficiency of MARS for selecting variables to serve as input in ANNs has been demonstrated by Kao and Chiu (2020). These researchers found that a network adjusted with prior selected variables can be up to three times faster than a network without variable selection while producing more accurate results.

Additionally, it was shown that marker selection significantly reduced the computational time required for fitting the ANN. The time to adjust the ANN for prediction was 11.57 times smaller for the two-step prediction using marker selection by MARS in SC1. While in SC2, the time was up to 10.87 times faster for ANN with the selection of the top 5% most important markers using Boosting. Moreover, an increased time required to fit the ANN model was observed as more markers were included as input variables. For instance, with the top 50% most important markers via Boosting, the time spent to fit the ANN model was only 1.03 times faster than the ANN fitting without marker selection in SC1 and 1.79 times faster than the ANN fitting without marker selection in SC2.

This study supports the hypothesis that marker selection efficiently reduces computational costs while increasing predictive accuracy. The explanation lies in the reduction of the marker pool, which, in turn, reduces the search space for the ANN, contributing to the improved predictive power of the model (Sant’Anna, Silva, Nascimento, & Cruz, 2020b).

Conclusion

Boosting and MARS proved to be efficient methodologies for marker selection. The proposed two-step prediction model emerged as an effective strategy for mitigating computational costs while enhancing the prediction accuracy of an ANN.

Acknowledgements

The authors gratefully acknowledge the Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG), Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq)

References

Publication Dates

History