R包-BGLR

BGLR 是一个使用贝叶斯方法进行基因组选择的一个 R 包,可用方法包括 BayesA , BayesB , BayesC(Pi), Bayes Lasso (对应参数为 BL)。

支持的表型

  • 连续性状(censored or not)
  • 分类性状:二分类(标记为01或者12),或者有序多分类(标记为 1到 n )

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

支持的模型

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

函数参数

使用 BGLR 函数具有以下参数,其中除了表型向量 其余参数具有默认值。

1
2
3
4
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 是使用的模型(FIXEDBRRBayesA, BLBayesB, BayesC

    1
    2
    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 中缺失元素的索引
  • varESD.varE : 残差及其标准误
  • fit :模型拟合情况

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

应用举例

BGLR 自带了两个数据集(小鼠和小麦)

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

连续表型

使用小鼠数据,其模型为

其中:

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

为相应的设计矩阵。

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#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,缺失用该列均值填充)

1
2
3
4
5
6
7
> 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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#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的概率(不清楚怎么计算的),如下

1
2
3
4
5
6
7
8
> 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 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
> 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$yHatfmBin$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.

  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!
  • Copyrights © 2019-2024 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信