Cloud Hosting

The FE-related traits and genomic information were obtained for 1,156 animals from an experimental breeding program at the Beef Cattle Research Center (Institute of Animal Science IZ).

Animals were from an experimental breeding program at the Beef Cattle Research Center at the Institute of Animal Science (IZ) in Sertozinho, So Paulo, Brazil. Since the 1980s, the experimental station has maintained three selection herds: Nellore control (NeC) with animals selected for yearling body weight (YBW) with a selection differential close to zero, within birth year and herd, while Nellore Selection (NeS) and Nellore Traditional (NeT) animals are selected for the YBW with a maximum selection differential, also within birth year and herd25. In the NeT herd, sires from commercial herds or NeS eventually were used in the breeding season, while the NeC and NeS were closed herds (only sires from the same herd were used in the breeding season), with controlled inbreeding rate by planned matings. In addition, the NeT herd has been selected for lower residual feed intake (RFI) since 2013. In the three herds, the animal selection is based on YBW measured at 378days of age in young bulls.

The FE-related traits were evaluated on 1156 animals born between 2004 and 2015 in a feeding efficiency trial, in which they were either housed in individual pens (683 animals) or group pens equipped with the GrowSafe feeding system (473 animals), with animals grouped by sex. From those, 146 animals were from the NeC herd (104 young bulls and 42 heifers), 300 from the NeS herd (214 young bulls and 86 heifers), and 710 from the NeT herd (483 young bulls and 227 heifers). Both feeding trials comprised at least 21 days for adaptation to the feedlot diet and management and at least 56 days for the data collection period. The young bull and heifers showed an average age at the end of the feeding trial was 36627.5 and 38445.4 days, respectively.

A total of 780 animals were genotyped with the Illumina BovineHD BeadChip assay (770k, Illumina Inc., San Diego, CA, USA), while 376 animals were genotyped with the GeneSeek Genomic Profiler (GGP Indicus HD, 77K). The animals genotyped with the GGP chip were imputed to the HD panel using FImpute v.326 with an expected accuracy higher than 0.97. Autosomal SNP markers with a minor allele frequency (MAF) lower than 0.10 and a significant deviation from HardyWeinberg equilibrium (P105) were removed, and markers and samples with call rate lower than 0.95 were also removed. An MAF lower than 10% was used to remove genetic markers with lower significance and noise information in a stratified population. After this quality control procedure, genotypes from 1,024 animals and 305,128 SNP markers remained for GS analyses. Population substructure was evaluated using a principal component analysis (PCA) based on the genomic relationship matrix using the ade4 R package (Supplementary Figure S1)27.

Animals were weighed without fasting at the beginning and end of the feeding trial, as well as every 14 days during the experimental period. The mixed ration (dry corn grain, corn silage, soybean, urea, and mineral salt) was offered ad libitum and formulated with 67% of total digestible nutrients (TDN) and 13% of crude protein (CP), aiming for an average daily gain (ADG) of 1.1kg.

The following feed efficiency-related traits were evaluated: ADG, dry matter intake (DMI), feed efficiency (FE), and RFI. In the individual pens, the orts were weighed daily in the morning before the feed delivery to calculate the daily dietary intake. In the group pens, the GrowSafe feeding system automatically recorded the feed intake. Thus, the DMI (expressed as kg/day) was estimated as the feed intake by each animal with subsequent adjustments for dry matter content. ADG was estimated as the slope of the linear regression of body weight (BW) on feeding trial days, and the FE was expressed as the ratio of ADG and DMI. Finally, RFI was calculated within each contemporary group (CG), as the difference between the observed and expected feed intake considering the average metabolic body weight (MBW) and ADG of each animal (Koch et al., 1963) as follows:

$$DMI=CG+ {beta }_{0}+{beta }_{1}ADG+{beta }_{2}MBW+varepsilon$$

where ({beta }_{0}) is the model intercept, ({beta }_{1}) and ({beta }_{2}) are the linear regression coefficients for (ADG) and ({MBW=BW}^{0.75}), respectively, and (varepsilon) is the residual of the equation representing the RFI estimate.

The contemporary groups (CG) were defined by sex, year of birth, type of feed trial pen (individual or collective) and selection herd. Phenotypic observations with values outside the interval of3.5 standard deviations below and above the mean of each CG for each trait were excluded, and the number of animals per CG ranged from 10 to 70.

The (co)variance components and heritability for FE-related traits were estimated considering a multi-trait GBLUP (MTGBLUP) as follows:

$$mathbf{y}=mathbf{X}{varvec{upbeta}}+mathbf{Z}mathbf{a}+mathbf{e},$$

Where ({varvec{y}}) is the matrix of phenotypic FE-related traits (ADG, FE, DMI, and RFI) of dimension Nx4 (N individuals andfour traits); ({varvec{upbeta}}) is the vector of fixed effects, linear and quadratic effects of cow age, and linear effect of animals age at the beginning of the test; (mathbf{a}) is the vector of additive genetic effects (breeding values) of animal, and (mathbf{e}) is a vector with the residual terms. The (mathbf{X}) and (mathbf{Z}) are the incidence matrices related to fixed (b) and random effects (a), respectively. It was assumed that the random effects of animals and residuals were normally distributed, as (mathbf{a}sim {text{N}}(0,mathbf{G}otimes {mathbf{S}}_{mathbf{a}})) and (mathbf{e}sim {text{N}}(0,mathbf{I}otimes {mathbf{S}}_{mathbf{e}})), where (mathbf{G}) is the additive genomic relationship matrix between genotyped individuals according to VanRaden28, (mathbf{I}) is an identity matrix,is the Kronecker product, and ({mathbf{S}}_{mathbf{a}}=left[begin{array}{ccc}{upsigma }_{{text{a}}1}^{2}& cdots & {upsigma }_{mathrm{a1,4}}\ vdots & ddots & vdots \ {upsigma }_{mathrm{a1,4}}& cdots & {upsigma }_{{text{a}}4}^{2}end{array}right]) and ({mathbf{S}}_{mathbf{e}}=left[begin{array}{ccc}{upsigma }_{{text{e}}1}^{2}& cdots & {upsigma }_{mathrm{e1,4}}\ vdots & ddots & vdots \ {upsigma }_{mathrm{e1,4}}& cdots & {upsigma }_{{text{e}}4}^{2}end{array}right]) are the additive genetic and residual (co)variance matrices, respectively. The G matrix was obtained according to VanRaden28: (mathbf{G}=frac{mathbf{M}{mathbf{M}}^{mathbf{^{prime}}}}{2sum_{{text{j}}=1}^{{text{m}}}{{text{p}}}_{{text{j}}}left(1-{{text{p}}}_{{text{j}}}right)}) where (mathbf{M}) is the SNP marker matrix with codes 0, 1, and 2 for genotypes AA, AB, and BB adjusted for allele frequency expressed as (2{{text{p}}}_{{text{j}}}), and ({{text{p}}}_{{text{j}}}) is the frequency of the second allele jth SNP marker.

The analyses were performed using the restricted maximum likelihood (REML) method through airemlf90 software29. The predictf90 software29 was used to obtain the phenotypes adjusted for the fixed effects and covariates (({{text{y}}}^{*}={text{y}}-{text{X}}widehat{upbeta })). The adjusted phenotypes were used as the response variable in the genomic predictions.

Tthe GEBVs accuracy (({{text{Acc}}}_{{text{GEBV}}})) in the whole population, was calculated based on prediction error variance (PEV) and the genetic variance for each FE-related trait (({upsigma }_{{text{a}}}^{2})) using the following equation30: ({text{Acc}}=1-sqrt{{text{PEV}}/{upsigma }_{{text{a}}}^{2}}) .

A forward validation scheme was applied for computing the prediction accuracies using machine learning and parametric methods, splitting the dataset based on year of birth, with animals born between 2004 and 2013 assigned as the reference population (n=836) and those born in 2014 and 2015 (n=188) as the validation set. For ML approaches, we randomly split the training dataset into fivefold to train the models.

Genomic prediction for FE-related traits considering the STGBLUP can be described as follows:

$${mathbf{y}}^{mathbf{*}}={varvec{upmu}}+mathbf{Z}mathbf{a}+mathbf{e}$$

where ({mathbf{y}}^{mathbf{*}}) is the Nx1 vector of adjusted phenotypic values for FE-related traits, (upmu) is the model intercept, (mathbf{Z}) is the incidence connecting observations; (mathbf{a}) is the vector of predicted values, assumed to follow a normal distribution given by ({text{N}}(0,{mathbf{G}}sigma_{a}^{2})) and (mathbf{e}) is the Nx1 vector of residual values considered normally distributed as ({text{N}}(0,mathbf{I}{upsigma }_{{text{e}}}^{2})), in which I is an identity matrix, ({upsigma }_{{text{e}}}^{2}) is the residual variance. The STGBLUP model was performed using blupf90+software29.

Genomic prediction for FE-related traits considering MTGBLUP can be described as follows:

$${mathbf{y}}^{mathbf{*}}={varvec{upmu}}+mathbf{Z}mathbf{a}+mathbf{e}$$

where ({mathbf{y}}^{mathbf{*}}) is the matrix of adjusted phenotypes of dimension Nx4, (upmu) is the trait-specific intercept vector, (mathbf{Z}) is the incidence matrix for the random effect; (mathbf{a}) is an Nx4 matrix of predicted values, assumed to follow a normal distribution given by ({text{MVN}}(0,{mathbf{G}} otimes {mathbf{S}}_{{mathbf{a}}})) where ({mathbf{S}}_{mathbf{a}}) represents genetic (co)variance matrix for the FE-related traits (44). The residual effects (e) were considered normally distributed as ({text{MVN}}(0,mathbf{I}otimes {mathbf{S}}_{mathbf{e}})) in which I is an identity matrix, and ({mathbf{S}}_{mathbf{e}}) is the residual (co)variance matrix for FE-related traits (44). The MTGBLUP was implemented in the BGLR R package14 considering a Bayesian GBLUP with a multivariate Gaussian model with an unstructured (co)variance matrix between traits (({mathbf{S}}_{mathbf{a}})) using Gibbs sampling with 200,000 iterations, including 20,000 samples as burn-in and thinning interval of 5 cycles. Convergence was checked by visual inspection of trace plots and distribution plots of the residual variance.

Five Bayesian regression models with different priors were used for GS analyses: Bayesian ridge regression (BRR), Bayesian Lasso (BL), BayesA, BayesB, and BayesC. The Bayesian algorithms for GS were implemented using the R package BGLR version 1.0914. The BGLR default priors were used for all models, with 5 degrees of freedom (dfu), a scale parameter (S), and . The Bayesian analyses were performed considering Gibbs sampling chains of 200,000 iterations, with the first 20,000 iterations excluded as burn-in and a sampling interval of 5 cycles. Convergence was checked by visual inspection of trace plots and distribution plots of the residual variance. For Bayesian regression methods, the general model can be described as follows:

$${mathbf{y}}^{mathbf{*}}=upmu +sum_{{text{w}}=1}^{{text{p}}}{{text{x}}}_{{text{iw}}}{{text{u}}}_{{text{w}}}+{{text{e}}}_{{text{i}}}$$

where (upmu) is the model intercept; ({{text{x}}}_{{text{iw}}}) is the genotype of the ith animal at locus w (coded as 0, 1, and 2); ({{text{u}}}_{{text{w}}}) is the SNP marker effect (additive) of the w-th SNP (p=305,128); and ({{text{e}}}_{{text{i}}}) is the residual effect associated with the observation of ith animal, assumed to be normally distributed as (mathbf{e}sim {text{N}}(0,{mathbf{I}upsigma }_{{text{e}}}^{2})).

The BRR method14 assumes a Gaussian prior distribution for the SNP markers (({{text{u}}}_{{text{w}}})), with a common variance ({(upsigma }_{{text{u}}}^{2})) across markers so that ({text{p}}left({{text{u}}}_{1},dots ,{{text{u}}}_{{text{w}}}|{upsigma }_{{text{u}}}^{2}right)=prod_{{text{w}}=1}^{{text{p}}}{text{N}}({{text{u}}}_{{text{w}}}{|0,upsigma }_{{text{u}}}^{2})). The variance of SNP marker effects is assigned a scaled-inverse Chi-squared distribution [({text{p}})(({upsigma }_{{text{u}}}^{2})={upchi }^{-2}({upsigma }_{{text{u}}}^{2}|{{text{df}}}_{{text{u}}},{{text{S}}}_{{text{u}}}))], and the residual variance is also assigned a scaled-inverse Chi-squared distribution with degrees of freedom (dfe)and scale parameters (Se).

Bayesian Lasso (BL) regression31 used an idea from Tibshirani32 to connect the LASSO (least absolute shrinkage and selection operator) method with the Bayesian analysis. In the BL, the source of variation is split intoresidual term(({upsigma }_{{text{e}}}^{2}))and variation due to SNP markers (({upsigma }_{{{text{u}}}_{{text{w}}}}^{2})). The prior distribution for the additive effect of the SNP marker (left[{text{p}}left({{text{u}}}_{{text{w}}}|{uptau }_{{text{j}}}^{2},{upsigma }_{{text{e}}}^{2}right)right]) follows a Gaussian distribution with marker-specific prior variance given by ({text{p}}left({{text{u}}}_{{text{w}}}|{uptau }_{{text{j}}}^{2},{upsigma }_{{text{e}}}^{2}right)=prod_{{text{w}}=1}^{{text{p}}}{text{N}}({{text{u}}}_{{text{w}}}left|0,{uptau }_{{text{j}}}^{2}{upsigma }_{{text{e}}}^{2}right)). This prior distribution leads to marker-specific shrinkage of their effect, whose their extent depends on the variance parameters (left({uptau }_{{text{j}}}^{2}right)). The variance parameters (left({uptau }_{{text{j}}}^{2}right)) is assigned as exponential independent and identically distributed prior,({text{p}}left( {{uptau }_{{text{j}}}^{2} left| {uplambda } right.} right) = mathop prod limits_{{{text{j}} = 1}}^{{text{p}}} {text{Exp}}left( {{uptau }_{{text{j}}}^{2} left| {{uplambda }^{2} } right.} right)) and the square lambda regularization parameter (({uplambda }^{2})) follows a Gamma distribution (({text{p}}left({uplambda }^{2}right)={text{Gamma}}({text{r}},uptheta ))), where r and (uptheta) are the rate and shape parameters, respectively31. Thus, the marginal prior for SNP markers is given by a double exponential (DE) distribution as follows: ({text{p}}left( {{text{u}}_{{text{w}}} left| {uplambda } right.} right) = int {{text{N}}left( {{text{u}}_{{text{w}}} left| {0,{uptau }_{{text{j}}}^{2} ,{upsigma }_{{text{e}}}^{2} } right.} right){text{Exp}}left( {{uptau }_{{text{j}}}^{2} left| {{uplambda }^{2} } right.} right)}), where the DE distribution places a higher density at zero and thicker tails, inducing stronger shrinkage of estimates for markers with relatively small effect and less shrinkage for markers with substantial effect. The residual variance (({upsigma }_{{text{e}}}^{2})) is specified as a scaled inverse chi-squared prior density, with degrees of freedom dfe and scale parameter Se.

BayesA method14,33 considers Gaussian distribution with null mean as prior for SNP marker effects (({{text{u}}}_{{text{w}}})), and a SNP marker-specific variance (({upsigma }_{{text{w}}}^{2})). The variance associated with each marker effect assumes a scaled inverse chi-square prior distribution, ({text{p}}left({upsigma }_{{text{w}}}^{2}right)={upchi }^{-2}left({upsigma }_{{text{w}}}^{2}|{{text{df}}}_{{text{u}}},{{text{S}}}_{{text{u}}}^{2}right)), with degrees of freedom (({{text{df}}}_{{text{u}}})) and scale parameter (({{text{S}}}_{{text{u}}}^{2})) treated as known14. Thus, BayesA places a t-distribution for the markers effects, i.e., ({text{p}}left({{text{u}}}_{{text{w}}}|{{text{df}}}_{{text{u}}},{{text{S}}}^{2}right)={text{t}}left(0,{{text{df}}}_{{text{u}}},{{text{S}}}_{{text{u}}}^{2}right)), providing a thicker-tail distribution compared to the Gaussian, allowing a higher probability of moderate to large SNP effects.

BayesB assumes that a known proportion of SNP markers have a null effect (i.e., a point of mass at zero), and a subset of markers with a non-null effect that follow univariate t-distributions3,12, as follows:

$${text{p}}left({{text{u}}}_{{text{w}}}|{text{df}},uppi ,{{text{df}}}_{{text{u}}},{S}_{B}^{2}right)=left{begin{array}{cc}0& mathrm{with probability pi }\ {text{t}}left({{text{u}}}_{{text{w}}}|{{text{df}}}_{{text{u}}},{S}_{B}^{2}right)& mathrm{with probability }left(1-uppi right)end{array}right.$$

where (uppi) is the proportion of SNP markers with null effect, and (1-uppi) is the probability of SNP markers with non-null effect contributing to the variability of the FE-related trait3. Thus, the prior distribution assigned to SNP with non-null effects is a scaled inverse chi-square distribution.

BayesC method34 assumes a spikeslab prior for marker effects, which refers to a mixture distribution comprising a fixed amount with probability (uppi) of SNP markers have a null effect, whereas a probability of 1 of markers have effects sampled from a normal distribution. The prior distribution is as follows:

$${text{p}}left({{text{u}}}_{{text{w}}},{upsigma }_{{text{w}}}^{2},uppi right)=left{prod_{{text{j}}=1}^{{text{w}}}left[uppi left({{text{u}}}_{{text{w}}}=0right)+left(1-uppi right){text{N}}(0,{upsigma }_{{text{w}}}^{2})right]*{upchi }^{-2}left({upsigma }_{{text{w}}}^{2}|{{{text{df}}}_{{text{u}}},mathrm{ S}}_{{text{B}}}^{2}right)*upbeta (uppi |{{text{p}}}_{0},{uppi }_{0}right},$$

Where ({upsigma }_{{text{w}}}^{2}) is the common variance for marker effect, ({{text{df}}}_{{text{u}}}) and ({{text{S}}}_{{text{B}}}^{2}) are the degrees of freedom and scale parameter, respectively, ({{text{p}}}_{0}) and ({uppi }_{0})[0,1] are the prior shape parameters of the beta distribution.

Two machine learning (ML) algorithms were applied for genomic prediction: Multi-layer Neural Network (MLNN) and support vector regression (SVR). The ML approaches were used to alleviate the standard assumption adopted in the linear methods, which restrict to additive genetic effects of markers without considering more complex gene action modes. Thus, ML methods are expected to improve predictive accuracy for different target traits. To identify the best combination of hyperparameters (i.e., parameters that must be tuned to control the learning process to obtain a model with optimal performance) in the supervised ML algorithms (MLNN and SVR), we performed a random grid search by splitting the reference population from the forward scheme into five-folds35.

In MLNN, handling a large genomic dataset, such as 305,128 SNPs, is difficult due to the large number of parameters that need to be estimated, leading to a significant increase in computational demand36. Therefore, an SNP pre-selection strategy based on GWAS results in the training population using an MTGBLUP method (Fig.1A) was used to reduce the number of markers to be considered as input on the MLNN. In addition, this strategy can remove noise information in the genomic data set. In this study, the traits displayed major regions explaining a large percentage of genetic variance, which makes using pre-selected markers useful37.

(A) Manhattan plot for percentage of genetic variance explained by SNP-marker estimated through multi-trait GWAS in training population to be used as pre-selection strategies for multi-layer neural network. (B) General representation of neural networks with two hidden layers used to model nonlinear dependencies between trait and SNP marker information. The input layer ((X={x}_{i,p})) considered in the neural network refers to the SNP marker information (coded as 0, 1, and 2) of the ith animal. The selected node represents the initial weight ((W={w}_{p})), assigned as random values between -0.5 and 0.5, connecting each input node to the first hidden layer and in the second layer the ({w}_{up}) refers to the output weight from the first hidden layer, b represents the bias which helps to control the values in the activation function. The output ((widehat{y})) layer represents a weighted sum of the input features mapped in the second layer.

The MLNN model can be described as a two-step regression38. The MLNN approach consists of three different layer types: input layer, hidden layer, and output layer. The input layer receives the input data, i.e., SNP markers. The hidden layer contains mapping processing units, commonly called neurons, where each neuron in the hidden layer computes a non-linear function (activation) of the weighted sum of nodes on the previous layer. Finally, the output layer provides the outcomes of the MLNN. Our proposed MLNN architecture comprises two fully connected hidden layers schematically represented in Fig.1B. The input layer in MLNN considered SNP markers that explained more than 0.125% of the genetic variance for FE-related traits (Fig.1A;~15k for ADG and DMI, and~16k for FE and RFI). The input covariate (X={{x}_{p}}) contains pre-selected SNP markers (p) with a dimension Nxp (N individuals and p markers). The pre-selected SNP markers are combined with each k neuron (with k=1, , Nr) through the weight vector ((W)) in the hidden layer and then summed with a neuron-specific bias (({b}_{k})) for computing the linear score for the neuron k as:({Z}_{i}^{[1]}=f({{b}_{k}}^{[1]}+X{W}^{[1]})) (Fig.1B). Subsequently, this linear score transformed using an activation function (fleft(.right)) to map k neuron-specific scores and produce the first hidden layer output ((fleft({z}_{1,i}right))). In the second-hidden layer, each neuron k receives a net input coming from hidden layer 1 as: ({Z}_{i}^{[2]}={{b}_{k}}^{left[2right]}+{Z}_{i}^{[1]}{W}^{[2]}), where ({W}^{[2]}) represents the weight matrix of dimension k x k (knumber of neurons) connecting the ({Z}_{i}^{[1]}) into the second hidden layer, and ({{b}_{k}}^{left[2right]}) is a bias term in hidden layer 2. Then, the activation function is applied to map the kth hidden neuron unit in the second hidden layer and generate the output layer as ({V}_{2,i}=fleft({z}_{2,i}right)). In the MLNN, a hyperbolic tangent activation function (({text{tanh}}left({text{x}}right)={{text{e}}}^{{text{x}}}-{{text{e}}}^{-{text{x}}}/{{text{e}}}^{{text{x}}}+{{text{e}}}^{-{text{x}}})) was adopted in the first and second layers, providing greater flexibility in the MLNN39.

The prediction of the adjusted FE-related trait was obtained as follows38:

$${mathbf{y}}^{mathbf{*}}=mathbf{f}left(mathbf{b}+{mathbf{V}}_{2,mathbf{i}}{mathbf{W}}_{0}right)+mathbf{e}$$

where ({mathbf{y}}^{mathbf{*}}) represents the target adjusted feed efficiency-related trait for the ith animal; (k) the number of neurons considered in the model and assumed the same in the first and second layer; ({mathbf{W}}_{0}) represents the weight from the k neuron in layer 2, (mathbf{b}) is related to the bias parameter. The optimal weights used in MLNN were obtained by minimizing the mean square error of prediction in the training subset40.

The MLNN model was implemented using the R package h2o (https://github.com/h2oai/h2o-3), with the random grid search using the h2o.grid function (https://cran.r-project.org/web/packages/h2o) to determine the number of neurons to maximize the prediction accuracy. We used the training population split into fivefold to assess the best neural network architecture and then apply it in the disjoint validation set41,42. We considered a total of 1000 epochs36, numbers of neurons ranging from 50 to 2500 with intervals of 100, and applied a dropout ratio of 0.2 and regularization L1 and L2 parameters as 0.0015 and 0.0005, respectively. In this framework, the MLNN was performed using two hidden layers of neural networks with the number of neurons (k) of 750 for ADG, 1035 for DMI, 710 for FE, and 935 for RFI obtained during the training process.

Support vector regression (SVR) is a kernel-based supervised learning technique used for regression analysis43. In the context of GS, the SVR uses linear models to implement nonlinear regression by mapping the predictor variables (i.e., SNP marker) in the feature space using different kernel functions (linear, polynomial, or radial basis function) to predict the target information, e.g., adjusted phenotype the GS44. SVR can map linear or nonlinear relationships between phenotypes and SNP markers depending on the kernel function. The best kernel function mapping genotype to phenotype (linear, polynomial, and radial basis) was determined using the training subset split into fivefold. The radial basis function (RBF) was chosen as it outperformed the linear and polynomial (degree equal 2) kernels in the training process, increasing 8.25% in predictive ability and showing the lowest MSE.

The general model for SVR using a RBF function can be described as38,45: ({mathbf{y}}_{mathbf{i}}^{mathbf{*}}=mathbf{b}+mathbf{h}{left(mathbf{m}right)}^{mathbf{T}}mathbf{w}+mathbf{e}), where (mathbf{h}{left(mathbf{m}right)}^{mathbf{T}}) represents the kernel radial basis function used to transform the original predictor variables, i.e. SNP marker information (({text{m}})), (b) denotes the model bias, and (w) represents the unknown regression weight vector. In the SVR, the learn function (mathbf{h}{left(mathbf{m}right)}^{mathbf{T}}) was given by minimizing the loss function. The SVR was fitted using an epsilon-support vector regression that ignores residual absolute value ((left|{y}_{i}^{*}-{widehat{y}}_{i}^{*}right|)) smaller than some constant () and penalize larger residuals46.

The kernel RBF function considered in the SVR follows the form: (mathbf{h}{left(mathbf{m}right)}^{mathbf{T}}=mathbf{exp}left(-{varvec{upgamma}}{Vert {mathbf{m}}_{mathbf{i}}-{mathbf{m}}_{mathbf{j}}Vert }^{2}right)), where the ({varvec{upgamma}}) is a gamma parameter to quantity the shapes of the kernel functions, (m)and({m}_{i}) are the vectors of predictor variables for labels i and j. The main parameters in SVR are the cost parameter (({text{C}})), gamma parameter (({varvec{upgamma}})), and epsilon ((upepsilon)). The parameters ({text{C}}) and (upepsilon) were defined using the training data set information as proposed by Cherkasky and Ma47: ({text{C}}={text{max}}left(left|overline{{{text{y}} }^{*}}+3{upsigma }_{{{text{y}}}^{*}}right|,left|overline{{{text{y}} }^{*}}-3{upsigma }_{{{text{y}}}^{*}}right|right)) and (upepsilon =3{upsigma }_{{{text{y}}}^{*}}left(sqrt{{text{ln}}left({text{n}}right)/{text{n}}}right)), in which the (overline{{{text{y}} }^{*}}) and ({upsigma }_{{{text{y}}}^{*}}) are the mean and the standard deviation of the adjusted FE-related traits on the training population, and n represents the number of animals in the training set. The gamma () was determined through a random search of values varying from 0 to 5, using the training folder split into fivefold. The better-trained SVR model considered the parameter of 2.097 for ADG, 0.3847 for DMI, 0.225 for FE, and 1.075 for RFI. The SVR was implemented using the e1071 R package48.

Prediction accuracy (acc) of the different statistical approaches was assessed by Pearsons correlation between adjusted phenotypes (({{text{y}}}^{*})) and their predicted values (({widehat{{text{y}}}}_{{text{i}}}^{*})) on the validation set, and root mean squared error (RMSE). The prediction bias was assessed using the slope of the linear regression of ({widehat{y}}_{i}^{*}) on ({{text{y}}}^{*}), for each model. The Hotelling-Williams test49 was used to assess the significance level of the difference in the predictive ability of Bayesian methods (BayesA, BayesB, BayesC, BL, and BRR), MTGBLUP, and machine learning (MLNN and SVR) against STGBLUP. The similarity between the predictive performance of the different models was assessed using Wards hierarchical clustering method with an Euclidian distance analysis. The relative difference (RD) in the predictive ability was measured as ({text{RD}}=frac{({{text{r}}}_{{text{m}}}-{{text{r}}}_{{text{STGBLUP}}})}{{{text{r}}}_{{text{STGBLUP}}}}times 100), where ({{text{r}}}_{{text{m}}}) represents the acc of each alternative approach (SVR, MLNN, and MTGBLUP, or Bayesian regression modelsBayesA, BayesB, BayesC, BL, and BRR), and ({{text{r}}}_{{text{STGBLUP}}}) is the predictive ability obtained using the STGBLUP method.

The animal procedures and data sampling presented in this study were approved and performed following the Animal Care and Ethical Committee recommendations of the So Paulo State University (UNESP), School of Agricultural and Veterinary Science (protocol number 18.340/16).

Read more here:
Benchmarking machine learning and parametric methods for genomic prediction of feed efficiency-related traits in ... - Nature.com