群體結(jié)構(gòu)是指材料的亞群分化情況，會(huì)導(dǎo)致標(biāo)記間的非連鎖關(guān)聯(lián)，進(jìn)而導(dǎo)致關(guān)聯(lián)分析結(jié)果出現(xiàn)假陽(yáng)性。
地理隔離、人工選擇、移民和遺傳漂變等都可能導(dǎo)致群體分化。
是指遺傳變異在物種或群體中的一種非隨機(jī)分布；
將各材料歸到每個(gè)亞群，計(jì)算每個(gè)材料基因組變異源于第K個(gè)亞群的可能性，用Q值表示，Q值越大，表明改材料來(lái)自這個(gè)亞群的可能性越大，一般可以用來(lái)推斷祖先群，個(gè)體血緣組成，還有雜交事件;
常用軟件：Admixture、Structure、Frappe等。
隨著技術(shù)的發(fā)展，Structure速度較慢，無(wú)法滿足大量分子標(biāo)記計(jì)算的需求，因此，admixture逐漸成為群體結(jié)構(gòu)分析的主流軟件。本文將介紹如何通過(guò)admixture進(jìn)行群體結(jié)構(gòu)計(jì)算。
1.下載及安裝
1.1 下載地址
http://dalexander.github.io/admixture/index.html
1.2 安裝
$ tar xvf admixture_linux-1.3.0.tar.gz
$ cd your/path/admixture_linux-1.3.0
# 調(diào)用：./admixture
# 查看幫助：./admixture --help
2. 群體結(jié)構(gòu)計(jì)算
2.1 整理成admixture所需的.ped(12recode)格式
在plink中將vcf文件轉(zhuǎn)換成admixture所需的.ped或.bed格式:
$ cd your/path/plink1.9
$ ./plink --vcf genotype.vcf --allow-extra-chr --recode12 --out genotype12 --autosome-num 27
--vcf 輸入文件名
--allow-extra-chr 允許其他格式染色體，如scaffold
--recode12 二進(jìn)制編碼
--out 輸出文件名
--autosome-num 設(shè)置染色體數(shù)目，默認(rèn)人類染色體數(shù)
2.2 Admixture
$ cd your/path/admixture_linux-1.3.0
# 創(chuàng)建任務(wù)文件
$ vim adm.sh
# vim 文件名
# i 輸入 左下角出現(xiàn)insert，可以輸入
for K in 2 3 4 5 6 7 8 9 10; do ./admixture --cv root12.ped $K | tee log${K}.out; done
# ESC鍵 insert消失
# 退出
$ :wq
# 提交任務(wù)
$ bsub -n 4 -o log sh adm.sh
#查看任務(wù)
$ bjobs
JOBID USER STAT QUEUE FROM_HOST EXEC_HOST JOB_NAME SUBMIT_TIME
913421 xxx RUN normal login 4*compute11 sh adm.sh Aug 24 01:14
每個(gè)K值都會(huì)生成兩個(gè)文件，.P和.Q
P：儲(chǔ)存推斷的祖先種群的等位基因頻率
Q：每個(gè)樣本中各個(gè)祖先種群所占的百分比。
3. 最佳分群數(shù)確定及可視化
3.1 確定最佳分群數(shù)
查看cv值，cv error最小的K值為最佳分群數(shù)。
$ grep -h CV log*.out
CV error (K=10): 0.65873
CV error (K=2): 0.71095
CV error (K=3): 0.63424
CV error (K=4): 0.68598
CV error (K=5): 0.67584
CV error (K=6): 0.66818
CV error (K=7): 0.66301
CV error (K=8): 0.66083
CV error (K=9): 0.65919
3.2 群體結(jié)構(gòu)可視化
將CV結(jié)果復(fù)制粘貼至Excel中，繪制折線圖。圖中可看出最佳分群數(shù)為K=3。
在R中繪制群體結(jié)構(gòu)圖
提供幾個(gè)我喜歡的配色：
K=3 "#FF4500","#9ACD32","#6495ED"
K=4 "#336666","darkred","steelblue","#CC9933"
K=5 "#FF4500","#5F7A61","#6495ED","#986D8E","#F6D167"
將K=3時(shí)的.Q文件拷貝至Windows中
> setwd("D:/數(shù)據(jù)/GWAS/群體結(jié)構(gòu)")
> library("ggplot2")
> install.packages(c("ggplot2","gridExtra","label.switching","tidyr","remotes"),repos="https://cloud.r-project.org")
> remotes::install_github('royfrancis/pophelper')
> library("pophelper")
> tbl=read.table("genotype.3.Q")
> pdf("admixture.pdf",width = 9,height = 3)
> colorpal =c("#FF4500","#9ACD32","#6495ED")
> cols=rep(colorpal,700)
> barplot(t(as.matrix(tbl)), col=cols, xlab="", ylab="Ancestry",border = NA)
> dev.off()
3.3 確定樣本屬于哪個(gè)亞群
當(dāng)確定最佳分群數(shù)是3時(shí)，打開(kāi)K=3時(shí)的.Q文件，文件共包含三列，每行為一個(gè)樣本，三列中哪一個(gè)數(shù)值最大，則這個(gè)樣本屬于哪一個(gè)亞群。
引用轉(zhuǎn)載請(qǐng)注明出處，如有錯(cuò)誤敬請(qǐng)指出。
