FASTAファイルをGREPするスクリプト

bioinformatics

FASTAファイルをアノテーションや配列でgrepする事は度々あるのですが，シェルのgrepだと合致した行だけが出てきて，そのアノテーションを含む配列が何であったか，あるいは，その配列を含む配列の注釈は何かを調べるには一手間必要です．以下のスクリプトは，FASTAを検索して，該当の（行ではなく）エントリを返すスクリプトになります．

  • 最後のエントリが表示されないバグ，修正しました．@hinaichigoさん，ありがとうございます！！！
  • 使い方

引数: FASTAフォーマットファイル
オプション(特に指定しなければ，すべての電車がおくれます):

-a: 検索範囲をannotationだけに絞ります
-s: 検索範囲を配列だけに絞ります．

例：

% ruby grep_fasta.rb -a ATNEK5 TAIR9_cds_20090619   
>AT3G20860.1 | Symbols: ATNEK5 | ATNEK5 (NIMA-RELATED KINASE5); ATP binding / kinase/ protein kinase/ protein serine/threonine kinase/ protein tyrosine kinase | chr3:7306147-7308434 FORWARD
ATGGCAAACAAGATAAGTGAAACTGCTTCCAAGATGGACGATTATGAGGTCGTGGAACAG
ATCGGACGAGGCGCTTTTGGTTCTGCTTTTCTTGTAATTCACAAGTCCGAGAGACGAAAG
TATGTGGTTAAGAAGATTCGTTTGGCTAAACAAACAGAGAGATGCAAGCTTGCTGCAATC
#!/usr/bin/env ruby
require 'optparse'

is_annot_search = true
is_seq_search = true

opt = OptionParser.new
opt.on("-a"){|v| is_seq_search = false}
opt.on("-s"){|v| is_annot_search = true}
opt.parse!(ARGV)


exp = ARGV.shift
file = ARGV.shift
fp = open(file, "r")
name = ""
seq = ""
cat_seq = ""
prev_name = ""

while line = fp.gets
  if line =~ />(.+)$/ || fp.eof
    if fp.eof
      seq += line
      cat_seq += line.chomp
    end
    name = $1
    if (is_annot_search == true && prev_name =~ /#{exp}/) || 
        (is_seq_search == true && cat_seq =~ /#{exp}/)
      puts ">#{prev_name}"
      puts seq
    end
    prev_name = name
    seq = ""
    cat_seq = ""
  else
    seq += line
    cat_seq += line.chomp
  end
end
fp.close()