Genotypes

The information used was from 300 animals (236 females, 64 males) born from 2001 to 2016 in eight herds located in Eastern, Central and Western Mexico. Hair samples were collected for genotyping by the company GeneSeek (Lincoln, https://www.neogen.com/, NE, USA), using the GeneSeek® Genomic Profiler Bovine LDv.4 panel, with 30,000 and 50,000 SNP markers, 150 animals with each Chip. Genotyping was performed on two separate occasions, initially 150 individuals with the 30K Chip and later another group of 150 individuals with the 50K Chip since the 30K Chip was not available at the time. The SNPs in common between the 30K and 50K chips (12,835 SNPs) were used. The proportions of missing values were calculated for each marker and for each individual. The average of missing values per individual was 2.09 % with a standard deviation of 7.50 %. The average call rate (not missing proportion for each marker) was 97.90 % with a standard deviation of 4.66 %. Markers with a call rate of less than 95 % were removed. The genotypes were recoded as AA= 0, AB= 1 and BB= 2, from which a matrix with 300 rows (individuals) and 12,835 columns (markers) was obtained, whose cells take values in the set {0,1,2,-}, where “-” denotes a missing value. For the 12,835 common markers of the two chips, the missing values were randomly imputed, generating samples of the Binomial(2,p^) distribution, where p^ is the frequency of the major allele, calculated from the observed marker genotypes. Monomorphic markers or those with minor allele frequency (MAF) less than 0.04 were removed. After quality control, 9,628 of the 12,835 SNPs in common between the two chips were available for further analyses.

Phenotypes

The phenotypic and pedigree information of the Braunvieh cattle population was obtained from the database of the Mexican Association of Breeders of Registered Swiss Cattle. Records of birth (BW), weaning (WW) and yearling (YW) weights were used for analysis. Phenotype editing was similar for BW, WW and YW, records of animals not genetically related to those genotyped or with missing information for herd, dam’s age and management were discarded. Contemporary groups (CG) were defined by removing animals in CG of 2 individuals or with variance equal to zero. For BW, the CGs were defined by combining the effects of the herd (8 herds), year (1998 to 2016) and season of birth; the seasons of birth were defined considering the Julian calendar, from 80 to 171d, spring; from 172 to 264 d, summer; from 265 to 354 d, autumn; from 355 to 366 d and from 1 to 79 d, winter. After editing data, for BW, 330 records were obtained. For WW and YW, the CGs were defined by combining the effects of the herd (6 herds), year (from 1998 to 2016), season of birth (same as BW) and management. In the case of WW, the management groups were defined in three ways: calves fed their mother’s milk; their mother’s milk plus balanced feed; and milk from their mother and nurse plus a balanced diet. For YW, the management groups were defined in three ways: grazing animals; grazing animals with feed concentrate; and housed animals with a balanced diet. The edition of WW and YW data ended with 267 and 232 records for further analyses. Table 1 shows a summary of the number of animals genotyped, and phenotyped for BW, WW and YW. Figure 1 shows the violin plots for BW, WW and YW, the sample mean is represented by the red dot and the sample median by the horizontal line within the box, from the plot, it is clear that the response variables have an asymmetric distribution and the circles with solid filling in it suggest the presence of outliers.

The sample mean is represented by the red dot and the sample median by the horizontal line inside the box

Figure 1 Violin plots of birth (PN=BW), weaning (PD=WW) and yearling (PA=YW) weights in a Braunvieh cattle population 

Models

Quantile regression model with markers (QRM)

The model for quantile regression is:

where yi is the value of the phenotype of the i-th animal; μ is an intercept; xit=(xi1,…,xip) represents the i-th row of the marker matrix, β=β1,…,βpt is the vector of regression coefficients associated with markers and wi are independent random variables such that their quantile θ∈(0,1) is zero. The estimation of the regression coefficients for a fixed interest quantile θ is obtained by solving the following minimization problem:

where ∑j=1pβj is the sum of the absolute values of the regression coefficients; λ is the penalty parameter that controls the intensity of regularization; and ρθ(⋅) is the function defined as¹:

where ti=yi-μ-xitβ. After estimating the parameters of the model, the breeding values estimated by markers (GEBV) are obtained by the following expression:

where β^j(τ) is the effect of the j-th marker, defined by the functional relationship obtained for the quantile of interest.

The QR model can be extended to include other terms, in particular for growth characteristics, the following model is used:

where yi is the value of the phenotype of the analyzed characteristic (BW, WW or YW) of the i-th animal, μ is an intercept; sit=(si1,...,sif) the i-th row of the incidence matrix for fixed effects (sex, dam’s age, management), ς=(ς1,...,ςf)t the regression coefficients for fixed effects, cit=(ci1,...,cit) the i-th row of the incidence matrix for random effects of contemporary group (54, 43 and 37 for BW, WW and YW), ϱ=(ϱ1,...,ϱt)t random effects of contemporary group, the rest of the terms as described above.

GBLUP

The model is given by:

where zit=(zi1,..,zin) is the i-th row of the matrix that connects phenotypes with genotypes, u=(u1,...,un)t is the vector of random effects for animals. Additive, contemporary group and residual genetic variances are assume Varu=Gσu2, Varc=Iσcg2, and Vare=Iσe2, respectively. The matrix of genomic relationships, G, is calculated as described by Lopez-Cruz et al¹⁴ and Pérez-Rodríguez et al¹⁵; briefly, G = WW’/p, where W is the standardized and centered marker matrix (each marker centered by subtracting the mean allele frequency and standardized by dividing by the standard deviation of the sample of the allele frequency), p is the total number of markers, ei normal and independent random variables with normal distribution with mean 0 and variance σe2.

Single-step quantile regression (QRH) model

This method is considered an extension of the quantile model for a relationship matrix constructed using matrices of relationships for genotyped and non-genotyped animals and of which a pedigree is available. The resulting matrix is known in the literature as matrix H^16,17, this matrix is given by:

where, A _gg is a submatrix of A for genotyped animals, G _a = β G + α;β and α are obtained by solving the system of equations:

The QRH model is given by:

where Varu=σH2H, the rest of the terms as described above.

Single-step GBLUP regression (ssGBLUP) model

The ssGBLUP model is equivalent to the GBLUP model described above with the difference that the genomic relationship matrix G is replaced with the extended genetic relationship matrix H, it is assumed that Varu=HσH2.

Cross-validation

The predictive capacity of the models was evaluated by cross-validation, which was performed as follows. The dataset was divided into five groups of the same size {S1,S2,…,S5}, 80 % of the data was used for training of the model, the remaining 20 % for validation. For example, {S2} is used as a validation group and the set {S1,S3,…,S5} for training of the model. The models were fitted using the training set, and the fitted model was used to obtain predictions for the validation set. This procedure was repeated five times and predictions were obtained for each group. Correlations between observed and predicted phenotypes were calculated and averaged for the test sets¹⁸. Note that because these are actual values, the true breeding values are not known, but only the observed phenotypes are available, the fitted model provides predictions for breeding values and predictions of other fixed and environmental effects, with which a prediction of the phenotype is obtained, which is contrasted with the true value of the phenotype.

Simulation

In order to evaluate the predictive power of the QR model against GBLUP, an asymmetric data simulation with the presence of outliers was also carried out; the simulation of the present work is analogous to that used by Pérez-Rodríguez et al¹⁰. The main idea is to highlight that the quantile regression model works adequately in the presence of atypical observations, inhomogeneous variances and response variables with responses with asymmetric distribution. In the context of selection, for example, it is not unusual to have asymmetric distributions for phenotypes due to the process itself, since, if one selects for some characteristic Y, and if there is in addition to this another characteristic of interest O, then the conditional distribution of Y |O>o¹⁹ is asymmetric. In the context of genomic selection, it is also common to find subsets of observations that differ significantly from the rest and these observations could be considered atypical. Montesinos-López et al²⁰ proposed a model with Laplace errors and showed that it predicts adequately even in the presence of outliers, the proposed model is a special case of the quantile regression model that is studied in the present work. The 9,628 SNPs resulting from the quality control described above for 300 animals were considered, the simulation of the data was carried out considering the linear model:

where i=1,…, 300, with μ=39 for BW, it was assumed that the errors come from a biased normal distribution (SNc) with mean 0, variance σ2 (scale parameter σ) and asymmetry index γ1, that is ei~SNC(0,σ,γ1), with σ=1-h2, h2 with a value of 0.35, γ1=2πρ34π-11-2ρ2π-3/2, ρ∈0.950, 0.975, 0.999 were considered, leading to different degrees of positive bias. Only positive values of γ1 were considered since the negative bias is obtained simply by changing the sign of the ei's and therefore the conclusions obtained for the case of positive bias will also be valid for the negative case^21,22. Fifty markers with non-zero effect were fixed, simulating them from a normal distribution with mean 0 and variance 1-h2/50, the rest of the markers were set at 0; the positions of the sampled markers were taken at random. To introduce outliers in the phenotypes, a certain proportion of the residues o ei~SNC(0,3,γ1) were randomly generated, two proportions were considered, 5 and 10 %, so samples from a mixture of two components of biased normal distributions were taken. Six datasets were generated, three different asymmetry coefficients 0.950, 0.975, 0.999 with their two alternatives of outlier proportion 5 % and 10 %. The asymmetric normal distribution has been used in genomic prediction²² and its use in channeled selection has also been suggested²³. Once the data were generated, the QRM model was fitted with θ=0.25, 0.50, 0.75 to compare it with GBLUP. The selection of quantiles was made according to Nascimento et al¹¹, who consider these three possibilities when the distribution of phenotypes is asymmetric θ∈{0.25,0.75} or when the distribution is symmetric θ 0.50, since our fundamental interest in this work focused on the modeling of possibly asymmetric data and with the presence of outliers. The selection of the parameters was also made for computational convenience since the fitting of the model is done by using intensive computational techniques based on Markov chain Monte Carlo, as mentioned in the section on software and fitting of the models. For each analysis, the correlation between true and estimated β's, the correlation between true Xβ and estimated Xβ^ signals and the component of variance associated with the residuals for each model, which is a way to evaluate the goodness of fit of the models, were calculated. The Deviation Information Criterion (DIC) was also considered, which can be used to select candidate models; models with lower DIC are preferred to models with higher DIC²⁴.

Software and model fitting

The quantile regression models were fitted using a computational strategy similar to that described by Pérez-Rodríguez et al¹⁰. Adaptations of algorithms to include fixed and random effects do not present great computational difficulty. The codes for the fitting of the models were developed in the programming languages R²⁵ and C. The codes for the fitting of the models were organized in such a way that they can be easily run from the statistical software R and are available by requesting them to the first author of the present study. In all cases, three quantiles were selected, θ={0.25, 0.50, 0.75}. The GBLUP and ssGBLUP models were fitted with the BGLR library of R²⁶.

Real data

Tables 2, 3, and 4 show the results of the experiment conducted with BW, WW, and YW data from a Braunvieh cattle population, evaluated under two scenarios 1) with marker information only, and 2) marker and pedigree information. In general, the highest correlations between observed and predicted values were obtained with QR, except for BW, where the correlations of GBLUP and ssGBLUP were higher than those obtained with QRM and QRH, however, the correlations of QRM θ=0.75 and QRH θ=0.75 were close to those obtained with GBLUP and ssGBLUP (0.7902 vs 0.7924), (0.6981 vs 0.7055), respectively. The lowest MSE values were obtained with QRM θ=0.75 and QRH θ=0.75 in the WW characteristic, while in the BW and YW characteristics, the lowest values were obtained with GBLUP and ssGBLUP. The variance components associated with the error obtained with QRM and QRH were lower than those obtained with GBLUP and ssGBLUP. In general, the lowest DIC values were obtained with QRM θ=0.75 and QRH θ=0.75, except for BW with the markers-only scenario, where the lowest DIC was obtained with QRM θ=0.25.

Simulated data

The results of the simulation exercise where QR is compared with GBLUP under different degrees of asymmetry and proportions of outliers are shown in Table 5. Column 2 records the correlations between the “true” marker effects and the estimated marker effects, the correlations obtained with QR were higher than those obtained with GBLUP. Column 3 shows the correlations between the “true signals” and the estimated ones, the highest correlations were obtained with QR. Column 4 records the estimation of the variance components associated with the error and column 5 the DIC values, the lowest values in both columns were obtained with QR θ=0.75.

Real data

The observed and predicted phenotype correlations obtained from cross-validation with actual data were higher when using QRM and QRH in the WW and YW characteristics. For BW, the highest correlations were obtained with GBLUP and ssGBLUP; however, in this study, only three quantiles 0.25, 0.50 and 0.75 were tested, there is evidence in other studies where QR is better than GBLUP, as in the case of the work of Nascimento et al⁴, who compared QR with models such as BLASSO, BayesB and RR-BLUP. These authors found a 15.15 % gain in the predictive capacity of QR compared to RR-BLUP, it should be noted that, mathematically, RR-is equivalent to GBLUP, in addition to the fact that the datasets used in this experiment presented asymmetry.

The values of the mean squared error (MSE) measure the average of the squared error, that is, the difference between the estimator and what is estimated, so low values are preferred; the MSE averages of QRM and QRH were lower than those obtained with GBLUP and ssGBLUP only for WW. The residual variance estimator is an indication that how well or poorly the model fits the observed data, low values are preferred; the smallest variance components of the error were obtained with QRM and QRH for the three characteristics analyzed. Finally, the DIC value is used to select candidate models and, like MSE and error variance components, low values are preferred. The lowest DIC values were obtained with QRM θ=0.75 and QRH θ=0.75, except in the marker-only scenario and BW, where the lowest DIC was obtained with QRM θ=0.25. The mean squared error, the residual variance and the DIC are values that help to choose the best fit model. When examining these values together, it is observed that QRM and QRH are better in some of them, while in others they are not, that is, QR has a performance equal to or greater than GBLUP and ssGBLUP; although it should be noted that only three quantiles were tested and that QR has advantages when used in asymmetric data and outliers, for this case there were only outliers and the distributions did not present asymmetry. Mendes et al²⁷ compared QR with the Bayesian method of LASSO (BLASSO), these authors reported a 6.7 % and 20.0 % increase in accuracy and considered quantiles 0.15 and 0.45 in the evaluation of carcass yield and bacon thickness, respectively, however, the characteristics evaluated in their study were asymmetric.

In the analysis of real data, a limitation of the present study is the sample size, which can impact the variability of the parameters estimated with the models and consequently the variability of the predictions, however, all the models were fitted using the same information and therefore the comparison of the predictive capacity of the models is considered reasonable, the ideal would be to have large sample sizes, but, due to economic limitations, this is not always possible. On the other hand, it is currently very common to use prediction models in which the number of phenotypic records is smaller than the number of predictors (SNPS), that is n≪p, even in this context, numerous studies have shown that Bayesian methods provide sophisticated tools that allow deriving reasonable predictions as long as the regularization parameters are selected properly, for example using cross-validation methods^28-30.

Simulated data

In the simulated data experiment, the correlations between “true” marker effects and estimated effects as well as correlations of “true” and estimated signals were higher when QR was used compared to GBLUP. These results are similar to those obtained by other researchers¹⁰, who simulated data with three different coefficients of asymmetry 0.75, 0.95, 0.999 with 5 % and 10 % of outliers and found that the correlations obtained with QR were higher than those obtained with Bayesian ridge regression (BRR), in addition, this pattern was more evident with a greater asymmetry and proportion of outliers. In this study, simulations with asymmetry coefficients of 0.950, 0.975, 0.999 were carried out and the quantiles with which higher correlations were obtained varied between 0.25 and 0.50; the advantage of QR is that different quantiles can be tested, obtaining better results depending on the quantile used, this advantage in the ability to predict the effects of markers and signals has been taken advantage of by other researchers⁴ , who found no trait association using the traditional GWAS model of single SNP, but, when using QR with extreme quantiles such as 0.1, the model was able to find up to 7 SNPs associated with the characteristics studied.

The coefficients of variance of the error indicate how well the proposed model fits the studied data, the smaller this value, the better the fit, the DIC is another value that is used to compare candidate models. Models with a smaller DIC are preferred to models with a larger DIC²⁴. The lowest residual variance estimators and DIC values were obtained with QR θ=0.75, perhaps this is because high asymmetry coefficients 0.950, 0.975, 0.999 were used in the simulation, so therefore a quantile that fits best is expected to be the highest, in this case 0.75.

QR performed equally or better than GBLUP and ssGBLUP to predict growth characteristics BW, WW and YW, the advantages of this method are more noticeable when the data are more biased and present a higher proportion of outliers, as in the case of the simulation experiment.