Admixture 进行遗传结构分析

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)

美化