（1）自交系预测
模型：
y = Xβ + Zu + ε代码：
# sommer的GS预测
library(sommer) # 加载sommer包，第一次安装请先运行install.packages("sommer")
data(DT_wheat)
DT <- DT_wheat
GT <- GT_wheat
colnames(DT) <- paste0("X",1:ncol(DT)) # 列名为X1~X4
DT <- as.data.frame(DT) # 转换为数据框
DT$id <- as.factor(rownames(DT)) # 把行号作为列并作为因子
# select environment 1
rownames(GT) <- rownames(DT) # 给基因型赋值相同的行号（名称）
K <- A.mat(GT) # additive relationship matrix
colnames(K) <- rownames(K) <- rownames(DT) # K的行列相同，都赋值名称
# GBLUP pedigree-based approach
set.seed(12345)
y.trn <- DT
vv <- sample(rownames(DT),round(nrow(DT)/5)) # 抽取20%作为验证群体（被预测）
y.trn[vv,"X1"] <- NA
head(y.trn)
## GBLUP
ans <- mmer(X1~1,
random=~vs(id,Gu=K),
rcov=~units,
data=y.trn, verbose = FALSE) # kinship based
ans$U$`u:id`$X1 <- as.data.frame(ans$U$`u:id`$X1)
rownames(ans$U$`u:id`$X1) <- gsub("id","",rownames(ans$U$`u:id`$X1))
cor(ans$U$`u:id`$X1[vv,],DT[vv,"X1"], use="complete")
## rrBLUP
ans2 <- mmer(X1~1,
random=~vs(list(GT)),
rcov=~units,
data=y.trn, verbose = FALSE) # kinship based
u <- GT %*% as.matrix(ans2$U$`u:GT`$X1) # BLUPs for individuals
rownames(u) <- rownames(GT)
cor(u[vv,],DT[vv,"X1"]) # same correlation
数据使用的是CIMMYT的全球小麦项目组提供的数据，共有599个个体。y.trn表示籽粒产量，每一列表示一个环境，共有4个环境，例子中只计算了80%预测20%的结果。
（2）杂交种的预测
杂交种的预测使用GCA效应和SCA效应。数据来自1254年测定的玉米单交试验。123个dent和86个flint，用测定的1254材料预测未测定的9324个材料。例子中只计算了第一个性状GY。
模型：
y = Xβ + Zu1 + Zu2 + ZuS + ε代码：
library(sommer)
data(DT_technow)
DT <- DT_technow # 表型数据
Md <- Md_technow # dent马齿基因型矩阵
Mf <- Mf_technow # flint硬粒基因型矩阵
Ad <- Ad_technow # dent的加性效应矩阵
Af <- Af_technow # flint的加性效应矩阵
y.trn <- DT
vv1 <- which(!is.na(DT$GY)) # 获得有数据的材料行号
vv2 <- sample(vv1, 100) # 在有数据的材料行号中随机选择100个材料
y.trn[vv2,"GY"] <- NA # 将抽取的100个材料设为NA
anss2 <- mmer(GY~1, # 固定效应
random=~vs(dent,Gu=Ad) + vs(flint,Gu=Af), # 随机效应，并指定方差协方差矩阵
rcov=~units,
data=y.trn, verbose = FALSE)
summary(anss2)$varcomp
zu1 <- model.matrix(~dent-1,y.trn) %*% anss2$U$`u:dent`$GY
zu2 <- model.matrix(~flint-1,y.trn) %*% anss2$U$`u:flint`$GY
u <- zu1+zu2+anss2$Beta[1,"Estimate"]
cor(u[vv2,], DT$GY[vv2])
VarComp VarCompSE Zratio Constraint
u:dent.GY-GY 16.06423 2.5737578 6.241548 Positive
u:flint.GY-GY 11.42070 2.1591718 5.289390 Positive
units.GY-GY 16.81801 0.7689509 21.871368 Positive
> cor(u[vv2,], DT$GY[vv2])
[1] 0.922328预测时，所有有数据的个体作为建模群体，这里随机减去100个个体的测定值用于验证，即建模群体是有数据的个体-100。zu1是模型中的Zu1，Z是model.matrix(~dent-1,y.trn)，即dent的设计矩阵；u1是anss2$U$u:dent$GY是GY性状dent的BLUP值。相应的，zu2是flint的设计矩阵×GY性状flint的BLUP。
anns2$Beta[1,”Estimate”]是固定效应的BLUP。
最后，100个个体的测定值和这100个个体的BLUP值的相关系数作为预测结果。
上面的例子只使用了GCA进行预测，因为SCA在这里面没有对预测精度起到多大的影响，相反，由于使用了很大的矩阵，导致计算时间明显增加。
发表评论