require 'open3'
require 'narray'
require 'maasha/fasta'
+require 'maasha/mum'
class AlignError < StandardError; end;
-CONSENSUS = 0.20
-FREQ_MIN = 2
-QUAL_MIN = 20
-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
# Class for aligning sequences.
class Align
- attr_reader :entries
+ attr_accessor :options
+ attr_reader :entries
# Class method to align sequences in a list of Seq objects and
# return these as a new list of Seq objects.
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
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 longest perfect 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)
+ q_space_beg = 0
+ q_space_end = q_entry.length
+ s_space_beg = 0
+ s_space_end = s_entry.length
+
+ q_entry.seq.downcase!
+ s_entry.seq.downcase!
+
+ ap = AlignPair.new(q_entry, s_entry)
+
+ ap.align(q_space_beg, q_space_end, s_space_beg, s_space_end, 15, [])
+
+ ap.upcase_match
+ ap.insert_gaps
+
+ #ap.dump_bp and exit
+
+ self.new([q_entry, s_entry])
+ end
# Method to initialize an Align object with a list of aligned Seq objects.
- def initialize(entries)
+ def initialize(entries, options = {})
@entries = entries
+ @options = options
end
# Method that returns the length of the alignment.
@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
- cons = Consensus.new(@entries)
- cons.to_seq
+ cons = Consensus.new(@entries, @options)
+ cons.consensus
end
# Method to pretty print an alignment from an Align object.
def to_s
- cons = self.consensus
- cons.mask_seq_soft!(QUAL_MIN) unless cons.qual.nil?
+ cons = Consensus.new(@entries, @options)
+ cons.mask_entries!
max_name = @entries.group_by { |entry| entry.seq_name.length }.max.first
+ output = ""
+
@entries.each do |entry|
- entry.mask_seq_soft!(QUAL_MIN) unless entry.qual.nil?
- puts entry.seq_name + (" " * (max_name + 3 - entry.seq_name.length )) + entry.seq
+ output << entry.seq_name + (" " * (max_name + 3 - entry.seq_name.length )) + entry.seq + $/
end
- output = ""
- output << " " * (max_name + 3) + cons.seq + $/
- output << " " * (max_name + 3) + cons.qual.tr("[@-h]", " ..........ooooooooooOOOOOOOOOO") unless cons.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 for creating a pairwise alignment of two specified sequences. The
+ # alignment is created by casting a search space the size of the sequences
+ # and save the longest perfect 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.
+ class AlignPair
+ attr_reader :matches
+
+ # Method to initialize an AlignPair object given two Seq objects.
+ def initialize(q_entry, s_entry)
+ @q_entry = q_entry
+ @s_entry = s_entry
+ @matches = []
+ end
+
+ # Method that locates the longest perfect match in a specified search space
+ # by comparing two sequence indeces. The first index is created by hashing
+ # for the first sequence the position of all non-overlapping unique oligos
+ # with the oligo as key an the position as value. The second index is
+ # created for the second sequence by hashing the position of all
+ # overlapping unique oligos in a similar way. The longest match is located
+ # by comparing these indeces and the longest match is expanded maximally
+ # within the specified search space. Finally, the longest match is used
+ # to define a left and a right search space which are recursed into
+ # until no more matches are found.
+ def align(q_space_beg, q_space_end, s_space_beg, s_space_end, kmer_size, matches)
+# puts "search space: #{q_space_beg}-#{q_space_end} x #{s_space_beg}-#{s_space_end}"
+ new_matches = []
+
+ matches.each do |match|
+ new_matches << match if match.q_beg >= q_space_beg and
+ match.q_end <= q_space_end and
+ match.s_beg >= s_space_beg and
+ match.s_end <= s_space_end
+ end
+
+ while new_matches.empty? and kmer_size > 0
+ if @q_entry.length >= kmer_size and @s_entry.length >= kmer_size
+ new_matches = MUMs.find(@q_entry, @s_entry, {"q_space_beg" => q_space_beg,
+ "q_space_end" => q_space_end,
+ "s_space_beg" => s_space_beg,
+ "s_space_end" => s_space_end,
+ "kmer_size" => kmer_size} )
+ end
+
+ unless new_matches.empty?
+ longest_match = new_matches.sort_by { |match| match.length }.pop
+# pp longest_match
+
+ @matches << longest_match
+
+ q_space_beg_left = (q_space_beg == 0) ? 0 : q_space_beg + 1
+ s_space_beg_left = (s_space_beg == 0) ? 0 : s_space_beg + 1
+ q_space_end_left = longest_match.q_beg - 1
+ s_space_end_left = longest_match.s_beg - 1
+
+ q_space_beg_right = longest_match.q_end + 1
+ s_space_beg_right = longest_match.s_end + 1
+ q_space_end_right = (q_space_end == @q_entry.length) ? @q_entry.length : q_space_end - 1
+ s_space_end_right = (s_space_end == @s_entry.length) ? @s_entry.length : s_space_end - 1
+
+ if q_space_end_left - q_space_beg_left > 0 and s_space_end_left - s_space_beg_left > 0
+ align(q_space_beg_left, q_space_end_left, s_space_beg_left, s_space_end_left, kmer_size, new_matches)
+ end
+
+ if q_space_end_right - q_space_beg_right > 0 and s_space_end_right - s_space_beg_right > 0
+ align(q_space_beg_right, q_space_end_right, s_space_beg_right, s_space_end_right, kmer_size, new_matches)
+ end
+ end
+
+ kmer_size /= 2
+ end
+ end
+
+ # Method for debugging purposes that upcase matching sequence while non-matches
+ # sequence is kept in lower case.
+ def upcase_match
+ @matches.each do |match|
+ @q_entry.seq[match.q_beg .. match.q_end] = @q_entry.seq[match.q_beg .. match.q_end].upcase
+ @s_entry.seq[match.s_beg .. match.s_end] = @s_entry.seq[match.s_beg .. match.s_end].upcase
+ end
+ end
+
+ # Method that insert gaps in sequences based on a list of matches and thus
+ # creating an alignment.
+ # TODO check boundaries!
+ def insert_gaps
+ @matches.sort_by! { |m| m.q_beg }
+
+ q_gaps = 0
+ s_gaps = 0
+
+ @matches.each do |match|
+ diff = (q_gaps + match.q_beg) - (s_gaps + match.s_beg)
+
+ #pp "q_gaps #{q_gaps} s_gaps #{s_gaps} diff #{diff}"
+
+ if diff < 0
+ @q_entry.seq.insert(match.q_beg + q_gaps, "-" * diff.abs)
+ q_gaps += diff.abs
+ elsif diff > 0
+ @s_entry.seq.insert(match.s_beg + s_gaps, "-" * diff.abs)
+ s_gaps += diff.abs
+ end
+ end
+
+ diff = @q_entry.length - @s_entry.length
+
+ if diff < 0
+ @q_entry.seq << ("-" * diff.abs)
+ else
+ @s_entry.seq << ("-" * diff.abs)
+ end
+ end
+
+ # Method that dumps all matches as Biopiece records for use with
+ # plot_matches.
+ def dump_bp
+ @matches.sort_by! { |m| m.q_beg }
+
+ @matches.each { |m|
+ puts "Q_BEG: #{m.q_beg}"
+ puts "Q_END: #{m.q_end}"
+ puts "S_BEG: #{m.s_beg}"
+ puts "S_END: #{m.s_end}"
+ puts "STRAND: +"
+ puts "---"
+ }
+ end
+ end
+
class Consensus
# Method to initialize a Consensus object given a list of aligned Seq object.
- def initialize(entries)
+ def initialize(entries, options)
@entries = entries
- @cols = entries.first.length
- @rows = entries.size
+ @options = options
+
+ @cols = entries.first.seq.length
+ @rows = entries.size
+
+ @has_qual = entries.first.qual.nil? ? false : true
+
@na_seq = NArray.byte(@cols, @rows)
- @na_qual = 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)
- m1 = mask_high_qual
- m2 = mask_conserved_columns
- m3 = mask_conserved_residues
+ @entries.each_with_index do |entry, i|
+ na_seq[true, i] = NArray.to_na(entry.seq.upcase, "byte")
+ end
- @na_seq *= m3
- @na_qual *= m3
+ na_seq += ((na_seq.ne('-'.ord) - @na_seq.ne(0)) * ' '.ord)
- @na_seq *= (m1 | m2)
- @na_qual *= (m1 | m2)
+ @entries.each_with_index do |entry, i|
+ entry.seq = na_seq[true, i].to_s
+ end
end
- def to_seq
+ # 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
+ new_seq.qual = consensus_qual if @has_qual
new_seq.type = "dna"
new_seq
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") - SCORE_ILLUMINA
+ @na_qual[true, i] = NArray.to_na(entry.qual, "byte") - SCORE_BASE if @has_qual
end
end
- def mask_high_qual
- @na_qual > QUAL_MIN
+ # 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
+ mask = mask_quality_min
+
+ @na_seq *= mask
+ @na_qual *= mask
+ end
+
+ if @options.has_key? :quality_mean
+ mask = mask_quality_mean
+
+ @na_seq *= mask
+ @na_qual *= mask
+ end
+ end
+
+ if @options.has_key? :sequence_min
+ mask = mask_sequence_min
+
+ @na_seq *= mask
+ @na_qual *= mask if @has_qual
+ end
+
+ if @options.has_key? :gap_max
+ mask = mask_gap_max
+
+ @na_seq *= mask
+ @na_qual *= mask if @has_qual
+ end
+
+ if @options.has_key? :residue_min
+ mask = mask_residue_min
+
+ @na_seq *= mask
+ @na_qual *= mask if @has_qual
+ end
end
- def mask_conserved_columns
- @na_seq * (@na_seq - @na_seq[true, 0]).sum(1).eq(0) > 0
+ # 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).to_type("int").sum(1) >= @options[:sequence_min])
+ mask
end
- def mask_conserved_residues
- factor = (1 / @rows.to_f)
- 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 that indicates which residue frequencies that are above the residue_min.
+ # The residue frequencies are calculated for each column and residue type as the
+ # number of each residue type over the sum of all non-gap residues in that column.
+ # 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).to_type("float").sum(1)
+
+ 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)).to_f * factor > CONSENSUS
- mask_T = (mask_T * mask_T.sum(1)).to_f * factor > CONSENSUS
- mask_C = (mask_C * mask_C.sum(1)).to_f * factor > CONSENSUS
- mask_G = (mask_G * mask_G.sum(1)).to_f * factor > CONSENSUS
+ mask_A = (mask_A * mask_A.sum(1)) * factor >= cons_min
+ mask_T = (mask_T * mask_T.sum(1)) * factor >= cons_min
+ mask_C = (mask_C * mask_C.sum(1)) * factor >= cons_min
+ mask_G = (mask_G * mask_G.sum(1)) * factor >= cons_min
mask_A | mask_T | mask_C | mask_G
end
+ # Mask that indicates which columns contain fewer gaps than max_gap.
+ # 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).to_type("float").sum(1) / @rows > @options[:gap_max]
+
+ mask
+ end
+
+ # Mask that indicates which residues in an alignment are above quality_min.
+ # Positions with subquality are 0'ed - all others are set to 1.
+ def mask_quality_min
+ @na_qual > @options[:quality_min]
+ end
+
+ # Mask that indicates which columns have a quality mean above quality_mean which
+ # is the mean of all non-gap quality residues in that column. Columns with less then
+ # quality_mean are 0'ed, otherwise set to 1.
+ def mask_quality_mean
+ mask = NArray.byte(@cols, @rows) + 1
+ residues = @na_seq.ne(0).to_type("int").sum(1)
+ quality = @na_qual.to_type("float").sum(1)
+
+ mask * (quality / residues).round > @options[:quality_mean]
+ end
+
# Method to calculate a consensus sequence from a Consensus object.
def consensus_seq
cons = NArray.byte(@cols)
# 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") + SCORE_BASE).to_s
end
end
end
-
-
__END__
-
-cons |= ((@na_seq & BIT_A > 0).sum(1).to_type("float") * (1 / @rows.to_f) > CONSENSUS) * BIT_A
-cons |= ((@na_seq & BIT_T > 0).sum(1).to_type("float") * (1 / @rows.to_f) > CONSENSUS) * BIT_T
-cons |= ((@na_seq & BIT_C > 0).sum(1).to_type("float") * (1 / @rows.to_f) > CONSENSUS) * BIT_C
-cons |= ((@na_seq & BIT_G > 0).sum(1).to_type("float") * (1 / @rows.to_f) > CONSENSUS) * BIT_G
projects / biopieces.git / blobdiff