[BioRuby-ja] FlatFile blat parser

Tomoaki NISHIYAMA tomoakin @ kenroku.kanazawa-u.ac.jp
2007年 8月 22日 (水) 23:10:56 EDT


西山@金沢です

blatはmultifastaをqueryにしたときも、blastと異 
なり、queryが変わる度に
headerを出力ことはせず、各行にqueryとtargetの 
idがでます。
それに対応して、Bio::Blat::Reportは
FlatFileとして処理させたとき1つのentryでたくさんの 
hitという
事になっているのだと思いますが、使う都合としては、
query毎にバラしたファイルで検索するのはやりたくないし、
結果はやっぱりquery毎に欲しいということがあります。
例えば、query毎にbest hitの場所を選びたいと言う場合 
です。

blatの出力に区切りが無いとはいえ、同じqueryの結果は続けて 
でて来ますので
FlatFileとしての処理時は同じqueryの名前が続く
範囲を以て1つのentryとして返してくれると便利ということで
flatfile_splitterをつけてみました。
splitする段で各行をparseするので、
Hit.newを各行について再度呼ばなくて済むようにtextではなくて
HitのArrayを返すようにしてみました。(速度はおそらく2割程 
度の差)

また、全部一括システムで全てをHashに押し込んで同じ 
queryのhit毎に処理すると、
百数十MBytesのpslファイルを処理するのに
数GBytes以上のメモリーが必要になっていたのですが、
ずっと少ないメモリーで処理出来るようになります。

いかがでしょう。

$ diff  -u report.rb reportex.rb
--- report.rb   2006-06-29 20:54:33.000000000 +0900
+++ reportex.rb 2007-08-23 11:24:34.000000000 +0900
@@ -41,16 +41,70 @@
      # If multiple query sequences are given, these values
      # will be incorrect.
      #
-    class Report #< DB
+    class Reportex #< DB
        # Delimiter of each entry. Bio::FlatFile uses it.
        # In Bio::Blat::Report, it it nil (1 entry 1 file).
        DELIMITER = RS = nil # 1 file 1 entry
+      def self.flatfile_splitter(dbclass, stream)
+#        puts "Reportex:flatfile_splitter called"
+        BlatSplitter.new(stream)
+      end
+
+      class BlatSplitter
+        def initialize(stream)
+          @stream = stream
+        end
+        def flatfile_splitter(dbclass, stream)
+          BlatSplitter.new(stream)
+        end
+        def skip_leader
+          s = @stream.gets
+          if s =~ /psLayout/
+            while s !~ /------/
+              s = @stream.gets
+            end
+          else
+            @stream.ungets(s)
+          end
+        end
+        def get_entry
+          p0 = @entry_pos_flag ? @stream.pos : nil
+          firstline = @stream.gets
+          if firstline == nil
+            return nil
+          end
+          e = firstline
+          hits = Array.new
+          first = Hit.new(firstline)
+          hits << first
+          query = first.query.name
+          nextline=@stream.gets
+          if nextline != nil
+            nexthit = Hit.new(nextline)
+            while nexthit.query.name == query
+              e += nextline
+              hits << nexthit
+              nextline = @stream.gets
+              break if nextline == nil
+              nexthit = Hit.new(nextline)
+            end
+          end
+          @stream.ungets(nextline) if nextline != nil
+          p1 = @entry_pos_flag ? @stream.pos : nil
+          @entry_start_pos = p0
+          @entry = e
+          @entry_ended_pos = p1
+          @hits = hits
+        end
+      end
+
        # Creates a new Bio::Blat::Report object from BLAT result  
text (String).
        # You can use Bio::FlatFile to read a file.
        # Currently, results created with options -out=psl (default) or
        # -out=pslx are supported.
        def initialize(text)
+        if text.is_a?(String)
          flag = false
          head = []
          @hits = []
@@ -72,6 +126,9 @@
            end
          end
          @columns = parse_header(head)
+        elsif text.is_a?(Array)
+          @hits=text
+        end
        end
        # hits of the result.



-- 
Tomoaki NISHIYAMA

Advanced Science Research Center,
Kanazawa University,
13-1 Takara-machi,
Kanazawa, 920-0934, Japan




BioRuby-ja メーリングリストの案内