【GS专栏】10-全基因组选择实战之BGLR

知识/政策规划

2021-04-06


  1. BGLR简介

  R包 BGLR于2014年最先发表在《GENETICS》[1]杂志上,至今引用已经500多次。BGLR全称“Bayesian Generalized Linear Regression”,顾名思义,主要为贝叶斯方法(包括参数和半参数)设计开发。BGLR的主函数是其同名函数,用法如下:

  BGLR(y, response_type = "gaussian", a=NULL, b=NULL,ETA = NULL, nIter = 1500,

  burnIn = 500, thin = 5, saveAt = "", S0 = NULL,

  df0 =5, R2 = 0.5, weights = NULL,

  verbose = TRUE, rmExistingFiles = TRUE, groups=NULL)

 

  主要参数说明:

  •y:数值向量,表型,可含NA;

  •ETA:一个两水平的列表,用于指定线性回归函数(默认一个截距)的协变量和其他类型的随机效应,如

  ETA=list(list(X=W,model="FIXED"),list(X=Z,model="BL"), list(K=G,model="RKHS"));

  •nIter,burnIN,thin:都为整数,分别是控制采样器迭代的次数、丢弃的样本数以及用于计算后验的稀疏方法。

  前面我们了解到,不同的贝叶斯方法,最大的区别在于标记的先验假设及密度分布不同。那么,对应到BGLR函数中,它的使用方法可总结如下表。

 

 

  更具体的用法请参考官方文档:BGLR: A Statistical Package for Whole Genome Regression and Prediction[2]。

  BGLR一直在更新维护,除了构建求解常见的Bayes模型外,BGLR还能做很多其他事。比如基因与环境互作效应(GxE)、多性状、以及多组学等复杂模型也是比较方便处理的。

  在Github中,BGLR提供了不少示例,尤其针对单性状模型,感兴趣的朋友可自行去查看。

  https://github.com/gdlc/BGLR-R

  若想试着体验一下BGLR包,最简单的方法,只需运行一个demo函数即可。

  library(BGLR)

  #BayesA

  demo(BA)

  #BayesB

  demo(BB)

  #Bayesian LASSO

  demo(BL)

  #Bayesian Ridge Regression

  demo(BRR)

  #BayesCpi

  demo(BayesCpi)

  #RKHS

  demo(RKHS)

  #Binary traits

  demo(Bernoulli)

  #Ordinal traits

  demo(ordinal)

  #Censored traits

  demo(censored)

 

  BGLR输出的结果非常丰富,它由21个元素的列表组成,包括了输入的参数(默认参数也输出)、求解的参数以及拟合的结果等,同时会生成一系列dat文件。其中,yHat就是表型预测值。

  bayes <- BGLR(phenotype,ETA = list(list(X=Markers,model="BayesA")),

  nIter = 1000,burnIn = 500,thin = 5)

  > str(bayes)

  List of 21

  $ y : num [1:242] 0.705 0.389 0.117 0.375 0.414 ...

  $ a : NULL

  $ b : NULL

  $ whichNa : int(0)

  $ saveAt : chr "F:/GS_test/"

  $ nIter : num 1000

  $ burnIn : num 500

  $ thin : num 5

  $ weights : num [1:242] 1 1 1 1 1 1 1 1 1 1 ...

  $ verbose : logi TRUE

  $ response_type: chr "gaussian"

  $ df0 : num 5

  $ S0 : num 0.295

  $ yHat : num [1:242] 0.659 0.425 0.313 0.47 0.493 ...

  $ SD.yHat : num [1:242] 0.0909 0.111 0.1379 0.1357 0.145 ...

  $ mu : num 0.587

  $ SD.mu : num 0.0797

  $ varE : num 0.053

  $ SD.varE : num 0.0081

  $ fit :List of 4

  ..$ logLikAtPostMean: num 51.6

  ..$ postMeanLogLik : num 15.8

  ..$ pD : num 71.7

  ..$ DIC : num 40.2

  $ ETA :List of 1

  ..$ :List of 18

  .. ..$ model : chr "BayesA"

  .. ..$ Name : chr "ETA_1"

  .. ..$ p : int 1000

  .. ..$ colNames : chr [1:1000] "PZB00451.1" "PZB00416.2" "PZB00579.1" "PZB00257.1" ...

  .. ..$ MSx : num 681

  .. ..$ df0 : num 5

  .. ..$ R2 : num 0.5

  .. ..$ S0 : num 0.000434

  .. ..$ shape0 : num 1.1

  .. ..$ rate0 : num 231

  .. ..$ S : num 0.000137

  .. ..$ b : Named num [1:1000] 0.000377 -0.000156 -0.000634 0.00397 0.000164 ...

  .. .. ..- attr(*, "names")= chr [1:1000] "PZB00451.1" "PZB00416.2" "PZB00579.1" "PZB00257.1" ...

  .. ..$ varB : num [1:1000] 4.35e-05 5.14e-05 4.15e-05 5.15e-05 5.07e-05 ...

  .. ..$ NamefileOut: chr "F:/GS_test/ETA_1_ScaleBayesA.dat"

  .. ..$ saveEffects: logi FALSE

  .. ..$ SD.b : Named num [1:1000] 0.00625 0.00558 0.00592 0.0074 0.00677 ...

  .. .. ..- attr(*, "names")= chr [1:1000] "PZB00451.1" "PZB00416.2" "PZB00579.1" "PZB00257.1" ...

  .. ..$ SD.varB : num [1:1000] 3.79e-05 5.66e-05 3.70e-05 8.43e-05 6.50e-05 ...

  .. ..$ SD.S : num 3.77e-05

  - attr(*, "class")= chr "BGLR"

 

  2. BGLR代码演示

  示例数据来自CIMMYT的小麦数据GYSS[3]。这里以242个小麦样本的1000个标记作为演示。

  数据加载查看

  data(GYSS)

  dim(Markers)

  Markers[1:5,1:5] #{0,1,2}

  length(phenotype)

  summary(phenotype)

 

  > dim(Markers)

  [1] 242 1000

  > Markers[1:5,1:5] #{0,1,2}

  PZB00451.1 PZB00416.2 PZB00579.1 PZB00257.1 PZB00236.2

  DT10 0 0 0 0 2

  DT100 0 2 0 0 0

  DT101 0 0 2 0 0

  DT102 0 0 0 0 0

  DT103 2 0 0 0 0

  > length(phenotype)

  [1] 242

  > summary(phenotype)

  Min. 1st Qu. Median Mean 3rd Qu. Max.

  0.003333 0.334914 0.529265 0.551283 0.725641 1.396603

 

  编写函数,导出目标结果

  library(BGLR)

  bayes_fun <- function(phenotype,Markers,model){

  bayes <- BGLR(phenotype,ETA = list(list(X=Markers,model=model)),

  nIter = 1000,burnIn = 500,thin = 5,saveAt = "",

  S0 = NULL,R2 = 0.5, weights = NULL,verbose = TRUE,

  rmExistingFiles = TRUE, groups=NULL)

  bayes_pred <- bayes$yHat

  return(list(predict_value=bayes_pred,

  r=cor(bayes_pred,phenotype)))

  }

 

  运行Bayes各类模型

  ## BayesA

  bayesA <- bayes_fun(phenotype,Markers,model="BayesA")

  ## BayesB

  bayesB <- bayes_fun(phenotype,Markers,model="BayesB")

  ## BayesC

  bayesC <- bayes_fun(phenotype,Markers,model="BayesC")

  ## Bayesian Lasso

  bl <- bayes_fun(phenotype,Markers,model="BL")

  bayesL <- bl$r

  ## Bayesian Ridge Regression

  brr <- bayes_fun(phenotype,Markers,model="BRR")

 

  Bayes模型结果比较

  > bayesA$r

  [1] 0.9141349

  > bayesB$r

  [1] 0.9110653

  > bayesC$r

  [1] 0.8829751

  > bl$r

  [1] 0.9110347

  > brr$r

  [1] 0.8815097

 

  总体而言,各类贝叶斯参数方法结果差异并不是很大。

 

  半参数方法RKHS示例

 

  # Computing the genomic relationship matrix

  X <-scale(x = Markers,center=TRUE,scale=TRUE)

  G <-tcrossprod(X)/ncol(X)

  rkhs <- BGLR(y=phenotype, ETA=list(list(K=G, model= "RKHS")),

  nIter = 1000, burnIn = 500, thin = 5)

  rkhs_pred <- rkhs$yHat

  rkhs_r <- cor(rkhs_pred,phenotype)

  rkhs_r

 

  > rkhs_r

  [1] 0.9018886

 

  RKHS的结果也很接近。当然这里只是一个粗略的比较,数据量也较小。贝叶斯类方法是非常耗时的,当数据量较大时,要酌情考虑。

 

  References

  [1] 《GENETICS》: https://www.genetics.org/content/198/2/483

  [2] BGLR: A Statistical Package for Whole Genome Regression and Prediction: http://genomics.cimmyt.org/BGLR-extdoc.pdf

  [3] CIMMYT的小麦数据GYSS: https://repository.cimmyt.org/xmlui/handle/10883/2976


相关推荐

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

百奥云基于二代/三代测序数据以及中间结果文件(如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