R包-sommer

数据分析 / R
字数统计: 3k   |   阅读时长≈ 12 分钟

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):id 或 random= ~ 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-2024 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信