科学网—标签

相关帖子	版块	作者	回复/查看	最后发表

yzhlinscau 2019-9-17 21:09

ASReml是遗传评估的先锋软件，但其R版本在R语言中的运行，受R语言的影响，有时速度上确实会有些慢。Python在数据量比较大时，运行速度确实快过R语言。所以，我也想试试，能否在Python中运行asreml-R程序。答案是可行。不过，R与Python两种语言，在语法和编程上有所区别，仍需继续摸索。具体示例请见链接： https://blog.csdn.net/yzhlinscau/article/details/100939096

个人分类: ASReml|2161 次阅读|0 个评论

2019 R与ASReml-R动植物遗传参数评估与全基因组选择培训会

vsnc 2019-3-7 09:47

ASReml是拟合线性混合模型的优秀数据分析软件，由英国洛桑试验站与 NSW Department of Primary Industries的Arthur Gilmour博士开发，可利用灵活的混合线性模型和广义线性模型来处理大规模的数据，实现大数据高效、快速的分析。ASReml涵盖育种等专业领域所需复杂模型，可对数量性状、阈值性状、分类性状、SNP标记等多个性状进行分析，并可实现固定效应、随机效应值的预测以及显著性检验、遗传参数的估计、全基因组选择等。目前，ASReml软件已在世界范围内广泛地应用于动物、水产、林业、作物以及医学等多个领域的研究。为了给ASReml用户以及对遗传评估、混合线性模型感兴趣的科研人员，提供一个交流的平台，ASReml软件的中国子公司VSNC特在北京举办“R与ASReml-R动植物遗传参数评估与全基因组选择培训会”，并现场举办ASReml国际认证考试，通过的学员将会获得ASReml国际资格认证证书。欢迎各领域的研究人员参会。会议地点：北京上园饭店（会议中心）会议时间： 2019年3月30日-3月31日培训软件： ASReml-R4 培训优势：理论知识与实际案例相结合，提高学员解决实际问题的能力提供ASReml-R4软件，提供国际认证考试与证书培训日程注册费：参会须知：会务组可协助安排住宿，费用自理；培训期间免费提供午餐；请各位老师同学自备笔记本电脑；报名付费后工作人员会为各位老师及同学申请并发送软件（邮件），同时协助参会人员在会前安装，请各位老师注意查收邮件。由于软件绑定电脑，请各位参会者携带安装好软件的电脑参加培训。报名及付费方式：下载附件《报名回执》，报名回执表请发送至邮箱：China@vsni.co.uk，并将报名费汇入以下账户（报名时间以收到汇款时间为准）：开户名：北京维斯恩思软件有限责任公司开户行：中国建设银行北京中关村分行帐号：1100 1007 3000 5301 7767 注：汇款时请务必注明单位、姓名（例如：中国农业大学张三注册费）；发票及通知（加盖公章）将于会议当天统一发放。会务联系人：张娟（13121623804 / 010-88400822 / 010-62680244）邮箱：China@vsni.co.uk 附件：第二轮培训通知及报名回执 R与ASReml-R动植物遗传参数评估与全基因组选择培训会第二轮通知-2019.03.docx

1487 次阅读|0 个评论

R语言中计算blue和blup值的包：lme4、nlme、MCMCglmm、asreml

yijiaobai 2016-1-14 13:39

数据： tree.rda 利用同一数据，演示不同软件包：lme4、nlme、MCMCglmm和asreml估计BLUE值和BLUP值的代码。固定效应：Block 随机效应：Family 代码： ###############lme4中运行blue和blup的方法############### load(tree.rda) library(lme4) df-lme - lmer(Height~1+Block+(1|Family),data=tree) print(df_lme) anova(df_lme) ####求方差 ranef(df_lme) ####求随机效应的BLUP值 fixef(df_lme) ####求固定效应的BLUE值 ####################nlme包中运行的blue和blup值####################### library(nlme) df_nlme - lme(Height~1+Block,random = ~1|Family,data=tree) print(df_nlme) anova(df_nlme) random.effects(df_nlme) fixed.effects(df_nlme) ###############以上是mcmc包的结果，结果好像不符合blup和blue值################### library(MCMCglmm) df_mcmc - MCMCglmm::MCMCglmm(Height~1+Block,random=~Family,pr = T,family = 'gaussian',data=tree) summary(df_mcmc) posterior.mode(df_mcmc$VCV) posterior.mode(df_mcmc$Sol) #######################以上是asreml软件包计算的结果########### library(asreml) df_asreml - asreml(Height ~1+Block,random =~ Family,data = tree) summary(df_asreml)$varcomp wald(df_asreml) coef(df_asreml) 参考资料：童春发. 林木遗传模型统计分析及R语言实现 . 科学出版社, 2014. 林元震. R与ASReml-R统计分析教程 . 中国林业出版社, 2014.

个人分类: 农学统计|23162 次阅读|0 个评论

ASReml V4.0的新特性

热度 1 yzhlinscau 2014-1-7 10:59

ASReml的软件开发者Arthur Gilmour一直在坚持软件的不断升级和更新，对于一位年过60的人来说，实属不易。ASReml现已开发至4.0版。自4.0版起，基本不用设置初始值，先前的版本看起来比较professional，现在的版本相比而言，似乎更傻瓜化，犹如单反相机与数码相机。ASReml4.0的代码格式，更趋于asreml-r的格式。举个简单的示范，比较下V3.0和V4.0的差异。 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 !PART 1 # Random regression, V3.0 !continue !maxit 30 !STEP 0.1 root ~ mu Age Age.Block !r pol(Age,2).Clone !f mv 1 2 1 Age 0 US !+10 !GP !S2==1 10*0 1400 pol(Age,2).Clone 2 pol(Age,2) 0 US !+6 !GP 6*0 Clone !PART 2 # V4.0 !continue !maxit 30 root ~ mu Age Age.Block !r us(pol(Age,2)).Clone !f mv residual us(Age).id(1400) 从代码中可知，V4.0的代码简洁多了，是否更像数码相机，对于想入门ASReml的人来说，不是更酷吗？此外，ASReml还真可分析非常复杂的模型，我发现了asreml-r的一个模型可以运行，但在ASReml win版本中，不知如何运行，特意请教了Arthur Gilmour，他非常及时地给予了反馈。相对而言，asreml-r的代码最为简洁，因此建议初学者，还是先学习asreml-r版本为好。现把代码粘贴如下： 1 2 3 4 5 6 7 8 9 10 11 12 13 # asreml-r code df.asr-asreml(cbind(root,branch)~trait+trait:Age+trait:Age:Block, random=~diag(trait):Age:Clone, rcov=~at(Age):ar1(Row):ar1(Col):corh(trait), maxit=100,data=df) # asreml v4.0 （新版本尚未发布，目前的版本无法运行！！） root, branch ~ Trait Trait.Age Trait.Age.Block mv, !r xfa1(Trait).Age.Clone residual at(Age).ar1(Row).ar1(Col).us(Trait)

个人分类: ASReml|4372 次阅读|1 个评论

双性状分析之ASReml-R篇

yzhlinscau 2013-10-23 22:26

数据集与单性状的一样，以心材密度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-052.209498e-05 2.893405 trait:Fam!trait.h5:dj 7.547191e-02 7.547191e-024.645840e-02 1.624505 trait:Fam!trait.h5:h5 4.560441e+02 4.560441e+021.897371e+02 2.403558 R!variance 1.000000e+00 1.000000e+00 NA NA R!trait.dj:dj 4.978785e-04 4.978785e-043.142438e-0515.843704 R!trait.h5:dj -2.129037e-01-2.129037e-017.361576e-02-2.892093 R!trait.h5:h5 5.325492e+03 5.325492e+033.364866e+0215.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 forfixed effects Response: y Terms added sequentially; adjusted forthose 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 #遗传力计算 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.45517110.1450985 pincalc(fm.2, h2_B ~ 4* V3/(V3 + V7)) Estimate SE h2_B 0.31551760.1251766 pincalc(fm.2, gCORR ~ V2/sqrt(V1 * V3)) Estimate SE gCORR 0.44200840.2573905 pincalc(fm.2, pCORR ~ (V2+V6)/sqrt((V1+V5)*(V3+V7))) Estimate SE pCORR -0.076255540.04494997 从运行的结果可知， A 性状 dj 的遗传力为 0.455 +- 0.145 ， B 性状 h5 的遗传力为 0.316 +- 0.125 ，遗传相关为 0.442 +- 0.257 ，表型相关为 - 0.076 +- 0.044 。相关的显著性，可以通过相关值和误差的大小进行判断。

个人分类: ASReml|4710 次阅读|0 个评论

单性状分析之ASReml-R篇

热度 1 yzhlinscau 2013-10-23 22:22

以 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 ) 。方法一、手动计算遗传力： h 2 i =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 # 提取方差分量 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 ], names(pframe)), pframe) X - as.vector(attr(tvalue, gradient)) X - 0 tname - if (length(transform) == 3) transform ] 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 * X * k)) data.frame(row.names = tname, Estimate = tvalue, SE = se) }

个人分类: ASReml|9598 次阅读|3 个评论

帐号		自动登录	找回密码
密码			注册

关闭安全验证

标签: ASReml

相关帖子

相关日志

关闭 安全验证

标签: ASReml

相关帖子

相关日志

关闭安全验证