1. GCTA計(jì)算單性狀遺傳力常用參數(shù)

1.1 --reml（必須）

這部分，是使用reml的方法進(jìn)行估計(jì)方差組分。默認(rèn)的是AI算法，可以使用EM算法。

1.2 --reml-alg 指定迭代方法（非必須）

--reml-alg 0 # AI算法，默認(rèn)算法 --reml-alg 1 # EM算法

1.3 --reml-priors 迭代初始值

指定迭代初始值，當(dāng)數(shù)據(jù)量較大時(shí)，指定初始值，會(huì)加快迭代速度。

--reml-priors 0.45 0.55

表示Va為0.45， Ve為0.55

1.4 --grm（必須）

指定GRM矩陣

--grm # 接二進(jìn)制文件GRM的前綴 --grm-bin # 同上 --grm-gz # 接文本的GRM文件前綴

推薦使用二進(jìn)制的文件，因?yàn)樗俣瓤?，占用空間少。

1.5 --covar（非必須）

這是接因子協(xié)變量的，第一列和第二列分別是FID和IID，后面接因子協(xié)變量，比如場(chǎng)年季

1.6 --qcovar（非必須）

接的是數(shù)字協(xié)變量，比如PCA，比如初生重等

1.7 --pheno（必須）

接的是表型數(shù)據(jù)，格式也是plink的格式，第一列FID，第二列IID，第三列是表型數(shù)據(jù)（缺失用NA表示）

如果是多個(gè)表型，想指定第四列為表型，可以用--mpheno 2，表示第四列為分析的性狀。

1.8 ----mpheno 2 表型數(shù)據(jù)第四列

如果表型數(shù)據(jù)中有多列，可以設(shè)置第四列為分析的性狀。

1.9 --reml-pred-rand 輸出育種值BLUP

加上這個(gè)代碼，會(huì)輸出BLUP值，GBLUP育種值。

1.10 --blup-snp 輸出SNP效應(yīng)值

再加上--reml-pred-rand和--bfile，加上plink二進(jìn)制文件和育種值信息，會(huì)計(jì)算SNP的效應(yīng)值，類(lèi)似rrblup的值。

1.11 --blup-snp 輸出SNP效應(yīng)值

再加上--reml-pred-rand和--bfile，加上plink二進(jìn)制文件和育種值信息，會(huì)計(jì)算SNP的效應(yīng)值，類(lèi)似rrblup的值。

2. 數(shù)據(jù)準(zhǔn)備

2.1 表型數(shù)據(jù)

三列，第一列是FID，第二列IID，第三列是表型數(shù)據(jù)y，沒(méi)有行頭，空格隔開(kāi)。

2.2 基因型數(shù)據(jù)

plink的二進(jìn)制文件

2.3 協(xié)變量

這里，示例數(shù)據(jù)中，沒(méi)有提供協(xié)變量信息。如果提供，可以按照第一列是FID，第二列是IID，其它是協(xié)變量的方法整理數(shù)據(jù)。協(xié)變量分為數(shù)字協(xié)變量和因子協(xié)變量，要分開(kāi)整理。

3. 構(gòu)建GRM矩陣

3.1 使用Yang的方法

這里，默認(rèn)的是Yang的方法。

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

結(jié)果文件：

3.2 使用Van的方法

這里，用Van的方法，類(lèi)似我們GBLUP估計(jì)所用的矩陣構(gòu)建形式。

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

結(jié)果文件：

4. 單性狀估算遺傳力和標(biāo)準(zhǔn)誤

這里，已經(jīng)構(gòu)建好了GRM矩陣，指定表型數(shù)據(jù)，進(jìn)行遺傳力的估計(jì)

4.1 使用Yang的GRM矩陣

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

結(jié)果可以看出：

• Vg：加性方差組分為0.022
• Ve：殘差方差組分為0。969
• Vp：表型方差組分（Vg+Ve），為0.991

遺傳力為0.02，標(biāo)準(zhǔn)誤是0.008

4.2 使用Van的GRM矩陣

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

結(jié)果可以看出：

• Vg：加性方差組分為0.022
• Ve：殘差方差組分為0。97
• Vp：表型方差組分（Vg+Ve），為0.991

遺傳力為0.02，標(biāo)準(zhǔn)誤是0.008

兩種方法，結(jié)果基本一致。

4.3 使用ASReml作為對(duì)比

將二進(jìn)制的GRM文件，變?yōu)閍sreml支持的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))

「方差組分和遺傳力結(jié)果：」

結(jié)果和GCTA一致。

5. GBLUP育種值計(jì)算和比較

5.1 GCTA計(jì)算GEBV值

「代碼：」

gcta64 --bfile ../test --make-grm --make-grm-alg 1 --out g2 
gcta64 --grm g2 --pheno ../test.phen --reml --out re2 --reml-pred-rand

計(jì)算的育種值：re2.indi.blp，第四列為GEBV。

5.2 ASReml計(jì)算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進(jìn)行比較：

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

可以看出，結(jié)果完全一樣。