大家好，我是邓飞。

基因与环境互作，在育种中也称为品种与环境互作，它是指某些品种在不同环境中表现的差异较大，存在互作效应。

在植物或者林木中，它是指同一个品种或者无性系在不同的环境中表现不同，有些环境适合特定的品种。

在动物中或者水产中，由于个体只有一个，没有克隆体，所以这里的个体通常指品种或者品系或者家系，也可以分析出某些家系适合的环境或者需要避免的环境。

基因与环境互作，弄清楚它的模式，我们可以利用（特定环境最好的品种，特定品种最好的环境），我们可以避免（避免一些环境），我们可以选择（选择特定环境下的品种，或者选择稳定性强的品种），很多手段可以用。

GWAS和GS分析中，也可以分析位点与环境的互作，人类PRS计算中也可以考虑基因与环境互作。

1. 环境

宏环境

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中的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），分解基因与环境互作。