关注我们：
【GS专栏】8-全基因组选择实战之RRBLUP
知识/政策规划
2021-03-15
1.rrBLUP包简介
R 包rrBLUP能方便快速地实现RRBLUP(Ridge Regression Best Linear Unbiased Prediction)模型，该包主要有A.mat和mixed.solve两个核心函数。
A.mat函数用来过滤和填充基因型，返回加性效应关系矩阵(即Kinship)。
A.mat(
X, #编码为{-1,0,1}的标记矩阵，可包含小数点(表明已填充过)和NA
min.MAF=NULL, #最小等位基因频率
max.missing=NULL, #最大缺失值比例
impute.method="mean", #填充方法，mean/EM
tol=0.02, #EM算法的收敛标准
n.core=1, #EM算法的并行核心数(仅unix)
shrink=FALSE, #收缩估计
return.imputed=FALSE #返回填充标记矩阵
)
mixed.solve函数求解模型参数。
mixed.solve(
y, #观测值
Z=NULL, #随机效应矩阵
K=NULL, #随机效应协变量矩阵
X=NULL, #固定效应矩阵
method="REML", #最大似然估计方法，ML/REML
bounds=c(1e-09, 1e+09), #岭回归参数两元素边界
SE=FALSE, #是否计算标准误
return.Hinv=FALSE #是否H取逆，一般在GWAS中用，忽略之
)
mixed.solve函数返回值(列表)：
• Vu：σ^2_u的估计值
• Ve：σ^2_e的估计值
• beta：BLUE(β)，即训练群表型均值
• u：BLUP(u)，即预测值
• LL：maximized log-likelihood(ML/REML值)
如果设了标准误参数，则会返回：
• beta.SE：β标准误
• u.SE：u^*-u
2. 实操
安装和加载包
if(!require(rrBLUP))install.packages("rrBLUP")
library(rrBLUP)
读入表型和基因型数据
示例数据来源：
表型 traits.txt[1]
基因型 snp.txt[2]
• 查看表型数据，共96个品种3个性状。
> Pheno <- as.matrix(read.table(file ="traits.txt", header=TRUE))
> dim(Pheno)
[1] 96 3
> head(Pheno)
grain_yield pht_height Heading_Date
[1,] 4869 78 57.5
[2,] 4621 74 56.0
[3,] 4557 66 58.0
[4,] 5484 77 59.0
[5,] 4641 69 55.5
[6,] 4516 76 58.5
• 查看基因型数据，96个品种共1178个SNP标记(已编码)。
> Markers <- as.matrix(read.table(file="snp.txt"), header=F)
> dim(Markers)
[1] 96 1178
> Markers[1:10,1:10]
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
[1,] -1 1 1 -1 -1 1 -1 1 1 1
[2,] -1 1 1 -1 -1 1 -1 1 1 1
[3,] -1 1 1 -1 -1 1 -1 -1 -1 1
[4,] -1 1 1 -1 -1 1 1 1 1 1
[5,] 1 1 1 1 -1 1 -1 -1 -1 1
[6,] -1 1 1 -1 -1 1 -1 -1 -1 1
[7,] -1 1 1 -1 -1 1 0 1 1 1
[8,] -1 1 1 -1 -1 1 -1 1 1 1
[9,] -1 1 1 -1 -1 1 -1 1 1 1
[10,] -1 1 1 -1 -1 1 -1 1 1 1
数据清洗
所谓数据清洗，是对基因型和表型数据进行质控，过滤掉那些不可信的数据，以及填充缺失的数据。在正式数据分析前，一般会对整个数据做一个评估，按照一定标准进行预处理。
在本示例数据中，暂未对表型数据做处理，我们使用A.mat函数来对基因型数据进行过滤和填充。
#remove markers with more than 50% missing data
> impute = A.mat(Markers,max.missing=0.5,impute.method="mean",return.imputed=T)
> Markers_impute2 = impute$imputed
> dim(Markers)
[1] 96 1178
> dim(Markers_impute2)
[1] 96 1176
按50%缺失值过滤，并按均值填充后，剩下1176个标记。
划分数据集
建模之前，需要对群体定义训练集和测试集。如这里我们将96个品种的60%划为训练集，40%划为测试集。
#training-60% validation-40%
train = as.matrix(sample(1:96, 58))
test <- setdiff(1:96,train)
# train set
Pheno_train = Pheno[train,]
m_train = Markers_impute2[train,]
# valid set
Pheno_valid = Pheno[test,]
m_valid = Markers_impute2[test,]
建模
我们以三个性状中的产量性状为例进行建模，其他两个同理。
> yield = Pheno_train[,1]
> yield_answer <- mixed.solve(yield, Z=m_train, K=NULL, SE = FALSE, return.Hinv=FALSE)
> str(yield_answer)
List of 5
$ Vu : num 148
$ Ve : num 324718
$ beta: num [1(1d)] 4357
$ u : num [1:1176(1d)] 0.1739 0.0373 0.319 1.9616 0.6987 ...
..- attr(*, "dimnames")=List of 1
.. ..$ : chr [1:1176] "V1" "V2" "V3" "V4" ...
$ LL : num -448
预测
基于模型对验证集的表型值进行预测。
> pred_yield_valid = m_valid %*% as.matrix(yield_answer$u)
> pred_yield = pred_yield_valid[,1] + yield_answer$beta
> pred_yield
[1] 4447.569 4468.462 4359.035 4377.164 4396.756 4408.629 4410.807 4404.205 4473.221
[10] 4400.442 4429.143 4474.376 4477.764 4470.311 4420.335 4408.489 4387.122 4383.179
[19] 4434.570 4459.401 4424.827 4521.290 4432.745 4399.275 4676.461 4643.213 4510.135
[28] 4567.752 4377.642 4408.139 4334.978 4371.179 4350.259 4335.628 4347.399 4420.289
[37] 4106.352 4156.124
模型验证
模型建立好之后，需要评价模型的准确性。最简单直接地，我们用预测的表型值和真实的表型值之间的相关性(如pearson相关性系数)来评估。也可以用均方误差(Mean-square Error, MSE)来度量估计量与被估计量之间的差异程度。
> yield_valid = Pheno_valid[,1]
> YLD_accuracy <- cor(pred_yield_valid, yield_valid, use="complete" )
> YLD_accuracy
[,1]
[1,] 0.5289435
从上结果看出，真实值和预测值的相关性约为0.53。
交叉验证
单次的建模预测具有随机性，并不能真实地反映模型的准确性，因此我们需要进行交叉验证(Cross-Validation, CV)。交叉验证也称为循环估计，是统计上将数据样本切割成较小集合的使用方法。它能从有限的学习数据中获取尽可能多的有效信息，从多个方向开始学习样本，可以有效地避免陷入局部最小值，同时也可以在一定程度上避免过拟合问题。
交叉验证一般有简单交叉验证、K-fold交叉验证、留一验证等形式，每种形式各有特点，其中K-fold交叉验证应用最广。在机器学习模型中经常会见到，后面有机会我们专门安排一期来介绍。
这里，我们使用简单交叉验证方法，循环500次来提高模型的泛化能力，将上面的代码整合一下。
## cross validation for many cycles for yield only
> traits=1
> cycles=500
> accuracy = matrix(nrow=cycles, ncol=traits)
> for(r in 1:cycles){
+ train= as.matrix(sample(1:96, 58))
+ test<-setdiff(1:96,train)
+ Pheno_train=Pheno[train,]
+ m_train=Markers_impute2[train,]
+ Pheno_valid=Pheno[test,]
+ m_valid=Markers_impute2[test,]
+
+ yield=Pheno_train[,1]
+ yield_answer<-mixed.solve(yield, Z=m_train, K=NULL, SE = FALSE, return.Hinv=FALSE)
+ pred_yield_valid = m_valid %*% as.matrix(yield_answer$u)
+ pred_yield=pred_yield_valid[,1]+yield_answer$beta
+ pred_yield
+ yield_valid = Pheno_valid[,1]
+ accuracy[r,1] <-cor(pred_yield_valid, yield_valid, use="complete" )
+ }
> mean(accuracy)
[1] 0.2797321
循环500次，模型的准确性变为0.28，可以看出单次的随机性影响还是很大的。大家感兴趣的可以用K折交叉验证的方法来试一下，看看结果有什么不同。当然，影响GS模型准确性的因素非常之多，比得到结果更重要的是如何解释结果。
本次的R包rrBLUP的实操就介绍到这里了，你学会了吗?我们下次再见。
下一篇：
相关推荐
百奥繁育 | 良种选育的信息管理系统
09-06
2023
8月25日，江西省水稻种质资源创新交流及智能育种技术研讨会在江西省农科院顺利召开。江西农业大学、江西师范大学、广东省农科院、江西省农科院、中国水稻所江西早稻中心、南昌市农科院、赣州市农科所有关专家代表以及省内种企代表共约40人参与了本次交流会议。
08-28
2023
育种4.0时代，比拼的是科技创新，关键是通过数字化、信息化、智能化让育种过程缩短周期、提升效率、降低风险，而这也是长沙百奥云数据科技有限公司（百奥云）的创业初衷。
08-22
2023
有人说，植物育种的过程就像工厂流水线。挑选的种质即原材料，根据市场需求进行设计与加工，送入田间生产线后，一代代优中选优，还得经过严格的产量、品质、抗性测试，过五关斩六将，最终拿到审定编号，成为可以推广的成熟产品。 虽然流程相似，但育种工作远不及工厂生产那样标准。田间环境气候复杂、真实性状判断困难、水肥条件难以统一……种种困难下，选育良种成了概率事件，育种家们都有种“尽人事听天命”的无奈感。
08-22
2023
1. 利用PopGenomics快速对基因组实现高质量组装 基因组组装是将全基因组测序的小片段（reads长度100 bp-100 kb）通过算法拼接成尽量长的片段（contig 和scaffold，长度几十kb 到Mb 不等）或者整条染色体的过程。获得包含基因组全序列的参考基因组是对动植物进行基因组学研究和育种利用的前提[1]。 由于植物基因组具有非常丰富的多样性，参考已发表的少数物种组装新的物种，有时却无法达到理想的组装效果。测序技术发展提供了短序列测序、单分子测序、光学图谱、Hi-C 图谱等多种测序技术及其组合的组装方案，到目前为止，已经有上千个植物基因组被组装（图1）。然而，如何以最低成本快速获得满足需求的高质量基因组，仍然是科研人员和育种数据分析人员普遍面临的一个问题。
08-22
2023
7月26日，顶着炎炎夏日 在公司CSO桂進老师的带领下 举行了一场别开生面的劳动体验。
08-04
2023