Blast2GO アノテーションが付けられた SIMAP データベースを使用して GO アノテーションを実行しようとしています。すべて問題ありませんが、エントリ番号が GO に関連付けられているファイルでアクセッション番号を見つけようとすると問題が発生します。問題は、実際には入力ファイルに数字があるのに、スクリプトがその数字を見つけられないことです。良い結果が得られずにいくつかのことを試しました(再一致、リストに挿入してから要素を抽出するなど)。GOがエントリ番号に関連付けられているファイルには、次の構造があります(アクセッション番号、GOターム、blats2goスコア):
1f0ba1d119f52ff28e907d2b5ea450db GO:0007154 79
1f0ba1d119f52ff28e907d2b5ea450db GO:0005605 99
Python コード:
import re
from Bio.Blast import NCBIXML
from Bio import SeqIO
input_file = open('/home/fpiston/Desktop/test_go/test2.fasta', 'rU')
result_handle = open('/home/fpiston/Desktop/test_go/test2.xml', 'rU')
save_file = open('/home/fpiston/Desktop/test_go/test2.out', 'w')
fh = open('/home/fpiston/Desktop/test_go/Os_Bd_Ta_blat2go_fake', 'rU')
q_dict = SeqIO.to_dict(SeqIO.parse(input_file, "fasta"))
blast_records = NCBIXML.parse(result_handle)
hits = []
for blast_record in blast_records:
if blast_record.alignments:
list = (blast_record.query).split()
if re.match('ENA|\w*|\w*', list[0]) != None:
list2 = list[0].split("|")
save_file.write('%s\t' % list2[1])
else:
save_file.write('%s\t' % list[0])
for alignment in blast_record.alignments:
for hsp in alignment.hsps:
h = alignment.hit_def
for l in fh:
ls = l.split() #at this point all right
if h in ls: #here, 'h' in not found in 'fh'
print h
print 'ok'
save_file.write('%s\t' % ls[1])
save_file.write('\n')
hits.append(blast_record.query.split()[0])
misses =set(q_dict.keys()) - set(hits)
for i in misses:
list = i.split("|")
if len(list) > 1:
save_file.write('%s\t' % list[1])
else:
save_file.write('%s\t' % list)
save_file.write('%s\n' % 'no_match')
save_file.close()
これはmartineau (fh.seek(0))を修正したコードです。
#!/usr/bin/env python
import sys
import re
from Bio.Blast import NCBIXML
from Bio import SeqIO
input_file = sys.argv[1] #queries sequences in fasta format
out_blast_file = sys.argv[2] #name of the blast results file
output_file = sys.argv[3] #name of the output file
result_handle = open(out_blast_file, 'rU')
fh = open('/home/fpiston/Desktop/test_go/Os_Bd_Ta_blat2go', 'rU')
q_dict = SeqIO.to_dict(SeqIO.parse(open(input_file), "fasta"))
blast_records = NCBIXML.parse(result_handle)
save_file = open(output_file, 'w')
hits = []
for blast_record in blast_records:
if blast_record.alignments:
list = (blast_record.query).split()
if re.match('ENA|\w*|\w*', list[0]) != None:
list2 = list[0].split("|")
save_file.write('\n%s\t' % list2[1])
else:
save_file.write('\n%s\t' % list[0])
for alignment in blast_record.alignments:
for hsp in alignment.hsps:
hit = alignment.hit_def
save_file.write('%s\t' % hit)
fh.seek(0)
for l in fh:
ls = l.split()
if ls[0] in hit:
save_file.write('%s\t' % ls[1])
hits.append(blast_record.query.split()[0])
misses =set(q_dict.keys()) - set(hits)
for i in misses:
list = i.split("|")
if len(list) > 1:
save_file.write('\n%s\t' % list[1])
else:
save_file.write('\n%s\t' % list)
save_file.write('%s' % 'no_match')
save_file.close()