Approximated prediction of genomic selection accuracy when reference and candidate populations are related
- Jean-Michel Elsen^{1, 2}Email author
Received: 23 July 2015
Accepted: 8 January 2016
Published: 3 March 2016
Abstract
Background
Genomic selection is still to be evaluated and optimized in many species. Mathematical modeling of selection schemes prior to their implementation is a classical and useful tool for that purpose. These models include formalization of a number of entities including the precision of the estimated breeding value. To model genomic selection schemes, equations that predict this reliability as a function of factors such as the size of the reference population, its diversity, its genetic distance from the group of selection candidates genotyped, number of markers and strength of linkage disequilibrium are needed. The present paper aims at exploring new approximations of this reliability.
Results
Two alternative approximations are proposed for the estimation of the reliability of genomic estimated breeding values (GEBV) in the case of non-independence between candidate and reference populations. Both were derived from the Taylor series heuristic approach suggested by Goddard in 2009. A numerical exploration of their properties showed that the series were not equivalent in terms of convergence to the exact reliability, that the approximations may overestimate the precision of GEBV and that they converged towards their theoretical expectations. Formulae derived for these approximations were simple to handle in the case of independent markers. A few parameters that describe the markers’ genotypic variability (allele frequencies, linkage disequilibrium) can be estimated from genomic data corresponding to the population of interest or after making assumptions about their distribution. When markers are not in linkage equilibrium, replacing the real number of markers and QTL by the “effective number of independent loci”, as proposed earlier is a practical solution. In this paper, we considered an alternative, i.e. an “equivalent number of independent loci” which would give a GEBV reliability for unrelated individuals by considering a sub-set of independent markers that is identical to the reliability obtained by considering the full set of markers.
Conclusions
This paper is a further step towards the development of deterministic models that describe breeding plans based on the use of genomic information. Such deterministic models carry low computational burden, which allows design optimization through intensive numerical exploration.
Background
The effectiveness of genomic selection comes from the possibility of predicting breeding values on un-phenotyped and young animals [1]. Genomic selection promised and proved to be extremely efficient and beneficial for dairy cattle (e.g. [2–7]), but debate continues for other species and production sectors (e.g. [8–12]). A key criterion to decide whether or not selection schemes (also referred to here as breeding plans) should include genomic information is the reliability of the genomic predictor. It was clearly shown that this reliability depends on the structure of the reference population and on the characteristics of the marker set used. The size of this reference population, its diversity, the genetic distance between the reference and the group of selection candidates genotyped, the number of markers, and the degree or strength of the linkage disequilibrium are the main factors that influence this reliability [13–23].
An extensive literature exists on the mathematical modeling of selection schemes prior to their implementation, in order, for instance, to optimize their design, or to evaluate the usefulness of new technologies such as embryo transfer, sperm selection, DNA markers and others (e.g. [24] for a review). These models account for factors such as selection intensities and maintenance or loss of genetic variability. Among these parameters, the precision of breeding value estimates is central. To model genomic selection schemes, equations that predict this reliability as a function of the factors cited above are needed (e.g. [6, 25, 26]).
The quantitative influence of these factors (size of the reference population, its diversity, etc.) was assessed by simulation studies [18–21, 27, 28]. An equation that predicts the reliability of genomic evaluation in the very simple situation of independent quantitative trait loci (QTL), that are perfectly marked by single nucleotide polymorphisms (SNPs) and populations (reference and candidates) of unrelated individuals was derived [13]. This approach was extended to the case when only a part of the genetic variability is imperfectly marked by SNPs [16, 19], and the situation of non-independence between reference and candidate populations was explored [17]. It was demonstrated that genomic information captures historical linkage disequilibrium, short-term linkage between QTL and markers and additive relationships between reference and candidate individuals, the equation of the reliability accounting for these three phenomena being derived in a very simple case of one QTL marked by a single SNP [22].
A Taylor expansion of a matrix inverse involved in the reliability formula was suggested [18], which led to the algebraic development of an approximation. This approximation seems to work well in the simple situation but lacks generality. In this paper, an alternative approximation is proposed, opening a way to include non-independence between reference and candidate populations, and between markers.
After a formalization of the genomic selection context, the principles that underlie these approximations are presented and their properties are compared by using a simple example. Then, the new approximation is derived when reference and candidate animals are related. This is illustrated by some numerical examples. Finally, the extension to the linkage disequilibrium situation is described.
Methods
General framework
Although the prediction equations derived below were based on a number of simplifying assumptions, it is important to first draw a complete description of the biological framework, as a basis to subsequently simplify the discussion.
The SNP effects are estimated in a reference population, P_{r}, comprising n _{ r } individuals. The genomic estimated breeding values (GEBV) are calculated for a population of candidates for selection and used in breeding, P_{c}, comprising n _{ c } individuals.
Let \( {\mathbf{\mathcal{P}}} = \left( {{\mathbf{\mathcal{P}}}_{r} ,{\mathbf{\mathcal{P}}}_{c} } \right) \) the population structure (including pedigree relationships between individuals and marker allele frequencies, but not including genotypes and phenotypes).
Individuals are characterized by their genotypes at n _{ M } markers (observed) and at n _{ Q } QTL (unknown). Alleles will be noted A _{ m } and B _{ m } for the marker m and A _{ q } and B _{ q } for the QTL q. Let a _{ tim } ∊ {0, 1, 2} and a _{ tiq } ∊ {0, 1, 2} be the numbers of B _{ m } (and respectively, B _{ q }) alleles that an individual i from population P_{t} (P_{r} or P_{c}) carries at marker m (respectively, QTL q). Let p _{ tm } and p _{ tq } be the frequencies of alleles B _{ m } and B _{ q } in P_{t}.
Genotypic values will be assigned to the different markers and QTL genotypes. Following [18], genotypes will be coded as x _{ tim } = a _{ tim } − 2p _{ tm } and w _{ tiq } = a _{ tiq } − 2p _{ tq }. Different codifications can be proposed [15]. In particular, as described for instance in [29], genotypic values may be standardized, i.e. x _{ tim } = (a _{ tim } − 2p _{ tm })/σ _{ tm } and w _{ tiq } = (a _{ tiq } −2p _{ tq })/σ _{ tq }, with variances σ _{ tm } ^{2} = 2p _{ tm }(1 − p _{ tm }) and σ _{ tq } ^{2} = 2p _{ tq }(1 − p _{ tq }). Most of the following developments are given with the first codification here, and the results with the second codification are described in a specific section.
These genotypic values are assembled in matrices X (dim (X) = (n _{ r } + n _{ c }) × n _{ M }) and W (dim (W) = (n _{ r } + n _{ c }) × n _{ Q }). Sub-matrices corresponding to sub-populations will be noted in the following way: X ^{′} = (X _{ r } ^{′} , X _{ c } ^{′} ) and W ^{′} = (W _{ r } ^{′} , W _{ c } ^{′} ).
The genetic model assumes additivity of QTL effects. The additive genetic value of an individual is described as \( g_{ti} = \sum\nolimits_{q = 1}^{{n_{Q} }} {w_{tiq} \alpha_{q} } \) and, in general, \( {\mathbf{ g = W\alpha }} \). The phenotypic values when observed are y = g + ɛ.
A statistical model describes the performances in the reference population as random variables for which the expectations are linear combinations of SNP effects: \( y_{ri} = \sum\nolimits_{m = 1}^{{n_{S} }} {X_{rim} \beta_{m} + e_{ri} } \) and, in general, y = Xβ + e.
In these models, the SNP (or QTL) effects may be considered as fixed, or random. Since the number of SNPs is much bigger than the number of individuals (n _{ M } ≫ n _{ r }) the second solution is generally chosen in the statistical model (but not always see [1, 13]).
In the random model, a distribution \( {\mathcal{L}}\left( {{\varvec{\uptheta}}_{{\upbeta }} ,{\mathbf{V}}_{{\upbeta }} } \right) \) (respectively \( {\mathcal{L}}\left( {{\varvec{\uptheta}}_{\upalpha} ,{\mathbf{V}}_{\upalpha} } \right) \)) of the SNP (respectively QTL) effects is assumed, with θ _{ β } (respectively, θ _{α}) being the vector of expectations and V _{β} (respectively, V _{α}) being the matrix of variances. For a full description of the variability, the V _{β} and V _{α} matrices are each subdivided into four blocks corresponding to the reference and candidate populations and their covariances. Covariances between the α and β vectors have also to be considered. Most generally, the SNP (QTL) effects are supposed i.i.d. giving V _{β} = Iσ _{β} ^{2} (V _{α} = Iσ _{α} ^{2} ). The interpretation of these QTL effects is nicely debated in Gianola et al. [30]. In the frequentist view, we simply have to imagine that QTL effects are randomly sampled from a distribution with a σ _{α} ^{2} variance. In the Bayesian context, the prior variability of the SNP effects was most generally described as heteroskedastic or even coming from mixtures of SNPs with or without an effect on the trait.
The expectations \( {\varvec{\uptheta}}_{{\upbeta }} \left( {{\varvec{\uptheta}}_{{\upalpha }} } \right) \) are generally assumed equal to zero, but when information about population history is available (in particular, when we know it is a mixed population), non-zero values should be considered.
The vector q = Xβ is a quantity similar but not equal to the genetic value g. Its element q _{ ti } is the molecular score of individual i in population t. This vector may be segmented in two parts: \( {\mathbf{q^{\prime}}} = \left( {{\mathbf{q^{\prime}_{{{r}}}}} ,{\mathbf{q^{\prime}_{{\mathbf{c}}} }}} \right) \).
Since the variances may be defined within a population, we have \( v\left( {{\mathbf{q}}_{{\mathbf{r}}} |{\mathbf{X}}} \right) = {\mathbf{X}}_{{\mathbf{r}}} {\mathbf{V}}_{{{\bf\upbeta r}}} {\mathbf{X}}^{\prime}_{{\mathbf{r}}} \), and v(q _{ c }|X) = X _{ c } V _{βc} X _{ c } ^{′} . The residual variance is v(e) = I σ _{ e } ^{2} . Assuming that the distribution of marker effects is centered (\( {\varvec{\uptheta}}_{{\varvec{\upbeta}}} = {\bf 0} \)) and i.i.d. (\( {\mathbf{V}}_{{{\upbeta r}}} = {\mathbf{I}}\sigma_{{{\upbeta r}}}^{2} \, {\text{and}} \,{\mathbf{V}}_{{{\upbeta c}}} = {\mathbf{I}}\sigma_{{{\upbeta c}}}^{2} \)), and extending Gianola et al. [30], we have \( v\left( {q_{ri} } \right) = {\text{E}}_{{\mathbf{X}}} \left[ {v\left( {q_{\text{ri}} |{\mathbf{X}}} \right)} \right] = \sigma_{\beta r}^{2} \sum\nolimits_{m} {2p_{mr} \left( {1 - p_{mr} } \right)} = \sigma_{\beta r}^{2} \sum\nolimits_{m} {\sigma_{mr}^{2} = \sigma_{\beta r}^{2} \tau_{r} } \) in the reference population, and v(q _{ ci }) = σ _{βc} ^{2} τ _{ c } in the candidate population. Assuming that the distribution of the marker effects and genotypes are the same in P_{r} and P_{c}, i.e. \( p_{rm} = p_{cm} = p_{m} , p_{rq} = p_{cq} = p_{q} \), thus τ _{ r } = τ _{ C } = τ and σ _{ βr } ^{2} = σ _{ βc } ^{2} = σ _{ β } ^{2} , we define σ _{ q } ^{2} = τσ _{β} ^{2} . Thus, \( v\left( {{\mathbf{q}}|{\mathbf{X}}} \right) = \frac{1}{\tau }{\mathbf{XX^{\prime}}}\sigma_{{q}}^{2} \). These equations hold even if the markers are in linkage disequilibrium (LD) as shown in Eq. A2 from Gianola et al. [30].
We note σ ^{2} as the total phenotypic variance, i.e. σ ^{2} = σ _{ q } ^{2} + σ _{ e } ^{2} , and ν ^{2} as the proportion of this variance explained by the molecular score \( \left( {\nu^{2} = \frac{{\sigma_{q}^{2} }}{{\sigma^{2} }}} \right) \). The ratio \( \frac{{\sigma_{q}^{2} }}{{\sigma_{e}^{2} }} \) will be noted γ.
The SNP effects β may be estimated in different ways. The genomic best linear unbiased prediction (BLUP) will only be considered here, with \( {\hat{\varvec{\upbeta} }} = {\text{cov}}\left( {{\varvec{\upbeta}},{\mathbf{y}}} \right){\text{var}}\left( {\mathbf{y}} \right)^{ - 1} {\mathbf{y}} \). Classically, this equation becomes \( {\hat{\varvec{\upbeta }}} = \sigma_{{{\upbeta }}}^{2} {\mathbf{X}}^{\prime}_{{\mathbf{r}}} \left[ {\sigma_{{{\upbeta }}}^{2} \left( {{\mathbf{X}}_{{\mathbf{r}}} {\mathbf{X}}^{\prime}_{{\mathbf{r}}} + {\mathbf{I}}{{\uplambda }}_{{{\upbeta }}} } \right)} \right]^{ - 1} {\mathbf{y}} = \left( {{\mathbf{X}}^{\prime}_{{\mathbf{r}}} {\mathbf{X}}_{{\mathbf{r}}} + {\mathbf{I}}{{\uplambda }}_{{{\upbeta }}} } \right)^{ - 1} {\mathbf{X}}^{\prime}_{{\mathbf{r}} }{\mathbf{y}} \) with λ_{β} = σ _{e} ^{2} /σ _{β} ^{2} . The linear combination \( \hat{\varvec{q}}_{\varvec{c}} = {\mathbf{X}}_{\varvec{c}} {\hat{\varvec{\upbeta}}} \) is the GBLUP vector for candidates in P_{c}. It must be emphasized that these estimations and predictions are conditional on the genotypic structures defined by X (X _{ r } and \( {\mathbf{X}}_{\varvec{c}} \)).
Given X, the reliability of the GBLUP is \( r^{2} \left( {g_{ci} ,\hat{q}_{ci} |{\mathbf{X}}} \right) = \frac{{cov^{2} \left( {g_{ci} ,\hat{q}_{ci} |{\mathbf{X}}} \right)}}{{v\left( {g_{ci} |{\mathbf{X}}} \right)v\left( {\hat{q}_{ci} |{\mathbf{X}}} \right)}} \).
In [16], the reliability is described (Eq. 6 in [16]) as \( r\left( {g_{ci} ,\hat{q}_{ci} } \right) = r\left( {g_{ci} ,q_{ci} } \right) \times r\left( {q_{ci} ,\hat{q}_{ci} } \right) \), by ignoring the conditioning on X. In Goddard et al. [18], the reliability is described as \( r_{{g_{ci} ,\hat{q}_{ci} }}^{2} = \frac{{v\left( {\hat{q}_{ci} } \right)}}{{v\left( {g_{ci} } \right)}} = \frac{{v\left( {q_{ci} } \right)}}{{v\left( {g_{ci} } \right)}}\frac{{v\left( {\hat{q}_{ci} } \right)}}{{v\left( {q_{ci} } \right)}} \). In this formulation, \( \frac{{v\left( {q_{ci} } \right)}}{{v\left( {g_{ci} } \right)}} \) is the proportion of the genetic variance explained by the markers and \( \frac{{v\left( {\hat{q}_{ci} } \right)}}{{v\left( {q_{ci} } \right)}} \) is the accuracy of estimated marker effects. This is similar to the \( {\text{qr}}_{{\hat{Q}}} \) reported by Dekkers et al. [25]. All these reliability formulae are approximations since \( cov^{2} \left( {g_{ci} ,\hat{q}_{ci} } \right) = cov^{2} \left( {\sum w_{ciq} \alpha_{q} ,\sum x_{cis} \hat{\beta }_{s} } \right) \ne v\left( {\hat{q}_{ci} } \right) = v\left( {\sum x_{cis} \hat{\beta }_{s} } \right) \), in general.
Situation analyzed in this paper
In the following, ignoring the difficulty that was mentioned above, we will assume \( r^{2} \left( {q_{\text{ci}} ,\hat{q}_{ci} |{\mathbf{X}}} \right) = \frac{{cov^{2} \left( {q_{ci} ,\hat{q}_{ci} |{\mathbf{X}}} \right)}}{{v\left( {q_{ci} |{\mathbf{X}}} \right)v\left( {\hat{q}_{ci} |{\mathbf{X}}} \right)}} = \frac{{v\left( {\hat{q}_{ci} |{\mathbf{X}}} \right)}}{{v\left( {q_{ci} |{\mathbf{X}} } \right)}} \). We are interested in a single candidate in \( \varvec{P}_{\varvec{c}} \) with a x _{ c } vector of marker genotypes.
Formulae were simplified in two ways. (1) the i index of the candidate was omitted in the following developments: the genetic value of the candidate is noted q _{ c }, estimated by \( \hat{q}_{c} = cov\left( {q_{c} ,{\mathbf{y}}} \right)v\left( {\mathbf{y}} \right)^{ - 1} {\mathbf{y}} \), and its precision is \( r^{2} \left( {q_{\text{c}} ,\hat{q}_{c} |{\mathbf{X}}} \right) = \frac{{{\text{v}}(\hat{q}_{c} |{\mathbf{X}})}}{{{\text{v}}\left( {q_{\text{c}} |{\mathbf{X}}} \right)}} \), with \( v\left( {q_{\text{c}} |{\mathbf{X}}} \right) = \sigma_{{{\upbeta }}}^{2} {\mathbf{x}}_{{\mathbf{c}}} {\mathbf{x}}^{\prime}_{{\mathbf{c}}} \) and \( v\left( {\hat{q}_{c} |{\mathbf{X}}} \right) = \sigma_{{{\upbeta }}}^{2} {\mathbf{x}}_{{\mathbf{c}}} {\mathbf{X}}^{\prime}_{{\mathbf{r}}} \left( {{\mathbf{X}}_{{\mathbf{r}}} {\mathbf{X}}^{\prime}_{{\mathbf{r}}} + {\mathbf{I}}{{\uplambda }}_{{{\upbeta }}} } \right)^{ - 1} {\mathbf{X}}_{{\mathbf{r}}} {\mathbf{x}}^{\prime}_{{\mathbf{c}}} \) (where \( {\mathbf{x}}_{{\mathbf{c}}} \) is a row vector); (2) the r indices of reference individuals were most often omitted, which resulted in y _{ i } for their phenotypes and q _{ i } for their molecular scores.
In fact, our objective was to estimate the expectation of this precision across the variation domain of X _{ r } and x _{ c } given the pedigree structure \( \left( {{\mathbf{\mathcal{P}}}_{\varvec{r}} ,{\mathbf{\mathcal{P}}}_{\varvec{c}} } \right):\;{\text{E}}_{{\mathbf{X}}} \left[ {r^{2} \left( {q_{\text{c}} ,\hat{q}_{c} |{\mathbf{X}}} \right)|{\mathbf{\mathcal{P}}}} \right] \). It will be noted \( E\left[ {r_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right] \).
The following approximation was made: \( {\text{E}}\left[ {r_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right] = \frac{{{\text{E}}_{{\mathbf{X}}} \left[ {{\text{v}}(\hat{q}_{c} |{\mathbf{X}})} \right]}}{{{\text{E}}_{{\mathbf{X}}} \left[ {{\text{v}}\left( {q_{\text{c}} |{\mathbf{X}}} \right)} \right]}} = \frac{{{\text{E}}\left[ {{\text{v}}\left( {\hat{q}_{c} } \right)} \right]}}{{{\text{E}}\left[ {{\text{v}}\left( {q_{\text{c}} } \right)} \right]}} \).
Let A be the pedigree relationship matrix between individuals in \( {\mathbf{\mathcal{P}}} \). Its blocks are \( {\mathbf{A}} = \left( {\begin{array}{*{20}c} {{\text{a}}_{\text{cc}} } & {{\mathbf{A}}_{{{\mathbf{cr}}}} } \\ {{\mathbf{A}}_{{{\mathbf{rc}}}} } & {{\mathbf{A}}_{{{\mathbf{rr}}}} } \\ \end{array} } \right) \). Let \( {\mathbf{G}}^{ *} = {\mathbf{XX^{\prime}}} = \left( {\begin{array}{*{20}c} {{\mathbf{x}}_{{\mathbf{c}}} {\mathbf{x}}^{\prime}} & {{\mathbf{x}}_{{\mathbf{c}}} {\mathbf{X}}^{\prime}_{{\mathbf{r}}} } \\ {{\mathbf{X}}_{{\mathbf{r}}} {\mathbf{x}}^{\prime}_{{\mathbf{c}}} } & {{\mathbf{X}}_{{\mathbf{r}}} {\mathbf{X}}^{\prime}_{{\mathbf{r}}} } \\ \end{array} } \right) \), which results in \( {\mathbf{V}} = \frac{1}{\tau }{\mathbf{G}}^{ *} \sigma_{q}^{2} + {\mathbf{I}}\sigma_{e}^{2} \). It must be noted that the σ _{ e } ^{2} term in the diagonal of the V submatrix corresponding to the candidate population is artificial since candidates are not phenotyped.
We have E[G*] = A τ. The limits of this equality will be discussed below. As indicated above, the denominator of the expected reliability \( {\text{E}}_{{\mathbf{X}}} \left[ {v\left( {q_{\text{c}} |{\mathbf{X}}} \right)} \right] \), is τσ _{β} ^{2} = σ _{ q } ^{2} . Approximating \( {\text{E}}\left[ {v\left( {\hat{q}_{c} } \right)} \right] \) by E[cov(q _{ c }, y)]E[v(y)]^{−1} E[cov(y, q _{ c })] is useless because it makes an oversimplification of the relationships between the reference and the candidate population: it considers separately the marginal distributions of x _{ c } X _{ r } ^{′} and (X _{ r } X _{ r } ^{′} + Iλ_{β})^{−1}, while these random matrices are correlated. Estimating directly E[cov(q _{ c }, y)v(y)^{−1} cov(y, q _{ c })] seems impossible in the general case. The approach of Goddard et al. [18] avoids this difficulty, i.e. the variance \( v\left( {\hat{q}_{c} |{\mathbf{X}}} \right) = \sigma_{{{\upbeta }}}^{2} {\mathbf{x}}_{{\mathbf{c}}} {\mathbf{x}}^{\prime}_{{\mathbf{c}}} + \sigma_{e}^{2} - \frac{1}{{\left\{ {{\mathbf{V}}^{ - 1} } \right\}_{{\varvec{cc}}} }} \), and V ^{−1} is approximated by a second degree Taylor expansion (\( {\mathbf{V}}^{ - 1} \sim {\varvec{\Lambda}}\left( {\mathbf{X}} \right) \)), giving \( v\left( {\hat{q}_{c} |{\mathbf{X}}} \right)\sim \sigma_{{{\upbeta }}}^{2} {\mathbf{x}}_{{\mathbf{c}}} {\mathbf{x}}^{\prime}_{{\mathbf{c}}} + \sigma_{e}^{2} - \frac{1}{{{\varvec{\Lambda}}_{\text{cc}} \left( {{\mathbf{x}}_{{\mathbf{c}}} ,{\mathbf{X}}_{{\mathbf{r}}} } \right)}} \).
Alternative approximations of the reliability
Extension of Goddard’s formula
In their “heuristic approximation for V ^{*−1}”, Goddard et al. [18] considered the situation where unrelated individuals are included in the reference and candidate populations, that is E[G ^{*}] = I τ and \( {\mathbf{G}}^{ *} = {\mathbf{I}}\tau + {\mathbf{E}} \), with E, a “noise” matrix centered on the null matrix \( {\bf 0}. \) A direct extension of their development would be the following. The matrix \( {\mathbf{V}} = \frac{1}{\tau }{\mathbf{G}}^{ *} \sigma_{q}^{2} + {\mathbf{I}}\sigma_{e}^{2} \) can be written as:
V = σ _{ e } ^{2} (I + A γ)[I + D γ],
with \( {\mathbf{D}} = \left( {{\mathbf{I}} + {\mathbf{A}}\gamma } \right)^{ - 1} \left( {\frac{1}{\tau }{\mathbf{G}}^{ *} - {\mathbf{A}}} \right) = {\mathbf{T}}\left( {\frac{1}{\tau }{\mathbf{G}}^{ *} - {\mathbf{A}}} \right) \),
and \( \gamma = \frac{{\sigma_{q}^{2} }}{{\sigma_{e}^{2} }} \). Thus, \( {\mathbf{V}}^{ - 1} = \frac{1}{{\sigma_{e}^{2} }}\left[ {{\mathbf{I}} + {\mathbf{D}}\gamma } \right]^{ - 1} {\mathbf{T}} \). The inverse matrix [I + D γ]^{−1} will be approximated using a Taylor series. It must be emphasized that the Taylor series I − D γ + (D γ)^{2} − (D γ)^{3} + ··· converges towards [I + D γ]^{−1} only if the highest Eigen value of D γ is smaller than 1, i.e. if \( \left( {{\mathbf{D}}\gamma } \right)^{\text{t}} \to {\bf 0} \) when t → ∞.
The second order approximation of V ^{−1} is equal to \( \frac{1}{{\sigma_{e}^{2} }}\left( {{\mathbf{I}} - {\mathbf{D}}\gamma + {\mathbf{D}}^{2} \gamma^{2} } \right){\mathbf{T}} \). As E[D] = 0 and \( {\text{E}}\left[ {{\mathbf{D}}^{2} } \right] = {\mathbf{T}}\left( {\frac{1}{{\tau^{2} }}{\text{E}}\left[ {{\mathbf{G}}^{ *} {\mathbf{TG}}^{ *} } \right] - {\mathbf{ATA}} } \right), \) its expectation \( {\text{E}}\left[ {\varvec{\Lambda}} \right] = \frac{1}{{\sigma_{e}^{2} }}\left( {{\mathbf{I}} - {\text{E}}\left[ {\mathbf{D}} \right]\gamma + {\text{E}}\left[ {{\mathbf{D}}^{2} } \right]\gamma^{2} } \right){\mathbf{T}} \) i.e. \( {\text{E}}\left[ {\varvec{\Lambda}} \right] = \frac{1}{{\sigma_{e}^{2} }}\left( {{\mathbf{I}} - \gamma^{2} {\mathbf{TATA}} + \frac{{\gamma^{2} }}{{\tau^{2} }}{\text{E}}\left[ {{\mathbf{TG}}^{ *} {\mathbf{TG}}^{ *} } \right]} \right){\mathbf{T}} \).
A difficulty with this approximation comes from the T term. As an example, consider a reference population composed of n _{ r } half-sibs of the candidate, \( {\mathbf{T}} = {{\upxi }}{\mathbf{I}} + {{\uppsi }}{\mathbf{J}} \) with \( {{\upxi }} = \frac{4}{4 + 3\gamma } \). As \( {\mathbf{T}}^{\varvec{t}} = {{\upxi }}^{\varvec{t}} {\mathbf{I}} + \left[ {n_{r}^{t} {{\upxi }}^{\varvec{t}} + \cdots } \right]{\mathbf{J}} \), the J coefficient will tend to ∞ as soon as \( n_{r} {{\upxi }} = \frac{{4n_{r} }}{4 + 3\gamma } > 1 \), a very realistic situation. Thus, the convergence of the Taylor series will be a balance between the increase of \( {\mathbf{T}}^{\varvec{t}} \) and decrease of \( \left[ {{\mathbf{D}}\gamma } \right]^{\varvec{t}} \).
Another approximation of the reliability
Principle
Using the classical matrix inversion lemma, the variance \( v\left( {\hat{q}_{{\text{c}}} |{\mathbf{x}}_{{\mathbf{c}}} ,{\mathbf{X}}_{{\mathbf{r}}} } \right) = \sigma _{\upbeta }^{2} {\mathbf{x}}_{{\mathbf{c}}} {\mathbf{X}}^{\prime}_{{\mathbf{r}}} \left( {{\mathbf{X}}_{{\mathbf{r}}} {\mathbf{X}}^{\prime}_{{\mathbf{r}}} + {\mathbf{I}}\uplambda _{\upbeta } } \right)^{{ - 1}} {\mathbf{X}}_{{\mathbf{r}}} {\mathbf{x}}^{\prime}_{{\text{c}}} \) may also be defined as \( {\boldsymbol{v}}\left( {\hat{q}_{c} |{\mathbf{x}}_{{\mathbf{c}}} ,{\mathbf{X}}_{{\mathbf{r}}} } \right) = \sigma_{{{\upbeta }}}^{2} {\mathbf{x}}_{{\mathbf{c}}} {\mathbf{x}}^{\prime}_{{\mathbf{c}}} - \sigma_{\text{e}}^{2} {\mathbf{x}}_{{\mathbf{c}}} \left( {{\mathbf{X}}^{\prime}_{{\mathbf{r}}} {\mathbf{X}}_{{\mathbf{r}}} + {\mathbf{I}}{{\uplambda }}_{{{\upbeta }}} } \right)^{ - 1} {\mathbf{x}}^{\prime}_{{\mathbf{c}}} \).
\( {\mathbf{X}}^{\prime}_{{\mathbf{r}}} {\mathbf{X}}_{{\mathbf{r}}} \) is a very large matrix \( \left( {n_{M} \times n_{M} } \right) \) that describes the LD between markers: its elements tend to be smaller when they are more distant from the diagonal.
Elements of E[X _{ r } ^{ ′ } X _{ r }] are the following: \( {\text{E}}\left[ {{\mathbf{X}}^{\prime}_{{\mathbf{r}}} {\mathbf{X}}_{{\mathbf{r}}} } \right]_{\text{ml}} = {\text{E}}\left[ {\mathop {\sum\nolimits_{\text{i}} {\left( {a_{im} - 2p_{m} } \right)\left( {a_{il} - 2p_{l} } \right)} }\limits_{{}} } \right] = 2n_{r} \Updelta_{ml}, \) with \( \Updelta_{ml} \) the LD between loci m and l.
E[X _{ r } ^{′} X _{ r }]_{mm} = E[ ∑ _{i}(a _{ im } − 2p _{ m })^{2}] = n _{ r } σ _{m} ^{2} . Let \( {\mathbf{C}} = {\mathbf{I}}{{\uplambda }}_{{{\upbeta }}} + n_{r} {\text{diag}}\left[ {\sigma_{1}^{2} , \ldots ,\sigma_{{n_{M} }}^{2} } \right] \), the X _{ r } ^{′} X _{ r } + Iλ_{β} matrix may be written as:
\( {\mathbf{X_r^\prime}} {\mathbf{X}}_{{\mathbf{r}}} + {\mathbf{I}}{{\uplambda }}_{{{\upbeta }}} = \left[ {\left( {{\mathbf{X}}^{\prime}_{{\mathbf{r}}} {\mathbf{X}}_{{\mathbf{r}}} - n_{r} {\text{diag}}\left[ {\sigma_{1}^{2} , \ldots ,\sigma_{{n_{M} }}^{2} } \right]} \right){\mathbf{C}}^{ - 1} + {\mathbf{I}}} \right]{\mathbf{C}} \), which results in:
X _{ r } ^{′} X _{ r } + Iλ_{β} = [I + B]C.
The convergence of the Taylor series I − B + B ^{2} − B ^{3} + ··· to (I + B)^{−1} depends on the structure of the B matrix, which varies depending on the sample. However, we can examine the case of its expectation E[B].
E[B]_{ mm } = 0 and \( {\text{E}}\left[ {\mathbf{B}} \right]_{ml} = \frac{{2n_{r} \Updelta_{ml} }}{{{{\uplambda }}_{{{\upbeta }}} + n_{r} \sigma_{\text{l}}^{2} }} \). The ratio λ_{β} is proportional to the number of markers (\( {{\uplambda }}_{{{\upbeta }}} = n_{M} \frac{{\bar{\sigma }_{\text{m}}^{2} \sigma_{e}^{2} }}{{\sigma_{\text{q}}^{2} }}) \) and dominates the denominator when n _{ M } ≫ n _{ r }. The (m, l) term in E[B]^{2}, i.e. \( {\text{E}}\left[ {\mathbf{B}} \right]_{{\varvec{ml}}}^{2} = \sum\nolimits_{k} {\frac{{4n_{r}^{2} \Updelta_{mk} \Updelta_{kl} }}{{\left( {{{\uplambda }}_{{{\upbeta }}} + n_{r} \sigma_{\text{k}}^{2} } \right)\left( {{{\uplambda }}_{{{\upbeta }}} + n_{r} \sigma_{\text{l}}^{2} } \right)}}} \), is of order \( \frac{1}{{n_{M} }} \). Thus, we expect the Taylor series to converge to (I + E[B])^{−1}.
First order approximation
Results
Application in the case of independent markers
This situation either assumes low density marker information, or corresponds to the idea of an effective number of loci that was developed by Goddard [16, 31]. In the first case, the proportion of the genetic variance explained by the markers \( \frac{{v\left( {q_{ci} } \right)}}{{v\left( {g_{ci} } \right)}} \) is small and this quantity should be considered when estimating the genomic precision.
First approximation
Coefficients describing the genotypes’ distributions moments when using the relation \( {\text{E}}\left[ {X_{im}^{{d_{i} }} X_{jm}^{{d_{j} }} \cdots X_{Km}^{{d_{K} }} } \right] = p_{m} \left( {1 - p_{m} } \right)\alpha_{ij \cdots K}^{{d_{i} d_{j} \cdots d_{K} }} - \left[ {p_{m} \left( {1 - p_{m} } \right)} \right]^{2} \gamma_{ij \cdots K}^{{d_{i} d_{j} \cdots d_{K} }} \) from Additional file 1
E[X _{ im }] | \( \alpha_{i}^{1} = 0 \) \( \gamma_{i}^{1} = 0 \) |
E[X _{ im } ^{2} ] | \( \alpha_{i}^{2} = 2 \) \( \gamma_{i}^{2} = 0 \) |
E[X _{ im } ^{4} ] | \( \alpha_{i}^{4} = 2 \) \( \gamma_{i}^{4} = 0 \) |
E[X _{ im } X _{ jm }] | \( \alpha_{ij}^{11} = 4\delta_{ 1} + 2\left[ {\delta_{ 2} + \delta_{ 3} + \delta_{ 4} + \delta_{ 5} + \delta_{ 9} + \delta_{ 1 2} } \right] + \delta_{ 10} + \delta_{ 1 1} + \delta_{ 1 3} + \delta_{ 1 4} = 4a_{ij} \) \( \gamma_{ij}^{ 1 1} = 0 \) |
E[X _{ im } X _{ jm } ^{3} ] | \( \alpha_{ij}^{13} = 16\delta_{1} + 2\left( {\delta_{2} + \delta_{3} ) + 8(\delta_{4} + \delta_{5} } \right) + 2\left( {\delta_{9} + \delta_{12} } \right) + \delta_{10} + \delta_{11} + \delta_{13} + \delta_{14} \) \( \gamma_{ij}^{13} = 2 4\delta_{ 1} + 1 2\left( {\delta_{ 4} + \delta_{ 5} } \right) \) |
E[X _{ im } ^{2} X _{ jm } ^{2} ] | \( \alpha_{ij}^{22} = 1 6\delta_{ 1} + 4\left( {\delta_{ 2} + \delta_{ 3} + \delta_{ 4} + \delta_{ 5} } \right) + 2\left( {\delta_{ 9} + \delta_{ 1 2} } \right) + \delta_{ 10} + \delta_{ 1 1} + \delta_{ 1 3} + \delta_{ 1 4} \) \( \gamma_{ij}^{22} = 4 8\delta_{ 1} + 8\left( {\delta_{ 2} + \delta_{ 3} + \delta_{ 4} + \delta_{ 5} } \right) - 4\delta_{ 1 5} - 1 6\delta_{ 6} - 8\left( {\delta_{ 7} + \delta_{ 8} } \right) \) |
Moments of genotypes’ distributions depending on genotype codification
Expectations | Genotype codification | |
---|---|---|
x _{ tim } = a _{ tim } − 2p _{ tm } | x _{ tim } = (a _{ tim } − 2p _{ tm })/σ _{ tm } | |
E[X _{ im }] | 0 | 0 |
E[X _{ im } ^{2} ] | σ _{m} ^{2} | 1 |
E[X _{ im } ^{4} ] | σ _{m} ^{2} | 1/σ _{m} ^{2} |
E[X _{ im } X _{ jm }] | 2a _{ ij } σ _{m} ^{2} | 2a _{ ij } |
E[X _{ im } X _{ jm } ^{3} ] | \( \frac{1}{2}\alpha_{ij}^{13} \sigma_{\text{m}}^{2} - \frac{1}{4}\gamma_{ij}^{13} \sigma_{\text{m}}^{4} \) | \( \frac{1}{2}\alpha_{ij}^{13} /\sigma_{\text{m}}^{2} - \frac{1}{4}\gamma_{ij}^{13} \) |
E[X _{ im } ^{2} X _{ jm } ^{2} ] | \( \frac{1}{2}\alpha_{ij}^{22} \sigma_{\text{m}}^{2} - \frac{1}{4}\gamma_{ij}^{22} \sigma_{\text{m}}^{4} \) | \( \frac{1}{2}\alpha_{ij}^{22} /\sigma_{\text{m}}^{2} - \frac{1}{4}\gamma_{ij}^{22} \) |
Second approximation
Parameter estimation
Expectation of elements involved in precision formulae when a uniform \( ( f ( p ) = 1 ) \) or a U shaped distribution of allelic frequencies is assumed \( \left( {f( p) = {k \mathord{\left/ {\vphantom {k {2p\left( {1\text{ - }p} \right)}}} \right. \kern-0pt} {2p\left( {1-p} \right)}}} \right) \)
Element | Expectation | |
---|---|---|
Uniform | U shaped | |
E[σ _{ m } ^{2} ] | 1/3 | k |
E[σ _{ m } ^{4} ] | 2/15 | k/3 |
E[ρ _{ m }] | \( 1 - 2\frac{h}{\omega }\theta \) | \( \frac{k}{\omega }\theta \) |
E[ρ _{ m } ^{2} ] | \( \left( {\frac{4\theta }{\omega } + \frac{2}{h}} \right)\left( {\frac{1 + h}{1 + 4h}} \right)^{2} - \frac{4\theta h}{\omega } \) | \( \frac{k}{{\omega^{2} }}\left[ {\theta \left( {\omega - \frac{2h}{\omega }} \right) - 1} \right] \) |
E[ρ _{ m } ^{2} /σ _{ m } ^{2} ] | \( \frac{1}{{\omega^{2} }}\left[ {\theta \left( {\omega - \frac{2h}{\omega }} \right) - 1} \right] \) | \( \frac{k}{{2\omega^{3} }}\left\{ {2\theta + \frac{\omega }{h}} \right\} \) |
The case of unrelated individuals
Non-independence between reference and candidate population, a simple example
Alternative genotypes codification
Based on Additional file 2, the expectation of ζ parameter is \( \frac{k}{4}\left[ {2\log \left( {N_{e} - 1} \right) + 2\frac{{N_{e} \left( {N_{e} - 2} \right)}}{{N_{e} - 1}}} \right] \) for a U-shaped distribution of alleles frequencies and log (N _{ e } − 1) for a uniform distribution.
The case of markers in linkage disequilibrium
So far, following Goddard [16], we considered the situation of n _{ M } independent segments that each carries a single QTL in LD with a single marker. More typically, the genomic information consists of a large number of non-independent markers. This non-independence comes from long-term effects due to bottlenecks, mutations, migrations, etc. and short-term effects due to family structure.
Effective and equivalent numbers of independent loci
Once marker allele frequencies and between-marker LD are estimated in a population of interest, the equivalent number of independent loci which can be estimated from formula (9) and this parameter can be used in models that predict the genetic gain expected from a genomic selection scheme applied to this population.
In the more general situation, prior to the observation of the X _{ r } matrix, a simple approximation for \( n_{{M_{i} }} \) is obtained assuming equal variances σ _{ m } ^{2} = s ^{2}, and using the relation between expected LD and effective population size N _{ e } as derived by Sved [32]: \( {\text{E}}\left[ {2\Updelta_{ml} } \right] = {{\sigma_{m} \sigma_{l} } \mathord{\left/ {\vphantom {{\sigma_{m} \sigma_{l} } {\sqrt {1 + 4N_{e} d_{lm} } }}} \right. \kern-0pt} {\sqrt {1 + 4N_{e} d_{lm} } }} \) with d _{ lm } the distance between ordered loci l and m, such that \( d_{lm} = {{\left| {l - m} \right|L} \mathord{\left/ {\vphantom {{\left| {l - m} \right|L} {n_{M} }}} \right. \kern-0pt} {n_{M} }} \) with L the genome length in Morgan. With those hypotheses, let U = tr[(γn _{ R } R + n _{ M } I)^{−1}] with \( {\mathbf{R}}_{\text{ml}} = \sqrt {{{n_{M} } \mathord{\left/ {\vphantom {{n_{M} } {\left( {n_{M} + 4N_{e} \left| {l - m} \right|L} \right)}}} \right. \kern-0pt} {\left( {n_{M} + 4N_{e} \left| {l - m} \right|L} \right)}}} \).
Towards an exact treatment of linkage disequilibrium
For the candidate c as for the reference individual i, the pair of genetic values may originate from the same parent (and coded on the same chromosome) or not, giving four types of \( \left( {g_{cls} ,g_{cmt} ,g_{ilu,} g_{imv} } \right) \) vectors. In type 1 (\( s = t \, and \, u = v \)), both alleles (belonging to loci \( m\; {\text{and}} \;l \)) of each pair of loci (one for c and one for i) are on the same chromosome (may be from the two fathers, the two dams, c’s father and i’s dam or i’s father and c’s dam). In type 2 (\( s = t \, and \, u \ne v \)), both alleles (belonging to loci \( m \;{\text{and}} \;l \)) of the candidate are on the same chromosome, while alleles of the reference individual i are not on the same chromosome.Type 3 (\( s \ne t \, and \, u = v \)) is the reverse from type 2. In type 4 (\( s \ne t \, and \, u \ne v \)), alleles of loci \( m\; {\text{and}} \;l \) of both individuals and i are on different chromosomes.
For each of these situations, the identity by descent (IBD) status between alleles at locus \( m \) on chromosomes ct and iv, and at locus \( l \) on chromosomes cs and iu are considered. There are four, as follows:
Expectations of products of four allelic values received by two individuals at two loci depending on the IBD status and parental origins of the alleles
\( {\mathcal{S}} \) | \( {\mathcal{T}} \) | \( {\text{E}}\left[ {g_{cls} g_{cmt} g_{ilu} g_{imv} |{\mathcal{S} }\,\& \text{ }\,T} \right] \) |
---|---|---|
\( {\mathcal{S}}_{ml} \) | \( s = t \, and \, u = v \) | \( \left( {1 - p_{m} } \right)\left( {1 - p_{l} } \right)p_{m} p_{l} + \Updelta_{lm} \left( {1 - 2p_{l} } \right)\left( {1 - 2p_{m} } \right) \) |
\( s = t \, \, and \, \, u \ne v \) | \( \left( {1 - p_{m} } \right)\left( {1 - p_{l} } \right)p_{m} p_{l} + \Updelta_{lm} \left( {1 - 2p_{l} } \right)\left( {1 - 2p_{m} } \right) \) | |
\( s \ne t \, and \, u = v \) | \( \left( {1 - p_{m} } \right)\left( {1 - p_{l} } \right)p_{m} p_{l} + \Updelta_{lm} \left( {1 - 2p_{l} } \right)\left( {1 - 2p_{m} } \right) \) | |
\( s \ne t \, and \, u \ne v \) | (1 − p _{ m }) (1 − p _{ l })p _{ m } p _{ l } | |
\( s = t \, and \, u = v \) | \( \Updelta_{lm}^{2} \times {{\left[ {p_{m}^{3} + \left( {1 - p_{m} } \right)^{3} } \right]} \mathord{\left/ {\vphantom {{\left[ {p_{m}^{3} + \left( {1 - p_{m} } \right)^{3} } \right]} {\left[ {p_{m} \left( {1 - p_{m} } \right)} \right]}}} \right. \kern-0pt} {\left[ {p_{m} \left( {1 - p_{m} } \right)} \right]}} \) | |
\( s = t \, and \, u = v \) | \( \Updelta_{lm}^{2} \times {{\left[ {p_{l}^{3} + \left( {1 - p_{l} } \right)^{3} } \right]} \mathord{\left/ {\vphantom {{\left[ {p_{l}^{3} + \left( {1 - p_{l} } \right)^{3} } \right]} {\left[ {p_{l} \left( {1 - p_{l} } \right)} \right]}}} \right. \kern-0pt} {\left[ {p_{l} \left( {1 - p_{l} } \right)} \right]}} \) | |
\( s = t \, and \, u = v \) | \( \Updelta_{lm}^{2} \times \left( {1 - 2p_{m} } \right)\left( {1 - 2p_{l} } \right) \) |
As an illustration, we consider again the situation of a candidate (c), that is the son of one of the n _{ r } individuals in P_{r} and assume that c ^{’s }dam is unrelated to the sire. In formula (2), the summation over the reference individuals i comprises a single term for the sire of the candidate and n _{ r } − 1 terms for the individual that are unrelated to the c members of this reference population.
Based on Additional files 1 and 5, expectations involved in the precision formulae (2) are:
E[x _{ cm } ^{2} X _{ rim } ^{2} ] = p _{ m }(1 − p _{ m }), and
\( \begin{aligned} {\text{E}}\left[ {x_{cl} x_{cm} X_{il} X_{im} } \right] \hfill \\ = \left( {1 - p_{m} } \right)\left( {1 - p_{l} } \right)p_{m} p_{l} + \Updelta_{lm} \left( {1 - 2p_{l} } \right)\left( {1 - 2p_{m} } \right) \hfill \\ + 2\Updelta_{lm}^{2} \left[ {r_{ml} \left( {\frac{{p_{m}^{3} + \left( {1 - p_{m} } \right)^{3} }}{{p_{m} \left( {1 - p_{m} } \right)}} + \frac{{p_{l}^{3} + \left( {1 - p_{l} } \right)^{3} }}{{p_{l} \left( {1 - p_{l} } \right)}}} \right) + \left( {1 - r_{ml} } \right)\left( {1 - 2p_{m} } \right)\left( {1 - 2p_{l} } \right)} \right], \hfill \\ \end{aligned} \) when i is the sire of c; and E[x _{ cm } ^{2} X _{ rim } ^{2} ] = 4[p _{ m }(1 − p _{ m })]^{2} and \( {\text{E}}\left[ {x_{cl} x_{cm} X_{il} X_{im} } \right] = 4\Updelta_{lm}^{2} \left( {1 - 2p_{m} } \right)\left( {1 - 2p_{l} } \right) \), when i and c are unrelated.
Numerical evaluation
Simulation of allele frequencies
In the following numerical evaluation of the formulae derived above, allele frequencies were simulated following an inverse transform sampling (e.g. [32]): \( n_{M} \) allele frequency cumulative distribution function values u _{ m } were simulated in a uniform \( {\mathcal{U}}\left( {0,1} \right) \), and corresponding allele frequencies p _{ m }, i.e. such as \( u_{m} = \int_{{1/2n_{r} }}^{{p_{m} }} {f\left( p \right)dp} \), computed by \( p_{m} = \frac{{\left( {2n_{r} - 1} \right)^{{\left( {2u_{m} - 1} \right)}} }}{{1 + \left( {2n_{r} - 1} \right)^{{\left( {2u_{m} - 1} \right)}} }} \).
Basic situation: no LD and unrelated individuals
Convergence of Taylor series and quality of expectation of the reliability approximations were tested for different population sizes (\( n_{r} = 500, 1000, 1500 \;{\text{and}}\; 2500 \)), numbers of markers (\( n_{M} = 50, 100, 250, 1000, 1500, 2000 \;{\text{and}}\; 2500 \)) and proportions of the phenotypic variance explained by the molecular score (\( \nu^{2} = 0.1, 0.4 \;{\text{and}} \;0.7 \)). Given the set of allele frequencies \( p_{m} \left( {m = 1 \ldots n_{M} } \right) \), genotypes X of n _{ r } + 1 individuals were generated and the G matrix was built. The reliability of the candidate individual GEBV, \( r^{2} = \frac{{{\text{v}}(\hat{q}_{c} |{\mathbf{X}})}}{{{\text{v}}\left( {q_{\text{c}} |{\mathbf{X}}} \right)}} \) was computed as described in the section «Situation analyzed» , as well as approximations considering 1–10 elements in the Taylor series I − D γ + D ^{2} γ ^{2} − D ^{3} γ ^{3}··· The convergence of the series as predicted by the value (lower or higher than 1) of the matrix’s largest eigenvalue was checked numerically, by estimating the mean of this largest eigenvalue from five simulations in each case studied (\( n_{r} = 200\; {\text{to}} \;1000 ;n_{M} = 100 \;{\text{to}}\; 2000 \;{\text{and}} \;\nu^{2} = 0.1, 0.4, 0.7). \)
This limited number of replications was chosen after observation of a very limited variance of this eigenvalue. Finally, the asymptotic values of the suggested approximations [formulae (5) and (6)] were computed using the number of independent segments as described by [4]. The process was repeated 50 times and the means of those exact or approximated reliabilities computed.
Performances of the first approximation \( \left( {\tilde{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right) \) for an unrelated reference population as a function of the number of markers (n _{ M }) and reference population size (n _{ R }), assuming ν ^{2} = 0.4
n _{ M } | n _{ R } | True value \( \left( {{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right) \) | Approximation \( \left( {\hat{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right) \) | Convergence criteria | \( {\text{E}}\left[ {\hat{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right] \) |
---|---|---|---|---|---|
50 | 500 | 0.92 | 2.14 | 2.74 | 2.05 |
50 | 1000 | 0.96 | 2.30 | 2.62 | 2.23 |
50 | 1500 | 0.96 | 2.31 | 2.57 | 2.28 |
50 | 2000 | 0.97 | 2.43 | 2.60 | 2.37 |
100 | 500 | 0.85 | 1.71 | 2.58 | 1.69 |
100 | 1000 | 0.90 | 2.01 | 2.49 | 2.05 |
100 | 1500 | 0.95 | 2.21 | 2.58 | 2.18 |
100 | 2000 | 0.96 | 2.31 | 2.60 | 2.26 |
250 | 500 | 0.67 | 1.09 | 2.53 | 1.09 |
250 | 1000 | 0.81 | 1.56 | 2.53 | 1.55 |
250 | 1500 | 0.85 | 1.72 | 2.54 | 1.72 |
250 | 2000 | 0.89 | 1.96 | 2.53 | 1.97 |
1000 | 500 | 0.32 | 0.39 | 0.61 | 0.38 |
1000 | 1000 | 0.52 | 0.72 | 2.45 | 0.73 |
1000 | 1500 | 0.64 | 0.99 | 2.51 | 0.99 |
1000 | 2000 | 0.71 | 1.18 | 2.51 | 1.18 |
1500 | 500 | 0.25 | 0.28 | 0.27 | 0.28 |
1500 | 1000 | 0.42 | 0.54 | 1.88 | 0.54 |
1500 | 1500 | 0.54 | 0.76 | 2.49 | 0.76 |
1500 | 2000 | 0.61 | 0.91 | 2.50 | 0.91 |
2000 | 500 | 0.20 | 0.22 | 0.20 | 0.22 |
2000 | 1000 | 0.35 | 0.43 | 0.95 | 0.43 |
2000 | 1500 | 0.46 | 0.61 | 2.28 | 0.61 |
2000 | 2000 | 0.54 | 0.76 | 2.48 | 0.76 |
2500 | 500 | 0.16 | 0.17 | 0.16 | 0.17 |
2500 | 1000 | 0.30 | 0.35 | 0.44 | 0.35 |
2500 | 1500 | 0.40 | 0.50 | 1.61 | 0.50 |
2500 | 2000 | 0.49 | 0.65 | 2.40 | 0.66 |
Performances of the second approximation \( \left( {\tilde{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right) \) for an unrelated reference population as a function of the number of markers (n _{ M }) and reference population size (n _{ R }), assuming ν ^{2} = 0.4
n _{ M } | n _{ R } | True value \( \left( {{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right) \) | Approximation \( \left( {\tilde{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right) \) | 10th order approximation | \( {\text{E}}\left[ {\tilde{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right] \) |
---|---|---|---|---|---|
50 | 500 | 0.92 | 0.91 | 0.91 | 0.91 |
50 | 1000 | 0.96 | 0.95 | 0.94 | 0.95 |
50 | 1500 | 0.96 | 0.97 | 0.97 | 0.96 |
50 | 2000 | 0.97 | 0.97 | 0.97 | 0.97 |
100 | 500 | 0.85 | 0.83 | 0.82 | 0.83 |
100 | 1000 | 0.90 | 0.91 | 0.90 | 0.91 |
100 | 1500 | 0.95 | 0.94 | 0.94 | 0.94 |
100 | 2000 | 0.96 | 0.95 | 0.95 | 0.95 |
250 | 500 | 0.67 | 0.71 | 0.67 | 0.71 |
250 | 1000 | 0.81 | 0.83 | 0.81 | 0.82 |
250 | 1500 | 0.85 | 0.88 | 0.87 | 0.88 |
250 | 2000 | 0.89 | 0.90 | 0.90 | 0.90 |
1000 | 500 | 0.32 | 0.41 | 0.31 | 0.40 |
1000 | 1000 | 0.52 | 0.59 | 0.52 | 0.57 |
1000 | 1500 | 0.62 | 0.68 | 0.64 | 0.67 |
1000 | 2000 | 0.69 | 0.73 | 0.70 | 0.73 |
1500 | 500 | 0.24 | 0.32 | 0.23 | 0.31 |
1500 | 1000 | 0.42 | 0.50 | 0.42 | 0.48 |
1500 | 1500 | 0.52 | 0.60 | 0.53 | 0.59 |
1500 | 2000 | 0.60 | 0.66 | 0.61 | 0.67 |
2000 | 500 | 0.19 | 0.26 | 0.17 | 0.28 |
2000 | 1000 | 0.34 | 0.43 | 0.33 | 0.44 |
2000 | 1500 | 0.46 | 0.53 | 0.46 | 0.53 |
2000 | 2000 | 0.53 | 0.60 | 0.54 | 0.61 |
2500 | 500 | 0.16 | 0.22 | 0.14 | 0.23 |
2500 | 1000 | 0.30 | 0.38 | 0.28 | 0.38 |
2500 | 1500 | 0.40 | 0.48 | 0.39 | 0.47 |
2500 | 2000 | 0.47 | 0.55 | 0.48 | 0.56 |
Table 6 shows that the second Taylor series converges always when \( \nu^{2} = 0.4 \). The proposed approximation was generally biased upwards. This over-estimation of the precision was generally limited but increased as the number of markers and the reference population size decreased. The maximum over-estimation observed was 37.5 % (0.22 instead of 0.16 with a standard error less than 0.02). Based on the results in Table 5, it appears that when the first Taylor series converges, the proposed approximation is also slightly over-estimated. The expectation of the approximations, as given in formulae (5) and (6) are very close to the observations.
No LD and non-independence between reference and candidate population
Performances of the first approximation \( \left( {\hat{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right) \) when the parents of candidate belong to the reference population as a function of the number of markers (n _{ M }) and reference population size (n _{ R }), assuming ν ^{2} = 0.4
n _{ M } | n _{ R } | True value \( \left( {{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right) \) | Approximation \( \left( {\hat{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right) \) | 10th order approximation | \( {\text{E}}\left[ {\hat{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right] \) |
---|---|---|---|---|---|
1000 | 500 | 0.37 | 0.42 | 0.58 | 0.47 |
1000 | 1000 | 0.56 | 0.73 | 2.47 | 0.82 |
1000 | 1500 | 0.65 | 0.95 | 2.50 | 1.04 |
1000 | 2000 | 0.72 | 1.17 | 2.52 | 1.26 |
1500 | 500 | 0.31 | 0.34 | 0.33 | 0.37 |
1500 | 1000 | 0.46 | 0.56 | 1.87 | 0.63 |
1500 | 1500 | 0.56 | 0.73 | 2.44 | 0.81 |
1500 | 2000 | 0.62 | 0.87 | 2.50 | 0.96 |
2000 | 500 | 0.27 | 0.29 | 0.27 | 0.32 |
2000 | 1000 | 0.40 | 0.46 | 0.89 | 0.52 |
2000 | 1500 | 0.50 | 0.62 | 2.24 | 0.69 |
2000 | 2000 | 0.57 | 0.76 | 2.48 | 0.84 |
2500 | 500 | 0.24 | 0.25 | 0.24 | 0.27 |
2500 | 1000 | 0.36 | 0.40 | 0.49 | 0.45 |
2500 | 1500 | 0.46 | 0.55 | 1.79 | 0.61 |
2500 | 2000 | 0.52 | 0.67 | 2.35 | 0.74 |
Performances of the second approximation \( \left( {\tilde{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right) \) when the parents of the candidates belong to the reference population as a function of the number of markers (n _{ M }) and reference population size (n _{ R }), assuming ν ^{2} = 0.4
n _{ M } | n _{ R } | True value \( \left( {{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right) \) | Approximation \( \left( {\tilde{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right) \) | 10th order approximation | \( {\text{E}}\left[ {\tilde{r}_{{q_{\text{c}} ,\hat{q}_{\text{c}} }}^{2} } \right] \) |
---|---|---|---|---|---|
1000 | 500 | 0.37 | 0.46 | 0.35 | 0.46 |
1000 | 1000 | 0.53 | 0.60 | 0.54 | 0.61 |
1000 | 1500 | 0.64 | 0.70 | 0.65 | 0.69 |
1000 | 2000 | 0.71 | 0.75 | 0.72 | 0.75 |
1500 | 500 | 0.30 | 0.39 | 0.26 | 0.40 |
1500 | 1000 | 0.47 | 0.55 | 0.46 | 0.51 |
1500 | 1500 | 0.56 | 0.63 | 0.56 | 0.61 |
1500 | 2000 | 0.63 | 0.69 | 0.64 | 0.68 |
2000 | 500 | 0.27 | 0.36 | 0.22 | 0.35 |
2000 | 1000 | 0.40 | 0.49 | 0.38 | 0.48 |
2000 | 1500 | 0.50 | 0.58 | 0.50 | 0.56 |
2000 | 2000 | 0.57 | 0.64 | 0.57 | 0.62 |
2500 | 500 | 0.24 | 0.33 | 0.20 | 0.32 |
2500 | 1000 | 0.34 | 0.44 | 0.31 | 0.45 |
2500 | 1500 | 0.44 | 0.53 | 0.43 | 0.53 |
2500 | 2000 | 0.37 | 0.46 | 0.35 | 0.46 |
Example of the use of the second approach
As an illustration of formula (3) different situations that differ in the relationships between the candidate and reference populations were compared. Coefficients of formula (3) were estimated using the elements in Table 3. An effective reference population size of 200, the genotyping of 10,000 markers and a heritability of 0.4 were assumed. Scenarios included no individuals related to the candidate in the reference population, its sire, both parents, 1–10 half-sibs (or uncles), and a combination of parental and half-sib information.
Equivalent number of independent loci
This number was computed using formula (8), for various effective population sizes (\( N_{e} = 100 \;{\text{to}} \;1000 \)), heritabilities (\( h^{2} = 0.1 \;{\text{to}}\; 0.5 \)), total numbers of loci (\( n_{M} = 1000 \;{\text{to}}\; 10,000) \) and reference population sizes (\( n_{R} = 1000 \;{\text{to}}\; 2500 \)).
- (1)
The trace T of (E[X _{ r } ^{′} X _{ r }] + λ_{β} I)^{−1} is a decreasing function of n_{ r }: as a consequence, the larger is the population size, the smaller is T, which is proportional to the marker effects conditional variances \( v\left(\varvec{\beta}\right) - cov\left( {\varvec{\beta},{\mathbf{y}}} \right)v\left( \varvec{y} \right)^{ - 1} cov( {{\mathbf{y}},\varvec{\beta}}) \)) and the higher is the variance of the estimated molecular score (\( {\text{v}}\left( {q_{c} |{\mathbf{y}}} \right) = {\mathbf{x}}_{{\mathbf{c}}} cov\left( {\varvec{\beta},{\mathbf{y}}} \right)v\left( \varvec{y} \right)^{ - 1} cov\left( {{\mathbf{y}},\varvec{\beta}} \right){\mathbf{x}}^{\prime}_{{\mathbf{c}}} ) \).
- (2)
The trace T is always higher in the situation of LD than for independent markers \( \left( {{\text{T}}_{LD} > {\text{T}}_{LE} } \right) \).
- (3)The rate of decrease is higher for T_{ LD } than for T_{ LE }. On the whole, the reliability for a given number of observed markers corresponds to the reliability that is reached with a larger number of independent loci when the size of the reference population is larger.
Discussion
The objective of this paper was to explore approximations of the precision of genomic selection when the selection candidate has relatives in the reference population. Two approximations were developed and numerically compared.
These approximations were based on Taylor expansions of a matrix inverse M ^{−1}. In both cases, the initial matrix is the sum of the identity matrix and a perturbation (M = I + E). Convergence of these series is not guaranteed and depends on the behavior of the perturbation (I − E + E ^{2} − E ^{3} → (I + E)^{−1} if \( {\mathbf{E}}^{\text{t}} \to 0 \) when t → ∞). With the first approximation, derived from the appendix in [18], this convergence failed when the number of makers was too small (less than 1500 in our example) or the heritability was greater than 0.1. This was only observed when ν ^{2} = 0.7 with the second approximation. This is fully consistent with the deviation to one of the largest eigenvalue of the E matrix.
The expectation of the proposed approximation, when data were simulated with the model corresponding to the hypotheses underlying their algebraic development, was very close to the mean value after 50 simulations. Thus, extremely fast estimation of the precision is possible, which allows intensive optimization and comparison of selection schemes.
When individuals are unrelated and markers are in linkage equilibrium, we obtain an estimation of the GEBV accuracy which differs from that of Goddard et al. [18]. This is surprising since that approach was said to be based on the Taylor approximation used here. Their formula may be obtained in a simpler way [see Additional file 6]. However, relaxing the assumption of “absence of between-individual relationship” is not straightforward using this approach.
A strong limit of our new approximation comes from the limitation to the first order term of the Taylor series. Deriving algebra was only possible at this stage. The side effect is that no genotypic covariance terms between reference individuals appear in this approximation. As a consequence, only the direct relationships between candidate and reference individuals play a role in the estimation, but not the structure within the reference population. This is unfortunate, because accuracies of genomic prediction are obviously affected by the construction of the reference population. Our last numerical example, in which there is a linear trend with the number of half-sibs, reveals this drawback: two half-sibs of the candidates are treated as unrelated and the information that they carry is just the double of that of a single half-sib. Future developments should focus on this limitation, for instance to derive the expectation of the x _{ c } C ^{−1} B ^{2} x _{ c } ^{′} term.
Our work focused on the BLUP precision of the molecular score \( r^{2} \left( {q_{ci} ,\hat{q}_{ci} |{\mathbf{X}}} \right) = \frac{{v\left( {\hat{q}_{ci} } \right)}}{{v\left( {q_{ci} } \right)}} \) but left aside the proportion of the genetic variance that is captured by the markers \( \left( {\frac{{v ( {q_{ci} } )}}{{v ( {g_{ci} } )}}} \right) \). This last term could be treated as in Goddard et al. [18]: \( \frac{{v\left( {q_{ci} } \right)}}{{v ( {g_{ci} })}} = b = \frac{{n_{M} }}{{n_{M} + M_{e} }} \) with M _{ e } the number of independent segments. As noted in the section on the general framework, the quantity \( \frac{{v ( {\hat{q}_{ci} } )}}{{v ( {g_{ci} } )}} = b \times r ( {q_{ci} ,\hat{q}_{ci} |{\mathbf{X}}} ) \) is only an approximation of these GEBV reliabilities i.e. \( r^{2} \left( {g_{ci} ,\hat{q}_{ci} |{\mathbf{X}}} \right) = \frac{{cov^{2} \left( {g_{ci} ,\hat{q}_{ci} |{\mathbf{X}}} \right)}}{{v\left( {g_{ci} |{\mathbf{X}}} \right)v\left( {\hat{q}_{ci} |{\mathbf{X}}} \right)}} \). Equality between those quantities is obtained when \( {\mathbf{X}} = {\mathbf{W}} \)(identity between statistical and genetical models), a condition assumed in Goddard [16] where markers and QTL are modeled as a series of uncorrelated pairs.
All the developments shown in this paper are based on the hypothesis that the reliability of GEBV based on non-independent markers for a trait controlled by \( n_{Q} \) QTL that are in incomplete LD with the markers can be approached by the reliability of GEBV when there are n _{ M } independent segments that carry a single QTL in LD with a single marker. A few difficulties arose when applying this approach proposed by Goddard [16]. How many independent markers should be considered? The reasoning in Goddard [16] was based on the idea of an effective number of loci (M _{ e }) corresponding to a given variance of realized relationships. Here, we proposed the alternative equivalent number of independent loci (M _{ i }) which corresponds to a given reliability. We showed that this M _{ i } number depends on the size of the reference population and on heritability, a dependence that does not occur with M _{ e }. If we invert the argument, controlling the level of realized relationships variance with the effective number of loci (M _{ e }) does not seem to be a good approach to control the estimated GEBV reliability.
As detailed by Hayes et al. [17], the effective number of independent chromosome segments depends on the population structure. The higher is the mean relationship level, the smaller is this effective number. However, we suggest the use of this number as estimated from a set of unrelated individuals, or of its expectation prior to any observation, assuming independence between individuals. Without formal proof, the idea was that long-term LD was considered by using an effective (or equivalent) number of independent loci while short-term non-independence was taken into account with our formalization of the matrix’s expectations that is developed in Additional file 1. A complete proof of the procedure is still needed.
Regardless of the definition of \( M_{e} \) or M _{ i }, there is no reason that the number of independent loci must equal the number of QTL, which is unknown, contrary to the hypothesis about pairs of marker-QTL (in practice, since the QTL effects are random variables, many segments will only have very small effects on the trait, thus simulating the more likely situation of a limited number of “real” QTL). Equating \( {\mathbf{X}} \) and W as well as σ _{β} ^{2} and σ _{α} ^{2} has no clear justification. The variance \( v\left( {\hat{q}_{c} |{\mathbf{X}}} \right) \) of the molecular score should not be σ _{β} ^{2} x _{ c } X _{ r } ^{′} (X _{ r } X _{ r } ^{′} + Iλ_{β})^{−1} X _{ r } x _{ c } ^{′} but σ _{α} ^{2} x _{ c } X _{ r } ^{′} (X _{ r } X _{ r } ^{′} + Iλ_{β})^{−1}(W _{ r } W _{ r } ^{′} + Iλ_{α})(X _{ r } X _{ r } ^{′} + Iλ_{β})^{−1} X _{ r } x _{ c } ^{′} . This other formula assembles two sets of unknown parameters: the variances σ _{α} ^{2} and σ _{β} ^{2} , and the genotypes X and W. It is often assumed that \( {{\sigma }}_{{{\upbeta }}}^{2} = \sigma_{g}^{2} /\left( {n_{M} \bar{\tau }} \right) \) (e.g. [1]), which results in an overestimation of the \( {{\uplambda }}_{{{\upbeta }}} \) parameter since LD is not considered. Working on the number of independent loci (\( M_{e}\, {\text{or}}\,M_{i} \)) apparently solves this difficulty. The QTL variance \( {{\sigma }}_{{{\alpha }}}^{2} = \sigma_{g}^{2} /\left( {n_{Q} \bar{\tau }} \right) \) could be derived based on a hypothesis about the number of QTL. The situation is more difficult for the genotype matrices since the W _{ r } matrix is not observed.
If the framework considered so far (n _{ M } markers-QTL pairs with strong LD within pairs and no LD between pairs) is partly retained, a slight improvement is possible considering the element b of the genetic variability explained by SNPs. The idea would be to replace, in the formulae used in this paper, σ _{q} ^{2} by b × σ _{g} ^{2} . Element b can be derived by considering that the markers’ (β) and QTL’ (α) effects are fixed in the genetic and statistical models. Leaving aside the singularity of X _{ r } ^{′} X _{ r } when the number of SNPs is large, the marker effects are now estimated by \( {\hat{\varvec{\upbeta }}} = \left( {{\mathbf{X}}^{\prime}_{{\mathbf{r}}} {\mathbf{X}}_{{\mathbf{r}}} } \right)^{ - 1} {\mathbf{X}}^{\prime}_{{\mathbf{r} }}{\mathbf{y}} \) and the molecular score defined as \( {\hat q}= {\mathbf{X}}_{{\mathbf{r}}} {\hat{\varvec{\upbeta }}} \), while the genetic value was g = W _{ r } α. Given the genotype matrices, the sample genetic variability is v _{ g } = α ^{ ′ } W _{ r } ^{′} W _{ r } α and the sample molecular score variability y ^{ ′ } X _{ r }(X _{ r } ^{′} X _{ r })^{−1} X _{ r } ^{′} y with an expectation v _{ q } = α ^{′} W _{ r } ^{′} X _{ r }(X _{ r } ^{′} X _{ r })^{−1} X _{ r } ^{′} W _{ r } α. The part of the genetic variability explained by the SNPs is the ratio b = v _{ q }/v _{ g }.
Expectations of the matrices’ product elements \( \left\{ {{\mathbf{X}}^{\prime}_{{\mathbf{r}}} {\mathbf{X}}_{{\mathbf{r}}} } \right\}_{{\varvec{ml}}} \) are \( 2n_{r} \Updelta_{ml} \) off diagonal and 2n _{ r } p _{ m }(1 − p _{ m }) = n _{ r } σ _{ m } ^{2} in the diagonal, with similar expressions for \( {\mathbf{W^{\prime}_{\mathbf{r}}}} {\mathbf{X}}_{{\mathbf{r}}} \) and W _{ r } ^{′} W _{ r } elements.
Following Goddard [16], approximating expectations of the matrices’ functions by the function of their expectation, and assuming that (1) markers are independent, (2) each QTL q is in LD with only one marker m(q), with a LD value \( \Updelta_{qm\left( q \right)} \), and (3) individuals are unrelated: v _{ g } = n _{ r } ∑ _{ q } α _{ q } ^{2} σ _{ q } ^{2} , we get \( v_{q} \sim 4n_{r} \sum\nolimits_{q} {\frac{{\Updelta_{qm\left( q \right)}^{2} }}{{\sigma_{m\left( q \right)}^{2} }}\alpha_{q}^{2} } = n_{r} \sum\nolimits_{q} {r_{qm\left( q \right)}^{2} \alpha_{q}^{2} \sigma_{q}^{2} } \), and \( {\text{b}} = \frac{{\mathop \sum \nolimits_{q} r_{qm\left( q \right)}^{2} \alpha_{q}^{2} \sigma_{q}^{2} }}{{\mathop \sum \nolimits_{q} \alpha_{q}^{2} \sigma_{q}^{2} }} \), corresponding to Eq. (4) in [16].
The ratio b is the weighted mean of LD r ^{2}. Unfortunately, neither α _{ q } ^{2} nor σ _{ q } ^{2} are known. The unweighted mean \( \frac{{\mathop \sum \nolimits_{q} r_{qm\left( q \right)}^{2} }}{{n_{q} }} = \bar{r}^{2} \) may be a fruitful approximation. Following Sved [33], the expectation of \( r_{qm\left( q \right)}^{2} \) is \( \frac{1}{{1 + 4N_{e} c}} \) with c being the distance, in Morgan, between the QTL and its marker. Let L be the total length of the genome, and assume an equal distance L/n _{ M } between each successive marker \( {\text{b}}\sim \int_{0}^{{L/2n_{M} }} {\frac{1}{{1 + 4N_{e} c}}\frac{1}{{{L \mathord{\left/ {\vphantom {L {2n_{M} }}} \right. \kern-0pt} {2n_{M} }}}}dc = \frac{{n_{M} }}{{2N_{e} L}}\left[ {{ \log }\left( {1 + {{2N_{e} L} \mathord{\left/ {\vphantom {{2N_{e} L} {n_{M} }}} \right. \kern-0pt} {n_{M} }}} \right)} \right]} \).
Expectation of the ratio of variances vs. the ratio of the variance expectations considering different reference population sizes and numbers of markers (ν ^{2} = 0.4, 50 simulations)
n _{ r } | n _{ M } | \( {\text{E}}\left[ {{\text{v}}\left( {\hat{q}_{\text{c}} } \right)/{\text{v}}\left( {q_{\text{c}} } \right)} \right] \) | \( {\text{E}}\left[ {{\text{v}}\left( {\hat{q}_{\text{c}} } \right)} \right]/{\text{E}}\left[ {{\text{v}}\left( {q_{\text{c}} } \right)} \right] \) |
500 | 1000 | 0.403 | 0.401 |
1000 | 1000 | 0.726 | 0.725 |
1500 | 1000 | 1.010 | 1.008 |
2000 | 1000 | 1.212 | 1.212 |
500 | 1500 | 0.270 | 0.269 |
1000 | 1500 | 0.535 | 0.534 |
1500 | 1500 | 0.753 | 0.753 |
2000 | 1500 | 0.944 | 0.944 |
500 | 2000 | 0.213 | 0.213 |
1000 | 2000 | 0.414 | 0.413 |
1500 | 2000 | 0.597 | 0.597 |
2000 | 2000 | 0.760 | 0.759 |
500 | 2500 | 0.175 | 0.175 |
1000 | 2500 | 0.349 | 0.348 |
1500 | 2500 | 0.515 | 0.514 |
2000 | 2500 | 0.670 | 0.669 |
The theory presented here was developed by considering a single selection candidate. When candidates are diversely related to the reference population, as suggested in Goddard et al. [18], the candidates should be examined one by one. Moreover, non-independence between candidates should be considered. A further step towards the modeling of genomic selection could be an approximation of the mean genetic values of selected individuals when GEBV reliabilities are heterogeneous.
A few other hypotheses were made in this paper, including additivity and i.i.d. of QTL effects, and the use of GBLUP. As long as the objective is to model and optimize breeding plans, only relative values are interesting and we assumed that these hypotheses were not critical.
Conclusions
The objective of this paper was to provide a further step towards the development of deterministic models that describe genomic breeding plans. Such deterministic models carry low computational burden and thus allow design optimization through intensive numerical exploration.
We proposed two alternative approximations of the estimation of GEBV reliability in the case of non-independence between candidate and reference populations. Both were derived from the Taylor series heuristic approach suggested by Goddard [16]. A numerical exploration of their properties showed that the series were not equivalent in terms of convergence to the exact reliability, that the approximations may overestimate GEBV precision and that they perfectly converged toward their theoretical expectations.
Formulae derived for these approximations were simple to handle in the case of independent markers. A few parameters that describe the markers’ genotypic variability (allele frequencies, linkage disequilibrium) can be estimated from genomic data corresponding to the population of interest or estimated after assumption about their distribution.
When markers are not in linkage equilibrium (i.e. there is LD), replacing the real number of markers and QTL by an effective or equivalent number of independent loci, as proposed by Goddard [16] and Hayes et al. [17], is a practical solution. Research efforts are still needed to overcome some strong limits of this approach.
Declarations
Acknowledgements
This work was partly done when the author was on sabbatical leave in the Animal Genetic and Breeding Unit (AGBU) in Armidale, Australia. This sabbatical was supported by a grant from AGBU and from INRA (métaprogramme Selgen). Andrew Swan, Julius van der Werf, Mike Goddard, Anne Ricard and Bruno Goffinet are thanked for their many useful comments. Rob Banks is particularly thanked for his help at many levels.
Competing interests
The author declares that he has no competing interests.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
- Meuwissen THE, Hayes BJ, Goddard ME. Prediction of total genetic value using genome-wide dense marker maps. Genetics. 2011;157:1819–29.Google Scholar
- Schaeffer LR. Strategy for applying genome-wide selection in dairy cattle. J Anim Breed Genet. 2006;123:218–23.View ArticlePubMedGoogle Scholar
- König S, Simianer H, Willam A. Economic evaluation of genomic breeding programs. J Dairy Sci. 2009;92:382–91.View ArticlePubMedGoogle Scholar
- McHugh N, Meuwissen THE, Cromie AR, Sonesson AK. Use of female information in dairy cattle genomic breeding programs. J Dairy Sci. 2011;94:4109–18.View ArticleGoogle Scholar
- de Roos P, Schrooten WC, Veerkamp RF, van Arendonk JAM. Effects of genomic selection on genetic improvement, inbreeding, and merit of young versus proven bulls. J Dairy Sci. 2011;94:1559–67.View ArticlePubMedGoogle Scholar
- Pryce JE, Goddard ME, Raadsma HW, Hayes BJ. Deterministic models of breeding scheme designs that incorporate genomic selection. J Dairy Sci. 2010;93:5455–66.View ArticlePubMedGoogle Scholar
- Meuwissen THE, Hayes BJ, Goddard ME. Accelerating improvement of livestock with genomic selection. Annu Rev Anim Biosci. 2013;1:221–37.View ArticlePubMedGoogle Scholar
- Sonesson AK, Meuwissen THE. Testing strategies for genomic selection in aquaculture breeding programs. Genet Sel Evol. 2009;41:37.View ArticlePubMedPubMed CentralGoogle Scholar
- Ibánẽz-Escriche N, Fernando RL, Toosi A, Dekkers JCM. Genomic selection of purebreds for crossbred performance. Genet Sel Evol. 2009;41:12.View ArticlePubMedPubMed CentralGoogle Scholar
- Wolc A, Arango J, Settar P, Fulton JE, O’Sullivan NP, Preisinger R, et al. Persistence of accuracy of genomic estimated breeding values over generations in layer chickens. Genet Sel Evol. 2011;43:33.View ArticleGoogle Scholar
- Tribout T, Larzul C, Phocas F. Efficiency of genomic selection in a purebred pig male line. J Anim Sci. 2012;45:4164–76.View ArticleGoogle Scholar
- Shumbusho F, Raoul J, Astruc JM, Palhiere I, Elsen JM. Potential benefits of genomic selection on genetic gain of small ruminant breeding programs. J Anim Sci. 2013;91:3644–57.View ArticlePubMedGoogle Scholar
- Daetwyler HD, Villanueva B, Woolliams JA. Accuracy of predicting the genetic risk of disease using a genome-wide approach. PLoS One. 2008;3:e3395.View ArticlePubMedPubMed CentralGoogle Scholar
- Legarra A, Robert-Granié C, Manfredi E, Elsen JM. Performance of genomic selection in mice. Genetics. 2008;180:611–8.View ArticlePubMedPubMed CentralGoogle Scholar
- Van Raden PM. Efficient methods to compute genomic predictions. J Dairy Sci. 2008;91:4414–23.View ArticleGoogle Scholar
- Goddard ME. Genomic selection: prediction of accuracy and maximisation of long term response. Genetica. 2009;136:245–57.View ArticlePubMedGoogle Scholar
- Hayes B, Visscher P, Goddard M. Increased accuracy of artificial selection by using the realized relationship matrix. Genet Res. 2009;91:47–60.View ArticleGoogle Scholar
- Goddard ME, Hayes BJ, Meuwissen THE. Using the genomic relationship matrix to predict the accuracy of genomic selection. J Anim Breed Genet. 2011;128:409–21.View ArticlePubMedGoogle Scholar
- Buch LH, Kargo M, Berg P, Lassen J, Sørensen AC. The value of cows in reference populations for genomic selection of new functional traits. Animal. 2011;6:880–6.View ArticleGoogle Scholar
- Clark SA, Hickey JM, Daetwyler HD, van der Werf JHJ. The importance of information on relatives for the prediction of genomic breeding values and the implications for the makeup of reference data sets in livestock breeding schemes. Genet Sel Evol. 2012;44:4.View ArticlePubMedPubMed CentralGoogle Scholar
- Wientjes YCJ, Veerkamp RF, Calus MPL. The effect of linkage disequilibrium and family relationships on the reliability of genomic prediction. Genetics. 2013;193:621–31.View ArticlePubMedPubMed CentralGoogle Scholar
- Habier D, Fernando RL, Garrick DJ. Genomic BLUP decoded: a look into the black box of genomic prediction. Genetics. 2013;194:597–607.View ArticlePubMedPubMed CentralGoogle Scholar
- Erbe M, Gredler B, Seefried FR, Bapst B, Simianer H. A function accounting for training set size and marker density to model the average accuracy of genomic prediction. PLoS One. 2013;8:e81046.View ArticlePubMedPubMed CentralGoogle Scholar
- Weller JI. Economic aspects of animal breeding. London: Chapman & Hall; 1994.Google Scholar
- Dekkers JCM. Prediction of response to marker-assisted and genomic selection using selection index theory. J Anim Breed Genet. 2007;124:331–41.View ArticlePubMedGoogle Scholar
- Daetwyler HD, Pong-Wong R, Villanueva B, Woolliams A. The impact of genetic architecture on genome-wide evaluation methods. Genetics. 2010;185:1021–31.View ArticlePubMedPubMed CentralGoogle Scholar
- Habier D, Fernando RL, Dekkers JCM. The impact of genetic relationship information on genome-assisted breeding values. Genetics. 2007;177:2389–97.PubMedPubMed CentralGoogle Scholar
- Pszczola M, Strabel T, Mulder HA, Calus MPL. Reliability of direct genomic values for animals with different relationships within and to the reference population. J Dairy Sci. 2012;95:389–400.View ArticlePubMedGoogle Scholar
- Wientjes YCJ, Veerkamp RF, Bijma P, Bovenhuis H, Schrooten C, Calus MPL. Empirical and deterministic accuracies of across-population genomic prediction. Genet Sel Evol. 2015;47:5.View ArticlePubMedPubMed CentralGoogle Scholar
- Gianola D, De Los Campos G, Hill WG, Manfredi E, Fernando R. Additive genetic variability and the bayesian alphabet. Genetics. 2009;183:347–63.View ArticlePubMedPubMed CentralGoogle Scholar
- Visscher PM, Medland SE, Ferreira MAR, Morley KI, Zhu G, Cornes BK, et al. Assumption-free estimation of heritability from genome-wide identity-by-descent sharing between full siblings. PLoS Genet. 2006;2:e41.View ArticlePubMedPubMed CentralGoogle Scholar
- Sigman K. Lecture notes introduction to discrete-time Markov chains. http://www.columbia.edu/~ks20/stochastic-I/stochastic-I-MCI.pdf. 2009.
- Sved JA. Linkage disequilibrium and homozygosity of chromosome segments in finite populations. Theor Popul Biol. 1971;2:125–41.View ArticlePubMedGoogle Scholar
- Wright S. The distribution of gene frequencies in populations. Proc Natl Acad Sci USA. 1937;23:307–20.View ArticlePubMedPubMed CentralGoogle Scholar
- La Gillois M. relation d’identité en génétique. Ann Inst Henri Poincaré. 1964;B2:1–94.Google Scholar
- Harris DL. Genotypic covariances between inbred relatives. Genetics. 1964;50:1319–48.PubMedPubMed CentralGoogle Scholar
- Jacquard A. Logique du calcul des coefficients d’identité entre deux individus. Populations. 1966;21:751–76.View ArticleGoogle Scholar