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     while result = track_C(self.seq, self.length, pattern, offset, max_mismatches, max_insertions, max_deletions)
  68       pos, length, mis, ins, del = result
  69
  70       match = Match.new(pos, length, mis, ins, del, self.seq[pos ... pos + length])
  71
  72       if block_given?
  73         yield match
  74       else
  75         matches << match
  76       end
  77
  78       offset = match.pos + 1
  79     end
  80
  81     return matches.empty? ? nil : matches unless block_given?
  82   end
  83
  84   private
  85
  86   inline do |builder|
  87     # Macro for matching nucleotides including ambiguity codes.
  88     builder.prefix %{
  89       #define MATCH(A,B) ((bitmap[A] & bitmap[B]) != 0)
  90     }
  91
  92     # Defining a struct for returning match information.
  93     builder.prefix %{
  94       typedef struct
  95       {
  96         unsigned int end;
  97         unsigned int mis;
  98         unsigned int ins;
  99         unsigned int del;
 100       } match;
 101     }
 102
 103     # Bitmap for matching nucleotides including ambiguity codes.
 104     # For each value bits are set from the left: bit pos 1 for A,
 105     # bit pos 2 for T, bit pos 3 for C, and bit pos 4 for G.
 106     builder.prefix %{
 107       char bitmap[256] = {
 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, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
 113           0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
 114           0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
 115           0, 0, 9,12, 2, 2,13, 3, 0, 6, 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, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 118           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 119           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 120           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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       };
 125     }
 126
 127     # Backtrack algorithm for matching a pattern (p) starting in a sequence (s) allowing for mis
 128     # mismatches, ins insertions and del deletions. ss is the start of the sequence, used only for
 129     # reporting the match endpoints.
 130     builder.prefix %{
 131       match* backtrack(
 132         char         *ss,   // Sequence start
 133         char         *s,    // Sequence
 134         char         *p,    // Pattern
 135         unsigned int  mis,  // Max mismatches
 136         unsigned int  ins,  // Max insertions
 137         unsigned int  del,  // Max deletions
 138         int state           // 1 deletion, -1 insertion, 0 match or mismatch
 139       )
 140       {
 141         match found = {0, 0, 0, 0};
 142         match *f = &found;
 143
 144         while (*s && MATCH(*s, *p)) ++s, ++p;    // OK to always match longest segment
 145
 146         if (!*p)
 147         {
 148           f->end = (unsigned int) (s - ss);
 149           f->mis = mis;
 150           f->ins = ins;
 151           f->del = del;
 152
 153           return f;
 154         }
 155         else
 156         {
 157           if (mis && *s && *p            && (f = backtrack(ss, s + 1, p + 1, mis - 1, ins, del, 0))) return f;
 158           if (ins && *p && (state != 1)  && (f = backtrack(ss, s, p + 1, mis, ins - 1, del, -1)))    return f;
 159           if (del && *s && (state != -1) && (f = backtrack(ss, s + 1, p, mis, ins, del - 1, 1)))     return f;
 160         }
 161
 162         return NULL;
 163       }
 164     }
 165
 166     # Find pattern (p) in a sequence (s) starting at pos, with at most mis mismatches, ins
 167     # insertions and del deletions.
 168     builder.c %{
 169       VALUE track_C(
 170         VALUE _s,         // Sequence
 171         VALUE _len,       // Sequence length
 172         VALUE _p,         // Pattern
 173         VALUE _pos,       // Start position
 174         VALUE _max_mis,   // Maximum mismatches
 175         VALUE _max_ins,   // Maximum insertions
 176         VALUE _max_del    // Maximum deletions
 177       )
 178       {
 179         char         *s       = StringValuePtr(_s);
 180         char         *p       = StringValuePtr(_p);
 181         unsigned int  len     = FIX2UINT(_len);
 182         unsigned int  pos     = FIX2UINT(_pos);
 183         unsigned int  max_mis = FIX2UINT(_max_mis);
 184         unsigned int  max_ins = FIX2UINT(_max_ins);
 185         unsigned int  max_del = FIX2UINT(_max_del);
 186
 187         char         *ss    = s;
 188         match        *f     = NULL;
 189         unsigned int  mbeg  = 0;
 190         unsigned int  mlen  = 0;
 191         unsigned int  mmis  = 0;
 192         unsigned int  mins  = 0;
 193         unsigned int  mdel  = 0;
 194         int           state = 0;
 195         VALUE         ary   = 0;
 196
 197         if (pos < len) // FIXME: pos < len - plen + 1
 198         {
 199           s += pos;
 200
 201           while (*s)
 202           {
 203             if ((f = backtrack(ss, s, p, max_mis, max_ins, max_del, state)))
 204             {
 205               mbeg = (s - ss);
 206               mlen = f->end - mbeg;
 207               mmis = max_mis - f->mis;
 208               mins = max_ins - f->ins;
 209               mdel = max_del - f->del;
 210
 211               ary = rb_ary_new();
 212               rb_ary_push(ary, INT2FIX(mbeg));
 213               rb_ary_push(ary, INT2FIX(mlen));
 214               rb_ary_push(ary, INT2FIX(mmis));
 215               rb_ary_push(ary, INT2FIX(mins));
 216               rb_ary_push(ary, INT2FIX(mdel));
 217               return ary;
 218             }
 219
 220             s++;
 221           }
 222         }
 223
 224         return Qnil;
 225       }
 226     }
 227   end
 228
 229   # Class containing match information.
 230   class Match
 231     attr_reader :pos, :length, :mis, :ins, :del, :match
 232
 233     def initialize(pos, length, mis, ins, del, match)
 234       @pos    = pos
 235       @length = length
 236       @mis    = mis
 237       @ins    = ins
 238       @del    = del
 239       @match  = match
 240     end
 241
 242     def to_s
 243       "#{@pos}:#{@length}:#{@mis}:#{@ins}:#{@del}:#{@match}"
 244     end
 245   end
 246 end
 247
 248
 249 __END__