code_ruby/lib/maasha/align/pair.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/matches'
  26 require 'maasha/math_aux'
  27
  28 FACTOR_SCORE_LENGTH =  1.0
  29 FACTOR_SCORE_DIAG   = -1.41
  30 KMER                = 32
  31
  32 # Module with stuff to create a pairwise aligment.
  33 module PairAlign
  34   # Class for creating a pairwise alignment.
  35   class AlignPair
  36     # Class method to create a pairwise alignment of two given Seq objects.
  37     def self.align(q_entry, s_entry)
  38       self.new(q_entry, s_entry)
  39     end
  40
  41     # Method to inialize a pairwise alignment given two Seq objects.
  42     def initialize(q_entry, s_entry)
  43       @q_entry = q_entry
  44       @s_entry = s_entry
  45       @matches = []
  46
  47       @q_entry.seq.downcase!
  48       @s_entry.seq.downcase!
  49
  50       space = Space.new(0, 0, @q_entry.length - 1, @s_entry.length - 1)
  51
  52       align_recurse(@q_entry.seq, @s_entry.seq, space, KMER)
  53       matches_upcase
  54       gaps_insert
  55     end
  56
  57     private
  58
  59     # Method that creates an alignment by chaining matches, which are
  60     # subsequences shared between two sequences. This recursive method
  61     # functions by considering only matches within a given search space. If no
  62     # matches are given these will be located and matches will be included
  63     # depending on a calculated score. New search spaces spanning the spaces
  64     # between the best scoring matches and the search space boundaries will be
  65     # cast and recursed into.
  66     def align_recurse(q_seq, s_seq, space, kmer, matches = [])
  67       matches = matches_select_by_space(matches, space)
  68       matches = matches_select_by_score(matches, space)
  69
  70       while (matches.size == 0 and kmer > 0)
  71         matches = Matches.find(q_seq, s_seq, space.q_min, space.s_min, space.q_max, space.s_max, kmer)
  72         #matches = Mem.find(q_seq, s_seq, kmer, space.q_min, space.s_min, space.q_max, space.s_max)
  73
  74         if @matches.empty?
  75           matches.sort_by! { |m| m.length }
  76         else
  77           matches = matches_select_by_score(matches, space)
  78         end
  79
  80         kmer /= 2
  81       end
  82
  83       if best_match = matches.pop
  84         @matches << best_match
  85
  86         space_left  = Space.new(space.q_min, space.s_min, best_match.q_beg - 1, best_match.s_beg - 1)
  87         space_right = Space.new(best_match.q_end + 1, best_match.s_end + 1, space.q_max, space.s_max)
  88
  89         align_recurse(q_seq, s_seq, space_left, kmer, matches)  unless space_left.empty?
  90         align_recurse(q_seq, s_seq, space_right, kmer, matches) unless space_right.empty?
  91       end
  92     end
  93
  94     # Method to select matches that lies within the search space.
  95     def matches_select_by_space(matches, space)
  96       new_matches = matches.select do |match|
  97         match.q_beg >= space.q_min and
  98         match.s_beg >= space.s_min and
  99         match.q_end <= space.q_max and
 100         match.s_end <= space.s_max
 101       end
 102
 103       new_matches
 104     end
 105
 106     # Method to select matches based on score.
 107     def matches_select_by_score(matches, space)
 108       matches_score(matches, space)
 109
 110       matches.select { |match| match.score > 0 }
 111     end
 112
 113     def matches_score(matches, space)
 114       matches.each do |match|
 115         score_length = match_score_length(match)
 116         score_diag   = match_score_diag(match, space)
 117
 118         match.score = score_length + score_diag
 119       end
 120
 121       matches.sort_by! { |match| match.score }
 122     end
 123
 124     def match_score_length(match)
 125       match.length * FACTOR_SCORE_LENGTH
 126     end
 127
 128     def match_score_diag(match, space)
 129       if space.q_dim > space.s_dim   # s_dim is the narrow end
 130         dist_beg = Math.dist_point2line(match.q_beg,
 131                                         match.s_beg,
 132                                         space.q_min,
 133                                         space.s_min,
 134                                         space.q_min + space.s_dim,
 135                                         space.s_min + space.s_dim)
 136
 137         dist_end = Math.dist_point2line( match.q_beg,
 138                                          match.s_beg,
 139                                          space.q_max - space.s_dim,
 140                                          space.s_max - space.s_dim,
 141                                          space.q_max,
 142                                          space.s_max)
 143       else
 144         dist_beg = Math.dist_point2line( match.q_beg,
 145                                          match.s_beg,
 146                                          space.q_min,
 147                                          space.s_min,
 148                                          space.q_min + space.q_dim,
 149                                          space.s_min + space.q_dim)
 150
 151         dist_end = Math.dist_point2line( match.q_beg,
 152                                          match.s_beg,
 153                                          space.q_max - space.q_dim,
 154                                          space.s_max - space.q_dim,
 155                                          space.q_max,
 156                                          space.s_max)
 157       end
 158
 159       dist_min = dist_beg < dist_end ? dist_beg : dist_end
 160
 161       dist_min * FACTOR_SCORE_DIAG
 162     end
 163
 164     # Method for debugging purposes that upcase matching sequence while non-matches
 165     # sequence is kept in lower case.
 166     def matches_upcase
 167       @matches.each do |match|
 168         @q_entry.seq[match.q_beg .. match.q_end] = @q_entry.seq[match.q_beg .. match.q_end].upcase
 169         @s_entry.seq[match.s_beg .. match.s_end] = @s_entry.seq[match.s_beg .. match.s_end].upcase
 170       end
 171     end
 172
 173     # Method that insert gaps in sequences based on a list of matches and thus
 174     # creating an alignment.
 175     def gaps_insert
 176       @matches.sort_by! { |m| m.q_beg }
 177
 178       q_gaps = 0
 179       s_gaps = 0
 180
 181       match = @matches.first
 182       diff  = (q_gaps + match.q_beg) - (s_gaps + match.s_beg)
 183
 184       if diff < 0
 185         @q_entry.seq.insert(0, "-" * diff.abs)
 186         q_gaps += diff.abs
 187       elsif diff > 0
 188         @s_entry.seq.insert(0, "-" * diff.abs)
 189         s_gaps += diff.abs
 190       end
 191
 192       @matches[1 .. -1].each do |m|
 193         diff = (q_gaps + m.q_beg) - (s_gaps + m.s_beg)
 194
 195         if diff < 0
 196           @q_entry.seq.insert(m.q_beg + q_gaps, "-" * diff.abs)
 197           q_gaps += diff.abs
 198         elsif diff > 0
 199           @s_entry.seq.insert(m.s_beg + s_gaps, "-" * diff.abs)
 200           s_gaps += diff.abs
 201         end
 202       end
 203
 204       diff = @q_entry.length - @s_entry.length
 205
 206       if diff < 0
 207         @q_entry.seq << ("-" * diff.abs)
 208       else
 209         @s_entry.seq << ("-" * diff.abs)
 210       end
 211     end
 212   end
 213
 214   # Class for containing a search space between two sequences q and s.
 215   class Space
 216     attr_reader :q_min, :s_min, :q_max, :s_max
 217
 218     def initialize(q_min, s_min, q_max, s_max)
 219       @q_min = q_min
 220       @s_min = s_min
 221       @q_max = q_max
 222       @s_max = s_max
 223     end
 224
 225     def q_dim
 226       @q_max - @q_min + 1
 227     end
 228
 229     def s_dim
 230       @s_max - @s_min + 1
 231     end
 232
 233     def empty?
 234       if @q_max - @q_min >= 0 and @s_max - @s_min >= 0
 235         return false
 236       end
 237
 238       true
 239     end
 240   end
 241 end
 242