code_ruby/lib/maasha/seq/digest.rb

   1 # Copyright (C) 2007-2011 Martin A. Hansen.
   2
   3 # This program is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU General Public License
   5 # as published by the Free Software Foundation; either version 2
   6 # of the License, or (at your option) any later version.
   7
   8 # This program is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 # GNU General Public License for more details.
  12
  13 # You should have received a copy of the GNU General Public License
  14 # along with this program; if not, write to the Free Software
  15 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  16
  17 # http://www.gnu.org/copyleft/gpl.html
  18
  19 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  20
  21 # This software is part of the Biopieces framework (www.biopieces.org).
  22
  23 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  24
  25 # Error class for all exceptions to do with Digest.
  26 class DigestError < StandardError; end
  27
  28 module Digest
  29   # Method to get the next digestion product from a sequence.
  30   def each_digest(pattern, cut_pos)
  31     pattern = disambiguate(pattern)
  32     results = []
  33     offset  = 0
  34
  35     self.seq.upcase.scan pattern do
  36       pos = $`.length + cut_pos
  37
  38       if pos >= 0 and pos < self.length - 2
  39         subseq = self[offset .. pos - 1]
  40         subseq.seq_name = "#{self.seq_name}[#{offset}-#{pos - offset - 1}]"
  41
  42         block_given? ? (yield subseq) : (results << subseq)
  43       end
  44
  45       offset = pos
  46     end
  47
  48     if offset < 0 or offset > self.length
  49       offset = 0
  50     end
  51
  52     subseq = self[offset .. -1]
  53     subseq.seq_name = "#{self.seq_name}[#{offset}-#{self.length - 1}]"
  54
  55     block_given? ? (yield subseq) : (results << subseq)
  56
  57     return results unless block_given?
  58   end
  59
  60   private
  61
  62   # Method that returns a regexp object with a restriction
  63   # enzyme pattern with ambiguity codes substituted to the
  64   # appropriate regexp.
  65   def disambiguate(pattern)
  66     ambiguity = {
  67       'A' => "A",
  68       'T' => "T",
  69       'U' => "T",
  70       'C' => "C",
  71       'G' => "G",
  72       'M' => "[AC]",
  73       'R' => "[AG]",
  74       'W' => "[AT]",
  75       'S' => "[CG]",
  76       'Y' => "[CT]",
  77       'K' => "[GT]",
  78       'V' => "[ACG]",
  79       'H' => "[ACT]",
  80       'D' => "[AGT]",
  81       'B' => "[CGT]",
  82       'N' => "[GATC]"
  83     }
  84
  85     new_pattern = ""
  86
  87     pattern.upcase.each_char do |char|
  88       if ambiguity[char]
  89         new_pattern << ambiguity[char]
  90       else
  91         raise DigestError, "Could not disambiguate residue: #{char}"
  92       end
  93     end
  94
  95     Regexp.new(new_pattern)
  96   end
  97 end