added Seq.shuffle method and tests

[biopieces.git] / code_ruby / lib / maasha / align.rb

diff --git a/code_ruby/lib/maasha/align.rb b/code_ruby/lib/maasha/align.rb
index 32532fe26be362277eb36670ffa400aea7df5066..ffa41c4e41e8107d4fe828ef6b81c529b320741f 100755 (executable)
--- a/code_ruby/lib/maasha/align.rb
+++ b/code_ruby/lib/maasha/align.rb
@@ -25,18 +25,19 @@
 require 'pp'
 require 'open3'
 require 'narray'
+require 'maasha/align/pair'
 require 'maasha/fasta'
 
 class AlignError < StandardError; end;
 
-ALPH_DNA    = %w{A T C G}
-ALPH_AMBI   = %w{A T C G M R W S Y K V H D B N}
+ALPH_DNA  = %w{A T C G}
+ALPH_AMBI = %w{A T C G M R W S Y K V H D B N}
 
 BIT_INDEL = 0
-# BIT_A = 1 << 0
-# BIT_T = 1 << 1
-# BIT_C = 1 << 2
-# BIT_G = 1 << 3
+BIT_A = 1 << 0
+BIT_T = 1 << 1
+BIT_C = 1 << 2
+BIT_G = 1 << 3
 
 BIT_M = BIT_A | BIT_C
 BIT_R = BIT_A | BIT_G
@@ -91,6 +92,8 @@ class Align
   attr_accessor :options
   attr_reader   :entries
 
+  include PairAlign
+
   # Class method to align sequences in a list of Seq objects and
   # return these as a new list of Seq objects.
   def self.muscle(entries, verbose = false)
@@ -101,7 +104,7 @@ class Align
       entries.each do |entry|
         raise AlignError, "Duplicate sequence name: #{entry.seq_name}" if index.has_key? entry.seq_name
 
-        index[entry.seq_name] = entry
+        index[entry.seq_name] = entry.dup
 
         stdin.puts entry.to_fasta
       end
@@ -126,6 +129,17 @@ class Align
 
     self.new(result)
   end
+  
+  # Class method to create a pairwise alignment of two given Seq objects. The
+  # alignment is created by casting a search space the size of the sequences
+  # and save the best scoring match between the sequences and recurse into
+  # the left and right search spaced around this match. When all search spaces
+  # are exhausted the saved matches are used to insert gaps in the sequences.
+  def self.pair(q_entry, s_entry)
+    AlignPair.align(q_entry, s_entry)
+
+    self.new([q_entry, s_entry])
+  end
 
   # Method to initialize an Align object with a list of aligned Seq objects.
   def initialize(entries, options = {})
@@ -143,6 +157,22 @@ class Align
     @entries.size
   end
 
+  # Method that returns the identity of an alignment with two members.
+  def identity
+    if self.members != 2
+      raise AlignError "Bad number of members for similarity calculation: #{self.members}"
+    end
+
+    na1 = NArray.to_na(@entries[0].seq.upcase, "byte")
+    na2 = NArray.to_na(@entries[1].seq.upcase, "byte")
+
+    shared   = (na1 - na2).count_false
+    total    = (@entries[0].length < @entries[1].length) ? @entries[0].length : @entries[1].length
+    identity = shared.to_f / total
+
+    identity
+  end
+
   # Method to create a consensus sequence from an Align object and
   # return a new Seq object with the consensus.
   def consensus
@@ -170,6 +200,15 @@ class Align
     output
   end
 
+  # Method for iterating each of the aligned sequences.
+  def each
+    if block_given?
+      @entries.each { |entry| yield entry }
+    else
+      return @entries
+    end
+  end
+
   private
 
   class Consensus
@@ -223,7 +262,7 @@ class Align
     def na_add_entries
       @entries.each_with_index do |entry, i|
         @na_seq[true, i]  = NArray.to_na(entry.seq.downcase.tr(TR_NUC, TR_HEX), "byte")
-        @na_qual[true, i] = NArray.to_na(entry.qual, "byte") - SCORE_BASE if @has_qual
+        @na_qual[true, i] = NArray.to_na(entry.qual, "byte") - Seq::SCORE_BASE if @has_qual
       end
     end
 
@@ -336,7 +375,7 @@ class Align
 
     # Method to calculate a consensus quality from a Consensus object.
     def consensus_qual
-      (@na_qual.mean(1).round.to_type("byte") + SCORE_BASE).to_s
+      (@na_qual.mean(1).round.to_type("byte") + Seq::SCORE_BASE).to_s
     end
   end
 end