code_ruby/lib/maasha/align/matches.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 require 'maasha/align/match'
  26
  27 # Class containing methods for located all non-redundant Maximally Extended Matches (MEMs)
  28 # between two strings.
  29 class Matches
  30   # Class method to located all non-redudant MEMs bewteen two sequences within a given
  31   # search space and seeded with kmers of a given size.
  32   def self.find(q_seq, s_seq, q_min, s_min, q_max, s_max, kmer)
  33     m = self.new
  34     m.matches_find(q_seq, s_seq, q_min, s_min, q_max, s_max, kmer)
  35   end
  36
  37   # Method that finds all maximally expanded non-redundant matches shared
  38   # between two sequences inside a given search space.
  39   def matches_find(q_seq, s_seq, q_min, s_min, q_max, s_max, kmer)
  40     matches   = []
  41     redundant = {}
  42
  43     s_index = index_seq(s_seq, s_min, s_max, kmer)
  44
  45     q_pos = q_min
  46
  47     while q_pos <= q_max - kmer + 1
  48       q_oligo = q_seq[q_pos ... q_pos + kmer]
  49
  50       s_index[q_oligo].each do |s_pos|
  51         match = Match.new(q_pos, s_pos, kmer)
  52
  53         unless redundant[match.q_beg * 1_000_000_000 + match.s_beg]
  54           match.expand(q_seq, s_seq, q_min, s_min, q_max, s_max)
  55           matches << match
  56
  57           (0 ... match.length).each do |j|
  58             redundant[(match.q_beg + j) * 1_000_000_000 + match.s_beg + j] = true
  59           end
  60         end
  61       end
  62
  63       q_pos += 1
  64     end
  65
  66     matches
  67   end
  68
  69   # Method that indexes a sequence within a given interval such that the
  70   # index contains all oligos of a given kmer size and the positions where
  71   # this oligo was located.
  72   def index_seq(seq, min, max, kmer, step = 1)
  73     index_hash = Hash.new { |h, k| h[k] = [] }
  74
  75     pos = min
  76
  77     while pos <= max - kmer + 1
  78       oligo = seq[pos ... pos + kmer]
  79       index_hash[oligo] << pos
  80
  81       pos += step
  82     end
  83
  84     index_hash
  85   end
  86 end