이번 스크립트는 중간에 gene list 파일을 만들 필요 없는 스크립트.
-중복이 되는 유전자가 있으면 제거하고 중복되지 않게 유전자 서열 저장
import glob,sys
from Bio import SeqIO
blast_result_file = sys.argv[1]
db_fasta_file = sys.argv[2]
output_file = sys.argv[3]
lists = []
for x in open(blast_result_file):
lists.append((x.split('\t')[1]).strip())
if len(lists) > 1:
lists.sort()
last = lists[-1]
for i in range(len(lists)-2,-1,-1):
if last == lists[i]:
del lists[i]
else:
last=lists[i]
db = SeqIO.parse(open(db_fasta_file), format='fasta')
ow = open(output+'.txt','w')
for rec in db:
name = rec.description
seq = rec.seq.tostring()
for x in lists:
if name.strip() == x.strip() or (name.strip()).find(x.strip()) != -1:
if name[:1] == '>':
ow.write(name+'\n'+seq+'\n')
else:
ow.write('>'+name+'\n'+seq+'\n')
ow.close()
사용법은 지난번 스크립트와 동일
ex) this.py blast_output_file db_file output_file
blast_output_file의 형식은 blast돌릴때 -m 8 형식입니다.
db_file은 fasta 형식이면 됩니다. 물론 blast_output_file에 gene 이름이 있어야
db_file에서도 찾아주겠죠..
여기서 나온 blast_output_file에서 두번째 컬럼의 gene 이름으로 db_file에서
서열을 추출하는 방식입니다.
Posted by gwlee
