code_ruby/lib/maasha/backtrack.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 require 'inline' # RubyInline
  26
  27 # Error class for all exceptions to do with BackTrack.
  28 class BackTrackError < StandardError; end
  29
  30 # Module containing code to locate nucleotide patterns in sequences allowing for
  31 # ambiguity codes and a given maximum mismatches, insertions, and deletions. The
  32 # pattern match engine is based on a backtrack algorithm.
  33 # Insertions are nucleotides found in the pattern but not in the sequence.
  34 # Deletions are nucleotides found in the sequence but not in the pattern.
  35 # Algorithm based on code kindly provided by j_random_hacker @ Stackoverflow.
  36 module BackTrack
  37   OK_PATTERN = Regexp.new('^[flsycwphqrimtnkvadegu]+$')
  38   MAX_MIS    = 5 # Maximum number of mismatches allowed
  39   MAX_INS    = 5 # Maximum number of insertions allowed
  40   MAX_DEL    = 5 # Maximum number of deletions allowed
  41
  42   # ------------------------------------------------------------------------------
  43   #   str.scan(pattern[, max_mismatches [, max_insertions [,max_deletions]]])
  44   #   -> Array
  45   #   str.scan(pattern[, max_mismatches [, max_insertions [,max_deletions]]]) { |match|
  46   #     block
  47   #   }
  48   #   -> Match
  49   #
  50   # ------------------------------------------------------------------------------
  51   # Method to iterate through a sequence to locate pattern matches starting at a
  52   # given offset allowing for a maximum number of mismatches, insertions, and
  53   # deletions. Insertions are nucleotides found in the pattern but not in the
  54   # sequence. Deletions are nucleotides found in the sequence but not in the
  55   # pattern. Matches found in block context return the Match object. Otherwise
  56   # matches are returned in an Array of Match objects.
  57   def patscan(pattern, offset = 0, max_mismatches = 0, max_insertions = 0, max_deletions = 0)
  58     raise BackTrackError, "Bad pattern: #{pattern}" unless pattern.downcase =~ OK_PATTERN
  59     raise BackTrackError, "offset: #{offset} out of range (0 ... #{self.length - 1})"        unless (0 ... self.length).include? offset
  60     raise BackTrackError, "max_mismatches: #{max_mismatches} out of range (0 .. #{MAX_MIS})" unless (0 .. MAX_MIS).include? max_mismatches
  61     raise BackTrackError, "max_insertions: #{max_insertions} out of range (0 .. #{MAX_INS})" unless (0 .. MAX_INS).include? max_insertions
  62     raise BackTrackError, "max_deletions: #{max_deletions} out of range (0 .. #{MAX_DEL})"   unless (0 .. MAX_DEL).include? max_deletions
  63
  64     matches = []
  65
  66     if block_given?
  67       while result = self.track(pattern, offset, max_mismatches, max_insertions, max_deletions)
  68         yield Match.new(result.first, result.last, self.seq[result.first ... result.first + result.last])
  69         offset = result.first + 1
  70       end
  71     else
  72       while result = self.track(pattern, offset, max_mismatches, max_insertions, max_deletions)
  73         matches << Match.new(result.first, result.last, self.seq[result.first ... result.first + result.last])
  74         offset = result.first + 1
  75       end
  76     end
  77
  78     return matches.empty? ? nil : matches unless block_given?
  79   end
  80
  81   private
  82
  83   inline do |builder|
  84     builder.add_static "id_seq", 'rb_intern("@seq")', "ID"
  85
  86     # Macro for matching nucleotides including ambiguity codes.
  87     builder.prefix %{
  88       #define MATCH(A,B) ((bitmap[(unsigned short) A] & bitmap[(unsigned short) B]) != 0)
  89     }
  90
  91     # Bitmap for matching nucleotides including ambiguity codes.
  92     # For each value bits are set from the left: bit pos 1 for A,
  93     # bit pos 2 for T, bit pos 3 for C, and bit pos 4 for G.
  94     builder.prefix %{
  95       unsigned short bitmap[256] = {
  96           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  97           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  98           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  99           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 100           0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
 101           0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
 102           0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
 103           0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
 104           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 105           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 106           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 107           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 108           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 109           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 110           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 111           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
 112       };
 113     }
 114
 115     # Backtrack algorithm for matching a pattern (p) starting in a sequence (s) allowing for mm
 116     # mismatches, ins insertions and del deletions. ss is the start of the sequence, used only for
 117     # reporting the match endpoints.
 118     builder.prefix %{
 119       unsigned int backtrack(char* ss, char* s, char* p, int mm, int ins, int del)
 120       {
 121           unsigned int r = 0;
 122
 123           while (*s && MATCH(*s, *p)) ++s, ++p;    // OK to always match longest segment
 124
 125           if (!*p)
 126               return (unsigned int) (s - ss);
 127           else
 128           {
 129               if (mm && *s && *p && (r = backtrack(ss, s + 1, p + 1, mm - 1, ins, del))) return r;
 130               if (ins && *p &&      (r = backtrack(ss, s, p + 1, mm, ins - 1, del)))     return r;
 131               if (del && *s &&      (r = backtrack(ss, s + 1, p, mm, ins, del - 1)))     return r;
 132           }
 133
 134           return 0;
 135       }
 136     }
 137
 138     # Find pattern (p) in a sequence (@seq) starting at pos, with at most mm mismatches, ins
 139     # insertions and del deletions.
 140     builder.c %{
 141       static VALUE track(char* p, unsigned int pos, int mm, int ins, int del)
 142       {
 143           VALUE        seq = rb_ivar_get(self, id_seq);
 144           char*        s   = StringValuePtr(seq);
 145           char*        ss  = s;
 146           unsigned int e;
 147           VALUE        tuple;
 148
 149           if (pos < strlen(s))
 150           {
 151               s += pos;
 152
 153               while (*s)
 154               {
 155                   e = backtrack(ss, s, p, mm, ins, del);
 156
 157                   if (e)
 158                   {
 159                       tuple = rb_ary_new();
 160                       rb_ary_push(tuple, INT2FIX((int) (s - ss)));
 161                       rb_ary_push(tuple, INT2FIX((int) e - (s - ss)));
 162                       return tuple;
 163                   }
 164
 165                   s++;
 166               }
 167           }
 168
 169           return Qnil;
 170        }
 171     }
 172   end
 173
 174   # Class containing match information.
 175   class Match
 176     attr_reader :pos, :length, :match
 177
 178     def initialize(pos, length, match)
 179       @pos    = pos
 180       @length = length
 181       @match  = match
 182     end
 183
 184     def to_s
 185       "#{pos}:#{length}:#{match}"
 186     end
 187   end
 188 end
 189
 190
 191 __END__