育种数据分析之放飞自我分享 http://blog.sciencenet.cn/u/yijiaobai 关注:生物统计,数量遗传,混合线性模型,生物信息,R,Perl,Python,GWAS,GS相关方法,文章及代码

博文

GCTA学习7 | 计算单性状遗传力和标准误

已有 5495 次阅读 2022-1-18 21:12 |个人分类:GWAS|系统分类:科研笔记

前面的几节中(请翻看博客之前的博文),我们介绍了GCTA计算G矩阵,本节我们介绍,如果使用GCTA进行遗传力的估计。


1. GCTA计算单性状遗传力常用参数

1.1 --reml(必须)

这部分,是使用reml的方法进行估计方差组分。默认的是AI算法,可以使用EM算法。

1.2 --reml-alg 指定迭代方法(非必须)

--reml-alg 0 # AI算法,默认算法 --reml-alg 1 # EM算法

1.3 --reml-priors 迭代初始值

指定迭代初始值,当数据量较大时,指定初始值,会加快迭代速度。

--reml-priors 0.45 0.55

表示Va为0.45, Ve为0.55

1.4 --grm(必须)

指定GRM矩阵

--grm # 接二进制文件GRM的前缀 --grm-bin # 同上 --grm-gz # 接文本的GRM文件前缀

推荐使用二进制的文件,因为速度快,暂用空间少。

1.5 --covar(非必须)

这是接因子协变量的,第一列和第二列分别是FID和IID,后面接因子协变量,比如场年季

1.6 --qcovar(非必须)

接的是数字协变量,比如PCA,比如初生重等

1.7 --pheno(必须)

接的是表型数据,格式也是plink的格式,第一列FID,第二列IID,第三列是表型数据(缺失用NA表示)

如果是多个表型,想指定第四列为表型,可以用--mpheno 2,表示第四列为分析的性状。

1.8 ----mpheno 2 表型数据第四列

如果表型数据中有多列,可以设置第四列为分析的性状。

1.9 --reml-pred-rand 输出育种值BLUP

加上这个代码,会输出BLUP值,GBLUP育种值。

1.10 --blup-snp 输出SNP效应值

再加上--reml-pred-rand--bfile,加上plink二进制文件和育种值信息,会计算SNP的效应值,类似rrblup的值。

1.11 --blup-snp 输出SNP效应值

再加上--reml-pred-rand--bfile,加上plink二进制文件和育种值信息,会计算SNP的效应值,类似rrblup的值。

2. 数据准备

2.1 表型数据

三列,第一列是FID,第二列IID,第三列是表型数据y,没有行头,空格隔开。

2.2 基因型数据

plink的二进制文件

2.3 协变量

这里,示例数据中,没有提供协变量信息。如果提供,可以按照第一列是FID,第二列是IID,其它是协变量的方法整理数据。协变量分为数字协变量和因子协变量,要分开整理。

3. 构建GRM矩阵

3.1 使用Yang的方法

这里,默认的是Yang的方法。

gcta64 --bfile ../test --make-grm --make-grm-alg 0 --out g1

结果文件:

3.2 使用Van的方法

这里,用Van的方法,类似我们GBLUP估计所用的矩阵构建形式。

gcta64 --bfile ../test --make-grm --make-grm-alg 1 --out g2

结果文件:

4. 单性状估算遗传力和标准误

这里,已经构建好了GRM矩阵,指定表型数据,进行遗传力的估计

4.1 使用Yang的GRM矩阵

gcta64 --grm g1 --pheno ../test.phen --reml --out re1

结果可以看出:

  • Vg:加性方差组分为0.022
  • Ve:残差方差组分为0。969
  • Vp:表型方差组分(Vg+Ve),为0.991

遗传力为0.02,标准误是0.008

4.2 使用Van的GRM矩阵

gcta64 --grm g2 --pheno ../test.phen --reml --out re2

结果可以看出:

  • Vg:加性方差组分为0.022
  • Ve:残差方差组分为0。97
  • Vp:表型方差组分(Vg+Ve),为0.991

遗传力为0.02,标准误是0.008

两种方法,结果基本一致。

4.3 使用ASReml作为对比

将二进制的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一致。

5. GBLUP育种值计算和比较

5.1 GCTA计算GEBV值

「代码:」

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。

5.2 ASReml计算BLUP值

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)

5.3 GCTA与ASReml的GEBV进行比较:

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相关的知识。




https://blog.sciencenet.cn/blog-2577109-1321626.html

上一篇:GCTA学习5 | GCTA计算PCA及可视化
下一篇:GCTA学习8 | 计算多性状遗传力和遗传相关
收藏 IP: 223.90.189.*| 热度|

0

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

数据加载中...

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

GMT+8, 2024-10-11 17:03

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部