在分析變異過濾得到SNP時(shí)，一般都會用PHYLIP構(gòu)建NJ進(jìn)化樹

但是phylip又不能直接把vcf文件作為輸入文件，它的輸入格式要求如下

image.png

第一行的兩個(gè)數(shù)字分別為樣本數(shù)和SNP數(shù)目
第二行的第一列為物種名稱，一定要是10個(gè)字符的長度，長度不足可以用空格填充.第二列則是要用于比對的SNP序列。兩列之間沒有分隔符，而且后面的SNP每10個(gè)字符有一個(gè)空格.
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"]
###參數(shù)傳遞
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的位點(diǎn)
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)
#將每個(gè)vcf中的snp連接起來，如果該位點(diǎn)存在突變，則輸出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文件