An effective hyper-parameter can increase the prediction accuracy in a single-step genetic evaluation

The H-matrix best linear unbiased prediction (HBLUP) method has been widely used in livestock breeding programs. It can integrate all information, including pedigree, genotypes, and phenotypes on both genotyped and non-genotyped individuals into one single evaluation that can provide reliable predictions of breeding values. The existing HBLUP method requires hyper-parameters that should be adequately optimised as otherwise the genomic prediction accuracy may decrease. In this study, we assess the performance of HBLUP using various hyper-parameters such as blending, tuning, and scale factor in simulated and real data on Hanwoo cattle. In both simulated and cattle data, we show that blending is not necessary, indicating that the prediction accuracy decreases when using a blending hyper-parameter <1. The tuning process (adjusting genomic relationships accounting for base allele frequencies) improves prediction accuracy in the simulated data, confirming previous studies, although the improvement is not statistically significant in the Hanwoo cattle data. We also demonstrate that a scale factor, $\alpha$, which determines the relationship between allele frequency and per-allele effect size, can improve the HBLUP accuracy in both simulated and real data. Our findings suggest that an optimal scale factor should be considered to increase prediction accuracy, in addition to blending and tuning processes, when using HBLUP.

1 Introduction

Genomic prediction can achieve a more accurate prediction of additive genetic values at an early life stage, compared to the conventional pedigree-based prediction. Genomic prediction has been applied to a broad range of disciplines, including animal breeding (Hayes et al., 2009) and human disease risk prediction (Abraham et al., 2016; Inouye et al., 2018; Khera et al., 2018). The accuracy of genomic prediction is important, which depends on several factors such as marker density, linkage disequilibrium (LD) between the quantitative trait loci (QTLs) and markers, the sample size of reference, the heritability of the trait, the number of QTLs, and the distribution of QTL effects. The prediction accuracy is also determined by the method used (Yan et al., 2022).

Genomic prediction requires genotypic information for both discovery and target samples. Genome-wide single nucleotide polymorphisms (SNPs) are typically used to estimate the genomic relationship matrix (GRM) for the genotyped samples so that breeding values (in livestock) can be estimated for the target samples, given the phenotypic information of discovery samples (De Los Campos et al., 2009; VanRaden et al., 2009). In many cases, we may have individuals with useful phenotypic information that are not genotyped, but they may be linked with genotyped samples through a pedigree, i.e., missing genotype data. To address this problem, a single-step genomic best linear unbiased prediction (ssGBLUP) method was introduced, in which phenotypic information on both genotyped and non-genotyped individuals in the pedigree can be used simultaneously to maximise the prediction accuracy of genotyped target individuals (Legarra et al., 2009; Christensen and Lund, 2010; Christensen et al., 2012; McWhorter et al., 2022).

SsGBLUP uses an H-matrix that is a harmonised matrix of a pedigree-based numerator relationship matrix (NRM) and a GRM; therefore, we will use the term H-matrix best linear unbiased prediction (HBLUP). The H-matrix allows us to use the information of non-genotyped individuals in genomic prediction using a data augmentation technique (see (Legarra et al., 2009; Misztal et al., 2009) and Legarra et al., 2014). HBLUP has been widely used in the genetic evaluation of livestock and has been employed in the national genetic evaluation program in many countries (Gao et al., 2012; McMillan and Swan, 2017; Brown et al., 2018; Chung et al., 2018; Johnston et al., 2018; Meyer et al., 2018; Teissier et al., 2018; Oliveira et al., 2019; Mäntysaari et al., 2020; Alkhoder and Liu, 2021). There are numerous studies reporting that HBLUP outperforms traditional GBLUP (Baloche et al., 2014; Gao et al., 2018; Gowane et al., 2019; Mancisidor et al., 2021).

In HBLUP, there are two main hyper-parameters that can determine its performance. First, blending is one of the hyper-parameters that can provide a weighted sum of genomic and numerator relationships, using an arbitrary weight typically ranging from 0.5 to 0.99 (Meyer et al., 2018). This process is essential because it ensures GRM, which is a positive definite matrix, avoids numerical problems in HBLUP (VanRaden, 2008; Legarra et al., 2009). Second, tuning is another important hyper-parameter that can adjust GRM, accounting for the allele frequencies in the base population that are inferred from the information of NRM (Legarra et al., 2009; Misztal et al., 2009; Chen et al., 2011; Vitezica et al., 2011). Note that GRM is typically based on genotyped samples in the last few generations, whereas NRM includes the information of founders in the base population through the pedigree. Third, a scale factor is a novel hyper-parameter for HBLUP to be introduced in this study, which can generate different kinds of GRMs, accounting for the relationship between allele frequency and per-allele effect size, i.e., per-allele effect sizes vary, depending on a function proportional to [p (1 − p)] ^α, where p is the allele frequency (Speed et al., 2012; Speed et al., 2017; Schoech et al., 2019; Momin et al., 2023). Negative $\alpha$ values indicate larger effect sizes for rare variants, and the choice of $\alpha$ may determine the HBLUP accuracy, i.e., an optimal $\alpha$ can increase the accuracy. Bouwman et al. (2017) considered alternative scale factors in GBLUP, which were applied to a dairy cattle dataset. However, they did not test the impact of alternative scale factors in HBLUP; therefore, it is unclear how the scale factor, as a hyper-parameter, determines the HBLUP accuracy.

In this study, we investigate the three hyper-parameters, blending, tuning, and $\alpha$, to assess how they affect HBLUP accuracy, using simulated and real data. There are several tuning methods (Legarra et al., 2009; Chen et al., 2011; Vitezica et al., 2011; Meyer et al., 2018) among which we test two most frequently used approaches, i.e., methods by Chen et al. (2011) (Chen et al., 2011) and Vitezica et al. (2011) (Vitezica et al., 2011), referred to as tune = 1 and 2 in this study. For blending, we investigate a wide range of weighting factors ($\theta$) to assess the performance of HBLUP. In the analyses, we use the direct Average Information algorithm (Lee and Van Der Werf, 2006; Yang et al., 2011) that is robust to the numerical problem caused by non-positive definite GRM so that we can assess all kinds of weighting factors in blending, including $\theta$ = 1. We also assess HBLUP performance, varying the scale factor, ranging from $\alpha$ = −1.5 to 1.5, in the estimation of GRM. We consider the three hyper-parameters simultaneously to obtain optimal values for blending, tuning, and $\alpha$, using a grid search method (Bergstra et al., 2011). Then, the performance of HBLUP with optimal values is compared to performances with less optimal values.

2 Material and methods

2.1 Simulated data

QMSim software (Sargolzaei and Schenkel, 2009) was used for simulation since it can efficiently generate a large-scale dataset including genotypic and pedigree information. We simulated three different scenarios that differed in terms of the effective population size, mating design, and family structure. Two different effective population sizes are determined at 100 and 1,000 individuals with 100 generations in order to mimic livestock (a half-sib design) and human (a full-sib design) populations.

I. The historical population consists of 100 generations. For the initial 95 generations, the effective population size ($N_e$) keeps fixed at 100 individuals, consisting of 50 female individuals and 50 male individuals. Two offspring are generated with random selection and random mating of parents. In the following five generations (95th–100th), the number of progenies gradually increased to 1,000. In the last generation of the historical population (the 100th generation), we randomly select 50 male individuals and 500 female individuals as the founders, and each male individual is mated with ten female individuals and each female individual produced two offspring (i.e., a half-sib design). The current population consists of five generations with 1,000 offspring in each generation (101–105th generations), which is used for the main analyses. The details of applied parameters in the simulation of genotypic and pedigree data are listed in Table 1. The steps to simulate the historical and current populations are illustrated in Supplementary Figure S1.

II. In the second simulation scenario, $N_e$ = 1,000 is used (500 female individuals and 500 male individuals) with a historical population of 100 generations. The population size for each generation in the historical population with 100 generations is constant (N = 1,000). In the subsequent five generations (101st–105th), each male individual is mated with one female individual and each female individual produced two offspring (i.e., a full-sib design), and 1,000 offspring were generated in total. Thus, the founder population size is 1,000.

III. In the third scenario, $N_e$ and the number of generations in the historical population are the same as in the first scenario ($N_e$ = 100 with 100 generations). However, In the last generation of the historical population (100th) and the subsequent five generations (101st–105th), the mating design and family structure are the same as the second scenario, i.e., one male individual is mated with one female individual to produce two progeny per mating (full-sib design), producing 1,000 offspring in total in each generation.

TABLE 1

TABLE 1. Parameters of historical population and genotyping data simulation in the first scenario using QMSim software.

In order to simulate the phenotypes of a complex trait based on the simulated genotyped data, we used a model,

$y_i={\mathbf{Z}}_{\mathit{i}}\mathbf{u}+e_i(1)$

where $y_i$ is the phenotypic value, ${\mathbf{Z}}_{\mathbf{i}}$ is the vector of SNP genotypes and $e_i$ is the residual effect for the $i^{th}$ individual, and u is the vector of SNP effects. In this phenotypic simulation, we randomly selected 1,000 SNPs as causal variants, and u was drawn from a normal distribution such that the mean and variance of the genetic effects are $mean\left({\mathbf{Z}}_{\mathit{i}}\mathbf{u}\right)=\mathbf{0}$ and $var\left({\mathbf{Z}}_{\mathit{i}}\mathbf{u}\right)=h^2$. The residual effects were generated from a normal distribution with mean = 0 and variance $=1-h^2$. In the phenotypic simulation, the SNP effects, u, are scaled by [2p (1 − p)]^α, considering a non-negligible relationship between allele frequency and per-allele effect size (Speed et al., 2012; Speed et al., 2017; Schoech et al., 2019; Momin et al., 2023), which is a function of alpha ranging from −1.5 to 1.5 in the simulation.

In the HBLUP analysis, for three simulation scenarios, it is assumed that the pedigree information is available for the last five generations (101–105th generations), and the genotypic information is available for the individuals from the last two generations (104–105th generations), noting that the sample size in each of the last five generations is 1,000. Furthermore, it is noted that the phenotypes are available for all individuals. We conducted 3,000 replicates of the simulations under three different scenarios with specified simulation parameters. By running multiple replicates, we were able to estimate the variance and uncertainty in the results and obtain a more accurate assessment of the effects of different factors on the population. Replicating the simulation multiple times is a common practice in simulation studies as it can increase the reliability and validity of the results by reducing the impact of chance events and providing a more robust assessment of the effects of the factors being studied.

2.2 Real data

2.2.1 Hanwoo cattle data

In this study, we applied statistical analyses to genotypic and phenotypic data from Hanwoo beef cattle. The total number of animals with pedigree information was 84,020, and among them, 13,800 animals were genotyped for 52,791 genome-wide SNPs, and 25,502 animals were recorded for their phenotypes. The number of animals available for both genotypic and phenotypic information was 9,072. The following criteria were applied for quality control (QC) using PLINK: minor allele frequency below 0.01 (MAF), filtering SNPs with a call rate lower than 95% (GENO = 0.05), individual missingness more than 5% (MIND = 0.05), and Hardy–Weinberg Equilibrium p-value threshold lower than 1e-04 (HWE). After QC, the number of individuals did not change, and the SNPs number was 42,795. The Hanwoo beef cattle data included five carcass traits: carcass weight, eye muscle area, back fat thickness, marbling score, and adjusted 12 months weight. The total number of animals with non-missing records for each carcass trait with and without genotypic information can be seen in Table 2.

TABLE 2

TABLE 2. The number of individuals available for phenotypes with and without genotypic information for five carcass traits in the Hanwoo cattle dataset.

In the HBLUP analysis for the Hanwoo cattle data, animals available for phenotypes and genotypes ($N_{g,p}$) (see Table 2) are randomly divided into five groups. In a 5-fold cross-validation, one of the five groups is selected as the target dataset, and the remaining groups are used as the discovery dataset, which is repeated five times, and the averaged phenotypic prediction accuracy is calculated. The technical details of the training and validation of HBLUP can be seen in Figure 1.

FIGURE 1

FIGURE 1. A diagram showing the experimental designs and how to select the target and discovery samples for simulated and Hanwoo cattle datasets. In the simulated dataset, the number of founders depends on the simulation scenarios ($f_n$ = 550, 1,000, and 550 for simulation scenarios 1, 2, and 3). The sample size in each generation ($G_i$) is 1,000. Therefore, the sample size in the whole population is $N_{all}$ = ${\displaystyle \sum _{i=1}^N}{G}_i+f_n$. The sample sizes of target and discovery samples are denoted as $N_t$ and $N_d$. In Hanwoo cattle data, the phenotypic and genotypic information is partly missing. The number of animals without genotype and phenotype ($N_{ng,np}$), animals without genotype but with phenotype ($N_{ng,p}$), animals with genotype but without phenotype ($N_{g,np}$), and animals with both genotype and phenotype ($N_{g,p}$) are shown in the diagram. $N_g$ is the total number of genotyped animals. In HBLUP, for the animals with both genotype and phenotype ($N_{g,p}$), 5-fold cross-validation is applied, and each fold is selected as the target dataset ($N_t$), and the remaining animals with phenotypes are used as the discovery samples ($N_d$). The best linear unbiased predictions for the phenotypes of the target samples are obtained. In order to calculate the prediction accuracy, we used Pearson’s correlation coefficients between the true and predicted phenotypes for the target samples. It is noted that the target dataset is selected from the last generations (offspring) and should be predicted by the previous generations (discovery population).

2.3 Estimating NRM, GRM, and HRM

2.3.1 Numerator relationship matrix

NRM is denoted as A which is estimated based on the pedigree and has been used in Henderson’s mixed model equation (Henderson, 1975) to obtain estimated breeding values. Following Legarra et al., 2014, A matrix can be formulated as follows.

$\mathit{A}=\left[\begin{array}{cc}{\mathbf{A}}_{11}& {\mathbf{A}}_{12}\\ {\mathbf{A}}_{21}& {\mathbf{A}}_{22}\end{array}\right](2)$

Where ${\mathbf{A}}_{11}$ and ${\mathbf{A}}_{22}$ denote the numerator relationships for the groups of non-genotyped and genotyped individuals, and ${\mathbf{A}}_{12}$ and ${\mathbf{A}}_{21}$ are the numerator relationships between non-genotyped and genotyped individuals.

2.3.2 Scale factor ($\mathit{\alpha }$) and GRM

Following Momin et al., 2023, the variance of the $i^{th}$ genetic variant (v_i) can be expressed as a function of the allele substation effect ($u$) and the allele frequency ($p_i$), which can be written as

$Var\left(v_i\right)=2p_i\left(1-p_i\right){\gamma }_i^2={\left[2p_i\left(1-p_i\right)\right]}^{1+2\alpha }\times {u_i}^2(3)$

where ${\gamma }_i=u_i\times {\left[2p_i\left(1-p_i\right)\right]}^{\alpha }$ is the allele effect size ($u_i$) that can vary, depending on the allele frequency and the scale factor, $\alpha$ (Speed et al., 2012; Speed et al., 2017), which can be explained by evolutionary forces such as selections, mutations, immigrations, and genetic drift. In the classical model [36], $\alpha$ is assumed to be zero for all traits. Another widely used $\alpha$ value is $\alpha$ = −0.5, assuming that the genetic variance of the causal variant has a uniform distribution across the minor allele frequency spectrum. However, there have been reports that optimal $\alpha$ values vary, depending on traits and populations (Speed et al., 2012; Speed et al., 2017; Momin et al., 2023). Following Speed et al., 2020, the genomic relationship matrix can be formulated as a function of $\alpha$, which can be written as

${\mathit{G}}_{ij}=\frac{1}{d}{\displaystyle \sum _{k=1}^L}\left[\left(x_{jk}-2p_k\right)\left(x_{ik}-2p_k\right)\right]{\left[2p_k\left(1-p_k\right)\right]}^{2\alpha }(4)$

where ${\mathit{G}}_{ij}$ is the genomic relationship between the $i^{th}$ and $j^{th}$ individuals and $L$ is the total number of SNPs, $p_k$ is the allele frequency of the k^th SNP, $x_{jk}$ is the SNP genotype coefficient of the j^th individual at the k^th SNP, and d is the expected diagonals computed as $d=L\bullet \mathbb{E}\left[{\left(x_{ik}-2p_k\right)}^{\mathbf{2}}{\left[2p_k\left(1-p_k\right)\right]}^{2\alpha }\right]$. Eq. 4 is implemented in LDAK software (Speed et al., 2012).

Note that Eq. 4 with $\alpha$ = −0.5 is equivalent to the genomic relationship estimation implemented in PLINK, GCTA, and option 2 in BLUPf90 (VanRaden, 2008; Yang et al., 2010; Misztal et al., 2018), and Eq. 4 with $\alpha$ = 0 is equivalent to option 1 in BLUPf90 (VanRaden, 2008; Misztal et al., 2018).

In the HBLUP analysis, we will vary $\alpha$ from −1.5 to 1.5 to find an optimal $\alpha$ value that can improve the phenotypic prediction accuracy and compare the performance with the conventional HBLUP (with $\alpha$ = −0.5 or 0). It is noted that $\alpha$ value is used to estimate GRM using the restricted maximum likelihood (REML) method.

2.3.3 H-matrix (HRM) best linear unbiased prediction

In the HBLUP analysis, GRM (G) is computed based on genotypic information, and NRM (A) is estimated using the pedigree information of the population. Following Legarra et al., 2009, given estimated G and A (from Eqs 3, 4), the H matrix can be derived as

$\mathbf{H}=\left[\begin{array}{cc}{\mathbf{A}}_{11}+{\mathbf{A}}_{12}{\mathbf{A}}_{22}^{-1}\left(\mathbf{G}-{\mathbf{A}}_{22}\right){\mathbf{A}}_{22}^{-1}{\mathbf{A}}_{21}& {\mathbf{A}}_{12}{\mathbf{A}}_{22}^{-1}\mathbf{G}\\ \mathbf{G}{\mathbf{A}}_{22}^{-1}{\mathbf{A}}_{21}& \mathbf{G}\end{array}\right](5)$

In the HBLUP analysis, the simulated data were divided into two groups; one group included the individuals in the first three generations, and the other group included individuals in the last two generations in the current population (101–105th generations). We used the genotypic information of the last two generations and the full pedigree information across the five generations to estimate the H matrix. In cattle data, animals available for phenotypes and genotypes were considered (see Table 2) to estimate GRM, and then the HRM was estimated using a combination of NRM estimated based on whole pedigree (84,020 individuals) and GRM.

2.3.4 Blending

GRM is typically a non-positive definite matrix. In the process of HBLUP, it is usually required to modify GRM to be positive definite so that it can be inverted without any numerical problem (VanRaden, 2008). This modification method is called “blending” which shrinks the genomic relationships toward the pedigree relationships, using an arbitrary weight, $\theta$, typically ranging from 0.5 to 0.99 (VanRaden, 2008; Meyer et al., 2018; McWhorter et al., 2022). The blended GRM can be written as

${\mathbf{G}}_{blended}=\theta \mathbf{G}+\left(1-\theta \right){\mathbf{A}}_{22}\forall 0\le \theta \le 1(6)$

2.3.5 Tuning

The tuning process adjusts GRM, accounting for the allele frequencies in the base population, using the information from NRM that includes the information of founders in the base population through the pedigree (Legarra et al., 2009; Misztal et al., 2009; Chen et al., 2011; Vitezica et al., 2011; Hsu et al., 2017). The tuned GRM (${\mathbf{G}}_{tuned}$) is computed as

${\mathbf{G}}_{tuned}=\beta {\mathbf{G}}_{blended}+\omega \mathbf{J}(7)$

where J is a matrix with the same size as GRM, all elements are equal to one, and $\omega$ and $\beta$ are tuning parameters that can be used to adjust GRM, accounting for base allele frequencies. In this study, we use the two most frequently used methods to obtain the tuning parameters, $\omega$ and $\beta$. Following Chen et al., 2011, the first method (referred to as tune = 1) computes $\omega ,and\beta$ as

$\omega =\frac{\left({\mathrm{I}}^{\prime }{\mathbf{A}}_{22}\mathrm{I}-{\mathrm{I}}^{\prime }\mathbf{G}\mathrm{I}\right)}{n_2^2}\beta =\frac{\frac{\left[{\displaystyle \sum _{i=1}^n}{\mathbf{A}}_{{22}_{i,i}}-{\mathrm{I}}^{\prime }{\mathbf{A}}_{22}\mathrm{I}\right]}{n_2}}{\frac{\left[{\displaystyle \sum _{i=1}^n}{\mathit{G}}_{i,i}-{\mathrm{I}}^{\prime }\mathit{G}\mathrm{I}\right]}{n_2}}(8)$

where I is an array with the size of $n\times 1$ and all values equal to one. Following Vitezica et al., 2011, the second method (referred to as tune = 2) can be written as

$\omega =\frac{\left({\mathrm{I}}^{\prime }{\mathbf{A}}_{22}\mathrm{I}-{\mathrm{I}}^{\prime }\mathbf{G}\mathrm{I}\right)}{n_2^2}\beta =1(9)$

Please note that Eqs 8, 9 have been implemented in BLUPf90 (Misztal et al., 2018) as the second and fourth tuning options (i.e., TunedG = 2 or 4).

2.4 Linear mixed model

In the analyses, we used a linear mixed model that can be written as

$\mathit{y}=\mathbf{X}\mathit{b}+\mathbf{Z}\mathit{g}+\mathit{e}(10)$

where $\mathit{y}$ denotes a vector of phenotypic value, $\mathit{b}$ is a vector of the (environmental) fixed effects, $\mathit{g}$ is a vector of random additive genetic effect that is distributed based on $N\left(0,\mathbf{H}{\sigma }_g^2\right)$, where $\mathbf{H}$ can be derived from Eq. 5 and ${\sigma }_g^2$ denotes the genetic variance. Both $\mathbf{X}$ and $\mathbf{Z}$ are the incidence matrixes. Finally, the residual effect vector is shown by $\mathit{e}$ distributed as $N\left(0,\mathbf{I}{\sigma }_e^2\right)$ where $\mathbf{I}$ is an identity matrix and ${\sigma }_e^2$ is the residual variance.

We employed the restricted maximum likelihood (REML) method, fitting the $\mathbf{H}$ matrix, to estimate genetic variance and heritability, which is referred to as HREML in this study. The Akaike Information Criterion (AIC) was used to assess the goodness of fitness of the model as $=2P-2\times \ln \left(L\right)$ , where $\ln \left(L\right)$ is the log likelihood from HREML, and $P$ is the number of parameters. Given the estimated variances and heritability from HREML, HBLUP was used to obtain individual genetic values. We used MTG2.22 (Lee and Van der Werf, 2016; Lee et al., 2017) genomic analysis software to perform the HREML and HBLUP methods.

2.5 Grid search to find optimal hyper-parameters

One of the well-known methods to find the best configuration of hyper-parameters is the grid search (LaValle et al., 2014). In the grid search, all possible combinations of hyper-parameters are considered to evaluate the performance of prediction models. We considered two tuning methods and without tuning (Tune = 0, 1, and 2). The blending step size in the grid search is 0.1 from 0 to 1 and 0.02 from 0.9 to 1.0. Meanwhile, the step size for $\alpha$ is 0.1 from −1 and 1.

2.6 Key performance metrics

This study uses critical performance metrics to evaluate the accuracy and effectiveness of the prediction and estimation methods. The specific performance metrics will depend on the specific research question, the data type, and the prediction method’s goals. Using multiple performance metrics to provide a more comprehensive assessment of the model’s performance is common.

2.6.1 Root Mean squared error (RMSE)

In genomic analysis studies, we often use a performance metric called RMSE to see how good a model is at making predictions. It's like a measuring stick to compare the model’s guesses with the real answers. We calculate RMSE by taking the differences between the model’s predictions and the actual values, squaring them, and then finding the average of those squared differences. Finally, we take the square root of that average. RMSE is useful because it's simple to understand and tells us how far off the model’s predictions are from the real answers, on average. The smaller the RMSE is related the better the model is at making accurate predictions, which means the guesses are closer to the real answers.

2.6.2 R-value

In the study of genes and their effects on physical traits, scientists often use a tool called the Pearson correlation coefficient (R-value). This helps them figure out if there’s a connection between the two things they’re studying. If the coefficient is high and positive, that means when one thing goes up, the other thing tends to go up too. If it’s high and negative, that means when one thing goes up, the other thing tends to go down.

2.6.3 Akaike information criteria (AIC)

In genomic analysis studies, researchers use statistical models to understand the relationship between genes and traits. The Akaike Information Criteria (AIC) is a metric that compares different models and determines the best one. It was developed by Hirotugu Akaike in 1974 (Akaike, 1974) and is based on the principle of maximum likelihood, which aims to estimate the parameters of a statistical model that is most likely to have produced the observed data. The AIC value represents the amount of information the model loses when it approximates the true underlying process. A lower AIC value indicates a better model fit and a higher likelihood of accurately predicting new data. AIC is a valuable metric for model selection because it takes into account both the goodness of fit and the complexity of the model. It penalizes models with more parameters, which can help prevent overfitting and improve the generalisability of the model to new data. Furthermore, AIC can be used in a wide range of statistical models, including linear regression, generalized linear models, and mixed effects models. It plays a crucial role in model selection, allowing us to choose the model that best fits the data while avoiding overfitting and ensuring that the model is generalizable to new data.

3 Results

3.1 Simulated data

Figure 2A shows that the tuning process can improve the phenotypic prediction accuracy (referred to as R-value) when using the simulated data, which is a Pearson correlation coefficient between the observed and predicted phenotypes in the target dataset, confirming previous studies. However, it should be noted that the improvement in prediction accuracy between Blend = 0.9 and Blend = 1 is only 0.003, which may be considered relatively small. The tuning process with the first option (tune = 1; Eq. 8) appears to better perform than the second option (tune = 2; Eq. 9) for this simulated data. However, this shows that tuning GRM before blending had a negligible impact on genomic predictions (McWhorter et al., 2022). Furthermore, blending ($\theta$ <1) does not significantly improve the HBLUP accuracy for this simulated data (Figure 2A; Supplementary Figure S2). Figure 2B represents the impact of $\alpha$ value on the HBLUP’s performance, showing that the prediction accuracy increases when $\alpha$ value used in estimating GRM is close to the true $\alpha$ value used in the phenotypic simulation. When varying simulation scenarios (e.g., a small or large effective population size with full-sib designs), a similar result is observed: the phenotypic prediction accuracy improves when applying the tunning process or when using optimal $\alpha$ (Supplementary Figures S3–S6). In addition, Figure 2C shows the importance of $\alpha$ value in decreasing the root mean square error (RMSE) prediction of the HBLUP, and Blend = 1 proposes less RMSE compared with Blend = 0.1 and Blend = 0.5.

FIGURE 2

FIGURE 2. HBLUP accuracy and hyper-parameters. (A) The HBLUP accuracy (R-value) improves when using tune = 1 (Eq. 8) or tune = 2 (Eq. 9). However, blending ($\theta$ < 1) would not increase the accuracy for this simulated dataset. (B) Optimal $\alpha$ values can increase the accuracy, and also (C) can decrease the RMSE indicating that the choice of $\alpha$ is important in HBLUP. We simulated genotypes and phenotypes in 3,000 replications in which simulation parameters of $h^2=$ 0.8, $N_e=100$ for 100 historical generations and a half-sib design (50 male individuals, 500 female individuals) were used. The true $\alpha$ values used in the phenotypic simulation were −0.5 or 0. The error bars are 95% CI over the 3,000 replications.

Mimicking a real dataset in which multiple replicates are not possible, we used a single simulation data to assess the HBLUP accuracy, varying hyper-parameters (Figure 3). All possible configurations of tuning, blending, and $\alpha$ values were evaluated using the grid search method where the prediction accuracy was measured using 5-fold cross-validation (see Methods and Supplementary Figures S7, S8). Figure 3 shows the HBLUP accuracy averaged over 5-fold cross-validation when varying hyper-parameters. The highest phenotypic prediction accuracy was achieved with tune = 1, blend = 1, and $\alpha$ = 0 when using the true $\alpha$ = 0, and with tune = 1, blend = 0.9, and $\alpha$ = −0.5 when using the true $\alpha$ = −0.5 in the simulations (See Figure 3 and Supplementary Figures S9, S10 for average RMSE). This shows that the optimal $\alpha$ values found in the grid search are approximately in agreement with the true simulated values.

FIGURE 3

FIGURE 3. HBLUP accuracy averaged over 5-fold cross validation in a grid search with various configurations of the hyper-parameters, using a single simulation dataset. The best configuration found in the grid search consists of (A) and tune = 1, blend = 0.9, and α = −0.5 when using α = −0.5 in the simulation, and (B) tune = 1, blend = 1 and α = 0 (in estimating GRM) when using α = 0 in the simulation. The population parameters used in the simulation are h² = 0.8, N_e = 100 for 100 historical generations, N_SNPs = 9000, chromosome number = 30 and α = 0 or −0.5. Mimicking livestock population, a half-sib design (50 sires, 10 dams per sire and 2 offspring per dam) was applied to the last 5 generations. Full pedigree across the 5 generations were used in HBLUP. Among 2000 offspring in the last 2 generations, 5 subsets each with a random 400 individuals were used as target datasets in the 5-fold cross validation. To predict for each target dataset, the remaining 5150 (across the 5 generations) were used as the discovery dataset.

3.2 Cattle data

We used pedigree, genotype, and phenotype data of Korean native cattle (Hanwoo), which is a unique and important breed in the beef industry (Kim et al., 2017; Srivastava et al., 2021), to assess the HBLUP accuracy with various hyper-parameters including $\alpha$. We first estimated optimal hyper-parameters that provided the lowest Akaike information criteria ($\mathrm{A}\mathrm{I}\mathrm{C}$) value based on the residual maximum log-likelihood for each trait, using HREML (Figure 4). We observed that $∆AIC$ was not uniformly distributed across different $\alpha$ values, and optimal $\alpha$ values were largely different across five carcass traits (Figure 4A). On the other hand, a blending parameter $\theta$ = 1 provided the lowest $∆AIC$ values for all traits except for EMA ($\theta$ = 0.86), indicating that a blended GRM with $\theta$ < 1 did not increase the goodness of fit when using HREML in general (Figure 4B). Finally, Figure 4C shows that tune = 2 could achieve better goodness of fit, compared with tune = 1 or tune = 0 (i.e., without tuning), in most cases. For BFT and MS traits, tune = 1 and 0 provided the lowest AIC (Figure 4C) although the AIC was not significantly lower than tune = 2 (difference in AIC less than 1). The best-performed hyper-parameters for five traits can be seen in Supplementary Table S1.

FIGURE 4

FIGURE 4. HREML estimation accuracy depends on $\alpha$ estimated in the genotyped samples and HRM. (A) Evaluating the impact of $\alpha$ values on the $∆AIC$ for five different traits of the Hanwoo cattle dataset using HREML in a univariate linear mixed model with different tuning methods and blending coefficients. The Akaike Information Criterion (AIC) was used to show the goodness of fitness of the model as $AIC=2P-2\times \ln \left(L\right)$, where $2\times \ln \left(L\right)$ is the HREML log likelihood, and $P$ is the number of parameters. $∆AIC=AIC-{AIC}_{optimal}$, where AIC is obtained with the corresponding $\alpha$ value at the x-axis and ${AIC}_{optimal}$ is the AIC for the optimal $\alpha$. It is observed that optimal $\alpha$ varies across traits. Whole individuals with available phenotypes were applied in estimating the heritability based on Table 2. (B) A performance comparison between two different blending coefficients (0.5–1) in order to estimate the HRM using HREML with optimal tuning method and optimal $\alpha$ value. (C) The performance of tune = 1 (Eq. 8) compared with the tune = 2 (Eq. 9), without considering tuning in estimating the HRM with the applied optimal blending and $\alpha$ values.

We also used a grid search to assess the performance of all hyper-parameters (Figure 5) in which HBLUP accuracies of all possible configurations of tuning, blending, and $\alpha$ values were evaluated in 5-fold cross-validation. Figure 5 shows the HBLUP accuracy averaged over 5-fold cross-validation when varying $\alpha$, tuning, and blending values for five carcass traits. In Figure 5A, we observed that the accuracy of HBLUP could be considerably increased or decreased, depending on the choice of $\alpha$ values. In contrast, Figure 5B shows that the highest HBLUP accuracy was achieved with a blending parameter $\theta$ = 1 for all traits except EMA ($\theta$ = 0.86), indicating that blended GRM would not improve the HBLUP accuracy in most cases. Finally, Figure 5C indicates that the tuning process would not substantially improve the HBLUP accuracy for all carcass traits in Hanwoo cattle data. In addition, the 2-D landscape of grid search results for HREML estimation accuracy depending on $\alpha$ estimated in the genotyped samples and making HRM can be seen in Supplementary Figure S11. The best configuration of the hyper-parameters for each trait is shown in Supplementary Table S1.

FIGURE 5

FIGURE 5. The performance of HBLUP when (A) varying $\alpha$, (B) blending, and (C) tuning hyper-parameters for five carcass traits. The five carcass traits include carcass weight (cwt), eye muscle area (ema), adjusted 12 months weight (adj-w12), marbling score (ms), and back fat thickness (bft). There are a total of 84,020 animals in the Hanwoo cattle pedigree, of which 9,072 animals have both phenotypic and genotypic records that are randomly divided into five validation groups (Table 2). Each set of the five groups is selected as the target samples, and all the phenotyped animals except the target samples were used as the discovery dataset. This 5-fold cross-validation was used to validate the performance of HBLUP.

4 Discussion

HBLUP or ssGBLUP has been widely used in livestock breeding programs (De Los Campos et al., 2009; VanRaden et al., 2009). The HBLUP method (e.g., BLUPf90) requires hyper-parameters to integrate the information of genomic and pedigree relationship matrices, which should be optimised to increase the accuracy of genomic prediction (Legarra et al., 2009; Chen et al., 2011; Vitezica et al., 2011; Meyer et al., 2018). In this study, we evaluated the performance of HBLUP with various hyper-parameters such as blending, tuning, and scale factor, using simulated and real Hanwoo cattle datasets.

In our simulation scenario, we employed random mating and random selection instead of artificial selection based on phenotypes or estimated breeding values because the purpose of this simulation study is to demonstrate how the novel hyper-parameter, alpha, works in a simplified simulation setting. Nevertheless, we have applied our approach to real cattle data that have been subjected to artificial selection. By doing so, we believe that we have verified the performance of the hyper-parameter in a realistic setting. In both simulation and real data, allele frequencies can be altered significantly due to genetic drift and selection (Falconer and Mackay, 1996; Hartl and Clark, 1997; Lynch and Walsh, 1998).

The scale factor, $\alpha$, can determine the relationship between allele frequency and per-allele effect size. In the simulation, HBLUP accuracy can be the highest when using GRM scaled by the true $\alpha$ value used in the phenotypic simulation, indicating that the choice of $\alpha$ value is important although this has never been considered as a hyper-parameter in HBLUP. In fact, the performance of HBLUP is shown to vary across the carcass traits in the cattle data used in this study, confirming previous studies that reported that optimal α values vary, depending on traits and populations (Speed et al., 2012; Speed et al., 2017; Momin et al., 2023). Importantly, using less optimal $\alpha$ values may decrease HBLUP accuracy significantly, which should be carefully checked before conducting genetic evaluations, emphasising that the scale factor is not less important, compared to other hyper-parameters such as blending and tuning.

In both simulated and cattle data, blending ($\theta <1$) would not really improve the phenotypic prediction accuracy except that of the cattle traits, EMA; the blending of ${\theta }_{optimal}=0.86$ could increase the accuracy. The accuracy would increase more when GRM was blended with higher weights, which is clearly shown in Supplementary Figure S2. This is not totally unexpected because richer information can come from GRM (e.g., Mendelian sampling variance within sibs), and blended GRM may lose some of such information. When the mixed model equation is used for HREML or HBLUP (Henderson, 1953; Misztal et al., 2018), a non-positive definite GRM may cause a numerical problem, for which the blending process is essential. This may be one of the reasons blending has been an important hyper-parameter in HBLUP. However, the direct Average Information algorithm can use a non-positive definite GRM without blending ( $\theta =1$) and there is a method that can provide a positive definite GRM (Momin et al., 2023). In any case, we recommend optimising the blending hyper-parameter as the optimal blending can vary, depending on data, in which $\theta =1$ should also be explicitly evaluated. A non-positive definite GRM means the matrix has one or more negative eigenvalues, which can cause problems in certain computations. For example, when using the GRM in a linear mixed model for genomic prediction, the non-positive definiteness can lead to negative variance estimates, which are not biologically meaningful (Legarra et al., 2009). In addition, non-positive definite matrices can cause numerical instability in various computations, such as matrix inversion and eigenvalue decomposition (Kang et al., 2008). One common cause of a non-positive definite GRM is the presence of genotyping errors, which can lead to negative pairwise genetic distances between individuals. Small sample sizes, genotyping errors, and mean bias from the current GRM method can contribute to this problem (Momin et al., 2023). However, there is a reliable algorithm available to address non-positive definite GRMs, which is the direct average information algorithm (Lee and Van Der Werf, 2006; Yang et al., 2011; Lee and Van Der Werf, 2016).

The tuning process adjusts GRM, accounting for the allele frequencies in the base population, assuming that the founders in the base population are not genotyped but are linked through the pedigree. As expected, the widely used tuning method (tune = 1 (Chen et al., 2011) implemented in BLUPf90 option 2) could improve the prediction accuracy in the simulated data, indicating that the base allele frequencies are correctly accounted for. However, the improvement caused by tune = 1 or 2 was not remarkable in the Hanwoo cattle data. This is probably due to the fact that the pedigree information in the real data is not accurate enough to trace the founders, or the genotypes may capture substantial information about the base allele frequencies.

The grid search benefits include being able to provide reproducible results, being fast to implement, being simple to develop for parallel computing, and being efficient in exploring a low-dimensional hyper-parameter space. Moreover, for the large-scale hyper-parameters search space, there are a large number of other hyper-parameters optimisation methods, such as genetic/evolutionary algorithms, swarm intelligence methods, stochastic/random search techniques, and co-evolutionary algorithms (Dudzik et al., 2021). These methods are able to provide robust performance in exploring the multi-modal search space (Kuyu and Vatansever, 2021).

In conclusion, existing hyper-parameters such as blending and tuning in HBLUP are important in general, and their optimal values or options should be properly sought to achieve a reliable genetic evaluation. Depending on the data, optimal values can vary, and unnecessary or over-parametrised blending or tuning can produce adverse effects on the prediction accuracy. The scale factor, a novel hyper-parameter to be introduced in HBLUP, should be explicitly optimised to increase the prediction accuracy, given that the impact of the scale factor is competitive with other hyper-parameters, blending and tuning. We suggest including the scale factor, $\alpha$, in HBLUP as a hyper-parameter.

Data availability statement

The data analyzed in this study is subject to the following licenses/restrictions: The SNP genotypic data and phenotypic data of Korean Hanwoo cattle used in this study are deposited and available at the digital repository of NIAS, South Korea (https://www.nias.go.kr/). Moreover, the simulated data generated and used in this study is available by https://github.com/a1708192/HBLUP_gridsearch.git. Requests to access these datasets should be directed to SoL, lhyungm@korea.kr (lhyungm@gmail.com https://www.nias.go.kr/.

Ethics statement

The current study was approved by the Animal Care and Use Committee of the National Institute of Animal Science (NIAS), Rural Development Administration (RDA), and Nonghyup Hanwoo Improvement Center, South Korea. SL (lhyungm@korea.kr).

Author contributions

Conceptualization: MN, SL, BT, JW, and SHL. methodology: MN and SHL. Data curation: MN, SL, and BT formal analysis: MN, SL, and MM investigation: MN and SHL. Software: MN, MM, and SHL. Validation: MN visualization: MN resources: SL and BT funding acquisition: SL and SHL. Writing–original draft preparation: MN and SHL. Writing–review and editing: MN, MM, BT, JW, and SHL. Project administration: SL and SHL. Supervision: SHL. All authors contributed to the article and approved the submitted version.

Funding

We acknowledge the support of the Cooperative Research Program for Agriculture Science and Technology Development (PJ01609901) from the Rural Development Administration, Republic of Korea, in sharing the Hanwoo cattle dataset and technical details. This research is also supported by the Australian Research Council (DP190100766).

Conflict of interest

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1104906/full#supplementary-material

References

Abraham, G., Havulinna, A. S., Bhalala, O. G., Byars, S. G., De Livera, A. M., Yetukuri, L., et al. (2016). Genomic prediction of coronary heart disease. Eur. heart J. 37 (43), 3267–3278. doi:10.1093/eurheartj/ehw450

PubMed Abstract | CrossRef Full Text | Google Scholar

Akaike, H. (1974). A new look at the statistical model identification. IEEE Trans. Automatic Control 19 (6), 716–723. doi:10.1109/TAC.1974.1100705

CrossRef Full Text | Google Scholar

Alkhoder, H., and Liu, Z. (2021). Application of a Single-Step SNP BLUP Model to Conformation Traits of German Holsteins. (INTERBULL BULLETIN: Leeuwarden, The Netherlands) 56, 30–40.

Google Scholar

Baloche, G., Legarra, A., Sallé, G., Larroque, H., Astruc, J. M., Robert-Granié, C., et al. (2014). Assessment of accuracy of genomic prediction for French Lacaune dairy sheep. J. Dairy Sci. 97 (2), 1107–1116. doi:10.3168/jds.2013-7135

PubMed Abstract | CrossRef Full Text | Google Scholar

Bergstra, J., Bardenet, R., Bengio, Y., and Kégl, B. (2011). Algorithms for hyper-parameter optimization. Advances in Neural Information Processing Systems 24 (NIPS 2011) (ACM: Granada, Spain).

Google Scholar

Bouwman, A. C., Hayes, B. J., and Calus, M. P. (2017). Estimated allele substitution effects underlying genomic evaluation models depend on the scaling of allele counts. Genet. Sel. Evol. 49 (1), 1–13. doi:10.1186/s12711-017-0355-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, D. J., Swan, A. A., Boerner, V., Li, L., Gurman, P. M., McMillan, A. J., et al. (2018). “Single-step genetic evaluations in the Australian sheep industry,” in Proceedings of the world congress on genetics applied to livestock production, (Armidale, Australia: University of New England).

Google Scholar

Chen, C. Y., Misztal, I., Aguilar, I., Legarra, A., and Muir, W. M. (2011). Effect of different genomic relationship matrices on accuracy and scale. J. animal Sci. 89 (9), 2673–2679. doi:10.2527/jas.2010-3555

PubMed Abstract | CrossRef Full Text | Google Scholar

Christensen, O. F., and Lund, M. S. (2010). Genomic prediction when some animals are not genotyped. Genet. Sel. Evol. 42 (1), 1–8. doi:10.1186/1297-9686-42-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Christensen, O. F., Madsen, P., Nielsen, B., Ostersen, T., and Su, G. (2012). Single-step methods for genomic evaluation in pigs. animal 6 (10), 1565–1571. doi:10.1017/S1751731112000742

PubMed Abstract | CrossRef Full Text | Google Scholar

Chung, K. Y., Lee, S. H., Cho, S. H., Kwon, E. G., and Lee, J. H. (2018). Current situation and future prospects for beef production in South Korea—a review. Asian-Australasian J. Animal Sci. 31 (7), 951. doi:10.5713/ajas.18.0187

PubMed Abstract | CrossRef Full Text | Google Scholar

De Los Campos, G., Naya, H., Gianola, D., Crossa, J., Legarra, A., Manfredi, E., et al. (2009). Predicting quantitative traits with regression models for dense molecular markers and pedigree. Genetics 182 (1), 375–385. doi:10.1534/genetics.109.101501

PubMed Abstract | CrossRef Full Text | Google Scholar

Dudzik, W., Nalepa, J., and Kawulok, M. (2021). Evolving data-adaptive support vector machines for binary classification. Knowledge-Based Syst. 227, 107221. doi:10.1016/j.knosys.2021.107221

CrossRef Full Text | Google Scholar

Falconer, D. S., and Mackay, T. F. C. (1996). Introduction to quantitative genetics. Essex, England: Longman, 254–256.

Google Scholar

Gao, H., Christensen, O. F., Madsen, P., Nielsen, U. S., Zhang, Y., Lund, M. S., et al. (2012). Comparison on genomic predictions using three GBLUP methods and two single-step blending methods in the Nordic Holstein population. Genet. Sel. Evol. 44 (1), 1–8. doi:10.1186/1297-9686-44-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, H., Koivula, M., Jensen, J., Strandén, I., Madsen, P., Pitkänen, T., et al. (2018). Genomic prediction using different single-step methods in the Finnish red dairy cattle population. J. dairy Sci. 101 (11), 10082–10088. doi:10.3168/jds.2018-14913

PubMed Abstract | CrossRef Full Text | Google Scholar

Gowane, G. R., Lee, S. H., Clark, S., Moghaddar, N., Al-Mamun, H. A., and van der Werf, J. H. (2019). Effect of selection and selective genotyping for creation of reference on bias and accuracy of genomic prediction. J. Animal Breed. Genet. 136 (5), 390–407. doi:10.1111/jbg.12420

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartl, D. L., and Clark, A. G. (1997). Principles of population genetics. Sunderland, MA, USA: Sinauer Associates.

Google Scholar

Hayes, B. J., Bowman, P. J., Chamberlain, A. J., and Goddard, M. E. (2009). Invited review: Genomic selection in dairy cattle: Progress and challenges. J. dairy Sci. 92 (2), 433–443. doi:10.3168/jds.2008-1646

PubMed Abstract | CrossRef Full Text | Google Scholar

Henderson, C. R. (1975). Best linear unbiased estimation and prediction under a selection model. Biometrics 31, 423–447. doi:10.2307/2529430

PubMed Abstract | CrossRef Full Text | Google Scholar

Henderson, C. R. (1953). Estimation of variance and covariance components. Biometrics 9 (2), 226–252. doi:10.2307/3001853

CrossRef Full Text | Google Scholar

Hsu, W. L., Garrick, D. J., and Fernando, R. L. (2017). The accuracy and bias of single-step genomic prediction for populations under selection. G3 Genes, Genomes, Genet. 7 (8), 2685–2694.

CrossRef Full Text | Google Scholar

Inouye, M., Abraham, G., Nelson, C. P., Wood, A. M., Sweeting, M. J., Dudbridge, F., et al. (2018). Genomic risk prediction of coronary artery disease in 480,000 adults: Implications for primary prevention. J. Am. Coll. Cardiol. 72 (16), 1883–1893. doi:10.1016/j.jacc.2018.07.079

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnston, D. J., Ferdosi, M. H., Connors, N. K., Boerner, V., Cook, J., Girard, C. J., et al. (2018). “Implementation of single-step genomic BREEDPLAN evaluations in Australian beef cattle,” in Proceedings of the world congress on genetics applied to livestock production (Armidale, Australia: University of New England).

Google Scholar

Kang, H. M., Zaitlen, N. A., Wade, C. M., Kirby, A., Heckerman, D., Daly, M. J., et al. (2008). Efficient control of population structure in model organism association mapping. Genetics 178 (3), 1709–1723. doi:10.1534/genetics.107.080101

PubMed Abstract | CrossRef Full Text | Google Scholar

Khera, A. V., Chaffin, M., Aragam, K. G., Haas, M. E., Roselli, C., Choi, S. H., et al. (2018). Genome-wide polygenic scores for common diseases identify individuals with risk equivalent to monogenic mutations. Nat. Genet. 50 (9), 1219–1224. doi:10.1038/s41588-018-0183-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, S., Alam, M., and Park, M. N. (2017). Breeding initiatives for Hanwoo cattle to thrive as a beef industry–A review study. J. Anim. Breed. Genom 1, 103. doi:10.12972/jabng.20170011

CrossRef Full Text | Google Scholar

Kuyu, Y. Ç., and Vatansever, F. (2021). Advanced metaheuristic algorithms on solving multimodal functions: Experimental analyses and performance evaluations. Archives Comput. Methods Eng., 28 1–3. doi:10.1007/s11831-021-09555-0

CrossRef Full Text | Google Scholar

LaValle, S. M., Branicky, M. S., and Lindemann, S. R. (2004). On the relationship between classical grid search and probabilistic roadmaps. Int. J. Robotics Res. 23 (7-8), 673–692. doi:10.1177/0278364904045481

CrossRef Full Text | Google Scholar

Lee, S. H., Clark, S., and Van Der Werf, J. H. (2017). Estimation of genomic prediction accuracy from reference populations with varying degrees of relationship. PloS one 12 (12), e0189775. doi:10.1371/journal.pone.0189775

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S. H., and Van Der Werf, J. H. (2006). An efficient variance component approach implementing an average information REML suitable for combined LD and linkage mapping with a general complex pedigree. Genet. Sel. Evol. 38 (1), 1–19. doi:10.1186/1297-9686-38-1-25

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, S. H., and Van der Werf, J. H. (2016). MTG2: An efficient algorithm for multivariate linear mixed model analysis based on genomic information. Bioinformatics 32 (9), 1420–1422. doi:10.1093/bioinformatics/btw012

PubMed Abstract | CrossRef Full Text | Google Scholar

Legarra, A., Christensen, O. F., Aguilar, I., and Misztal, I. (2014). Single Step, a general approach for genomic selection. Livest. Sci. 166, 54–65. doi:10.1016/j.livsci.2014.04.029

CrossRef Full Text | Google Scholar

Legarra, A., Aguilar, I., and Misztal, I. (2009). A relationship matrix including full pedigree and genomic information. J. dairy Sci. 92 (9), 4656–4663. doi:10.3168/jds.2009-2061

PubMed Abstract | CrossRef Full Text | Google Scholar

Lynch, M., and Walsh, B. (1998). Genetics and analysis of quantitative traits. Sunderland, MA: Sinauer Associates.

Google Scholar

Mancisidor, B., Cruz, A., Gutiérrez, G., Burgos, A., Morón, J. A., Wurzinger, M., et al. (2021). ssGBLUP method improves the accuracy of breeding value prediction in huacaya alpaca. Animals 11 (11), 3052. doi:10.3390/ani11113052

PubMed Abstract | CrossRef Full Text | Google Scholar

Mäntysaari, E. A., Koivula, M., and Strandén, I. (2020). Symposium review: Single-step genomic evaluations in dairy cattle. J. dairy Sci. 103 (6), 5314–5326. doi:10.3168/jds.2019-17754

PubMed Abstract | CrossRef Full Text | Google Scholar

McMillan, A. J., and Swan, A. A. (2017). Weighting of genomic and pedigree relationships in single step evaluation of carcass traits in Australian sheep. Proceedings of the 22nd Conference of the Association for the Advancement of Animal Breeding and Genetics (AAABG), Townsville, Queensland, Australia, July 2017 22, 557–560.

Google Scholar

McWhorter, T. M., Bermann, M., Garcia, A. L., Legarra, A., Aguilar, I., Misztal, I., et al. (2022). Implication of the order of blending and tuning when computing the genomic relationship matrix in single-step GBLUP. J. Animal Breed. Genet. 140, doi:10.1111/jbg.12734

CrossRef Full Text | Google Scholar

Meyer, K., Tier, B., and Swan, A. (2018). Estimates of genetic trend for single-step genomic evaluations. Genet. Sel. Evol. 50 (1), 1–11.

PubMed Abstract | CrossRef Full Text | Google Scholar

Misztal, I., Legarra, A., and Aguilar, I. (2009). Computing procedures for genetic evaluation including phenotypic, full pedigree, and genomic information. J. dairy Sci. 92 (9), 4648–4655.

PubMed Abstract | CrossRef Full Text | Google Scholar

Misztal, I., Tsuruta, S., Lourenco, D. A. L., Masuda, Y., Aguilar, I., Legarra, A., et al. (2018). Manual for blupf90 family of programs. Athens, Georgia: University of Georgia.

Google Scholar

Momin, M. M., Shin, J., Lee, S., Truong, B., Benyamin, B., and Lee, S. H. (2023). A method for an unbiased estimate of cross-ancestry genetic correlation using individual-level data. Nat. Commun. 14, 722. doi:10.1038/s41467-023-36281-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Oliveira, H. R., Lourenco, D. A. L., Masuda, Y., Misztal, I., Tsuruta, S., Jamrozik, J., et al. (2019). Application of single-step genomic evaluation using multiple-trait random regression test-day models in dairy cattle. J. dairy Sci. 102 (3), 2365–2377. doi:10.3168/jds.2018-15466

PubMed Abstract | CrossRef Full Text | Google Scholar

Sargolzaei, M., and Schenkel, F. S. (2009). QMSim: A large-scale genome simulator for livestock. Bioinformatics 25 (5), 680–681. doi:10.1093/bioinformatics/btp045

PubMed Abstract | CrossRef Full Text | Google Scholar

Schoech, A. P., Jordan, D. M., Loh, P. R., Gazal, S., O’Connor, L. J., Balick, D. J., et al. (2019). Quantification of frequency-dependent genetic architectures in 25 UK Biobank traits reveals action of negative selection. Nat. Commun. 10 (1), 1–10. doi:10.1038/s41467-019-08424-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Speed, D., Cai, N., Johnson, M. R., Nejentsev, S., and Balding, D. J. (2017). Reevaluation of SNP heritability in complex human traits. Nat. Genet. 49 (7), 986–992. doi:10.1038/ng.3865

PubMed Abstract | CrossRef Full Text | Google Scholar

Speed, D., Hemani, G., Johnson, M. R., and Balding, D. J. (2012). Improved heritability estimation from genome-wide SNPs. Am. J. Hum. Genet. 91 (6), 1011–1021. doi:10.1016/j.ajhg.2012.10.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Speed, D., Holmes, J., and Balding, D. J. (2020). Evaluating and improving heritability models using summary statistics. Nat. Genet. 52 (4), 458–462. doi:10.1038/s41588-020-0600-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Srivastava, S., Lopez, B. I., Kumar, H., Jang, M., Chai, H. H., Park, W., et al. (2021). Prediction of Hanwoo cattle phenotypes from genotypes using machine learning methods. Animals 11 (7), 2066. doi:10.3390/ani11072066

PubMed Abstract | CrossRef Full Text | Google Scholar

Teissier, M., Larroque, H., and Robert-Granié, C. (2018). Weighted single-step genomic BLUP improves accuracy of genomic breeding values for protein content in French dairy goats: A quantitative trait influenced by a major gene. Genet. Sel. Evol. 50 (1), 1–12. doi:10.1186/s12711-018-0400-3

PubMed Abstract | CrossRef Full Text | Google Scholar

VanRaden, P. M. (2008). Efficient methods to compute genomic predictions. J. dairy Sci. 91 (11), 4414–4423. doi:10.3168/jds.2007-0980

PubMed Abstract | CrossRef Full Text | Google Scholar

VanRaden, P. M., Van Tassell, C. P., Wiggans, G. R., Sonstegard, T. S., Schnabel, R. D., Taylor, J. F., et al. (2009). Invited review: Reliability of genomic predictions for North American Holstein bulls. J. dairy Sci. 92 (1), 16–24. doi:10.3168/jds.2008-1514

PubMed Abstract | CrossRef Full Text | Google Scholar

Vitezica, Z. G., Aguilar, I., Misztal, I., and Legarra, A. (2011). Bias in genomic predictions for populations under selection. Genet. Res. 93 (5), 357–366. doi:10.1017/S001667231100022X

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, X., Zhang, T., Liu, L., Yu, Y., Yang, G., Han, Y., et al. (2022). Accuracy of genomic selection for important economic traits of cashmere and meat goats assessed by simulation study. Front. Veterinary Sci. 9. doi:10.3389/fvets.2022.770539

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J., Benyamin, B., McEvoy, B. P., Gordon, S., Henders, A. K., Nyholt, D. R., et al. (2010). Common SNPs explain a large proportion of the heritability for human height. Nat. Genet. 42 (7), 565–569. doi:10.1038/ng.608

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J., Lee, S. H., Goddard, M. E., and Visscher, P. M. (2011). Gcta: A tool for genome-wide complex trait analysis. Am. J. Hum. Genet. 88 (1), 76–82. doi:10.1016/j.ajhg.2010.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar