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的个数。