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