Solved – Calculated breeding values using markers using animal model in R

blupmixed modelrridge regression

Animal model (frequently used in animal science and sometime in human or plants) is mixed model with:

$y$ = $X$$b$ + $Z$$u$ + $e$

y is observed values for any quantitative variable, $Xb$ is fixed effects terms while $Zu$ is random effects, and $e$ is residual. Covariance between terms of $u$ can be defined as:
$G$ = $A$$\sigma$$^2$$u$ . $A$ is additive relationship matrix whereas $\sigma$$^2$$u$ is additive genetic variance. $A$ can estimated from pedigree or from markers most commonly as $M$$M'$ where $M$ is marker scores (-1,0,1) x observation matrix.

enter image description here

If there are no fixed effects in the model above model can be re-written as:

$y$ = $Z$$u$ + $e$

As $A$ is estimated from large number of markers with multi-colinarity it is logical to shrink the $u$ toward 0 using ridge shrinkage. From the above information we can calculate Best Linear Predictors (BLUPs).

enter image description here

Here for individual i, the BLUP estimate is $yi$, which is summation of $ai$ is additive genetic effect and $ei$ is residual.

Is there any R software package that can do this ? Examples appreciated.

Best Answer

The mixed.solve or kin.blup functions from R package rrBLUP can do the type of analysis you are talking. In addition to reference manual in CRAN, the official website consists of vignettes. Here is an example from the package:

require(rrBLUP)

#random population of 200 individuals/lines with 1000 markers
M <- matrix(rep(0,200*1000),200,1000)
for (i in 1:200) {
  M[i,] <- ifelse(runif(1000)<0.5,-1,1)
}

#random phenotypes
u <- rnorm(1000)
g <- as.vector(crossprod(t(M),u))
h2 <- 0.5  #heritability
y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))


#predict breeding values, here we are using MM' as relatedness measure 
# using function A.mat 
ans <- mixed.solve(y,K=A.mat(M), method="REML", bounds=c(1e-09, 1e+09))
# y is vector of y observations, K is kinship or relatedness matrix
# the defualt method is REML, the bounds specifies the upper and lower bound for 
# ridge parameter 

# BLUP estimation 
ans$u 

# if you are interested in prediction marker effects 
ans <- mixed.solve(y,Z=M)  #By default K = I

# MM' matrix, measure of relatedness, use function A.mat  
A <- A.mat(M)

# The kin.blup is wrapper function for mixed.solve, can also be used.
Related Question