# 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
oligo = seq[pos ... pos + kmer]
index_hash[oligo] << pos
- pos += 1
+ pos += step
end
index_hash
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
@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
# 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
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 }
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,
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