关注我们：
【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
相关推荐
百奥繁育 | 良种选育的信息管理系统
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