Structrue分析

安装软件 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)