背景
群体遗传学中测出很多个个体,得到了最终的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是一个进行单倍型分析的一个软件。具体的使用我后续会单独再具体介绍。
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)
image
转载请注明出处
微信公众号:oddxix
简书作者:oddxix
网友评论