A model with metafounders is similar to an animal model but the scale of the genetic variances is shifted, in particular if there are several metafounders. The model with metafounders accounts for the fact that crossing two animals with different origins may result in a higher variance than either of the origins.
Much of the theory comes from Legarra et al. 2015 (Genetics). The genetic variance component associated with metafounders is called \sigma_{\text{related}}^{2} in the paper, not the same as \sigma_{\text{unrelated}}^{2} used in "regular" evaluations and to make it simple, in that paper we propose to obtain one as a function of the other:
\sigma_{\text{related}}^{2} = \sigma_{\text{unrelated}}^{2}/k
Scale k is lower than 1, resulting in
\sigma_{\text{related}}^{2} > \sigma_{\text{unrelated}}^{2} .
For the case of a single metafounder this is very simple and it is actually exact:
k = 1 - \frac{\gamma}{2}
And it works very well in the little experiences that we made.
For multiple origins we proposed in the paper that, in practice, the population under study is a mixture in equal proportion of all origins (this is approximately true in single-breed sheep ) resulting in the expression
k = 1 + \frac{\overline{\text{diag}\left( \mathbf{\Gamma} \right)}}{2} - \overline{\mathbf{\Gamma}}
Practicalities in blupf90
These two expressions are the default in blupf90test and blupf90iod2, and the way it works is that k is computed internally and genetic variances modified internally. This can be avoided using OPTION lscale false and then the variance declared in the parameter file remains unchanged.
However, in some populations one metafounder (say number 1) is much more represented than the other one. Then it makes sense to use e.g.
k = 1 - \frac{\Gamma_{\left( 1,1 \right)}}{2}
, but this possibility has not been programmed. The workaround for blupf90test and blup90iod2 is to scale genetic variances manually, e.g. if the genetic covariance matrix is
G_{0} = \begin{pmatrix} 10 & 2 \\ 2 & 4 \\ \end{pmatrix}
and the scale factor is, to give an example, k = 0.85 (e.g. see later) then we must do two things:
1. put in the parameter file
\begin{pmatrix} 11.76 & 2.35 \\ 2.35 & 4.71 \\ \end{pmatrix}
, which is equal to
G_{0}/k = \begin{pmatrix} 10 & 2 \\ 2 & 4 \\ \end{pmatrix}/0.85
2. put OPTION lscale false
An example. Assume a breed where most all animals come from in-country (original) and a few are from another country. We have
\Gamma = \begin{pmatrix} 0.49 & 0.38 \\ 0.38 & 0.57 \\ \end{pmatrix}
and using the default scaling that we programmed in blupf90 results in k = 0.81 . However, if we consider that all individuals are in practice "original" we haven k = 1 - \frac{0.49}{2} = 0.76 . So I would use this number instead.
A possibility is a weighted mean of elements of \Gamma . But the problem is, first, that we need to weight the values in the diagonal and in the off-diagonal of \mathbf{\Gamma} using and argument like in page 461 of the metafounders paper (right column, top), which assumes that all metafounders contribute equally (the algebra needs to be derived but it should be easy). Still this is not completely correct, because the number of founders related to each metafounder is not necessarily equal to their contribution to the overall population. It could happen that "external" are only 1% of the founders but they contribute to 10% of the overall populations.
How to obtain scaling factor in complex crosses
Now we put a sketch of a more complete theory. If we have a whole relationship matrix across animals (excluding metafounders) say D, then the variance of an animal sampled at random is on expectation s^{2} = \overline{\text{diag}\left( D \right)} - \overline{D} (Legarra, Theor. Pop. Biol. 2016).
In the regular relationship matrix \mathbf{A} , this variance is
s^{2} = \left( \overline{\text{diag}\left( A \right)} - \overline{A} \right)\sigma_{\text{unrelated}}^{2} = \left( 1 + \overline{F} - \overline{A} \right)\sigma_{\text{unrelated}}^{2} \approx \left( 1 - \overline{F} \right)\sigma_{\text{unrelated}}^{2}
The last approximation is because inbreeding is half parents' relationships and assumes slow increase of inbreeding.
In the relationship matrix \mathbf{A}^{\left( \mathbf{\Gamma} \right)} this is
s^{2} = \left( \overline{\text{diag}\left( \mathbf{A}^{\left( \mathbf{\Gamma} \right)} \right)} - \overline{\mathbf{A}^{\left( \mathbf{\Gamma} \right)}} \right)\sigma_{\text{related}}^{2} = \left( 1 + \overline{F^{\left( \Gamma \right)}} - \overline{\mathbf{A}^{\left( \mathbf{\Gamma} \right)}} \right)\sigma_{\text{related}}^{2} \approx \left( 1 - \overline{F^{\left( \Gamma \right)}} \right)\sigma_{\text{related}}^{2}
The approximation here \overline{F^{\left( \Gamma \right)}} \approx \overline{A^{\left( \Gamma \right)}}/2 is less clear !!
So then if we want to equate both s^{2} and we use both approximations. We force
\left( 1 - \overline{F^{\left( \Gamma \right)}} \right)\sigma_{\text{related}}^{2} = \left( 1 - \overline{F} \right)\sigma_{\text{unrelated}}^{2}
which results in
\frac{\sigma_{\text{unrelated}}^{2}}{\sigma_{\text{related}}^{2}} = \frac{\left( 1 - \overline{F^{\left( \Gamma \right)}} \right)}{\left( 1 - \overline{F} \right)} and
therefore
k = \frac{\left( 1 - \overline{F^{\left( \Gamma \right)}} \right)}{\left( 1 - \overline{F} \right)} .
We could equally use the less approximate expressions involving \overline{\mathbf{A}} .
Potentially, this can be programmed in blupf90: the inbreeding F^{\left( \Gamma \right)} is computed internally, and we could also compute F by the same subroutine by calling it twice: one with \mathbf{\Gamma} = \mathbf{0} and another with the true \mathbf{\Gamma} .
I show now that the equivalence that I propose here is approximately the same as k = 1 - \frac{\gamma}{2} for a single metafounder. We have that
A^{\left( \Gamma \right)} = A\left( 1 - \frac{\gamma}{2} \right) + \gamma
which results in
F^{\left( \Gamma \right)} = F + \frac{\gamma}{2}
and therefore
k = \frac{\left( 1 - \overline{F^{\left( \Gamma \right)}} \right)}{\left( 1 - \overline{F} \right)} = \frac{1 - \overline{F} - \frac{\gamma}{2}}{1 - \overline{F}} = 1 - \frac{\gamma}{2\left( 1 - \overline{F} \right)}
Which is almost the same as k = 1 - \frac{\gamma}{2} .
All these ideas are *approximations* because at the end, the "classical" A and the A^{\Gamma} are too different models.
I think that
k = 1 + \frac{\overline{\text{diag}\left( \mathbf{\Gamma} \right)}}{2} - \overline{\mathbf{\Gamma}}
is good for UPGs in cattle/sheep for pure breed if they are well defined, and I believe that
k = \frac{\left( 1 - \overline{F^{\left( \Gamma \right)}} \right)}{\left( 1 - \overline{F} \right)}
should work well in cases of strong unbalance among UPGs. But this needs to be verified with real data.
