fixed seq qual length check

[biopieces.git] / bp_bin / digest_seq

diff --git a/bp_bin/digest_seq b/bp_bin/digest_seq
index 3fab531b4249cd144a5092cb58eb15c4d2b29d87..2aa9ff9837fded6ad832d44fffacf99b08f96742 100755 (executable)
--- a/bp_bin/digest_seq
+++ b/bp_bin/digest_seq
@@ -24,39 +24,44 @@
 
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-# Splits sequences in the stream at a given restriction enzyme's cleavage points.
+# Split sequences in the stream at a given restriction enzyme's cleavage points.
 
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-
-require 'biopieces'
-require 'fasta'
-require 'seq'
+require 'maasha/biopieces'
+require 'maasha/seq'
 
 casts = []
 casts << {:long=>'pattern', :short=>'p', :type=>'string', :mandatory=>true, :default=>nil, :allowed=>nil, :disallowed=>nil}
 casts << {:long=>'cut_pos', :short=>'c', :type=>'int',    :mandatory=>true, :default=>nil, :allowed=>nil, :disallowed=>nil}
 
-bp = Biopieces.new
+options = Biopieces.options_parse(ARGV, casts)
+
+Biopieces.open(options[:stream_in], options[:stream_out]) do |input, output|
+  input.each_record do |record|
+    if record[:SEQ_NAME] and record[:SEQ]
+      seq = Seq.new_bp(record)
 
-options = bp.parse(ARGV, casts)
+      seq.each_digest(options[:pattern].to_s, options[:cut_pos]) do |digest|
+        new_record = digest.to_bp
 
-bp.each_record do |record|
-  bp.puts record
+        if new_record[:SEQ_NAME] =~ /\[(\d+)-(\d+)\]$/
+          s_beg = $1
+          s_end = $2
+        end
 
-  if record.has_key? :SEQ_NAME and record.has_key? :SEQ
-    seq    = Seq.new(record[:SEQ_NAME], record[:SEQ])
-    digest = Digest.new(seq, options[:pattern].to_s, options[:cut_pos])
+        new_record[:S_BEG]    = s_beg 
+        new_record[:S_END]    = s_end
+        new_record[:REC_TYPE] = "DIGEST"
 
-    digest.each do |subseq|
-      new_record            = subseq.to_bp
-      new_record[:REC_TYPE] = "DIGEST"
-      bp.puts new_record
+        output.puts new_record
+      end
+    else
+      output.puts record
     end
   end
 end
 
-
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<