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 | random = ~ vsr(re, 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 andy
is the response variable. -
Gtc
:方差组分限制,是一个 t × t 矩阵(t 为性状数目),元素含义如下1
2
3
40: 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 | fixed= y1~x |
这里的 dsr()
和usr()
函数分别定义对角化的协方差结构和无结构的协方差结构。在上面的例子中,假如 env 的因子数目为3,则 dsr(env)
则创建了一个 3×3 的对角矩阵。
似然比检验
似然比检验有两个用途:检查某个随机因子的显著性,检查某个特定的方差组分元素的显著性。
第一种方式,假如我们想要检查加入 the effect of a spatial kerne
是否有用,首先我们看不加入的模型
1 | data(DT_cpdata) |
然后我们加入这个因子试试
1 | mix2 <- mmer(Yield~1, |
最后我们进行似然比检验
1 | lrt <- anova(mix1, mix2) |
结果如下
1 | > lrt <- anova(mix1, mix2) |
第二种方式,检验一个协方差组分的显著性,例如在多性状模型中,我们想看性状之间的遗传协方差是否显著,同样地,我们可以构建两个不同的模型(加入遗传协方差和不加入遗传协方差)
这里第二个模型的写法看不懂,这里用到的模型和上面不一样,应该用一样的模型才对。
1 | data(DT_example) |
GBLUP 和 RRBLUP
GBLUP 示例代码如下,这里 DT 是表型数据框,GT 是基因型数据框(-101格式,行名称需要修改为个体号),使用 sommer 自带的 A.mat() 函数构建G阵
这里表型人为制造了 1/5 的缺失值,然后计算完 GBLUP 值之后,计算 GEBV 和表型的相关,做为准确性。
1 | import sommer |
rrBLUP 方法如下,就是将 random 部分做相应的修改,其中 buildGu = FALSE
指定不创建协方差矩阵(因为水平数目太多,并且亲缘关系矩阵只是一个单位矩阵,这样可以节约内存和时间;如果你需要指定特定的亲缘关系矩阵,则不能设置这个)。
1 | ## rrBLUP |
GBLUP 实操脚本
如果自己处理基因型数据,那么首先是将芯片号替换为个体号,然后生成 012 格式的基因型数据,然后调用以下脚本
1 | library(sommer) |
然后处理表型数据,将相应转为因子
1 | # 使用 sommer 进行 GS 计算 |
最后调用 sommer 包进行 GBLUP 计算,得到方差组分和育种值
1 | # 双性状模型 |
显性效应和上位效应
sommer 包可以通过 D.mat()
和 E.mat()
来构建显性效应和上位效应的亲缘关系矩阵。
这里我们需要先复制两次个体号所在列,创建 idd 和 ide 列,分别做为显性和上位效应指定的列。
1 |
|
永久环境效应
类似于显性效应和上位效应,永久环境效应也需要复制 id 列行成一列新列,不然结果会有问题(sommer应该不能指定同一列作为多个因子),举例如下
1 | # 单性状重复力模型 |