【原】環(huán)境遺傳相關(guān) | 育種中的基因與環(huán)境互作

育種數(shù)據(jù)分析 2021-11-18

展開全文

基因與環(huán)境互作

基因與環(huán)境互作，植物中同一個品種多年多點種植，評價基因與環(huán)境互作，找到品種最適合推廣的區(qū)域。

動物育種中，由于動物基因型不一樣，一般用同一個家系在不同的環(huán)境中養(yǎng)殖，雖然進行不一樣，但是個體間有親緣關(guān)系，也可以計算基因與環(huán)境互作（或者叫做環(huán)境遺傳相關(guān)）。比如種豬培育時，生長條件盡量和商品種條件類似，可以避免選擇種豬表現(xiàn)優(yōu)良而商品種表現(xiàn)不好的品系（遺傳相關(guān)較低，環(huán)境互作較高）。

基因組時代的到來，可以構(gòu)建G矩陣或者H矩陣來表示個體間的親緣關(guān)系，這樣也可以計算基因與環(huán)境互作，選擇特定環(huán)境下合適的品系，也可以實現(xiàn)了。

育種中，數(shù)據(jù)分析的不同模型，都是為育種服務(wù)的，通過模型了解應(yīng)用場景，然后會更深刻的了解育種的選擇。

1. 環(huán)境

宏環(huán)境

Macro-environment

環(huán)境 (Environment): 環(huán)境定義為影響一個基因型表現(xiàn)的一組非遺傳因素。

非生物因素, 如土壤的物理和化學特性、氣候因子（如光照，降雨量和溫度）等
生物因素, 包含害蟲、病原體、線蟲和雜草等

微環(huán)境

微環(huán)境 (Micro-environment): 微環(huán)境定義為單個植株或小區(qū)的生長環(huán)境，兩個不同的植株或小區(qū)生長在同樣的微環(huán)境中生長的可能性幾乎是0。

「基因與環(huán)境互作中環(huán)境一般是指宏環(huán)境，微環(huán)境一般視為隨機誤差效應(yīng)」

2. 基因型與環(huán)境互作的幾種模式

模式一：一個基因型在兩個環(huán)境下都優(yōu)于另一個基因型，差異相等，兩條線平行，說明不存在基因型與環(huán)境互作
模式二：一個基因型在兩個環(huán)境下都優(yōu)于另一個基因型，基因型間的差異在兩個環(huán) 境下不相等，基因型 2 隨著環(huán)境效應(yīng)的增加表現(xiàn)出更大的優(yōu)勢。在這種模式下，基型之間的差異因環(huán)境而變，也就是說存在基因型和環(huán)境間的互作
模式三：是交叉互作，基因型間差異的絕對值在兩個環(huán)境下是相等的，這時的基因型效應(yīng)為 0，只存在環(huán)境效應(yīng)和互作效應(yīng)
模式四：是交叉互作，基因型間差異的絕對值在兩個環(huán)境下不相等

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

3. 基因型與環(huán)境互作的利用方式

忽略它，選擇平均值高的，適應(yīng)性廣的品種
降低它，將環(huán)境分組，分區(qū)域推廣品種
利用它，強調(diào)品種對特殊環(huán)境的適應(yīng)性，并利用它推廣特定品種最優(yōu)的地點

4. 植物中的基因與環(huán)境互作

「特點：」同一個基因型（ID）可以有多個個體，分別種植于不同環(huán)境中，用一般線性模型就可以分析。

這里的i為基因型，j為地點，k為重復(fù)

這個很簡單，也是最常見的用法。更復(fù)雜的模型，比如各地點方差異質(zhì)，我們后面介紹。這里我們重點介紹動物中基因與環(huán)境互作的方法。

5. 動物中的基因與環(huán)境互作

「特點：」

動物同一基因型（不考慮雙胞胎或者克隆體）只能在一個環(huán)境中，要想計算基因與環(huán)境互作，需要利用半同胞、全同胞或者系譜的關(guān)系，放到動物模型的框架下進行分析。

「常用的數(shù)據(jù)格式：」

利用系譜構(gòu)成的A矩陣，用多性狀動物模型計算，比如半同胞的個體在不同的環(huán)境中，親代子代的個體在不同的環(huán)境中
利用全基因組SNP信息構(gòu)建G矩陣（或者H矩陣），利用多性狀GBLUP，計算環(huán)境的遺傳相關(guān)

「系譜AbLUP計算環(huán)境遺傳相關(guān)數(shù)據(jù)量要求：」

至少要50~100個家系在不同的環(huán)境中
每個家系至少要包含50~100個個體

「計算方法：」

將兩個環(huán)境的數(shù)據(jù)，變?yōu)閮闪袛?shù)據(jù)，利用雙性狀動物模型進行分析
性狀的遺傳相關(guān)，即為環(huán)境的遺傳相關(guān)
因為同一個個體只能在一個環(huán)境匯總，所以殘差的協(xié)方差組分為0，矩陣為diag矩陣（對角矩陣）

6. 代碼演示

原始數(shù)據(jù)，包括系譜數(shù)據(jù)和表型數(shù)據(jù)，表型數(shù)據(jù)觀測值為phe，有兩個環(huán)境場地（A和B），現(xiàn)在要計算A和B的環(huán)境遺傳相關(guān)：

> 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

「整理數(shù)據(jù)，變?yōu)槎嘈誀钅Ｐ偷臄?shù)據(jù)格式：」

> 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矩陣，構(gòu)建雙性狀動物模型，加性用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. 非常大的坑

因為我的數(shù)據(jù)是模擬的數(shù)據(jù)，兩個場的相關(guān)系數(shù)幾乎為1，用us矩陣時，發(fā)現(xiàn)相關(guān)系數(shù)為0，反復(fù)測試一直失敗，后來我用asreml中的corgh函數(shù)，試了一下，果然成功了。

corgh函數(shù)，輸出的直接是遺傳相關(guān)及其標準誤，原來我以為它和變換后的us矩陣等價，看起來我還是太年輕了。

「us錯誤的結(jié)果：」

> 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

可以看到，遺傳相關(guān)竟然為0.。。。，這是錯誤的！

「corgh正確的結(jié)果：」

> 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