R包-sommer

sommer 包可以进行 GBLUP 或 RRBLUP 分析,不过只能用于小规模数据。

sommer 包功能和注意事项

  • sommer 包是通过对 V 矩阵(mmer函数) 或 系数矩阵(mmec函数) 直接求逆从而求解 MME 的 R包,因为需要直接求逆,直接求逆的矩阵维度最好小于10000。

    通过实际测试,当矩阵维度超过 1万时,sommer 计算速度非常慢,几乎卡死。

  • 表型不能有缺失,sommer 的表型缺失填充使用中位数,有点怪异,默认会剔除含缺失的表型行。

  • 支持计算显性效应、上位效应等。

  • 支持一步法 (通过 H.mat() 构建H矩阵,但是这里 A 阵需要自己构建,文档中是建议使用 pedigreemm 的 getA 函数)

  • 支持多倍体生物的GS

安装和帮助文档

首先安装 sommer 包

1
install.packages("sommer")

安装完之后,运行下面命令打开文档

1
vignette("v1.sommer.quick.start")

之后还有 v2 到 v6 的文档

协方差结构基础

方差组分的结果有下面几种:

  • 无结构的协方差矩阵(unstructured covariance matrices, US),最常用的方差组分结构。

  • 对角化的协方差矩阵(Diagonal covariance matrices, DIAG), 非对角线元素为 0,如下

  • 复合对称协方差矩阵(Compound simmetry covariance structures, CS),如下

  • 一阶自我回归协方差矩阵(First order autoregressive covariance structures, AR1),如下

mmer 函数

sommer 包中有两个核心的函数 mmec 和 mmer 函数,它们都是利用 AI 算法和 MME 来计算方差组分和估计育种值的。二者区别在于,当需要估计的系数 (coefficients to estimate, c) 小于表型数目(records, r)时,我们调用 mmec 函数;反过来,当 c > r 时,我们使用 mmer 函数,此时我们使用直接求逆算法更快(对系数矩阵直接求逆)

由于一般来说,系数矩阵维度都是超过表型数目的,因此我们一般只会调用 mmer 函数。

mmer(fixed, random, rcov, data, weights, W, nIters=20, tolParConvLL = 1e-03,
tolParInv = 1e-06, init=NULL, constraints=NULL,method=“NR”, getPEV=TRUE,
naMethodX=“exclude”, naMethodY=“exclude”,returnParam=FALSE,
dateWarning=TRUE,date.warning=TRUE,verbose=TRUE, reshapeOutput=TRUE, stepWeight=NULL,
emWeight=NULL)

这里我们重点介绍一些常用参数:

  • fixed : 指定性状和固定效应,固定格式为 fixed = 反应变量 ~ 固定因子1 + 固定因子2 ,其中 fixed=可以省略。

  • random : 指定随机因子,固定格式为 random = ~ 随机因子1 + 随机因子2

  • rcov: 指定残差,例如 rcov = ~ units

  • data : 指定数据集

  • nIters : 最大迭代次数,默认为 20

  • tolParConvLL : 对数似然值变化的收敛阈值

  • init :方差组分初始值,默认由软件自己提供。

  • mothod :估计方差组分算法,NR 或 AI

  • naMethodX: 自变量是否可以含有缺失值,include (将缺失值填充为中位数?)或 exclude(剔除含有缺失值的行)

  • naMethodY: 自变量是否可以含有缺失值,include (将缺失值填充为中位数?)或include2 (只对多性状模型中至少有一个表型不为缺失的行进行填充,应该也是填充为中位数),或 exclude(剔除含有缺失值的行)

  • verbose :TRUE/FALSE, 是否反馈迭代过程

输出结果为一个列表,包括以下内容:

  • Vi :表型的协方差矩阵的逆矩阵
  • P :投影矩阵

多性状模型

固定因子

多性状模型固定因子设置如下,就是性状部分改为 cbind(y1,y2) , 后面和单性状一样。

1
fixed = cbind(y1,y2) ~ x

随机因子

设置随机因子如下,加入我们需要设置一个无结构化的方差组分的随机因子 re ,则设置如下(unsm(2) 中的2 为性状数目)。

对于残差项,其永远对应着 units

1
2
random = ~ vsr(re, Gtc=unsm(2))
rcov= ~ vsr(units, Gtc=unsm(2))

这里 vsr 就是指定随机效应的方差组分结构的函数,完整参数如下:

1
vsr(..., Gu=NULL, Gti=NULL, Gtc=NULL, reorderGu=TRUE, buildGu=TRUE)

其中:

  • Gu :随机因子的不同水平之间的已知的协方差矩阵(如亲缘关系矩阵,矩阵格式,行名称和列名称改为个体号),默认为单位矩阵,矩阵维度可以大于随机因子水平数目。

  • Gti :方差组分初始值,是一个 t × t 矩阵(t 为性状数目)。

    Gti 初始值需要经过 scaled,要么使用之前跑的 mmer 函数的 sigma_scaled 结果, 例如

    1
    initialVal <- mix1$sigma_scaled$`u:id`

    要么可以通过下面的公式计算。

    1
    m = x/var(y)

    where x is the desired initial value and y is the response variable.

  • Gtc :方差组分限制,是一个 t × t 矩阵(t 为性状数目),元素含义如下

    1
    2
    3
    4
    0: not to be estimated
    1: estimated and constrained to be positive (i.e. variance component)
    2: estimated and unconstrained (can be negative or positive, i.e. covariance component)
    3: not to be estimated but fixed (value has to be provided in the Gti argument)

    默认就是创建一个无结构的矩阵,即使用 unsm(2) 函数(2为性状数目,创建一个对角线元素为1,非对角线元素为2的矩阵)。也可以使用 diag() 函数指定为对角矩阵(即将非对角线元素设为0,对角线元素为 1);也可以使用 fixm() 函数指定使用方差组分初始值作为真实值(创建一个元素均为3的矩阵)。

  • reorderGu : TRUE/FALSE ,没看懂

此时的 re 具有以下的方差结构(即sommer是按照因子排序的,不是我们一般的按照个体排序)。

如果存在额外的未知的协方差结构,例如

举个例子,随机因子 id 是嵌套在因子 env 中的,那么按照以前的写法,模型应该写为 random= ~ diag(env):idrandom= ~ usr(env):id 格式,但是现在我们应该写为

1
2
fixed= y1~x
random= ~ vsr(dsr(env),id) or random= ~ vsr(usr(env),id)

这里的 dsr()usr()函数分别定义对角化的协方差结构无结构的协方差结构。在上面的例子中,假如 env 的因子数目为3,则 dsr(env) 则创建了一个 3×3 的对角矩阵。

似然比检验

似然比检验有两个用途:检查某个随机因子的显著性,检查某个特定的方差组分元素的显著性。

第一种方式,假如我们想要检查加入 the effect of a spatial kerne 是否有用,首先我们看不加入的模型

1
2
3
4
5
6
7
8
9
10
11
12
data(DT_cpdata)
DT <- DT_cpdata
GT <- GT_cpdata
MP <- MP_cpdata
### mimic two fields
A <- A.mat(GT)
mix1 <- mmer(Yield~1,
random=~vsr(id, Gu=A) +
vsr(Rowf) +
vsr(Colf),
rcov=~vsr(units), nIters=3,
data=DT, verbose = TRUE)

然后我们加入这个因子试试

1
2
3
4
5
6
7
mix2 <- mmer(Yield~1,
random=~vsr(id, Gu=A) +
vsr(Rowf) +
vsr(Colf) +
spl2Da(Row,Col),
rcov=~vsr(units), nIters=3,
data=DT,verbose = TRUE)

最后我们进行似然比检验

1
lrt <- anova(mix1, mix2)

结果如下

1
2
3
4
5
6
7
8
> lrt <- anova(mix1, mix2)
Likelihood ratio test for mixed models
==============================================================
Df AIC BIC loLik Chisq ChiDf PrChisq
mod2 8 304.5293 308.4210 -151.2647
mod1 7 305.1825 309.0741 -151.5912 0.65315 1 0.41899
==============================================================
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

第二种方式,检验一个协方差组分的显著性,例如在多性状模型中,我们想看性状之间的遗传协方差是否显著,同样地,我们可以构建两个不同的模型(加入遗传协方差和不加入遗传协方差)

这里第二个模型的写法看不懂,这里用到的模型和上面不一样,应该用一样的模型才对。

1
2
3
4
5
6
7
8
9
10
11
12
data(DT_example)
DT <- DT_example
DT$EnvName <- paste(DT$Env,DT$Name)
modelBase <- mmer(cbind(Yield, Weight) ~ Env,
random= ~ vsr(Name, Gtc=diag(2)), # here is diag()
rcov= ~ vsr(units, Gtc=unsm(2)), nIters=3,
data=DT,verbose = TRUE)
modelCov <- mmer(cbind(Yield, Weight) ~ Env,
random= ~ vsr(usr(Env),Name, Gtc=unsm(2)), # here is unsm()
rcov= ~ vsr(dsr(Env),units, Gtc=unsm(2)), nIters=3,
data=DT,verbose = TRUE)
lrt <- anova(modelBase, modelCov)

GBLUP 和 RRBLUP

GBLUP 示例代码如下,这里 DT 是表型数据框,GT 是基因型数据框(-101格式,行名称需要修改为个体号),使用 sommer 自带的 A.mat() 函数构建G阵

这里表型人为制造了 1/5 的缺失值,然后计算完 GBLUP 值之后,计算 GEBV 和表型的相关,做为准确性。

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
import sommer
data(DT_wheat)
DT <- DT_wheat
GT <- GT_wheat
colnames(DT) <- paste0("X",1:ncol(DT))
DT <- as.data.frame(DT);DT$id <- as.factor(rownames(DT))
# select environment 1
rownames(GT) <- rownames(DT)
K <- A.mat(GT) # additive relationship matrix
# GBLUP pedigree-based approach
set.seed(12345)
y.trn <- DT
vv <- sample(rownames(DT),round(nrow(DT)/5))
y.trn[vv,"X1"] <- NA
## GBLUP
ans <- mmer(X1~1,
random=~vsr(id,Gu=K),
rcov=~units, nIters=3,
data=y.trn,verbose = TRUE) # kinship based
ans$U$`u:id`$X1 <- as.data.frame(ans$U$`u:id`$X1)
# 这一步是将育种值结果的行名称中的id替换为0,但是我发现原本的行名称就没有id,因此不需要运行
#rownames(ans$U$`u:id`$X1) <- gsub("id","",rownames(ans$U$`u:id`$X1))
cor(ans$U$`u:id`$X1[vv,],DT[vv,"X1"], use="complete")

## varcomp
summary(ans)
summary(ans)$varcomp
vpredict(ans, h2 ~V1/(V1+V2)) # 遗传力

rrBLUP 方法如下,就是将 random 部分做相应的修改,其中 buildGu = FALSE 指定不创建协方差矩阵(因为水平数目太多,并且亲缘关系矩阵只是一个单位矩阵,这样可以节约内存和时间;如果你需要指定特定的亲缘关系矩阵,则不能设置这个)。

1
2
3
4
5
6
7
8
## rrBLUP
ans2 <- mmer(X1~1,
random=~vsr(list(GT), buildGu = FALSE),
rcov=~units, getPEV = FALSE, nIters=3,
data=y.trn,verbose = FALSE) # kinship based
u <- GT %*% as.matrix(ans2$U$`u:GT`$X1) # BLUPs for individuals
rownames(u) <- rownames(GT)
cor(u[vv,],DT[vv,"X1"]) # same correlation

GBLUP 实操脚本

如果自己处理基因型数据,那么首先是将芯片号替换为个体号,然后生成 012 格式的基因型数据,然后调用以下脚本

1
2
3
4
5
6
7
8
9
10
library(sommer)
library(data.table)

# 处理基因型,构建 G阵
M012 = fread("renew_chipid.raw", header=T)

G = as.matrix(M012[,!c(1:6)])
G = G - 1
rownames(G) = M012$IID
G <- A.mat(G)

然后处理表型数据,将相应转为因子

1
2
3
# 使用 sommer 进行 GS 计算
dd = read.table("phe_uncode", header=F)
for(i in 1:3) dd[,i] = as.factor(dd[,i])

最后调用 sommer 包进行 GBLUP 计算,得到方差组分和育种值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 双性状模型
mod = mmer(cbind(V4,V5) ~ V2,
random = ~ vsr(V1,Gu=G,Gtc=unsm(2)) + vsr(V3,Gtc=unsm(2)),
rcov= ~ vsr(units, Gtc=unsm(2)),
data=dd)


# 保存方差组分
write.table(summary(mod)$varcomp, "summary_mmer.txt",quote=F)

# 保存育种值
mod$U$`u:V1`$V4 <- as.data.frame(mod$U$`u:V1`$V4)
mod$U$`u:V1`$V5 <- as.data.frame(mod$U$`u:V1`$V5)

write.table(mod$U$`u:V1`$V4, "mmer_gblup1.txt",quote=F,col.names=F)
write.table(mod$U$`u:V1`$V5, "mmer_gblup2.txt",quote=F,col.names=F)

显性效应和上位效应

sommer 包可以通过 D.mat()E.mat() 来构建显性效应和上位效应的亲缘关系矩阵。

这里我们需要先复制两次个体号所在列,创建 idd 和 ide 列,分别做为显性和上位效应指定的列。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

# 显性效应和上位效应
data(DT_cpdata)
DT <- DT_cpdata
GT <- GT_cpdata
MP <- MP_cpdata
DT$idd <-DT$id; DT$ide <-DT$id
### look at the data
A <- A.mat(GT) # additive relationship matrix
D <- D.mat(GT) # dominance relationship matrix
E <- E.mat(GT) # epistatic relationship matrix

ans.ADE <- mmer(color~1,
random=~vsr(id,Gu=A) + vsr(idd,Gu=D)+ vsr(ide,Gu=E),
rcov=~units, nIters=20,
data=DT,verbose = TRUE)
summary(ans.ADE)$varcomp
vpredict(ans.ADE, h2 ~ (V1) / ( V1+V2+V3+V4) ) # narrow sense
vpredict(ans.ADE, h2 ~ (V1+V2+V3) / ( V1+V2+V3+V4) ) # broad-sense

永久环境效应

类似于显性效应和上位效应,永久环境效应也需要复制 id 列行成一列新列,不然结果会有问题(sommer应该不能指定同一列作为多个因子),举例如下

1
2
3
4
5
6
7
# 单性状重复力模型
# 需要复制id新增一列
dd$id2 = dd$V1
mod = mmer(V4 ~ V2+V3,
random = ~ vsr(V1,Gu=G,Gtc=unsm(1)) + vsr(id2,Gtc=unsm(1)),
rcov= ~ vsr(units, Gtc=unsm(1)),
data=dd)
  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!
  • Copyrights © 2019-2023 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信