Pythonでvcfファイルのコードを処理する
3099 ワード
問題の背景:
私は今何百もの同じ菌株のゲノムを持っていますが、私は差の大きいサンプルだけで後続の分析をする必要があります.私は似たような冗長サンプルを削除します.では、私はどうすればいいですか.直接megaで進化木を構築することはできない.ゲノムが大きすぎて分析できないからだ.私はこれらのサンプルの1つを参照ゲノムとして、残りのサンプルをこのサンプルに比較して、call SNP、最終的に数百のvcfファイルを得て、これらのvcfファイルはこれらのサンプルの参照ゲノムに対するSNPを含んで、1つのスクリプトを書いてこれらのSNPをそれぞれ1つのシーケンスとして接続して、すべてのシーケンスは1つのfastaファイルに出力して、megaで進化木を構築すればいいです!
以下はpythonコードです
コードは参考にして、具体的な問題は具体的に分析します.の
私は今何百もの同じ菌株のゲノムを持っていますが、私は差の大きいサンプルだけで後続の分析をする必要があります.私は似たような冗長サンプルを削除します.では、私はどうすればいいですか.直接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()
コードは参考にして、具体的な問題は具体的に分析します.の