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

Taskinen et al. Genet Sel Evol (2017) 49:36 DOI 10.1186/s12711-017-0310-9 RESEARCH ARTICLE Single-step SNP-BLUP withon-the-fly imputed genotypes andresidual polygenic effects Matti Taskinen* , Esa A. Mäntysaari and Ismo Strandén 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 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 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. © The Author(s) 2017. 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. 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 G−1 g . In addition, matrix Gg 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 G−1 APY approximation to the G−1 g 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 Open Access G enetics S election Evolution *Correspondence: [email protected] Natural Resources Institute Finland (Luke), Myllytie 1, Jokioinen, Finland Page 2 of 15 Taskinen et al. Genet Sel Evol (2017) 49:36 in APY will give different evaluations, which may affect selection decisions. An alternative formulation called hereinafter singlestep single nucleotide polymorphism BLUP (ssSNPBLUP) [6] overcomes some of the major computational challenges in ssGBLUP. In particular, there is no need to construct or invert the genomic relationship matrix Gg . 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 ssSNPBLUP. 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 Two mixed linear effects models, i.e. with the same observations y but different fixed ( b and  b ) and random effects ( u and  u ), residual errors ( e and  e ), model matrices ( X , Z ,  X , and  Z ) and variance structures ( G , R ,  G , and  R ), are said to be 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(1) can be solved from separate equations for fixed and random effects [14] as: where matrix V needs to be invertible. The size of matrix V is the number of observations, which can be very (1) y=Xb +Zu +e,Var(u)=G,Var(e)=R (2) y=  X  b+  Z u+  e,Var( u)=  G,Var( e)=  R (3)  Xb =  X  b, V=ZGZ′+R=  Z  G  Z′+  R=  V . (4)  b =(X′V−1X) −1 X′V−1y  u=GZ′V− 1 (y−X  b) , Page 3 of 15 Taskinen et al. Genet Sel Evol (2017) 49:36 large, and, therefore, solving the mixed model using this method is seldom feasible in practice. Mixed model equations In practice, Eq.(4) can be solved by Henderson’s MME [14] as follows: where the variance matrix of the random effects G and the residual variance matrix R need to be invertible. MME(5) usually lead to more sparse matrix systems than Eq.(4). Equivalent models bysplitting the variance matrix Suppose the variance matrix G can be expressed as a matrix product as follows: where M is rectangular and  G is an invertible square matrix. Here, the matrix  G could be, for example, an identity matrix and could have different dimensions than the G matrix. Matrices G and Z are always together in Eq.(4). This allows us to re-parametrize the model as: and thereafter: where the equivalent model has the same residual variance matrix (  R=R ), the new model matrix  Z=ZM , and the random effects  u are defined as: Now, according to Eq. (3), the quantities with a tilde together with the same observation vector y , the original fixed effects (  b=  b ), and design matrix (  X=X ) form a linearly equivalent model(2). Linearly equivalent MME Original fixed effects  b and the new random effects  u of the linearly equivalent model can be solved, similarly as in Eq.(4), from: or from the corresponding MME(5): (5)  X ′ R −1 XX ′ R −1 Z Z ′ R − 1XZ ′ R − 1Z + G − 1  b  u = X ′ R −1 y Z ′ R − 1y , (6) G=M  GM′, (7) V =ZGZ ′ +R =ZM  GM′Z′+R =  Z  G  Z′+  R=  V, (8)  u =M  GM′Z′V−1(y−X  b) =M  G  Z′ V− 1 (y−X  b)=M u, (9)  u=  G  Z ′ V − 1(y − X  b) . (10)  b =(X′  V−1X) −1 X′  V−1y  u=  G  Z′ V− 1 (y−X  b), Because of the equivalence  u=M u in Eq.(8), the original random effects  u can be obtained from the solution of the linearly equivalent MME (11) by pre-multiplying the new random effects  u with matrix M , i.e. where identity matrix Ib has dimension of  b . Note that the inverse of the MME matrix in Eq.(12) is usually not evaluated explicitly, but rather, the corresponding linear matrix equation is solved using either direct or iterative solution methods. The original MME (5) can, thus, be solved from a modified linearly equivalent MME (11) where the number of random effects in  u could be smaller or larger than in  u , the variance structure  G could be easier to obtain or invert than the original G , or the new matrix system could be otherwise numerically more efficient. Presentation forGBLUP andSNP‑BLUP As an example, consider single-trait GBLUP MME [15]. The variance matrix of random effects u is based on a genomic relationship matrix ( Gg ), which describes the genomic relationships between individuals, i.e. Var (u) = σ 2 u G g . The genomic relationship matrix is usually fully dense and increases in order as the number of genotyped animals increases. The inverted genomic relationship matrix Gg −1 is needed in the solution of the GBLUP MME: where a single trait case is assumed for simplicity, R= σ 2 eI , and = σ 2 e σ 2 u . There are many ways to construct Gg . Assume the genomic relationship matrix Gg can be expressed, in simplified form, using (centered and scaled) marker matrix Zm [15] so that: where and Im is an identity matrix of size equal to the number of markers. Now, a linearly equivalent MME(11), alternative to GBLUP MME(13), can be derived and original effects solved similarly as in(12) so that: (11)  X′R−1XX ′R−1  Z  Z′R− 1 X Z′R− 1  Z+ G− 1  b  u = X′R−1y  Z ′ R − 1y . (12)  b  u  =  Ib0 0M  X′R−1XX ′R−1  X  Z′R−1X  Z′R−1  Z+  G−1 −1 X′R−1y  Z′R−1y , (13)  b  u = X′XX ′Z Z′XZ ′Z + G−1 g−1 X′y Z′y , (14) Gg= Z m Z ′ m= M  GM ′, (15) M=Zmand  G=Im, Page 4 of 15 Taskinen et al. Genet Sel Evol (2017) 49:36 where  Z=ZZm and  is the same variance ratio as in Eq. (13). This equivalent MME system, known as the SNP-BLUP [16], has markers as random effects instead of individuals. Random effects are also “orthogonalized”, and, thus, inversion of the dense genomic relationship matrix Gg is avoided. Note that, if the marker effects have unequal variances, the relationship matrix can be build as: where and the matrix B is a diagonal covariance matrix describing the variances of different marker effects. Now the solution of SNP-BLUP becomes: 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 H [1, 2]. The relationships in H are described by the pedigree-based relationship matrix A , and the genomic relationship matrix among genotyped animals by Gg . Relationships among non-genotyped individuals are constructed from the pedigree but modified according to the relationships among genotyped animals. Assuming that the non-genotyped individuals are denoted with suband super-scripts 1 and the genotyped individuals by suband super-scripts 2, the pedigree relationship matrix and its inverse are the following: Assume, as before, that the genomic relationship matrix has the form Gg= Z m Z ′ m . If the genotyped population contains identical twins, i.e. clones, or if there are more individuals than markers, the genomic relationship matrix Gg becomes singular and the usual MME(5) cannot be constructed. Hence, Gg is commonly adjusted by regressing it towards the pedigree relationship matrix A22 with:  b  u = Ib0 0Zm  X′XX ′  Z  Z′X Z′ Z+  Im−1 X′y  Z ′ y , (16) Gg= Z m BZ ′ m= M  GM ′, (17) M=Zmand  G=B,  b  u = Ib0 0Zm  X′XX ′  Z  Z′X Z′ Z+  B− 1 −1 X′y  Z ′ y . (18) A= A11 A12 A21 A22  and A−1 = A 11 A 12 A21 A22 . (19) Gw=wA22 +(1−w)Gg, where w is a scalar weight between 0 and 1 and can be interpreted as the relative weight on the polygenic effect [2, 15]. Single‑step relationship matrix The inverse of the ssGBLUP variance matrix H is [1, 2]: Using block LDL decomposition, it is equal to: where I1 and I2 are identity matrices of size equal to the number of non-genotyped and genotyped animals, respectively. For imputation of genotypes: is a regression prediction or an 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 H has a similar decomposition: By using block matrix inversion identities (23) and to the “imputed” A22 matrix of Gw (19) the variance matrix H (24) can be alternatively expressed as: where Gimp is an imputed genomic relationship matrix as follows: Note that this operates on Gg instead of Gw . Also, note that Gimp has a size equal to the number of all animals. Genotyped animals have observed marker data but nongenotyped animals have imputed marker data. (20) H −1 = A−1 +00 0G w− 1 − (A 22 ) − 1  (21) = A 11 A 12 A 21 A 22 +Gw− 1 −(A22)− 1 . (22) H −1 = I10 − A′ imp I2  A 11 0 0G w − 1  I1−Aimp 0I 2 , (23) Aimp = A12(A22) − 1 =− (A11) −1 A 12 (24) H= I1Aimp 0I 2  (A11) −1 0 0G w I10 A′ imp I2 . (25) (A 11 ) −1 =A11 −A12(A22)− 1 A21 (26)  Aimp I2  A22  A′ imp I2 = A − (A11) −1 0 00 , (27) H= (1 − w)  (A11) −1 0 00 + wA + (1 − w)Gimp , (28) G imp = Aimp I2  Gg  A′ imp I2 . Page 5 of 15 Taskinen et al. Genet Sel Evol (2017) 49:36 Equivalent single‑step MME In the following, the single-step relationship matrix H will be expressed as six different decompositions equivalent to Eq.(6): for different matrices Mi and  Gi . All of these lead to a linearly equivalent ssGBLUP MME system of new sets of random effects with, potentially, different numerical properties. The linearly equivalent MME of these modified sets of random effects are similar to Eq.(11) and the original effects can be solved similarly as in Eq.(12): where  Zi=ZMi and a single trait case is assumed. Note that in these equivalent MME, square matrix  Gi represents the covariance structure for the reparametrized random effects, Mi hasa row for each original random effect in u to change the model matrix Z , and  Zi is the redefined model matrix. In order to derive linearly equivalent MME for multiple trait cases: where ⊗ is the Kronecker product, the genetic (co)variance matrix G0 of size number of traits is assumed to have a decomposition: For example, a simple case would have M0 as identity matrix and  G0 as G0 . Now the variance matrix G can be expressed using the decompositions of the single-step relationship matrix(29) similarly as in Eq.(6): where Effects can then be solved from linearly equivalent multiple trait MME (12): where  Zi=Z(M0⊗Mi) . (29) H=Mi  GiM′ i ,i = 1, ..., 6,  b  u  =  Ib0 0M i  X′XX ′  Zi  Z′ iX Z′ i Zi+   G− 1 i−1 X′y  Z ′ i y , (30) G=G0⊗H, (31) G0=M0  G0M′ 0. (32) G=G0⊗H=(M0  G0M′ 0)⊗(Mi  GiM′ i) (33) =(M0⊗Mi)(  G0⊗  Gi)(M0⊗Mi)′ (34) =M  GM′, (35) M=M0⊗Miand  G=  G0⊗  Gi.   b  u  =  Ib0 0M 0⊗Mi  ×  X′R−1XX ′R−1  Zi  Z′ iR−1X  Z′ iR−1  Zi+  G−1 0⊗  G−1 i  −1  X′R−1y  Z′ iR−1y , Basic equivalent ssGBLUP MME The LDL decomposition(24) can be used directly to build the first linearly equivalent ssGBLUP MME of this paper using H=M1  G1M′ 1 where: From the modified relationship matrix  G1 (36), it can be seen that this basic equivalent ssGBLUP MME has random effects for non-genotyped animals with variances (A 11 ) −1 and for genotyped animals with variances of the adjusted genomic relationship matrix Gw . The number of effects in the system is, thus, the same as in the original ssGBLUP. Basic RPG ssSNP‑BLUP MME Other linearly equivalent ssGBLUP MME of the form(29) can be derived as well. Note that the adjusted genomic relationship matrix Gw in Eq.(19) can be expressed by matrix products as follows: where Gg= Z m Z ′ m . The second linearly equivalent ssGBLUP MME can be built by substituting Gw in Eq. (38) to Eq. (36): In this form, we avoid the inverse of Gg in the MME (11). The coefficients w and (1−w) were also split using square roots to matrix M2 so that the new variance matrix  G2 can be inverted even when w is 0 or 1. The first group of new random effects in Eq.(39), for the non-genotyped animals, is the same as in the first, basic equivalent ssGBLUP in Eq.(36). However, the genotyped animals now have random effects related through the variance matrix A22 . Effects in this second effect group can be seen as 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 SNPBLUP(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. (36) M 1 = I1Aimp 0I 2  and  G1 = (A11) −1 0 0G w. (37) G w = I2I2  wA22 0 0(1 − w)Gg  I2 I2  (38) = I2Zm  wA22 0 0(1 − w)Im  I2 Z′ m, (39) M 2= � I1 √ wAimp √ 1−wAimpZm 0√wI2√1−wZm � � G2=  (A11)−100 0A 22 0 0 0I m   . Page 6 of 15 Taskinen et al. Genet Sel Evol (2017) 49:36 Expanded RPG ssSNP‑BLUP MME A third linearly equivalent MME of form(29) can be derived from the alternative expression of matrix H in Eq. (27) and by splitting Gg= Z m Z ′ m in Eq.(28) as: Here, matrices E1 and E2 are rectangular sparse incidence matrices that select the subsets of non-genotyped and genotyped animals, respectively, from the A matrix. Both E1 and E2 have the same number of columns, i.e. number of all animals. Matrix E1 has a row for each nongenotyped and matrix E2 for each genotyped animal corresponding to animal’s column in matrices Z1 and Z2 of Each row of both E1 and E2 has only one non-zero element, a value one at the column corresponding to that animal’s location among all of the animals. Hence, when rows and columns of the matrices of all animals are in the same order as in matrix A in Eq.(18), matrix  E1 E2  is an identity matrix of the size of all animals. 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 √1−w . Also, the second group, related through the pedigree relationship matrix 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→1 . At the limit ( w=1 ), the expanded RPG ssSNP-BLUP MME 3 (40) has clearly the recognizable covariance structure of A . The basic equivalent ssGBLUP MME1(36) and the basic RPG ssSNPBLUP MME 2(39) are models where the genotyped animals act as base animals and the non-genotyped animals are regressed on them. For the other direction of w→0 , the basic equivalent ssGBLUP MME 1(36) converges to an alternative presentation of the standard ssGBLUP. However, it divides (40) M 3= �√ 1−wI1 √ wE1 √ 1−wAimpZm 0√wE2√1−wZm � � G3=  (A11)−100 0 A0 0 0I m   . (41) Z=[Z1Z2]. the breeding values of non-genotyped animals into regressions on genotyped animals and into non-imputed breeding values that are not conditional on them. Similarly, at the limit w=0 the basic and the expanded RPG ssSNP-BLUP MME, 2(39) and 3(40), coincide with the simple ssSNP-BLUP without residual polygenic effects. For example, in the single trait case, the MME coefficient matrix of the basic RPG ssSNP-BLUP MME2 (39) is where W=(Z 1 Aimp +Z 2 )Zm . 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 Zm 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 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 The inverse of the pedigree relationship matrix A−1 can be expressed efficiently [18] as: where animals are sorted in (reversed) age order from the youngest to the oldest using sparse permutation matrix Q′ , so that matrix L= (I − 1 2 P) ′ D 1 2 becomes a lower triangular matrix in Q′A−1Q=LL′ . The diagonal matrix D has values 4/(4−k−Fs) where k is the number of known parents and Fs is the sum of parent inbreeding coefficients. In the “parental matrix” P on row i , there are 1s in columns corresponding to parents of animal i . The parental matrix can be interpreted, together with identity (42)     X ′ XX ′ Z1 0X ′ W Z′ 1XZ ′ 1Z1+A11 0Z ′ 1W 00A22−10 W′XW ′Z10W ′W + Im     , (43) A −1 = Q  I − 1 2P ′ D  I − 1 2P  Q′ = QLL′Q′ , Page 7 of 15 Taskinen et al. Genet Sel Evol (2017) 49:36 matrix I , as a very sparse lower triangular “Cholesky” matrix( L )[18]. The pedigree relationship matrix A can be expressed as the inverse of its inverse, and so the submatrices of the pedigree relationship matrix and its inverse(18) can be obtained by selecting the appropriate rows and columns as follows: where i,j=1, 2 . Matrix-vector products Aijx and Aijx can be efficiently computed using these decompositions [19, 20]. The submatrices Aij are very sparse, so they could alternatively be expressed as separate sparse matrices. The inverse of A−1 or particular parts of it (e.g. (A 11 ) −1 ) are, however, in general non-sparse and, thus, these inverse matrix terms should never be computed explicitly. The sparse submatrix A11 can be expressed using sparsity preserving Cholesky factorization so that: where matrix L1 is sparse lower triangular and Q1 sparse permutation matrix. Note that matrix L1 has to be computed explicitly, as opposed to matrix L in Eqs. (43)to(46). Matrix L1 has some more fill-ins compared to A11 but is still very sparse and efficient in use. In [4], it was demonstrated that the computations remain affordable even when the dataset size grows. Computations involving matrix inversions of A−1 or parts of it can be transformed into solutions of sparse matrix equation systems [20]: where v and v1 are vectors of appropriate sizes into which the inverse matrix operations are performed. Here the backslash ( \ ) is an operator indicating forward or back‑ ward substitutions and emphasizes the importance of avoiding inverting matrices. In other words, L\y is the solution x of equation Lx =y or can be expressed as solving ( L , y ). Note that the matrix products are carefully (44) A=(A− 1 )−1 =Q(L′)− 1 L− 1 Q′, (45) A ij = EiQ(L′) −1 L−1Q′E ′ j (46) Aij = EiQLL ′ Q ′ E ′ j, (47) A11 =Q1L1L′ 1Q′ 1, (48) Av =Q(L′) −1 L−1Q′v = Q  L ′\ L \ (Q ′ v)  (49) ( A11) −1 v1 = Q1  L ′ 1\ L1 \ (Q ′ 1 v1) , nested with parenthesis so that only matrix-vector operations are performed. Furthermore, following [4, 9] the inverse of the matrix A22 , needed in the inversion of the second modified relationship matrix  G2 in Eq.(39), can be expressed efficiently using block matrix inversion identity similar to Eq.(25) as: where all terms can be computed using Eqs. (45)to(49). On‑the‑fly imputation operation The derived equivalent formulations in Eqs. (36), (39), and (40) contain imputation operator Aimp =− (A11) −1 A 12 (23) in matrices Mi that are needed when operating with the modified model matrix  Zi=ZMi and its transpose  Z′ i= M ′ i Z ′ , 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 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. Aimpv2 ) or to a vector of the marker effects vm through the marker matrix (i.e. AimpZmvm ). Thus, in the transpose  Z′ i side, the imputation term operates on a vector of size equal to the number of non-genotyped animals (i.e. A′ imp v 1 ). In all cases, the size of the vector term that operates on (A 11 ) −1 equals the number of non-genotyped animals, i.e. size of v1 . For example, the imputed genomic marker data term, −(A 11 ) −1A 12 Zmvm , that expands the genomic information from genotyped to non-genotyped animals, can be calculated using: and Eq.(49) so that: Vector v1 is calculated from Eq.(46), or without constructing any of the matrices by using rules for A−1 by pedigree information [18], or, alternatively, as a sparse matrix-vector product of separate sparse matrix A12 . 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 (50) (A22)− 1 =A 22 −A 21 (A 11 ) −1A 12 , (51)  v2=Zmvm (52) v1=A12 v2 (53) − (A11) −1 A12Z m v m=− Q1  L ′ 1\ L1 \ (Q ′ 1 v1) . Page 8 of 15 Taskinen et al. Genet Sel Evol (2017) 49:36 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 Zm . 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  Zi and another time for the transpose  Z′ 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 (  Z′ i R − 1 y ) in Eq.(11) and once at the end of the iteration in order to retrieve the original random effects (  u=Mi  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 M as in Eqs.(6) and(29). The term (A 11 ) −1 in Eqs.(39) and (40) can be orthogonalized by using the sparse Cholesky factorization in Eq.(47) as follows: where Note that the permutation operator Q1 can be performed outside the inverse operator. However, the term A22 in Eq.(39) seems to be much more difficult to decompose. Still, it can be expressed using the sparse decomposition of the full matrix A in Eq.(44) as: (54) (A 11 )−1 =M11  G11M′ 11, (55) M11 =Q1 ( L′ 1 ) −1 and  G11 =I1. (56) A22 =E2AE′ 2=M22  G22M′ 22, where and The rectangular matrix  A 1 2 2 has dimensions number of genotyped individuals times total number of individuals, hence the trade-off here is that Eq.(56) will expand the second random effect group of Eq.(39) from genotyped to all individuals (size of I ). Using Eqs.(54), (56), and (58), the linearly equivalent MME (39) and (40) can be “orthogonalized” into fourth: and fifth equivalent MME. Both of these linearly equivalent(29) ssSNP-BLUPs share the same orthogonal variance structure: and, thus, both also have the same number of new random effects: random effects for the genotyped individuals, two sets of random effects for the non-genotyped individuals, and random effects for the markers. The difference in equivalent MME 4 (59) and 5 (60) is on how they divide the RPG on non-genotyped animals. The fourth equivalent MME (59) can be called orthog‑ onal 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 A22 , as the covariance structure for the genotyped animals, was reparametrized in Eq.(56) using the full pedigree relationship matrix A . This reparametrization increases the number of corresponding new random effects from genotyped animals to all animals in the pedigree. However, computations involving A22 require only the genotyped individuals and their ancestors. Thus, to reduce the number of extra new effects, A22 can be expressed using a smaller pedigree and relationship (57) M 22 =  A 1 2 2 and  G22 = I , (58)  A 1 2 k= E k Q(L′)−1,k =1, 2. (59) M 4=   M11 √wAimp � A 1 2 2 √ 1−wAimpZm 0√w � A 1 2 2 √1 − wZm  , (60) M 5=  √ 1−wM11 √w � A 1 2 1 √ 1−wAimpZm 0√w � A 1 2 2 √1 − wZm  , (61) � G 4= � G5=   I1 00 0I 0 0 0I m  , Page 9 of 15 Taskinen et al. Genet Sel Evol (2017) 49:36 matrix  A containing the genotyped animals and their ancestors [4]. Let the inverse of the pedigree relationship matrix(43) of this smaller pedigree be: where  Q and  L are as before in Eq.(44) but involve genotyped animals and their ancestors only. Matrix A22 in Eq.(56) can then be represented using the smaller pedigree as: where  E2 selects the genotyped individuals from the smaller pedigree, and the size of identity matrix Iganc is the number of genotyped animals and their ancestors. The sixth linearly equivalent, reduced orthogonal ssSNP‑BLUP MME, can be derived from Eqs.(59) and (61) as: Here only the non-genotyped ancestors of the genotyped have two sets of random effects, all the other nongenotyped and all genotyped animals have single sets of effects, in addition to the marker effects. 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 (62)  A− 1 =  Q  L  L′  Q′, (63) A22 =  E 2  A  E ′ 2=  M 22  G 22  M ′ 22, (64)  M 22 =  E2  Q(  L ′ ) −1 and  G22 = I ganc, (65) M 6= � M11 √ wAimp � M22 √ 1−wAimpZm 0√w� M22 √1−wZm � � G6=  I100 0I ganc 0 0 0I m   . 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 Gg without A22 adjustment, i.e. in the case w=0 . Comparison statistics The original ssGBLUP with inverse variance matrix 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 Gg −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. Accuracy of the formulations was tested by solving the MME with the PCG method using Octave’s PCG routine (pcg) with the diagonal of the MME coefficient matrix as the preconditioner, or without preconditioning. Convergence tolerance in pcg was relative residual norm. The tolerance was chosen to be small ( 10−12 ) so that all solved effects, without doubt, converged. Accuracies, or rather the differences from the “exact solution”, were calculated as relative residual errors (  ei ) between iteratively obtained MME solutions ( si ) and the direct solution ( sdirect ) of the original ssGBLUP: where subscript i=1, ...,6 is the formulation number. 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 (66)  ei = � s direct − s i� �sdirect �,