大家好，我是邓飞，之前介绍了：丢失的遗传力是个什么鬼？

里面的遗传力分为：family遗传力，SNP遗传力和GWAS遗传力，他们的关系如下：

family遗传力可以通过同卵双胞胎、全同胞、半同胞、亲子、甚至F2，BC群体估计。

SNP遗传力是全部SNP的遗传力，混合线性模型中和GBLUP估计遗传力等价，这里我们介绍一下计算方法。

GWAS遗传力是显著SNP的解释百分比，具体可以参考我写的系列博客：GWAS分析中SNP解释百分比PVE | 第四篇，MLM模型中如何手动计算PVE？

这里介绍一下SNP的遗传力如何计算，我们使用的是GCTA和ASReml软件，通过计算比较，可以得到SNP的遗传力就是GBLUP的遗传力的结论。所以，医学和动植物育种在这个概念上，是一致的。

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

可以看出，结果完全一样。