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