# Single-step SNP-BLUP with on-the-fly imputed genotypes and residual polygenic effects

- Matti Taskinen
^{1}Email authorView ORCID ID profile, - Esa A. Mäntysaari
^{1}and - Ismo Strandén
^{1}

**49**:36

**DOI: **10.1186/s12711-017-0310-9

© The Author(s) 2017

**Received: **9 September 2016

**Accepted: **17 March 2017

**Published: **30 March 2017

## Abstract

### Background

Single-step genomic best linear unbiased prediction (BLUP) evaluation combines relationship information from pedigree and genomic marker data. The inclusion of the genomic information into mixed model equations requires the inverse of the combined relationship matrix \({\mathbf{H}}\), which has a dense matrix block for genotyped animals.

### Methods

To avoid inversion of dense matrices, single-step genomic BLUP can be transformed to single-step single nucleotide polymorphism BLUP (SNP-BLUP) which have observed and imputed marker coefficients. Simple block LDL type decompositions of the single-step relationship matrix \({\mathbf{H}}\) were derived to obtain different types of linearly equivalent single-step genomic mixed model equations with different sets of reparametrized random effects. For non-genotyped animals, the imputed marker coefficient terms in the single-step SNP-BLUP were calculated on-the-fly during the iterative solution using sparse matrix decompositions without storing the imputed genotypes. Residual polygenic effects were added to genotyped animals and transmitted to non-genotyped animals using relationship coefficients that are similar to imputed genotypes. The relationships were further orthogonalized to improve convergence of iterative methods.

### Results

All presented single-step SNP-BLUP models can be solved efficiently using iterative methods that rely on iteration on data and sparse matrix approaches. The efficiency, accuracy and iteration convergence of the derived mixed model equations were tested with a small dataset that included 73,579 animals of which 2885 were genotyped with 37,526 SNPs.

### Conclusions

Inversion of the large and dense genomic relationship matrix was avoided in single-step evaluation by using fully orthogonalized single-step SNP-BLUP formulations. The number of iterations until convergence was smaller in single-step SNP-BLUP formulations than in the original single-step GBLUP when heritability was low, but increased above that of the original single-step when heritability was high.

## Background

The first model to simultaneously combine genomic information with non-genotyped animal information was single-step best linear unbiased prediction (BLUP) [1, 2] or ssGBLUP. When the number of genotyped animals is large, ssGBLUP may become computationally infeasible for practical purposes because it requires the inverses of dense matrices of size equal to the number of genotyped animals, particularly the inverse of the genomic relationship matrix \({{\mathbf{G}}_g^{-1}}\). In addition, matrix \({\mathbf{G}}_g\) can be singular when the number of genotyped individuals exceeds the number of markers. Computational challenges may have been a reason for the slow adoption of ssGBLUP instead of a multi-step approach. A computationally scalable alternative, the algorithm for proven and young (APY), has been suggested [3]. In APY, a sparse \({{\mathbf{G}}_{APY}^{-1}}\) approximation to the \({{\mathbf{G}}_g^{-1}}\) matrix is created by setting a diagonal matrix for a group of individuals. In practice, APY has been able to reduce computational costs significantly when the number of genotyped animals is very large [4, 5]. However, different sets of core animals in APY will give different evaluations, which may affect selection decisions.

An alternative formulation called hereinafter single-step single nucleotide polymorphism BLUP (ssSNP-BLUP) [6] overcomes some of the major computational challenges in ssGBLUP. In particular, there is no need to construct or invert the genomic relationship matrix \({\mathbf{G}}_g\). The original idea in ssSNP-BLUP circumvents the problems in ssGBLUP by predicting or imputing genotypes for non-genotyped animals, and relying on computationally less demanding SNP-based prediction instead of breeding value based prediction. An additional advantage is that the marker effect solutions are easier to use for interim predictions.

The ssSNP-BLUP has some computational challenges as well. A simple implementation for ssSNP-BLUP generates and stores genotypes of all SNPs for all animals. This will lead to very large disk storage and fast reading requirements that will prohibit use of the approach for large populations. An alternative is to make the required genotype imputations “on-the-fly” instead of storing the very large amount of imputed genotypes to file. For practical purposes, the on-the-fly imputation requires a fast computing approach for the imputation step and/or fast convergence of the iterative method.

The original description for ssSNP-BLUP approach presents a wider range of models than ssGBLUP such as the use of a number of different Bayesian SNP model formulations [6]. In ssGBLUP, in contrast to ssSNP-BLUP, it is typical to include residual polygenic (RPG) information to enhance genomic information by including pedigree-based relationships into the genomic relationship matrix. Thus, the genomic relationship matrix is usually “adjusted” with part of the pedigree relationship matrix either to supply more additive relationship information or simply to make the genomic relationship matrix invertible. In the original ssSNP-BLUP formulation, this information has not been included. The computational challenges of ssGBLUP has led to the introduction of several equivalent models (e.g., [7–9]). However, these alternative approaches have had poor convergence by iterative methods [7, 9]. One reason is that the covariance structures have a poorer condition number which is a ratio of the largest and smallest eigenvalues and is used to measure numerical stability [9]. Some of the alternative versions of ssGBLUP have had SNP effects. Because ssSNP-BLUP is equivalent to ssGBLUP and ssGBLUP has many equivalent forms, alternative formulations of mixed model equations (MME) can be derived for ssSNP-BLUP as well.

In this paper, simple block LDL type decompositions [10] of the ssGBLUP relationship matrix are derived to obtain several linearly equivalent MME. This allows derivation and testing of several equivalent ssSNP-BLUP MME that avoid making and storing the imputed genotypes. In this paper, an explicit imputation of genotypes is not needed but instead we apply sparse matrix decompositions to attain pedigree-based regressions of evaluations of genotyped animals on non-genotyped animals using sparse matrix decompositions. Accuracy and iteration convergence of the derived equivalent MME are tested on a small Nordic dairy cattle dataset.

## Methods

For any model, an infinite number of equivalent models exist. In the following, first we recall the concept of equivalent models and how equivalent models can be made by attaching covariance information to the model design matrix. As an example, equivalence of genomic BLUP (GBLUP) and SNP-BLUP is presented. Then, ssGBLUP is recalled, and its covariance structure formulated using LDL decomposition. Equivalent MME are derived where information from the genomic relationship matrix are attached to the model design matrix similarly as shown for GBLUP and SNP-BLUP, resulting in ssSNP-BLUP. Finally, MME having orthogonalized random effects or diagonal covariance structures are derived in order to improve the convergence in iterative methods. A small dataset is used to illustrate performance of the derived equivalent models.

### Linearly equivalent models

*linearly equivalent models*[11–13] if the expected values and the variances of the observations are equal. Thus, models (1) and (2) are equivalent if:

#### Separate equations for fixed and random effects

#### Mixed model equations

#### Equivalent models by splitting the variance matrix

#### Linearly equivalent MME

*linearly equivalent MME*(11) by pre-multiplying the new random effects \(\widetilde{{\mathbf{u}}}\) with matrix \({\mathbf{M}}\), i.e.

The original MME (5) can, thus, be solved from a modified linearly equivalent MME (11) where the number of random effects in \(\widetilde{{\mathbf{u}}}\) could be smaller or larger than in \(\widehat{{\mathbf{u}}}\), the variance structure \(\widetilde{{\mathbf{G}}}\) could be easier to obtain or invert than the original \({\mathbf{G}}\), or the new matrix system could be otherwise numerically more efficient.

#### Presentation for GBLUP and SNP-BLUP

### Single-step SNP-BLUP

In ssGBLUP [1, 2], some individuals have genomic information while some have only pedigree information. The model for ssGBLUP is a special case of mixed effect models, where the between-animal relationships are modeled via the aggregated relationship matrix \({\mathbf{H}}\) [1, 2]. The relationships in \({\mathbf{H}}\) are described by the pedigree-based relationship matrix \({\mathbf{A}}\), and the genomic relationship matrix among genotyped animals by \({\mathbf{G}}_g\). Relationships among non-genotyped individuals are constructed from the pedigree but modified according to the relationships among genotyped animals.

#### Single-step relationship matrix

*imputation operator*that expands the genomic relationship information from the genotyped to the non-genotyped individuals [2, 6, 17]. By inverting the LDL decomposition of Eq. (22), the variance matrix \({\mathbf{H}}\) has a similar decomposition:

#### Equivalent single-step MME

Note that in these equivalent MME, square matrix \(\widetilde{{\mathbf{G}}}_i\) represents the covariance structure for the reparametrized random effects, \({\mathbf{M}}_i\) has a row for each original random effect in \({\mathbf{u}}\) to change the model matrix \({\mathbf{Z}}\), and \(\widetilde{{\mathbf{Z}}}_i\) is the redefined model matrix.

###
*Basic equivalent ssGBLUP MME*

From the modified relationship matrix \(\widetilde{{\mathbf{G}}}_1\) (36), it can be seen that this *basic equivalent ssGBLUP* MME has random effects for non-genotyped animals with variances \({({\mathbf{A}}^{11})}^{-1}\) and for genotyped animals with variances of the adjusted genomic relationship matrix \({\mathbf{G}}_w\). The number of effects in the system is, thus, the same as in the original ssGBLUP.

###
*Basic RPG ssSNP-BLUP MME*

*residual polygenic effects*that can describe effects that the marker effects are unable to model [8]. The third group of random effects are the marker effects as in SNP-BLUP (15) and so, this decomposition (39) can be called

*basic RPG ssSNP-BLUP*MME. Compared to the original ssGBLUP, the second equivalent MME has marker effects in addition to the animal effects.

###
*Expanded RPG ssSNP-BLUP MME*

The third equivalent MME (40) has three groups of effects similar to the second MME (39). The third effect group has, again, the orthogonal marker effects, and so this formulation is a ssSNP-BLUP as well. The first effect group, for the non-genotyped animals has, however, a constant multiplier \(\sqrt{1-w}\). Also, the second group, related through the pedigree relationship matrix \({\mathbf{A}}\), has now effects for all animals, and not just for the genotyped animals. Thus, in this third *expanded RPG ssSNP-BLUP* MME, the non-genotyped animals have two sets of random effects.

###
*Special cases of equivalent ssGBLUP MME*

These three equivalent MME, (36), (39) and (40), will approach the usual animal model when \(w \rightarrow 1\). At the limit (\(w = 1\)), the expanded RPG ssSNP-BLUP MME 3 (40) has clearly the recognizable covariance structure of \({\mathbf{A}}\). The basic equivalent ssGBLUP MME 1 (36) and the basic RPG ssSNP-BLUP MME 2 (39) are models where the genotyped animals act as base animals and the non-genotyped animals are regressed on them.

In the case where \(w = 0\), i.e. there is no adjustment of the genomic relationship matrix and, therefore, no residual polygenic effects, the basic and the expanded RPG ssSNP-BLUP MME, 2 (39) and 3 (40), are essentially the same MME as was derived by Fernando et al. [6]. The main differences are that they have moved the centering term of the marker matrix \({\mathbf{Z}}_m\) into an additional fixed effect, and they proposed to solve the MME using Bayesian regression.

### Efficient implementation

The three linearly equivalent MME based on Eqs. (36), (39), and (40) contain inverted and non-inverted terms of the pedigree relationship matrix \({\mathbf{A}}\). In an efficient set up to solve these MME, all these terms can be expressed, or modified into a form that can be expressed, with sparse matrices or sparse decompositions of sparse matrices. This is expected to give three efficient implementations of these MME. In practice, matrix equations of these MME are assumed to be solved iteratively by the preconditioned conjugate gradient (PCG) algorithm. Then, only a matrix-vector product of the MME coefficient matrix times a vector is performed once every iteration.

#### Sparse matrices and decompositions

*sparsity preserving Cholesky factorization*so that:

*forward or backward substitutions*and emphasizes the importance of avoiding inverting matrices. In other words, \({\mathbf{L}}{\backslash }{\mathbf{y}}\) is the solution \({\mathbf{x}}\) of equation \({\mathbf{L}}{\mathbf{x}}= {\mathbf{y}}\) or can be expressed as solving (\({\mathbf{L}}\),\({\mathbf{y}}\)). Note that the matrix products are carefully nested with parenthesis so that only matrix-vector operations are performed.

#### On-the-fly imputation operation

The derived equivalent formulations in Eqs. (36), (39), and (40) contain imputation operator \({\mathbf{A}}_{imp} = -{({\mathbf{A}}^{11})}^{-1}{\mathbf{A}}^{12}\) (23) in matrices \({\mathbf{M}}_i\) that are needed when operating with the modified model matrix \(\widetilde{{\mathbf{Z}}}_i = {\mathbf{Z}}{\mathbf{M}}_i\) and its transpose \(\widetilde{{\mathbf{Z}}}^{\prime}_i = {\mathbf{M}}^{\prime}_i{\mathbf{Z}}^{\prime}\), and when calculating the original random effects in Eq. (12). When the MME are solved by the PCG iteration algorithm, the core of the algorithm is a multiplication of the so-called direction vector \({\mathbf{v}}\) by the left hand side of the MME (11). In this multiplication, the imputation operator, as part of the MME coefficient matrix, operates either with a part of the vector pertaining to random effects of the genotyped animals (i.e. \({\mathbf{A}}_{imp}{\mathbf{v}}_2\)) or to a vector of the marker effects \({\mathbf{v}}_m\) through the marker matrix (i.e. \({\mathbf{A}}_{imp}{\mathbf{Z}}_m{\mathbf{v}}_m\)). Thus, in the transpose \(\widetilde{{\mathbf{Z}}}^{\prime}_i\) side, the imputation term operates on a vector of size equal to the number of non-genotyped animals (i.e. \({\mathbf{A}}_{imp}^{\prime}{\mathbf{v}}_1\)).

*imputed genomic marker data*term, \(-{({\mathbf{A}}^{11})}^{-1}{\mathbf{A}}^{12}{\mathbf{Z}}_m{\mathbf{v}}_m\), that expands the genomic information from genotyped to non-genotyped animals, can be calculated using:

Note that the actual imputation of the genomic marker information is not needed. The imputation operation is performed only implicitly, “on-the-fly” during the iterative solution without the need to use, for example, disk storage. In the normal imputation process [6], the marker information needs to be calculated for thousands, or even hundreds of thousands marker vectors of genotyped animals, i.e. columns of marker matrix \({\mathbf{Z}}_m\). The predicted marker data matrix contains real numbers and can be very large, and, thus, takes a lot of time and disk space to generate and use.

In the on-the-fly imputation process of genetic effects, however, imputation is an operation on a “projection vector” of the genotyped animals, i.e. a linear combination of the marker vectors. It needs to be performed only twice within each iteration round for each trait. Once for matrix \(\widetilde{{\mathbf{Z}}}_i\) and another time for the transpose \(\widetilde{{\mathbf{Z}}}^{\prime}_i\) in matrix multiplication of MME coefficient matrix of Eq. (11). The imputation operation is also needed once before the iteration when calculating the right-hand-side for the new random effects (\(\widetilde{{\mathbf{Z}}}_i^{\prime}{{\mathbf{R}}}^{-1}{\mathbf{y}}\)) in Eq. (11) and once at the end of the iteration in order to retrieve the original random effects (\(\widehat{{\mathbf{u}}}= {\mathbf{M}}_i\widetilde{{\mathbf{u}}}\)) in Eq. (12).

#### Orthogonalization of random effects

In the SNP-BLUP versions of the derived equivalent ssGBLUP, in Eqs. (39) and (40), the marker effects are orthogonal, i.e. their covariance matrix is diagonal. It turns out that the PCG iteration numbers of these two equivalent ssSNP-BLUP MME are considerably larger than the original ssGBLUP. In Eq. (39), the RPG effects and genomic values predicted by SNPs have colinearity, and in Eq. (40), the RPG and the animal effects for non-genotyped animals are difficult to separate.

The key for maintaining good numerical properties of the original ssGBLUP seems to be to “orthogonalize” the remaining new random effects, too. The remaining variance structures can be orthogonalized by splitting the variance matrices and attaching the two “halfs” into the coefficient matrices \({\mathbf{M}}\) as in Eqs. (6) and (29).

The fourth equivalent MME (59) can be called *orthogonal ssSNP-BLUP* MME, and the fifth MME (60), originating from the expanded RPG ssSNP-BLUP (40), *orthogonal expanded ssSNP-BLUP* MME.

#### Reduction of the number of effects by using ancestors of genotyped animals

Matrix \({\mathbf{A}}_{22}\), as the covariance structure for the genotyped animals, was reparametrized in Eq. (56) using the full pedigree relationship matrix \({\mathbf{A}}\). This reparametrization increases the number of corresponding new random effects from genotyped animals to all animals in the pedigree. However, computations involving \({\mathbf{A}}_{22}\) require only the genotyped individuals and their ancestors. Thus, to reduce the number of extra new effects, \({\mathbf{A}}_{22}\) can be expressed using a smaller pedigree and relationship matrix \(\widehat{{\mathbf{A}}}\) containing the genotyped animals and their ancestors [4].

*reduced orthogonal ssSNP-BLUP*MME, can be derived from Eqs. (59) and (61) as:

### Data

The derived MME were tested using a small Nordic Red dairy cattle dataset and a simple model. The small dataset and the model were partially chosen in order to be able to use direct sparse matrix solutions of the original ssGBLUP to obtain accurate “correct solutions”. The data were deregressed proofs of milk yield that were based on estimated breeding values from the Nordic production trait evaluations by NAV (Nordic Evaluations, Denmark). There were 73,579 animals in the pedigree of which 2885 were genotyped. Genotyped animals together with their ancestors form a smaller pedigree of 6833 animals. The animals had been genotyped with the Illumina Bovine SNP50 Bead Chip (Illumina, San Diego, USA). The analysis used 37,526 SNPs that passed quality control. There were 66,426 non-genotyped and 1222 genotyped animals with phenotypes. Hence, 1663 animals had a genotype but no phenotype. We considered a single trait model and assumed a heritability of 0.5. The genomic data contained one pair of animals with identical genomic marker data and a couple of more near identical pairs that led to problems for the inversion of the genomic relationship matrix \({\mathbf{G}}_g\) without \({\mathbf{A}}_{22}\) adjustment, i.e. in the case \(w = 0\).

### Comparison statistics

The original ssGBLUP with inverse variance matrix \({{\mathbf{H}}}^{-1}\) of Eq. (20) and the six linearly equivalent formulations of the form (29), from Eqs. (36), (39), (40), (59), (60), (61), and (65), were implemented and tested in an Octave [21] environment. Sparse matrix factorizations were based on CHOLMOD routines [22]. Six different weights \(w\) were tested: 0.00, 0.01, 0.10, 0.20, 0.30, and 1.00. Because of the singularity in the inverse of the genomic relationship matrix \({{\mathbf{G}}_g}^{-1}\) with the test data, the original ssGBLUP matrix (20) and the first equivalent formulation (36) were not calculated when \(w = 0\). The derived new formulations were compared against the original ssGBLUP, mainly focusing on efficiency, accuracy, and number of iterations.

Efficiency of the derived equivalent ssGBLUP formulations relies on the sparsity of the pedigree relationship matrices and their decompositions. Inverse matrix operations of these sparse matrices were transformed into solving sparse lower triangular matrix systems. The efficiency of these solving operations depend on the sparsity structure of the matrices, i.e. number of non-zero elements.

Implementations of the derived equivalent ssGBLUP formulation were not yet streamlined for speed and, thus, the execution times were neither optimal, nor comparable. The performance of the formulations is, therefore, tested by comparing the number of iterations of the iterative solution. The purpose was to demonstrate that the iteration counts are comparable to those obtained by the original ssGBLUP. With a larger number of genotyped individuals, the inversion of the genomic relationship matrix in the original ssGBLUP becomes a bottleneck and ssSNP-BLUP formulations, avoiding the inverse, become relatively faster.

Two versions of the genomic marker matrix encodings were tested. The first was VanRaden method 1 matrix (\({\mathbf{Z}}_m^{vR1}\)) where the marker information is centered around the mean of the observed genotyped animals, and scaled [15]. The second was “−1,0,1 encoding” (\({\mathbf{Z}}_m^{-1,0,1}\)) where the 0,1,2 genotypes were assigned values of −1,0,1, then scaled. In both cases, the centered \({\mathbf{Z}}_m\) was scaled by dividing by \(\sqrt{\sum _{i=1}^{m} 2p_i (1-p_i)}\) where \(m\) is the number of markers, and \(p_i\) is the allele frequency of marker \(i\). In case of −1,0,1 coding, the allele frequencies were all \(p_i=0.5\), and the scaling factor was \(\sqrt{\frac{m}{2}}\). In VanRaden 1, the allele frequencies were those of the observed genotypes.

## Results and discussion

### Efficiency

Number of rows (\(N_r\)), columns (\(N_c\)), non-zeros (\(N_z\)), and mean number of non-zeros on row or column (\(M_z\)) by matrix used in MME coefficient matrices in the test case with 73,579 animals

Matrix | \(N_r\) | \(N_c\) | \(N_z\) | \(M_z\) |
---|---|---|---|---|

\({\mathbf{A}}^{-1}\) | 73,579 | 73,579 | 235,669 | 3.20 |

\({\mathbf{A}}^{11}\) | 70,694 | 70,694 | 168,818 | 2.39 |

\({\mathbf{A}}^{12}\) | 70,694 | 2885 | 62,489 | 0.88/21.66 |

\({\mathbf{A}}^{22}\) | 2885 | 2885 | 4362 | 1.51 |

\({\mathbf{L}}\) | 73,579 | 73,579 | 191,026 | 2.60 |

\({\mathbf{L}}_1\) | 70,694 | 70,694 | 181,536 | 2.57 |

\({\widehat{{\mathbf{A}}}}^{-1}\) | 6833 | 6833 | 22,612 | 3.31 |

\(\widehat{{\mathbf{L}}}\) | 6833 | 6833 | 18,173 | 2.66 |

PCG iteration was implemented such that all operations were matrix by vector products. Thus, there was no need to build matrix \({\mathbf{Z}}_{imp}^{\prime}{\mathbf{Z}}_1^{\prime}{{\mathbf{R}}}^{-1}{\mathbf{Z}}_1{\mathbf{Z}}_{imp}\) but instead the required computations were performed stepwise from right to left in \({\mathbf{Z}}_{imp}^{\prime}({\mathbf{Z}}_1^{\prime}({{\mathbf{R}}}^{-1}({\mathbf{Z}}_1({\mathbf{Z}}_{imp}{\mathbf{v}}))))\) where \({\mathbf{v}}\) is a vector. This approach saves memory and allows fast computations [19].

Because of the very sparse pedigree relationship (\({{\mathbf{A}}}^{-1}\)) and factorization (\({\mathbf{L}}\)) matrices, the operations on each PCG iteration step that involve inverses of the sparse matrices, e.g. the imputation operations, can be calculated in linear time with respect to the number of individuals. The matrix multiplication of the marker matrix \({\mathbf{Z}}_m\) is also linear if the number of markers is assumed to be constant. Hence, the cubic complexity of inverting the genomic relationship matrix can be replaced by linear computational complexity of ssSNP-BLUP formulations. This holds if the iteration counts remain low.

### Accuracy

All the derived ssGBLUP MME formulations in this paper were linearly equivalent. Consequently, the final iterative solutions by all formulations were numerically equal and dictated only by the convergence tolerance. All formulations had “real” relative residual errors \(\widetilde{{\mathbf{e}}}_i\) in Eq. (66) of \(10^{-10}\) scale when using the residual convergence tolerance \(10^{-12}\) in pcg. Thus, all formulations indicated convergence to the same solutions, and no divergence of the iterative method was observed.

### Speed of iterative convergence

#### Genetic relationship: VanRaden 1

Number of iterations in PCG of linearly equivalent ssGBLUP MME using marker matrix \({\mathbf{Z}}_m^{vR1}\) (VanRaden 1), convergence tolerance \(10^{-12}\), heritability 0.5, and no preconditioning under different polygenic proportions* w*

MME | Size | Weight \({w}\) | |||||
---|---|---|---|---|---|---|---|

0.00 | 0.01 | 0.10 | 0.20 | 0.30 | 1.00 | ||

Orig. | 73,580 | – | 362 | 358 | 358 | 355 | 357 |

1 | 73,580 | – | 536 | 512 | 536 | 523 | 536 |

2 | 111,106 | 646 | 1203 | 1085 | 987 | 886 | 536 |

3 | 181,800 | 646 | 4023 | 3660 | 3307 | 3040 | 355 |

4 | 181,800 | 193 | 196 | 193 | 191 | 190 | 182 |

5 | 181,800 | 193 | 197 | 194 | 192 | 190 | 181 |

6 | 115,054 | 193 | 196 | 193 | 191 | 190 | 182 |

When \(w = 1\), all equivalent MME coincide with the traditional animal model with no genomic information and the numbers of iterations are smaller, even for the badly behaving third formulation (40). When \(w = 0\), the original ssGBLUP and the first equivalent MME (36), basic equivalent ssGBLUP, cannot be calculated because of the singularity in the inverse of the genomic relationship matrix \({{\mathbf{G}}_g}^{-1}\) in the test case. The numbers of iterations of the second (39) and third (40) MME are much lower in this case, which implies that the higher iteration counts are due to the RPG terms.

All fully orthogonalized equivalent ssSNP-BLUP MME 4 to 6 needed about half the number of iterations of the original ssGBLUP in all cases of \(w\). Note that the fifth formulation (60), i.e. orthogonal expanded ssSNP-BLUP, is an orthogonalized version of the badly behaving third formulation (40) while the fourth (59), i.e. orthogonal ssSNP-BLUP, and the sixth (65), i.e. reduced orthogonal ssSNP-BLUP, are orthogonalized versions of the second MME (39). The applied orthogonalization clearly reduced the number of iterations. Of the three orthogonalized ssSNP-BLUP formulations 4 to 6, the last, number 6 (65) is preferred because it has the smallest number of unknowns (115,054).

It should be noted that the convergence results apply only on the data used in the example. The equivalent MME formulation 2 (39), i.e. basic RPG ssSNP-BLUP, had 37,526 equations more than the MME formulation 1 (36), i.e. basic equivalent ssGBLUP. These extra equations were for the SNP solutions. Furthermore, equivalent MME 3 (40), i.e. expanded RPG ssSNP-BLUP, had again 70,694 new equations associated with RPG of non-genotyped animals. It remains to be tested whether MME 2 and 3 are more competitive when the number of genotyped animals is either close to or larger than the number of SNPs.

Number of iterations in PCG of linearly equivalent ssGBLUP MME using marker matrix \({\mathbf{Z}}_m^{vR1}\) (VanRaden 1), convergence tolerance \(10^{-12}\), heritability 0.5, and diagonal preconditioning under different polygenic proportions *w*

MME | Size | Weight \({w}\) | |||||
---|---|---|---|---|---|---|---|

0.00 | 0.01 | 0.10 | 0.20 | 0.30 | 1.00 | ||

Orig. | 73,580 | – | 93 | 62 | 59 | 57 | 57 |

1 | 73,580 | – | 152 | 153 | 158 | 158 | 171 |

2 | 111,106 | 847 | 1176 | 1099 | 1007 | 950 | 171 |

3 | 181,800 | 847 | 4058 | 3658 | 3195 | 2873 | 57 |

4 | 181,800 | 442 | 443 | 702 | 853 | 920 | 283 |

5 | 181,800 | 442 | 455 | 720 | 868 | 917 | 124 |

6 | 115,054 | 442 | 443 | 702 | 853 | 920 | 283 |

However, all fully orthogonalized ssSNP-BLUP formulations 4 to 6 converged much more slowly with diagonal preconditioning. For these formulations, the inverse of the diagonal of the MME matrix is not a good approximation of the inverse MME coefficient matrix. This could mean that the inverse MME coefficient matrix is not diagonally dominant or that the off-diagonal parts of the MME matrix contribute to the diagonal of the inverse MME matrix.

Number of iterations in PCG of linearly equivalent ssGBLUP MME using marker matrix \({\mathbf{Z}}_m^{vR1}\) (VanRaden 1), convergence tolerance \(10^{-12}\), heritability 0.1, and no preconditioning under different polygenic proportions *w*

MME | Size | Weight \({w}\) | |||||
---|---|---|---|---|---|---|---|

0.00 | 0.01 | 0.10 | 0.20 | 0.30 | 1.00 | ||

Orig. | 73,580 | – | 619 | 622 | 624 | 624 | 623 |

1 | 73,580 | – | 758 | 733 | 701 | 728 | 744 |

2 | 111,106 | 736 | 1005 | 970 | 907 | 836 | 748 |

3 | 181,800 | 736 | 3422 | 3315 | 3154 | 2934 | 623 |

4 | 181,800 | 74 | 74 | 72 | 72 | 71 | 71 |

5 | 181,800 | 74 | 73 | 72 | 70 | 71 | 70 |

6 | 115,054 | 74 | 74 | 72 | 72 | 71 | 71 |

Number of iterations in PCG of linearly equivalent ssGBLUP MME using marker matrix \({\mathbf{Z}}_m^{vR1}\) (VanRaden 1), convergence tolerance \(10^{-12}\), heritability 0.1, and diagonal preconditioning under different polygenic proportions *w*

MME | Size | Weight \({w}\) | |||||
---|---|---|---|---|---|---|---|

0.00 | 0.01 | 0.10 | 0.20 | 0.30 | 1.00 | ||

Orig. | 73,580 | – | 178 | 131 | 126 | 122 | 119 |

1 | 73,580 | – | 191 | 156 | 148 | 145 | 139 |

2 | 111,106 | 438 | 491 | 459 | 429 | 398 | 139 |

3 | 181,800 | 438 | 1732 | 1592 | 1416 | 1247 | 119 |

4 | 181,800 | 107 | 105 | 130 | 151 | 160 | 82 |

5 | 181,800 | 107 | 107 | 132 | 151 | 162 | 49 |

6 | 115,054 | 107 | 105 | 130 | 151 | 160 | 82 |

When the heritability was equal to 0.2, the diagonally preconditioned original ssGBLUP and non-preconditioned orthogonalized ssSNP-BLUP formulations all achieved relative convergence of \(10^{-12}\) in about the same 100 iterations (results not shown).

#### Genetic relationship: −1,0,1 encoding

Number of iterations in PCG of linearly equivalent MME using marker matrix \({\mathbf{Z}}_m^{-1,0,1}\) (−1,0,1 encoding), convergence tolerance \(10^{-12}\), heritability 0.5, and no preconditioning under different polygenic proportions *w*

MME | Size | Weight \({w}\) | |||||
---|---|---|---|---|---|---|---|

0.00 | 0.01 | 0.10 | 0.20 | 0.30 | 1.00 | ||

Orig. | 73,580 | – | 372 | 366 | 362 | 364 | 357 |

1 | 73,580 | – | 518 | 512 | 512 | 515 | 536 |

2 | 111,106 | 620 | 1141 | 1082 | 996 | 896 | 536 |

3 | 181,800 | 620 | 3794 | 3586 | 3390 | 3096 | 355 |

4 | 181,800 | 189 | 193 | 192 | 190 | 190 | 182 |

5 | 181,800 | 189 | 193 | 191 | 190 | 192 | 181 |

6 | 115,054 | 189 | 193 | 192 | 190 | 190 | 182 |

Number of iterations in PCG of linearly equivalent MME using marker matrix \({\mathbf{Z}}_m^{-1,0,1}\) (−1,0,1 encoding), convergence tolerance \(10^{-12}\), heritability 0.5, and diagonal preconditioning under different polygenic proportions *w*

MME | Size | Weight \({w}\) | |||||
---|---|---|---|---|---|---|---|

0.00 | 0.01 | 0.10 | 0.20 | 0.30 | 1.00 | ||

Orig. | 73,580 | – | 96 | 64 | 60 | 58 | 57 |

1 | 73,580 | – | 144 | 144 | 147 | 152 | 171 |

2 | 111,106 | 785 | 1059 | 1068 | 1010 | 953 | 171 |

3 | 181,800 | 785 | 3570 | 3411 | 3080 | 2826 | 57 |

4 | 181,800 | 418 | 413 | 636 | 762 | 817 | 283 |

5 | 181,800 | 418 | 423 | 673 | 805 | 857 | 124 |

6 | 115,054 | 418 | 413 | 636 | 762 | 817 | 283 |

### Compatibility of the genomic relationship matrix \({\mathbf{G}}_w\)

Convergence properties of the PCG method depend on the model used, data, and parameters. The model in our statistical analysis was very simple but we used different genomic relationship matrices. When genomic data is involved, several approaches are available to construct the genomic relationship matrix and the marker matrix \({\mathbf{Z}}_m\). Previously, Strandén and Christensen [23] had shown that differences in \({\mathbf{Z}}_m\) marker matrix in SNP-BLUP type models can give different mixing properties in Markov chain Monte Carlo computations. In ssGBLUP, the genomic and pedigree relationship matrices need to be constructed properly in order to avoid bias in the breeding values. We used different genomic relationship matrices by changing the \(w\) parameter without trying to maximize prediction ability in our data. Our choice for the family of genomic relationship matrices showed differences in convergence that may be partly due to, a potentially suboptimal, \({{\mathbf{H}}}^{-1}\) matrix.

In practice, the genomic relationship matrix in ssGBLUP or ssSNP-BLUP can be built such that it is “compatible” with the pedigree relationship matrix. There are several strategies for this e.g., [24, 25]. For each strategy an equivalent ssSNP-BLUP can be derived similarly as done for the family of relationship matrices in our study. For example, when the genomic relationship matrix has the form \(a {\mathbf{1}}{\mathbf{1}}^{\prime} + {\mathbf{G}}_g\) instead of \({\mathbf{G}}_w\), there will be one effect due to \(a {\mathbf{1}}{\mathbf{1}}^{\prime}\) instead of the residual polygenic effects. The equations for the marker effects remain the same.

Some studies have indicated that the \({({\mathbf{A}}_{22})}^{-1}\) matrix in ssGBLUP should be scaled by a factor less than 1, e.g., [26, 27]. This suggests that not only the genomic relationship matrix should be carefully constructed but also the pedigree relationship matrix should be adjusted. This has been properly formulated in the metafounders approach [28] where the pedigree relationship matrix is modified.

### Comparison to other approaches

The model by Fernando et al. [6] did not include RPG effects. The existence of RPG effects can be justified because not all the additive genetic variation can be explained by marker genotypes [2, 29]. Moreover Goddard et al. [30] suggested that with a finite number of markers, estimates of the genomic relationships are subject to error. The error can be due to sampling variation or inaccuracies in the analysis, and is inversely related to the number of markers in the analysis. In practice, the RPG has been included in most genomic evaluations because a moderate \(w\) is known to reduce the prediction bias in young selection candidates [27]. However, \(w\) might not have notable effects in single-step models [27].

The RPG effect could be included in the models of Fernando et al. [6] as a general random effect with \(Var({\mathbf{u}}) \sim N({\mathbf{0}}, w {\mathbf{A}})\). However, when the same observations of non-genotyped animals in [6] are modelled by RPG, the iterative approaches might face problems for separating the RPG and the “imputation residual” [6] (i.e. \({\mathbf{a}}_1\)) effects from each other. This was clearly visible in our equivalent MME 3 (40), i.e. expanded RPG ssSNP-BLUP, which showed very poor convergence.

When RPG was not included (\(w = 0\)) in equivalent MME 2 and 3, basic and expanded RPG ssSNP-BLUP, they converged much faster. However, all the other equivalent MME alternatives reached convergence with a smaller number of iterations.

Liu et al. [8] proposed a different approach for the use of marker effects in single-step models. Their model for the observations did not include the SNP effects, but instead the SNP effects were introduced into the MME as augmented effects correlated with the aggregated genomic breeding values, i.e. the sum of the genomic breeding value and the RPG effect. In this way, the marker matrix generated the dependencies among genotyped animals similar to genomic relationships in the \({\mathbf{H}}\) matrix. However, the MME by Liu et al. [8] was reported to have problematic convergence properties [8] (Liu Z: personal communication). This may have been due to poor condition number of MME coefficient matrix which resulted from large off-diagonal elements in blocks connecting marker effects and aggregated genomic values.

## Conclusions

A procedure was presented to derive linearly equivalent MME formulations for ssGBLUP. Six ssGBLUP based MME were derived, of which five were ssSNP-BLUP. In ssSNP-BLUP, inversion of the genomic relationship matrix is avoided. Three of the derived formulations were fully orthogonalized such that all random effects had diagonal covariance matrices.

All matrix operations and matrices, except the marker matrix \({\mathbf{Z}}_m\), were expressed using sparse matrices and sparse decompositions. During the iteration, the on-the-fly imputation from the genotyped to the non-genotyped animals was used on the genomic breeding values and residual polygenic effects without any need to explicitly predict and store genotypes of relatives of genotyped animals. All implementations used efficient matrix times vector operations where very sparse matrices were saved in memory. This enables efficient iteration on data-based implementations for large datasets and models.

All the derived MME gave exactly the same breeding value estimates at convergence with the small test data. Using the expected heritability of 0.5, the fully orthogonalized non-preconditioned ssSNP-BLUP formulations needed more iterations than the diagonally preconditioned original ssGBLUP. However, the number of iterations depended on the heritability value used. With a smaller heritability (0.1), the non-preconditioned ssSNP-BLUP formulations needed less iterations. The ssSNP-BLUP formulations did not benefit from diagonal preconditioner in PCG iteration.

The blending of the genomic relationship matrix \({\mathbf{G}}_g\) with the pedigree relationship matrix \({\mathbf{A}}_{22}\) was transformed into an additional residual polygenic effect in MME. The number of estimated additional random effects from all animals was reduced to include genotyped and their non-genotyped ancestors only. However, this did not improve convergence, although the number of unknowns to solve was smaller. When no residual polygenic effect was included, the ssSNP-BLUP models converged with a smaller number of iterations.

In conclusion, inversion of the large and dense genomic relationship matrix in the ssGBLUP can be avoided by using fully orthogonalized ssSNP-BLUP formulations. Although the new algorithms are more complicated than the original ssGBLUP, numerical efficiency is better when the number of genotyped individuals is large. The number of iterations until convergence by PCG was smaller in orthogonalized ssSNP-BLUP than in the original ssGBLUP when the heritability was low, but increased above that of the original ssGBLUP when heritability was higher.

The performance of the new MME should be further tested by analyzing data with more genotyped and non-genotyped animals.

## Declarations

### Authors' contributions

MT derived the formulae together with IS and EAM. MT did the main data analysis with IS supplying the data and performing additional comparisons. MT wrote the first drafts of the manuscript, and EAM and IS helped to revise and finalize it. All authors read and approved the final manuscript.

### Acknowledgements

This work was a part of the Genomics in Herds project which is a joint effort of Luke and Aarhus University. Breeding organizations Viking Genetics, Faba, and Breed4Food, the MAKERA foundation (Ministry of Agriculture and Forestry) and Valio Ltd were all participating funding of the work.

### Competing interests

The authors declare that they have no competing interests.

**Open Access**This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## Authors’ Affiliations

## References

- 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.View ArticlePubMedGoogle Scholar
- Christensen OF, Lund MS. Genomic prediction when some animals are not genotyped. Genet Sel Evol. 2010;42:2.View ArticlePubMedPubMed CentralGoogle Scholar
- Misztal I, Legarra A, Aguilar I. Using recursion to compute the inverse of the genomic relationship matrix. J Dairy Sci. 2014;97:3943–52.View ArticlePubMedGoogle Scholar
- 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 genotyped animals. J Dairy Sci. 2016;99:1968–74.View ArticlePubMedGoogle Scholar
- Masuda Y, Misztal I, VanRaden PM. Single-step GBLUP using APY inverse for protein yield in U.S. Holstein with a large number of genotyped animals. In: Proceedings of the 2016 ASAS-ADSA-CSAS-WSASAS joint annual meeting; Salt Lake City. 2016.Google Scholar
- Fernando RL, Dekkers JCM, Garrick DJ. A class of Bayesian methods to combine large numbers of genotyped and non-genotyped animals for whole-genome analyses. Genet Sel Evol. 2014;46:50.View ArticlePubMedPubMed CentralGoogle Scholar
- Legarra A, Ducrocq V. Computational strategies for national integration of phenotypic, genomic, and pedigree data in a single-step best linear unbiased prediction. J Dairy Sci. 2012;95:4629–45.View ArticlePubMedGoogle Scholar
- Liu Z, Goddard ME, Reinhardt F, Reents R. A single-step genomic model with direct estimation of marker effects. J Dairy Sci. 2014;97:5833–50.View ArticlePubMedGoogle Scholar
- Strandén I, Mäntysaari EA. Comparison of some equivalent equations to solve single-step GBLUP. In: Proceedings of the 10th World Congress on genetics applied to Livestock production. Vancouver; 2014. p. 22.Google Scholar
- Golub GH, Van Loan CF. Matrix computations. Baltimore: Johns Hopkins University Press; 1985.Google Scholar
- Henderson CR. Applications of linear models in animal breeding. Guelph: University of Guelph Press; 1984.Google Scholar
- Mrode R, Thompson R. An alternative algorithm for incorporating the relationships between animals in estimating variance components. J Anim Breed Genet. 1989;106:89–95.View ArticleGoogle Scholar
- Quaas RL. Linear prediction. In: BLUP School Handbook. Armidale: A.G.B.U., University of New England; 1984. https://www.cabdirect.org/cabdirect/abstract/19850186676.
- Henderson CR. Sire evaluation and genetic trends. In: Proceedings of the animal breeding and genetics symposium in honor of Dr. Jay L. Lush: 29 July 1972. Blacksburg; 1973. p. 10–41.Google Scholar
- VanRaden PM. Efficient methods to compute genomic predictions. J Dairy Sci. 2008;91:4414–23.View ArticlePubMedGoogle Scholar
- 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.View ArticlePubMedGoogle Scholar
- Legarra A, Aguilar I, Misztal I. A relationship matrix including full pedigree and genomic information. J Dairy Sci. 2009;92:4656–63.View ArticlePubMedGoogle Scholar
- Henderson CR. A Simple method for computing the inverse of a numerator relationship matrix used in prediction of breeding values. Biometrics. 1976;32:69–83.View ArticleGoogle Scholar
- Strandén I, Lidauer M. Solving large mixed linear models using preconditioned conjugate gradient iteration. J Dairy Sci. 1999;82:2779–87.View ArticlePubMedGoogle Scholar
- Colleau JJ. An indirect approach to the extensive calculation of relationship coefficients. Genet Sel Evol. 2002;34:409–21.View ArticlePubMedPubMed CentralGoogle Scholar
- Eaton JW, et al. GNU Octave. www.gnu.org/software/octave. Accessed 15 Mar 2017.
- Chen Y, Davis TA, Hager W, Rajamanickam S. Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate. ACM Trans Math Softw. 2008;35:22.View ArticleGoogle Scholar
- Strandén I, Christensen OF. Allele coding in genomic evaluation. Genet Sel Evol. 2011;43:25.View ArticlePubMedPubMed CentralGoogle Scholar
- Vitezica ZG, Aguilar I, Misztal I, Legarra A. Bias in genomic predictions for populations under selection. Genet Res (Camb). 2011;93:357–66.View ArticlePubMedGoogle Scholar
- Gao H, Christensen OF, Madsen P, Nielsen US, Zhang Y, Lund MS, et al. Comparison on genomic predictions using three GBLUP methods and two single-step blending methods in the Nordic Holstein population. Genet Sel Evol. 2012;44:8.View ArticlePubMedPubMed CentralGoogle Scholar
- Tsuruta S, Misztal I, Aguilar I, Lawlor T. Multiple-trait genomic evaluation of linear type traits using genomic and phenotypic data in US Holsteins. J Dairy Sci. 2011;94:4198–204.View ArticlePubMedGoogle Scholar
- Koivula M, Strandén I, Pösö J, Aamand GP, Mäntysaari EA. Single-step genomic evaluation using multitrait random regression model and test-day data. J Dairy Sci. 2015;98:2775–84.View ArticlePubMedGoogle Scholar
- Legarra A, Christensen OF, Vitezica ZG, Aguilar I, Misztal I. Ancestral relationships using metafounders: finite ancestral populations and across population relationships. Genetics. 2015;200:455–68.View ArticlePubMedPubMed CentralGoogle Scholar
- Erbe M, Gredler B, Seefried FR, Bapst B, Simianer H. A function accounting for training set size and marker density to model the average accuracy of genomic prediction. PLoS One. 2013;8:e81046.View ArticlePubMedPubMed CentralGoogle Scholar
- Goddard ME, Hayes BJ, Meuwissen THE. Using the genomic relationship matrix to predict the accuracy of genomic selection. J Anim Breed Genet. 2011;128:409–21.View ArticlePubMedGoogle Scholar