code_ruby/lib/maasha/seq/patternmatcher.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 edit distance.
  29 # Insertions are nucleotides found in the pattern but not in the sequence.
  30 # Deletions are nucleotides found in the sequence but not in the pattern.
  31 #
  32 # Inspired by the paper by Bruno Woltzenlogel Paleo (page 197):
  33 # http://www.logic.at/people/bruno/Papers/2007-GATE-ESSLLI.pdf
  34 module PatternMatcher
  35   # ------------------------------------------------------------------------------
  36   #   str.patmatch(pattern[, pos[, max_edit_distance]])
  37   #   -> Match or nil
  38   #   str.patscan(pattern[, pos[, max_edit_distance]]) { |match|
  39   #     block
  40   #   }
  41   #   -> Match
  42   #
  43   # ------------------------------------------------------------------------------
  44   # Method to iterate through a sequence to locate pattern matches starting from a
  45   # given position and allowing for a maximum edit distance. Matches found in
  46   # block context return the Match object. Otherwise matches are returned in an
  47   # Array.
  48   def patmatch(pattern, pos = 0, max_edit_distance = 0)
  49     self.patscan(pattern, pos, max_edit_distance) do |m|
  50       return m
  51     end
  52   end
  53
  54   # ------------------------------------------------------------------------------
  55   #   str.patscan(pattern[, pos[, max_edit_distance]])
  56   #   -> Array or nil
  57   #   str.patscan(pattern[, pos[, max_edit_distance]]) { |match|
  58   #     block
  59   #   }
  60   #   -> Match
  61   #
  62   # ------------------------------------------------------------------------------
  63   # Method to iterate through a sequence to locate pattern matches starting from a
  64   # given position and allowing for a maximum edit distance. Matches found in
  65   # block context return the Match object. Otherwise matches are returned in an
  66   # Array.
  67   def patscan(pattern, pos = 0, max_edit_distance = 0)
  68     matches = []
  69
  70     while result = match_C(self.seq, self.length, pattern, pattern.length, pos, max_edit_distance)
  71       match = Match.new(*result, self.seq[result[0] ... result[0] + result[1]]);
  72
  73       if block_given?
  74         yield match
  75       else
  76         matches << match
  77       end
  78
  79       pos = match.beg + 1
  80     end
  81
  82     return matches unless block_given?
  83   end
  84
  85   private
  86
  87   inline do |builder|
  88     # Macro for matching nucleotides including ambiguity codes.
  89     builder.prefix %{
  90       #define MAX_PAT 1024
  91     }
  92
  93     builder.prefix %{
  94       #define MATCH(A,B) ((bitmap[A] & bitmap[B]) != 0)
  95     }
  96
  97     builder.prefix %{
  98       typedef struct
  99       {
 100         unsigned int mis;
 101         unsigned int ins;
 102         unsigned int del;
 103         unsigned int ed;
 104       } score;
 105     }
 106
 107     # Bitmap for matching nucleotides including ambiguity codes.
 108     # For each value bits are set from the left: bit pos 1 for A,
 109     # bit pos 2 for T, bit pos 3 for C, and bit pos 4 for G.
 110     builder.prefix %{
 111       char bitmap[256] = {
 112           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 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, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
 117           0, 0, 9,12, 2, 2,13, 3, 0, 6, 0, 0, 0, 0, 0, 0,
 118           0, 1,14, 4,11, 0, 0, 8, 7, 0, 0,10, 0, 5,15, 0,
 119           0, 0, 9,12, 2, 2,13, 3, 0, 6, 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           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       };
 129     }
 130
 131     builder.prefix %{
 132       void vector_init(score *vec, unsigned int vec_len)
 133       {
 134         unsigned int i = 0;
 135
 136         for (i = 1; i < vec_len; i++)
 137         {
 138           vec[i].ins = i;
 139           vec[i].ed  = i;
 140         }
 141       }
 142     }
 143
 144     builder.prefix %{
 145       void vector_print(score *vec, unsigned int vec_len)
 146       {
 147         unsigned int i = 0;
 148
 149         for (i = 0; i < vec_len; i++)
 150         {
 151           printf("i: %d   mis: %d   ins: %d   del: %d   ed: %d\\n", i, vec[i].mis, vec[i].ins, vec[i].del, vec[i].ed);
 152         }
 153
 154         printf("---\\n");
 155       }
 156     }
 157
 158     builder.prefix %{
 159       int match_found(score *vec, unsigned int pat_len, unsigned int max_ed)
 160       {
 161         return (vec[pat_len].ed <= max_ed);
 162       }
 163     }
 164
 165     builder.prefix %{
 166       void vector_update(score *vec, char *seq, char *pat, unsigned int pat_len, unsigned int pos)
 167       {
 168         score diag = vec[0];
 169         score up   = {0, 0, 0, 0};   // insertion
 170         score left = vec[1];            // deletion
 171         score new  = {0, 0, 0, 0};
 172
 173         unsigned int i = 0;
 174
 175         for (i = 0; i < pat_len; i++)
 176         {
 177           if (MATCH(seq[pos], pat[i]))                        // match
 178           {
 179             new = diag;
 180           }
 181           else
 182           {
 183             if (left.ed <= diag.ed && left.ed <= up.ed)       // deletion
 184             {
 185               new = left;
 186               new.del++;
 187             }
 188             else if (diag.ed <= up.ed && diag.ed <= left.ed)   // mismatch
 189             {
 190               new = diag;
 191               new.mis++;
 192             }
 193             else if (up.ed <= diag.ed && up.ed <= left.ed)      // insertion
 194             {
 195               new = up;
 196               new.ins++;
 197             }
 198             else
 199             {
 200               printf("This should not happen\\n");
 201               exit(1);
 202             }
 203
 204             new.ed++;
 205           }
 206
 207           diag = vec[i + 1];
 208           up   = new;
 209           left = vec[i + 2];
 210
 211           vec[i + 1] = new;
 212         }
 213       }
 214     }
 215
 216     builder.c %{
 217       VALUE match_C(
 218         VALUE _seq,       // Sequence
 219         VALUE _seq_len,   // Sequence length
 220         VALUE _pat,       // Pattern
 221         VALUE _pat_len,   // Pattern length
 222         VALUE _pos,       // Offset position
 223         VALUE _max_ed     // Maximum edit distance
 224       )
 225       {
 226         char         *seq     = (char *) StringValuePtr(_seq);
 227         char         *pat     = (char *) StringValuePtr(_pat);
 228         unsigned int  seq_len = FIX2UINT(_seq_len);
 229         unsigned int  pat_len = FIX2UINT(_pat_len);
 230         unsigned int  pos     = FIX2UINT(_pos);
 231         unsigned int  max_ed  = FIX2UINT(_max_ed);
 232
 233         score         vec[MAX_PAT] = {0};
 234         unsigned int  vec_len      = pat_len + 1;
 235         unsigned int  match_beg    = 0;
 236         unsigned int  match_len    = 0;
 237
 238         VALUE         match_ary;
 239
 240         vector_init(vec, vec_len);
 241
 242         while (pos < seq_len)
 243         {
 244           vector_update(vec, seq, pat, pat_len, pos);
 245
 246           if (match_found(vec, pat_len, max_ed))
 247           {
 248             match_len = pat_len - vec[pat_len].ins + vec[pat_len].del;
 249             match_beg = pos - match_len + 1;
 250
 251             match_ary = rb_ary_new();
 252             rb_ary_push(match_ary, INT2FIX(match_beg));
 253             rb_ary_push(match_ary, INT2FIX(match_len));
 254             rb_ary_push(match_ary, INT2FIX(vec[pat_len].mis));
 255             rb_ary_push(match_ary, INT2FIX(vec[pat_len].ins));
 256             rb_ary_push(match_ary, INT2FIX(vec[pat_len].del));
 257
 258             return match_ary;
 259           }
 260
 261           pos++;
 262         }
 263
 264         return Qfalse;  // no match
 265       }
 266     }
 267   end
 268
 269   class Match
 270     attr_accessor :beg, :length, :mis, :ins, :del, :match
 271
 272     def initialize(beg, length, mis, ins, del, match)
 273       @beg    = beg
 274       @length = length
 275       @mis    = mis
 276       @ins    = ins
 277       @del    = del
 278       @match  = match
 279     end
 280   end
 281 end
 282
 283
 284 __END__