code_ruby/lib/maasha/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 class Digest
  29   include Enumerable
  30
  31   # Initialize a digest object with the following arguments:
  32   # - seq: A sequence object.
  33   # - pattern: A restriction enzyme recognition pattern.
  34   # - cut_pos: Offset from match position where the enzyme cuts.
  35   def initialize(seq, pattern, cut_pos)
  36     @seq     = seq
  37     @pattern = disambiguate(pattern)
  38     @cut_pos = cut_pos
  39     @offset  = 0
  40   end
  41
  42   # Method to get the next digestion product from a sequence.
  43   def each
  44     @seq.seq.upcase.scan @pattern do
  45       pos = $`.length + @cut_pos - 1
  46
  47       if pos >= 0 and pos < @seq.length - 2
  48         seq = Seq.new("#{@seq.seq_name}[#{@offset}-#{pos}]", @seq.seq[@offset .. pos], @seq.type)
  49
  50         yield seq
  51       end
  52
  53       @offset = pos + 1
  54     end
  55
  56     if @offset < 0
  57       @offset = 0
  58     elsif @offset > @seq.length
  59       @offset = 0
  60     end
  61
  62     seq = Seq.new("#{@seq.seq_name}[#{@offset}-#{@seq.length - 1}]", @seq.seq[@offset .. @seq.length], @seq.type)
  63
  64     yield seq
  65
  66     self # conventionally
  67   end
  68
  69   private
  70
  71   # Method that returns a regexp object with a restriction
  72   # enzyme pattern with ambiguity codes substituted to the
  73   # appropriate regexp.
  74   def disambiguate(pattern)
  75     ambiguity = {
  76       'A' => "A",
  77       'T' => "T",
  78       'U' => "T",
  79       'C' => "C",
  80       'G' => "G",
  81       'M' => "[AC]",
  82       'R' => "[AG]",
  83       'W' => "[AT]",
  84       'S' => "[CG]",
  85       'Y' => "[CT]",
  86       'K' => "[GT]",
  87       'V' => "[ACG]",
  88       'H' => "[ACT]",
  89       'D' => "[AGT]",
  90       'B' => "[CGT]",
  91       'N' => "[GATC]"
  92     }
  93
  94     new_pattern = ""
  95
  96     pattern.upcase.each_char do |char|
  97       if ambiguity.has_key? char
  98         new_pattern << ambiguity[char]
  99       else
 100         raise DigestError, "Could not disambiguate residue: #{char}"
 101       end
 102     end
 103
 104     Regexp.new(new_pattern)
 105   end
 106 end