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

参考：https://cran./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)

gid	y
1	-52.34296
2	-10.69953
3	21.22609
4	64.13197
5	57.48252
6	71.74167

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

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

	1	2	3	4	5
1	1.97080231	0.01052722	-0.04917231	-0.04563456	-0.06999679
2	0.01052722	1.99001872	0.04887964	0.07653841	-0.03626756
3	-0.04917231	0.04887964	2.01534579	-0.03944350	-0.09596709
4	-0.04563456	0.07653841	-0.03944350	1.99830027	0.02817576
5	-0.06999679	-0.03626756	-0.09596709	0.02817576	1.99781785

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)$varcompASReml: 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

	gamma	component	std.error	z.ratio	constraint
giv(gid).giv	1.12492	521.7146	176.7487	2.951731	Positive
R!variance	1.00000	463.7794	309.5252	1.498358	Positive

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

	Estimate	SE
h2	0.529394	0.2450095

pre <- predict(mod,"gid")$prediction$pvalsASReml: 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灵活（这个结论有推广之嫌）