安装软件 admixtrue
文件 snp
##不能有multiallel
bcftools view -M2 -v snps SNP.filtered.vcf > zy-bi.vcf
##挑选每个基因上的一个位点，需要不连锁
perl randSnps.pl < zy-bi.vcf > zy_un.vcf
## 染色体名称修改为全数字的，不能有字母
awk '{gsub(/gene/,""); print}' zy_un.vcf > zy1.vcf
awk '{gsub(/.fasta/,""); print}' zy1.vcf> zy_name.vcf
##将vcf转换为ped文件
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
####上面是思琪姐的方法，下面是admixture——manual方法
SNP.filtered.vcf
bcftools view -M2 -v snps SNP.filtered.vcf > aa.vcf
awk '{gsub(/gene/,""); print}' aa.vcf > bb.vcf
awk '{gsub(/.fasta/,""); print}' bb.vcf> cc.vcf
plink2 --vcf cc.vcf --make-bed --out fjs --allow-extra-chr --threads 12
plink2 --bfile fjs --allow-extra-chr --indep-pairwise 50 10 0.1
plink2 --bfile fjs --extract plink2.prune.in --make-bed --allow-extra-chr --out prunedData
for K in {2..6};
do admixture --cv prunedData.bed $K |tee log${K}.out;
done
grep -h CV log*.out
## R作图
setwd("E:/yangling02/SNP")
tbl=read.table("ZY.2.Q")
barplot(t(as.matrix(tbl)), col=rainbow(2),xlab="Individual #", ylab="Ancestry", border=NA)
tbl=read.table("ZY.12.Q")
barplot(t(as.matrix(tbl)), col=rainbow(12),xlab="Individual #", ylab="Ancestry", border=NA)
tbl=read.table("ZY.9.Q")
barplot(t(as.matrix(tbl)), col=rainbow(9),xlab="Individual #", ylab="Ancestry", border=NA)
