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