一千零一技|如何计算环境遗传相关：育种中的基因与环境互作

基因与环境互作

1. 环境

「参考：」

❝

http://www./common/UploadFiles/file/teaching/%E6%95%B0%E9%87%8F%E9%81%97%E4%BC%A0%E6%95%99%E5%AD%A62013/Ch7.pdf

❞

宏环境

Macro-environment

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

• 非生物因素, 如土壤的物理和化学特性、气候因子（如光照，降雨量和温度）等
• 生物因素, 包含害虫、病原体、线虫和杂草等

微环境

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

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

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

• 模式一：一个基因型在两个环境下都优于另一个基因型，差异相等，两条线平行，说明不存在基因型与环境互作
• 模式二：一个基因型在两个环境下都优于另一个基因型，基因型间的差异在两个环 境下不相等，基因型 2 随着环境效应的增加表现出更大的优势。在这种模式下，基型之间的差异因环境而变，也就是说存在基因型和环境间的互作
• 模式三：是交叉互作，基因型间差异的绝对值在两个环境下是相等的，这时的基因型效应为 0，只存在环境效应和互作效应
• 模式四：是交叉互作，基因型间差异的绝对值在两个环境下不相等

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

3. 基因型与环境互作的利用方式

• 忽略它，选择平均值高的，适应性广的品种
• 降低它，将环境分组，分区域推广品种
• 利用它，强调品种对特殊环境的适应性，并利用它推广特定品种最优的地点

4. 植物中的基因与环境互作

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

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

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

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

「特点：」

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

「常用的数据格式：」

• 利用系谱构成的A矩阵，用多性状动物模型计算，比如半同胞的个体在不同的环境中，亲代子代的个体在不同的环境中
• 利用全基因组SNP信息构建G矩阵（或者H矩阵），利用多性状GBLUP，计算环境的遗传相关

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

• 至少要50~100个家系在不同的环境中
• 每个家系至少要包含50~100个个体

「计算方法：」

• 将两个环境的数据，变为两列数据，利用双性状动物模型进行分析
• 性状的遗传相关，即为环境的遗传相关
• 因为同一个个体只能在一个环境汇总，所以残差的协方差组分为0，矩阵为diag矩阵（对角矩阵）

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，反复测试一直失败，后来问了大神（asreml的大神「栾生老师」），他建议我用asreml中的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个性状的多性状模型，运算量大而且不容易收敛。可以考虑 利用随机回归模型，勒让德多项式（Legendre polynomials），分解基因与环境互作。

「示例代码：」

10. 飞哥感言

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