单性状动物模型矩阵形式计算BLUP值

动物育种
邓飞    2020-09-08    404

  公众号叫:“育种数据分析之放飞自我”,但是竟然没有一个相关的专辑,因此,决定将之前写过的关于育种数据分析相关的博文,重新整理一下,包括软件包有些更新后语法变化了,有些软件包不在维护有了更好的替代品。经常有人给我发过来很久以前的博客,告诉我运行报错,然后我发现有些排版时错误的,有些R包已经不再维护退出CRAN了,因此重新整理一下,确保发出去的东西能够运行,整理一下之前的博文。

  专业性的博文增加生活的段子,其乐无穷。

  段子:

  同事去大骨头, 吃过之后进行了抽奖, 然后向我炫耀, 看我抽到了一个小猪佩奇。 照片如下:

  我看过之后, 冷静的说, 首先, 我们去吃饭, 奖品都是送的, 不需要抽奖的。

  其次, 这是乔治, 不是佩奇。

  1, 数据

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

  下面是生成数据的代码

  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

  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, 计算亲缘关系逆矩阵

  library(nadiv)

  提取系谱信息

  ped <- dat[,2:4]

  ped

  ID    Sire    Dam

  1    0    0

  2    0    0

  3    1    0

  4    1    2

  5    3    2

  计算亲缘关系逆矩阵

  pped = prepPed(ped)

  pped

  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

  Ainv = makeAinv(pped)$Ainv

  Ainv

  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

  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构建矩阵, 比较方便

  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

  构建单元矩阵

  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的矩阵

  y <- as。matrix(dat$weight)

  y

  140

  152

  135

  143

  160

  混合线性方程组

  XpZ <- crossprod(X,Z);XpZ

  Chang1    1    1    1    0    0

  Chang2    0    0    0    1    1

  X’X

  XpX <- crossprod(X) ;XpX

  Chang1    Chang2

  Chang1    3    0

  Chang2    0    2

  Z’X

  ZpX <- crossprod(Z,X);ZpX

  Chang1    Chang2

  1    0

  1    0

  1    0

  0    1

  0    1

  Z’Z

  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

  Xpy <- crossprod(X,y);Xpy

  Chang1    427

  Chang2    303

  Z’y

  Zpy <- crossprod(Z,y);Zpy

  140

  152

  135

  143

  160

  K

  K <- 2;K

  2

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

  LHS

  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左手矩阵和上图结果一致。

  RHS <- rbind(Xpy,Zpy)

  RHS

图片.png  求解BLUP值

  solve(LHS)%*%RHS

  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

图片.png  可以看到, 结果虽然结果不一致, 但是PPT里面的结果是错误的…

  所以说,PPT里面的内容也不一定是正确的,现场演示之后,发现PPT里面的结果是错误的,这该如何圆场???现场翻车记!!!

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

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

  Ainv

  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

本文来自微信公众号【育种数据分析之放飞自我】公众号ID:R-breeding;未经许可谢绝二次转载至其他网站。


本文为专栏作者授权科易网发表,版权归原作者所有。文章系作者个人观点,不代表科易网立场,转载请联系原作者。如有任何疑问,请联系ky@1633.com。

熟练使用统计软件和遗传评估软件,负责开发了具有数据管理、常规遗传评估、MAS、GWAS和GS的育种云平台,具有丰富的全基因组选择分析及平台建设经验。 [查看更多]

13164阅读数 9238访问数 924订阅数
热门观点

向邓飞提问