X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=code_ruby%2Flib%2Fmaasha%2Fseq.rb;h=cf4b6e6802883d5447c7bbad7157ae65bc0c95b1;hb=80706855913e8607239241fd2938f47d1fb76799;hp=b919a7f31cff1a32555960e4fa60a03f4668eb34;hpb=93cdd1c3b5b94db85707fbffe84be90095cce1b4;p=biopieces.git diff --git a/code_ruby/lib/maasha/seq.rb b/code_ruby/lib/maasha/seq.rb index b919a7f..cf4b6e6 100644 --- a/code_ruby/lib/maasha/seq.rb +++ b/code_ruby/lib/maasha/seq.rb @@ -23,13 +23,13 @@ # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< require 'maasha/bits' -require 'maasha/seq/backtrack' require 'maasha/seq/digest' -#require 'maasha/seq/patscan' -require 'maasha/seq/patternmatcher' require 'maasha/seq/trim' require 'narray' +autoload :BackTrack, 'maasha/seq/backtrack.rb' +autoload :Dynamic, 'maasha/seq/dynamic.rb' + # Residue alphabets DNA = %w[a t c g] RNA = %w[a u c g] @@ -67,18 +67,16 @@ 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 - #include Patscan - include PatternMatcher - include BackTrack + # Quality scores bases + SCORE_BASE = 33 + SCORE_MIN = 0 + SCORE_MAX = 40 + include Digest include Trim @@ -152,6 +150,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. @@ -336,17 +335,13 @@ class Seq key end - # Method to reverse complement sequence. - def reverse_complement - self.reverse - self.complement - self + # Method to reverse the sequence. + def reverse + Seq.new(self.seq_name, self.seq.reverse, self.type, self.qual ? self.qual.reverse : self.qual) end - alias :revcomp :reverse_complement - # Method to reverse the sequence. - def reverse + def reverse! self.seq.reverse! self.qual.reverse! if self.qual self @@ -356,6 +351,26 @@ class Seq def complement raise SeqError, "Cannot complement 0 length sequence" if self.length == 0 + entry = Seq.new + entry.seq_name = self.seq_name + entry.type = self.type + entry.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? @@ -363,6 +378,8 @@ class Seq else raise SeqError, "Cannot complement sequence type: #{self.type}" end + + self end # Method to determine the Hamming Distance between @@ -392,12 +409,28 @@ class Seq seq_new end - # Method to shuffle a sequence readomly inline. + # Method to return a new Seq object with shuffled sequence. + def shuffle + Seq.new(self.seq_name, self.seq.split('').shuffle!.join, self.type, self.qual) + end + + # Method to shuffle a sequence randomly inline. def shuffle! self.seq = self.seq.split('').shuffle!.join self end + # 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 + + self.seq << entry.seq + self.qual << entry.qual unless entry.qual.nil? + + 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) @@ -415,25 +448,18 @@ class Seq qual = self.qual[start .. stop] unless self.qual.nil? end - Seq.new(self.seq_name.dup, seq, self.type, qual) + Seq.new(self.seq_name, seq, self.type, qual) # TODO changed self.seq_name.dup to self.seq_name -> consequence? 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) - 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 - - if length == 0 - self.seq = "" - self.qual = "" unless self.qual.nil? - else - stop = start + length - 1 + s = subseq(start, length) - self.seq = self.seq[start .. stop] - self.qual = self.qual[start .. stop] unless self.qual.nil? - end + self.seq_name = s.seq_name + self.seq = s.seq + self.type = s.type + self.qual = s.qual self end @@ -492,21 +518,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) @@ -545,34 +556,56 @@ class Seq self end - # Method to convert quality scores inbetween formats. - # Sanger base 33, range 0-40 - # Solexa base 64, range -5-40 - # Illumina13 base 64, range 0-40 - # Illumina15 base 64, range 3-40 - # Illumina18 base 33, range 0-41 - def convert_scores!(from, to) - unless from == to - na_qual = NArray.to_na(self.qual, "byte") + # Method that determines if a quality score string can be + # absolutely identified as base 33. + def qual_base33? + self.qual.match(/[!-:]/) ? true : false + end + + # Method that determines if a quality score string may be base 64. + def qual_base64? + self.qual.match(/[K-h]/) ? true : false + end - case from.downcase - when "sanger" then na_qual -= 33 - 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 + # 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 to.downcase - when "sanger" then na_qual += 33 - 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 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 coerce quality scores to be within the 0-40 range. + def qual_coerce!(encoding) + raise SeqError, "Missing qual" if self.qual.nil? + 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 + + 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 @@ -585,6 +618,7 @@ class Seq na_qual = NArray.to_na(self.qual, "byte") na_qual -= SCORE_BASE + na_qual.mean end