R包-BGLR

BGLR 是一个使用贝叶斯方法进行基因组选择的一个 R 包。

支持的表型

• 连续性状（censored or not）
• 分类性状：二分类（标记为01或者12），或者有序多分类（标记为 1到 n ）

表型可以有缺失，但是自变量不能有缺失。

支持的模型

支持的回归模型如下，这里 Flat 就是固定因子，

函数参数

使用 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 。

• response_type : 表型向量的类型，“gaussian” 或 “ordinal”

• ETA : 一个双层列表(list)，用于指定自变量，默认是只有截距，具体用法见下面的例子。

  里面的每一个子列表表示一种类型的自变量，其中 X 是相应的设计矩阵（可以用 model.matrix 语法，见下面应用中的例子），model 是使用的模型（FIXED ，BRR，BayesA, BL，BayesB, BayesC）

  ETA=list(list(X=W, model="FIXED"),
  list(X=Z,model="BL")),

• nIter,burnIn, thin: (integer) the number of iterations, burn-in and thinning .

• saveAt : 字符串，输出文件前缀。

输出

BGLR 函数会返回一个列表，其中有后验均值估计值及其标准误估计值，其元素包括:

• y ：表型
• whichNa ：y 中缺失元素的索引
• varE 和 SD.varE ： 残差及其标准误
• fit ：模型拟合情况

同时也会输出一些文件，为某些数值的采样结果，用于用户评估收敛情况。举例而言，mu.dat 就是截距的采样结果。

应用举例

BGLR 自带了两个数据集（小鼠和小麦）

• 小鼠数据中包含 1814 个个体， 10346 个 SNP
• 小麦数据中包含 599 个个体

连续表型

使用小鼠数据，其模型为

其中：

• 为性别和窝大小的效应，为固定因子
• 为笼子的效应，视为随机因子（先验分布为正态分布）
• 为标记的效应，先验分布为 double-exponential

为相应的设计矩阵。

使用 BGLR 的脚本如下，注意这里的 ETA 部分使用了 model.matrix 函数的语法来生成设计矩阵。

#1# Loading and preparing the input data
library(BGLR); data(mice);
Y<-mice.pheno; X<-mice.X; A=mice.A;
y<-Y$Obesity.BMI; y<-(y-mean(y))/sd(y)

#2# Setting the linear predictor
ETA<-list( list(~factor(GENDER)+factor(Litter),
data=Y,model='FIXED'),
list(~factor(cage),data=Y, model='BRR'),
list(X=X, model='BL')
)

#3# Fitting the model
fm<-BGLR(y=y,ETA=ETA, nIter=12000, burnIn=2000)
save(fm,file='fm.rda')

这里基因型数据就是一个矩阵，一行是一个样本（与 y 保持一致），一列是一个位点（基因型为 012，缺失用该列均值填充）

> X[0:5, 0:5]
           rs3683945_G rs3707673_G rs6269442_G rs6336442_G rs13475700_A
A048005080           1           1           1           1            0
A048006063           1           1           2           1            1
A048006555           2           0           2           2            0
A048007096           1           1           1           1            1
A048010273           2           0           2           2            0

然后我们可以从输出的列表中提取结果，其中

• #1# 部分是提取所有位点的效应估计值及其标准误，这里因为标记效应是第 3 个自变量，因此使用 fm$ETA[[3]] 。
• #2# 部分是个体的基因型值（所有位点加性效应的总和）与表型的散点图。
• #3# 部分是提取模型的拟合情况（DIC）和残差
• #4# 部分是提取后验分布的采样结果。

#1# Estimated Marker Effects & posterior SDs
par(mfrow=c(2,2))
bHat<- fm$ETA[[3]]$b
SD.bHat<- fm$ETA[[3]]$SD.b
plot(bHat^2, ylab='Estimated Squared-Marker Effect',
type='o',cex=.5,col=4,main='Marker Effects')

#2# Predictions
# Total prediction
# yHat<-fm$yHat
# tmp<-range(c(y,yHat))
# plot(yHat~y,xlab='Observed',ylab='Predicted',col=2,
# xlim=tmp,ylim=tmp); abline(a=0,b=1,col=4,lwd=2)
# Just the genomic part
gHat<-X%*%fm$ETA[[3]]$b
tmp<-range(c(y,gHat))
plot(gHat~y,xlab='Phenotype',
ylab='Predicted Genomic Value',col=2,
xlim=tmp,ylim=tmp); abline(a=0,b=1,col=4,lwd=2)

#3# Godness of fit and related statistics
fm$fit
fm$varE # compare to var(y)

#4# Trace plots
list.files()
# Residual variance
varE<-scan('varE.dat')
plot(varE,type='o',col=2,cex=.5,ylab=expression(var[e]));
abline(h=fm$varE,col=4,lwd=2);
abline(v=fm$burnIn/fm$thin,col=4)
# lambda (regularization parameter of the Bayesian Lasso)
lambda<-scan('ETA_3_lambda.dat')
plot(lambda,type='o',col=2,cex=.5,ylab=expression(lambda));
abline(h=fm$ETA[[3]]$lambda,col=4,lwd=2);
abline(v=fm$burnIn/fm$thin,col=4)

画图结果见下面，这里图3以及图4和文献中的结果不一样，文献中总共有 24000 个点，但是我得到的 ‘varE.dat’ 和 ‘ETA_3_lambda.dat’ 均只有 2400 个结果，暂时不清楚原因。

分类表型

对于分类变量，我们需要设置 response_type='ordinal' ，我们使用小麦的数据进行分析如下。

#1# Loading and preparing the input data
library(BGLR); data(wheat);
Y<-wheat.Y; X<-wheat.X; A<-wheat.A;
y<-Y[,1]
tst<-sample(1:nrow(X),size=150)

#2# Binary outcome
yBin<-ifelse(y>0,1,0)
yBinNA<-yBin ; yBinNA[tst]<-NA
ETA<-list(list(X=X,model='BL'))
fmBin<-BGLR(y=yBinNA,response_type='ordinal', ETA=ETA,
nIter=1200,burnIn=200)
head(fmBin$probs)
par(mfrow=c(1,2))
boxplot(fmBin$probs[-tst,2]~yBin[-tst],main='Training',ylab='Estimated prob.')
boxplot(fmBin$probs[tst,2]~yBin[tst],main='Testing', ylab='Estimated prob.')

#3# Ordinal outcome
yOrd<-ifelse(y<quantile(y,1/4),1,ifelse(y<quantile(y,3/4),2,3))
yOrdNA<-yOrd ; yOrdNA[tst]<-NA
ETA<-list(list(X=X,model='BL'))
fmOrd<-BGLR(y=yOrdNA,response_type='ordinal', ETA=ETA,
nIter=1200,burnIn=200)
head(fmOrd$probs)

分析过程和上面一样，输出结果中 fmBin$ETA[[1]]$b 为标记效应，fmBin$prob 给出每个样本估计表型是0或者1的概率（不清楚怎么计算的），如下

> head(fmBin$prob)
             0         1
[1,] 0.4488623 0.5511377
[2,] 0.6177665 0.3822335
[3,] 0.6119654 0.3880346
[4,] 0.3914834 0.6085166
[5,] 0.3705525 0.6294475
[6,] 0.4130708 0.5869292

上面脚本中给出了参考群和验证群的真实表型和预测概率的箱线图，如下

我想看一下所有样本的加性效应值结果，发现它和 fmBin$yHat 结果是一样的，因为 fmBin$mu 等于 0 。

> gHat<-X%*%fmBin$ETA[[1]]$b
> head(gHat)
          [,1]
[1,] -1.359382
[2,] -1.828429
[3,] -1.810299
[4,] -1.204968
[5,] -1.152511
[6,] -1.272660
> head(fmBin$yHat)
      775      2166      2167      2465      3881      3889
-1.359382 -1.828429 -1.810299 -1.204968 -1.152511 -1.272660
> fmBin$mu
[1] 0

自己查看了一下 fmBin$yHat 和 fmBin$prob 的相关系数在 99.9% 以上，二者应该是一致的。而且，我发现有一个 fmBin$threshold 为 -1.5 ，计算 （这里 是标准正态分布的累积分布函数）和 fmBin$prob 很相似，但是不完全一样。

改成使用 BayesC 方法，就是将脚本中的 model='BL' 改为 model='BayesC' ，结果见下图，比上面感觉好一点。

参考文献

1. Pérez P, de Los Campos G. Genome-wide regression and prediction with the BGLR statistical package[J]. Genetics, 2014, 198(2): 483-495.

2. Pérez P, de los Campos G. BGLR: a statistical package for whole genome regression and prediction[J]. Genetics, 2014, 198(2): 483-495.