背景
群體遺傳學中測出很多個個體纺裁,得到了最終的SNP vcf文件挥下,需要將其分成群體汰具,看哪幾個物種聚在一起。一般使用的軟件就是STRUCTURE粥烁,但是STREUTURE運行速度極慢贤笆,后面frappe軟件提升了速度,但是也不是很快页徐;admixture憑借其運算速度,成為了主流的分析軟件银萍。
admixture
1. admixture的輸入文件為:PLINK (.bed), ordinary PLINK (.ped), or EIGENSTRAT (.geno)
所以需要使用vcftools進行格式轉換
vcftools --vcf my.vcf --plink --out plink
輸出的結果為:plink.ped與plin.map文件
2.過濾SNP文件
使用Plink軟件变勇。Plink是一個開放的,免費的全基因組關聯分析工具贴唇。其分析的基礎是基因型和表型數據搀绣。通過整合gplink 和Haploview 使得結果變得可視化。 gplink 是基于Java的圖形用戶界面戳气,許多plink 命令可以通過其實現链患。而且易于與Haploview進行整合。Haploview是一個進行單倍型分析的一個軟件瓶您。具體的使用我后續(xù)會單獨再具體介紹麻捻。
plink --noweb --file plink --geno 0.05 --maf 0.05 --hwe 0.0001 --make-bed --out QC
輸出的結果為:QC.bed(binary file,genotype information)纲仍、QC.fam(first six columns of plink.bed)、QC.bim(extended MAP file:two extra cols=allele names)
得到bed文件后就可以使用admixture進行分析了
3.群體結構分析
如果不知道理想的K值可以設定1-5進行分析
for K in 1 2 3 4 5; do admixture --cv hapmap3.bed $K | tee log${K}.out; done
4.提取出CV值
grep -h CV log*.out
CV error (K=1): 0.55248
CV error (K=2): 0.48190
CV error (K=3): 0.47835
CV error (K=4): 0.48236
CV error (K=5): 0.48985
5.使用R畫圖
tbl=read.table("hapmap3.3.Q")
barplot(t(as.matrix(tbl)), col=rainbow(3),xlab="Individual #", ylab="Ancestry", border=NA)
轉載請注明出處
微信公眾號:oddxix
簡書作者:oddxix