code_ruby/lib/maasha/mum.rb

   1 #!/usr/bin/env ruby
   2
   3 require 'inline'
   4 require 'maasha/seq'
   5 require 'pp'
   6 #require 'profile'
   7
   8 BYTES_IN_BIN   = 4
   9 BYTES_IN_POS   = 4
  10 BYTES_IN_OLIGO = BYTES_IN_BIN + BYTES_IN_POS
  11
  12 class MUMs
  13   include Enumerable
  14
  15   def self.find(q_seq, s_seq, space, kmer_size)
  16     m = self.new(q_seq, s_seq, space, kmer_size)
  17     m.to_a
  18   end
  19
  20   def initialize(q_seq, s_seq, space, kmer_size)
  21     q_space_beg = space.q_min
  22     q_space_end = space.q_max
  23     s_space_beg = space.s_min
  24     s_space_end = space.s_max
  25
  26     @mums = []
  27
  28     q_index = Index.new(q_space_beg, q_space_end, q_seq.seq, kmer_size, kmer_size)
  29     s_index = Index.new(s_space_beg, s_space_end, s_seq.seq, kmer_size, 1)
  30
  31     find_mums(q_index, s_index, kmer_size)
  32
  33     @mums.each do |match|
  34       match.expand(q_seq.seq, s_seq.seq, q_space_beg, q_space_end, s_space_beg, s_space_end)
  35     end
  36   end
  37
  38   def each
  39     @mums.each do |match|
  40       yield match
  41     end
  42   end
  43
  44   private
  45
  46   def find_mums(q_index, s_index, kmer_size)
  47     q_index.each do |q_oligo|
  48       if s_pos = s_index[q_oligo]
  49         q_pos = q_oligo.unpack("II").last
  50
  51         last_match = @mums.last
  52         new_match  = Match.new(q_pos, s_pos, kmer_size)
  53
  54         if last_match and
  55            last_match.q_beg + last_match.length == new_match.q_beg and
  56            last_match.s_beg + last_match.length == new_match.s_beg
  57
  58            last_match.length += kmer_size
  59         else
  60           @mums << new_match
  61         end
  62       end
  63     end
  64   end
  65
  66   class Index
  67     def initialize(space_beg, space_end, seq, kmer_size, step_size)
  68       index_size = (space_end - space_beg) * BYTES_IN_OLIGO
  69       index      = "\0" * index_size
  70       index_size = create_index_C(space_beg, space_end, seq, kmer_size, step_size, index)
  71       sort_by_oligo_C(index, index_size)
  72       @index_size = uniq_oligo_C(index, index_size)
  73
  74       @index = index[0 ... @index_size * BYTES_IN_OLIGO]
  75     end
  76
  77     def each
  78       (0 ... @index_size - 1).each do |i|
  79         yield @index[BYTES_IN_OLIGO * i ... BYTES_IN_OLIGO * i + BYTES_IN_OLIGO]
  80       end
  81     end
  82
  83     def [](oligo)
  84       search_index_C(@index, @index_size, oligo)
  85     end
  86
  87     private
  88
  89     inline do |builder|
  90       builder.prefix %{
  91         unsigned char nuc2bin[256] = {
  92           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  93           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  94           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  95           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  96           0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0,
  97           0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
  98           0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0,
  99           0, 0, 0, 0, 1, 1, 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, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 102           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 103           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
 104           0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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         };
 109       }
 110
 111       builder.prefix %{
 112         typedef struct
 113         {
 114           unsigned int bin;
 115           unsigned int pos;
 116         } oligo;
 117       }
 118
 119       builder.prefix %{
 120         int oligo_cmp(const void *a, const void *b)
 121         {
 122           oligo *ia = (oligo *) a;
 123           oligo *ib = (oligo *) b;
 124
 125           return (int) (ia->bin - ib->bin);
 126         }
 127       }
 128
 129       builder.c %{
 130         VALUE create_index_C(
 131           VALUE _space_beg,
 132           VALUE _space_end,
 133           VALUE _seq,
 134           VALUE _kmer_size,
 135           VALUE _step_size,
 136           VALUE _index
 137         )
 138         {
 139           unsigned int   space_beg = NUM2UINT(_space_beg);
 140           unsigned int   space_end = NUM2UINT(_space_end);
 141           unsigned char *seq       = (unsigned char *) StringValuePtr(_seq);
 142           unsigned int   kmer_size = NUM2UINT(_kmer_size);
 143           unsigned int   step_size = NUM2UINT(_step_size);
 144           oligo         *index     = (oligo *) StringValuePtr(_index);
 145
 146           oligo        new  = {0, 0};
 147           unsigned int mask = (1 << (2 * kmer_size)) - 1;
 148           unsigned int bin  = 0;
 149           unsigned int size = 0;
 150           unsigned int i    = 0;
 151
 152           for (i = space_beg; i < space_beg + kmer_size; i++)
 153           {
 154             bin <<= 2;
 155             bin |= nuc2bin[seq[i]];
 156           }
 157
 158           new.bin       = bin;
 159           new.pos       = i - kmer_size;
 160           index[size++] = new;
 161
 162           while (i <= space_end)
 163           {
 164             bin <<= 2;
 165             bin |= nuc2bin[seq[i]];
 166
 167             i++;
 168
 169             if ((i % step_size) == 0) {
 170               new.bin       = (bin & mask);
 171               new.pos       = i - kmer_size;
 172               index[size++] = new;
 173             }
 174           }
 175
 176           return UINT2NUM(size);
 177         }
 178       }
 179
 180       builder.c %{
 181         void sort_by_oligo_C(
 182           VALUE _index,
 183           VALUE _index_size
 184         )
 185         {
 186           oligo        *index      = (oligo *) StringValuePtr(_index);
 187           unsigned int  index_size = NUM2UINT(_index_size);
 188
 189           qsort(index, index_size, sizeof(oligo), oligo_cmp);
 190         }
 191       }
 192
 193       builder.c %{
 194         VALUE uniq_oligo_C(
 195           VALUE _index,
 196           VALUE _index_size
 197         )
 198         {
 199           oligo        *index      = (oligo *) StringValuePtr(_index);
 200           unsigned int  index_size = NUM2UINT(_index_size);
 201
 202           unsigned int  new_index_size = 0;
 203           unsigned int  i              = 0;
 204           unsigned int  j              = 0;
 205
 206           for (i = 1; i < index_size; i++)
 207           {
 208             if (index[i].bin != index[j].bin)
 209             {
 210               j++;
 211               index[j] = index[i]; // Move it to the front
 212             }
 213           }
 214
 215           new_index_size = j + 1;
 216
 217           return UINT2NUM(new_index_size);
 218         }
 219       }
 220
 221       builder.c %{
 222         VALUE search_index_C(
 223           VALUE _index,
 224           VALUE _index_size,
 225           VALUE _item
 226         )
 227         {
 228           oligo        *index      = (oligo *) StringValuePtr(_index);
 229           unsigned int  index_size = NUM2UINT(_index_size);
 230           oligo        *item       = (oligo *) StringValuePtr(_item);
 231           oligo        *result     = NULL;
 232
 233           result = (oligo *) bsearch(item, index, index_size, sizeof(oligo), oligo_cmp);
 234
 235           if (result == NULL)
 236             return Qfalse;
 237           else
 238             return UINT2NUM(result->pos);
 239         }
 240       }
 241     end
 242   end
 243
 244   class Match
 245     attr_accessor :length, :score
 246     attr_reader :q_beg, :s_beg
 247
 248     def initialize(q_beg, s_beg, length)
 249       @q_beg  = q_beg
 250       @s_beg  = s_beg
 251       @length = length
 252       @score  = 0.0
 253     end
 254
 255     def q_end
 256       @q_beg + @length - 1
 257     end
 258
 259     def s_end
 260       @s_beg + @length - 1
 261     end
 262
 263     # Method to expand a match as far as possible to the left and right
 264     # within a given search space.
 265     def expand(q_seq, s_seq, q_space_beg, q_space_end, s_space_beg, s_space_end)
 266       length = expand_left_C(q_seq, s_seq, @q_beg, @s_beg, q_space_beg, s_space_beg)
 267       @length += length
 268       @q_beg  -= length
 269       @s_beg  -= length
 270
 271       length = expand_right_C(q_seq, s_seq, q_end, s_end, q_space_end, s_space_end)
 272       @length += length
 273     end
 274
 275     private
 276
 277     # Method to expand a match as far as possible to the right within a given
 278     # search space.
 279     def expand_right(q_seq, s_seq, q_space_end, s_space_end)
 280       while q_end + 1 <= q_space_end and s_end + 1 <= s_space_end and q_seq[q_end + 1] == s_seq[s_end + 1]
 281         @length += 1
 282       end
 283     end
 284
 285     inline do |builder|
 286       # Method to expand a match as far as possible to the left within a given
 287       # search space.
 288       builder.c %{
 289         VALUE expand_left_C(
 290           VALUE _q_seq,
 291           VALUE _s_seq,
 292           VALUE _q_beg,
 293           VALUE _s_beg,
 294           VALUE _q_space_beg,
 295           VALUE _s_space_beg
 296         )
 297         {
 298           unsigned char *q_seq       = (unsigned char *) StringValuePtr(_q_seq);
 299           unsigned char *s_seq       = (unsigned char *) StringValuePtr(_s_seq);
 300           unsigned int   q_beg       = NUM2UINT(_q_beg);
 301           unsigned int   s_beg       = NUM2UINT(_s_beg);
 302           unsigned int   q_space_beg = NUM2UINT(_q_space_beg);
 303           unsigned int   s_space_beg = NUM2UINT(_s_space_beg);
 304
 305           unsigned int   len = 0;
 306
 307           while (q_beg > q_space_beg && s_beg > s_space_beg && q_seq[q_beg - 1] == s_seq[s_beg - 1])
 308           {
 309             q_beg--;
 310             s_beg--;
 311             len++;
 312           }
 313
 314           return UINT2NUM(len);
 315         }
 316       }
 317
 318       builder.c %{
 319         VALUE expand_right_C(
 320           VALUE _q_seq,
 321           VALUE _s_seq,
 322           VALUE _q_end,
 323           VALUE _s_end,
 324           VALUE _q_space_end,
 325           VALUE _s_space_end
 326         )
 327         {
 328           unsigned char *q_seq       = (unsigned char *) StringValuePtr(_q_seq);
 329           unsigned char *s_seq       = (unsigned char *) StringValuePtr(_s_seq);
 330           unsigned int   q_end       = NUM2UINT(_q_end);
 331           unsigned int   s_end       = NUM2UINT(_s_end);
 332           unsigned int   q_space_end = NUM2UINT(_q_space_end);
 333           unsigned int   s_space_end = NUM2UINT(_s_space_end);
 334
 335           unsigned int   len = 0;
 336
 337           while (q_end + 1 <= q_space_end && s_end + 1 <= s_space_end && q_seq[q_end + 1] == s_seq[s_end + 1])
 338           {
 339             q_end++;
 340             s_end++;
 341             len++;
 342           }
 343
 344           return UINT2NUM(len);
 345         }
 346       }
 347     end
 348   end
 349 end