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'
  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 module BackTrack
  34   OK_PATTERN = Regexp.new('^[flsycwphqrimtnkvadegu]+$')
  35
  36   # ------------------------------------------------------------------------------
  37   #   str.scan(pattern[, max_mismatches [, max_insertions [,max_deletions]]])
  38   #   -> Array
  39   #   str.scan(pattern[, max_mismatches [, max_insertions [,max_deletions]]]) { |match|
  40   #     block
  41   #   }
  42   #   -> Match
  43   #
  44   # ------------------------------------------------------------------------------
  45   # Method to iterate through a sequence to locate pattern matches starting at a
  46   # given offset allowing for a maximum number of mismatches, insertions, and
  47   # deletions. Matches found in block context return the Match object. Otherwise
  48   # matches are returned in an Array of Match objects.
  49   def patscan(pattern, offset = 0, max_mismatches = 0, max_insertions = 0, max_deletions = 0)
  50     raise BackTrackError, "Bad pattern: #{pattern}" unless pattern.downcase =~ OK_PATTERN
  51     matches = []
  52
  53     if block_given?
  54       while result = self.track(pattern, offset, max_mismatches, max_insertions, max_deletions)
  55         yield Match.new(result.first, result.last, self.seq[result.first ... result.first + result.last])
  56         offset = result.first + 1
  57       end
  58     else
  59       while result = self.track(pattern, offset, max_mismatches, max_insertions, max_deletions)
  60         matches << Match.new(result.first, result.last, self.seq[result.first ... result.first + result.last])
  61         offset = result.first + 1
  62       end
  63     end
  64
  65     return matches.empty? ? nil : matches unless block_given?
  66   end
  67
  68   private
  69
  70   inline do |builder|
  71     builder.add_static "id_seq", 'rb_intern("@seq")', "ID"
  72
  73     builder.prefix %{
  74       #define MATCH(A,B) ((bitmap[(unsigned short) A] & bitmap[(unsigned short) B]) != 0)
  75     }
  76
  77     builder.prefix %{
  78       unsigned short bitmap[256] = {
  79           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  80           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  81           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  82           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  83           0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
  84           0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
  85           0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
  86           0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
  87           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  88           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  89           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  90           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  91           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  92           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  93           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  94           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
  95       };
  96     }
  97
  98     # ss is the start of the string, used only for reporting the match endpoints.
  99     builder.prefix %{
 100       unsigned int backtrack(char* ss, char* s, char* p, int mm, int ins, int del)
 101       {
 102           unsigned int r = 0;
 103
 104           while (*s && MATCH(*s, *p)) ++s, ++p;    // OK to always match longest segment
 105
 106           if (!*p)
 107               return (unsigned int) (s - ss);
 108           else
 109           {
 110               if (mm && *s && *p && (r = backtrack(ss, s + 1, p + 1, mm - 1, ins, del))) return r;
 111               if (ins && *s &&      (r = backtrack(ss, s + 1, p, mm, ins - 1, del)))     return r;
 112               if (del && *p &&      (r = backtrack(ss, s, p + 1, mm, ins, del - 1)))     return r;
 113           }
 114
 115           return 0;
 116       }
 117     }
 118
 119     # Find all occurrences of p starting at pos in s, with at most
 120     # mm mismatches, ins insertions and del deletions.
 121     builder.c %{
 122       static VALUE track(char* p, unsigned int pos, int mm, int ins, int del)
 123       {
 124           VALUE        seq = rb_ivar_get(self, id_seq);
 125           char*        s   = StringValuePtr(seq);
 126           char*        ss  = s;
 127           unsigned int e;
 128           VALUE        tuple;
 129
 130           if (pos < strlen(s))
 131           {
 132               s += pos;
 133
 134               while (*s)
 135               {
 136                   e = backtrack(ss, s, p, mm, ins, del);
 137
 138                   if (e)
 139                   {
 140                       tuple = rb_ary_new();
 141                       rb_ary_push(tuple, INT2FIX((int) (s - ss)));
 142                       rb_ary_push(tuple, INT2FIX((int) e - (s - ss)));
 143                       return tuple;
 144                   }
 145
 146                   s++;
 147               }
 148           }
 149
 150           return Qnil;
 151        }
 152     }
 153   end
 154
 155   class Match
 156     attr_reader :pos, :length, :match
 157
 158     def initialize(pos, length, match)
 159       @pos    = pos
 160       @length = length
 161       @match  = match
 162     end
 163
 164     def to_s
 165       "#{pos}:#{length}:#{match}"
 166     end
 167   end
 168 end
 169
 170
 171 __END__