N3 | 动物模型计算BLUP值

1, 数据

代码语言：javascript复制

这次使用一个PPT里面的数据, 用R语言演示一下如何做BLUP值计算.

下面是生成数据的代码

代码语言：javascript复制

Chang <- c(1,1,1,2,2)
ID <- c(1,2,3,4,5)
Sire <- c(0,0,1,1,3)
Dam <- c(0,0,0,2,2)
weight <- c(140,152,135,143,160)
dat <- data.frame(Chang,ID,Sire,Dam,weight)

代码语言：javascript复制

dat

Chang	ID	Sire	Dam	weight
1	1	0	0	140
1	2	0	0	152
1	3	1	0	135
2	4	1	2	143
2	5	3	2	160

2, 计算亲缘关系逆矩阵

代码语言：javascript复制

library(nadiv)

提取系谱信息

代码语言：javascript复制

ped <- dat[,2:4]

ped

ID	Sire	Dam
1	0	0
2	0	0
3	1	0
4	1	2
5	3	2

计算亲缘关系逆矩阵

代码语言：javascript复制

pped = prepPed(ped)

pped

代码语言：javascript复制

Warning message in prepPed(ped):
"Zero in the dam column interpreted as a missing parent"Warning message in prepPed(ped):
"Zero in the sire column interpreted as a missing parent"

首先, 将系谱进行一下转换, 使用nadiv的prepPed函数, 预处理. 它会自动不齐没有亲本的个体, 变为NA.

ID	Sire	Dam
1	NA	NA
2	NA	NA
3	1	NA
4	1	2
5	3	2

如果是计算逆矩阵的矩阵形式, 可以使用makeAinv(pped)$Ainv

代码语言：javascript复制

Ainv = makeAinv(pped)$Ainv

Ainv

代码语言：javascript复制

5 x 5 sparse Matrix of class "dgCMatrix"

1  1.8333333  0.5 -0.6666667 -1  .
2  0.5000000  2.0  0.5000000 -1 -1
3 -0.6666667  0.5  1.8333333  . -1
4 -1.0000000 -1.0  .          2  .
5  .         -1.0 -1.0000000  .  2

如果是计算逆矩阵的行列形式, 可以使用makeAinv(pped)$listAinv

代码语言：javascript复制

makeAinv(pped)$listAinv

	row	column	Ainv
1	1	1	1.8333333
5	2	1	0.5000000
6	2	2	2.0000000
10	3	1	-0.6666667
11	3	2	0.5000000
12	3	3	1.8333333
14	4	1	-1.0000000
15	4	2	-1.0000000
16	4	4	2.0000000
17	5	2	-1.0000000
18	5	3	-1.0000000
19	5	5	2.0000000

教科书的结果, 两者一样

3, 构建模型

$$ y = Xb Zu e $$

构建固定因子矩阵

这里使用函数model.matrix构建矩阵, 比较方便

代码语言：javascript复制

for(i in 1:4) dat[,i] <- as.factor(dat[,i])

X <- model.matrix(~Chang-1,dat)

X

	Chang1	Chang2
1	1	0
2	1	0
3	1	0
4	0	1
5	0	1

构建单元矩阵

代码语言：javascript复制

Z <- diag(length(unique(dat$ID)))
Z

1	0	0	0	0
0	1	0	0	0
0	0	1	0	0
0	0	0	1	0
0	0	0	0	1

构建y的矩阵

代码语言：javascript复制

y <- as.matrix(dat$weight)
y

140
152
135
143
160

混合线性方程组

代码语言：javascript复制

XpZ <- crossprod(X,Z);XpZ

Chang1	1	1	1	0	0
Chang2	0	0	0	1	1

X’X

代码语言：javascript复制

XpX <- crossprod(X) ;XpX

	Chang1	Chang2
Chang1	3	0
Chang2	0	2

Z’X

代码语言：javascript复制

ZpX <- crossprod(Z,X);ZpX

Chang1	Chang2
1	0
1	0
1	0
0	1
0	1

Z’Z

代码语言：javascript复制

ZpZ <- crossprod(Z);ZpZ

1	0	0	0	0
0	1	0	0	0
0	0	1	0	0
0	0	0	1	0
0	0	0	0	1

X’y

代码语言：javascript复制

Xpy <- crossprod(X,y);Xpy

Chang1	427
Chang2	303

Z’y

代码语言：javascript复制

Zpy <- crossprod(Z,y);Zpy

140
152
135
143
160

K

代码语言：javascript复制

K <- 2;K

2

代码语言：javascript复制

LHS <- rbind(cbind(XpX,XpZ),cbind(ZpX,ZpZ Ainv*K))
LHS

代码语言：javascript复制

7 x 7 sparse Matrix of class "dgCMatrix"
       Chang1 Chang2                            
Chang1      3      .  1.000000  1  1.000000  .  .
Chang2      .      2  .         .  .         1  1
1           1      .  4.666667  1 -1.333333 -2  .
2           1      .  1.000000  5  1.000000 -2 -2
3           1      . -1.333333  1  4.666667  . -2
4           .      1 -2.000000 -2  .         5  .
5           .      1  .        -2 -2.000000  .  5

可以看到, 里面的LHS左手矩阵和上图结果一致.

代码语言：javascript复制

RHS <- rbind(Xpy,Zpy)
RHS

求解BLUP值

代码语言：javascript复制

solve(LHS)%*%RHS

代码语言：javascript复制

7 x 1 Matrix of class "dgeMatrix"
           [,1]
[1,] 142.842105
[2,] 151.118421
[3,]  -2.462551
[4,]   3.052632
[5,]  -2.116397
[6,]  -1.387652
[7,]   2.150810

最后，为了方便大家重演，我将相关的生产数据代码，运行代码汇总如下：

代码语言：javascript复制

Chang <- c(1,1,1,2,2)ID <- c(1,2,3,4,5)Sire <- c(0,0,1,1,3)Dam <- c(0,0,0,2,2)weight <- c(140,152,135,143,160)dat <- data.frame(Chang,ID,Sire,Dam,weight)dat
library(nadiv)ped <- dat[,2:4]ped
pped = prepPed(ped)pped
Ainv = makeAinv(pped)$AinvAinv
makeAinv(pped)$listAinv
for(i in 1:4) dat[,i] <- as.factor(dat[,i])X <- model.matrix(~Chang-1,dat)X
Z <- diag(length(unique(dat$ID)))Z
y <- as.matrix(dat$weight)y
XpZ <- crossprod(X,Z);XpZ
XpX <- crossprod(X) ;XpX
ZpX <- crossprod(Z,X);ZpX
ZpZ <- crossprod(Z);ZpZ
Xpy <- crossprod(X,y);Xpy
Zpy <- crossprod(Z,y);Zpy
K <- 2;K
LHS <- rbind(cbind(XpX,XpZ),cbind(ZpX,ZpZ Ainv*K))LHS
RHS <- rbind(Xpy,Zpy)RHS
solve(LHS)%*%RHS

na

0 人点赞