ASR012. LMM with pedigree information - Atlantic salmon

The complete script for this example can be downloaded here:

Dataset

In this example, we will use the D008 and the D008ped datasets, which include the phenotypic data and pedigree file, respectively. A few lines of each of these data frames are shown below:

d008

indiv	mean_gill_score	amoebic_load
Fish_00001	3.00	28.72
Fish_00002	3.50	26.47
Fish_00003	2.00	32.78
Fish_00004	3.25	26.90
Fish_00005	3.00	30.34
Fish_00007	2.50	31.36

d008ped

indiv	sire	dam
Sire1	0	0
Sire2	0	0
Sire3	0	0
Sire4	0	0
Sire5	0	0
Sire6	0	0

Model

In this example, we will fit a LMM with pedigree information (i.e., Animal Model) based on the following equation: \[ y = \mu + gen+ e\ \]

where,

\(y\) is the amoebic load,

\(\mu\) is the overall mean,

\(gen\) is the random effect of genotype, with \(gen \sim \mathcal{N}(0,\,\sigma^{2}_{g})\),

\(e\) is the random residual effect, with \(e \sim \mathcal{N}(0,\,\sigma^{2}_{e})\).

The factor \(gen\) is associated with the numerator relationship matrix (NRM) \(\mathbf{A}\).

Now, let’s take a look at how to write the model with ASReml-R. Note that before fitting the model, the inverse of the NRM needs to be created using the pedigree information (see ASR013 for more details), and the column indiv from the phenotypic dataset needs to be set as factor.

asr012 <- asreml(
  fixed = amoebic_load ~ 1,
  random = ~vm(indiv, ainv),
  residual = ~units,
  data = d008
)

Note that the above model is including the inverse of the numerator relationship matrix by the use of the expression vm() that indicates variance model and that links a factor to the inverse matrix.

Exploring output

The residuals diagnostic plots are:

plot(asr012)

The variance components are:

summary(asr012)$varcomp

                component std.error   z.ratio bound %ch
vm(indiv, ainv)  3.805643 0.8391101  4.535333     P   0
units!R          6.721801 0.5879412 11.432776     P   0

The heritability is calculated as follow:

vpredict(asr012, h2 ~ V1/(V1+V2))

    Estimate         SE
h2 0.3614974 0.06850024

This result reports a heritability estimate of 0.36, which indicates that the studied trait is under important additive genetic control.

The overall mean, which is the only BLUE effect we have can be requested with:

summary(asr012, coef = TRUE)$coef.fixed

            solution std error  z.ratio
(Intercept) 31.61313 0.1635798 193.2582

The additive genetic effects (BLUP/breeding values) are:

BLUP <- as.data.frame(summary(asr012, coef = TRUE)$coef.random)
head(BLUP)

                        solution std.error    z.ratio
vm(indiv, ainv)_Sire1  0.1768296 0.8396683  0.2105945
vm(indiv, ainv)_Sire2 -1.1805100 1.2934350 -0.9126937
vm(indiv, ainv)_Sire3 -2.3679553 0.9298205 -2.5466801
vm(indiv, ainv)_Sire4  1.1015471 0.8689376  1.2676942
vm(indiv, ainv)_Sire5  0.7081991 0.8875566  0.7979200
vm(indiv, ainv)_Sire6 -2.9925001 1.2173379 -2.4582328

We can also get the predictions (adjusted means) as follow:

preds <- predict(asr012, classify = 'vm(indiv,ainv)')$pvals
head(preds)

  indiv predicted.value std.error    status
1 Sire1        31.78999 0.8121392 Estimable
2 Sire2        30.43264 1.2852457 Estimable
3 Sire3        29.24521 0.9072839 Estimable
4 Sire4        32.71468 0.8463674 Estimable
5 Sire5        32.32133 0.8664946 Estimable
6 Sire6        28.62068 1.2051464 Estimable

In this case, the predictions are simply corresponding to the sum of the overall mean and the breeding value, and they can be used to rank and select future parents.