From 08f02c04012e697b74c4da03bcda0f468d3a3dfe Mon Sep 17 00:00:00 2001
From: martinahansen <martinahansen@74ccb610-7750-0410-82ae-013aeee3265d>
Date: Wed, 12 Sep 2012 19:43:13 +0000
Subject: [PATCH] worked on pair aligner

git-svn-id: http://biopieces.googlecode.com/svn/trunk@1929 74ccb610-7750-0410-82ae-013aeee3265d
---
 code_ruby/lib/maasha/align/match.rb |  4 ++--
 code_ruby/lib/maasha/align/pair.rb  | 34 +++++++++++++++++++----------
 2 files changed, 25 insertions(+), 13 deletions(-)

diff --git a/code_ruby/lib/maasha/align/match.rb b/code_ruby/lib/maasha/align/match.rb
index e4d5fe0..c17ec7f 100644
--- a/code_ruby/lib/maasha/align/match.rb
+++ b/code_ruby/lib/maasha/align/match.rb
@@ -60,7 +60,7 @@ class Matches
   # Method that indexes a sequence within a given interval such that the
   # index contains all oligos of a given kmer size and the positions where
   # this oligo was located.
-  def index_seq(seq, min, max, kmer)
+  def index_seq(seq, min, max, kmer, step = 1)
     index_hash = Hash.new { |h, k| h[k] = [] }
 
     pos = min
@@ -69,7 +69,7 @@ class Matches
       oligo = seq[pos ... pos + kmer]
       index_hash[oligo] << pos
 
-      pos += 1
+      pos += step
     end
 
     index_hash
diff --git a/code_ruby/lib/maasha/align/pair.rb b/code_ruby/lib/maasha/align/pair.rb
index 6e2dcad..c17dba6 100644
--- a/code_ruby/lib/maasha/align/pair.rb
+++ b/code_ruby/lib/maasha/align/pair.rb
@@ -26,7 +26,7 @@ require 'maasha/align/match'
 require 'maasha/math_aux'
 
 FACTOR_SCORE_LENGTH =  1.0
-FACTOR_SCORE_DIAG   = -1.42
+FACTOR_SCORE_DIAG   = -1.41
 
 # Module with stuff to create a pairwise aligment.
 module PairAlign
@@ -47,7 +47,7 @@ module PairAlign
       @s_entry.seq.downcase!
 
       space = Space.new( 0, 0, @q_entry.length - 1, @s_entry.length - 1)
-      kmer  = 16
+      kmer  = 64
 
       align_recurse(@q_entry.seq, @s_entry.seq, space, kmer)
       matches_upcase
@@ -65,18 +65,19 @@ module PairAlign
     # cast and recursed into.
     def align_recurse(q_seq, s_seq, space, kmer, matches = [])
       matches = matches_select_by_space(matches, space)
+      matches = matches_select_by_score(matches, space)
 
       while (matches.size == 0 and kmer > 0)
         matches = Matches.find(q_seq, s_seq, space.q_min, space.s_min, space.q_max, space.s_max, kmer)
-        kmer /= 2
-      end
-
-      matches_score(matches, space)
 
-#      matches.each { |m| puts m.to_s(q_seq) }
+        if @matches.empty?
+          matches.sort_by! { |m| m.length }
+        else
+          matches = matches_select_by_score(matches, space)
+        end
 
-      unless @matches.empty?
-        matches = matches.select { |match| match.score > 0 }
+#      $stderr.puts "#{kmer}   #{matches.size}   #{space.q_dim} x #{space.s_dim}"
+        kmer /= 2
       end
 
       if best_match = matches.pop
@@ -102,12 +103,23 @@ module PairAlign
       new_matches
     end
 
+    # Method to select matches based on score.
+    def matches_select_by_score(matches, space)
+      matches_score(matches, space)
+
+      matches.select { |match| match.score > 0 }
+    end
+
     def matches_score(matches, space)
       matches.each do |match|
         score_length = match_score_length(match)
         score_diag   = match_score_diag(match, space)
 
         match.score = score_length + score_diag
+
+
+#        $stderr.puts "score_length: #{score_length}   score_diag: #{score_diag}   score: #{match.score}"   # DEBUG
+#        $stderr.puts match
       end
 
       matches.sort_by! { |match| match.score }
@@ -118,7 +130,7 @@ module PairAlign
     end
 
     def match_score_diag(match, space)
-      if space.q_dim >= space.s_dim   # s_dim is the narrow end
+      if space.q_dim > space.s_dim   # s_dim is the narrow end
         dist_beg = Math.dist_point2line(match.q_beg,
                                         match.s_beg,
                                         space.q_min,
@@ -321,7 +333,7 @@ module PairAlign
     end
 
     def empty?
-      if @q_max - @q_min > 0 and @s_max - @s_min > 0
+      if @q_max - @q_min >= 0 and @s_max - @s_min >= 0
         return false
       end
 
-- 
2.39.5