Estimation of individual animal SNP-BLUP reliability using full Monte Carlo sampling

Graphical Abstract Summary: Reliabilities of estimated breeding values from a SNP-BLUP model can be calculated using elements of the inverse coefficient matrix of the mixed model equations. Computation of the inverse is not feasible for large data sets when the model has a residual polygenic (RPG) effect. We developed a full Monte Carlo (MC) samplingbased method for approximating reliabilities in the SNP-BLUP model with an RPG effect. Reliabilities obtained by the full MC approach were compared with the corresponding exact reliabilities obtained by the GBLUP model. The full MC approach provided good approximations to the exact values with only a small upward bias. The full MC approach is computationally efficient even for large data sets.

I n the last decade, genomic selection has become the main source of genetic progress in dairy cattle breeding (Mäntysaari et al., 2020). Routine genomic evaluations in animal breeding are usually performed using genomic relationship-based BLUP (GBLUP) or random regression-based SNP-BLUP models. These 2 models are equivalent and, thus, yield equal EBV and prediction error variances (PEV) at the animal level, regardless of whether or not a residual polygenic (RPG) effect is fitted (Strandén and Garrick, 2009;Liu et al., 2016;Ben Zaabza et al., 2020a). Genetic variation cannot be completely explained by SNP markers because of the incomplete linkage disequilibrium between the SNP markers used and the QTL. The RPG effect models the proportion of genetic variance not captured by the SNP markers.
The PEV can be used to obtain the accuracy of EBV or its square-the reliability (Schaeffer, 2019). The calculation of EBV reliability requires elements from the inverse of the coefficient matrix of the mixed model equations (MME). In GBLUP models, the inverse of the genomic relationship matrix (G) has to be computed as well. The G matrix is a dense matrix of the size of the number of genotyped animals (n) without a known simple structure (Meyer et al., 2013). Thus, the MME and G inverse matrices need to be explicitly computed with a cubic computational cost O(n 3 ) of the number of genotyped animals. This is not tenable for very large values of n. The SNP-BLUP model requires no G −1 matrix and the MME size is bounded by the number of markers (m). When n > m, the SNP-BLUP model yields fewer equations to solve than the GBLUP model. However, although low-cost SNP arrays with approximately 50,000 genome-wide SNPs are often used today for most livestock species, even a 700,000-SNP chip array can be used (Meuwissen et al., 2013). With such a high-density SNP chips available, calculating the inversion of the MME in SNP-BLUP can therefore present, at some point, a considerable computational burden.
When the RPG effect is included in the SNP-BLUP model, the computational cost of inverting the MME in SNP-BLUP increases from O(m 3 ) to O[(n + m) 3 ]. Reducing the size of the MME can help decrease some of the computational costs. Ben Zaabza et al. (2020a,b) presented a Monte Carlo (MC) approach for the RPG effect in an SNP-BLUP model, in which the number of RPG effects was reduced from n to the number of MC samples. In this article, we describe a full MC-based SNP-BLUP approach that extends the approach in Ben Zaabza et al. (2020a) to reduce the computational requirements in the calculation of SNP-BLUP genomic reliabilities. We compare the performance of the fully MC-based SNP-BLUP approach to GBLUP where the model has both SNP marker effects and RPG effects.
The SNP-BLUP model with an RPG effect can be expressed as (Ben Zaabza et al., 2020a) where y is an n × 1 vector of observations; n is the number of phenotypic records; 1 n is a vector of n ones; μ is the unknown over- where w is the proportion of RPG, Z c is a centered and scaled genotype marker matrix of order n × m, the square matrix L of size n is the Cholesky decomposition of the pedigree-based relationship matrix of genotyped animals A 22 ; that is, A 22 = LL′; and g R is an m n × 1 vector of unknown genetic effects, m n = m + n.
where λ = σ σ e u 2 2 . Denote the coefficient matrix of this MME by C u and its inverse matrix elements as Approximate reliability (r 2 ) for EBV of animal j can be calculated as: where PEVj is the prediction error variance and σ j 2 the genetic variance for the jth individual. Here, PEVj can be calculated as where t j is equal to the jth row in U. The animal genetic variance σ j 2 is the diagonal element j in the MC approximate genomic relationship matrix UU′. Note that instead of using the diagonal elements of UU′ for σ j 2 , it is possible to compute and use the exact value. The exact value is where F j is the inbreeding coefficient of animal j. We applied the method to 2 data sets to demonstrate the efficiency of the proposed approximation. Our first data set (Data1) comprised genomic and pedigree information in Finnish Red dairy cattle. The pedigree included 231,168 animals. Genotypes were available for 19,757 animals. Because SNP-BLUP is computationally less demanding than GBLUP when the number of genotyped animals is greater than the number of SNP markers, the number of markers in Data1 was reduced to 11,729 SNPs by taking every fourth marker. The second data set (Data2) comprised genomic and pedigree information from a multibreed Irish beef cattle carcass conformation evaluation. Genotypes for 50,240 SNP markers were available for 222,619 animals. The pedigree included up to 13.35 million animals. A summary of Data2 is given in Mäntysaari et al. (2017).
Performance examination of the full-MC-SNP-BLUP reliability approximation involved 2 steps. First, a standard GBLUP model was used to calculate exact GBLUP model reliabilities. Second, the full-MC-SNP-BLUP model was used to calculate approximated reliabilities under 10 different MC sample sizes (10,000, 15,000, 20,000, 30,000, 40,000, 50,000, 60,000, 70,000, and 90,000). Reliabilities for the 2 approaches were computed using 3 RPG proportions (20, 50, and 80%). Performance statistics included correlation, maximum difference, and mean-squared error (MSE) between the exact GBLUP reliability and the full-MC-SNP-BLUP approximated reliability. We fitted the linear regression from the true GBLUP reliability on that from the full-MC-SNP-BLUP to investigate inflation (or deflation) in the approximated reliability. Computational performance was based on wall clock time at selected computing steps. Table 1 gives statistics on GBLUP versus full-MC-SNP-BLUP for Data1 and Data2. For Data1, correlations between the reliabilities from the true and approximation methods ranged from 0.997 to 1.000. The correlation was lowest with a small MC sample size and high RPG weight. For example, the lowest correlation (0.997) was observed with an MC sample size of 10,000 and an RPG proportion of 80%. This scenario had the largest MSE (675 × 10 −5 ) and the largest maximum difference (0.14).
The largest relative decrease in MSE was observed when the number of MC samples was increased from 10,000 to 20,000 (Table 1). Increasing the MC sample size beyond 20,000 led to only a slight decrease in MSE and the maximum difference between full-MC-SNP-BLUP and GBLUP model reliabilities. The MSE values, which reflect accuracy or closeness, were consistent with the correlations, which reflect the association. In fact, the higher the observed correlation, the smaller the MSE. The results reaffirmed the association of the MC sampling sizes with the MSE and maximum difference estimates previously reported in Ben Zaabza et al. (2020a). In their study, MC sampling was used only for the RPG effect, and so it was logical that the higher the RPG weight, the greater the inflation in reliability.
We compared the MSE estimates obtained in this study to those reported in Ben Zaabza et al. (2020a). Note that in our full-MC-SNP-BLUP approach, the size of the MME is determined by the number of MC samples, but in the MC-SNP-BLUP by Ben Zaabza et al. (2020a), the MME size was determined by the number of markers plus the number of MC samples. For example, in full-MC-SNP-BLUP with 20,000 MC samples, we generated and made an MME of size 20,000, when the full model size would have been 11,729 SNP markers + 19,757 genotyped animals. In the MC-SNP-BLUP approach, 20,000 MC samples were used to generate only the RPG effect (19,757 genotyped animals), which made an MME of size 31,729. Given that the computing time to invert the MME coefficient matrix increases in cubic terms of its size, we compared the performance of these approaches at equal MME sizes.
The MSE (66 × 10 −5 vs. 81 × 10 −5 ) reached by full-MC-SNP-BLUP (20,000 MC samples) was similar to that obtained by MC-SNP-BLUP (10,000 MC samples) with an RPG weight of 20% for Data1 (Table 1). When the RPG proportion was increased to 50%, the MSE achieved by the full-MC-SNP-BLUP approach (20,000 MC samples) was much lower, 124 × 10 −5 , than that achieved by MC-SNP-BLUP, 1,822 × 10 −5 (10,000 MC samples). The relative difference was even larger when the RPG weight was 80% (188 × 10 −5 vs. 8,280 × 10 −5 ). These results suggest that the full-MC-SNP-BLUP model approximation is less sensitive, in terms of MSE, to changes from low to high RPG weight than the MC-SNP-BLUP approach at the same MME sizes. Table 1 shows the intercept (b 0 ) and slope (b 1 ) of the regression of correct reliability from GBLUP on approximation by full-MC-SNP-BLUP for Data1. An unbiased approximation method would be expected to give a null intercept and a slope equal to 1. Although the regression of GBLUP on full-MC-SNP-BLUP gave estimates of b 0 close to 0 and b 1 close to 1, all scenarios led to inflated reliability because the slopes were slightly greater than 1. The inflation decreased as the MC sample size increased or the RPG weight decreased, which was also observed by MC-SNP-BLUP using the same data (Ben Zaabza et al., 2020a). However, the full-MC-SNP-BLUP approach appeared to produce intercept and slope estimates smaller than those obtained by MC-SNP-BLUP, especially with an RPG weight larger than 20%. For example, in the scenario with an 80% RPG weight, the slope was 1.70 with MC-SNP-BLUP (10,000 MC samples) but reduced to 1.08 with full-MC-SNP-BLUP (20,000 MC samples). This was partly a consequence of comparing equal-size MME results where the number of MC samples used for approximating the RPG effect was larger in full-MC-SNP-BLUP than in MC-SNP-BLUP. On the other hand, the SNP marker effects were accounted for without MC sampling in MC-SNP-BLUP but they were approximated with MC samples in full-MC-SNP-BLUP. However, the results suggest that when the same number of MC samples is used to approximate covariance structure between animals, the simulated genomic-based covariance structure is closer to the true one than the simulated pedigreebased covariance structure because the number of SNP markers is less than the number of RPG effects.
Correlations between the exact GBLUP and approximate full-MC-SNP-BLUP model reliabilities for Data2 are in Table 1. Increasing the number of MC samples increased the correlation between the model reliabilities, which became clearer as the RPG proportion increased. When the number of MC samples was 10,000 or more, the full-MC-SNP-BLUP approach gave high cor- We compared the analysis of Data2 between the full-MC-SNP-BLUP and MC-SNP-BLUP results at about equal MME sizes, as was done for Data1. Because Data2 had 50,240 markers, about 50,000 more MC samples were allowed in the full-MC-SNP-BLUP than in the MC-SNP-BLUP model when comparing the models. The MSE values with full-MC-SNP-BLUP were smaller than those obtained by MC-SNP-BLUP in all scenarios when the RPG proportion was greater than 20%, which was observed for Data1 as well. For example, with Data2, the MSE associated with full-MC-SNP-BLUP (70,000 MC samples) was smaller than observed with MC-SNP-BLUP (20,000 MC samples) using an RPG weight of 80% (757 × 10 −5 vs. 3,142 × 10 −5 ).
The regression of the reliabilities obtained from GBLUP on those approximated by full-MC-SNP-BLUP indicated a higher degree of correspondence than that between GBLUP and MC-SNP-BLUP. This became clearer, the larger the RPG weight. For example, for an MC sample size of 60,000, the slope was 1.16 (intercept −0.14) and 1.18 (intercept −0.15) for RPG weights of 50% and 80%, respectively, when analyzing Data2 using the full-MC-SNP-BLUP method. However, for an MC sample size of 10,000, the slopes (intercepts) corresponding to these RPG weights were 1.26 (−0.24) and 1.70 (−0.63), respectively, for Data2 using the MC-SNP-BLUP method. Nevertheless, the full-MC-SNP-BLUP approach tended to overestimate the reliabilities (data not shown). Furthermore, the lower the exact GBLUP reliability, the larger the overestimation.
The full-MC-SNP-BLUP computations used MC sampling to approximate the denominator σ j 2 . However, it is possible to compute the exact value as stated in the section Full MC-SNP-BLUP Derivations. We tested the exact σ j 2 values in the full-MC-SNP-BLUP model calculations. However, the use of the exact σ j 2 decreased the accuracy (higher maximum difference and MSE), especially for RPG higher than 20%. Similar results were obtained with the MC-SNP-BLUP approach (Ben Zaabza et al. 2020a). Table 2 gives the computing times for different steps in the calculation of reliability using full-MC-SNP-BLUP for Data1 and Data2, respectively. The computing times are averages over the used RPG weights. The full-MC-SNP-BLUP method involved 4 steps: making the MC samples for the marker effect, making the MC samples for the RPG effect, making the MME coefficient matrix, and inverting the MME matrix. Table 2 clearly shows that the computing times increased along with an increase in the number of MC samples, which would be expected, because the number of MC samples increases the number of computations in all the steps.
For Data1, the time required for MC sampling of the marker effects was often lowest among the computing steps ( Table 2). The computing times required for MC sampling of the RPG effect, on the other hand, were more expensive than those of the marker effects. However, both of these steps showed a linear increase in computing time with the number of MC samples. In contrast, the computing time for making the MME is roughly proportional to n n MC ( ) 2 . For instance, increasing the number of MC sample from 10,000 to 50,000 increased the computing time by a factor of 24, which is close to the squared ratio of the MME sizes; that is, computationally to invert the MME than to make the MME when n MC exceeded the number of genotyped animals (19,757). For Data2, the number of genotyped animals n was always higher than the number of MC samples n MC (Table 2). Consequently, making the MME took more time than inverting the MME in all scenarios, accounting for 3% (1.5 min) to 19% (109 min) of the total computing time. Thus, it can be expected that, for large data sets, making the MME can be computationally more demanding than inverting the MME coefficient matrix. Note that the relative differences between the total computing time and the 4 steps shown in Table 2 Table 3. In the analysis of Data1, the total computing time required by GB-LUP was less than by full-MC-SNP-BLUP when the number of MC samples was more than 30,000. For Data2, the full-MC-SNP-BLUP approach always needed less time than GBLUP, and the reduction in total computing time was 81% (2,504 min) even with the largest MC sample size (90,000 MC samples). In GBLUP, most of the computing time was spent in inverting the 2 matrices (G and MME coefficient matrix). These 2 matrix inversions accounted for 28% (71%) of the total computing time for Data1 (Data2). With full-MC-SNP-BLUP, the computing time for matrix inversion depended on the number of MC samples and took at most 36% (12%) for Data1 (Data2). Thus, it is likely that with larger numbers of genotyped animals, the GBLUP approach will become computationally unfeasible due to the time needed for matrix inversion; in full-MC-SNP-BLUP, this time can be kept acceptable by controlling the number of MC samples.
The reliabilities of the SNP-BLUP models can be approximated satisfactorily by using a full MC-based sampling method. The correlation between the correct reliabilities calculated by GBLUP and those approximated by the full-MC-SNP-BLUP model approached unity in all studied scenarios. The computing time required to invert the MME coefficient matrix was roughly in the proportion There were no external sources of funding for this study.
Viking Genetics (Randers, Denmark) and Nordic Cattle Genetic Evaluation (Aarhus, Denmark) are acknowledged for providing the Finnish dairy cattle genotype data. The authors acknowledge the Irish Cattle Breeding Federation (ICBF, Bandon, Cork, Ireland) for providing the beef cattle data. We gratefully acknowledge the very helpful comments by two anonymous reviewers and insightful suggestions by Journal of Dairy Science section editor Andrés Legarra.
The authors state that this research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.  Reported times are averages over 3 analyses using different residual polygenic weights.
2 Making G = during making the genomic relationship matrix; Inverting G = inversion of the G matrix; Making MME = making of the mixed model equations (MME); Inverting MME = inversion of the MME; and Total = total computing time.