例えばP972-blast.txtというファイルにBlast検索の結果がはいっているとき、
ruby blast_to_fasta_imp.rb P972-blast.txt > P972.fasta
のようにするとP972.fastaにfasta形式のデータが記録される。以下にblast_to_fasta_imp.rbの内容を解説する。
01: lines = ARGF.readlines.collect{|a| a.chomp}
02: seqs = {}
03: num_seqs = Hash.new(0)
04: seq = ''
05: name = nil
06: query = ''
07: q_name = 'Query '
08: seqs[q_name] = ''
09: index = [' ',' ']+('2'..'9').to_a + ('A'..'Z').to_a
10: lines.each do |line|
11: if /\A>/ =~ line
12: if name
13: seqs[name] = seq
14: seqs[q_name] = query if seqs[q_name].size < query.size
15: end
16: sections = line.split(/\| */)
17: if sections.size == 2
18: words = sections[1].split(" ")[1..-1]
19: else
20: words = sections[2].split(" ")
21: end
22: if words[1] == 'x'
23: base_name = words[0][0,4]+'_'+words[2][0,4]
24: else
25: base_name = words[0][0,4]+'_'+words[1][0,4]
26: end
27: num_seqs[base_name] += 1
28: name = base_name + index[num_seqs[base_name]] + ' '*(9-base_name.length)
29: seq = ''
30: query = ''
31: elsif /Sbjct/ =~ line
32: seq << line.split(/ +/)[2]+"\n"
33: elsif /Query / =~ line
34: query << line.split(/ +/)[2]+"\n"
35: end
36: end
37: seqs[name] = seq
38: seqs[q_name] = query if seqs[q_name].size < query.size
39: seqs.each do |name,seq|
40: print ">#{name}\n"
41: print seq
42: end
1-8行 入力ファイル（ARGF、P972-blast.txtを想定）の内容を改行を取り除いた後に１行ずつの配列としてlinesに格納する。seqsが空のハッシュ、num_seqsがデフォルト値が0の空のハッシュ、seqが空の文字列、nameがnilに設定される。また、Query配列も書き出すためにハッシュのキーとしてq_name、塩基配列としてqueryを用意する。q_nameの代わりに直接'Query 'という文字列を使ってもよいが、こうしておくと'Query 'という文字列を別のものに変えたい場合（例えば、コマンドラインから名前を指定する）、１か所の変更で済む。
1行目の表現は授業では扱わなかったが、rubyで多用される、イテレータと呼ばれる機能で、collectメソッドは配列、ハッシュなどの各要素に対して{}内の処理を行った結果を新たな配列として返す。例えば
[1, 2, 3].collect{|x| x*2}
は、[2, 4, 6]を返す。ほかに、条件に当てはまる要素だけを抜き出すselectメソッド、条件に当てはまるものだけを取り除くrejectメソッドなどがある。例えば、
[18, 0, 14, 1, 19, 12].select{|x| x % 2 ==0}
は[18, 0, 14, 12]を返す。%は整数の割り算の余りである。
num_seqsは塩基配列の定義から作られる名前が重複した時に番号で区別するための準備。また、9行では、そのために使うindex配列を[' ', ' ', '2', '3', ... , '9', 'A', 'B', ... , 'Z']として用意している。
10行 linesの内容を１行ずつlineという変数に入れて、対応するend (36行) までの処理を繰り返す。
11行の正規表現は\Aが行の先頭を示すので、先頭に>がある行が一致する。Blast検索の結果ではこのような行は塩基配列の比較部分で検索結果の配列の定義を示している。この行から次の定義の行までが一つの塩基配列である。これに当てはまらない場合は31行の条件によって、"Sbjct"という文字を含んでいれば、32行の処理を行い、そうでなければさらに33行の条件によって、"Query "という文字を含んでいれば34行の処理を行う。Queryの後ろのスペース一つがポイントで、これがないとファイルの先頭近くの"Query="という行も当てはまり、処理でエラーが出る。これらの条件に当てはまらなければ何もせずに次の行がlineとなって11行からの処理を繰り返す。
定義の行であったときはまず12行でnameに何か（定義の行から取り出した配列の名前）が入っているかどうかを確かめ、入っていればseqsハッシュにその名前をキーとしてseqの塩基配列を格納する。この名前と塩基配列は現在の定義の行の直前までのデータであることが分かるだろう。また、ファイルの最初の定義の行でここを実行している時はnameがnilになっているため、この処理は行わない。ハッシュに直前までのSbjctの配列を格納したら（13行）、14行で今までにseqsに入っているQuery配列の長さよりも今回のQuery配列の方が長ければそれをseqsに保存する。14行のように行の後ろに条件を書く形式を後置if文といい、条件が成立した時の処理が１つだけで、成立しない時の処理がない場合に使え、通常のif文よりも簡潔に書ける。
16-28行 nameを新しい名前にする。まず16行のline.split(/\| */)で"|"（パイプ文字）とその後に続く任意の数の空白で定義行を切り分けた配列が得られる。例えば、次の行は
>gb|DQ226511.1| Morus indica chloroplast, complete genome
[">gb", "DQ226511.1", "Morus indica chloroplast, complete genome"]となる。この３番目の文字列（[2]）をさらにスペースで切り分ける(split(' '))と["Morus", "indica", "chloroplast,", "complete", "genome"]となり、これをwordsとする。
17-21行で２種類のファイルフォーマットを判断し、それぞれ適切な部分を取り出している。18行の[1..-1]は２番目の要素から最後の要素までを意味する。Rubyの負の添字は後ろからという意味。
22行の条件は雑種かどうかの判断で、たとえばソメイヨシノならwordsは["Prunus", "x", "yedoensis", ...]、オランダイチゴは["Fragaria", "x", "ananasa", ...]となっていて、先頭から２つの文字列のそれぞれ先頭４文字をアンダースコアでつなぐと"Prun_x"、"Frag_x"となってしまう。できれば"Prun_yedo" や"Frag_anan"としたいので、この条件を入れた。２番目の単語が"x"でないときは普通に先頭の２語をつないでいる。この名前をbase_nameとする。
27-28行 さて、こうして作った名前だが、Blast検索結果に同じ種でも違う配列がヒットしているばあい、複数の配列が同じ名前になってしまう。そこで、num_seqsハッシュを用意しbase_nameをキーとして、この名前の出現回数を入れておく。27行のnum_seqs[base_name] += 1はnum_seqs[base_name] = num_seqs[base_name]+1と同じ。この回数が1回ならばindex[1]は空白、2回以上なら数字、10回以上ならアルファベットが順に使われる。このように、28行でnameを構築し、さらに、base_nameが9文字より短い場合は空白で埋めてnameを10文字にする。別の解決法としてアクセション番号を名前に含めることが考えられるが、Phylipの使用が前提の場合、名前が10文字に限られることとのトレードオフとなる。
37-38行 改良前のプログラムでは最後の配列が出力されていなかったため、13, 14行と同じ処理を行って最後の配列をハッシュに格納する。もしこの処理がもう少し複雑な場合、メソッドを定義してそれぞれの場所から呼ぶようにする。これらの処理に変更を加える場合に１箇所だけ変更すれば良いからである。
39行でseqsハッシュからキーと値の組をひとつずつname, seqに入れて40, 41行の処理を行っている。配列の部分は32, 34行ですでに改行記号("\n")を追加しているので、printのみで改行もされる。
2011年7月12日作成、2012年7月30日更新
國分 尚 (Hisashi Kokubun)