code_ruby/lib/maasha/align/match.rb

   1 # Copyright (C) 2007-2012 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 # Class containing methods for located all non-redundant Maximally Extended Matches (MEMs)
  26 # between two strings.
  27 class Matches
  28   # Class method to located all non-redudant MEMs bewteen two sequences within a given
  29   # search space and seeded with kmers of a given size.
  30   def self.find(q_seq, s_seq, q_min, s_min, q_max, s_max, kmer)
  31     m = self.new
  32     m.matches_find(q_seq, s_seq, q_min, s_min, q_max, s_max, kmer)
  33   end
  34
  35   # Method that finds all maximally expanded non-redundant matches shared
  36   # between two sequences inside a given search space.
  37   def matches_find(q_seq, s_seq, q_min, s_min, q_max, s_max, kmer)
  38     matches   = []
  39     redundant = Hash.new { |h, k| h[k] = [] }
  40
  41     s_index = index_seq(s_seq, s_min, s_max, kmer)
  42
  43     q_pos = q_min
  44
  45     while q_pos <= q_max - kmer + 1
  46       q_oligo = q_seq[q_pos ... q_pos + kmer]
  47
  48       s_index[q_oligo].each do |s_pos|
  49         match = Match.new(q_pos, s_pos, kmer)
  50
  51         unless match_redundant?(redundant, match)
  52           match_expand(match, q_seq, s_seq, q_min, s_min, q_max, s_max)
  53           matches << match
  54
  55           match_redundant_add(redundant, match)
  56         end
  57       end
  58
  59       q_pos += 1
  60     end
  61
  62     matches
  63   end
  64
  65   # Method that indexes a sequence within a given interval such that the
  66   # index contains all oligos of a given kmer size and the positions where
  67   # this oligo was located.
  68   def index_seq(seq, min, max, kmer, step = 1)
  69     index_hash = Hash.new { |h, k| h[k] = [] }
  70
  71     pos = min
  72
  73     while pos <= max - kmer + 1
  74       oligo = seq[pos ... pos + kmer]
  75       index_hash[oligo] << pos
  76
  77       pos += step
  78     end
  79
  80     index_hash
  81   end
  82
  83   # Method to check if a match is redundant.
  84   def match_redundant?(redundant, match)
  85     redundant[match.q_beg].each do |s_interval|
  86       if s_interval.include? match.s_beg and s_interval.include? match.s_end   # TODO test if include? is slow
  87         return true
  88       end
  89     end
  90
  91     false
  92   end
  93
  94   # Method that adds a match to the redundancy index.
  95   def match_redundant_add(redundant, match)
  96     (match.q_beg .. match.q_end).each do |q|
  97       redundant[q] << (match.s_beg .. match.s_end)
  98     end
  99   end
 100
 101   # Method that expands a match as far as possible to the left and right.
 102   def match_expand(match, q_seq, s_seq, q_min, s_min, q_max, s_max)
 103     match_expand_left(match, q_seq, s_seq, q_min, s_min, q_max, s_max)
 104     match_expand_right(match, q_seq, s_seq, q_min, s_min, q_max, s_max)
 105
 106     match
 107   end
 108
 109   # Method that expands a match as far as possible to the left.
 110   def match_expand_left(match, q_seq, s_seq, q_min, s_min, q_max, s_max)
 111     while match.q_beg > q_min and
 112           match.s_beg > s_min and
 113           q_seq[match.q_beg - 1] == s_seq[match.s_beg - 1]
 114       match.q_beg  -= 1
 115       match.s_beg  -= 1
 116       match.length += 1
 117     end
 118
 119     match
 120   end
 121
 122   # Method that expands a match as far as possible to the right.
 123   def match_expand_right(match, q_seq, s_seq, q_min, s_min, q_max, s_max)
 124     while match.q_end < q_max and
 125           match.s_end < s_max and
 126           q_seq[match.q_end + 1] == s_seq[match.s_end + 1]
 127       match.length += 1
 128     end
 129
 130     match
 131   end
 132
 133   # Class for containing a match between two sequences q and s.
 134   class Match
 135     attr_accessor :q_beg, :s_beg, :length, :score
 136
 137     def initialize(q_beg, s_beg, length, score = 0.0)
 138       @q_beg  = q_beg
 139       @s_beg  = s_beg
 140       @length = length
 141       @score  = score
 142     end
 143
 144     def q_end
 145       @q_beg + @length - 1
 146     end
 147
 148     def s_end
 149       @s_beg + @length - 1
 150     end
 151
 152     def to_s(seq = nil)
 153       s = "q: #{@q_beg} #{q_end} s: #{@s_beg} #{s_end} l: #{@length} s: #{@score}"
 154       s << " seq: #{seq[@q_beg .. q_end]}" if seq
 155       s
 156     end
 157   end
 158 end