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