読者です 読者をやめる 読者になる 読者になる

ピエール瀧になりたい

備忘録として

blastの出力からヒットした配列を抜き出す

blastをかけて、ヒットした部分の配列をfastaで出力するコード(biopythonとかまともに使えれば絶対こんなのあるんだろうな・・・)

makeblastdbは事前にやっておく

ゲノムからあるタンパク質周辺のシンテニーを知りたい場合などは、周辺領域の長さを指定することでごっそり抜き出したり

 >|python|
#!/usr/bin/env python
#coding: UTF-8

#よく使うモジュールは何も考えずに書く

import sys
import re
import os
import glob
import copy
import commands
from bisect import bisect
from Bio import SeqIO
import urllib2
import datetime


"""
使い方
blastEX blasyn query.fasta db.fasta num(default:0)

で走らせると以下のコマンドを実行し、ヒットした配列を取得

blastn -query query -db db.fasta -outfmt 6

numは取得する周辺配列の長さ。
num = 100 にすればヒットした配列の前後100aaを取得。
前後が100aa以下の場合は取れるだけの長さを取る。

"""

now = datetime.datetime.now().strftime("%Y%m%d_%H%M%S")

argvs = sys.argv # コマンドライン引数を格納したリストの取得
argc = len(argvs) # 引数の個数

blast = argvs[1]

#num変数の設定

if argc > 4:
num = int(argvs[4])
else:
num = 0

#blastの実行

os.system("%s -query %s -db %s -outfmt 6 > BlastResult_%s"%(blast,argvs[2],argvs[3],now))
outputFasta = open("Output_%s.fasta"%now,"w")

#配列の取得

for line in open("BlastResult_%s"%now) :
#print line
if len(line) > 2:
ID = line.split("\t")[1]
seqstart = min(int(line.split("\t")[8]),int(line.split("\t")[9]))
seqend = max(int(line.split("\t")[8]),int(line.split("\t")[9]))
index = 0
handle = open(argvs[3], "rU")
for record in SeqIO.parse(handle, "fasta"):
if record.id == ID:
if seqstart < num and seqend > len(record.seq)-num:
outputFasta.write( ">"+record.id +"\n" )
outputFasta.write( str(record.seq ) +"\n" )
elif seqstart < num:
outputFasta.write( ">"+record.id +"\n" )
outputFasta.write( str(record.seq[0:seqend + num] ) +"\n" )
elif seqend > len(record.seq)-num:
outputFasta.write( ">"+record.id +"\n" )
outputFasta.write( str(record.seq[seqstart-num:] ) +"\n" )
else:
outputFasta.write( ">"+record.id +"\n" )
outputFasta.write( str(record.seq[seqstart-num:seqend + num]) +"\n" )
handle.close()
outputFasta.close() ||<