我們上一節(jié)講到了如何處理vcf文件，那我們這節(jié)課給大家講一下，如何運用處理好的vcf文件進(jìn)行群體結(jié)構(gòu)計算。

首先給大家講一下為什么要做群體結(jié)構(gòu)呢，舉一個簡單的例子，黃種人的頭發(fā)是黑的，白種人的頭發(fā)是非黑的，你拿這個性狀在這倆人種里關(guān)聯(lián)，所有一個人種內(nèi)特異的SNP都會起來峰，那這個峰還很好看，那這個峰是不是你想要的呢？

答案是NO，因為這是群體結(jié)構(gòu)造成的影響。具體的原理我們這里不過多的陳述，如果想了解原理，可以看這里的群體結(jié)構(gòu)與親緣關(guān)系原理這一章節(jié)，里面有很詳盡的講解

http://www./class/home/index/series?id=41

我們這期主要是講重測序數(shù)據(jù)的GWAS分析，那重測序數(shù)據(jù)的SNP數(shù)目少則幾十萬，多則幾百萬上千萬，傳統(tǒng)的的計算群體結(jié)構(gòu)的軟件是structure，計算速度十分感人，如果您的項目是幾百個品種，上百萬甚至上千萬個SNP，可能計算到您畢業(yè)都不一定都算完，這里給大家推薦一個算法和structure一樣的，計算速度更快的軟件，那就是今年來引用率更高的admixture。

運用admixture需要準(zhǔn)備的文件，為處理好的vcf文件，我們建議您用admixture進(jìn)行群體結(jié)構(gòu)分析的時候，能將LD近的標(biāo)記過濾掉，只保留這個，這樣會提升您的運算速度，Admixture的美中不足的是不接受vcf文件，你需要將vcf文件轉(zhuǎn)換為admixture所接受的bed格式，這里推薦大家用plink軟件進(jìn)行轉(zhuǎn)換。

##按照LD過濾并轉(zhuǎn)換為bed格式

plink--vcf snp.int0.8.maf0.05.vcf --indep-pairwise 100 50 0.2 --outsnp.int0.8.maf0.05 --allow-extra-chr --make-bed

##轉(zhuǎn)換為admixture可以接受的格式

plink--bfile snp.int0.8.maf0.05 --extract snp.int0.8.maf0.05.prune.in --out prunData--recode 12 --allow-extra-chr

生成的文件就可以做admixture啦。

##做K=2時候的admixture

admixture--cv prunData.ped 2 >> log.txt

運行結(jié)束后會生成響應(yīng)的Q文件：

這個Q文件，稍加修改就可以進(jìn)行GWAS分析了，當(dāng)然您也可以用已經(jīng)做出來的結(jié)果，繪制發(fā)表級別的圖片（如下圖）。

以上部分就是admixture的全部講解了，當(dāng)然我們建議您能充分解群體結(jié)構(gòu)原理，會用admixture結(jié)果繪制發(fā)表級的圖片，一個高分的文章，圖片一定是精美的。

上面給大家講了如何用admixture去計算群體結(jié)構(gòu)，在做GWAS分析的時候，只矯正群體結(jié)構(gòu)是不夠的，親緣關(guān)系也會回GWAS結(jié)果造成一定的假陽性。這一節(jié)為大家講一下如何計算kinship。

先給大家介紹幾款計算親緣關(guān)系的軟件：目前常用的幾款軟件有GCTA、LDAK、SPAGeDi、TASSEL。

首先SPAGeDi是一款引用率非常高的軟件，在那個還是用芯片做GWAS的時代，它還是可以勝任親緣關(guān)系計算這項工作的，但是由于測序技術(shù)的不斷提高，做GWAS的標(biāo)記量越來越高，漸漸就形成了日益增長的需求和落后的計算能力之間的矛盾。舉個例子，如果你有幾百萬標(biāo)記，幾百個樣本，這款軟件會浪費你的青春。

那我們用什么如做大標(biāo)記量的親緣關(guān)系計算呢，GCTA、LDAK、TASSEL適合你，今天我們給您講解一下如何用TASSEL去計算kinship。

提問，用TASSEL計算kinship需要幾步。

答案：4步

第一步，您要在您的服務(wù)器上安裝一個TASSEL，Windows界面版的不能夠勝任這一項工作，具體的方法可以參考我們第一期。

第二步，準(zhǔn)備好您的vcf文件，當(dāng)然我們默認(rèn)您的vcf文件是處理好的那樣。

第三步，給vcf文件排序，排成tassel認(rèn)可的序列，如果您不排序，運行任何命令都會報錯，畢竟人家軟件牛，還是有點脾氣的哈。

命令格式：run_pipeline.pl -Xmx1536m-Xms512m -SortGenotypeFilePlugin -inputFile 你的vcf文件 -outputFile 輸出vcf文件的名字 -fileType VCF

run_pipeline.pl-Xmx1536m -Xms512m -SortGenotypeFilePlugin -inputFile lecture06_genotype.vcf-outputFile lecture06cp -fileType VCF

我們ls -t看一下：

生成了個重新排序tassel可以接受的vcf文件。

第四步，開始振奮人心的親緣關(guān)系分析，同樣也是一條命令：

run_pipeline.pl-Xmx1536m -Xms512m -importGuess lecture06cp.vcf -KinshipPlugin -methodCentered_IBS -endPlugin -export tassel_kinship.txt -exportType SqrMatrix

我們再ls -t看一下：

就生成了你所需要的文件啦：

這個文件就可以直接拿去做GWAS分析了。

當(dāng)然，我們今天只講了用tassel去算kinship，如果您想用別的方法比如GCTA、LDAK、SPAGeDi去計算，沒關(guān)系，基迪奧的GWAS課程有詳細(xì)的講解，從原理到方法再到技術(shù)，你想要的這里都有，同時也會教你，如何用您算出來的親緣關(guān)系文件去繪制發(fā)表級別的圖片。

教程鏈接：http://www./class/。