基因与环境互作

基因与环境互作，植物中同一个品种多年多点种植，评价基因与环境互作，找到品种最适合推广的区域。

动物育种中，由于动物基因型不一样，一般用同一个家系在不同的环境中养殖，虽然进行不一样，但是个体间有亲缘关系，也可以计算基因与环境互作（或者叫做环境遗传相关）。比如种猪培育时，生长条件尽量和商品种条件类似，可以避免选择种猪表现优良而商品种表现不好的品系（遗传相关较低，环境互作较高）。

基因组时代的到来，可以构建G矩阵或者H矩阵来表示个体间的亲缘关系，这样也可以计算基因与环境互作，选择特定环境下合适的品系，也可以实现了。

育种中，数据分析的不同模型，都是为育种服务的，通过模型了解应用场景，然后会更深刻的了解育种的选择。

1. 环境

宏环境

Macro-

环境 (): 环境定义为影响一个基因型表现的一组非遗传因素。

微环境

微环境 (Micro-): 微环境定义为单个植株或小区的生长环境，两个不同的植株或小区生长在同样的微环境中生长的可能性几乎是0。

「基因与环境互作中环境一般是指宏环境，微环境一般视为随机误差效应」

2. 基因型与环境互作的几种模式

最常见的互作是模式二和模式四。

3. 基因型与环境互作的利用方式4. 植物中的基因与环境互作

「特点：」同一个基因型（ID）可以有多个个体，分别种植于不同环境中，用一般线性模型就可以分析。

y_{ijk} = \mu.. + G_i + E_j + GE_{ij} + \{ijk}

这里的i为基因型，j为地点，k为重复

这个很简单，也是最常见的用法。更复杂的模型，比如各地点方差异质，我们后面介绍。这里我们重点介绍动物中基因与环境互作的方法。

5. 动物中的基因与环境互作

「特点：」

动物同一基因型（不考虑双胞胎或者克隆体）只能在一个环境中，要想计算基因与环境互作，需要利用半同胞、全同胞或者系谱的关系，放到动物模型的框架下进行分析。

「常用的数据格式：」

「系谱AbLUP计算环境遗传相关数据量要求：」

「计算方法：」

6. 代码演示

原始数据，包括系谱数据和表型数据，表型数据观测值为phe，有两个环境场地（A和B），现在要计算A和B的环境遗传相关：

> ped = asreml.read.table("ped.csv",header=T,sep=",")
> head(ped)
  ID Sire Dam
1  1    0   0
2  2    0   0
3  3    0   0
4  4    0   0
5  5    0   0
6  6    0   0
> dat = asreml.read.table("phe.csv",sep=",",header=T)
> head(dat)
  ID Sire Dam  Fam Changdi Sex      phe
1 61    1  54 1_54       A   M 196.8497
2 62    1  54 1_54       A   M 178.1221
3 63    1  54 1_54       A   M 163.6030
4 64    1  54 1_54       A   M 226.2328
5 65    1  54 1_54       A   F 215.7228
6 66    1  54 1_54       A   M 153.7567

复制

「整理数据，变为多性状模型的数据格式：」

> dat$phe_A = dat$phe
> dat$phe_B = dat$phe
> dat[dat$Changdi == "A",]$phe_B = NA
> dat[dat$Changdi == "B",]$phe_A = NA
> head(dat)
  ID Sire Dam  Fam Changdi Sex      phe    phe_A phe_B
1 61    1  54 1_54       A   M 196.8497 196.8497    NA
2 62    1  54 1_54       A   M 178.1221 178.1221    NA
3 63    1  54 1_54       A   M 163.6030 163.6030    NA
4 64    1  54 1_54       A   M 226.2328 226.2328    NA
5 65    1  54 1_54       A   F 215.7228 215.7228    NA
6 66    1  54 1_54       A   M 153.7567 153.7567    NA
> tail(dat)
         ID Sire  Dam       Fam Changdi Sex      phe phe_A    phe_B
14995 15055 5609 5207 5609_5207       B   M 253.1178    NA 253.1178
14996 15056 5609 5207 5609_5207       B   M 229.0524    NA 229.0524
14997 15057 5609 5207 5609_5207       B   F 247.3232    NA 247.3232
14998 15058 5609 5207 5609_5207       B   M 285.1402    NA 285.1402
14999 15059 5609 5207 5609_5207       B   M 243.7538    NA 243.7538
15000 15060 5609 5207 5609_5207       B   F 243.6527    NA 243.6527

复制

利用双性状的us矩阵，构建双性状动物模型，加性用us矩阵，残差用diag矩阵：

mod3 = asreml(cbind(phe_A,phe_B) ~ trait + Sex, random = ~ us(trait):vm(ID,ainv), 
              na.action = na.method(y = "include",x = "include"),
              residual = ~ units:diag(trait),data=dat)
summary(mod3)$varcomp
vpredict(mod3,rg ~ V2/sqrt(V1*V3))

复制

7. 非常大的坑

因为我的数据是模拟的数据，两个场的相关系数几乎为1，用us矩阵时，发现相关系数为0，反复测试一直失败，后来我用中的corgh函数，试了一下，果然成功了。

corgh函数，输出的直接是遗传相关及其标准误，原来我以为它和变换后的us矩阵等价，看起来我还是太年轻了。

「us错误的结果：」

> summary(mod3)$varcomp
                                     component std.error   z.ratio bound %ch
trait:vm(ID, ainv)!trait_phe_A:phe_A 327.89751  50.53725  6.488234     P 0.1
trait:vm(ID, ainv)!trait_phe_B:phe_A   0.00001        NA        NA     F 0.0
trait:vm(ID, ainv)!trait_phe_B:phe_B 330.40559  49.33407  6.697311     P 0.0
units:trait!R                          1.00000        NA        NA     F 0.0
units:trait!trait_phe_A              524.65267  27.71968 18.927086     P 0.0
units:trait!trait_phe_B              607.18007  27.66161 21.950282     P 0.0
> vpredict(mod3,rg ~ V2/sqrt(V1*V3))
       Estimate           SE
rg 3.038136e-08 3.259792e-09

复制

可以看到，遗传相关竟然为0.。。。，这是错误的！

「corgh正确的结果：」

> mod4 = asreml(cbind(phe_A,phe_B) ~ trait + Sex, random = ~ corgh(trait):vm(ID,ainv), 
+               na.action = na.method(y = "include",x = "include"),
+               residual = ~ units:diag(trait),data=dat)
Model fitted using the sigma parameterization.
ASReml 4.1.0 Mon May 10 21:15:17 2021
          LogLik        Sigma2     DF     wall    cpu
 1     -57872.83           1.0  14997 21:15:18    0.5 (1 restrained)
 2     -57525.98           1.0  14997 21:15:18    0.3 (1 restrained)
 3     -57236.46           1.0  14997 21:15:18    0.4 (1 restrained)
 4     -57122.92           1.0  14997 21:15:19    0.4 (1 restrained)
 5     -57099.16           1.0  14997 21:15:19    0.4 (1 restrained)
 6     -57096.26           1.0  14997 21:15:19    0.4
 7     -57096.15           1.0  14997 21:15:20    0.4
 8     -57096.15           1.0  14997 21:15:20    0.4
> summary(mod4)$varcomp
                                                 component    std.error    z.ratio bound %ch
trait:vm(ID, ainv)!trait!phe_B:!trait!phe_A.cor   0.990679  0.008861388 111.797266     U   0
trait:vm(ID, ainv)!trait_phe_A                  325.464194 46.679516413   6.972313     P   0
trait:vm(ID, ainv)!trait_phe_B                  339.460777 48.562234893   6.990221     P   0
units:trait!R                                     1.000000           NA         NA     F   0
units:trait!trait_phe_A                         525.666184 25.881161866  20.310765     P   0
units:trait!trait_phe_B                         602.430876 27.227647834  22.125704     P   0

复制

可以看到遗传相关为0.99，标准误为0.008。

A场和B场的遗传相关为0.99，说明基本没有基因与环境互作。

9. 更高级的模型

两个地点的遗传相关可以用双性状模型，三个地点的遗传相关可以用三性状模型，如果是22个地点呢？？？

22个性状的多性状模型，运算量大而且不容易收敛。可以考虑 利用随机回归模型，勒让德多项式（ ），分解基因与环境互作。

「示例代码：」

10. 飞哥感言

基因与环境互作，分析环境遗传相关的原理，代码及结果解析就是这样，有什么问题欢迎沟通。