例えばP972-blast.txtというファイルにBlast検索の結果がはいっているとき、
ruby blast_to_fasta.rb P972-blast.txt > P972.fasta
のようにするとP972.fastaにfasta形式のデータが記録される。以下にblast_to_fasta.rbの内容を解説する。
01: lines = ARGF.readlines
02: seqs = {}
03: seq = ""
04: name = nil
05:
ここまででP972-blast.txtの内容が１行ずつの配列としてlinesに格納され、seqsが空のハッシュ、seqが空の文字列、nameがnilに設定される。次の行でこのlinesの内容を１行ずつlineという変数に入れて、対応するend (22行) までの処理を繰り返す。
06: lines.each do |line|
07: line = line.chomp
08: if /\A>/ =~ line
09: if name
10: seqs[name] = seq
11: end
12: sections = line.split(/\| */)
13: if sections.size == 2
14: name = sections[1].split(" ")[1,2].join("_")
15: else
16: name = sections[2].split(" ")[0,2].join("_")
17: end
18: seq = ""
19: elsif /Sbjct/ =~ line
20: seq << line.split(/ +/)[2]+"\n"
21: end
22: end
23:
8行の正規表現は\Aが行の先頭を示すので、先頭に>がある行が一致する。Blast検索の結果ではこのような行は塩基配列の比較部分で検索結果の配列の定義を示している。この行から次の定義の行までが一つの塩基配列である。これに当てはまらない場合は19行の条件によって、"Sbjct"という文字を含んでいれば、20行の処理を行い、そうでなければ何もせずに次の行がlineとなって7行からの処理を繰り返す。
定義の行であったときはまず9行でnameに何か（定義の行から取り出した配列の名前）が入っているかどうかを確かめ、入っていればseqsハッシュにその名前をキーとしてseqの配列を格納する。この名前と配列は現在の定義の行の直前までのデータであることが分かるだろう。また、ファイルの最初の定義の行でここを実行している時はnameがnilになっているため、この処理は行わない。ハッシュに直前までの配列を格納したら、12〜17行でnameを新しい名前にし、18行でseqを空にする。
12行では"|"（パイプ文字）とその後に続く任意の数の空白で定義行を切り分けた配列をsectionsに入れる。例えば、NCBIの場合、次の行は
>gb|DQ226511.1| Morus indica chloroplast, complete genome
[">gb", "DQ226511.1", "Morus indica chloroplast, complete genome"]となる。DDBJとNCBIではこの行が異なるので、13行で判断し、NCBIのときは３番目の文字列（sections[2]）をさらにスペースで切り分ける(split(' '))と["Morus", "indica", "chloroplast,", "complete", "genome"]となり、先頭から２つの単語を取り出して([0,2])アンダースコア"_"でつないでいる（join('_')）。DDBJのときは２番目の文字列（sections[1]）の２番目から２つの単語が名前になっている。名前にアクセション番号を含める場合、この部分を変更すればよい。
24: seqs.each do |name,seq|
25: print ">#{name}\n"
26: print seq
27: end
24行でseqsハッシュからキーと値の組をひとつずつname, seqに入れて25,26行の処理を行っている。配列の部分は20行ですでに改行記号("\n")を追加しているので、printのみで改行もされる。
このプログラムには致命的な欠陥があり、最後の配列が出力されない。また、blast検索で問い合わせに使った配列 (Queryで表される) も出力されない。このあたりが改良すべき点である。
2011年5月30日作成、2012年6月18日更新
國分 尚 (Hisashi Kokubun)