用R进行G×E互作的AMMI模型分析

本博文最早发在netlify个人博客

在多地点试验中,突出的一个问题是林木基因型与环境之间往往存在显著的交互作用(Genotype by environment interaction, GEI),因此如何准确评估GEI对于后续林木良种的选育和推广至关重要。目前,农业上大多使用联合回归法 1、主效可加互作可乘模型(Additive main effects multiplicative interaction, AMMI) 2和基因型主效加基因型-环境互作效应双标图(Genotype Main Effect plus Genotype-by-Environment Interaction biplot, GGE双标图) 3分析多点试验,并据此来进行品种评价、试验点评价和品种生态区划分,其中GGE双标图越来越受关注。迄今,GGE双标图在林木上的应用 4仍然很少,仅在少数树种如辐射松、杨树和乐昌含笑中有过报道。

00 示例数据集

示例采用agricolae包的数据集plrv。

    data(plrv,package='agricolae')

    names(plrv)
    str(plrv)

01 使用AMMI()建模

使用agricolae包的AMMI()建模,具体代码如下:

    model<- with(plrv,AMMI(ENV=Locality,GEN=Genotype,REP=Rep,
                               Y=Yield,console=F,PC=T))

现在可查看方差分析结果与GEI部分的主成分结果如下:

    model$ANOVA # for aov
## Analysis of Variance Table
## 
## Response: Y
##            Df Sum Sq Mean Sq  F value    Pr(>F)    
## ENV         5 122284 24456.9 257.0382  9.08e-12 ***
## REP(ENV)   12   1142    95.1   2.5694  0.002889 ** 
## GEN        27  17533   649.4  17.5359 < 2.2e-16 ***
## ENV:GEN   135  23762   176.0   4.7531 < 2.2e-16 ***
## Residuals 324  11998    37.0                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    
    model$analysis # for pc
##     percent  acum Df     Sum.Sq   Mean.Sq F.value   Pr.F
## PC1    56.3  56.3 31 13368.5954 431.24501   11.65 0.0000
## PC2    27.1  83.3 29  6427.5799 221.64069    5.99 0.0000
## PC3     9.4  92.7 27  2241.9398  83.03481    2.24 0.0005
## PC4     4.3  97.1 25  1027.5785  41.10314    1.11 0.3286
## PC5     2.9 100.0 23   696.1012  30.26527    0.82 0.7059

如果只保留前2个主成分,则可通过自编函数来操作:

   ## only keep the first two PCs
    PC.res <- function(model, pcN=2) {
      maov<-model$ANOVA # for aov
      df.res<-maov['Residuals','Df'] # residual df
      MS.res<-maov['Residuals','Mean Sq']
      
      cat('GEI aov results:\n')
      print(maov['ENV:GEN',])
      cat('\n\n')
      
      GEI<-model$analysis # for pc
      row.names(GEI)
      nr<-nrow(GEI)
      for(i in 1:7) GEI[nr+1,i]<-sum(GEI[(1+pcN):nr,i])
      row.names(GEI)[nr+1]<-'PCres'
      
      GEI['PCres','acum']=GEI[nr,'acum']
      GEI['PCres','Mean.Sq']=round(GEI['PCres','Sum.Sq']/GEI['PCres','Df'],4)
      GEI['PCres','F.value']=round(GEI['PCres','Mean.Sq']/MS.res,4)
      GEI['PCres','Pr.F']=round(pf(GEI['PCres','F.value'], GEI['PCres','Df'], df.res,lower.tail=F),4)
      GEI$GEIp<-round(100*GEI$Sum.Sq/sum(GEI$Sum.Sq[1:nr]),4)

      
      for(i in c(1,3:4,8)) GEI[nr+2,i]<-sum(GEI[1:nr,i])
      row.names(GEI)[nr+2]<-'PCtl'
      GEI['PCtl','Mean.Sq']=round(GEI['PCtl','Sum.Sq']/GEI['PCtl','Df'],4)
      GEI['PCtl','F.value']=round(GEI['PCtl','Mean.Sq']/MS.res,4)
      nr1<-nrow(GEI)
      return(GEI[c(1:pcN,nr1-1,nr1),])
    } 

    PCs.aov<-PC.res(model) # for GEI aov
    PCs.aov
##       percent  acum  Df     Sum.Sq   Mean.Sq F.value   Pr.F     GEIp
## PC1      56.3  56.3  31 13368.5954 431.24501 11.6500 0.0000  56.2609
## PC2      27.1  83.3  29  6427.5799 221.64069  5.9900 0.0000  27.0501
## PCres    16.6 100.0  75  3965.6195  52.87490  1.4278 0.0192  16.6891
## PCtl    100.0    NA 135 23761.7949 176.01330  4.7531     NA 100.0001

查看ASV和YSI结果:

    (Idx2<-index.AMMI(model)) # for ASV
##           ASV YSI rASV rYSI    means
## G1  2.5131813  42   19   23 26.31947
## G10 1.1739456  28    6   22 26.34039
## G11 2.4615460  32   18   14 30.58975
## G12 2.0368354  29   11   18 28.17335
## G13 1.5167528  18    9    9 35.32583
## G14 4.8741897  30   27    3 38.75767
## G15 2.3623790  37   16   21 26.34808
## G16 2.0954103  36   12   24 26.01336
## G17 3.6050812  51   26   25 23.84175
## G18 2.3436592  23   15    8 36.11581
## G19 0.5966506  13    3   10 34.05974
## G2  1.3792025  21    8   13 31.28887
## G20 0.2026430  20    1   19 27.47748
## G21 2.7709324  36   20   16 28.98663
## G22 2.1722949  26   14   12 32.68323
## G23 0.9507170  11    4    7 36.19020
## G24 2.3663500  23   17    6 36.19602
## G25 0.5646275  13    2   11 33.26623
## G26 2.1652861  33   13   20 27.00126
## G27 5.5374138  56   28   28 16.15569
## G28 3.3470545  27   25    2 39.10400
## G3  1.7912464  25   10   15 30.10174
## G4  3.1531170  24   23    1 39.75624
## G5  2.8907699  26   21    5 36.95181
## G6  3.0764673  49   22   27 21.41747
## G7  1.2740344  33    7   26 22.98480
## G8  1.0521529  22    5   17 28.66655
## G9  3.3065468  28   24    4 38.63477

查看各基因型的PCs得分:

    head(model$biplot[,1:5],8)
##     type    Yield       PC1          PC2        PC3
## G1   GEN 26.31947 -1.508289 -1.258765244  0.1922031
## G10  GEN 26.34039 -0.799485 -0.220768053  0.9853880
## G11  GEN 30.58975 -1.495438  1.186549449 -0.9255252
## G12  GEN 28.17335  1.393354  0.332786322  0.7322688
## G13  GEN 35.32583  1.051708 -0.002555823  0.8156191
## G14  GEN 38.75767  3.083381 -1.995946966 -0.8797167
## G15  GEN 26.34808 -1.557371 -0.732314249  0.4143257
## G16  GEN 26.01336 -1.358809  0.741980068 -0.874801102 

结果可视化可以使用plot()和AMMI.contour()来生成双标图:

    plot(model)
    AMMI.contour(model,distance=0.7,shape=8,col="red",lwd=2,lty=5) 

在这里插入图片描述
图1 AMMI模型的双标图

对于GEI的结果可视化,不如GGE方法来得直观。后一篇博文将介绍如何使用R来可视化GGE方法。

03 方法局限性

3.1 AMMI 模型的局限性

只限定于固定效应模型;假定各地点误差同质;要求数据是平衡的。

3.2 AMMI 模型可视化的局限性

未能充分展示基因型的稳定性和高产性;未能充分展示试验地点间的关系。

04 参考文献


  1. Finlay K W, Wilkinson G N. The analysis of adaptation in a plant breeding programme [J]. Aust J Agric Res, 1963, 14: 742-754 ↩︎

  2. Gauch H G, Zobel R W. Identifying mega-environments and targeting genotypes [J]. Crop Sci, 1997, 37: 311-326 ↩︎

  3. Yan W. GGEbiplot—a Windows application for graphical analysis of multi-environment trial data and other types of two-way data [J]. Agron J, 2001, 93:1111-1118 ↩︎

  4. 程玲,张心菲, 张鑫鑫,等. 基于BLUP和GGE双标图的林木多地点试验分析[J]. 西北农林科技大学学报(自然科学版), 2018,46(3):87-93 [return] ↩︎

  • 0
    点赞
  • 11
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值