code_ruby/lib/maasha/align/mem.rb

   1 #!/usr/bin/env ruby
   2
   3 require 'maasha/align/match'
   4
   5 class Mem
   6   def self.find(q_seq, s_seq, kmer, q_space_beg, s_space_beg, q_space_end, s_space_end)
   7     m = self.new(q_seq, s_seq, kmer, q_space_beg, s_space_beg, q_space_end, s_space_end)
   8     m.find_mems
   9   end
  10
  11   def initialize(q_seq, s_seq, kmer, q_space_beg, s_space_beg, q_space_end, s_space_end)
  12     @q_seq       = q_seq.downcase
  13     @s_seq       = s_seq.downcase
  14     @kmer        = kmer
  15     @q_space_beg = q_space_beg
  16     @s_space_beg = s_space_beg
  17     @q_space_end = q_space_end
  18     @s_space_end = s_space_end
  19   end
  20
  21   def find_mems
  22     mask      = (1 << (2 * @kmer)) - 1;
  23     matches   = []
  24     pos_ary   = []
  25     map_ary   = []
  26     redundant = {}
  27
  28     index_seq(@q_seq, @q_space_beg, @q_space_end, @kmer, pos_ary, map_ary)
  29
  30     i   = @s_space_beg
  31     bin = 0
  32
  33     while i < @s_space_beg + @kmer
  34       bin <<= 2
  35
  36       case @s_seq[i]
  37       when 'a' then bin |= 0
  38       when 't' then bin |= 1
  39       when 'u' then bin |= 1
  40       when 'c' then bin |= 2
  41       when 'g' then bin |= 3
  42       else
  43         raise
  44       end
  45
  46       i += 1
  47     end
  48
  49     i =  @s_space_beg + @kmer
  50
  51     while i <= @s_space_end
  52       pos = pos_ary[bin & mask]
  53
  54       while pos
  55         match = Match.new(i - @kmer, pos, @kmer)
  56
  57         unless redundant[match.q_beg * 1_000_000_000 + match.s_beg]
  58           match.expand(@q_seq, @s_seq, 0, 0, @q_seq.length - 1, @s_seq.length - 1)
  59
  60           (0 ... match.length).each do |j|
  61             redundant[(match.q_beg + j) * 1_000_000_000 + match.s_beg + j] = true
  62           end
  63
  64           matches << match
  65         end
  66
  67         pos = map_ary[pos]
  68       end
  69
  70       bin <<= 2
  71
  72       case @s_seq[i]
  73       when 'a' then bin |= 0
  74       when 't' then bin |= 1
  75       when 'u' then bin |= 1
  76       when 'c' then bin |= 2
  77       when 'g' then bin |= 3
  78       else
  79         raise
  80       end
  81
  82       i += 1
  83     end
  84
  85     pos = pos_ary[bin & mask]
  86
  87     while pos
  88       match = Match.new(i - @kmer, pos, @kmer)
  89
  90       unless redundant[match.q_beg * 1_000_000_000 + match.s_beg]
  91         match.expand(@q_seq, @s_seq, 0, 0, @q_seq.length - 1, @s_seq.length - 1)
  92
  93         (0 ... match.length).each do |j|
  94           redundant[(match.q_beg + j) * 1_000_000_000 + match.s_beg + j] = true
  95         end
  96
  97         matches << match
  98       end
  99
 100       pos = map_ary[pos]
 101     end
 102
 103     matches
 104   end
 105
 106   def index_seq(seq, start, stop, kmer, pos_ary, map_ary)
 107     bin  = 0
 108     mask = (1 << (2 * kmer)) - 1;
 109
 110     i = start
 111
 112     while i < start + kmer
 113       bin <<= 2
 114
 115       case seq[i]
 116       when 'a' then bin |= 0
 117       when 't' then bin |= 1
 118       when 'u' then bin |= 1
 119       when 'c' then bin |= 2
 120       when 'g' then bin |= 3
 121       end
 122
 123       i += 1
 124     end
 125
 126     i = start + kmer
 127
 128     while i <= stop
 129       key = bin & mask
 130
 131       map_ary[i - kmer] = pos_ary[key]
 132       pos_ary[key]      = i - kmer
 133
 134       bin <<= 2
 135
 136       case seq[i]
 137       when 'a' then bin |= 0
 138       when 't' then bin |= 1
 139       when 'u' then bin |= 1
 140       when 'c' then bin |= 2
 141       when 'g' then bin |= 3
 142       else
 143         raise
 144       end
 145
 146       i += 1
 147     end
 148
 149     key = bin & mask
 150
 151     map_ary[i - kmer] = pos_ary[key]
 152     pos_ary[key]      = i - kmer
 153   end
 154 end
 155
 156
 157 __END__