博文
双性状分析之ASReml-R篇
|||
数据集与单性状的一样，以心材密度dj和5年生树高h5为目标性状，进行双性状的分析。
在做双性状分析前，一般先对各性状进行独自的单性状分析，获得各自的最佳模型，得到各性状的加性遗传方差和误差方差，对于后续双性状分析的模型，可以用于R结构和G结构的初始值设置。
asreml-r双性状的分析模型如下:
1 2 3 4 5 6 7 8 9 10 11 12 | ############ 代码清单 ########## library(asreml) fm.2 <- asreml( cbind( dj, h5 ) ~ trait + trait : Rep, # 固定效应 random =~ us(trait) : Fam， # G结构 rcov =~ units : us(trait), # R结构 maxit = 100, # 最大迭代次数 data = df, # 数据集 ) summary(fm.2)$varcomp wald(fm.2) coef(fm.2)$random |
分析结果如下:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 | > summary(fm.2)$varcomp # 方差分量 gamma component std.error z.ratio trait:Fam!trait.dj:dj 6.392972e-05 6.392972e-05 2.209498e-05 2.893405 trait:Fam!trait.h5:dj 7.547191e-02 7.547191e-02 4.645840e-02 1.624505 trait:Fam!trait.h5:h5 4.560441e+02 4.560441e+02 1.897371e+02 2.403558 R!variance 1.000000e+00 1.000000e+00 NA NA R!trait.dj:dj 4.978785e-04 4.978785e-04 3.142438e-05 15.843704 R!trait.h5:dj -2.129037e-01 -2.129037e-01 7.361576e-02 -2.892093 R!trait.h5:h5 5.325492e+03 5.325492e+03 3.364866e+02 15.826759 constraint trait:Fam!trait.dj:dj Positive trait:Fam!trait.h5:dj Positive trait:Fam!trait.h5:h5 Positive R!variance Fixed R!trait.dj:dj Positive R!trait.h5:dj Positive R!trait.h5:h5 Positive > wald(fm.2) # 固定效应显著性分析 Wald tests for fixed effects Response: y Terms added sequentially; adjusted for those above Df Sum of Sq Wald statistic Pr(Chisq) trait 2 70667 70667 < 2.2e-16 *** trait:Rep 8 373 373 < 2.2e-16 *** residual (MS) 1 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > coef(fm.2)$random # 随机效应值 effect trait_dj:Fam_70001 -8.291279e-03 trait_dj:Fam_70002 -1.895423e-03 trait_dj:Fam_70003 -8.473590e-04 trait_dj:Fam_70004 -3.198284e-03 trait_dj:Fam_70005 3.894355e-04 trait_dj:Fam_70006 1.567450e-03 trait_dj:Fam_70007 4.936646e-03 trait_dj:Fam_70008 3.264833e-03 trait_dj:Fam_70009 -3.500462e-03 trait_dj:Fam_70010 -8.938504e-03 …… trait_h5:Fam_70001 -1.180198e+01 trait_h5:Fam_70002 4.697574e+00 trait_h5:Fam_70003 -3.599515e+00 trait_h5:Fam_70004 -1.443880e+01 trait_h5:Fam_70005 -2.290498e+00 trait_h5:Fam_70006 1.206186e+01 trait_h5:Fam_70007 -1.117458e+01 trait_h5:Fam_70008 -1.570330e+00 trait_h5:Fam_70009 -1.546446e+01 trait_h5:Fam_70010 -1.006793e+01 …… |
有了性状的加性遗传方差与协方差，以及误差方差与协方差，就很容易计算各性状的遗传力，以及性状间的遗传相关与表型相关。由于表型相关在实际应用中指导价值不大，所以一般都只计算遗传相关。
遗传力和遗传相关计算的代码如下:
1 2 3 4 5 6 7 8 9 10 11 12 | pincalc<-dget( "d:/pin.R" ) # 载入pin()函数 summary(fm.2)$varcomp[ ,1:3] #遗传力计算 pincalc( fm.2, h2_A ~ 4 * V1 / ( V1 + V5 )) # A性状的个体遗传力 pincalc( fm.2, h2_B ~ 4 * V3 / ( V3 + V7 )) # B性状的个体遗传力 #遗传相关计算 pincalc( fm.2, gCORR ~ V2 / sqrt( V1 * V3 )) #A、B性状的遗传相关 #A、B性状的表型相关 pincalc( fm.2, pCORR ~ ( V2 + V6 ) / sqrt(( V1 + V5 )* ( V3 + V7 ))) |
运行结果如下:
1 2 3 4 5 6 7 8 9 10 11 12 | > pincalc(fm.2, h2_A ~ 4 * V1/(V1 + V5)) Estimate SE h2_A 0.4551711 0.1450985 > pincalc(fm.2, h2_B ~ 4 * V3/(V3 + V7)) Estimate SE h2_B 0.3155176 0.1251766 > pincalc(fm.2, gCORR ~ V2/sqrt(V1 * V3)) Estimate SE gCORR 0.4420084 0.2573905 > pincalc(fm.2, pCORR ~ (V2+V6)/sqrt((V1+V5)*(V3+V7))) Estimate SE pCORR -0.07625554 0.04494997 |
从运行的结果可知，A性状dj的遗传力为0.455 +-0.145，B性状h5的遗传力为0.316 +- 0.125，遗传相关为0.442 +- 0.257，表型相关为-0.076+- 0.044。相关的显著性，可以通过相关值和误差的大小进行判断。
https://blog.sciencenet.cn/blog-1114360-735461.html
上一篇：单性状分析之ASReml-R篇
下一篇：R高级绘图之lattice篇
