在分析变异过滤得到SNP时,一般都会用PHYLIP构建NJ进化树
但是phylip又不能直接把vcf文件作为输入文件,它的输入格式要求如下
第一行的两个数字分别为样本数和SNP数目
第二行的第一列为物种名称,一定要是10个字符的长度,长度不足可以用空格填充.第二列则是要用于比对的SNP序列。两列之间没有分隔符,而且后面的SNP每10个字符有一个空格.
python脚本如下:
###author qiujunhui
import os
import sys
import argparse
import vcf
pars=argparse.ArgumentParser()
pars.add_argument("-file_path",required=True,help="the vcf file directory which contain the vcf files need to handed")
pars.add_argument("-out_type",required=True,help="The output snp file type,it should be fasta or phy")
pars.add_argument("-out_path",required=True,help="The output file directory")
args=pars.parse_args()
type_list=["fasta","phy"]
###参数传递
if args.file_path:
file_index=sys.argv.index("-file_path")
file_path=sys.argv[file_index+1]
else:
raise Exception("Please input the file path!")
if args.out_type:
out_type_index=sys.argv.index("-out_type")
out_type=sys.argv[out_type_index+1]
if out_type not in type_list:
raise Exception("Please input the right out put type:fasta or phy")
else:
raise Exception("Please input the out put type")
if args.out_path:
out_path_index=sys.argv.index("-out_path")
out_path=sys.argv[out_path_index+1]
else:
raise Exception("Please input the out put directory")
###保存文件名
filelist=os.listdir(file_path)
file_num=len(filelist)
outname_list=[]
for file in filelist:
outname_list.append(file.split(".")[0])
dicREF_list=[]
dicALT_list=[]
for file in filelist:
ech_vcf=vcf.Reader(filename=r'%s/%s' % (file_path,file))
dicREF,dicALT={},{}
for SNP in ech_vcf:
if SNP.is_snp == 1:
dicREF[SNP.CHROM + '_' + str(SNP.POS)] = SNP.REF
if len(SNP.ALT) > 1:
dicALT[SNP.CHROM + '_' + str(SNP.POS)] = SNP.ALT[0]
else:
dicALT[SNP.CHROM + '_' + str(SNP.POS)] = SNP.ALT
dicREF_list.append(dicREF)
dicALT_list.append(dicALT)
#求出所有vcf文件中snp的位点
SNP_REF={}
for dic in dicREF_list:
SNP_REF=dict(SNP_REF,**dic)
pos_list=list(SNP_REF.keys())
pos_list.sort()
snp_num=len(pos_list)
#将每个vcf中的snp连接起来,如果该位点存在突变,则输出ALT,否则输出REF,保证每条SNP序列长度相同
all_list=[]
for ech_dic in dicALT_list:
ech_str=''
for pos in pos_list:
if pos in ech_dic:
theSNP = str(ech_dic[pos]).replace('[', '').replace(']', '')
ech_str+=theSNP
else:
theSNP = str(SNP_REF[pos]).replace('[', '').replace(']', '')
ech_str+=theSNP
all_list.append(ech_str)
if out_type=="fasta":
with open(r"%s/output.fasta" % out_path,"w") as f:
for i in range(0,len(outname_list)):
out_name=">"+outname_list[i]
print(out_name,file=f)
print(all_list[i],file=f)
if out_type=="phy":
with open(r"%s/output.phy" % out_path,"w") as f:
print("%d\t%d" %(file_num,snp_num),file=f)
for i in range(0,len(outname_list)):
out_str=outname_list[i]
if(len(out_str)<=10):
out_str=out_str+" "*(10-len(out_str))
else:
out_str=out_str[:10]
for index in range(10,len(all_list[i]),10):
out_str+=all_list[i][index-10:index]+" "
out_str+=all_list[i][index:]
print("%s" %(out_str),file=f)
可以选择输出fasta文件或者phy文件