code_ruby/lib/maasha/seq/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 # http://stackoverflow.com/questions/7557017/approximate-string-matching-using-backtracking/
  37 module BackTrack
  38   OK_PATTERN = Regexp.new('^[bflsycwphqrimtnkvadegu]+$')
  39   MAX_MIS    = 5 # Maximum number of mismatches allowed
  40   MAX_INS    = 5 # Maximum number of insertions allowed
  41   MAX_DEL    = 5 # Maximum number of deletions allowed
  42
  43   def patmatch(pattern, offset = 0, max_mismatches = 0, max_insertions = 0, max_deletions = 0)
  44     self.patscan(pattern, offset, max_mismatches, max_insertions, max_deletions) do |m|
  45       if block_given?
  46         yield m
  47       else
  48         return m
  49       end
  50
  51       break
  52     end
  53
  54     if block_given?
  55       yield nil
  56     else
  57       return nil
  58     end
  59   end
  60
  61   # ------------------------------------------------------------------------------
  62   #   str.scan(pattern[, max_mismatches[, max_insertions[, max_deletions]]])
  63   #   -> Array
  64   #   str.scan(pattern[, max_mismatches[, max_insertions[, max_deletions]]]) { |match|
  65   #     block
  66   #   }
  67   #   -> Match
  68   #
  69   # ------------------------------------------------------------------------------
  70   # Method to iterate through a sequence to locate pattern matches starting at a
  71   # given offset allowing for a maximum number of mismatches, insertions, and
  72   # deletions. Insertions are nucleotides found in the pattern but not in the
  73   # sequence. Deletions are nucleotides found in the sequence but not in the
  74   # pattern. Matches found in block context return the Match object. Otherwise
  75   # matches are returned in an Array of Match objects.
  76   def patscan(pattern, offset = 0, max_mismatches = 0, max_insertions = 0, max_deletions = 0)
  77     raise BackTrackError, "Bad pattern: #{pattern}"                                          unless pattern.downcase =~ OK_PATTERN
  78     raise BackTrackError, "offset: #{offset} out of range (0 ... #{self.length})"            unless (0 ... self.length).include? offset
  79     raise BackTrackError, "max_mismatches: #{max_mismatches} out of range (0 .. #{MAX_MIS})" unless (0 .. MAX_MIS).include? max_mismatches
  80     raise BackTrackError, "max_insertions: #{max_insertions} out of range (0 .. #{MAX_INS})" unless (0 .. MAX_INS).include? max_insertions
  81     raise BackTrackError, "max_deletions: #{max_deletions} out of range (0 .. #{MAX_DEL})"   unless (0 .. MAX_DEL).include? max_deletions
  82
  83     matches = []
  84
  85     if block_given?
  86       while result = scan_C(self.seq, self.length, pattern, offset, max_mismatches, max_insertions, max_deletions)
  87         yield Match.new(result.first, result.last, self.seq[result.first ... result.first + result.last])
  88         offset = result.first + 1
  89       end
  90     else
  91       while result = scan_C(self.seq, self.length, pattern, offset, max_mismatches, max_insertions, max_deletions)
  92         matches << Match.new(result.first, result.last, self.seq[result.first ... result.first + result.last])
  93         offset = result.first + 1
  94       end
  95     end
  96
  97     return matches.empty? ? nil : matches unless block_given?
  98   end
  99
 100   private
 101
 102   inline do |builder|
 103     # Macro for matching nucleotides including ambiguity codes.
 104     builder.prefix %{
 105       #define MATCH(A,B) ((bitmap[A] & bitmap[B]) != 0)
 106     }
 107
 108     # Bitmap for matching nucleotides including ambiguity codes.
 109     # For each value bits are set from the left: bit pos 1 for A,
 110     # bit pos 2 for T, bit pos 3 for C, and bit pos 4 for G.
 111     builder.prefix %{
 112       char bitmap[256] = {
 113           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 114           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 115           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 116           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 117           0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
 118           0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
 119           0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
 120           0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
 121           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 122           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 123           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 124           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 125           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 126           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 127           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 128           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
 129       };
 130     }
 131
 132     # Backtrack algorithm for matching a pattern (p) starting in a sequence (s) allowing for mis
 133     # mismatches, ins insertions and del deletions. ss is the start of the sequence, used only for
 134     # reporting the match endpoints. State is used to avoid ins followed by del and visa versa which
 135     # are nonsense.
 136     builder.prefix %{
 137       unsigned int backtrack(
 138         char         *ss,     // Sequence start
 139         char         *s,      // Sequence
 140         char         *p,      // Pattern
 141         unsigned int  mis,    // Max mismatches
 142         unsigned int  ins,    // Max insertions
 143         unsigned int  del,    // Max deletions
 144         int           state   // Last event: mis, ins or del
 145       )
 146       {
 147         unsigned int r = 0;
 148
 149         while (*s && MATCH(*s, *p)) ++s, ++p;    // OK to always match longest segment
 150
 151         if (!*p)
 152           return (unsigned int) (s - ss);
 153         else
 154         {
 155           if (mis && *s && *p &&            (r = backtrack(ss, s + 1, p + 1, mis - 1, ins, del, 0))) return r;
 156           if (ins && *p && (state != -1) && (r = backtrack(ss, s, p + 1, mis, ins - 1, del, 1)))     return r;
 157           if (del && *s && (state != 1)  && (r = backtrack(ss, s + 1, p, mis, ins, del - 1, -1)))    return r;
 158         }
 159
 160         return 0;
 161       }
 162     }
 163
 164     # Find pattern (p) in a sequence (s) starting at pos, with at most mis mismatches, ins
 165     # insertions and del deletions.
 166     builder.c %{
 167       VALUE scan_C(
 168         VALUE _s,     // Sequence
 169         VALUE _len,   // Sequence length
 170         VALUE _p,     // Pattern
 171         VALUE _pos,   // Start position
 172         VALUE _mis,   // Maximum mismatches
 173         VALUE _ins,   // Maximum insertions
 174         VALUE _del    // Maximum deletions
 175       )
 176       {
 177         char         *s   = StringValuePtr(_s);
 178         char         *p   = StringValuePtr(_p);
 179         unsigned int  len = FIX2UINT(_len);
 180         unsigned int  pos = FIX2UINT(_pos);
 181         unsigned int  mis = FIX2UINT(_mis);
 182         unsigned int  ins = FIX2UINT(_ins);
 183         unsigned int  del = FIX2UINT(_del);
 184
 185         char         *ss    = s;
 186         int           state = 0;
 187         unsigned int  e     = 0;
 188         VALUE         tuple;
 189
 190         if (pos < len)
 191         {
 192           s += pos;
 193
 194           while (*s)
 195           {
 196             if ((e = backtrack(ss, s, p, mis, ins, del, state)))
 197             {
 198               tuple = rb_ary_new();
 199               rb_ary_push(tuple, INT2FIX((int) (s - ss)));
 200               rb_ary_push(tuple, INT2FIX((int) e - (s - ss)));
 201               return tuple;
 202             }
 203
 204             s++;
 205           }
 206         }
 207
 208         return Qnil;
 209       }
 210     }
 211   end
 212
 213   # Class containing match information.
 214   class Match
 215     attr_reader :pos, :length, :match
 216
 217     def initialize(pos, length, match)
 218       @pos    = pos
 219       @length = length
 220       @match  = match
 221     end
 222
 223     def to_s
 224       "#{pos}:#{length}:#{match}"
 225     end
 226   end
 227 end
 228
 229
 230 __END__