yzhlinscau的个人博客分享 http://blog.sciencenet.cn/u/yzhlinscau

博文

G×E分析之ASReml-R篇

已有 5794 次阅读 2013-10-25 13:59 |个人分类:ASReml|系统分类:科研笔记

     为研究种植密度对辐射松生长性状的影响,以55个半同胞家系为试验材料,设置3种种植密度( S1:1×1m,S2:1×2m 和 S3:2×3m),每种植密度下设置5个区组,4株行式小区。试验设计为随机完全区组设计。在28年生的该试验林分,对所有存活树钻取1.3m处的木芯(15mm),共1456个木芯。测量数据分两套:一套是早期生长性状(树高h和胸径dbh),另一套是存活树的生长性状(年轮宽度rw和胸径dbh)。后者的胸径dbh数据是由年轮宽度rw(从软件WinDENDRO分析获取)转换而来的。因本试验林初始种植密度较小,树木个体之间的竞争比较激烈,导致多数木芯中年轮的可数数量不超过21个,因此,本试验的第二套数据仅分析从第2个年轮到年第21个年轮之间的数据。


      本试验的通用线性模型为:

      Y  = μ+ S + R+ F + SF + e

      式中,Y是树木性状的观测值,μ是总体均值,S是种植密度效应,R是内嵌于种植密度的区组效应,F是家系效应,SF是家系与种植密度互作效应,e是残差。其中,μ、S、R为固定效应,其余为随机效应。      

      家系稳定系数的计算公式如下:

     DBHik = ak + bk× DBHi

      式中,DBHik是每个种植密度下的每个家系dbh均值,DBHi 是每个种植密度下的所有家系dbh均值,ak是回归常数,bk是回归系数,即家系稳定系数。稳定系数b用于分析家系与种植密度互作效应。

   
首先,进行各家系年轮宽度rw的分析。代码如下:

1
2
3
4
5
6
7
############# 代码清单  #############
library(asreml)
RW <- asreml.read.table( file = 'RW.csv', header = T, sep = ',' )
RW.asr <- asreml(rw~ca, random =~ diag(pol(ca)):Fam + spl(ca) + spl(ca):Fam,
                      Knots = 12, data = RW, subset = Spacing == 2 )
rw.pv <- predict( RW.asr, classify = "spl(ca):Fam" )
xyplot(predicted.value~ca,data=rw.pv$predictions$pvals, groups=Fam, type="l")

上述代码是用于预测55个家系在第2种植密度下的年轮宽度变化趋势。运行后生成的图形如下:



第二,估算家系单株遗传力。分析代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
############# 代码清单  #############
 
library(asreml)
 
H<-asreml.read.table(file='earlygrowth.csv',header=T, sep = ',') # 读入第一套数据
H2<-asreml.read.table(file='lategrowth.csv',header=T, sep=',') # 读入第二套数据
 
H.asr<-asreml(h2~1+Rep, random=~Fam ,data =H,subset=Spacing==2)  # 树高模型
H.asr2<-asreml(dbh2~1+ Rep,random=~Fam, data=H,subset=Spacing==2) # 胸径模型
 
pincalc <- dget("pin.R")
summary(H.asr)$varcomp[ ,1:3]
pincalc( H.asr, hsq ~ 4 * V1/(V1 + V2))  # 计算树高或胸径单株遗传力

通过上述的程序代码,容易得到各年份树高的单株遗传力,并可绘制出辐射松在3种种植密度下,早期测定的树高单株遗传力的变化趋势图。结果表明,树高早期单株遗传力都比较高,从第二年的0.5下降到第五年的0.4,而且结合遗传力估算值的误差,可知这3种种植密度对树高遗传力没有显著的影响。


      同理,可以分析早期和晚期辐射松胸径的单株遗传力随着树龄的变化趋势(6-11)。辐射松胸径早期和晚期的家系单株遗传力变化要比树高的复杂。早期测定时,由于S1S2种植较密,树木之间的竞争导致遗传力下降。而在晚期测定时,S1由于竞争过于激烈而致使多数树木死亡,进而导致遗传力的估算出现严重偏差。S2S3估算的遗传力相对合理。各种植密度具体对家系遗传力的影响,感兴趣的读者请查看笔者的论文(Lin et al.,For Ecol Manage,2013,304:204-211)



第三,计算胸径的年年遗传相关。

1
2
3
4
5
6
7
8
9
10
########## 代码清单  age-age corr #############
library(asreml)
df <- asreml.read.table( file = 'fm.csv', header = T, sep = ',' )
h.asr <- asreml( cbind(h1,h5) ~ trait + trait:Rep, random =~ us(trait):Fam,
             rcov =~ units:diag(trait), maxit = 100, data = df )
summary(h.asr)$varcomp
wald(h.asr)
pincalc <- dget("pin.R")  # 载入pin()函数
summary(h.asr)$varcomp[ ,1:3]
pincalc( h.asr, r.aa ~ V2/sqrt(V1 * V3) )



家系胸径早晚相关的结果显示,早期测定时,不论哪种种植密度,与第10年生胸径的年年相关都比较高,而且随着树龄增加而加强,表明种植密度对早期胸径的年年相关没有影响。而在晚期测定中,除了S3胸径的年年相关变化曲线正常外,其余的S1S2,胸径的年年相关变化曲线均出现异常的偏差。这说明较激烈的树木竞争和较高的致死率会影响胸径年年相关的变化趋势。


第四,分析家系与种植密度之间是否存在显著的G×E效应。综合3种种植密度,统一分析Vf×s以及Vf×s/Vf比值随着树龄的变化趋势,以判断是否存在G×E效应。

1
2
3
4
5
6
############# 代码清单  #############
H2<-asreml.read.table(file='lategrowth.csv', header=T,sep=',')   # 读入第二套数据
H.asr3<- asreml(dbh2~1+Spa+Spa/Rep,random =~ Fam+Spa:Fam, data=H2 ) # 胸径模型
 
summary(H.asr3)$varcomp # 获取随机效应各变量的方差分量
wald(H.asr3)            # 获取固定效应的F统计检验量

对上述模型运行的结果,计算Vf×s/Vf比值,Vf×s/Vf分别是家系与种植密度互作的方差、家系的方差分量。根据模型的运行结果,Vf×s从第5年开始,一直到第14年,均显著,而家系效应一直不显著。在固定效应中,种植密度效应一直极其显著。现将Vf×s/Vf比值与树龄绘成图。从图中可知,Vf×s/Vf比值先随着树龄的增加而变大,在第10年时,比值达到峰值,而后随着树龄的增加而呈下降的趋势。Vf×s/Vf比值在第10年时最大,而且Vf×s效应显著,这表明家系与种植密度互作效应(F×S)在第10年时最强烈。


家系dbh均值在不同种植密度下的排名发生变化,这是G×E中的秩次改变效应。


但本例中的F×S效应是否由少数家系引起的,还需计算家系的稳定系数b。当b=1时,表明该家系比较稳定,在各种植密度下表现属于平均水平;当b>1时,表明该家系极不稳定,在相对好的环境当中表现优异;当b<1时,表明该家系表现较差但比较稳定,适合于相对差的环境中。



通过上述的表和图可知,本例中的F×S效应是少数表现优异的家系引起的。例如,家系70055的F×S效应比较强烈,其稳定系数b显著高于1,说明其在比较好的环境当中表现优异;而家系70037,其稳定系数b显著小于1,说明其表现虽然差,但比较稳定,适合于比较恶劣的环境。


上述这个案例,就是通过家系的稳定系数来研究G×E效应。




https://blog.sciencenet.cn/blog-1114360-736005.html

上一篇:空间分析之ASReml-R篇
下一篇:双性状分析之MCMCglmm篇
收藏 IP: 58.254.92.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (1 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-11-22 09:11

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部