X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=code_ruby%2Flib%2Fmaasha%2Fseq.rb;h=92ce19edaee4a9176487272a8f1a89c32c7ee806;hb=a945a21356f87f679f433279e502676f425ac8ee;hp=4a0af47ed0f5fde9d5ffe87ca19a3caec087f6d6;hpb=b982cc677363d7458963b610ec53bf2c4f7476a9;p=biopieces.git diff --git a/code_ruby/lib/maasha/seq.rb b/code_ruby/lib/maasha/seq.rb index 4a0af47..92ce19e 100644 --- a/code_ruby/lib/maasha/seq.rb +++ b/code_ruby/lib/maasha/seq.rb @@ -73,7 +73,7 @@ class SeqError < StandardError; end class Seq # Quality scores bases - SCORE_BASE = 64 + SCORE_BASE = 33 SCORE_MIN = 0 SCORE_MAX = 40 @@ -87,7 +87,7 @@ 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) @@ -98,9 +98,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 @@ -140,9 +140,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 @@ -190,24 +190,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 @@ -215,14 +215,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 @@ -257,7 +256,7 @@ class Seq self.seq = protein self.qual = nil - self.type = "protein" + self.type = :protein self end @@ -392,20 +391,17 @@ class Seq 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 @@ -420,6 +416,17 @@ class Seq 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) @@ -550,61 +557,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 @@ -617,6 +614,7 @@ class Seq na_qual = NArray.to_na(self.qual, "byte") na_qual -= SCORE_BASE + na_qual.mean end