Implications

Tools used in genomic selection range from those based on estimation of SNP effects to animal models using a genomic relationship matrix (GRM). When many animals are ungenotyped, modeling options include multi-step or single-step methods. Computations involving mixed models with SNP effects or the inverse of GRM have good computing properties, whereas those including their combination or plain GRM do not. Recent discovery of efficient inverse of GRM removed computing limits from genomic analyses involving GRM. Same discovery enables GREML for very large data sets.

Definition of terms and basic models

Let A denote a NRM based on pedigrees. Let Z be a matrix of gene content with z _ij being the number of occurrences of the minor allele in SNP j of animal i. The values of ‘raw’ z _ij are 0, 1 and 2, whereas that of centered z _ij is −2p _j , 1−p _j and 2−p _j , where p _j is gene frequency of SNP j. The simplest SNP BLUP model is

$$y_{i} \,{\equals}\,\mu {\plus}\mathop{\sum}\limits_j^{} {z_{{ij}} a_{j} {\plus}e_{i} } $$

where y _i is a (pseudo)phenotype of animal i, $${\bf a}\,\sim\,N\left( {0,{\bf I}\sigma _{a}^{2} } \right)$$ the vector of SNP effects, $${\bf e}\,\sim\,N\left( {0,{\bf I}\sigma _{e}^{2} } \right)$$ the vector of residuals, and $$\sigma _{a}^{2} $$ and $$\sigma _{e}^{2} $$ are SNP and residual variances, respectively. GBLUP can be defined as

$$y_{i} \,{\equals}\,\mu \,{\plus}\,u_{i} \,{\plus}\,e_{i} $$

where $${\bf u}\,\sim\,N\left( {0,{\bf G}\sigma _{u}^{2} } \right)$$ is a vector of animal effects, G the GRM and $$\sigma _{u}^{2} $$ the additive variance. G can be derived by a transformation from SNP BLUP:

$${\bf{G}} \,{\equals}\,{\bf{ZZ'}} {{\sigma _{a}^{{\rm 2}} } \over {\sigma _{u}^{{\rm 2}} }}$$

possibly with an alternative scaling factor (VanRaden, Reference VanRaden2008):

$${\bf{G}} \,{\equals}\,{{{\bf{ZZ'}} } \over {\mathop{\sum}\limits_j^{} {\left[ {p_{j} \left( {{\rm 1}{\minus}p_{j} } \right)} \right]} }}$$

If SNP BLUP involves heterogeneous variance for SNP: $${\bf a}\,\sim\,N(0,{\bf D}\sigma _{a}^{2} )$$ , where D is a diagonal matrix of weights, the equivalent ‘weighted’ G is

$${\bf{G}} \,{\equals}\,{{{\bf{ZDZ'}} } \over {\mathop{\sum}\limits_j^{} {\left[ {p_{j} \left( {{\rm 1}{\minus}p_{j} } \right)} \right]} }}$$

SNP BLUP and GBLUP are equivalent models (Strandén and Christensen, Reference Strandén and Lidauer2011).

When only a fraction of animals is genotyped, the numerator and genomic relationships can be combined (Legarra et al., Reference Legarra, Christensen, Aguilar and Misztal2009; Christensen and Lund, Reference Christensen and Lund2010):

$${\bf H}\,{\equals}\;{\bf A}{\plus}\left[ {\matrix{ {{\bf A}_{{12}} {\bf A}_{{22}}^{{{\minus}1}} } & {\bf 0} \cr {\bf 0} & {\bf I} \cr } } \right]\left[ {\matrix{ {\bf I} \cr {\bf I} \cr } } \right]({\bf G}{\minus}{\bf A}_{{22}} )\left[ {\matrix{ {\bf I} & {\bf I} \cr } } \right]\left[ {\matrix{ {{\bf A}_{{22}}^{{{\minus}1}} {\bf A}_{{21}} } & {\bf 0} \cr {\bf 0} & {\bf I} \cr } } \right]$$

where H is a combined matrix and indices in A refer to ungenotyped (1) and genotyped (2) animals. The inverse of matrix H is (Aguilar et al., Reference Aguilar, Legarra and Misztal2010; Christensen and Lund, Reference Christensen and Lund2010)

$${\bf H}^{{{\minus}1}} \,{\equals}\;{\bf A}^{{{\minus}1}} {\plus}\left[ {\matrix{ {\bf 0} & {\bf 0} \cr {\bf 0} & {{\bf G}^{{{\minus}1}} {\minus}{\bf A}_{{22}}^{{{\minus}1}} } \cr } } \right]$$

BLUP with matrix H is called ssGBLUP.

Computing operations in animals breeding

Common operations in animal breeding include solving of mixed model equations and variance component estimation either via REML or Bayesian methods with Gibbs sampling. Such operations in the animal breeding context are described in Mrode (Reference Misztal and Perez-Enciso2014) and Misztal (Reference Misztal, Legarra and Aguilar2014), and are typically based on Henderson’s mixed model equations (MME). Let this linear system of equations be

$${\bf{B}} \,{\bf{x}} \,{\bf{{\equals}}} \;{\bf{y}} $$

where B is left-hand-side (LHS), y the right-hand side, and x the vector of solutions. Solving mixed model equations can be accomplished by finite and iterative methods. In finite methods, exact solutions (except for numerical errors) are obtained in a finite number of steps. Popular methods for general matrices are based on the LU decomposition (B=LU, where L and U are lower and upper triangular, respectively) and for symmetric matrices the Cholesky decomposition (B=LL', where L is lower triangular). As often mixed model equations are not full rank, finite methods used for mixed models need to compute generalized solutions rather than cause numerical exceptions. For example, the Cholesky decomposition can be modified to

$${\bf{B\,{\equals}\,LDL}} \prime$$

where D is a diagonal matrix with zero elements for redundant equations and non-zeros otherwise. Although finite methods usually compute solutions with high precision, their cost is high, usually close to cubic. Finite methods can be used to calculate inverses, which are used in REML or for prediction error variance (PEV) calculations. In general, the cost of the inverse is a few times that of a solution by decomposition.

When matrices are sparse, with nearly all elements of LHS equal to 0, a ‘sparse’ Cholesky decomposition exists that avoids operations with zeros, leading to about quadratic rather than cubic costs (George and Liu, Reference George and Liu1981). Also, an algorithm exists for computing a ‘sparse’ inverse, where elements computed are those corresponding to non-zeroes in the original matrix (Takahashi et al., Reference Takahashi, Fagan and Chen1973). Such an inverse is sufficient to calculate traces in REML (Misztal and Perez-Enciso, Reference Misztal1993) or PEV. However, equations involving SNP effects are usually dense. When the system of equations is sparse but contains dense blocks, a sparse matrix package can be modified for greatly improved performance in REML and other applications (Masuda et al., Reference Masuda, Misztal, Tsuruta, Legarra, Aguilar, Lourenco, Fragomeni and Lawlor2015). This is quite new as most matrix packages are efficient for either dense or sparse matrices, but not for mixtures of both.

Iterative methods are such that progressively more accurate solutions are obtained every next round. Two such methods are popular in animal breeding. The first one is Gauss–Seidel iteration (GS), which is also a backbone of Gibbs sampling. GS is very simple, but requires access to equations by rows. The second one is preconditioned conjugate gradient (PCG), a method difficult to understand, but easy to implement and usually with superior convergence rate to GS (Tsuruta et al., Reference Tsuruta, Misztal and Stranden2001).

In animal breeding application, the size of LHS can be much larger than the size of data to generate LHS. For example, one line of data in a 10 trait model with 10 effects per trait includes 20 numbers (one for each trait and effect) but generates 10 000 contributions to LHS. A special implementation of iteration algorithms is matrix-free or by iteration on data (Schaeffer and Kennedy, Reference Schaeffer1986; Misztal and Gianola, Reference Misztal1987), where coefficients of LHS are regenerated from the data every round of iteration. Although the implementation of iteration on data by GS can be complicated as the coefficients need to be recreated row by row, such an implementation in PCG is trivial as PCG only requires a product of B by a vector (Bq, with q a given vector), with no individual elements necessary (except for diagonals).

Tricks and rules in mixed model and genomic computations

Convergence properties with various systems of equations

When mixed models are solved by iteration (GS or PCG), the convergence rate is better with fewer equations per phenotype. Experiences indicate that the convergence rate is very good with sire models, much slower with animal models, and could be especially slow with young animals not tied to phenotypes or progeny.

The regular MME equations involve an inverse of A or G, e.g.

$$\left( {{\bf{I}} {\plus}{\bf{A}} ^{{{\minus}1}} \alpha } \right) {\bf{ {\vskip -3pt \hskip 2pt {\hat }}{\hskip -2.5pt a}\,{\equals}\;y}} $$

An equivalent system of equation avoids creating an inverse (Strandén and Garrick, Reference Strandén and Christensen2009):

$$\left( {{\bf{A{\plus}I}} \alpha } \right){\bf{\bf{ {\vskip -3pt \hskip 2pt {\hat }}{\hskip -2.5pt a}\;{\equals}\;Ay}} $$

The convergence rate in such models is poor and may not occur for larger systems of equations. This is because A ⁻¹ ‘connects’ only parents and progeny resulting in ‘sparse’ equations, whereas all animals in A are connected, resulting in dense equations. As G differs little from A (VanRaden, Reference VanRaden2008; Wang et al., Reference Wang, Misztal and Legarra2014), the same thinking applies to models with G and G ⁻¹.

Computations in SNP BLUP (and BayesX)

SNP BLUP and its various forms (including BayesA, BayesB, etc.) are very popular models in genomic computations. When the number of SNP is small, the MME for SNP BLUP can be created explicitly. As the memory and compute requirements are quadratic, explicit storage of SNP BLUP equations is unfeasible with large number of SNP. For example, memory requirements for MME with 1 million SNP would be 10¹² elements also with 10¹² elements of MME created from each genotyped animal. Various options in such a case have been explored by Legarra and Misztal (Reference Golden, Fernando and Garrick2008) with a data set on about 2 k individuals and 20 k SNP markers. Solution by Cholesky decomposition took 2 h, and computations would increase cubically with the number of SNP, whereas memory would increase quadratically. With solution by GS and straightforward iteration on data, computations took about 4 days, with linear memory and quadratic memory costs. With a special GS implementation called GSRU, where currently computed solutions are used to adjust residuals for all individuals, computing took about 1 min and both computing and memory costs were linear. Finally, with PCG iteration on data, computing took only 20 s, with approximately linear costs for memory and computations. Although GSRU is the only realistic algorithm with millions of SNP for Bayesian models, the PCG approach can be useful for models where SNP variances are derived from SNP solutions (Zhang et al., Reference Zhang, Liu, Ding, Bijma, de Koning and Zhang2010; Sun et al., Reference Sun, Qu, Garrick, Dekkers and Fernando2012; Wang et al., Reference Wang, Misztal, Aguilar, Legarra and Muir2012). A hybrid method would use PCG iteration with inner (non-Monte Carlo-based) updates of weights (VanRaden, Reference VanRaden2008).

Variations of single-step equations

Due to the prohibitive cubic cost of inversion of G with many genotyped animals, alternative equations were proposed for ssGBLUP. Equations in Misztal et al. (Reference Misztal and Gianola2009) used H rather than H ⁻¹. Convergence was achieved only with a small number of animals. Equations by Legarra and Ducrocq (Reference Legarra and Ducrocq2012) used G and A ₂₂, which are easy to use, instead of G ⁻¹ and A ₂₂ ⁻¹. The convergence rate was slow with medium models and not achieved with large models (Aguilar et al., Reference Aguilar, Misztal, Legarra and Tsuruta2013). Liu et al. (Reference Liu, Goddard, Reinhardt and Reents2014) proposed equations with SNP effects for genotyped animals. They reported lack of convergence with real data. All these models illustrate computing advantages with inverses of dispersion matrices (A ⁻¹ or G ⁻¹). As SNP BLUP models have good computing properties, Fernando et al. (Reference Fernando, Dekkers and Garrick2014) looked into ssGBLUP where ungenotyped animals were imputed and only SNP effects (plus imputation errors for ungenotyped animals) were estimated. Although convergence rates were good, imputations required a supercomputer and imputed genotypes required large memory for big populations. Better computing times were reported using a computer with fast graphic cards (Golden et al., Reference Legarra, Aguilar and Misztal2016). Meuwissen et al. (Reference Meuwissen, Svendsen, Solberg and Ødegård2015) presented a new type of H matrix based on (possibly fractional) imputation of ancestors by segregation analyses, followed by a sophisticated tailoring of G to A to compensate for errors in imputation. Although the cost of segregation analysis is high, it is not yet clear whether such an approach offers accuracy benefits when most ancestors are genotyped, and the process of tailoring is rather complex. All non-traditional ssGBLUP would require new research and investment in code for implementation of more complex models, such as multiple-trait, random regression or maternal models.

Sparse genomic relationship matrix

Computations with ssGBLUP would be easy for large number of animals if G ⁻¹ was sparse, like A ⁻¹. Past studies suggested that G had a limited rank, as inversion is unstable with >5 to 10 k animals, and Gs as used in genetic evaluations are blended with A ₂₂ (VanRaden, Reference VanRaden2008; Aguilar et al. Reference Aguilar, Legarra and Misztal2010) to achieve full rank. The eigenvalue decomposition of G is:

$${\bf{G}} \,\,{\bf{{\equals}}} \;{\bf{UDU}} '$$

where U is a matrix of eigenvectors and D the matrix of eigenvalues. If all eigenvalues are positive, the inverse of G is

$${\bf{G}} ^{{{\minus}1}} \,\,{\bf{{\equals}}} \;{\bf{U}} '{\bf{D}} ^{{{\minus}1}} {\bf{U}} $$

If some eigenvalues are 0, the inverse does not exist. If, due to noise, many eigenvalues are close to 0 and carry no information, they can cause instability of the inverse through large elements of D ⁻¹.

Let D _t indicate a fraction of D with non-negligible eigenvalues, and let U _t be corresponding eigenvalues. Then

$${\bf{G}} ^{{\minus}} \,{\bf{{\equals}\;U}} _{t}^{'} {\bf{D}} _{t}^{{{\rm {\minus}1}}} {\bf{U}} _{t} $$

If G ⁻¹ is to be created explicitly, the computing cost is cubic with quadratic storage. If G ⁻ is used by a PCG algorithm as

$${\bf{G}} ^{{\minus}} \,{\bf{q\;{\equals}}}\, \left( {{\bf{U}} _{t}^{'} \left( {{\bf{D}} _{t}^{{{\rm {\minus}1}}} \left( {{\bf U}_{t}^{{}} {\bf q}} \right)} \right)} \right)$$

and the number of eigenvalues in D _t is small, G ⁻ q can be calculated and stored at linear cost with respect to the number of genotyped animals. However, computing U and D is expensive and requires storage of G in full.

When G has a limited dimensionality, Misztal et al. (Reference Misztal, Legarra and Aguilar2014) proposed to calculate G ⁻¹ using only a fraction of G and without eigenvalue decomposition. Let

$${\bf{G\,{\equals}}}\, \left[ {\matrix{ {{\bf{G}} _{{cc}} } & {{\bf{G}} _{{cn}} } \cr {{\bf{G}} _{{cn}} } & {{\bf{G}} _{{nn}} } \cr } } \right]$$

where index c denotes core animals and index n non-core animals. In the algorithm for proven and young animals (APY)

$${\bf G}^{{{\minus}1}} \,\approx\,\left[ {\matrix{ {{\bf G}_{{cc}}^{{{\minus}1}} } & {\bf 0} \cr {\bf 0} & {\bf 0} \cr } } \right]{\plus}\left[ {\matrix{ {{\minus}{\bf G}_{{cc}}^{{{\minus}1}} {\bf G}_{{cn}} } \cr {\bf I} \cr } } \right]{\bf M}^{{{\minus}1}} \left[ {\matrix{ {{\minus}{\bf G}_{{nc}} {\bf G}_{{cc}}^{{{\minus}1}} } & {\bf I} \cr } } \right],$$

where M is a diagonal matrix with elements $m_{i} \,{\equals}\,g_{{ii}} {\minus}{\bf G}_{{ic}} {\bf G}_{{cc}}^{{{\minus}1}} {\bf G}_{{ci}} $ for individual i in the core group. Inversion of G by APY has a cubic cost (and quadratic memory) for core individuals and a linear cost (and linear memory) for non-core individuals. Figure 1 shows an example non-zero structure of a regular inverse (after blending to make it positive-definite) and different APY inverses. The inverses are diagonal for non-core animals. As the regular G is singular for a large population, an infinite number of generalized inverses exist, including those by the APY algorithm with different sets of core animals.

Figure 1 Non-zero pattern of a regular (a) and APY (b to d) inverses of the genomic relationship matrix with different choices of core animals.

Core animals can be randomly selected (Fragomeni et al., Reference Fragomeni, Lourenco, Tsuruta, Masuda, Aguilar, Legarra, Lawlor and Misztal2015), with an optimal number of about 12 000 for Holsteins (Masuda et al., Reference Masuda, Tsuruta, Aguilar and Misztal2016) and 10 000 for Angus (Lourenco et al., Reference Lourenco, Tsuruta, Fragomeni, Masuda, Aguilar, Legarra, Bertrand, Amen, Wang, Moser and Misztal2015). The number of core animals is linked to effective population size (Misztal, Reference Mrode2016). Optimal number of core animals is equal to the number of the largest eigenvalues in D explaining 98% of the variance in G (Pocrnic et al., Reference Pocrnic, Lourenco, Masuda, Legarra and Misztal2016a). With such a number, prediction of breeding values with APY was more accurate than with a regular inverse. Computation of G ⁻¹ for 570 k genotyped Holsteins by APY took only about 2 h, whereas a month computing with a large memory would be required for a regular inverse (Masuda et al., Reference Masuda, Tsuruta, Aguilar and Misztal2016). The APY algorithm removes computing limits from inversion of G and subsequently ssGBLUP.

Computing of A₂₂

ssGBLUP as in Aguilar et al. (Reference Aguilar, Legarra and Misztal2010) requires $${\bf A}_{{22}}^{{{\minus}1}} $$ . Although A ⁻¹ is sparse, $${\bf A}_{{22}}^{{{\minus}1}} $$ may be relatively dense (Faux and Gengler, Reference Faux and Gengler2013), and its storage for potentially millions of genotyped animals impossible. This matrix can be calculated indirectly (Strandén and Mäntysaari, Reference Strandén and Mäntysaari2014).

$${\bf{A}} _{{22}}^{{{\minus}1}} \,{\bf{{\equals}\,A}} ^{{22}} {\minus}\left( {{\bf{A}} ^{{12}} } \right)'\left( {{\bf{A}} ^{{11}} } \right)^{{{\minus}1}} {\bf{A}} ^{{12}} $$

where all submatrices are sparse and can be stored explicitly. With PCG iteration, explicit $$A_{{22}}^{{{\minus}1}} $$ is not required, and its product with a vector q can be calculated every round as follows (Masuda et al., Reference Masuda, Tsuruta, Aguilar and Misztal2016):

$${\bf{A}} _{{22}}^{{{\minus}1}} {\bf{q\,{\equals}}}\, \left[ {{\bf{A}} ^{{22}} {\minus}\left( {{\bf{A}} ^{{12}} } \right)'\left( {{\bf{A}} ^{{11}} } \right)^{{{\minus}1}} {\bf{A}} ^{{12}} } \right]{\bf{q}} $$

In particular

$${\mib{s}} \,{\bf{{\equals}\, }} \left( {{\bf{A}} ^{{11}} } \right)^{{{\minus}1}} {\bf{A}} ^{{12}} {\bf{q}} $$

is computed as a solution to:

$${\bf{A}} ^{{11}} {\mib{s}} \,\,{\bf{{\equals}\,\, A}} ^{{12}} {\bf{q}} $$

As all matrices are sparse, the cost of above calculations is small (Masuda et al., Reference Masuda, Tsuruta, Aguilar and Misztal2016).

GREML and APY

GREML is REML using a GRM. Similarly, REML using matrix H can be called single-step GREML (ssGREML). In the past, the success with REML in animal populations was due to sparse inverse of A ⁻¹ and sparse matrix factorization/inversion (Misztal and Perez-Enciso, Reference Misztal1993; Pérez-Enciso et al., Reference Pérez-Enciso, Misztal and Elzo1994). If G is dense and large, Masuda et al. (Reference Masuda, Misztal, Tsuruta, Legarra, Aguilar, Lourenco, Fragomeni and Lawlor2015) showed that GREML with the old sparse matrix package is inefficient and unstable. A new package called YAMS recognizes dense blocks in MME, rearranges computations accordingly, allowing also for parallel computations for large dense blocks. With YAMS, GREML was much faster and more reliable.

Additional issues

This paper focused on algorithmic issues relevant to current mixed model computations in genomics in animal breeding. Many other issues exist and many will become relevant. For example, genotyping of animals with small amount of information leads to greater importance of avoidance of double counting in multi-step evaluations, or of more accurate match between genomic and additive relationships in ssGBLUP. With efforts to use sequence data, important operations are aligning short fragments, imputation, genomewide association studies, and incorporations of the detected variants in the evaluation. In dairy, a de facto global breeding scheme for Holsteins creates a need for a comprehensive genomic evaluation by Interbull that has a high accuracy, yet does not compromise privacy of the member countries. In general, while brute-force approaches may be formidable, as or more accurate yet inexpensive approaches may exist that exploit peculiarities of data including a limited effective population in farm populations. For instance, the sequence includes three G base pairs, but farm populations can be described by about 4000 to 20 000 chromosomal segments (Pocrnic et al., Reference Pocrnic, Lourenco, Masuda and Misztal2016b). This means a block size of about 150 to 750 kbases with hard to differentiate SNPs, and ability to reduce costs by working with blocks rather than individual SNP. This also suggests that increases of genomic accuracies past 4000 to 20 000 of high accuracy genotyped animals are limited (ignoring G×E and decays of predictivity over time). Similarly, if Interbull countries are unwilling to share genotypes on individual bulls or SNP effects derived from their population, they can contribute GEBV and genotypes of artificially generated animals (Maantysaari, 2015, personal communication) with their number and accuracies relevant to each country.