2020-07-21
群體進化前面兩期(qi)帶大家實(shi)戰了基于和（點(dian)擊(ji)查看），我們今天(tian)帶(dai)大(da)家進(jin)行(xing)實戰群(qun)體遺傳結構分析(xi)，依(yi)然(ran)含分析(xi)和(he)繪圖(tu)代(dai)碼哦(e)。
所(suo)謂群(qun)(qun)體(ti)遺(yi)傳(chuan)(chuan)結構(gou)(Structure)是指遺(yi)傳(chuan)(chuan)變異在物種(zhong)或群(qun)(qun)體(ti)中的一(yi)種(zhong)非隨機分(fen)布。按照(zhao)地理(li)分(fen)布或其他劃分(fen)標準可將一(yi)個群(qun)(qun)體(ti)劃分(fen)為若干亞(ya)(ya)群(qun)(qun)，處于同(tong)一(yi)亞(ya)(ya)群(qun)(qun)內的不同(tong)個體(ti)親緣關系(xi)較高，而亞(ya)(ya)群(qun)(qun)與亞(ya)(ya)群(qun)(qun)之(zhi)間的個體(ti)則親緣關系(xi)稍遠(yuan)。群(qun)(qun)體(ti)結構(gou)分(fen)析有助于理(li)解進化過(guo)程，并且可以(yi)通過(guo)基因型和(he)表型的關聯研(yan)究確定個體(ti)所(suo)屬的亞(ya)(ya)群(qun)(qun)。
群體(ti)進化研究的(de)“三駕馬車”中通過PCA和進化樹(shu)分析可(ke)以直(zhi)觀了解(jie)群體(ti)中個體(ti)間(jian)的(de)分類關系。但如果(guo)我(wo)們要(yao)探究：
目標(biao)群體劃分為幾個亞(ya)群更合理(li)？
直接(jie)按(an)照地域劃分(fen)的群體分(fen)組是否合理(li)？
群體間(jian)是(shi)否存(cun)在基(ji)因交流？
以及每個(ge)個(ge)體混血(xue)程(cheng)度(du)是多少(shao)？
這(zhe)(zhe)些疑問PCA圖(tu)和(he)進化樹樹形圖(tu)都無法告(gao)訴我(wo)們答案(an)，這(zhe)(zhe)時候(hou)就該高大上的Structure分(fen)析的遺傳(chuan)結構圖(tu)出(chu)場啦！
圖1 遺傳(chuan)結構分析結果圖
雖然(ran)圖(tu)形很高大上(shang)但看到(dao)花花綠綠的(de)條形柱圖(tu)你是(shi)不(bu)是(shi)也心存疑惑(huo)，這(zhe)張圖(tu)代表什么生(sheng)物學意(yi)義呢？怎么解(jie)答前面的(de)提出來(lai)疑問(wen)呢？接下來(lai)我們先(xian)幫大家一(yi)一(yi)解(jie)惑(huo)。
在進行Structure分(fen)析時我們只是獲得(de)了樣本的基因型，并不(bu)知(zhi)道這(zhe)個(ge)群體(ti)(ti)實際包含幾個(ge)亞群。把群體(ti)(ti)的亞群數(shu)(shu)稱為K值，可(ke)以先預設群體(ti)(ti)亞群數(shu)(shu)等于1~n，即K=1~n，然后(hou)模擬(ni)在K=x的情況下，通過貝葉斯算法推算群體(ti)(ti)是如何分(fen)群以及(ji)每個(ge)個(ge)體(ti)(ti)的祖(zu)先來源。最(zui)后(hou)再根據每次模擬(ni)的最(zui)大似(si)然值，找出劃分(fen)這(zhe)個(ge)群體(ti)(ti)的最(zui)佳K值。
具體解讀以下面圖(tu)2為例做詳細的(de)說明。
圖(tu)2 K=3時(shi)群體(ti)遺傳結構(gou)分布(bu)圖(tu)
圖2就是(shi)假設群體結(jie)構數(shu)的(de)先驗值K=3的(de)模擬結(jie)果(guo)。每(mei)個(ge)(ge)直(zhi)方柱(zhu)子(zi)(zi)代表(biao)1個(ge)(ge)樣(yang)本(ben)(ben)（這里popA包含10個(ge)(ge)樣(yang)本(ben)(ben)，樣(yang)本(ben)(ben)名a1~a10，popB包含17個(ge)(ge)樣(yang)本(ben)(ben)，樣(yang)本(ben)(ben)名b1~b17，popC包含20個(ge)(ge)樣(yang)本(ben)(ben)，樣(yang)本(ben)(ben)名c1~c20），柱(zhu)子(zi)(zi)的(de)顏(yan)(yan)色(se)以(yi)及(ji)顏(yan)(yan)色(se)比例(li)代表(biao)這個(ge)(ge)樣(yang)本(ben)(ben)屬于(yu)哪個(ge)(ge)亞群以(yi)及(ji)祖(zu)先來(lai)源(yuan)(yuan)構成(cheng)比例(li)是(shi)多(duo)少，可(ke)以(yi)很直(zhi)觀(guan)的(de)通過顏(yan)(yan)色(se)就可(ke)以(yi)判斷這47個(ge)(ge)個(ge)(ge)體是(shi)如(ru)何劃分為3個(ge)(ge)亞群體（橙(cheng)、藍、粉(fen)）。同(tong)時(shi)也(ye)可(ke)以(yi)看出(chu)有些樣(yang)本(ben)(ben)的(de)祖(zu)先來(lai)源(yuan)(yuan)單一（僅有1種(zhong)顏(yan)(yan)色(se)），有些樣(yang)本(ben)(ben)卻(que)由2個(ge)(ge)顏(yan)(yan)色(se)組(zu)成(cheng)。例(li)如(ru)樣(yang)本(ben)(ben)c1，它對應的(de)直(zhi)方柱(zhu)子(zi)(zi)由兩(liang)(liang)種(zhong)顏(yan)(yan)色(se)構成(cheng)，大(da)概是(shi)37%的(de)藍色(se)和63%的(de)粉(fen)色(se)。說明這個(ge)(ge)個(ge)(ge)體很有可(ke)能(neng)是(shi)從(cong)兩(liang)(liang)個(ge)(ge)祖(zu)先亞群雜交而來(lai)，且兩(liang)(liang)個(ge)(ge)祖(zu)先亞群各占了37%和63%的(de)血統。
但(dan)是每一(yi)個(ge)K值(zhi)都對應著一(yi)個(ge)結果圖，哪一(yi)個(ge)K值(zhi)才是最準確(que)的呢(ni)？這里就要(yao)引出(chu)下(xia)面圖3的這張圖啦。
圖(tu)3 LK分布(bu)圖(tu)
每個K值(zhi)(zhi)基于貝葉斯模型的(de)計算方法模擬的(de)結果，都會產生對(dui)應的(de)最大似然(ran)值(zhi)(zhi)（likelihood），LK值(zhi)(zhi)是(shi)取了自然(ran)對(dui)數(shu)后輸出(chu)的(de)（ln likelihood）。ln likelihood越大，說明K值(zhi)(zhi)越接近于真實情況，但一(yi)般隨著(zhu)K值(zhi)(zhi)升高(gao)，ln likelihood值(zhi)(zhi)也會進入平臺期。最優K值(zhi)(zhi)就是(shi)進入平臺期的(de)那個拐點K值(zhi)(zhi)。
看到這里，大(da)家是不(bu)是都(dou)想知道怎(zen)么才(cai)能實現(xian)群(qun)體(ti)遺傳(chuan)結構分析呢(ni)？接下來我(wo)們就帶(dai)大(da)家一起實戰(zhan)遺傳(chuan)結構分析。實戰(zhan)示例見(jian)以下鏈接：
鏈(lian)接：//pan.baidu.com/s/1o1gNmyrQE2AA94pa-xC-8w
提取碼：56bm
軟(ruan)(ruan)件(jian)(jian)Structure、FastStructure、Admixture都可以(yi)進(jin)(jin)行遺(yi)傳(chuan)結構分析(xi)(xi)，Structure的運(yun)行速度較慢，Admixture使(shi)用(yong)的模型(xing)與Structure相同，但(dan)其憑借高速的運(yun)算速度逐漸成為群(qun)體遺(yi)傳(chuan)結構分析(xi)(xi)的主流(liu)軟(ruan)(ruan)件(jian)(jian)，接(jie)下來就以(yi)軟(ruan)(ruan)件(jian)(jian)Admixture的使(shi)用(yong)方法為例(li)帶大家(jia)逐步進(jin)(jin)行分析(xi)(xi)。
VCF文件(jian)格式轉換
利用軟(ruan)件vcftools將群(qun)體(ti)(ti)SNP的vcf文件轉(zhuan)換成Plink可以分析的格(ge)式，具體(ti)(ti)轉(zhuan)換命令如下：
usage: vcftools --vcf test.vcf --plink --out test
##--vcf 輸入vcf文件
##--plink 將數據轉換為plink格式
##--out 輸出文件前綴
得(de)到的結果文件test.map和test.bed。
SNP過濾
在(zai)Structure分析過(guo)程中需要滿足兩點，其一，各個(ge)亞群(qun)內(nei)部(bu)個(ge)體應該(gai)符合哈代-溫伯格平衡；其二，每個(ge)SNP位點是獨立的(de)(de)不(bu)連鎖的(de)(de)。一般基于重測(ce)序得到(dao)所有SNP作為Structure分析的(de)(de)輸入文件是不(bu)對的(de)(de)，需要利用Plink進行SNP的(de)(de)過(guo)濾，生成.bed文件，命令如下：
plink --noweb --file test --maf 0.05 --hwe 0.0001 --make-bed --out filter
##--noweb 不(bu)連接網絡
##--file ; 指定輸入文件
##--maf 過濾掉(diao)低于(yu)次等位基因頻(pin)率閾值的位點
##--hwe 過濾(lv)掉低(di)于(yu)Hardy-Weinberg平衡檢驗(yan)p值低(di)于(yu)閾值的(de)位點
##--make-bed 數據轉換為二(er)進(jin)制格式(shi)
##--out ; 輸(shu)出文件前綴
Admixture計算和K值的確定
K是假設的群體亞群或者(zhe)祖(zu)先數，為確定(ding)最理(li)想的K值(zhi)，可以(yi)設定(ding)K=1~10，用軟(ruan)件Admixture進行計算，命(ming)令如下：
for i in {1..10};do admixture --cv filter.bed $i|tee out${i};done
從結果(guo)文(wen)件中提取出(chu)Cross-validation error值(zhi)(zhi)即不同K值(zhi)(zhi)的(de)錯誤率，一般認為(wei)CV error最(zui)(zui)小值(zhi)(zhi)對應的(de)K值(zhi)(zhi)為(wei)最(zui)(zui)佳K值(zhi)(zhi)。提取CV error命令如下(xia)：
grep -h CV out|awk -F ':' '{print NR"\t"$2}'|sed '1i\\K\tCV_error' > CV_for_plot.txt
然(ran)后將結果在(zai)R中可視化，最終得(de)到CV error的分布圖圖4。具體代碼(ma)如下：
library(ggplot2)
mydata<-read.table("CV_for_plot.txt",header=T,sep="\t")
qplot(x=K,y=CV_error,color=I('black'),data=mydata)+geom_line(color="red",lwd=1)
圖4 CV error的(de)分(fen)布(bu)圖
根據結果顯示K=3對應的CV error最小，故為最佳K值。
遺傳結構圖形繪制
根據K=3時軟(ruan)件(jian)(jian)Admixture計(ji)算(suan)的(de)結果(guo)文件(jian)(jian)filter.3.Q利(li)用R軟(ruan)件(jian)(jian)繪(hui)制遺傳結構(gou)分布(bu)圖。
具體畫圖代碼如下(xia)：
mydata<-read.table("filter.3.Q")
barplot(t(as.matrix(mydata)),col=c("#E69F00","#56B4E9","#CC79A7"),border=NA)
然后(hou)就得(de)到我們(men)的目標(biao)結果啦！
圖(tu)5 K=3時(shi)遺傳結構(gou)分布(bu)圖(tu)
以(yi)上就是(shi)基(ji)于群(qun)體SNP進行(xing)遺傳結構分析和圖(tu)形(xing)繪制的方法，是(shi)不是(shi)沒有想象的復雜(za)呢，趕快用(yong)自己的數(shu)據試試吧(ba)！沒有服務器資(zi)源的老師也可(ke)以(yi)用(yong)Structure軟件(jian)(jian)的windows的Java圖(tu)形(xing)界(jie)面(mian)版本試試，該(gai)軟件(jian)(jian)包含三個(ge)主要參數(shu)length of burn-in period、Number of MCMC Reps after burn-in、admixture model，前(qian)兩個(ge)參數(shu)一般設(she)置為10000，第三個(ge)參數(shu)是(shi)該(gai)軟件(jian)(jian)涉及(ji)的兩種模(mo)型(xing) no admixture model 和admixture model，前(qian)者假設(she)亞群(qun)間(jian)不存在(zai)雜(za)交，后者假設(she)亞群(qun)間(jian)存在(zai)雜(za)交，一般情況下選擇(ze)admixture model更合理哦~
不管基于哪(na)種軟件(jian)如果大家實戰分析(xi)的過程中有(you)任何疑問，歡迎(ying)和(he)我們一起探討，可在文末留言或者郵件(jian)交(jiao)流（genome_support@doudin.cn）。
下面(mian)是(shi)動植物產品線的產品類型，如果對其(qi)他(ta)產品比較感興趣的，也可以郵箱反饋我們，后(hou)續我們一(yi)一(yi)為(wei)大(da)家安排。
