GCTA(全基因組復(fù)雜性狀分析)工具開發(fā)目的是針對(duì)復(fù)雜性狀的全基因組關(guān)聯(lián)分析，評(píng)估SNP解釋的表型方差所占的比例。
(官網(wǎng)：http://cnsgenomics.com/software/gcta/)。
目前GCTA工具可實(shí)現(xiàn)以下功能：

1 評(píng)估全基因組SNP的親緣關(guān)系（遺傳關(guān)系）
2 評(píng)估全基因組SNP的近交系數(shù)
3 評(píng)估所有的常染色體SNP對(duì)于變異的解釋度（遺傳度）
4 評(píng)估遺傳方差與X-染色體的關(guān)聯(lián)
5 檢測(cè)遺傳方差對(duì)X-染色體的劑量補(bǔ)償效應(yīng)
6 預(yù)測(cè)單個(gè)個(gè)體和單個(gè)SNP的全基因組加性遺傳效應(yīng)
7 估計(jì)包含LD結(jié)構(gòu)的目標(biāo)SNP
8 根據(jù)觀察到的基因型數(shù)據(jù)模擬GWAS數(shù)據(jù)
9 轉(zhuǎn)化Illumina原始基因型數(shù)據(jù)為PLINK格式
10 在沒有個(gè)體層面的基因型數(shù)據(jù)下，條件與聯(lián)合分析GWAS的概括統(tǒng)計(jì)量
11 使用SNP數(shù)據(jù)估計(jì)兩個(gè)特性（疾?。┑倪z傳相關(guān)性
12 混合線性模型關(guān)聯(lián)分析

GCTA: a tool for genome-wide complex trait analysis
2011 被引量

下載、安裝軟件

先從下載開始吧，進(jìn)入官網(wǎng)，找到download，如圖

找到適合自己的版本下載吧。解壓完是這樣的。

根據(jù)所有的常染色體snp計(jì)算遺傳關(guān)系矩陣（GRM）

所需要的文件為上圖中的bim，fam，bed文件

./gcta64 --bfile test --autosome --maf 0.01 --make-grm --out test --thread-num 10

結(jié)果生成：

如果要是數(shù)據(jù)量比較大，可以分開染色體計(jì)算，代碼如下：

./gcta64 --bfile test --chr 1 --maf 0.01 --make-grm --out test_chr1 --thread-num 10
./gcta64 --bfile test --chr 2 --maf 0.01 --make-grm --out test_chr2 --thread-num 10
...
./gcta64 --bfile test --chr 22 --maf 0.01 --make-grm --out test_chr22 --thread-num 10

就會(huì)生成這樣的文件：

image.png

然后再把染色體merge到一起

./gcta64 --mgrm grm_chrs.txt --make-grm --out test

去掉隱性相關(guān)的代碼：

./gcta64 --grm test --grm-cutoff 0.025 --make-grm --out test_rm025

cutoff這個(gè)值可以根據(jù)自己的數(shù)據(jù)調(diào)整

計(jì)算遺傳度

--reml
輸入文件：上一步生成的文件、表型文件
表型文件后綴為.txt格式，不需要表頭，第一列為family ID, 第二列為individual ID 第三列為 phenotypes ，類似于PLINK的表型文件格式

./gcta64 --grm test --pheno test.phen --reml --out test --thread-num 10

生成結(jié)果文件：

這一步計(jì)算時(shí)也可以回歸協(xié)變量，自己準(zhǔn)備協(xié)變量文件，示例為PCA前10個(gè)主成分：

./gcta64 --grm test --pheno test.phen --reml --qcovar test_10PCs.txt --out test --thread-num 10

計(jì)算遺傳度也可以一條染色體一條染色體計(jì)算，代碼如下：

./gcta64 --grm test_chr1 --pheno test.phen --reml --out test_chr1 --thread-num 10
./gcta64 --grm test_chr2 --pheno test.phen --reml --out test_chr2 --thread-num 10
......
./gcta64 --grm test_chr22 --pheno test.phen --reml --out test_chr22 --thread-num 10

就會(huì)生成這樣的文件：

或者是把上一步merger之后的一起跑：

./gcta64 --mgrm grm_chrs.txt --pheno test.phen --reml --out test_all_chrs --thread-num 10

test.hsp文件內(nèi)容如下：

V(G)/V(P)就是我們要的遺傳度啦，為0.022347，遺傳度比較低啊，p值為1.3948e-03。

計(jì)算兩個(gè)形狀之間的遺傳相關(guān)性

--reml-bivar
準(zhǔn)備表型文件，后綴為.txt格式，不需要表頭，第一列為family ID, 第二列為individual ID 第三列和第四列為 phenotypes，類似于PLINK的表型文件格式。

代碼：

./gcta64 --reml-bivar --reml-bivar-nocove --grm test --pheno pheno.txt --reml-bivar-lrt-rg 0 --out test

Source Variance SE
V(G)_tr1 0.479647 0.179078 #trait 1 的遺傳方差和標(biāo)準(zhǔn)誤
V(G)_tr2 0.286330 0.181329 #trait 2 的遺傳方差和標(biāo)準(zhǔn)誤
C(G)_tr12 0.230828 0.147958 #trait 1 和 2 之間的遺傳協(xié)方差和標(biāo)準(zhǔn)誤 V(e)_tr1 0.524264 0.176650 #trait 1 的剩余方差和標(biāo)準(zhǔn)誤
V(e)_tr2 0.734654 0.181146 #trait 2 的剩余方差和標(biāo)準(zhǔn)誤
C(e)_tr12 0.404298 0.146863 #trait 1 和 2 的剩余協(xié)方差和標(biāo)準(zhǔn)誤
Vp_tr1 1.003911 0.033202
Vp_tr2 1.020984 0.033800
V(G)/Vp_tr1 0.477779 0.176457
V(G)/Vp_tr2 0.280445 0.176928
rG 0.622864 0.217458 # 遺傳相關(guān)性和標(biāo)準(zhǔn)誤
n 3669 # 樣本量

其中，rG即為我們想要的遺傳相關(guān)性，0.622864 和 0.217458分別代表兩個(gè)性狀/表型間的遺傳相關(guān)性（genetic correlation）和標(biāo)準(zhǔn)誤（Stand error）

————————————————————————
一般遺傳度LDSC也可以算，但兩者的不同是LDSC用summary數(shù)據(jù)計(jì)算，而GCTA用基因型數(shù)據(jù)算。貌似準(zhǔn)點(diǎn)吧。

其他推文：
1.https://mp.weixin.qq.com/s/z4TAoNZjmDSKS6mKY6FgaQ
2.https://www.cnblogs.com/chenwenyan/p/6219749.html