Abstract
Background
Singlestep methods provide a coherent and conceptually simple approach to incorporate genomic information into genetic evaluations. An issue with singlestep methods is compatibility between the markerbased relationship matrix for genotyped animals and the pedigreebased relationship matrix. Therefore, it is necessary to adjust the markerbased relationship matrix to the pedigreebased relationship matrix. Moreover, with data from routine evaluations, this adjustment should in principle be based on both observed marker genotypes and observed phenotypes, but until now this has been overlooked. In this paper, I propose a new method to address this issue by 1) adjusting the pedigreebased relationship matrix to be compatible with the markerbased relationship matrix instead of the reverse and 2) extending the singlestep genetic evaluation using a joint likelihood of observed phenotypes and observed marker genotypes. The performance of this method is then evaluated using two simulated datasets.
Results
The method derived here is a singlestep method in which the markerbased relationship matrix is constructed assuming all allele frequencies equal to 0.5 and the pedigreebased relationship matrix is constructed using the unusual assumption that animals in the base population are related and inbred with a relationship coefficient γand an inbreeding coefficient γ/2. Taken together, this γparameter and a parameter that scales the markerbased relationship matrix can handle the issue of compatibility between markerbased and pedigreebased relationship matrices. The full loglikelihood function used for parameter inference contains two terms. The first term is the REMLloglikelihood for the phenotypes conditional on the observed marker genotypes, whereas the second term is the loglikelihood for the observed marker genotypes. Analyses of the two simulated datasets with this new method showed that 1) the parameters involved in adjusting markerbased and pedigreebased relationship matrices can depend on both observed phenotypes and observed marker genotypes and 2) a strong association between these two parameters exists. Finally, this method performed at least as well as a method based on adjusting the markerbased relationship matrix.
Conclusions
Using the full loglikelihood and adjusting the pedigreebased relationship matrix to be compatible with the markerbased relationship matrix provides a new and interesting approach to handle the issue of compatibility between the two matrices in singlestep genetic evaluation.
Introduction
Singlestep methods for genetic evaluation [13] have recently become popular because they provide an approach to incorporate genomic information into genetic evaluations that is both coherent and conceptually simple. A singlestep method extends the usual pedigreebased method by replacing the additive relationship matrix constructed from pedigree by an additive relationship matrix that combines the markerbased relationship matrix for genotyped animals with the pedigreebased relationship matrix.
An issue with a singlestep method is compatibility between the markerbased relationship matrix for genotyped animals and the pedigreebased relationship matrix [46]. To handle this problem, it is necessary to determine which allele frequencies should be used in the markerbased relationship matrix and to adjust this matrix to the pedigreebased relationship matrix. In theory, one should use the allele frequencies in the founder population of the pedigree (base animals) for the markerbased and pedigreebased relationships to be compatible, but these allele frequencies are rarely available in practice since base animals are not genotyped. Chen et al. [5] and Forni et al. [7] concluded that using the observed allele frequencies improved accuracy of prediction compared to using allele frequencies equal to 0.5. Studies on how to adjust the markerbased relationship matrix to be compatible with the submatrix of the pedigreebased relationship matrix for genotyped animals have been reported [46]. These adjustments consisted of scaling and adding a number to all elements in the markerbased relationship matrix based on equating means of diagonal and offdiagonal elements in the two matrices, and it was demonstrated that accuracy of prediction increased and bias decreased. In relation to this problem of compatibility between markerbased and pedigreebased relationship matrices, with data from routine evaluations, selection affects the allele frequencies over time, and in principle both observed marker genotypes and observed phenotypes contain information about allele frequencies in the base population. Therefore, studies on the adjustment of the markerbased relationship matrix to the pedigreebased relationship matrix have overlooked the fact that it should in principle incorporate information on observed phenotypes.
This work explores two possibilities to solve the problem, i.e. 1) in which the pedigreebased relationship matrix is adjusted to the markerbased relationship matrix and 2) in which the singlestep genetic evaluation is extended by using a joint likelihood of observed phenotypes and observed marker genotypes. This results in a singlestep method in which the markerbased relationship matrix is constructed assuming all allele frequencies are equal to 0.5 and the pedigreebased relationship matrix is constructed using the unusual assumption that animals in the base population are related and inbreed with a relationship coefficient γ and an inbreeding coefficient γ / 2. The loglikelihood function used for parameter inference contains two terms. The first term is the REMLloglikelihood for the phenotypes conditional on the observed marker genotypes and the second term is the loglikelihood for the observed marker genotypes. The performance of the proposed method was evaluated using two simulated datasets.
Methods
This section presents, first, the statistical model on which the singlestep methods are based along the lines of Christensen and Lund [3], then the proposal on how to adjust the pedigreebased relationship matrix and finally the adjustment of the markerbased relationship matrix as previously used.
Singlestep model
The marker genotypes are summarised into a marker matrix m, where m_{ij} = −1, 0 or 1 if SNP j of individual i is 11, 12, or 22, respectively. In the following, capital M and lowercase m indicate whether marker genotypes are considered as random variables or as nonrandom variables (observed variables or integration variables), respectively.
Let us consider the simple model
where y is the vector of phenotypes, μ is the general mean, 1 is a vector of ones, a is the vector of breeding values, Z is an incidence matrix, and e is the vector of residual errors. The breeding value may be decomposed into a = g + a_{r}, where g is the vector of genomic effects and a_{r} = (a − g) is the vector of residual polygenic effects. The residual polygenic effects
where v_{j} = v_{j,j} is a parameter and
where
By rearranging terms and using that
where f(y∣m^{o}) is defined implicitly. The first term f(y∣m^{o}) is the density for the phenotypes given the observed marker genotypes, whereas the second term f(m^{o}) is the multivariate Gaussian density for the observed marker genotypes.
Singlestep methods are based on the density f(y∣m^{o}), which may be written (see [3]) as
where the vector a has a mean zero and a variancecovariance matrix H_{ω} with inverse
with G_{ω} = (1 − ω)G(m^{o}) + ωA_{11}. In the above formula, the subdivision is made according to genotyped and nongenotyped
animals. The matrices
Issues raised by the singlestep method are 1) how are the allele frequencies ρ that are used in G(m^{o}) obtained and 2) how is the required compatibility between G(m^{o}) and A_{11} that is evident in equation (1) reached. To investigate these issues, the joint density of observed phenotypes and marker genotypes, f(y, m^{o}) = f(y∣m^{o})f(m^{o}), is taken as the starting point. From this, the full marginal loglikelihood becomes
where
with n_{1} denoting the number of genotyped animals, is the loglikelihood of the observed marker
genotypes. Thus, the parameter ρ enters into both terms of the full loglikelihood in equation (2), and this implies
that estimation of allele frequencies ρ should be in principle based on this loglikelihood. In particular, if selection
has been performed the phenotypes will contain information about the allele frequencies
which may cause bias if ignored. However, estimating ρ and s = s(v) by maximising
Furthermore, when computing BLUP breeding values
Adjusting the pedigreebased relationship matrix
The derivation of the proposed adjustment of the pedigreebased relationship matrix is based on assigning priors on the highdimensional parameters ρ and v in the previously described singlestep model, and then considers the first and second order moments of the marginal distribution of M(integrating ρ and v). Appendix A shows that the resulting marker distribution satisfies
with
with
Based on the derivation above (see also Appendix A), the full marginal loglikelihood function for parameter inference becomes
where
is the loglikelihood of the observed marker genotypes. Using the loglikelihood in
equation (4) instead of the loglikelihood in equation (2) makes the estimation of
parameters feasible in practice by numeric maximization methods since the highdimensional
parameters ρ and v are replaced by two parameters, a parameter γ that determines the relationship and inbreeding of individuals in the base population
and a parameter
The computations require algorithms that compute A(γ)^{−1} and A_{11}(γ) efficiently. Computations of inbreeding coefficients with related base animals have been considered by [10,11] in the context of incomplete pedigrees. In Appendix B, algorithms are presented that extend the approaches of Quaas [12] for computing A^{−1} and of Colleau [13] for computing A_{11} to the case where base animals are related and inbred.
Maximisation of the full loglikelihood in equation (4) is done by first specifying
a discrete threedimensional grid of values for parameters
For faster computing, an alternative to using the full loglikelihood is to determine
parameters γ and
which has the solution
Substituting
which has to be maximised numerically to estimate γ.
Adjusting the markerbased relationship matrix (Gadjust)
Alternatively, an adjustment of the form G_{a} = Gβ + α is used, where G is the markerbased relationship matrix with allele frequencies
for α and β, where
Simulated example 1
This example is deliberately simple and very extreme, and constructed for the purpose
of showing that parameter estimates of γ and
The base population consists of two individuals, one sire and one dam with 15 biallelic markers, and their genotypes are simulated assuming independence between markers and with equal allele frequencies. Furthermore, it is assumed that the markers are independently inherited, are all QTL with allele substitution effect equal to 1, and heritability of the phenotype is equal to 1.
The two base individuals produce 100 offspring (generation 1) that all have observed phenotypes. The two individuals in generation 1 with the largest own phenotype value are selected as parents and produce 100 offspring (generation 2) that are all genotyped.
Two different approaches to estimate parameters are compared. In the first approach,
all parameters are estimated using the full loglikelihood in equation (4). In the
second approach, γ and
Simulated example 2
This example is inspired by a pig nucleus breeding scheme and consists of five generations in which all animals have recorded phenotypes. In each generation, 150 boars are mated to 1 500 sows to produce 15 000 offspring (50/50 males/females). For the next generation, boars with a high own phenotype value are chosen and sows are selected at random. The last three generations of selected boars are genotyped and a sixth generation of 300 candidate boars are also genotyped. The breeding value is the sum of 500 independent QTL effects simulated from a Gamma(5.4,0.42) distribution, and the heritability of the phenotype is 0.22. This dataset is described in more detail in Christensen and Lund [3].
Three different approaches are compared. In the first approach, all parameters are
estimated using the full loglikelihood in equation (4). In the second approach, parameters
γ and
Results
Simulated example 1
Table 3 shows that both γand
Table 3. Parameter estimates obtained with simulated dataset 1
Simulated example 2
Table 4 shows that γ and
Figure 1. Profile loglikelihood for γ and
Figure 2. Profile loglikelihood for ω. The profile loglikelihood function for parameter ω based on the full loglikelihood in equation (5); the plot is constructed with values in a discrete grid that have been standardised such that the maximum value is zero.
Table 4. Parameter estimates and prediction performance obtained with simulated dataset 2
Estimated parameters obtained with the Gadjust approach are also shown in Table 4, where it should be noted that
In terms of prediction performance, the correlation between predicted and true breeding value and the regression coefficient were 0.536 and 1.17, respectively, in both cases when the pedigreebased matrix was adjusted, but 0.493 and 1.34, respectively, when the Gadjust approach was used.
Discussion
This paper provides a coherent approach to handle the issue of compatibility between
pedigreebased and markerbased relationship matrices in singlestep genetic evaluation.
The approach is computationally fast and feasible for large datasets, as is the case
for previously developed singlestep methods. Parameter γ can adjust the pedigreebased relationship matrix to the markerbased relationship
matrix that is scaled by parameter
Using the approach developed in this paper may also provide insight on performance of other approaches for making markerbased and pedigreebased relationship matrices compatible in singlestep genetic evaluation. Examples of such approaches are the adjustments of the markerbased relationship matrix reported in [46,14,15] and also investigated here, and the approach by Meuwissen et al. [16]. In [16], first, an average of position specific identical by descent matrices based on linkage information [17] (hereafter denoted G_{FG}) were used instead of the combined relationship matrix (1), and second, the final method consisted of replacing the pedigreebased relationship matrix A by the matrix G_{FG} in matrix (1) and then adjusting the markerbased relationship matrix to the G_{FG} matrix. Based on simple simulated datasets, without selection, the method resulted in high accuracy and low bias. The computation of G_{FG} is computationally burdensome and therefore this approach was not considered here. It should be noted that for these other approaches, the combined relationship matrix is constructed based on observed marker genotypes only, and the inference on remaining parameters is based on the loglikelihood conditional on the observed marker genotypes. The approach in this paper may provide guidelines on when it is safe to base inference on such procedures in which, first the observed marker genotypes are used to construct a combined relationship matrix, and second the remaining inference is based on the loglikelihood for the phenotypes conditional on observed marker genotypes.
Evaluation of accuracy and bias of prediction with the simulated dataset 2 shows that
the approaches based on adjusting the pedigreebased relationship matrix, where
This paper suggests that the pedigreebased relationship matrix should be adjusted
instead of the markerbased relationship matrix. On the one hand, from a practical
point of view this would be rather inconvenient since standard software used for REML
estimation and BLUP would have to be modified, but on the other hand, it is conceptually
simpler and may be easier to extend. For example, when the population consists of
a mixture of breeds, it may be simpler to extend the approach in this paper and specify
a parametric structure on the relationships of the animals in the base population
and estimate those parameters, instead of developing an appropriate way of adjusting
the markerbased relationship matrix of the genotyped animals across breeds. An additional
issue is the interpretation of the genetic variance parameter
For the simulated dataset 1, the estimate of the relative polygenic weight ω was about 0, as was expected since in the simulation the markers were assumed to
capture all the genetic variation. However, for the simulated dataset 2,
Conclusions
Using the full loglikelihood and adjusting the pedigreebased relationship matrix to be compatible with the markerbased relationship matrix provides a new and interesting approach to handle the issue of compatibility between the two matrices in singlestep genetic evaluation.
Appendix A
The derivation of the proposed adjustment of the pedigreebased relationship matrix is based on assigning priors on the highdimensional parameters ρ and v in the previously described singlestep model. This derivation is presented here.
The model for the marker genotypes, M, given the allele frequencies ρ and variances v, is as in [3],
where v_{j} = v_{j,j} is a parameter and
With these prior distributions, a marginal distribution of M may be obtained by integrating ρ and v. Using wellknown formulas for conditional expectations and covariances, the mean and covariances become
Defining
with
Note that the random labelling of alleles as 1 or 2 is not important because the resulting expression does not depend on how alleles are labelled.
Using
This loglikelihood may also be viewed as the loglikelihood for
Appendix B
The aim here is to present algorithms for computing
Let us partition matrix
Similar to the usual A matrix (see [18]) a decomposition exists
where T is a lower triangular matrix with entries T_{ii} = 1, T_{ik} = (T_{f(i)k} + T_{m(i)k}) / 2 when i has known parents f(i) and m(i) and i > k, and T_{ik} = 0 otherwise (k > i or i has unknown parents), and
The inverse matrix becomes
Here T^{−1} is a lower triangular matrix with ones on the diagonal and the only nonzero elements
being −1 / 2 for offspringparent entries, and matrix
A procedure to obtain the diagonal of
Matrix
Matrix
Therefore, an algorithm (with notation similar to that in Appendix A in Misztal et
al. [19]) consists of first computing r = T^{T}x by solving the sparse system
and finally computing
Competing interests
The author declares that he has no competing interests.
Authors’ contributions
The author has read and approved the final manuscript.
Acknowledgements
The work was performed in a project funded through the Green Development and Demonstration Programme (grant no. 3405110279) by the Danish Ministry of Food, Agriculture and Fisheries, the Pig Research Centre and Aarhus University. The author thanks two anonymous reviewers, an associate editor and editor Jack Dekkers for their comments, and in particular editor Hélène Hayes for her careful and competent language revisions.
References

Legarra A, Aguilar I, Misztal I: A relationship matrix including full pedigree and genomic information.
J Dairy Sci 2009, 92:46564663. PubMed Abstract  Publisher Full Text

Aguilar I, Misztal I, Johnson DL, Legarra A, Tsuruta S, Lawlor TJ: Hot topic: A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluations of Holstein final score.
J Dairy Sci 2010, 93:743752. PubMed Abstract  Publisher Full Text

Christensen OF, Lund MS: Genomic prediction when some animals are not genotyped.
Genet Sel Evol 2010, 42:2. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Vitezica ZG, Aguilar I, Misztal I, Legarra A: Bias in genomic predictions for populations under selection.
Genet Res 2011, 93:357366. Publisher Full Text

Chen CY, Misztal I, Aguilar I, Legarra A, Muir WM: Effect of different genomic relationship matrices on accuracy and scale.
J Anim Sci 2011, 89:26732679. PubMed Abstract  Publisher Full Text

Christensen OF, Madsen P, Nielsen B, Ostersen T, Su G: Singlestep methods for genomic evaluation in pigs.
Animal 2012, 6:15651571. PubMed Abstract  Publisher Full Text

Forni S, Aguilar I, Misztal I: Different genomic relationship matrices for singlestep analysis using phenotypic, pedigree and genomic information.
Genet Sel Evol 2011, 43:1. BioMed Central Full Text

Gilmour AR, Thompson R, Cullis BR: Average information REML: an efficient algorithm for parameter estimation in linear mixed models.
Biometrics 1995, 51:14401450. Publisher Full Text

Gengler N, Mayeres P, Szydlowski M: A simple method to approximate gene content in large pedigree populations: application to the myostation gene in dualpurpose Belgian Blue cattle.
Animal 2007, 1:2128. PubMed Abstract  Publisher Full Text

VanRaden PM: Accounting for inbreeding and crossbreeding in genetic evaluation of large populations.
J Dairy Sci 1992, 75:31363144. Publisher Full Text

Colleau JJ, Sargolzaei M: MIM: an indirect method to assess inbreeding and coancestry in large incomplete pedigrees of selected dairy cattle.
J Anim Breed Genet 2011, 128:163173. PubMed Abstract  Publisher Full Text

Quaas RL: Computing the diagonal elements and inverse of a large numerator relationship matrix.
Biometrics 1976, 32:949953. Publisher Full Text

Colleau JJ: An indirect approach to the extensive calculation of relationship coefficients.
Genet Sel Evol 2002, 34:409421. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

VanRaden PM: Efficient methods to compute genomic predictions.
J Dairy Sci 2008, 91:44144423. PubMed Abstract  Publisher Full Text

Gao H, Christensen OF, Madsen P, Nielsen US, Zhang Y, Lund MS, Su G: Comparison on genomic predictions using GBLUP models and two singlestep blending methods with different relationship matrices in the Nordic Holstein population.
Genet Sel Evol 2012, 44:8. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Meuwissen THE, Luan T, Woolliams JA: The unified approach to the use of genomic and pedigree information in genomic evaluations revisited.
J Anim Breed Genet 2011, 128:429439. PubMed Abstract  Publisher Full Text

Fernando RL, Grossman M: Marker assisted selection using best linear unbiased prediction.
Genet Sel Evol 1989, 21:467477. BioMed Central Full Text

Mrode RA: Linear models for the prediction of animal breeding values. Wallingford: CABI Publishing; 2005.

Misztal I, Legarra A, Aguilar I: Computing procedures for genetic evaluation including phenotypic, full pedigree and genomic information.
J Dairy Sci 2009, 92:46484655. PubMed Abstract  Publisher Full Text

Aguilar I, Misztal I, Legarra A, Tsuruta S: Efficient computation of the genomic relationship matrix and other matrices used in singlestep evaluation.
J Anim Breed Genet 2011, 128:422428. PubMed Abstract  Publisher Full Text