【原】rrBLUP軟件包和asreml-r計算GBLUP比較

育種數(shù)據(jù)分析 2021-11-18

展開全文

參考：https://cran./web/packages/rrBLUP/rrBLUP.pdf

思路：

1，模擬200個體，每個個體1000SNP的基因型數(shù)據(jù)

2，模擬一個表型數(shù)據(jù)，首先模擬每個SNP的效應值，然后根據(jù)遺傳力，模擬表型值

3，使用rrBLUP計算A矩陣（G矩陣）

4，使用rrBLUP計算GBLUP，因為有真值（True breeding value），計算準確性

5，使用asreml軟件，計算GBLUP，這里G矩陣的對角線加0.01，防止奇異。同時使用attr將G的ID賦予rowNames，用于asreml的定義。

6，rrBLUP和asreml結(jié)果對比。

總體來說，這篇微信文是代碼的分享，然后炫耀我會使用這么多軟件。中午吃飯時，討論為什么要戒煙以及戒游戲，得到這個結(jié)論：人生的樂趣這么少，為什么還要戒。。。

1，生成模擬基因組數(shù)據(jù)

生成200*1000個SNP數(shù)據(jù)，其中有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，生成模擬表型數(shù)據(jù)，遺傳力為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)$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

	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$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

結(jié)論

rrBLUP和asreml結(jié)果基本一致（遺傳力定了，G矩陣也定了，求解方程結(jié)果還不一樣的話，are you kiding???)
注意，由于模擬數(shù)據(jù)都在變化，結(jié)果不一定和我的一模一樣（這個結(jié)論不專業(yè)，我已經(jīng)設置了set.seed）
rrBLUP定義固定因子以及殘差矩陣結(jié)構(gòu)沒有asreml靈活（這個結(jié)論有推廣之嫌）

轉(zhuǎn)藏 分享

QQ空間 QQ好友新浪微博微信

獻花（0） +1

來自：育種數(shù)據(jù)分析 > 《待分類》

舉報/認領

0條評論

發(fā)表

請遵守用戶評論公約

類似文章 更多

育種數(shù)據(jù)分析

關注對話

TA的最新館藏

如逆水行舟不進則退
單倍型分析：個體所對應的單倍型是？
如何利用系譜或者SNP數(shù)據(jù)劃分家系？
如何利用DeepSeek的API搭建本地知識庫
群體遺傳三劍客第三篇：megacc和ggtree進化樹分析
如何計算群體中的單倍型頻率

喜歡該文的人也喜歡更多

熱門閱讀換一換

乡下人产国偷v产偷v自拍,国产午夜片在线观看,婷婷成人亚洲综合国产麻豆,久久综合给合久久狠狠狠9

【原】rrBLUP軟件包和asreml-r計算GBLUP比較

1，生成模擬基因組數(shù)據(jù)

2，生成模擬表型數(shù)據(jù)，遺傳力為0.5

3，使用rrBLUP計算計算G矩陣

4，使用rrBLUP計算GBLUP，并比較準確性