added shred_seq biopiece

[biopieces.git] / bp_bin / shred_seq

diff --git a/bp_bin/shred_seq b/bp_bin/shred_seq
new file mode 100755 (executable)
index 0000000..a59b688
--- /dev/null
+++ b/bp_bin/shred_seq
@@ -0,0 +1,94 @@
+#!/usr/bin/env ruby
+
+# Copyright (C) 2007-2011 Martin A. Hansen.
+
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU General Public License
+# as published by the Free Software Foundation; either version 2
+# of the License, or (at your option) any later version.
+
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+
+# http://www.gnu.org/copyleft/gpl.html
+
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+# This program is part of the Biopieces framework (www.biopieces.org).
+
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+# Shred sequences in the stream into subsequences.
+
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+
+require 'biopieces'
+require 'seq'
+require 'pp'
+
+class Seq
+  # Method that shreds the sequence of a sequence object into
+  # random subsequences of a specified size. Subsequences are
+  # alternatly picked from the + or - strand (i.e. reverse complemented)
+  # until the specified coverage is reached.
+  #
+  # seq.shred { |subseq| } -> Seq
+  # seq.shred              -> []
+  def shred(size, max_cov)
+    entries = []
+    sum     = 0
+    strand  = '+'
+
+    while (sum / self.seq.length) < max_cov
+      start    = rand(self.seq.size - size)
+      stop     = start + size - 1
+      seq_name = self.seq_name + "[#{start + 1}-#{stop + 1}:#{strand}]"
+      seq      = self.seq[start .. stop]
+      entry    = Seq.new(seq_name, seq, 'dna')
+
+      entry.revcomp if strand == '-'
+
+      if block_given?
+        yield entry
+      else
+        entries << entry
+      end
+      
+      strand = (strand == '+') ? '-' : '+'
+
+      sum += size
+    end
+
+    entries
+  end
+end
+
+casts = []
+casts << {:long=>'size',     :short=>'s', :type=>'uint', :mandatory=>true, :default=>500, :allowed=>nil, :disallowed=>'0'}
+casts << {:long=>'coverage', :short=>'c', :type=>'uint', :mandatory=>true, :default=>100, :allowed=>nil, :disallowed=>'0'}
+
+bp = Biopieces.new
+
+options = bp.parse(ARGV, casts)
+
+bp.each_record do |record|
+  if record.has_key? :SEQ
+    entry = Seq.new(record[:SEQ_NAME], record[:SEQ], record[:SCORES])
+
+    entry.shred(options[:size], options[:coverage]) do |subentry|
+      bp.puts subentry.to_bp
+    end
+  end
+end
+
+# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+
+__END__