Treemix检测杂交事件

treemix由Joseph K. Pickrell和Jonathan K. Pritchard开发，文章Inference of population splits and mixtures from genome-wide allele frequency data. 通过从多个种群中获得等位基因频率，返回该种群的最大似然树，并推断可能发生的杂交事件。基于等位基因频率构建的进化树不能代表物种树。

1. 准备输入文件

（1）过滤缺失率和LD

#过滤缺失率
FILE=test
vcftools --gzvcf $FILE.vcf.gz --max-missing 1 --recode --stdout | gzip > $FILE.noN.vcf.gz
#下载网上过滤LD的脚本
wget https://github.com/joanam/scripts/raw/master/ldPruning.sh
chmod +x ldPruning.sh
./ldPruning.sh $FILE.noN.vcf.gz
gzip $FILE.noN.LDpruned.vcf

（2）转换格式

文件格式转换脚本集：conversion script
plink2treemix.py下载 here

输入文件通过vcf和一个clust文件产生。clust文件提供有关哪个样品属于哪个分类群的信息。 包含三列，其中第一列和第二列指示样本名称，第三列指示分类单元名称。

#生成clust文件
bcftools query -l $FILE.vcf.gz | awk '{split($1,pop,"."); print $1"\t"$1"\t"pop[2]}' > test.clust
#使用文件格式转换脚本集中的代码生成treemix输入文件
vcf2treemix.sh $FILE.LDpruned.vcf.gz $FILE.clust

2. 运行Treemix

for i in {0..5}
do
 treemix -i $FILE.treemix.frq.gz -m $i -o $FILE.$i -root GoldenJackal -bootstrap -k 500 -noss > treemix_${i}_log &
done

3. 结果可视化

R环境中

prefix="test.LDpruned"
library(RColorBrewer)
library(R.utils)
#加载treemix提供的R脚本，下载位置[`treemix`](https://speciationgenomics.github.io/usr/local/apps/treemix/1.12/bin/plotting_funcs.R)
source("plotting_funcs.R")
#绘制6个混合树
par(mfrow=c(2,3))
for(edge in 0:5){
  plot_tree(cex=0.8,paste0(prefix,".",edge))
  title(paste(edge,"edges"))
}
#绘制残差图
for(edge in 0:5){
 plot_resid(stem=paste0(prefix,".",edge),pop_order="dogs.list")
}

参考：
Treemix说明书
文件格式大全