VCF格式转换
利用vcftools将群体SNP的vcf文件转换成plink可以分析的格式
vcftools --vcf snp.vcf --plink --out snp
## 参数
--vcf 输入vcf文件
--plink 将数据转成plink格式
--out 输出文件后缀
##输出的结果： .map 和.ped 文件
SNP过滤
plink --noweb --file snp --maf 0.05 --hwe 0.0001 --make-bed --out filter.snp
##参数
--noweb 不联网
--file 输入文件（上一步的map和ped文件）
--maf 0.05 maf<0.05,则过滤
--hwe 0.0001 不符合哈迪温伯格则过滤
--make-bed 转成bed格式
--out 输出文件前缀
##输出结果：.bed .bim .fam文件
Admixture计算
##创建bash脚本
vi admixture.sh
##写入以下脚本
for i in {1..10};
do admixture --cv filter.snp.bed $i >> log.txt
done
##增加可执行权限
chmod u+x admixture.sh
##运行脚本
./admixture.sh
##参数
--cv .bed文件 K值
>>log.txt 从1到10运行的out结果都整合到log.txt中，以便grep提取cv行
提取CV error
grep CV log.txt | awk -F ':' '{print NR"\t"$2}' | sed '1i\K\tCV_error' >> CV_for_plot.txt
## 参数
awk -F ':' 以：分隔
NR 当前行号
"\t" 空格
$2 第二列
sed '1i\K\tCV_error' 表示给第一行添加K 空格 CV_error
sed '5i\内容' 表示给第五行添加内容
在R中进行结果可视化，确定最佳K值
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)
##CV error最小的为最佳K值
画图(假设最佳K=3，filter.snp.3.Q结果文件作为画群体结构图的输入源文件)
tbl=read.table("filter.snp.3.Q")
barplot(t(as.matrix(tbl)), col=rainbow(3),xlab="Individual #", ylab="Ancestry", border=NA)
