使用vcf文件构建进化树

发表于 | 更新于 | 分类于

首先，过滤 vcf 文件

1
2
3
4
5
6
7
1.
vcftools --vcf final.vcf --min-alleles 2.0 --max-alleles 2.0 --max-missing 0.95 --non-ref-af 0.05 --max-non-ref-af 0.95 --recode -c >filter.vcf

java -jar SnpSift.jar filter "(countHom() > X) & (QUAL >= 30 )" filter.vcf > filtered.vcf #（可选）X可以指定数字，如10，100 意义为纯合变异的样本数
2.
vcftools --vcf filter.vcf --plink --out filter #转plink ped格式
plink --file filter --recode transpose 12 --out filtered #转plink tped格式

然后，构建距离矩阵

1
2
3
4
5
6
7
8
9
1.
VCF2Dis -InPut filtered.vcf -OutPut p_dis.mat
2.
emmax-kin filtered -h -s #-s 构建IBS 亲缘矩阵
2.1
python hIBS2PHYLIP.py filtered.hIBS.kinf
2.2
cat filtered.hIBS.kinf|awk -v OFS='\t' '{for(i=1;i<=NF;i++){$i=sqrt(($i-1)^2)}}{print $0}' >PHYLIP_need.txt
cut -f1 -d ' ' filtered.tfam |paste - PHYLIP_need.txt >p_dis.mat #最后在结果第一行加样本数

最后，使用PHYLIPNEW构建进化树

1
2
PHYLIPNEW-3.69.650/bin/fneighbor -datafile p_dis.mat  -outfile tree.out.txt -matrixtype s -treetype n -outtreefile tree.out.tre
#-treetype n 指构建Neighbor-joining树，可选u 指构建UPGMA树 -outgrno 参数可指定外群的数量

hIBS2PHYLIP.py

1
2
3
4
5
6
7
8
9
10
import sys
import numpy as np
if (len(sys.argv)==2):
    with open(sys.argv[1],'r') as IN, open("PHYLIP_need.txt",'w') as OUT:
        f_np = np.loadtxt(IN,delimiter='\t')
        result = abs((f_np-1)) #减一取绝对值
        np.savetxt(OUT,result,fmt="%.6f",delimiter="\t")
        print("Done!")
else:
    print("Usage: python hIBS2PHYLIP.py test.hIBS.kinf")