||
前面的几节中(请翻看博客之前的博文),我们介绍了GCTA计算G矩阵,本节我们介绍,如果使用GCTA进行遗传力的估计。
这部分,是使用reml的方法进行估计方差组分。默认的是AI算法,可以使用EM算法。
--reml-alg 0 # AI算法,默认算法 --reml-alg 1 # EM算法
指定迭代初始值,当数据量较大时,指定初始值,会加快迭代速度。
--reml-priors 0.45 0.55
表示Va为0.45, Ve为0.55
指定GRM矩阵
--grm # 接二进制文件GRM的前缀 --grm-bin # 同上 --grm-gz # 接文本的GRM文件前缀
推荐使用二进制的文件,因为速度快,暂用空间少。
这是接因子协变量的,第一列和第二列分别是FID和IID,后面接因子协变量,比如场年季
接的是数字协变量,比如PCA,比如初生重等
接的是表型数据,格式也是plink的格式,第一列FID,第二列IID,第三列是表型数据(缺失用NA表示)
如果是多个表型,想指定第四列为表型,可以用--mpheno 2
,表示第四列为分析的性状。
如果表型数据中有多列,可以设置第四列为分析的性状。
加上这个代码,会输出BLUP值,GBLUP育种值。
再加上--reml-pred-rand
和--bfile
,加上plink二进制文件和育种值信息,会计算SNP的效应值,类似rrblup的值。
再加上--reml-pred-rand
和--bfile
,加上plink二进制文件和育种值信息,会计算SNP的效应值,类似rrblup的值。
三列,第一列是FID,第二列IID,第三列是表型数据y,没有行头,空格隔开。
plink的二进制文件
这里,示例数据中,没有提供协变量信息。如果提供,可以按照第一列是FID,第二列是IID,其它是协变量的方法整理数据。协变量分为数字协变量和因子协变量,要分开整理。
Yang
的方法这里,默认的是Yang
的方法。
gcta64 --bfile ../test --make-grm --make-grm-alg 0 --out g1
结果文件:
Van
的方法这里,用Van
的方法,类似我们GBLUP估计所用的矩阵构建形式。
gcta64 --bfile ../test --make-grm --make-grm-alg 1 --out g2
结果文件:
这里,已经构建好了GRM矩阵,指定表型数据,进行遗传力的估计
Yang
的GRM矩阵gcta64 --grm g1 --pheno ../test.phen --reml --out re1
结果可以看出:
遗传力为0.02,标准误是0.008
Van
的GRM矩阵gcta64 --grm g2 --pheno ../test.phen --reml --out re2
结果可以看出:
遗传力为0.02,标准误是0.008
两种方法,结果基本一致。
将二进制的GRM文件,变为asreml支持的ginv格式。
「asreml代码:」
mod = asreml(V3 ~ 1, random = ~ vm(V2,ginv), residual = ~ idv(units), workspace = "10Gb", dense = ~ vm(V2,ginv), data=dat) summary(mod)$varcomp vpredict(mod,h2 ~ V1/(V1+V2))
「方差组分和遗传力结果:」
结果和GCTA一致。
「代码:」
gcta64 --bfile ../test --make-grm --make-grm-alg 1 --out g2 gcta64 --grm g2 --pheno ../test.phen --reml --out re2 --reml-pred-rand
计算的育种值:re2.indi.blp
,第四列为GEBV。
mod = asreml(V3 ~ 1, random = ~ vm(V2,ginv), residual = ~ idv(units), workspace = "10Gb", dense = ~ vm(V2,ginv), data=dat) summary(mod)$varcomp vpredict(mod,h2 ~ V1/(V1+V2)) blup = coef(mod)$random %>% tiqu_blup() head(blup)
blup = coef(mod)$random %>% tiqu_blup() head(blup) aa = fread("../10_reml_van_uni/re2.indi.blp") head(aa) re = merge(blup,aa,by.x = "ID",by.y = "V2") head(re) re %>% select(effect,V4) %>% cor
可以看出,结果完全一样。
❝欢迎关注我的公众号:
❞育种数据分析之放飞自我
。主要分享R语言,Python,育种数据分析,生物统计,数量遗传学,混合线性模型,GWAS和GS相关的知识。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-10-11 17:03
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社