# Theoretical accuracy for indirect predictions based on SNP effects from single-step GBLUP

Andre Garcia, Ignacio Aguilar, Andres Legarra, Shogo Tsuruta, Ignacy Misztal, Daniela Lourenco

► To cite this version:

Andre Garcia, Ignacio Aguilar, Andres Legarra, Shogo Tsuruta, Ignacy Misztal, et al.. Theoretical accuracy for indirect predictions based on SNP effects from single-step GBLUP. Genetics Selection Evolution, 2022, 54 (1), 11 p. 10.1186/s12711-022-00752-4. http://hal-03800048

HAL Id: hal-03800048

https://hal.inrae.fr/hal-03800048

Submitted on 6 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Theoretical accuracy for indirect predictions based on SNP effects from single-step GBLUP

Andre Garcia<fnref n="1" />* ID, Ignacio Aguilar², Andres Legarra³, Shogo Tsuruta<fnref n="1" />, Ignacy Misztal<fnref n="1" /> and Daniela Lourenco<fnref n="1" />

**Abstract**

**Background:** Although single-step GBLUP (ssGBLUP) is an animal model, SNP effects can be backsolved from genomic estimated breeding values (GEBV). Predicted SNP effects allow to compute indirect prediction (IP) per individual as the sum of the SNP effects multiplied by its gene content, which is helpful when the number of genotyped animals is large, for genotyped animals not in the official evaluations, and when interim evaluations are needed. Typically, IP are obtained for new batches of genotyped individuals, all of them young and without phenotypes. Individual (theoretical) accuracies for IP are rarely reported, but they are nevertheless of interest. Our first objective was to present equations to compute individual accuracy of IP, based on prediction error covariance (PEC) of SNP effects, and in turn, are obtained from PEC of GEBV in ssGBLUP. The second objective was to test the algorithm for proven and young (APY) in PEC computations. With large datasets, it is impossible to handle the full PEC matrix, thus the third objective was to examine the minimum number of genotyped animals needed in PEC computations to achieve IP accuracies that are equivalent to GEBV accuracies.

**Results:** Correlations between GEBV and IP for the validation animals using SNP effects from ssGBLUP evaluations were ≥ 0.99. When all available genotyped animals were used for PEC computations, correlations between GEBV and IP accuracy were ≥ 0.99. In addition, IP accuracies were compatible with GEBV accuracies either with direct inversion of the genomic relationship matrix (G) or using the algorithm for proven and young (APY) to obtain the inverse of G. As the number of genotyped animals included in the PEC computations decreased from around 55,000 to 15,000, correlations were still ≥ 0.96, but IP accuracies were biased downwards.

**Conclusions:** Theoretical accuracy of indirect prediction can be successfully obtained by computing SNP PEC out of GEBV PEC from ssGBLUP equations using direct or APY G inverse. It is possible to reduce the number of genotyped animals in PEC computations, but accuracies may be underestimated. Further research is needed to approximate SNP PEC from ssGBLUP to limit the computational requirements with many genotyped animals.

**Background**

One of the ways to deal with the ever-increasing number of genotyped animals in single-step genomic best linear unbiased prediction (ssGBLUP) evaluations is to include only animals with valuable information (own and progeny records) in the evaluations, and then compute

indirect predictions (IP) for the remaining young, genotyped animals [1-3]. In future evaluations, when these animals (for instance, young heifers) have a record or progeny, they could be considered in ssGBLUP; and if they are culled, they would not be considered in ssGBLUP. In addition, IP can be a useful tool to provide fast, interim evaluations for young, genotyped animals, and can also serve as a genomic prediction for animals not included in official evaluations (for instance, genotypes sent from foreign countries). Such predictions reduce the time necessary between collecting a DNA sample and

<footnote>*Correspondence: andre.garcia@uga.edu</footnote>

<footnote n="1">¹ Department of Animal and Dairy Science, University of Georgia, Athens, GA 30602, USA</footnote>

<footnote>Full list of author information is available at the end of the article</footnote>
getting predictions on young animals, which allows farmers and artificial insemination (AI) studs to make faster management decisions and thus to reduce rearing costs by culling animals earlier [4, 5]. When genomic BLUP (GBLUP) or ssGBLUP is used for genomic evaluations, effects of single nucleotide polymorphisms (SNPs) are not readily available but they can be easily backsolved from genomic estimated breeding values (GEBV) using formulas as shown by VanRaden [6], Strandén and Garrick [7] and Wang et al. [8]. Once SNP effects are calculated, IP are obtained for young animals as the sum of the SNP effects weighted by the gene content. Most national dairy cattle genomic predictions use this procedure to run periodical evaluations, obtain estimates of SNP effects (although usually by multi-step procedures), and release fast interim predictions based on IP. In the following, and to avoid confusion, we will call GEBV the direct estimate of the genomic breeding value obtained through ssGBLUP, whereas we use the name IP for the estimate of the genomic breeding value obtained as the sum of SNP effects weighted by the gene content.

Typically, in animal breeding programs, the accuracy of predicted breeding values is calculated to help make selection decisions. Henderson [9] showed that accuracies of EBV can be obtained based on the prediction error variance (PEV), and the latter may be obtained by directly inverting the coefficient matrix of the BLUP mixed model equations (MME). When the system of equations becomes too big, it is impossible to invert the coefficient matrix to obtain PEV even with current computing resources. To overcome this limitation, approximations have been proposed and implemented for pedigree-based evaluations [10] and when genomic information is included [11–15]. Thus, the problem of calculating genomic accuracies of GEBV obtained through GBLUP or ssGBLUP has already been addressed. Similarly, it is interesting for producers to have a measure of the accuracies of IP, to make early selection decisions with more confidence.

Strandén and Christensen [16] showed how to calculate accuracies for IP based on the prediction error covariance (PEC) of SNP effects. These authors and, Tier et al. [17] as well, pointed out that the reliability of GEBV depends on allele coding, and that by back-solving SNP PEC from the same model (ssGBLUP) the accuracy of both GEBV and IP are correctly aligned.

Liu et al. [12] explained that the cost of obtaining such reliabilities from SNP-BLUP is smaller because the size of the left-hand side (LHS) of the MME depends mainly on the number of SNPs rather than the number of genotyped animals. Because of the equivalence between SNP-BLUP and GBLUP, it is also possible to obtain the PEC

for SNP effects when using (ss)GBLUP. However, for (ss)GBLUP, the computational cost depends on the number of genotyped animals. Derivations to obtain SNP PEC under the (ss)GBLUP model were described by Gualdron Duarte et al. [18] and Aguilar et al. [19]. In principle, exact computation of SNP PEC requires that the whole matrix of PEC across all genotyped animals is obtained from the inverse of the MME. This can be very costly in time and memory.

Pocric et al. [14] investigated the accuracy of genomic selection under a GBLUP model using the algorithm for proven and young (APY) and showed that when only a small number of eigenvalues from the genomic relationship matrix (GRM) was used, it was sufficient to account for a large portion of the genetic variation. Because the dimensionality of the genomic information is limited [20, 21], it is possible to reduce the number of animals needed to calculate SNP effects and IP [2, 22]. Likewise, the limited dimensionality could also allow for a reduction in the number of animals needed to obtain SNP PEC and accuracies for IP under (ss)GBLUP.

The objectives of this study were to: (1) present equations to compute individual accuracy of IP by back-solving PEC of GEBV from ssGBLUP into PEC of SNP effects, and to investigate the feasibility of this method; (2) test the algorithm for APY in PEC computations; (3) investigate the minimum number of genotyped animals for which the complete PEC matrix needs to be computed, to obtain IP accuracy in large genotyped populations.

# Methods

## Data and model

Data for this study were provided by the American Angus Association (Saint Joseph, MO) and included 230,639 animals in the pedigree and 38,000 post-weaning gain (PWG) phenotypes. Genotypes for 39,774 SNPs after quality control, were available for 60,000 animals born up to 2018. To mimic a real situation, genotyped animals were split into “old” (N=54,533) and “validation” (N=5467), i.e., young, genotyped animals predicted through IP. Validation animals were genotyped animals that were born in 2018 and for which their pedigree, progeny, and own records were omitted from the data. Validation animals were first excluded from ssGBLUP, and IP and their accuracies obtained; then, these IP accuracies were compared to accuracies when these same validation animals were included in ssGBLUP.

The statistical model for PWG was $y = Xb + Wu + e$, where y is a vector of post-weaning gain phenotypes and b is a vector of fixed contemporary group effects; u is the vector of random additive genetic effects, and e is the vector
of random residuals; X and W are the incidence matrices relating y with the effects in b and u, respectively. Genomic evaluations were implemented using single-step GBLUP.

In ssGBLUP, the inverse of the relationship matrix that combines pedigree and genomic information (H⁻¹) was constructed as in Aguilar et al. [23]:

$$H^{-1} = A^{-1} + \begin{bmatrix} 0 & 0 \\ 0 & G^{-1} - A_{22}^{-1} \end{bmatrix}, \quad (1)$$

where G⁻¹ is the inverse of the genomic relationship matrix, A⁻¹ and A₂₂⁻¹ are the inverses of the pedigree relationship matrix for all and genotyped animals, respectively. All three matrices considered inbreeding. The initial genomic relationship matrix was constructed as the type 1 matrix in VanRaden [6]:

$$G_0 = \frac{ZZ'}{2 \sum p_i (1 - p_i)}, \quad (2)$$

where Z is a matrix of gene content centered for twice the allele frequency of SNP i (pᵢ). Allele frequencies were calculated based on the current genotyped population and recalculated for each evaluation. In this study, G was constructed as:

$$G = b((1 - \alpha)G_0 + \alpha A_{22}) + 11'\delta, \quad (3)$$

where α = 0.05 is a blending parameter [6], and δ and b are tuning parameters calculated as in Vitezica et al. [24]:

$$\delta = \frac{1}{n^2} \left( \sum_i \sum_j A_{22i,j} - \sum_i \sum_j G_{i,j} \right) \text{ and } b = 1 - \frac{1}{2}\delta. \quad (4)$$

After tuning and blending steps, G is invertible and compatible with the pedigree relationships.

Since for large-scale genomic evaluations, it becomes unfeasible to directly invert G, the algorithm for proven and young (APY) was proposed by Misztal et al. [25] and Misztal [26] to overcome this limitation. In APY, the genotyped animals are divided into core (c) and non-core (n) animals:

$$G = \begin{bmatrix} G_{cc} & G_{cn} \\ G_{nc} & G_{nn} \end{bmatrix}. \quad (5)$$

And G⁻¹ is calculated as follows:

$$G_{APY}^{-1} = \begin{bmatrix} G_{cc}^{-1} & 0 \\ 0 & 0 \end{bmatrix} + \begin{bmatrix} -G_{cc}^{-1} & G_{cn} \end{bmatrix} M_{nn}^{-1} \begin{bmatrix} -G_{nc} & G_{cc}^{-1} \end{bmatrix} I. \quad (6)$$

With elements of Mnn, the Mendelian error diagonal matrix, obtained for the i th non-core animal as:

$$m_{nn,i} = g_{ii} - G_{ic} G_{cc}^{-1} G_{ci}. \quad (7)$$

The number of core animals in APY can be obtained as the number of largest eigenvalues explaining 98 to 99% of the variance in G, which can be found by the eigenvalue decomposition of G or the singular value decomposition of Z [21]. In this study, the number of eigenvalues explaining 98% and 99% of the variance in G was 11,413 and 15,242, respectively; therefore 15,000 core animals were randomly selected to be used in APY.

Once H⁻¹ is built, the ssGBLUP MME for PWG are:

$$\begin{bmatrix} X'R^{-1}X & X'R^{-1}W \\ W'R^{-1}X & W'R^{-1}W + H^{-1}\sigma_u^{-2} \end{bmatrix} \begin{bmatrix} \hat{\beta} \\ \hat{u} \end{bmatrix} = \begin{bmatrix} X'R^{-1}y \\ W'R^{-1}y \end{bmatrix}, \quad (8)$$

where R = Iσ² is the residual variance and σ²u the additive genetic variance; β and u are the estimates of fixed effects and GEBV, respectively.

## Benchmark GEBV and accuracy

A ssGBLUP evaluation using the complete data (i.e., including validation animals) was run to obtain benchmark GEBV accuracy (ACCgEBV) for validation animals. The ACCgEBV for animal j was calculated based on PEV from the inverse of the LHS of MME (8) as follows:

Let the inverse of the coefficient matrix of the MME (8) be:

$$C^{-1} = \begin{bmatrix} C^{\beta\beta} & C^{\beta u} \\ C^{\beta u} & C^{uu} \end{bmatrix}. \quad (9)$$

$$\text{Then, } \text{ACC}_{\text{GEBV}_j} = \sqrt{1 - \frac{\text{PEV}_j}{\sigma_u^2}}, \quad (10)$$

where PEVj is the diagonal element j in the prediction error variance matrix Cuu.

## Indirect predictions and accuracy

Before calculating IP, SNP effects from ssGBLUP were obtained as described in Wang et al. [8], using the POSTGSF90 program [27]. Recently, Legarra et al. [28] showed that under ssGBLUP, blending and tuning parameters need to be taken into account when backsolving GEBV into SNP effects (a):

$$\hat{a}|\hat{u} = (1 - \alpha)b \frac{1}{2 \sum p_i (1 - p_i)} Z'G^{-1}\hat{u}, \quad (11)$$

where u is a vector of GEBV from an ssGBLUP evaluation with the reduced data that does not include data for
validation animals. Once SNP effects are available, IP can be calculated as IP = Z validation ̂a, which reflects the marker-based predictions [28].

Liu et al. [12] showed how to compute accuracies for IP from a SNP-BLUP model using SNP PEC as follows. Let the inverse coefficient matrix of the SNP-BLUP MME be:

$$
\mathbf{C}_{\mathrm{g}}^{-1} = \begin{bmatrix} \mathbf{C}^{\beta\beta} \mathbf{C}^{\beta g} \\ \mathbf{C}^{\beta g} \mathbf{C}^{gg} \end{bmatrix}. \qquad (12)
$$

$$
\text{Then, } \mathrm{ACC}_{\mathrm{IP}_j} = \sqrt{1 - \frac{\mathrm{PEV}_j}{\sigma_u^2}}, \qquad (13)
$$

where ACCIPj is the accuracy of IP for animal j;
PEVj = zj Cgg zj'; Cgg is the SNP PEC matrix and zj is the
row vector from the Z matrix, that contains the centered
genotypes for animal j. Since SNP-BLUP and GBLUP
are equivalent models, SNP PEC from SNP-BLUP (or ss-
SNP-BLUP) or from (ss)GBLUP are the same. Gualdron
Duarte et al. [18] and Aguilar et al. [19], showed that PEC
of SNP effects can be calculated as follows:

$$
\mathrm{var}(\hat{\mathbf{a}}) = \mathrm{PEC} = \mathrm{var} \left( (1 - \alpha)b \frac{1}{2 \sum \mathbf{p}_i (1 - \mathbf{p}_i)} \mathbf{Z}' \mathbf{G}^{-1} \hat{\mathbf{u}} \right) \quad (14)
$$

$$
= (1 - \alpha)b \frac{1}{2 \sum p_i (1 - p_i)} \mathbf{Z}' \mathbf{G}^{-1} \mathrm{var}(\hat{\mathbf{u}}) \left( (1 - \alpha)b \frac{1}{2 \sum p_i (1 - p_i)} \mathbf{Z}' \mathbf{G}^{-1} \right), \quad (15)
$$

$$
= (1 - \alpha)b \frac{1}{2 \sum p_i (1 - p_i)} \mathbf{Z}' \mathbf{G}^{-1} (\mathrm{var}(\mathbf{u}) - \mathrm{var}(\hat{\mathbf{u}} - \mathbf{u})) \mathbf{G}^{-1} \mathbf{Z} \frac{1}{2 \sum p_i (1 - p_i)} b(1 - \alpha). \quad (16)
$$

Then,

$$
\mathrm{var}(\hat{\mathbf{a}}) = \mathrm{PEC} = (1 - \alpha)b \frac{1}{2 \sum p_i (1 - p_i)} \mathbf{Z}' \mathbf{G}^{-1} (\mathbf{G} \sigma_u^2 - \mathbf{C}^{u_2 u_2}) \mathbf{G}^{-1} \mathbf{Z} \frac{1}{2 \sum p_i (1 - p_i)} b(1 - \alpha). \quad (17)
$$

Therefore,

$$
\mathrm{var}(\hat{\mathbf{a}}) = \mathrm{PEC} = (1 - \alpha)b \frac{1}{2 \sum p_i (1 - p_i)} (\mathbf{Z}'\mathbf{G}^{-1}\mathbf{Z}\sigma_u^2 - \mathbf{Z}'\mathbf{G}^{-1}\mathbf{C}^{u_2u_2}\mathbf{G}^{-1}\mathbf{Z}) \frac{1}{2 \sum p_i (1 - p_i)} b(1 - \alpha). \quad (18)
$$

Note that α and b are blending and tuning parameters, accounted for in PEC computations, and C<sup>u<sub>2</sub>u<sub>2</sub></sup> is part of the inverse of the LHS of MME (8) corresponding to geno-typed animals.

Once SNP PEC is available, accuracy for IP for animal j
(ACC<sub>IP<sub>j</sub></sub>) can be computed as:

$$
ACC_{IP_j} = \sqrt{1 - \frac{(1 - \alpha)b \mathbf{z}_j \mathrm{var}(\hat{\mathbf{a}}) \mathbf{z}'_j}{\sigma_u^2}}. \quad (19)
$$

While the accuracy of IP can be easily obtained with
small datasets, obtaining C<sup>u<sub>2</sub>u<sub>2</sub></sup> becomes impractical in
large-scale evaluations because of the number of geno-
typed animals. To overcome this limitation, the dimension-
ality of genomic information was exploited by using the
APY algorithm to compute a sparser G<sup>-1</sup>. Lourenco et al.
[22] and Garcia et al. [2] showed that correlations between
IP obtained based on SNP effects from all genotyped ani-
mals or only core animals from APY under ssGBLUP were
higher than 0.98, with greatly reduced computing cost
when using only core animals.

**Implementation**

The implementation required changes in three exist-
ing programs from the BLUPF90 software suite [27].
Most of the changes followed those for the computa-
tion of p-values for SNP effects as described in [19]. In
a nutshell, the modifications in BLUPF90 allows stor-
ing the inverse of the LHS of the ssGBLUP MME in
binary format. For that, OPTION snp_var is required.
When this option is also used in POSTGSF90, the
<table><caption>**Table 4** Peak memory requirements for each scenario</caption><thead><tr><th rowspan="2">Scenarios</th><th colspan="3">Peak memory requirement (GB)<sup>a</sup></th></tr><tr><th>BLUPF90</th><th>POSTGSF90</th><th>PREDF90</th></tr></thead><tbody><tr><td>direct</td><td>195.60</td><td>228</td><td>11.6</td></tr><tr><td>apy</td><td>208.50</td><td>238</td><td>11.6</td></tr><tr><td>50K</td><td>182.75</td><td>211</td><td>11.6</td></tr><tr><td>40K</td><td>103.60</td><td>157</td><td>11.6</td></tr><tr><td>30K</td><td>69.40</td><td>113</td><td>11.6</td></tr><tr><td>20K</td><td>26.69</td><td>78</td><td>11.6</td></tr><tr><td>10K</td><td>5.65</td><td>52</td><td>11.6</td></tr><tr><td>5K</td><td>2.50</td><td>42</td><td>11.6</td></tr><tr><td>2K</td><td>0.63</td><td>38</td><td>11.6</td></tr><tr><td>core (15K)</td><td>18.50</td><td>63</td><td>11.6</td></tr><tr><td>hacc (15K)</td><td>17.90</td><td>63</td><td>11.6</td></tr><tr><td>coreula (15K)</td><td>17.90</td><td>63</td><td>11.6</td></tr><tr><td>hacc_ula (15K)</td><td>18.20</td><td>63</td><td>11.6</td></tr></tbody></table>

scenarios are in Table 4. The total time for PREDF90 to calculate IP and ACCIP for the 5467 validation animals was approximately 30 min of which more than 99% was to calculate the accuracies (Eq. (19)). Although the memory for PREDF90 was the same in direct and apy because of the equal number of validation animals to compute IP, the memory needed for BLUPF90 and POSTGSF90 was about 10 GB more in the apy scenario. This is because a few extra temporary matrices and vectors are needed in APY. One expects APY to considerably reduce memory usage in BLUPF90 and POSTGSF90, which is true with the block implementation of APY adapted to the preconditioned conjugate gradient algorithm [30, 31]. When the inverse of the MME is involved such as in the computation of PEC, or variance components estimation, a full matrix with the dimension of the number of genotyped animals should be allocated to receive the elements of APY $G^{-1}$, making the sparsity not exploitable, as shown previously by Junqueira et al. [32]. In addition, the authors point out that the creation of “fill in” elements may increase the amount of computing time necessary in the sparse inversion.

When trying to reduce computing resources by cutting down the number of genotyped animals used in the computation of PEC from 50K to 2K, correlations between ACCGEBV and ACCIP were still 0.99 with as few as 20K

genotyped animals (Fig. 1). Even when correlations between accuracies are high, we need to make sure that the ACCIP is unbiased, and that it is on the same scale as the ACCGEBV. This will ensure that IP and its accuracy can be used as interim evaluations or replacements for GEBV if the number of genotyped animals becomes extremely large and it is desirable not to have young animals in the evaluation. For all the scenarios, the regression coefficient ($b_1$) of ACCGEBV on ACCIP was used to evaluate dispersion and the intercept ($b_0$) was used to check the scale. If there is no dispersion, $b_1$ equals 1, and deviations from one indicate either under or overestimation of ACCIP. Regression coefficient and intercept for each scenario are presented in Figs. 2 and 3. No bias or scale differences were found when all genotyped animals (except for the validation) were used to calculate PEC in the scenarios direct and apy. However, reducing the number of genotyped animals to 20K increased inflation, as $b_1$ dropped from 0.98 to 0.74 (Fig. 2). At the same time, bias increased from 0.02 to 0.25 (Fig. 3), with a shift towards underestimation. Although an ad-hoc scaling factor based on tests with smaller datasets can correct the overall underestimation or overestimation, it is more difficult to reduce the inflation.

Using the number of eigenvalues explaining 99% of the variance of G (i.e., 15K) as the number of genotyped animals to obtain PEC resulted in a correlation of 0.98. Scenario hacc also had 15K core animals, but these were chosen based on high BLUP accuracy, and the correlation still reached 0.97. Correlations of 0.97 and 0.96 were found for the hacc ula and core ula scenarios, respectively, although the pedigree information was halved, and phenotypes reduced. As those 15K core animals are expected to represent most of the chromosome segments segregating in this Angus population, correlations remained high although only 25% of the available genotyped animals were used. The fact that correlations were not higher than 0.99 when using 15K core animals can be explained by the possible collinearity existing among those animals, thus slightly more animals would be required to reach a correlation higher than 0.99 [3]. In spite of the high correlations, using less data resulted in underestimation of ACCIP, which was expected because less data leads to higher prediction error and consequently lower accuracy. In the current study, average ACCIP for the scenarios with 15K genotyped animals to compute PEC ranged from 0.57 to 0.62 and the regression coefficient ranged from 0.60 to 0.69. The memory requirement for BLUPF90 and POSTGSF90 when using 15K genotyped animals was on average 10 and 28% of that when using all available genotyped animals because of the much smaller number of elements in the coefficient matrix of sGBLUP between the two scenarios.

<footnote><sup>a</sup> Real/resident memory (RSS); Linux server (x86_64) equipped with Intel Xeon E5-2470 2.30 GHz processors with 16 cores</footnote>

<footnote>direct: all genotyped animals (N = 54,533) and phenotypes with direct $G^{-1}$;</footnote>

<footnote>apy: all genotyped animals (N = 54,533) and phenotypes with APY $G^{-1}$; 50K-2K: all phenotypes and decreasing the number of genotyped animals from 50K to 2K; core: genotypes for core animals only (N = 15K) and all phenotypes;</footnote>

<footnote>hacc: genotypes for high accuracy animals only (N = 15K) and all phenotypes;</footnote>

<footnote>core ula: genotypes and phenotypes for core animals plus their progeny</footnote>

<footnote>phenotypes; hacc ula: genotypes and phenotypes for high accuracy animals plus their progeny phenotypes</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Correlation between ACC<sub>GEBV</sub> and ACC<sub>IP</sub>. *direct*: all genotyped animals (N = 54,533) and phenotypes with direct **G**⁻¹; *apy*: all genotyped animals (N = 54,533) and phenotypes with APY **G**⁻¹; *50K-2K*: all phenotypes and decreasing the number of genotyped animals from 50K to 2K; *core*: genotypes for core animals only (N = 15K) and all phenotypes; *hacc*: genotypes for high accuracy animals only (N = 15K) and all phenotypes; *coreula*: genotypes and phenotypes for core animals plus their progeny phenotypes; *hacc_era*: genotypes and phenotypes for high accuracy animals plus their progeny phenotypes</figcaption></figure>

As the number of genotyped animals in the PEC computations decreased further, ACC<sub>IP</sub> were underestimated and the difference in scale between ACC<sub>IP</sub> and ACC<sub>GEBV</sub> increased. For instance, *b*₁ was as low as 0.42 and *b*₀ as high as 0.56 when using 2K genotyped animals. Typically, when *b*₁ is lower than 1, the conclusion is that the

predictions are overestimated; however, this is true when *b*₀ is close to 0. When 30K or less genotyped animals were used to compute PEC, the intercept was not 0 and although *b*₁ was lower than 1, ACC<sub>IP</sub> were underestimated rather than overestimated (Table 2). Differences between ACC<sub>IP</sub> and ACC<sub>GEBV</sub> can also be seen in

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Regression coefficient (*b*₁) of the regression of ACC<sub>GEBV</sub> on ACC<sub>IP</sub>. *direct*: all genotyped animals (N = 54,533) and phenotypes with direct **G**⁻¹; *apy*: all genotyped animals (N = 54,533) and phenotypes with APY **G**⁻¹; *50K-2K*: all phenotypes and decreasing the number of genotyped animals from 50K to 2K; *core*: genotypes for core animals only (N = 15K) and all phenotypes; *hacc*: genotypes for high accuracy animals only (N = 15K) and all phenotypes; *core_era*: genotypes and phenotypes for core animals plus their progeny phenotypes; *hacc_era*: genotypes and phenotypes for high accuracy animals plus their progeny phenotypes</figcaption></figure>
number of genotyped animals is of interest for routine applications. Compared to the approach presented by Liu et al. [12] for a SNP-BLUP model, obtaining SNP PEC from ssGBLUP may be challenging because it depends on the number of animals rather than the number of SNPs included in the system of equations, therefore reducing the number of genotyped animals for PEC computations is critical. Methods to approximate PEC could be likewise helpful.

More research is needed to investigate whether SNP PEC computed from a smaller subset of genotyped animals can be used to approximate ACCIP with a larger number of genotyped animals and to account for large proportions of blending or the residual polygenic effect. Such tests could be hard to accomplish because obtaining ACCGEBV based on PEV as a benchmark is not feasible for large datasets, although approximations could be used. In addition, to be able to use smaller subsets of animals, fine tuning is still needed to refine the methods and to define which animals should be used in PEC computations to avoid biases on IP accuracy. Finally, although it is outside of the scope of this study, since SNP PEC accounts for the genomic contributions from ssGBLUP MME, a combination of our approach with existing PEV approximations may be useful to obtain GEBV accuracies for large-scale ssGBLUP evaluations.

## Conclusions

Indirect prediction accuracy can be successfully obtained by computing SNP PEC from the single-step mixed model equations using direct inversion of the genomic relationship matrix or by the APY algorithm. With at least 40K out of 60K genotyped animals included in PEC calculations, robust indirect prediction accuracies can be obtained without dispersion issues. To reduce the computational costs of inverting the left-hand-side of the mixed model equations, SNP PEC can be approximated by using a smaller subset of the genotyped animals. This yields high correlations, but fine tuning is still required to scale accuracies of indirect predictions up to accuracies of GEBV. Using genotyped sires with phenotyped progeny could help mitigate this issue. Further studies are needed to develop SNP PEC approximations and extend it to large-scale genomic data.

## Acknowledgements

The authors thank the American Angus Association and Angus Genetics Inc. for providing the data and financial support for the study.

## Author contributions

AG, DL, and ST designed the study and planned the analysis. IA, AL and DL refined the theory and scaling/tuning parameters. IA programmed the algorithms. AG performed the analysis. The manuscript was initially written by AG

and then completed by DL, IA, AL, ST, and IM. All authors read and approved the final manuscript.

## Funding

This study was supported by the American Angus Association and its subsidiary Angus Genetics Inc. (St. Joseph, MO) and by the Agriculture and Food Research Initiative Competitive Grant no. 2020-67015-31030 from the US Department of Agriculture's National Institute of Food and Agriculture.

## Availability of data and materials

The data used in this study were provided by the American Angus Association and are used under license for the current study, so they are not publicly available, and restrictions apply to their availability. Modified BLUPF90, POSTGSF90 (to save PEC of SNP effects), and PREDF90 (to compute accuracy of IP) are freely available for research purposes on up to 25,000 genotyped animals at http://nce.ads.uga.edu/html/projects/programs/.

## Declarations

### Ethics approval and consent to participate

Data used in this study were obtained from an existing database.

### Consent for publication

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

### Author details

¹Department of Animal and Dairy Science, University of Georgia, Athens, GA 30602, USA. ²Instituto Nacional de Investigación Agropecuaria (INIA), 11500 Montevideo, Uruguay. ³UMR GenPhySE, INRA Toulouse, BP52626, 31326 Castanet Tolosan, France.

Received: 22 March 2022 Accepted: 23 August 2022

Published online: 27 September 2022

## References

1. Lourenco DA, Tsuruta S, Fragomeni BO, Masuda Y, Aguilar I, Legarra A, et al. Genetic evaluation using single-step genomic best linear unbiased predictor in American Angus. J Anim Sci. 2015;93:2653-62.

2. Garcia ALS, Masuda Y, Tsuruta S, Miller S, Misztal I, Lourenco D. Indirect predictions with a large number of genotyped animals using the algorithm for proven and young. J Anim Sci. 2020;98:skaa15.

3. Tsuruta S, Lourenco DAL, Masuda Y, Lawlor TJ, Misztal I. Reducing computational cost of large-scale genomic evaluation by using indirect genomic prediction. JDS Commun. 2021;2:356-60.

4. Wiggans GR, VanRaden PM, Cooper TA. Technical note: Rapid calculation of genomic evaluations for new animals. J Dairy Sci. 2015;98:2039-42.

5. Nicolazzi EL, Durr JW, Wiggans GR. Genomics in the US dairy industry: current and future challenges. Interbull Bull. 2018;53:54-6.

6. VanRaden PM. Efficient methods to compute genomic predictions. J Dairy Sci. 2008;91:4414-23.

7. Strandén I, Garrick DJ. Technical note: derivation of equivalent computing algorithms for genomic predictions and reliabilities of animal merit. J Dairy Sci. 2009;92:2971-5.

8. Wang H, Misztal I, Aguilar I, Legarra A, Muir WM. Genome-wide association mapping including phenotypes from relatives without genotypes. Genet Res (Camb). 2012;94:73-83.

9. Henderson CR. Applications of linear models in animal breeding. Guelph: University of Guelph; 1984.

10. Misztal I, Wiggans GR. Approximation of prediction error variance in large-scale animal models. J Dairy Sci. 1988;71:27-32.

11. Misztal I, Tsuruta S, Aguilar I, Legarra A, VanRaden P, Lawlor T. Methods to approximate reliabilities in single-step genomic evaluation. J Dairy Sci. 2013;96:647-54.
12. Liu Z, VanRaden PM, Lidauer MH, Calus MP, Benhajali H, Jorjani H, et al. Approximating genomic reliabilities for national genomic evaluation. Interbull Bull. 2017;51:75–85.

13. Erbe M, Edel C, Pimentel ECG, Dodenhoff J, Götz KU. Approximation of reliability in single step models using the interbull standardized genomic reliability method. Interbull Bull. 2018;54:1–8.

14. Pocrnic I, Lourenco DAL, Masuda Y, Misztal I. Accuracy of genomic BLUP when considering a genomic relationship matrix based on the number of the largest eigenvalues: a simulation study. Genet Sel Evol. 2019;51:75.

15. Bermann M, Lourenco D, Misztal I. Efficient approximation of reliabilities for single-step genomic best linear unbiased predictor models with the Algorithm for Proven and Young. J Anim Sci. 2021;100:skab353.

16. Strandén I, Christensen OF. Allele coding in genomic evaluation. Genet Sel Evol. 2011;43:25.

17. Tier B, Meyer K, Swan A. On implied genetic effects, relationships and alternate allele coding. In: Proceedings of the 11th world congress on genetics applied to livestock production: 11–16 February 2018; Auckland. 2018.

18. Gualdron Duarte JL, Cantet RJ, Bates RO, Ernst CW, Raney NE, Steibel JP. Rapid screening for phenotype-genotype associations by linear transformations of genomic evaluations. BMC Bioinformatics. 2014;15:246.

19. Aguilar I, Legarra A, Cardoso F, Masuda Y, Lourenco D, Misztal I. Frequen-tist p-values for large-scale-single step genome-wide association, with an application to birth weight in American Angus cattle. Genet Sel Evol. 2019;51:28.

20. Pocrnic I, Lourenco DAL, Masuda Y, Legarra A, Misztal I. The dimensional-ity of genomic information and its effect on genomic prediction. Genet-ics. 2016;203:573–81.

21. Pocrnic I, Lourenco DAL, Masuda Y, Misztal I. Dimensionality of genomic information and performance of the Algorithm for Proven and Young for different livestock species. Genet Sel Evol. 2016;48:82.

22. Lourenco DAL, Legarra A, Tsuruta S, Moser D, Miller S, Misztal I. Tuning indirect predictions based on SNP effects from single-step GBLUP. Inter-bull Bull. 2018;52:48–53.

23. Aguilar I, Misztal I, Johnson DL, Legarra A, Tsuruta S, Lawlor TJ. Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score. J Dairy Sci. 2010;93:743–52.

24. Vitezica ZG, Aguilar I, Misztal I, Legarra A. Bias in genomic predictions for populations under selection. Genet Res (Camb). 2011;93:357–66.

25. Misztal I, Legarra A, Aguilar I. Using recursion to compute the inverse of the genomic relationship matrix. J Dairy Sci. 2014;97:3943–52.

26. Misztal I. Inexpensive computation of the inverse of the genomic relationship matrix in populations with small effective population size. Genetics. 2016;202:401–9.

27. Misztal I, Tsuruta S, Lourenco DAL, Masuda Y, Aguilar I, Legarra A, et al. Manual for BLUPF90 family of programs. 2014. http://nce.ads.uga.edu/wiki/lib/exe/fetch.php?media=blupf90_all2.pdf/ Accessed 16 Aug 2022.

28. Legarra A, Lourenco DA, Vitezica Z. Bases for genomic prediction. 2021. http://genoweb.toulouse.inra.fr/~alegarra/GSIP.pdf/. Accessed 05 Aug 21

29. Ben Zaabza H, Mäntysaari EA, Strandén I. Using Monte Carlo method to include polygenic effects in calculation of SNP-BLUP model reliability. J Dairy Sci. 2020;103:5170–82.

30. Fragomeni BO, Lourenco DAL, Tsuruta S, Masuda Y, Aguilar I, Legarra A, et al. Hot topic: Use of genomic recursions in single-step genomic best linear unbiased predictor (BLUP) with a large number of genotypes. J Dairy Sci. 2015;98:4090–4.

31. Masuda Y, Misztal I, Tsuruta S, Legarra A, Aguilar I, Lourenco DAL, et al. Implementation of genomic recursions in single-step genomic best linear unbiased predictor for US Holsteins with a large number of geno-typed animals. J Dairy Sci. 2016;99:1968–74.

32. Junqueira VS, Lourenco D, Masuda Y, Cardoso FF, Lopes PS, Silva FF, et al. Is single-step genomic REML with the algorithm for proven and young more computationally efficient when less generations of data are pre-sent? J Anim Sci. 2022;100:skac082.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in pub-lished maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

<figure><img src="image_5.png" /></figure>