changed upper/lower case output in assemble_pairs

[biopieces.git] / code_ruby / lib / maasha / seq / assemble.rb

diff --git a/code_ruby/lib/maasha/seq/assemble.rb b/code_ruby/lib/maasha/seq/assemble.rb
index e59e178b73df8ba45689836db02c852d5632f637..016322727cf15eb824f9e3f0518cca9b2ecb0241 100644 (file)
--- a/code_ruby/lib/maasha/seq/assemble.rb
+++ b/code_ruby/lib/maasha/seq/assemble.rb
@@ -22,12 +22,15 @@
 
 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
+# Class containing methods to assemble two overlapping sequences into a single.
 class Assemble
+  # Class method to assemble two Seq objects.
   def self.pair(entry1, entry2, options = {})
     assemble = self.new(entry1, entry2, options)
     assemble.match
   end
 
+  # Method to initialize an Assembly object.
   def initialize(entry1, entry2, options)
     @entry1  = entry1
     @entry2  = entry2
@@ -38,6 +41,7 @@ class Assemble
     @options[:overlap_max]      = [@options[:overlap_max], entry1.length, entry2.length].min
   end
 
+  # Method to locate overlapping matche between two sequences.
   def match
     overlap = @options[:overlap_max]
 
@@ -47,31 +51,52 @@ class Assemble
     while overlap >= @options[:overlap_min]
       hamming_dist = (na_seq1[-1 * overlap .. -1] ^ na_seq2[0 ... overlap]).count_true
 
-      if hamming_dist <= percent2real(overlap, @options[:mismatches_max])
-        merged = @entry1 + @entry2[overlap .. -1]
+      if hamming_dist <= (overlap * @options[:mismatches_max] * 0.01).round
+        entry_left  = @entry1[0 ... @entry1.length - overlap]
+        entry_right = @entry2[overlap .. -1]
 
         if @entry1.qual and @entry2.qual
-          qual1 = @entry1.qual[@entry1.length - overlap .. -1]
-          qual2 = @entry2.qual[0 ... overlap]
+          entry_overlap1 = @entry1[-1 * overlap .. -1]
+          entry_overlap2 = @entry2[0 ... overlap]
 
-          na_qual1 = NArray.to_na(qual1, "byte")
-          na_qual2 = NArray.to_na(qual2, "byte")
-
-          qual = ((na_qual1 + na_qual2) / 2).to_s
-
-          merged.seq_name = @entry1.seq_name + ":overlap=#{overlap}:hamming=#{hamming_dist}"
-          merged.qual[@entry1.length - overlap ... @entry1.length]  = qual
+          entry_overlap = merge_overlap(entry_overlap1, entry_overlap2)
+        else
+          entry_overlap = @entry1[-1 * overlap .. -1]
         end
 
-        return merged
+        entry_left.seq.downcase!
+        entry_overlap.seq.upcase!
+        entry_right.seq.downcase!
+        entry_merged          = entry_left + entry_overlap + entry_right
+        entry_merged.seq_name = @entry1.seq_name + ":overlap=#{overlap}:hamming=#{hamming_dist}"
+
+        return entry_merged
       end
 
       overlap -= 1
     end
   end
 
-  def percent2real(length, percent)
-      (length * percent * 0.01).round
+  # Method to merge sequence and quality scores in an overlap.
+  # The residue with the highest score at mismatch positions is selected.
+  # The quality scores of the overlap are the mean of the two sequences.
+  def merge_overlap(entry_overlap1, entry_overlap2)
+    na_seq = NArray.byte(entry_overlap1.length, 2)
+    na_seq[true, 0] = NArray.to_na(entry_overlap1.seq.downcase, "byte")
+    na_seq[true, 1] = NArray.to_na(entry_overlap2.seq.downcase, "byte")
+
+    na_qual = NArray.byte(entry_overlap1.length, 2)
+    na_qual[true, 0] = NArray.to_na(entry_overlap1.qual, "byte")
+    na_qual[true, 1] = NArray.to_na(entry_overlap2.qual, "byte")
+
+    mask_xor = na_seq[true, 0] ^ na_seq[true, 1] > 0
+    mask_seq = ((na_qual * mask_xor).eq( (na_qual * mask_xor).max(1)))
+
+    merged      = Seq.new()
+    merged.seq  = (na_seq * mask_seq).max(1).to_s
+    merged.qual = na_qual.mean(1).round.to_type("byte").to_s
+
+    merged
   end
 end