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