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
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)
Open3.popen3("muscle", "-quiet") do |stdin, stdout, stderr|
entries.each do |entry|
- raise AlignError, "Duplicate sequence name: #{entry.seq_name}" if index.has_key? entry.seq_name
+ raise AlignError, "Duplicate sequence name: #{entry.seq_name}" if index[entry.seq_name]
- index[entry.seq_name] = entry
+ index[entry.seq_name] = entry.dup
stdin.puts entry.to_fasta
end
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 = {})
@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
# Method to pretty print an alignment from an Align object.
def to_s
cons = Consensus.new(@entries, @options)
+ cons.mask_entries!
- entries = cons.each
- cons_entry = cons.consensus
-
- max_name = entries.group_by { |entry| entry.seq_name.length }.max.first
+ max_name = @entries.group_by { |entry| entry.seq_name.length }.max.first
output = ""
- entries.each do |entry|
+ @entries.each do |entry|
output << entry.seq_name + (" " * (max_name + 3 - entry.seq_name.length )) + entry.seq + $/
end
- output << " " * (max_name + 3) + cons_entry.seq + $/
- output << " " * (max_name + 3) + cons_entry.qual.tr("[@-h]", " ..........ooooooooooOOOOOOOOOO") unless cons_entry.qual.nil?
+ cons_entry = cons.consensus
+
+ output << " " * (max_name + 3) + cons_entry.seq
+ output << $/ + " " * (max_name + 3) + cons_entry.qual.tr("[@-h]", " ..........ooooooooooOOOOOOOOOO") unless cons_entry.qual.nil?
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
@has_qual = entries.first.qual.nil? ? false : true
- @na_seq = NArray.byte(entries.first.length, entries.size)
- @na_qual = NArray.byte(entries.first.length, entries.size) if @has_qual
+ @na_seq = NArray.byte(@cols, @rows)
+ @na_qual = NArray.byte(@cols, @rows) if @has_qual
na_add_entries
+ consensus_calc
+ end
+
+ # Method that lowercase residues that have been removed during
+ # the determination of the consensus sequence.
+ def mask_entries!
+ na_seq = NArray.byte(@cols, @rows)
+
+ @entries.each_with_index do |entry, i|
+ na_seq[true, i] = NArray.to_na(entry.seq.upcase, "byte")
+ end
+
+ na_seq += ((na_seq.ne('-'.ord) - @na_seq.ne(0)) * ' '.ord)
+
+ @entries.each_with_index do |entry, i|
+ entry.seq = na_seq[true, i].to_s
+ end
+ end
+
+ # Method that returns a Sequence object with a consensus sequence
+ # for the entries in an Align object.
+ def consensus
+ new_seq = Seq.new
+ new_seq.seq = consensus_seq
+ new_seq.qual = consensus_qual if @has_qual
+ new_seq.type = "dna"
+
+ new_seq
+ end
+ private
+
+ # Method to add a Seq entry object to the two NArrays; @na_seq and @na_qual
+ 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") - Seq::SCORE_BASE if @has_qual
+ end
+ end
+
+ # Method to calculate a consensus sequence from a list of sequenced stored in two
+ # NArrays.
+ def consensus_calc
if @has_qual
- if @options.has_key? :quality_min
+ if @options[:quality_min]
mask = mask_quality_min
@na_seq *= mask
@na_qual *= mask
end
- if @options.has_key? :quality_mean
+ if @options[:quality_mean]
mask = mask_quality_mean
@na_seq *= mask
end
end
- if @options.has_key? :sequence_min
+ if @options[:sequence_min]
mask = mask_sequence_min
@na_seq *= mask
@na_qual *= mask if @has_qual
end
- if @options.has_key? :gap_max
+ if @options[:gap_max]
mask = mask_gap_max
@na_seq *= mask
@na_qual *= mask if @has_qual
end
- if @options.has_key? :residue_min
+ if @options[:residue_min]
mask = mask_residue_min
@na_seq *= mask
end
end
- # Method for iterating over the rows of an Align object and
- # return Seq objects for each sequence.
- def each
- entries = []
-
- (0 ... @rows).each do |i|
- entry = Seq.new
- entry.seq_name = @entries[i].seq_name
- entry.seq = @na_seq[true, i].to_s.tr!(TR_HEX, TR_NUC).upcase
- entry.qual = (@na_qual[true, i].to_type("byte") + SCORE_ILLUMINA).to_s if @has_qual
- entry.type = @entries[i].type
-
- if block_given?
- yield entry
- else
- entries << entry
- end
- end
-
- return entries unless block_given?
- end
-
- # Method that returns a Sequence object with a consensus sequence
- # for the entries in an Align object.
- def consensus
- new_seq = Seq.new
- new_seq.seq = consensus_seq
- new_seq.qual = consensus_qual if @has_qual
- new_seq.type = "dna"
-
- new_seq
- end
-
- private
-
- 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_ILLUMINA if @has_qual
- end
- end
-
# Mask that indicates which columns have more than sequence_min sequences.
# Columns with less than sequence_min are 0'ed, otherwise set to 1.
def mask_sequence_min
mask = NArray.byte(@cols, @rows) + 1
- mask *= ((@na_seq > 0).sum(1) >= @options[:sequence_min])
+ mask *= ((@na_seq > 0).to_type("int").sum(1) >= @options[:sequence_min])
mask
end
# Positions with residues above the residue_min are indicated with 1.
def mask_residue_min
cons_min = @options[:residue_min]
- factor = 1 / @na_seq.ne(0).sum(1).to_type("float")
+ factor = 1 / @na_seq.ne(0).to_type("float").sum(1)
- mask_A = @na_seq & BIT_A > 0
- mask_T = @na_seq & BIT_T > 0
- mask_C = @na_seq & BIT_C > 0
- mask_G = @na_seq & BIT_G > 0
+ mask_A = (@na_seq & BIT_A > 0).to_type("int")
+ mask_T = (@na_seq & BIT_T > 0).to_type("int")
+ mask_C = (@na_seq & BIT_C > 0).to_type("int")
+ mask_G = (@na_seq & BIT_G > 0).to_type("int")
mask_A = (mask_A * mask_A.sum(1)) * factor >= cons_min
mask_T = (mask_T * mask_T.sum(1)) * factor >= cons_min
# Columns with more gaps are 0'ed, others are set to 1.
def mask_gap_max
mask = NArray.byte(@cols, @rows) + 1
- mask *= @na_seq.ne(0).sum(1).to_type("float") / @rows > @options[:gap_max]
+ mask *= @na_seq.ne(0).to_type("float").sum(1) / @rows > @options[:gap_max]
mask
end
# Method to calculate a consensus quality from a Consensus object.
def consensus_qual
- (@na_qual.mean(1).round.to_type("byte") + SCORE_ILLUMINA).to_s
+ (@na_qual.mean(1).round.to_type("byte") + Seq::SCORE_BASE).to_s
end
end
end
projects / biopieces.git / blobdiff