|||
以5年生的树高h5为目标性状,进行单性状分析。将区组Rep视为固定效应,家系Fam和小区Plot视为随机效应,分析代码如下:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 | ############ 代码清单 ########## setwd( “d:\asreml_data\data” ) # 指定文件路径 library(asreml) # 载入asreml程序包 df <- asreml.read.table( file = ’fm.csv’, header = T, sep = ’,’ ) #读入数据 # names(df); head(df); str(df); summary(df) # 数据集结构 # 分析模型如下: fm <- asreml( h5 ~ 1 + Rep, # 固定效应 random = ~ Fam + Plot, # 随机效应 data = df, # 目标数据集 subset = Spacing == '3' , # 目标数据选择 maxit = 30 # 最大迭代次数 ) # 结果提取命令: plot(fm) # 查看数据是否合理 wald(fm) # 查看固定效应中因子的显著性 summary(fm)$varcomp # 查看方差分量 coef(fm)$random # 查看随机效应值 coef(fm)$fixed # 查看固定效应值 |
ASReml-r基本分析流程如下:
一、判断模型运行后,似然值是否收敛?如未收敛,加大maxit值,或者修改模型。
1 2 3 4 5 6 7 8 9 10 11 12 13 | > fm<-asreml(h5 ~1+ Rep, random =~ Fam + Plot, subset = Spacing == '3', data = df ) asreml 3.0-1 (31 August 2012), Library: 3.0hj (15 November 2011), IA32 LogLik S2 DF wall cpu -2671.0284 5274.2040 551 15:45:43 0.0 (1 restrained) -2668.2939 5295.5663 551 15:45:43 0.0 (1 restrained) -2667.7550 5319.8774 551 15:45:43 0.0 (1 restrained) -2667.7138 5330.6550 551 15:45:43 0.0 (1 restrained) -2667.7114 5332.7691 551 15:45:43 0.0 (1 restrained) -2667.7112 5333.2179 551 15:45:43 0.0 -2667.7112 5333.3478 551 15:45:43 0.0 -2667.7112 5333.3422 551 15:45:43 0.0 Finished on: Sun Aug 11 15:45:43 2013 LogLikelihood Converged |
运行的迭代结果显示,进行了9次迭代后,似然值(LogLikelihood)收敛。
二、判断试验因变量数据是否合理?
对于线性模型,因变量数据应当满足正态性、等方差性、线性、独立性的原则。Asreml-r通过plot()函数可以作出四幅图(图1),用于正态性、线性的判断。上半部分的2幅图是判断因变量数据是否成正态分布,左上角的柱形图应呈正态分布,右上角的图形是QQ图,QQ图中的点应落在45oC角度的直线。下半部分的2幅图是判断因变量与自变量是否成线性关系,图中的点应随机分布于直线两旁。上述的4幅图,均表明,试验数据符合正态性和线性的原则。对于等方差性、独立性的判断,读者可自行采用前文线性回归诊断的方法进行验证。
> plot(fm)
图1 残差图
三、判断固定效应中的因子是否显著(F检验)?如不显著,剔除后,重新运行模型。
1 2 3 4 5 6 7 8 9 10 | > wald(fm) Wald tests for fixed effects Response: h5 Terms added sequentially; adjusted for those above Df Sum of Sq Wald statisticPr(Chisq) (Intercept) 1 103078397 19327.2 < 2.2e-16 *** Rep 4 1497380 280.8 < 2.2e-16 *** residual (MS) 5333 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1 |
对于固定效应,F统计结果表明,因子Rep的F值= 280.8,P值= 2.2e-16 < 0.001,显示区组Rep的效应极显著。
四、判断随机效应中的因子是否显著?
简单判断:z.ratio >= 1.5,就认定该因子效应显著。如有不显著的因子时,去掉后,重新运行模型。更为精确的检验方法是LRT检验。
1 2 3 4 5 | > summary(fm)$varcomp gamma component std.error z.ratio constraint Fam!Fam.var 8.287372e-02 4.419939e+02 1.871792e+02 2.361341 Positive Plot!Plot.var 1.011929e-07 5.396963e-043.412361e-05 15.815922 Boundary R!variance 1.000000e+00 5.333342e+03 3.372135e+0215.815922 Positive |
方差分量提取结果表明,因子Plot为Boundary,即其方差分量值非常小或负值,需舍弃。而家系Fam和误差R均为positive,且z.ration均大于1.5,显示家系和误差的方差分量显著。
因此,去掉不显著的因子Plot,重新运行模型(修改前后模型的迭代情况类似),并获取方差分量,结果如下:
1 2 3 4 | >summary(fm)$varcomp gamma component std.error z.ratio constraint Fam!Fam.var0.08287373 441.9941 187.1784 2.361353 Positive R!variance 1.00000000 5333.3441 337.2136 15.815920 Positive |
五、计算单株遗传力
遗传力(heritability)是重要的遗传参数之一,可简单理解为亲本性状遗传给子代的能力。在本例中,家系为半同胞家系,考虑到单株遗传力比家系遗传力更有实际应用价值,因此,只计算单株遗传力(individual heritability)。
方法一、手动计算遗传力:
h2i=4*Vf/(Vf+Ve)=4*441.9941/(441.9941+5333.3441)= 0.306
方法二、通过编程计算遗传力及其误差
1 2 3 | pincalc<-dget("d:/pin.R") # 载入pin函数 summary(fm)$varcomp[,1:3] # 提取方差分量 pincalc(fm, h2 ~ 4 * V1/(V1 + V2)) # 计算遗传力 |
运行结果如下:
>pincalc(fm, h2 ~ 4 * V1/(V1 + V2))
Estimate SE
h2 0.3061252 0.1239232
六、提取育种值
育种值(breeding value)是另一重要的遗传参数,是决定数量性状的基因加性效应值。从理论上讲,育种值是能100%的遗传给下代的,但毕竟是根据表型值进行间接估计推导出来的,所以育种值也称为估计育种值。计算育种值的目的, 是预测选择育种的效果。ASReml利用BLUP方法可以获得较精确的育种值。
ASReml-r提取育种值的部分结果如下:
> coef(fm)$random
effect
Fam_70001 -4.8627845
Fam_70002 7.0716240
Fam_70003 -3.2965300
Fam_70004-13.4163219
Fam_70005 -2.9936053
Fam_70006 12.2379264
Fam_70007-17.8974327
Fam_70008 -5.3405846
Fam_70009-14.2814921
Fam_70010 -4.2185406
Fam_70011 14.2496541
Fam_70012 -8.5259029
Fam_70015 -3.2158516
Fam_70016 -1.4967663
Fam_70017 8.9747632
Fam_70018 -4.4080512
……
最后,pin()函数代码如下:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 | ###### pin()函数代码 ###### function (object, transform) { pframe <- as.list(object$gammas) names(pframe) <- paste("V", seq(1, length(pframe)), sep = "") tvalue <- eval(deriv(transform[[length(transform)]], names(pframe)), pframe) X <- as.vector(attr(tvalue, "gradient")) X[object$gammas.type == 1] <- 0 tname <- if (length(transform) == 3) transform[[2]] else "" n <- length(pframe) i <- rep(1:n, 1:n) j <- sequence(1:n) k <- 1 + (i > j) Vmat <- object$ai se <- sqrt(sum(Vmat * X[i] * X[j] * k)) data.frame(row.names = tname, Estimate = tvalue, SE = se) } |
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-10-19 21:40
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社