#删除多等位基因的位点
bcftools view -M2 -v snps - exclude SNP.filtered.vcf > BSZ_ZY_bi.vcf
###plink过滤
#1 转换格式
plink --vcf file.vcf --recode --out file
#2 位点过滤
--mind 0.10 去除基因型丢失率大于10%的个体样本：
--geno 0.05 大于95%的个体都具有的变异位点才保留，其他去除：
--maf 0.01 次等位基因频率，频率较低的第二等位基因的频率（防止假阳性）。
--hwe 0.0001 保留符合Hardy—Weinbery 的变异位点
### LD过滤生成 vcf 格式（服务器上有LD.sh脚本）
plink2 --vcf zy_bsz-bi.vcf --mind 0.10 --maf 0.05 --geno 0.05 --hwe 0.0001 \
--out zy_bsz_fi --recode vcf-iid \
--allow-extra-chr --set-missing-var-ids @:# \
--keep-allele-order
plink2 --vcf zy_bsz_fi.vcf \
--indep-pairwise 50 10 0.2 \
--out tmp.ld --allow-extra-chr \
--set-missing-var-ids @:#
plink2 --vcf zy_bsz_fi.vcf \
--extract tmp.ld.prune.in \
--out all.LDfilter --recode vcf-iid \
--keep-allele-order --allow-extra-chr \
--set-missing-var-ids @:#
#3 stuctrue作图
bcftools view -M2 -v snps snp.vcf > zy-bi.vcf
perl randSnps.pl < zy-bi.vcf > zy_un.vcf
awk '{gsub(/gene/,""); print}' zy_un.vcf > zy1.vcf
awk '{gsub(/.a/,""); print}' zy1.vcf> zy_name.vcf
plink2 --vcf zy_name.vcf --make-bed --out ZY --allow-extra-chr --threads 12
for K in {2..12};
do admixture --cv ZY.bed $K |tee log${K}.out;
done
grep -h CV log*.out
#CV error (K=2): 0.40129
#V error (K=3): 0.48130
#CV error (K=4): 0.47680
#CV error (K=5): 0.58563
#CV error (K=6): 0.60299
#CV error (K=7): 0.67999
#CV error (K=8): 0.69322
#CV error (K=9): 0.53627
