【GS专栏】8-全基因组选择实战之RRBLUP

知识/政策规划

2021-03-15


  1.rrBLUP包简介

  RrrBLUP能方便快速地实现RRBLUP(Ridge Regression Best Linear Unbiased Prediction)模型,该包主要有A.matmixed.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个品种共1178SNP标记(已编码)。

  > 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模型准确性的因素非常之多,比得到结果更重要的是如何解释结果。

  本次的RrrBLUP的实操就介绍到这里了,你学会了吗?我们下次再见。


相关推荐

新品发布 | 百奥云推出基因组局部组装服务

百奥云基于二代/三代测序数据以及中间结果文件(如fq/bam),自研局部组装算法与流程,并结合了大数据产品基因云湖(GenoLake)的海量数据管理、查询与分析等功能,为用户提供个性化的局部组装服务。通过百奥云的局部组装方案,用户可以快速获取基因/蛋白全长序列信息,避免传统实验所带来的时间和资源浪费,加快决策周期,提高研究和生产效率。同时,还可有效地解决大规模测序数据存储、查询和分析等问题,提高数据利用率,获取更全面的基因组信息,助力动植物研究和育种生产的顺利进行。

04-07

2024

AI大数据 | 百奥云育种数据科学家顾林林发表一种高效精准的集成学习基因组选择方法ELPGV

是否存在一种方法可以整合各种模型的结果以期望获得更加精准的预测呢?

04-03

2024

百奥云为南繁种业发展注入新动力,水稻基因体检服务首次参展

百奥云本次大会携水稻基因体检服务首次参展,诚邀各位专家、老师及同仁共同探讨探讨未来农业发展的新趋势和技术创新。

03-21

2024

喜讯 | 百奥云助力湖南农大解析水稻耐盐杂种优势之谜

近日,湖南农业大学刘次桃/段美娟团队在国际著名期刊《Journal of Integrative Plant Biology》上发表研究论文:The OsWRKY72–OsAAT30/OsGSTU26 module mediates reactive oxygen species scavenging to drive heterosis for salt tolerance in hybrid rice。百奥云大数据部门负责人彭建祥以共同作者身份为本研究提供了个性化的数据分析支持。

03-21

2024

水稻智能育种联盟 | 水稻种质资源、品种观摩暨育种创新提升工程恳谈会圆满举行

2024年3月21日上午,水稻智能育种联盟2024水稻种质资源、品种观摩暨育种创新提升工程恳谈会圆满举行,本次会议共有三十余名水稻种业同仁参会,分别来自十余家种企及科研单位,就如何充分利用数字化能力,将我国的种质资源优势进一步转化为育种创新优势的问题齐聚一起,共商水稻种质资源交流与合作共赢,探索水稻品种自主创新能力提升路径

03-21

2024