rrBLUP软件包和asreml-r计算GBLUP比较

这篇微信文，主要介绍GBLUP计算的两个软件rrBLUP和asreml-r的计算比较，通过模拟的数据，比较两者的结果。

参考：https://cran.r-project.org/web/packages/rrBLUP/rrBLUP.pdf

思路：

1，模拟200个体，每个个体1000SNP的基因型数据

2，模拟一个表型数据，首先模拟每个SNP的效应值，然后根据遗传力，模拟表型值

3，使用rrBLUP计算A矩阵（G矩阵）

4，使用rrBLUP计算GBLUP，因为有真值（True breeding value），计算准确性

5，使用asreml软件，计算GBLUP，这里G矩阵的对角线加0.01，防止奇异。同时使用attr将G的ID赋予rowNames，用于asreml的定义。

6，rrBLUP和asreml结果对比。

总体来说，这篇微信文是代码的分享，然后炫耀我会使用这么多软件。中午吃饭时，讨论为什么要戒烟以及戒游戏，得到这个结论：人生的乐趣这么少，为什么还要戒。。。

1，生成模拟基因组数据

生成200*1000个SNP数据，其中有200个个体，每个个体有1000个SNP

rm(list=ls())
set.seed(100)
M <- matrix(rep(0,200*1000),200,1000)
for (i in 1:200) {
  M[i,] <- ifelse(runif(1000)<0.5,-1,1)
}
rownames(M) <- 1:200

2，生成模拟表型数据，遗传力为0.5

#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)))
data <- data.frame(y=y,gid=1:200)
data$gid <- as.factor(data$gid)
data<- data[,c(2,1)]
head(data)

3，使用rrBLUP计算计算G矩阵

if(!requireNamespace("rrBLUP")) install.packages("rrBLUP")
library(rrBLUP) 
rownames(M) <- 1:200
A <- A.mat(M)
A[1:5,1:5]

4，使用rrBLUP计算GBLUP，并比较准确性

ans <- kin.blup(data=data,geno="gid",pheno="y",K=A)
accuracy <- cor(g,ans$g)
accuracy
plot(g,ans$g)

0.73764143220772

5，使用asreml-r计算GBLUP

if (!requireNamespace("learnasreml")) devtools::install_github("dengfei2013/learnasreml")
library(asreml)
library(learnasreml)

diag(A) = diag(A) + 0.01

ginv <- write_relation_matrix(A,type = "ginv")
attr(ginv,"rowNames") <- 1:200
mod <- asreml(y ~ 1, random=~ giv(gid),
              ginverse = list(gid=ginv),data=data)
summary(mod)$varcomp

ASReml: Fri Nov 09 21:10:14 2018

     LogLik         S2      DF      wall     cpu
   -831.0248   1261.3325   199  21:10:14     0.0
   -829.6764   1027.8914   199  21:10:14     0.0
   -828.6137    764.9393   199  21:10:14     0.0
   -828.1743    535.8671   199  21:10:14     0.0
   -828.1459    468.9919   199  21:10:14     0.0
   -828.1457    463.5712   199  21:10:14     0.0
   -828.1457    463.7794   199  21:10:14     0.0

Finished on: Fri Nov 09 21:10:14 2018

LogLikelihood Converged 

pin(mod,h2 ~V1/(V1+V2))

pre <- predict(mod,"gid")$prediction$pvals

ASReml: Fri Nov 09 21:10:18 2018

     LogLik         S2      DF      wall     cpu
   -828.1457    463.7690   199  21:10:18     0.0
   -828.1457    463.7692   199  21:10:18     0.0
   -828.1457    463.7694   199  21:10:18     0.0
   -828.1457    463.7695   199  21:10:18     0.0

Finished on: Fri Nov 09 21:10:18 2018

LogLikelihood Converged 

计算GBLUP和TBV的相关性

cor(g,pre$predicted.value)
plot(g,pre$predicted.value)

0.736187643745443

计算asreml和rrBLUP的GBLUP相关性

cor(pre$predicted.value,ans$g)
plot(pre$predicted.value,ans$g)

0.999913776260597

结论

• rrBLUP和asreml结果基本一致（遗传力定了，G矩阵也定了，求解方程结果还不一样的话，are you kiding???)
  

• 注意，由于模拟数据都在变化，结果不一定和我的一模一样（这个结论不专业，我已经设置了set.seed）

• rrBLUP定义固定因子以及残差矩阵结构没有asreml灵活（这个结论有推广之嫌）