Inexpensive Computation of the Inverse of the Genomic Relationship Matrix in Populations with Small Effective Population Size.

Many computations with SNP data including genomic evaluation, parameter estimation, and genome-wide association studies use an inverse of the genomic relationship matrix. The cost of a regular inversion is cubic and is prohibitively expensive for large matrices. Recent studies in cattle demonstrated that the inverse can be computed in almost linear time by recursion on any subset of ∼10,000 individuals. The purpose of this study is to present a theory of why such a recursion works and its implication for other populations. Assume that, because of a small effective population size, the additive information in a genotyped population has a small dimensionality, even with a very large number of SNP markers. That dimensionality is visible as a limited number of effective SNP effects, independent chromosome segments, or the rank of the genomic relationship matrix. Decompose a population arbitrarily into core and noncore individuals, with the number of core individuals equal to that dimensionality. Then, breeding values of noncore individuals can be derived by recursions on breeding values of core individuals, with coefficients of the recursion computed from the genomic relationship matrix. A resulting algorithm for the inversion called "algorithm for proven and young" (APY) has a linear computing and memory cost for noncore animals. Noninfinitesimal genetic architecture can be accommodated through a trait-specific genomic relationship matrix, possibly derived from Bayesian regressions. For populations with small effective population size, the inverse of the genomic relationship matrix can be computed inexpensively for a very large number of genotyped individuals.

FOR animals and plants, many genomic analyses with SNP data use one of two approaches. Either effects of SNP markers are estimated with best linear unbiased prediction (SNP-BLUP) (Meuwissen ; VanRaden 2008; Gianola ; Piepho 2009) or a genomic relationship matrix (GRM) is used in genomic BLUP (GBLUP) (VanRaden 2008). Estimation of SNP effects makes SNP selection and estimation of SNP variance easy, leading to straightforward single-trait prediction and genome-wide association study (GWAS). GBLUP is easier to use in more complex models (e.g., multiple traits) and for parameter estimation because existing BLUP including parameter estimation methodology can be used, although the use of GBLUP for GWAS is more complex (Zhang ). For prediction, SNP-BLUP (possibly with SNP weighting) and GBLUP are equivalent models (VanRaden 2008) but they differ in computing cost. SNP estimation includes the same number of SNPs independent of the number of individuals. Adding extra individuals incurs linear computing costs and no additional storage. GBLUP usually requires an inverse of GRM, and explicit inversion requires quadratic memory and cubic computations.

When only a small fraction of the population is genotyped, GBLUP can be extended to single-step GBLUP (ssGBLUP) (Aguilar ; Christensen and Lund 2010). In this method, a numerator relationship matrix (NRM) for all individuals and a GRM are combined and then applied to BLUP. Benefits of ssGBLUP include simplicity of application (another BLUP), avoidance of double counting, and accounting for preselection on Mendelian sampling (Legarra ). However, ssGBLUP also requires an inverse of GRM.

The inverse of GRM can be computed with general algorithms only for up to perhaps 150,000 individuals because of memory and computing time limitations. However, the number of genotyped individuals across animal populations is expanding. In dairy cattle, >950,000 Holsteins have been genotyped in the United States as of November 2015 (Council on Dairy Cattle Breeding; https://www.cdcb.us/Genotype/cur_density.html). Several algorithms were proposed to lower the cost of ssGBLUP (Legarra and Ducrocq 2012; Fernando ; Liu ), but either they are computationally expensive or the algorithms do not converge with a large number of genotypes.

Past progress in animal breeding resulted to a large degree from a fast algorithm to invert the NRM (Henderson 1976). Although the cost of explicit inversion of the NRM is cubic with the number of animals, the cost of creating that inverse directly by recursion is very low (Henderson 1976; Quaas 1988). When animals are ordered from oldest to youngest, a recursion for each animal includes at most only two terms (one for each parent). Consequently, the inverse of the NRM can be created at a linear cost.

Faux applied recursions on relatives of genotyped individuals to the GRM; however, the inverse was not accurate. Misztal postulated recursions on “proven” animals (with their own phenotypes or their progeny phenotyped) and called the methodology an algorithm for proven and young (APY) animals. The APY was tested in a population of Holsteins with a total of 100,000 genotyped animals and different groups of animals in recursion (Fragomeni ). When recursions were on proven bulls only, the correlation of genomic estimated breeding values (GEBVs) for selection candidates with APY G−1 with GEBVs from a complete inverse was >0.99. When only cows were in the recursion, the correlation remained >0.99. When the recursion included random subsets of 5000, 10,000, and 15,000 animals, the correlations were 0.97, 0.98, and 0.99, respectively, with minimal variability among replicates. Moreover, the convergence rates with random subsets were superior, indicating better numerical conditioning. The APY was also applied to a commercial population of Angus cattle (Lourenco ). With recursions on 4000, 8000, and 33,000 animals, the APY accounted for 84%, 97%, and 100% of accuracy gains of ssGBLUP over BLUP, respectively. The APY computing and storage costs when applied to cattle are almost linear, which allows for inverting practically any GRM size. However, why the APY works and its possible internal limitations have not been addressed. The first purpose of this study was to develop a theory explaining why recursion on a limited number of individuals results in an accurate GRM inverse. The second purpose is to determine implications of that theory for populations other than cattle.

Methods

Recursions and the inverse of the numerator relationship matrix

Henderson’s (1976) inverse of the relationship matrix can be derived by recursion. Let u be a vector of additive effects or breeding values (BVs) distributed as where A is a numerator relationship matrix and, for simplicity, the additive variance is set to 1. Following developments in Misztal , the joint distribution of can be written asA notation below n:m denotes an index from n to m. Assuming normality, the conditional distributions arewith part of the ith row of A, and with the recursion equationwhere In matrix notationIn matrix notation, the recursions can be written aswhere M is a diagonal matrix. Because the inverse of A can be computed aswhich is a form of a Cholesky decomposition with diagonal M = diag{m) and lower diagonal p = {p} with entries defined for elements with indexes j ≤ i. Each different ordering of individuals will lead to a different matrix P but identical A−1. For random ordering P is likely dense, and the cost of the Cholesky decomposition with dense coefficients is cubic. Additional costs of the inverse will include calculation of A and P.

Henderson (1976) discovered rules to create A−1 at a low cost. Indirectly, the rules are based on a recursion (Quaas 1988),where s and d refer to the sire and dam of animal i, and φ is the Mendelian sampling. Subsequently, when animals are ordered from the oldest to the youngest,P has at most only two nonzero elements per row corresponding to parents and equal to 0.5, and the inverse can be calculated at a low and linear cost. Additionally, computing A is no longer needed as the nonzero elements of P are known and diagonal elements of M are easy to compute. Even though the inverse by Henderson (1976) is very simple, it is not an approximation. In practice, it may be more accurate than an inverse derived from inverting A explicitly because of fewer computations and thus lower rounding errors.

Recursions and the inverse of the genomic relationship matrix

Let u be distributed as where G is a genomic relationship matrix. The recursion equations are the same as previously,but with different coefficients,and a similar inverseThis inverse can be computed at a low cost only when P is sparse and only a small fraction of G needs to be computed. The next sections show that both conditions can be met for populations with limited effective population size.

Recursion and SNP model

Assume availability of an optimal number of SNP markers (called effective SNP) such that increasing that number would not increase the accuracy of prediction. Let with and be a SNP BLUP, where y is a vector of phenotypes (or phenotype equivalents), µ is population mean, a is a vector of SNP marker effects, Z is a centered matrix of gene content, e is a vector of residual effects, I is an identity matrix, is SNP marker variance, and is residual variance. If u is a vector of BVs, then where ε is the fraction of BVs unexplained by SNP effects with Applications to farm animals using medium-size SNP chips usually assume (VanRaden 2008; Goddard ).

Divide individuals arbitrarily into two groups: core individuals denoted as c and other (noncore) individuals denoted as n. Then and The conditional expectation of a|uc is BLUP prediction of SNP effects calculated from BVs of core animalswhere Let where is prediction error, is prediction error variance (PEV), and are independent. ThenIf SNP effects nearly fully explain BVs (),where is a matrix that relates breeding values of noncore to core individuals and Note that the combined error term has a nondiagonal variance but can be approximated as diagonal especially when the number of core individuals is equal to or greater than the number of SNPs ().

In the formula that relates noncore to core animals, using fewer core animals than the number of SNPs necessary would lead to increased PEV, and using more core animals than the number of SNPs should not affect PEV. Note that when the number of core animals is the same as the number of effective SNPs and Z is invertible, BVs of core animals contain almost the same information as these SNP effects orConsequently, BVs of core individuals act as linear combinations of effective SNP effects and BVs of noncore individuals depend approximately only on BVs of core individuals. The above formula is useful only for presentation as in practice the differences in GEBVs with slightly different numbers of core animals are minuscule (see Figure 1).

Figure 1

Correlations between genomic estimated breeding values (GEBVs) for selection candidates using regular and the APY inverse of the genomic relationship matrix (GRM) with various numbers of base individuals (Fragomeni ). Correlations are based on analysis of 10,102,702 final scores on 6,930,618 Holstein cows, with genotypes available on 100,000 animals; and correlations are based on GEBVs for 49,611 selection candidates.

Recursions and independent chromosome segments

For populations with limited effective population size (Ne), the genome is broken into a small number called homogenic or independent chromosome segments (ICS), with the number of segments inversely proportional to Ne (Stam 1980; Daetwyler ). If the number of ICS in a population is Me and each segment has an additive effect, the BV of each individual is a sum of effects of chromosome segments present in that individual. The following derivations use the previous derivations while substituting effective SNP effects by ICS effects.

Let s be a vector of additive effects of ICS. Assume that these effects explain nearly all the additive variance. Let t be a fraction of segment j in individual i, and assume that the value of t is t. Assuming a gametic model, t would take values of 0, 1, or 2, although this assumption is not critical here. Then where T is a matrix that relates u to chromosome segments, var(s) = is segment variance, and ε is the fraction of BVs unexplained by ICS effects.

Following the previous derivations, divide the individuals into two groups: core and noncore, and If the number of core animals is equal to Me, T is full rank, and or nearly all the additive information present in ICS is also present in BVs of core animals. Then,and like before, BVs of noncore individuals depend approximately only on BVs of core individuals.

Recursions and the APY formula

Based on previous derivations, when the number of core animals is sufficiently large, BVs of noncore animals depend only on BVs of core animals. However, obtaining the values of P and corresponding errors from the formula with effective SNPs or ICS is hard. If the GRM is available, the inverse can be obtained indirectly, by applying Equations 1 and 2. In (1), decompose the component matrices P and M into sections due to core and noncore animalsSubsequently, the complete recursion isas the term relating BVs of noncore animals to BVs of noncore animals (P) is null. Following Equation 2, the inverse of G isIf the inverse corresponding to core animals is available, and the complete inverse can be simplified toComputing P directly is impossible especially when the effective SNPs or ICS are not known; however, P can be computed indirectly from the GRM. Denotewhere G is a GRM and is the additive variance. Using conditional distributions, for individual i in the noncore group, and the inverse of G can be calculated asThe above equation is the same as that reported by Misztal for the APY with proven animals as the core group and young animals as the noncore group. The formula above requires that only G be full rank. When G is not of full rank, the APY inverse may in fact be a generalized inverse.

The above derivations can be simplified if the recursion includes only the noncore animals. ThenSubsequentlyandwhich leads to the same formula as derived previously.

Computing costs

Inversion of G by APY has a cubic cost (and quadratic memory) for core individuals and a linear cost (and linear memory) for noncore individuals. Savings in memory and computing are due to ignoring storage and computations for the blocks of noncore × noncore animals (except diagonal elements) for both G and APY (see Figure 2). Assume n core and p noncore individuals. Although regular requires ∼ memory and computations, APY requires only ∼2np memory and computations. If is the fraction of individuals that compose the core group and APY would require only ∼2β memory and computations of the regular algorithm. If n = 10,000 and n + P = 600,000, this is equivalent to 3% of the storage and 0.05% of the computations of the regular algorithm. When the number of core animals is limited, the APY effectively removes limits from computing G−1.

Figure 2

Sparsity pattern of a regular genomic relationship matrix (G) and its inverse (G−1) and elements of the genomic relationship matrix needed for construction of the APY (G for the APY) and the APY inverse (APY G−1).

Determination of effective number of SNP markers

Assume an SNP model with a very large number of possibly redundant SNP markers. The real dimensionality of the SNP information and subsequently the number of core animals required to account for all information in these markers can be determined by singular value decomposition (e.g., Wall ). In a formula that relates BVs to SNP effects (ignoring the error term), apply singular value decomposition Z = UΔV, where U and V are unitarian matrices (UU′ = I and VV′ = I) and Δ ={δ} is a diagonal matrix with singular values on the diagonal. ThenLet Δ be a matrix of rows of Δ where small singular values (say < ϕ) are zeroed. A well-chosen ϕ would retain the accuracy of BVswhile reducing the dimensionality of the SNP information; the number of nonzero values in can be called the effective number of SNPs and the nonzero fraction of could be called effective SNP markers. The parameter ϕ and the effective number of SNPs can be obtained from eigenvalues of the GRM. On the variance scalewhere the formula above shows eigendecomposition of G with ΔΔ′ as a diagonal matrix of eigenvalues. Following Janss , the average fraction of retained variance over all individuals with small eigenvalues set to 0 is proportional to the sum of retained eigenvalues. Subsequently φ2 can be chosen to retain a high fraction of the variance (say 0.98) bySummarizing, the number of effective independent SNPs and subsequently the minimum number of core animals can be derived from eigenvalue analysis of the GRM.

Genetic architecture and the GRM

While the derivations for the APY included SNP effects or effects of ICS, these effects are absent from final APY derivations, which depend on GRM only. Therefore, any information on specific architecture of a trait, if present, needs to be included in the GRM. In the GBLUP case (same variance of each SNP effect) the GRM can be derived from SNP BLUP as where q is a scaling factor (VanRaden 2008). For weighted SNP BLUP, where and D is a diagonal matrix of weights, the GRM becomes Weights can be computed with the SNP model (e.g., Gianola ), the GBLUP model (Zhang , 2015; Sun ), or ssGBLUP (Wang ).

Discussion

Derivations using SNP effects

Because adjacent SNP markers carry limited information due to linkage disequilibrium, the required number of core individuals in the APY may be equal to some number of “effective independent” SNP markers, as discussed; such a number would be almost independent of the actual number of markers if the number of actual markers is large enough. In Holsteins, predictions using SNP-BLUP were considerably more accurate with 40,000 than with 10,000 SNPs (VanRaden ). However, improvements in accuracy using SNP-BLUP with higher-density SNP chips (VanRaden ) or sequencing (Druet ) are very small. Predictions with APY for the national evaluation of ∼7 million Holsteins were accurate when the number of core animals was >10,000 (Figure 1), with little improvement beyond that number (Fragomeni ). This suggests the effective number of SNPs for Holstein cattle to be ∼10,000. For GWAS, Li presented methods for estimating the effective number of independent markers based on eigenvalues of a SNP correlation matrix. For a human HapMap CEU population, they estimated the number of independent SNP markers at <620,000. The effective population size for that population was estimated at ∼3100 (Tenesa ). Assuming that the number of such markers is proportional to an effective population size, a simple extrapolation for Holsteins (effective population size ∼100) leads to <20,000 effective independent SNP markers. Pintus found that ∼15,000 eigenvalues extracted from a matrix similar to the correlation matrix based on 40,000 SNPs explained 99% of variance for various traits of Holstein bulls. Marginally higher realized accuracies with higher SNP density could partially be due to lower sampling errors in the GRM (VanRaden 2008; Goddard ).

There is a question of whether the number of core animals is trait dependent. In this study, APY was derived as an extension of the theory for the inverse of the numerator relationship matrix, which is not trait specific. On the other hand, when all QTL are SNPs and identified, the number of effective SNPs is equal to the number of QTL.

ICS

The main advantage of derivations with ICS is linking the number of core animals to an effective population size. Many definitions exist for the average number of ICS (Me) (Daetwyler ), with the upper bound Me = 4N, where Ne is effective population size and L is the length of the genome in morgans (Stam 1980). Assuming that most commercial populations of animals have Ne ≤ 100 and L ≤ 30, Me ≤ 12,000.

The concept of ICS is abstract, with segments not directly identifiable. Because ICS are associated with linkage-disequilibrium blocks (Cuppen 2005), the number of ICS and the number of effective independent SNPs may be similar under the polygenic model. Estimates of the number of ICS in commercial livestock populations vary greatly when estimated from realized accuracies (Brard and Ricard 2015), probably because those accuracies are dependent on selection intensity, with smaller accuracies under stronger selection (Bijma 2012; Lourenco ). Another reason could be an implicit assumption of a constant size of each chromosome segment. In fact, Stam (1980) gave formulas for a distribution of segment size. Assuming that the proportion of variance explained is approximately proportional to segment size, a relatively small number of the largest segments may explain a large fraction of the variance while the remaining majority of segments may explain a small fraction of the variance. This is analogous to eigenvalue distribution. In a study by Fragomeni , the correlations of GEBVs obtained with regular and APY G−1 were very high, 0.94 with only 2000 core animals, increasing to 0.99 with 15,000 animals. In a study by Lourenco , 84% of accuracy gains due to genomic information were obtained using 4000 animals, with an increase to 97% with 8000 animals and 100% with 33,000 animals. As the reason for both a small number of ICS and reduced dimensionality of SNP information is linkage disequilibrium, the number of the largest effective SNPs and the number of largest ICS explaining the same amount of variance may be similar.

If the number of effective SNPs and the number of ICS are similar, these numbers can be estimated from eigenvalue decomposition of the GRM, as presented. This requires availability of a genotyped population with the size a few times larger than the number of ICS. For very large populations, computations of eigenvalues may be a limiting factor.

A relationship between Ne and Me allows one to determine applicability of the APY for different populations. In general, the APY can lead to computational savings when the number of genotyped individuals is large and at least twice the number of ICS. Assuming (from Holsteins) Ne ≈ 100, the APY is useful for ≥20,000 individuals. Extrapolating, with Ne ≈ 3000 (a smaller estimate in some human populations), the APY would be useful for ≥600,000 individuals. Therefore, the APY appears to be more useful for animal populations with small Ne.

Determining the number of core individuals

Too few core individuals would reduce accuracy of the inverse, whereas too many would increase cost and possibly decrease numerical stability because of unnecessary computations. The optimal number of core individuals could be defined such that increasing that number brings no notable improvement in accuracy of GEBVs calculated using APY This number of individuals can be determined by performing many analyses with different numbers of core individuals and comparing GEBVs or realized accuracies. Alternatively, the number of effective SNPs can be calculated from eigenvalue decomposition of the GRM as the number of the largest eigenvalues that explain 98–99% of the variation in the GRM, assuming that the remaining variation is due to sampling noise. The choice of the number of core animals is not critical.

Choice of animals for recursion

Based on the provided theory, the choice of animals for recursion as core animals is not critical as long as appropriate matrices are full rank, and this excludes only multiple copies of clones. In practice, the choice may have some impact. The APY is sparse, with the location of dense blocks dependent on the definition of core animals. When such an inverse is used in mixed-model equations solved iteratively, the convergence rate may vary with the choice of core animals. In Holsteins, the best convergence rate was found with core animals selected randomly whereas the slowest was with cows as core animals (Fragomeni ). Another factor is quality of genotypes. In commercial populations, animals may be genotyped with SNP chips of different density followed by, sometimes multiple, imputation to a standard SNP density. In genetic evaluation including ∼570,000 genotyped animals, use of popular sires as core resulted in slightly greater accuracy at the same size of recursion than random choices (Masuda et al. 2016). Therefore, the selection process may include the quality of genotypes. In particular, more “valuable” animals, e.g., popular sires, are more likely to have more accurate genotypes.

Which inverse is more accurate?

A priori is not clear which of regular and APY is superior. Masuda et al. (2016) looked at realized accuracies of Holstein bulls and found that those obtained with APY were marginally (<0.01) greater. If the proposed theory for the APY is applicable, the APY is accurate if the number of core individuals reaches the number of ICS or effective number of independent SNPs. Therefore, for large populations most computations with regular G−1 are redundant and may lead to numerical problems. As argued before, the real rank of GRM (disregarding very small eigenvalues) is likely equal to the number of ICS (or independent SNPs). The GRM with a larger number of individuals is singular and standard procedures to invert it include blending, a weighted mean of the original GRM with a NRM (VanRaden 2008; Aguilar ); the primary purpose of blending is numerical stability of inversion as GEBVs obtained with blending 1–10% of the NRM were nearly identical (Misztal ). Another reason for higher accuracy of APY could be less influence from sampling errors in the GRM, as the APY does not use the block of GRMs due to noncore individuals (except for diagonals); the standard error of an element of the GRM assuming 0.5 allele frequencies is where m is the number of SNP markers (VanRaden 2008).

Concluding remarks

The presented theory explains why and when recursions on a small subset of animals lead to an efficient computation of inverse of the GRM. Such an inverse is accurate when the subset is as large as the number of ICS or the rank of the GRM; for Ne = 100 the number of ICS is ∼10,000. When the number of genotyped individuals is much larger than the number of ICS, the APY inverse is sparse and facilitates genomic evaluation, parameter estimation, and GWAS at greatly reduced cost and potentially higher accuracy than a regular inverse.