2023年12月7日发(作者：)
将多个样本的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=ntParser()
_argument("-file_path",required=True,help="the vcf file directory which contain the vcf files need to handed")
_argument("-out_type",required=True,help="The output snp file type,it should be fasta or phy")
_argument("-out_path",required=True,help="The output file directory")
args=_args()
type_list=["fasta","phy"]
###参数传递
if _path:
file_index=("-file_path")
file_path=[file_index+1]
else:
else:
raise Exception("Please input the file path!")
if _type:
out_type_index=("-out_type")
out_type=[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 _path:
out_path_index=("-out_path")
out_path=[out_path_index+1]
else:
raise Exception("Please input the out put directory")
###保存文件名
filelist=r(file_path)
file_num=len(filelist)
outname_list=[]
for file in filelist:
outname_((".")[0])
dicREF_list=[]
dicALT_list=[]
for file in filelist:
ech_vcf=(filename=r'%s/%s' % (file_path,file))
dicREF,dicALT={},{}
for SNP in ech_vcf:
if _snp == 1:
dicREF[ + '_' + str()] =
if len() > 1:
dicALT[ + '_' + str()] = [0]
else:
dicALT[ + '_' + str()] =
dicREF_(dicREF)
dicALT_(dicALT)
#求出所有vcf文件中snp的位点
SNP_REF={}
for dic in dicREF_list:
SNP_REF=dict(SNP_REF,**dic)
pos_list=list(SNP_())
pos_()
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_(ech_str)
if out_type=="fasta":
with open(r"%s/" % 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/" % out_path,"w") as f:
print("%dt%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]
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文件
发布评论