An efficient exact method to obtain GBLUP and single-step GBLUP when the genomic relationship matrix is singular.

BACKGROUND: The mixed linear model employed for genomic best linear unbiased prediction (GBLUP) includes the breeding value for each animal as a random effect that has a mean of zero and a covariance matrix proportional to the genomic relationship matrix ([Formula: see text]), where the inverse of [Formula: see text] is required to set up the usual mixed model equations (MME). When only some animals have genomic information, genomic predictions can be obtained by an extension known as single-step GBLUP, where the covariance matrix of breeding values is constructed by combining the pedigree-based additive relationship matrix with [Formula: see text]. The inverse of the combined relationship matrix can be obtained efficiently, provided [Formula: see text] can be inverted. In some livestock species, however, the number [Formula: see text] of animals with genomic information exceeds the number of marker covariates used to compute [Formula: see text], and this results in a singular [Formula: see text]. For such a case, an efficient and exact method to obtain GBLUP and single-step GBLUP is presented here.
RESULTS: Exact methods are already available to obtain GBLUP when [Formula: see text] is singular, but these require working with large dense matrices. Another approach is to modify [Formula: see text] to make it nonsingular by adding a small value to all its diagonals or regressing it towards the pedigree-based relationship matrix. This, however, results in the inverse of [Formula: see text] being dense and difficult to compute as [Formula: see text] grows. The approach presented here recognizes that the number r of linearly independent genomic breeding values cannot exceed the number of marker covariates, and the mixed linear model used here for genomic prediction only fits these r linearly independent breeding values as random effects.
CONCLUSIONS: The exact method presented here was compared to Apy-GBLUP and to Apy single-step GBLUP, both of which are approximate methods that use a modified [Formula: see text] that has a sparse inverse which can be computed efficiently. In a small numerical example, predictions from the exact approach and Apy were almost identical, but the MME from Apy had a condition number about 1000 times larger than that from the exact approach, indicating ill-conditioning of the MME from Apy. The practical application of exact SSGBLUP is not more difficult than implementation of Apy.

Background

In animal breeding, two equivalent mixed linear models have been used for genomic prediction using phenotypes on genotyped individuals [1]. In the first, random effects of markers are explicitly included in the model [2, 3]. We will refer to this model as the marker effects model (MEM). In the second, the breeding value of each animal, which is a linear combination of the random marker effects, is included as a random effect [1, 2, 4, 5]. We will refer to this model as the breeding value model (BVM). The mixed model equations (MME) that corresponds to the MEM has order , where p is the number of non-genetic effects and k is the number of marker covariates, and the MME that correspond to the BVM has order , where is the number of animals. When genomic data were first available, the number of animals with genotypic and phenotypic records was much smaller than the number k of marker effects. Thus, genomic prediction with the BVM was more efficient than using the MEM [1, 5], and prediction using this approach is now known as GBLUP.

However, at present, in some livestock species such as dairy cattle, has increased to over 100,000 if not 1 million. When exceeds k, the matrix of genomic relationships will have at least eigen values that are zero, and therefore is guaranteed to be singular. In practice, depending on the effective population size, some of the smallest of the k largest eigen values may be very near to zero if not zero. In either event, the MME that require the inverse of cannot be employed to obtain GBLUP. In that situation, an alternate form of the MME [4, 6–8] that can accommodate a singular can be employed, but this results in a completely dense set of MME of order . Thus, when is large, this formulation of the MME is not useful for computing GBLUP. An alternative is to use a modified matrix obtained from by adding a small value to all its diagonals or by regressing it towards the pedigree-based relationship matrix, , so that it retains full rank, but this is no longer an exact representation of the model if the markers completely explain the breeding values. Furthermore, this modified relationship matrix still has a dense inverse, which may be impossible to compute when is large.

Suppose the rank of is . Then, we will show here how to obtain exact GBLUP without approximation from a set of MME that has order , which can be much lower than We also show how this approach can be used to obtain exact single-step GBLUP without approximation when some animals have not been genotyped. These formulations are useful to better understand predictions that are obtained by using the recursive algorithm for “parents (core)” and “young (noncore)” animals i.e. Apy, which is gaining popularity [9-13] as an approach to approximate the inverse of [9] or [13]. The exact inverse of the nonsingular matrix will be dense whereas Apy approximates this with a sparse inverse [9, 10]. We will show here that when a full-rank is obtained by adding a small number to the diagonals of only noncore animals, the inverse calculated in Apy for a suitable choice of core animals will be sparse and an exact inverse of , but the inverse may be ill conditioned. The approximate inverse calculated in Apy cannot ever be that of which is singular when . The Apy algorithm will never yield exact GBLUP predictions contrary to the claims in [9, 11], but it has been demonstrated to be a useful approximation for some choices of [11-13].

Theory

Let denote the centered marker genotype covariate matrix of order with which is the case when the number of genotyped animals is larger than the number k of marker covariates. Then, the row rank r of is [14]. Suppose is ordered such that its first r rows are linearly independent and are denoted . It follows that the remaining dependent rows of , denoted can be written as a linear combination:so thatNow, a commonly-used form of the genomic relationship matrix [5] becomeswhere it can be seen that the last rows are a linear combination of the first r rows. The last columns of are similarly a linear combination of the first r columns. Thus, in this case, is singular and its inverse does not exist. It can be seen from (2) that can be written as:where and

GBLUP when G is singular

In the following, we will assume that the vector of breeding values of animals can be adequately modeled as:where the vector of marker effects is assumed to have zero mean and covariance matrix It follows that the covariance matrix of the breeding values is:where To proceed, we further assume the following mixed linear model for the vector of phenotypic values:where is a vector of non-genetic fixed effects, and are incidence matrices relating and to and is a vector of residuals with zero mean and covariance matrix Here, we have assumed that the markers fully explain the breeding values. If this is not the case, a random polygenic residual effect with zero mean and covariance matrix that is proportional to can be included in the model.

Strategy I

When is singular, one strategy to get the BLUP of is to use the formula:where is a solution to the systemand [7]. When is large, this strategy is not computationally feasible because the matrix is dense, has order , and its inverse is needed in (7) and (8).

Strategy II

Another strategy is to get the solution to the following MME as proposed by Harville [6]:where These MME are dense and have order . Thus, although the above approaches do not require inverting , explicitly using these MME do not provide a feasible approach as the number of genotyped animals approaches or exceeds a million, because storing and solving such large and dense system of equations would exceed the capacity of the typical computer used for genetic evaluation. An implementation with iteration on data using the PCG algorithm may be feasible by computing matrix products like in parts as [15]. However, Aguilar et al. [16] reported these asymmetric equations do not scale up well and suffer convergence problems.

Strategy III

We show here that it is possible to obtain BLUP of by solving a set of MME that has order , which can be much lower than . To do so, the breeding values of the r animals with genotypes is denoted and the breeding values of the animals with genotypes is denoted The model for the breeding values in (4) can be written as:and writing as in (1), this becomes:Note that the vector of breeding values given by (11) is identical to (4), and thus these two vectors have the same covariance matrix that is given by (5).

Now, using (11) for in (6), the mixed linear model for the phenotypic values can be written in terms of as:The random effect of this model has order r and can be much lower than the order of Furthermore, as the models given by (6) and (12) have the same first and second moments, and thus they are equivalent models and yield the same BLUP for [7]. The MME for the model (12) arewhere is the inverse of the non-singular matrix , and The BLUP of is obtained as

Strategy IV

A key assumption in Strategy III is that the matrix of marker covariates can be reordered such that the first r rows are linearly independent and the remaining dependent rows can be expressed as a linear combination of the first set of r linearly independent rows. Determining the precise rank of may be inexact as the eigen values of decay slowly [17]. On the one hand, if the chosen contains less rows than the rank of , it would not be possible to express as . On the other hand, if contains more rows than the rank of will be singular. Even when the number of rows in is equal to the rank of may be ill conditioned if the smallest eigen value of is close to zero. The condition number of a matrix is represented by the ratio of the largest to the smallest eigen value, and it is 1 for a perfectly conditioned matrix and a large number for an ill-conditioned matrix. There are many combinations of individuals that can be placed in , but the condition number of the resultant may vary greatly according to the chosen combination. The condition number of will impact the condition number of the resultant MME, and poorly-conditioned equations take longer to solve iteratively than well-conditioned equations. In comparing the choice of core used in Apy in a pig evaluation, Ostersen et al. [18] reported similar numbers of PCG iterations for non-genomic analyses and 8 choices of core, but the correlation between the Apy-SSGBLUP and SSGBLUP ranged from 0.93 to more than 0.99 for genotyped animals. That paper did not report the criterion used to determine PCG convergence.

One way to improve the condition number of the MME is to fit an equivalent model obtained by orthonormalizing the rows of . Suppose , where . Then, the transformed vector of breeding values will have a genomic covariance matrix:Then, as in [17], formulating the model in terms of , which has a well-conditioned covariance matrix, will result in a well-conditioned MME.

Another way to improve the condition of the MME without explicitly reordering is by using an RQ decomposition [19] that involves expressing as where is a lower triangular matrix and is a orthogonal matrix. The RQ decomposition applies to the rows of a matrix in the same manner that the QR decomposition is applied to the columns. Exploiting the decomposition, the model equation for the phenotypic values can be written in terms of as:where now The MME for this model are:and predictions for all individuals on the original scale can be obtained as . Also, the marker effects can be obtained as . Note that this factorization does not require us to know or determine the rank of . Furthermore, the orthogonal matrix can be obtained by applying the RQ factorization to just the first k rows of , for which the number of operations is proportional to [19]. The matrix can be obtained as .

Comparison to Apy-GBLUP

The efficient algorithm to obtain the inverse of the additive relationship matrix is based on the property that the additive relationships between an animal and any non-descendant (an individual that is not a descendant) can be written as a linear combination of the relationships between the non-descendant and the parents of the animal [20, 21]. This property of additive relationships also allows construction of the additive relationship matrix by the tabular method [22]. The so-called Apy algorithm [9, 10] attempts to extend this idea to genomic relationships by classifying animals into two groups: “core” and “noncore” animals. The Apy algorithm seems to imply that the relationship between a noncore animal and any other animal can be written as a linear combination of relationships between the other animal and the animals in the core group. We will refer to this property of the genomic relationships that is required for Apy as the Apy property. Provided this property holds, it is claimed that Apy results in an efficient inverse of that leads to exact calculations of GBLUP [9, 11]. However, when is singular and cannot have an inverse. Thus, Apy-GBLUP cannot be exact.

To better understand the matrix portrayed as an inverse by the Apy algorithm, the genomic relationship matrix is partitioned into the core and noncore animals as follows:where the subscripts c and n denote the core and noncore animals. The Apy algorithm implies that is nonsingular and that can be written as where . Similarly, . Now, using these results, can be written as:whereAssuming is nonsingular, the inverse of can be obtained as follows. We start by expressing as:Then, the inverse of can be written as:which is identical to the formula given in Misztal et al. [9] provided is diagonal.

To examine the situation when and is singular, suppose the animals with genotypes in are considered as the core animals and those with genotypes in are considered as the noncore animals. Then, , and , and where from (2) can also be written as . Given this definition of the core and noncore animals, it can be seen from (2) that the Apy property holds, because the rows of for the noncore animals is a linear function of those in the core. Furthermore, it can be seen from (2) that:which shows that is null, and thus (16) cannot be computed. In this situation, can be replaced by for a positive scalar s. Then, (16) gives the exact inverse for a matrix of modified genomic relationships that is obtained by adding s to only the diagonals of the noncore group. If the scalar s is chosen to be small, will be close to . Regardless of the size of s, the resulting inverse is sparse because the sub-matrix corresponding to in the inverse has non-zero elements only on the diagonal. If the core group is chosen such that has rank less than r the rank of , the matrix will not be null, but as can be seen by examining Eq. (2) and demonstrated in the numerical example, it is not likely to be diagonal as assumed in the Apy algorithm. In this case the inverse computed by the Apy algorithm is the inverse of:where sets all the off-diagonal elements of to zero, thus always leading to an approximation for that choice of core. Also, when is blended with will generally not be diagonal (see Additional file 3), and the inverse obtained in Apy is of given by (17), where the off-diagonal elements of have been set to zero.

Numerical example

A small example with seven animals is used to illustrate the calculation of GBLUP and Apy-GBLUP. The pedigree for the seven animals is in Table 1. Genotype covariates coded as −1, 0, 1 at four loci are in Table 2. Julia scripts and results for GBLUP by strategies I to IV and for Apy-GBLUP are in Additional file 1. Only the calculations by strategy III and by Apy-GBLUP are described below.

Table 1

Pedigree for numerical example

Animal	Sire	Dam	PV	BV	EBV
1	0	0	99.25	−0.25	0.14
2	0	0	97.92	−0.94	−0.95
3	0	0	103.2	1.12	1.09
4	1	2	99.39	−1.01	−0.69
5	1	2	102.03	0.79	0.25
6	1	3	100.59	0.18	0.14
7	1	3	101.7	1.55	1.08

PV, BV and EBV are the phenotypic values, breeding values and the BLUPs of the BV

Table 2

Genotype covariates at four loci

Animal	Locus 1	Locus 2	Locus 3	Locus 4
1	0.0	0.0	−1.0	0.0
2	−1.0	1.0	0.0	0.0
3	1.0	0.0	−1.0	0.0
4	−1.0	0.0	0.0	1.0
5	0.0	1.0	0.0	1.0
6	0.0	1.0	−1.0	0.0
7	1.0	1.0	−1.0	0.0

The first step in this approach is to reorder the rows of such that the first r rows are linearly independent, where r is the rank of . As described below, this can be done using Gaussian elimination with pivoting on to transform it to row echelon form, where all elements below the diagonal are zero. Starting in row , zeros are obtained below the diagonal by subtracting a multiple of row i from each subsequent row. Before doing these row operations to obtain zeros under the diagonal, the element with the largest absolute value is located in the sub-matrix comprising all rows below row and all columns to the right of column . Then by swapping rows and columns, this element is moved to the diagonal. If the element with the largest absolute value is zero, Gaussian elimination is terminated. The rank of the matrix is the number of non-zero diagonals in the transformed matrix, and the rows used for Gaussian elimination provide a maximal set of linearly independent rows.

When Gaussian elimination was applied to genotype covariates in Table 2, the resulting matrix is in Table 3. All four diagonals of this matrix are non-zero, and so has a rank equal to four. As a result of swapping rows, the rows were ordered as 2, 7, 1, 4, 5, 6, 3, where rows 2, 7, 1 and 4 were used for Gaussian elimination. Thus, these four rows are a linearly independent set, and they were taken to form and rows 5, 6, and 3 were taken to form The genomic relationship matrix that was constructed using the reordered genotype covariates is in Table 4. The upper-left, sub matrix of this relationship matrix, denoted , gives the relationships for individuals 2, 7, 1, and 4, and it is non-singular because the genotype covariates for these four individuals are linearly independent. Now, the matrix can be calculated using (3), and is in Table 5.

Table 3

Genotype matrix transformed to row echelon form by Gaussian elimination with pivoting

−1.0	1.0	0.0	0.0
0.0	2.0	−1.0	0.0
0.0	0.0	−1.0	0.0
0.0	0.0	0.0	1.0
0.0	0.0	0.0	0.0
0.0	0.0	0.0	0.0
0.0	0.0	0.0	0.0

Table 4

Genomic relationship matrix

0.5	0.0	0.0	0.25	0.25	0.25	−0.25
0.0	0.75	0.25	−0.25	0.25	0.5	0.5
0.0	0.25	0.25	0.0	0.0	0.25	0.25
0.25	−0.25	0.0	0.5	0.25	0.0	−0.25
0.25	0.25	0.0	0.25	0.5	0.25	0.0
0.25	0.5	0.25	0.0	0.25	0.5	0.25
−0.25	0.5	0.25	−0.25	0.0	0.25	0.5

Table 5

The matrix that relates to as

0.0	1.0	−1.0	1.0
0.5	0.5	0.5	0.0
−0.5	0.5	0.5	0.0

The last three rows of the matrix of genomic relationships in Table 4 can be written as a linear combination of the first four rows as shown in (2), by using the matrix in Table 5. Thus, breeding values for individuals 5, 6, and 3 can be written as:where is the vector of breeding values for individuals 2, 7, 1, and 4. Now, the phenotypes for these seven individuals are modeled in terms of , which has a non-singular covariance matrix proportional to All seven individuals in this example have one phenotypic value, and so assuming that the vector of fixed effects contains a single element for the overall mean, the matrix for this example is equal to a vector of seven 1s and is equal to an identity matrix of order seven. It follows that for in Table 5. The MME to fit the overall mean ( and the breeding values are in Table 6, where a value of 1.0 was used for BLUP of is obtained as . Results for strategies I through IV are in Additional file 1, and they are all identical, as expected. The condition numbers of the left-hand-side of the MME for strategies II through IV were 10.9, 11.3, and 6.8, demonstrating the improved condition of the MME obtained by fitting a RQ transformed vector of breeding values.

Table 6

Mixed model equations for and

	μ	u_2	u_7	u_1	u_4
μ	7.0	1.0	3.0	1.0	2.0
u_1	1.0	4.5	−1.0	1.0	−2.0
u_7	3.0	−1.0	5.5	−3.5	3.0
u_1	1.0	1.0	−3.5	9.5	−3.0
u_4	2.0	−2.0	3.0	−3.0	6.0
rhs	704.08	96.62	305.62	99.12	201.42
sol	100.43	−0.95	1.08	0.14	−0.69

The last two rows give the right-hand-side and the solutions of the equations

Apy-GBLUP

Here, we can see that if animals 2, 7, 1, and 4 are used as the core group, the Apy property is met because the last three rows of , which correspond to the animals in the noncore group, can be written as a linear combination of the first four rows, which correspond to the animals in the core group, using the matrix in Table 5 (see Additional file 1). Equation (2) shows that this property also holds for the columns of , where the last three columns of can be written as a linear combination of the first four columns. In this case, the matrix , the inverse of which is needed in the Apy algorithm, is null (Additional file 1). In order to proceed with the Apy algorithm, we set for a small value of s such as 0.0001. The inverse that is obtained from equation (16) will now be sparse because the sub-matrix corresponding to in the inverse is diagonal. Inverting a modified matrix, , by adding s to the diagonals of corresponding to the animals in the noncore group gives the same result (Additional file 1). Setting up and solving the MME for and assuming give results that are approximate but very close in this instance to the exact BLUP results obtained by strategies I through IV (Additional file 1), but the condition number of these MME was 56,548, which indicates that they are ill-conditioned relative to those for strategies II through IV. However, if individuals 2, 7, and 1 are chosen as the core animals, the Apy property does not hold. In that case, the last four rows of cannot be written as a linear combination of the first three rows (Additional file 2). Furthermore, the matrix computed by using equation (15) is not diagonal (Additional file 2). Now, the matrix that is inverted in the Apy algorithm deviates substantially from , and as a result, solving the MME for and assuming gives results that are substantially different from the exact BLUP (Additional file 2).

Recent publications [12, 13, 18] in which the Apy algorithm was applied to obtain a matrix portrayed as the inverse of the genomic relationship matrix use rather than the singular . This approach applied to the example gives a solution that is neither the same as the exact solution obtained using any of the strategies I to IV (Additional file 1), nor the exact solution to the MME constructed with the blended genomic relationship matrix. However, the condition number of these equations was 62.1, which is much better than that obtained without blending but poorer than with strategies II through IV.

Exact single-step GBLUP when Ggg is singular

Single-step GBLUP (SS-GBLUP) was proposed [23, 24] to obtain genomic evaluations when genotypes are not available on all animals.

Let denote the breeding values of animals with genotypes and denote the breeding values of those without genotypes. Now, the mixed linear model for SSGBLUP can be written as:It is convenient to similarly partition the vector of phenotypic values as , where are phenotypic values from animals that were not genotyped and are from animals that were genotyped. However, because is singular, is written as in Eq. (11) in terms of , and then the model becomes:where , and . The MME that correspond to (18) are:whereTo obtain , as in [23], is written as:where, in the last line, we have used the identity: . Now, is written in terms of and as:where and , and is written asFollowing [23], , and . Unlike in [23], two variance components and are used here, where is the additive genetic variance and stems from the prior used for the marker effects (, and its relationship to the genetic variance may not be straightforward [25]. Finally, is as follows:whereFollowing [26], the MME given in equation (19) can be augmented to avoid the expression involving in . To do so, equation (19) is first rewritten to show the partitions for and . Note that the matrix has the following form:where relates to and relates to . Then, can be partitioned as:where . Now the MME that show the partitions for and are:where and are partitions of corresponding to and . Consider now the following augmented MME:These equations do not have in them, and so they may be easier to construct. However, the left-hand-side is not positive definite and it has been reported that these equations are poorly conditioned [26]. Elimination of from Eq. (21) results in equation (20), and thus, solutions for and for from Eq. (21) are identical to those from Eq. (20).

The model for SSGBLUP can also be formulated in terms of as:and the MME corresponding to model (22) are:where and .

Comparison to Apy-SSGBLUP

The SSGBLUP method given in [23] requires computing the inverse of the matrix of genomic relationships and of the matrix of additive relationships for the genotyped animals. At the time those papers were published, was typically smaller than the number of markers so that was relatively small and of full rank. Since then has greatly increased in most livestock applications. Computational effort in matrix manipulation is determined by the number of non-zero coefficients and these increase as increases. To fully store a dense matrix of order one million in single precision requires about 4 TB. Therefore, it would be advantageous to have a sparse representation of all the large matrices involved in the MME.

Furthermore, the matrix is singular when and thus cannot be inverted when more animals than the number of SNPs have been genotyped. This suggests that there should be a sparse representation of . Suppose has rank r and it is ordered such that the first r rows are linearly independent. Then, the sub-matrix of the first r rows and columns of denoted gives the genomic relationships among the r core animals of the Apy algorithm, and that sub-matrix is nonsingular. The remaining animals are referred to as noncore and their genomic relationship matrix is denoted . When the genomic relationship matrix has not been blended, a nonsingular matrix can be obtained by adding a small value to the diagonals of for the animals in the noncore group, and in the inverse of , the sub-matrix corresponding to will be diagonal. This exact inverse of that particular can be obtained efficiently using the Apy algorithm given in [9]. When the genomic relationship matrix is blended with the matrix resulting from the Apy algorithm is the exact inverse of the where all the off-diagonals of have been ignored. That matrix may or may not be a close approximation of the blended genomic relationship matrix depending on the size of the core, the particular animals chosen for the core, the relationship among noncore animals and the relationship between core and noncore animals. Regardless of the form of the matrix the Apy algorithm is applied to, the resultant inverse is sparse.

The MME for Apy-SSGBLUP, which includes the inverse of , are:In addition to the inverse , SSGBLUP requires the inverse of . However, is a dense matrix, and so subtracting it from the inverse of will make the resultant matrix dense.

Part of the appeal of the Apy algorithm was to obtain a sparse representation of the MME for SSGBLUP. Accordingly, Misztal et al. [9] proposed that Apy could also be used to approximate the inverse of the nonsingular . However, the nature of depends on the genotyping strategy such that genotyping unrelated individuals results in a diagonal whereas genotyping relatives results in non-zero off-diagonals between each related pair. If off-diagonal elements in the noncore sub matrix of are not well predicted by , the Apy inverse can significantly depart from its true inverse as easily demonstrated by using an example (see Additional file 1). This means the adequacy of Apy applied to will depend on the pedigree structure, the nature of the genotyping strategy, and the choice of core group. Presumably, this inadequacy of Apy for inverting has been recognized because recent implementations [13] have adopted an alternative approach that is computationally more demanding than applying Apy to approximate the inverse of . Rather than forming prior to solving the MME, a partitioned matrix inverse result is used to calculate products such as as , where is the solution to . This requires storing the sparse matrices and the sparse Cholesky factors of . Each PCG iteration involves a matrix product for a different vector , which requires one forward and one backward triangular solve to obtain , two sparse matrix vector multiplications, and one vector subtraction.

Both the MME for Apy-SSGBLUP and that for SSGBLUP using strategy IV (SIV-SSGBLUP), include equations for the same fixed effects and the random effects corresponding to the breeding values of animals that were not genotyped. In the MME for Apy-SSGBLUP, there is an additional vector of random effects corresponding to the breeding values for animals that were genotyped, which comprises sub-vectors representing core and noncore animals. In contrast, the MME for SSGBLUP using strategy IV contains a vector of random effects that is not larger than k regardless of the number of animals genotyped. If the core size in Apy-SSGBLUP was chosen to be k, Eq. (24) would contain an additional random effect of order equal to the number of noncore animals compared to Eq. (23), and this number increases with the number of animals genotyped.

Given a core of k animals, both MME contain a dense matrix on the diagonal. Both MME contain the same sparse block on the diagonal for non-genotyped animals. Comparing the upper off-diagonals of the two sets of symmetric MME, that for Apy-SSGBLUP has the sparse matrix whereas SIV-SSGBLUP has the product of that matrix with the mostly dense matrix . Rather than forming the dense product, matrix computations involving that matrix can be done more efficiently when in parts (e.g. ) storing only and in memory. The Apy-SSGBLUP MME contain on the upper diagonal a dense block that does not appear in SIV-SSGBLUP and which increases in size as more animals are genotyped. The computation required to form the diagonal block of SIV-SSGBLUP involves computing , where is column i of . This calculation is virtually identical to the computation of in Apy-SSGBLUP, but the former needs to be done for each genotyped animal once whereas the latter needs to be done for each PCG iteration.

Discussion

When the number of genotyped animals is larger than the number k of marker covariates, the matrix of genomic relationships becomes singular. In this situation, we have shown here how to obtain exact GBLUP without any approximation from either Eqs. (13) or (14) of order or , where is the rank of . The MME given by Eq. (9) can also be used to obtain GBLUP without approximation, but these asymmetric MME are of order . When more individuals are genotyped and grows, the order of those MME (9) also grows. In contrast, the order of the MME presented here (Eqs. 13 and 14) will remain constant even as grows.

An alternative to these exact GBLUP calculations is used in Apy-GBLUP. Here, the pedigree is divided into two groups of animals: the core group and the noncore group. We have shown here that the inverse computed in the Apy algorithm is for a modified genomic relationship matrix, where the sub-matrix of genomic relationships among the noncore group of animals is replaced by . If the core group is chosen such that the rank of is equal to the rank of will be null and the Apy algorithm will fail. In that case, the diagonals of can be set to some small value, but this can result in ill-conditioned MME as shown by the example in Additional file 1. The MME can be ill-conditioned even when is not null but contains very small values on the diagonal. Although the MME for Apy-GBLUP will also grow with , it contains a block that is diagonal, and thus is very sparse.

The approach presented here can also be used to obtain exact SSGBLUP when some animals are not genotyped. In contrast to the Apy algorithm, the method presented here is never an approximation. In agreement with [10], “BVs of core individuals can all be written as linear combinations of effective SNP effects” when SNP effects fully explain the BV. In contrast to the claim in [10] that “BVs of noncore individuals depend approximately only on the BVs of the core individuals” we have shown that the BVs of noncore individuals are an exact linear function of the BVs of the core individuals when the rank of is equal to the rank of This requires the core group to contain at least as many animals as the rank of When the number of genotyped animals exceeds the number of markers, will be singular and its rank cannot be greater than the number of markers. Only when the rank of is less than the rank of , will the “BVs of noncore individuals depend approximately only on the BVs of the core individuals”.

The Apy algorithm when applied to may or may not be a good approximation depending on the particular . It will be exact for any core if genotyped animals are all unrelated as the matrix is strictly diagonal.The quality of the approximation will erode with increases in the number of large-magnitude off-diagonal elements in . Demonstrating with real data that the Apy gives good results in one or more field data sets is no guarantee that it will perform well for all applications. This raises concerns that the same could be true for the application of Apy to the genomic relationship matrix. When the number of genotyped individuals increases and the number of core animals remains constant, there may be a large increase in the number of off-diagonal coefficients in . Those coefficients are ignored in the Apy algorithm, and the predictions approximated by Apy are expected to deviate further from the exact predictions as more coefficients are ignored. Thus, inference that the Apy algorithm based on 100,000 or 500,000 genotyped animals is appropriate cannot be extrapolated to similar data structures with a million or more animals genotyped.

If SNP effects do not fully explain the BV, an additional polygenic effect for all animals can be readily fitted in addition to , the breeding values explained by the markers for a subset of genotyped animals. Lourenco et al. [12] used default options of BLUP90IOD2, which means they blended with , and included competitive results from Apy compared to exact predictions obtained by direct inversion for various analyses with that were smaller than 52,000. Fragomeni et al. [11] limited their analyses to in order to allow direct inversion of a matrix based on 42,503 SNPs but do not mention whether blending was used. In the absence of blending, the rank of matrix could not exceed 42,503 and a direct inverse of does not exist. Pocrnic et al. [27] simulated genotypes on 75,000 individuals and blended with . They showed that Apy exceeded the accuracy of exact ssGBLUP by direct inversion. Their QTL effects were simulated from a Gamma distribution, which creates a few loci with large effects. In those circumstances, methods such as BayesB typically outperform GBLUP [3], and Apy with a small set of core individuals may similarly benefit from reducing the dimension of the model. Masuda et al. [13] blended with to guarantee nonsingularity of the blended matrix with greater than 500,000. The exact inverse of that blended matrix will be a dense matrix of order , which will make exact calculations computationally infeasible when exceeds about 150,000. This makes it impossible to compare the accuracy of Apy approximations to exact predictions using that approach. They show high correlations between approximations for different core definitions but the correlations between their approximations and the exact predictions are not known.

There are no published results demonstrating the comparative accuracy of Apy and the exact approach when is too large for direct inversion of However, using the exact SSGBLUP calculations presented in this paper such a comparison is feasible, requiring only special computation for and in the MME (23). Computation of the matrix involves the same calculations as required to impute genotypes for non-genotyped animals as presented in Fernando et al. [28]. The computation of R is straightforward and analogous to matrix that is fundamental to computations in the Apy algorithm. Accordingly, we do not consider that the practical application of exact SSGBLUP will be any more difficult than implementation of Apy.

Conclusions

When the number of genotyped animals exceeds the number of marker loci, the genomic relationship matrix cannot be full rank. We introduce an approach that partitions the genotyped animals into two sets, one of which can be referred to as core animals, and the other as non-core animals whose breeding values can be written as a linear function of the breeding values of core animals. The MME used for genomic prediction are then constructed with only the breeding values of the core animals, and with phenotpyes of the non-core animals contributing to the predictions for core animals through their linear relationships to the core animals. The estimated breeding values of the non-core animals are obtained as a linear function of the estimates of the breeding values of the core animals. This gives exact solutions for all animals. Another approach is to blend the genomic relationship matrix with a numerator relationship matrix or a scaled identity matrix to ensure the blended genomic relationship matrix is full rank. In that case, standard mixed model computing procedures can be used, but the increase in computing effort will be proportional to the cube of the number of animals genotyped. That effort can be reduced by approximating the inverse of the blended genomic relationship matrix using the Apy algorithm. That approximation also partitions the animals into core and non-core groups, but explicitly fits both sets of animals in the MME. In some cases, it has been reported that this approach gives useful approximations. However, the computing effort for that approximate approach is similar to that of the exact approach introduced here.