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

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

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

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

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

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

##按照LD過(guò)濾并轉(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時(shí)候的admixture

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

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

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

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

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

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

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

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

提問(wèn)，用TASSEL計(jì)算kinship需要幾步。

答案：4步

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

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

第三步，給vcf文件排序，排成tassel認(rèn)可的序列，如果您不排序，運(yùn)行任何命令都會(huì)報(bào)錯(cuò)，畢竟人家軟件牛，還是有點(diǎ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看一下：

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

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

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

我們?cè)賚s -t看一下：

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

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

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

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