将多个样本的vcf文件转化为Phylip输入格式的python脚

在分析变异过滤得到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文件