View topic - GBLUP for multi-trait multi-environment data | Forum | VSN International
ZachHu Posted: Wed May 17, 2017 12:05 am
我想要构建一个来运行多性状多环境数据的gblup模型。
数据如下：
resp tr gid env
-0.157 y1 1 E1
-0.162 y1 2 E1
0.558 y2 1 E1
-0.254 y2 2 E1
0.799 y1 1 E2
-0.775 y1 2 E2
0.451 y2 1 E2
-0.083 y2 2 E2
....
列resp是在不同环境中生长的每个系的表型观测值。我还用SNP数据推断了系间的现实遗传相关性矩阵（K）。
我希望在分析模型中包含以下所有假设，但不确定如何实现它。
- 每个性状都有自己的遗传方差，且在两个环境中均不同。
- 每种性状的残差方差在两个环境中相同。
- 对这两个性状间的遗传相关很感兴趣。
但是，如果我修改假设 (1) 令Vg (y1 | E1 ) == Vg (y1 | E2)、Vg (y2 | E1 ) == Vg (y2 | E2)，我可以使用下面的代码：
fit <- asreml(fixed = resp ~ trait,
as.multivariate = tr,
random = ~ ped(gid):us(tr),
rcov=~ units:idh(tr),
data=data,
ginverse=list(gid=invG)
)
得到
gamma component std.error z.ratio constraint
ped(gid):tr!tr.y1:y1 0.8361243 0.8361243 0.18913556 4.420767 Positive
ped(gid):tr!tr.y2:y1 -0.4060155 -0.4060155 0.13895987 -2.921818 Positive
ped(gid):tr!tr.y2:y2 0.9651585 0.9651585 0.18877820 5.112659 Positive
R!variance 1.0000000 1.0000000 NA NA Fixed
R!tr.y1 0.3576103 0.3576103 0.03688209 9.696044 Positive
R!tr.y2 0.2876700 0.2876700 0.02930997 9.814749 Positive
我认为随机效应应该拟合成
random = ~ ped(gid):us(tr):idv(env)
但报错：
Warning message:
Abnormal termination
Singularity in Average Information Matrix
Results may be erroneous
到目前为止，我能想到的是将随机效应定义为
random = ~ ped(gid):us(tr):id(env)
这使我能够预测每个环境中每个系的育种值。估计的方差分量是
gamma component std.error z.ratio constraint
ped(gid):tr:env!tr.y1:y1 0.9753176 0.9753176 0.18580430 5.249166 Positive
ped(gid):tr:env!tr.y2:y1 -0.3653416 -0.3653416 0.11255951 -3.245764 Positive
ped(gid):tr:env!tr.y2:y2 0.8362916 0.8362916 0.13807872 6.056630 Positive
R!variance 1.0000000 1.0000000 NA NA Fixed
R!tr.y1 0.3550098 0.3550098 0.04438868 7.997754 Positive
R!tr.y2 0.2385684 0.2385684 0.03416360 6.983115 Positive
似乎，ped(gid):tr:env!tr.y1:y1项是y1env x gid的互作。那么，如何解释ped(gid):tr:env!tr.y2:y1？是环境与遗传相关的相互作用吗？
ZachHu Posted: Wed May 17, 2017 11:00 pm
我认为随机效应的以下设置可能满足所有这三个假设：
random = ~ ped(gid):us(tr):at(env)
得到：
gamma component std.error z.ratio constraint
ped(gid):tr:at(env, E1)!tr.y1:y1 0.7794293 0.7794293 0.19148003 4.070552 Positive
ped(gid):tr:at(env, E1)!tr.y2:y1 -0.5044509 -0.5044509 0.14821215 -3.403573 Positive
ped(gid):tr:at(env, E1)!tr.y2:y2 0.8792044 0.8792044 0.19605386 4.484504 Positive
ped(gid):tr:at(env, E2)!tr.y1:y1 1.4917692 1.4917692 0.36251091 4.115102 Positive
ped(gid):tr:at(env, E2)!tr.y2:y1 -0.2372390 -0.2372390 0.18256586 -1.299471 Positive
ped(gid):tr:at(env, E2)!tr.y2:y2 0.8042773 0.8042773 0.18445587 4.360269 Positive
R!variance 1.0000000 1.0000000 NA NA Fixed
R!tr.y1 0.3431727 0.3431727 0.04314069 7.954733 Positive
R!tr.y2 0.2218639 0.2218639 0.03253456 6.819332 Positive
上面，行1-3是环境E1中的方差分量，行4-6是E2的方差分量。在生物学上它对我有意义。
Arthur Posted: Wed May 17, 2017 11:12 pm
根据我对你模型的理解，我会按照env然后gid对数据记录性状进行排序，以便将成对的记录作为相关观察值（来自同一实验单元（env和gid）的两个性状值相关））
然后（在ASReml中）我会写：
resp ~ mu tr*env !r !{ tr.env.gid us(4).grm1(gid) !}
residual id(??).us(2)
其中??是成对的数量（实验单位）。
对于asreml-r，你需要将其转换为R语法，需要用函数str()。
你说“我认为随机效应应该这样拟合random = ~ ped(gid):us(tr):idv(env)”，而你应该写成random = ~ ped(gid):us(tr):id(env)。因为在矩阵中不能有两个缩放因子（1个隐含在us(tr)，另一个明确在idv(env)）。但是这假设env间没有协方差且每个env有相同的遗传方差。你可以写成比如random = ~ ped(gid):us(tr):coru(env)来加入相关，但这不能给每个env不同的遗传方差。我上面写的模型可以扩展为每个env给出单独的误差矩阵，通过排序env然后gid然后tr。残差部分写成residual id(#1).us(tr) id(#2).us(tr)，#替换成每个env的个数。
