X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=code_ruby%2Flib%2Fmaasha%2Fseq.rb;h=7377da533d50ea34fb4dc12c3e2bfb7f05a7e944;hb=d71b5469d06eeaa8de7d83b501eb17a0dedd7999;hp=8833acba2f54ca3eb22407912f8421d324e6ff52;hpb=c54a2cc53ee2a7db8028ef593be593e527dc50c2;p=biopieces.git diff --git a/code_ruby/lib/maasha/seq.rb b/code_ruby/lib/maasha/seq.rb index 8833acb..7377da5 100644 --- a/code_ruby/lib/maasha/seq.rb +++ b/code_ruby/lib/maasha/seq.rb @@ -27,8 +27,12 @@ require 'maasha/seq/digest' require 'maasha/seq/trim' require 'narray' -autoload :BackTrack, 'maasha/seq/backtrack.rb' -autoload :Dynamic, 'maasha/seq/dynamic.rb' +autoload :BackTrack, 'maasha/seq/backtrack' +autoload :Dynamic, 'maasha/seq/dynamic' +autoload :Homopolymer, 'maasha/seq/homopolymer' +autoload :Hamming, 'maasha/seq/hamming' +autoload :Levenshtein, 'maasha/seq/levenshtein' +autoload :Ambiguity, 'maasha/seq/ambiguity' # Residue alphabets DNA = %w[a t c g] @@ -67,15 +71,15 @@ TRANS_TAB11 = { "GTG" => "V", "GCG" => "A", "GAG" => "E", "GGG" => "G" } -# Quality scores bases -SCORE_BASE = 64 -SCORE_MIN = 0 -SCORE_MAX = 40 - # Error class for all exceptions to do with Seq. class SeqError < StandardError; end class Seq + # Quality scores bases + SCORE_BASE = 33 + SCORE_MIN = 0 + SCORE_MAX = 40 + include Digest include Trim @@ -86,10 +90,10 @@ class Seq def self.new_bp(record) seq_name = record[:SEQ_NAME] seq = record[:SEQ] - type = record[:SEQ_TYPE] + type = record[:SEQ_TYPE].to_sym if record[:SEQ_TYPE] qual = record[:SCORES] - self.new(seq_name, seq, type, qual) + self.new(seq_name: seq_name, seq: seq, type: type, qual: qual) end # Class method that generates all possible oligos of a specifed length and type. @@ -97,9 +101,9 @@ class Seq raise SeqError, "Cannot generate negative oligo length: #{length}" if length <= 0 case type.downcase - when /dna/ then alph = DNA - when /rna/ then alph = RNA - when /protein/ then alph = PROTEIN + when :dna then alph = DNA + when :rna then alph = RNA + when :protein then alph = PROTEIN else raise SeqError, "Unknown sequence type: #{type}" end @@ -121,16 +125,20 @@ class Seq oligos end - # Initialize a sequence object with the following arguments: - # - seq_name: Name of the sequence. - # - seq: The sequence. - # - type: The sequence type - DNA, RNA, or protein - # - qual: An Illumina type quality scores string. - def initialize(seq_name = nil, seq = nil, type = nil, qual = nil) - @seq_name = seq_name - @seq = seq - @type = type - @qual = qual + # Initialize a sequence object with the following options: + # - :seq_name Name of the sequence. + # - :seq The sequence. + # - :type The sequence type - DNA, RNA, or protein + # - :qual An Illumina type quality scores string. + def initialize(options = {}) + @seq_name = options[:seq_name] + @seq = options[:seq] + @type = options[:type] + @qual = options[:qual] + + if @seq and @qual and @seq.length != @qual.length + raise SeqError, "Sequence length and score length mismatch: #{@seq.length} != #{@qual.length}" + end end # Method that guesses and returns the sequence type @@ -139,9 +147,9 @@ class Seq raise SeqError, "Guess failed: sequence is nil" if self.seq.nil? case self.seq[0 ... 100].downcase - when /[flpqie]/ then return "protein" - when /[u]/ then return "rna" - else return "dna" + when /[flpqie]/ then return :protein + when /[u]/ then return :rna + else return :dna end end @@ -149,6 +157,7 @@ class Seq # by inspecting the first 100 residues. def type_guess! self.type = self.type_guess + self end # Returns the length of a sequence. @@ -188,24 +197,24 @@ class Seq # Method that returns true is a given sequence type is DNA. def is_dna? - self.type == 'dna' + self.type == :dna end # Method that returns true is a given sequence type is RNA. def is_rna? - self.type == 'rna' + self.type == :rna end # Method that returns true is a given sequence type is protein. def is_protein? - self.type == 'protein' + self.type == :protein end # Method to transcribe DNA to RNA. def to_rna raise SeqError, "Cannot transcribe 0 length sequence" if self.length == 0 raise SeqError, "Cannot transcribe sequence type: #{self.type}" unless self.is_dna? - self.type = 'rna' + self.type = :rna self.seq.tr!('Tt','Uu') end @@ -213,14 +222,13 @@ class Seq def to_dna raise SeqError, "Cannot reverse-transcribe 0 length sequence" if self.length == 0 raise SeqError, "Cannot reverse-transcribe sequence type: #{self.type}" unless self.is_rna? - - self.type = 'dna' + self.type = :dna self.seq.tr!('Uu','Tt') end # Method to translate a DNA sequence to protein. def translate!(trans_tab = 11) - raise SeqError, "Sequence type must be 'dna' - not #{self.type}" unless self.type == 'dna' + raise SeqError, "Sequence type must be 'dna' - not #{self.type}" unless self.type == :dna raise SeqError, "Sequence length must be a multiplum of 3 - was: #{self.length}" unless (self.length % 3) == 0 case trans_tab @@ -255,7 +263,7 @@ class Seq self.seq = protein self.qual = nil - self.type = "protein" + self.type = :protein self end @@ -286,8 +294,8 @@ class Seq raise SeqError, "Missing seq_name" if self.seq_name.nil? or self.seq_name == '' raise SeqError, "Missing seq" if self.seq.nil? or self.seq.empty? - seq_name = self.seq_name.to_s - seq = self.seq.to_s + seq_name = self.seq_name + seq = self.seq.dup unless wrap.nil? seq.gsub!(/(.{#{wrap}})/) do |match| @@ -297,7 +305,7 @@ class Seq seq.chomp! end - ">" + seq_name + $/ + seq + $/ + ">" + seq_name.to_s + $/ + seq + $/ end # Method that given a Seq entry returns a FASTQ entry (a string). @@ -333,17 +341,20 @@ class Seq key end - # Method to reverse complement sequence. - def reverse_complement - self.reverse - self.complement - self - end + # Method to reverse the sequence. + def reverse + entry = Seq.new( + seq_name: self.seq_name, + seq: self.seq.reverse, + type: self.type, + qual: (self.qual ? self.qual.reverse : self.qual) + ) - alias :revcomp :reverse_complement + entry + end # Method to reverse the sequence. - def reverse + def reverse! self.seq.reverse! self.qual.reverse! if self.qual self @@ -353,6 +364,27 @@ class Seq def complement raise SeqError, "Cannot complement 0 length sequence" if self.length == 0 + entry = Seq.new( + seq_name: self.seq_name, + type: self.type, + qual: self.qual + ) + + if self.is_dna? + entry.seq = self.seq.tr('AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'TCGAAYRWSKMDHBVNtcgaayrwskmdhbvn') + elsif self.is_rna? + entry.seq = self.seq.tr('AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'UCGAAYRWSKMDHBVNucgaayrwskmdhbvn') + else + raise SeqError, "Cannot complement sequence type: #{self.type}" + end + + entry + end + + # Method that complements sequence including ambiguity codes. + def complement! + raise SeqError, "Cannot complement 0 length sequence" if self.length == 0 + if self.is_dna? self.seq.tr!('AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'TCGAAYRWSKMDHBVNtcgaayrwskmdhbvn') elsif self.is_rna? @@ -360,70 +392,95 @@ class Seq else raise SeqError, "Cannot complement sequence type: #{self.type}" end + + self end # Method to determine the Hamming Distance between # two Sequence objects (case insensitive). - def hamming_distance(seq) - self.seq.upcase.hamming_distance(seq.seq.upcase) + def hamming_distance(entry, options = nil) + if options and options[:ambiguity] + Hamming.distance(self.seq, entry.seq) + else + self.seq.upcase.hamming_distance(entry.seq.upcase) + end + end + + # Method to determine the Edit Distance between + # two Sequence objects (case insensitive). + def edit_distance(entry) + Levenshtein.distance(self.seq, entry.seq) end # Method that generates a random sequence of a given length and type. def generate(length, type) raise SeqError, "Cannot generate sequence length < 1: #{length}" if length <= 0 - case type.downcase - when "dna" - alph = DNA - when "rna" - alph = RNA - when "protein" - alph = PROTEIN + case type + when :dna then alph = DNA + when :rna then alph = RNA + when :protein then alph = PROTEIN else raise SeqError, "Unknown sequence type: #{type}" end seq_new = Array.new(length) { alph[rand(alph.size)] }.join("") self.seq = seq_new - self.type = type.downcase + self.type = type seq_new end - # Method to shuffle a sequence readomly inline. + # Method to return a new Seq object with shuffled sequence. + def shuffle + Seq.new( + seq_name: self.seq_name, + seq: self.seq.split('').shuffle!.join, + type: self.type, + qual: self.qual + ) + end + + # Method to shuffle a sequence randomly inline. def shuffle! self.seq = self.seq.split('').shuffle!.join self end - # Method that returns a subsequence of from a given start position - # and of a given length. - def subseq(start, length = self.length - start) - raise SeqError, "subsequence start: #{start} < 0" if start < 0 - raise SeqError, "subsequence length: #{length} < 0" if length < 0 - raise SeqError, "subsequence start + length > Seq.length: #{start} + #{length} > #{self.length}" if start + length > self.length + # Method to add two Seq objects. + def +(entry) + new_entry = Seq.new() + new_entry.seq = self.seq + entry.seq + new_entry.type = self.type if self.type == entry.type + new_entry.qual = self.qual + entry.qual if self.qual and entry.qual + new_entry + end - if length == 0 - seq = "" - qual = "" unless self.qual.nil? - else - stop = start + length - 1 + # Method to concatenate sequence entries. + def <<(entry) + raise SeqError, "sequences of different types" unless self.type == entry.type + raise SeqError, "qual is missing in one entry" unless self.qual.class == entry.qual.class - seq = self.seq[start .. stop] - qual = self.qual[start .. stop] unless self.qual.nil? - end + self.seq << entry.seq + self.qual << entry.qual unless entry.qual.nil? - Seq.new(self.seq_name, seq, self.type, qual) # TODO changed self.seq_name.dup to self.seq_name -> consequence? + self end - # Method that replaces a sequence with a subsequence from a given start position - # and of a given length. - def subseq!(start, length = self.length - start) - s = subseq(start, length) + # Index method for Seq objects. + def [](*args) + entry = Seq.new + entry.seq_name = self.seq_name.dup if self.seq_name + entry.seq = self.seq[*args] + entry.type = self.type + entry.qual = self.qual[*args] unless self.qual.nil? - self.seq_name = s.seq_name - self.seq = s.seq - self.type = s.type - self.qual = s.qual + entry + end + + # Index assignment method for Seq objects. + def []=(*args, entry) + self.seq[*args] = entry.seq[*args] + self.qual[*args] = entry.qual[*args] unless self.qual.nil? self end @@ -437,7 +494,7 @@ class Seq start = rand(self.length - length + 1) end - self.subseq(start, length) + self[start .. start + length] end # Method that returns the residue compositions of a sequence in @@ -453,23 +510,6 @@ class Seq comp end - # Method that returns the length of the longest homopolymeric stretch - # found in a sequence. - def homopol_max(min = 1) - return 0 if self.seq.nil? or self.seq.empty? - - found = false - - self.seq.upcase.scan(/A{#{min},}|T{#{min},}|G{#{min},}|C{#{min},}|N{#{min},}/) do |match| - found = true - min = match.size > min ? match.size : min - end - - return 0 unless found - - min - end - # Method that returns the percentage of hard masked residues # or N's in a sequence. def hard_mask @@ -482,21 +522,6 @@ class Seq ((self.seq.scan(/[a-z]/).size.to_f / (self.len - self.indels).to_f) * 100).round(2) end - # Hard masks sequence residues where the corresponding quality score - # is below a given cutoff. - def mask_seq_hard_old(cutoff) - seq = self.seq.upcase - scores = self.qual - i = 0 - - scores.each_char do |score| - seq[i] = 'N' if score.ord - SCORE_BASE < cutoff - i += 1 - end - - self.seq = seq - end - # Hard masks sequence residues where the corresponding quality score # is below a given cutoff. def mask_seq_hard!(cutoff) @@ -540,61 +565,51 @@ class Seq def qual_base33? self.qual.match(/[!-:]/) ? true : false end - - # Method that determines if a quality score string can be - # absolutely identified as base 64. + + # Method that determines if a quality score string may be base 64. def qual_base64? self.qual.match(/[K-h]/) ? true : false end - # Method to determine if a quality score is valid. + # Method to determine if a quality score is valid accepting only 0-40 range. def qual_valid?(encoding) raise SeqError, "Missing qual" if self.qual.nil? - case encoding.downcase - when "sanger" then return true if self.qual.match(/^[!-~]*$/) - when "454" then return true if self.qual.match(/^[@-~]*$/) - when "solexa" then return true if self.qual.match(/^[;-~]*$/) - when "illumina13" then return true if self.qual.match(/^[@-~]*$/) - when "illumina15" then return true if self.qual.match(/^[@-~]*$/) - when "illumina18" then return true if self.qual.match(/^[!-~]*$/) + case encoding + when :base_33 then return true if self.qual.match(/^[!-I]*$/) + when :base_64 then return true if self.qual.match(/^[@-h]*$/) else raise SeqError, "unknown quality score encoding: #{encoding}" end false end - # Method to convert quality scores inbetween formats. - # Sanger base 33, range 0-40 - # 454 base 64, range 0-40 - # Solexa base 64, range -5-40 - # Illumina13 base 64, range 0-40 - # Illumina15 base 64, range 0-40 - # Illumina18 base 33, range 0-41 - def convert_scores!(from, to) - unless from == to - na_qual = NArray.to_na(self.qual, "byte") + # Method to coerce quality scores to be within the 0-40 range. + def qual_coerce!(encoding) + raise SeqError, "Missing qual" if self.qual.nil? - case from.downcase - when "sanger" then na_qual -= 33 - when "454" then na_qual -= 64 - when "solexa" then na_qual -= 64 - when "illumina13" then na_qual -= 64 - when "illumina15" then na_qual -= 64 - when "illumina18" then na_qual -= 33 - else raise SeqError, "unknown quality score encoding: #{from}" - end + case encoding + when :base_33 then self.qual.tr!("[J-~]", "I") + when :base_64 then self.qual.tr!("[i-~]", "h") + else raise SeqError, "unknown quality score encoding: #{encoding}" + end - case to.downcase - when "sanger" then na_qual += 33 - when "454" then na_qual += 64 - when "solexa" then na_qual += 64 - when "illumina13" then na_qual += 64 - when "illumina15" then na_qual += 64 - when "illumina18" then na_qual += 33 - else raise SeqError, "unknown quality score encoding: #{to}" - end + self + end + + # Method to convert quality scores. + def qual_convert!(from, to) + raise SeqError, "unknown quality score encoding: #{from}" unless from == :base_33 or from == :base_64 + raise SeqError, "unknown quality score encoding: #{to}" unless to == :base_33 or to == :base_64 + if from == :base_33 and to == :base_64 + na_qual = NArray.to_na(self.qual, "byte") + na_qual += 64 - 33 + self.qual = na_qual.to_s + elsif from == :base_64 and to == :base_33 + self.qual.tr!("[;-?]", "@") # Handle negative Solexa values from -5 to -1 (set these to 0). + na_qual = NArray.to_na(self.qual, "byte") + na_qual -= 64 - 33 self.qual = na_qual.to_s end @@ -607,6 +622,7 @@ class Seq na_qual = NArray.to_na(self.qual, "byte") na_qual -= SCORE_BASE + na_qual.mean end @@ -618,22 +634,26 @@ class Seq regex_start = Regexp.new(start_codons.join('|'), true) regex_stop = Regexp.new(stop_codons.join('|'), true) - while pos_beg and pos_beg < self.length - size_min - if pos_beg = self.seq.index(regex_start, pos_beg) - if pos_end = self.seq.index(regex_stop, pos_beg) - length = (pos_end - pos_beg) + 3 + while pos_beg = self.seq.index(regex_start, pos_beg) + pos_end = pos_beg + 3 + + while pos_end = self.seq.index(regex_stop, pos_end) + length = (pos_end - pos_beg) + 3 - if (length % 3) == 0 - if size_min <= length and length <= size_max - subseq = self.subseq(pos_beg, length) + if (length % 3) == 0 + if size_min <= length and length <= size_max + subseq = self[pos_beg ... pos_beg + length] - orfs << [subseq, pos_beg, pos_end + 3] - end + orfs << [subseq, pos_beg, pos_end + 3] end + + break end - pos_beg += 1 + pos_end += 1 end + + pos_beg += 1 end if pick_longest