Pythonでvcfファイルのコードを処理する

3099 ワード

問題の背景:
私は今何百もの同じ菌株のゲノムを持っていますが、私は差の大きいサンプルだけで後続の分析をする必要があります.私は似たような冗長サンプルを削除します.では、私はどうすればいいですか.直接megaで進化木を構築することはできない.ゲノムが大きすぎて分析できないからだ.私はこれらのサンプルの1つを参照ゲノムとして、残りのサンプルをこのサンプルに比較して、call SNP、最終的に数百のvcfファイルを得て、これらのvcfファイルはこれらのサンプルの参照ゲノムに対するSNPを含んで、1つのスクリプトを書いてこれらのSNPをそれぞれ1つのシーケンスとして接続して、すべてのシーケンスは1つのfastaファイルに出力して、megaで進化木を構築すればいいです!
以下はpythonコードです
#Pyvcf python    vcf      
import vcf
import os
import numpy
import collections
#vcf       
filepath=r'C:\Users\18019\Desktop\VCF\Vibrio.VCF'
filelist=os.listdir(filepath)
output_name=[]
for ech in filelist:
    output_name.append(ech.replace('_',''))
output_name+=['VibrioFF75']
dicREF_list=[]
dicALT_list=[]
for ech in filelist:
    ech_vcf=vcf.Reader(filename=r'C:\Users\18019\Desktop\VCF\Vibrio.VCF\%s' % ech)
    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 i in dicREF_list:
    SNP_REF=dict(SNP_REF,**i)
pos_list=list(SNP_REF.keys())
pos_list.sort()
print(len(pos_list))
#   vcf  snp    ,         ,   ALT,    REF,    SNP      
all_list=[]
for ech_dic in dicALT_list:
    ech_dic_index = dicALT_list.index(ech_dic)
    ech_name=output_name[ech_dic_index]
    ech_list=[]
    for pos in pos_list:
        if pos in ech_dic:
            theSNP = str(ech_dic[pos]).replace('[', '').replace(']', '')
            ech_list.append(theSNP)
        else:
            theSNP = str(SNP_REF[pos]).replace('[', '').replace(']', '')
            ech_list.append(theSNP)
    all_list.append(ech_list)
#   SNP           ,         
snp_array=[]
for pos in pos_list:
    theSNP=str(SNP_REF[pos]).replace('[','').replace(']','')
    snp_array.append(theSNP)

for i in all_list:
    snp_array=numpy.row_stack((snp_array,i))
snp_array=numpy.transpose(snp_array)
array_len=len(snp_array)
print(array_len)
snp_least=[]
i=0
while i'+output_name[-1]+'
' for pos in pos_list: theSNP=str(SNP_REF[pos]).replace('[','').replace(']','') allStr+=theSNP output=open(r'C:\Users\18019\Desktop\VCF\Vibrio.VCF\total_snp.txt','w') output.write(allStr) output.close()

コードは参考にして、具体的な問題は具体的に分析します.の