code_ruby/lib/maasha/align/match.rb

   1 # Copyright (C) 2007-2013 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 require 'maasha/math_aux'
  27
  28 # Class for describing a match between two sequences q and s.
  29 class Match
  30   attr_accessor :q_beg, :s_beg, :length, :score
  31
  32   def initialize(q_beg, s_beg, length)
  33     @q_beg  = q_beg
  34     @s_beg  = s_beg
  35     @length = length
  36     @score  = 0.0
  37   end
  38
  39   def q_end
  40     @q_beg + @length - 1
  41   end
  42
  43   def s_end
  44     @s_beg + @length - 1
  45   end
  46
  47   # Method to expand a match as far as possible to the left and right
  48   # within a given search space.
  49   def expand(q_seq, s_seq, q_space_beg, s_space_beg, q_space_end, s_space_end)
  50     length = expand_left_C(q_seq, s_seq, @q_beg, @s_beg, q_space_beg, s_space_beg)
  51     @length += length
  52     @q_beg  -= length
  53     @s_beg  -= length
  54
  55     length = expand_right_C(q_seq, s_seq, q_end, s_end, q_space_end, s_space_end)
  56     @length += length
  57
  58     self
  59   end
  60
  61   def redundant?(matches)
  62     # Speed-up with binary search
  63
  64     matches.each do |m|
  65       if Math.dist_point2line(self.q_beg, self.s_beg, m.q_beg, m.s_beg, m.q_end, m.s_end) == 0
  66         return true
  67       end
  68     end
  69
  70     false
  71   end
  72
  73   def to_s(seq = nil)
  74     s = "q: #{@q_beg}, s: #{@s_beg}, l: #{@length}, s: #{@score}"
  75     s << " #{seq[@q_beg .. q_end]}" if seq
  76
  77     s
  78   end
  79
  80   private
  81
  82   inline do |builder|
  83     # Method to expand a match as far as possible to the left within a given
  84     # search space.
  85     builder.c %{
  86       VALUE expand_left_C(
  87         VALUE _q_seq,
  88         VALUE _s_seq,
  89         VALUE _q_beg,
  90         VALUE _s_beg,
  91         VALUE _q_space_beg,
  92         VALUE _s_space_beg
  93       )
  94       {
  95         unsigned char *q_seq       = (unsigned char *) StringValuePtr(_q_seq);
  96         unsigned char *s_seq       = (unsigned char *) StringValuePtr(_s_seq);
  97         unsigned int   q_beg       = NUM2UINT(_q_beg);
  98         unsigned int   s_beg       = NUM2UINT(_s_beg);
  99         unsigned int   q_space_beg = NUM2UINT(_q_space_beg);
 100         unsigned int   s_space_beg = NUM2UINT(_s_space_beg);
 101
 102         unsigned int   len = 0;
 103
 104         while (q_beg > q_space_beg && s_beg > s_space_beg && q_seq[q_beg - 1] == s_seq[s_beg - 1])
 105         {
 106           q_beg--;
 107           s_beg--;
 108           len++;
 109         }
 110
 111         return UINT2NUM(len);
 112       }
 113     }
 114
 115     builder.c %{
 116       VALUE expand_right_C(
 117         VALUE _q_seq,
 118         VALUE _s_seq,
 119         VALUE _q_end,
 120         VALUE _s_end,
 121         VALUE _q_space_end,
 122         VALUE _s_space_end
 123       )
 124       {
 125         unsigned char *q_seq       = (unsigned char *) StringValuePtr(_q_seq);
 126         unsigned char *s_seq       = (unsigned char *) StringValuePtr(_s_seq);
 127         unsigned int   q_end       = NUM2UINT(_q_end);
 128         unsigned int   s_end       = NUM2UINT(_s_end);
 129         unsigned int   q_space_end = NUM2UINT(_q_space_end);
 130         unsigned int   s_space_end = NUM2UINT(_s_space_end);
 131
 132         unsigned int   len = 0;
 133
 134         while (q_end + 1 <= q_space_end && s_end + 1 <= s_space_end && q_seq[q_end + 1] == s_seq[s_end + 1])
 135         {
 136           q_end++;
 137           s_end++;
 138           len++;
 139         }
 140
 141         return UINT2NUM(len);
 142       }
 143     }
 144   end
 145 end