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

博文

单性状分析之ASReml-R篇

已有 9297 次阅读 2013-10-23 22:22 |个人分类:ASReml|系统分类:科研笔记| ASReml, ASReml-R, 单性状分析

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统计结果表明,因子RepF= 280.8P= 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

方差分量提取结果表明,因子PlotBoundary,即其方差分量值非常小或负值,需舍弃。而家系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)
}




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


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

0

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

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

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

GMT+8, 2024-3-29 15:51

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部