X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=code_ruby%2Flib%2Fmaasha%2Fseq.rb;h=e4b6bfd7c27d5933d0bd9c6b65210bd9a7b32ea0;hb=b2ea0b5a51a558478af60e9df4c643dd58552086;hp=394d6fb3ca205b9ee5b86757efdc88ff4ca9ab30;hpb=8d629b585e06ef165de5447895a95fc9230815fd;p=biopieces.git diff --git a/code_ruby/lib/maasha/seq.rb b/code_ruby/lib/maasha/seq.rb index 394d6fb..e4b6bfd 100644 --- a/code_ruby/lib/maasha/seq.rb +++ b/code_ruby/lib/maasha/seq.rb @@ -1,4 +1,4 @@ -# Copyright (C) 2007-2011 Martin A. Hansen. +# Copyright (C) 2007-2012 Martin A. Hansen. # This program is free software; you can redistribute it and/or # modify it under the terms of the GNU General Public License @@ -22,10 +22,15 @@ # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< -require 'maasha/digest' -require 'maasha/patternmatcher' require 'maasha/bits' -#require 'maasha/patscan' +require 'maasha/seq/digest' +require 'maasha/seq/trim' +require 'narray' + +autoload :BackTrack, 'maasha/seq/backtrack' +autoload :Dynamic, 'maasha/seq/dynamic' +autoload :Homopolymer, 'maasha/seq/homopolymer' +autoload :Levenshtein, 'maasha/seq/levenshtein' # Residue alphabets DNA = %w[a t c g] @@ -33,16 +38,48 @@ RNA = %w[a u c g] PROTEIN = %w[f l s y c w p h q r i m t n k v a d e g] INDELS = %w[. - _ ~] -# Quality scores bases -SCORE_PHRED = 33 -SCORE_ILLUMINA = 64 +# Translation table 11 +# (http://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=cgencodes#SG11) +# AAs = FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG +# Starts = ---M---------------M------------MMMM---------------M------------ +# Base1 = TTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGG +# Base2 = TTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGG +# Base3 = TCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAG +TRANS_TAB11_START = { + "TTG" => "M", "CTG" => "M", "ATT" => "M", "ATC" => "M", + "ATA" => "M", "ATG" => "M", "GTG" => "M" +} + +TRANS_TAB11 = { + "TTT" => "F", "TCT" => "S", "TAT" => "Y", "TGT" => "C", + "TTC" => "F", "TCC" => "S", "TAC" => "Y", "TGC" => "C", + "TTA" => "L", "TCA" => "S", "TAA" => "*", "TGA" => "*", + "TTG" => "L", "TCG" => "S", "TAG" => "*", "TGG" => "W", + "CTT" => "L", "CCT" => "P", "CAT" => "H", "CGT" => "R", + "CTC" => "L", "CCC" => "P", "CAC" => "H", "CGC" => "R", + "CTA" => "L", "CCA" => "P", "CAA" => "Q", "CGA" => "R", + "CTG" => "L", "CCG" => "P", "CAG" => "Q", "CGG" => "R", + "ATT" => "I", "ACT" => "T", "AAT" => "N", "AGT" => "S", + "ATC" => "I", "ACC" => "T", "AAC" => "N", "AGC" => "S", + "ATA" => "I", "ACA" => "T", "AAA" => "K", "AGA" => "R", + "ATG" => "M", "ACG" => "T", "AAG" => "K", "AGG" => "R", + "GTT" => "V", "GCT" => "A", "GAT" => "D", "GGT" => "G", + "GTC" => "V", "GCC" => "A", "GAC" => "D", "GGC" => "G", + "GTA" => "V", "GCA" => "A", "GAA" => "E", "GGA" => "G", + "GTG" => "V", "GCG" => "A", "GAG" => "E", "GGG" => "G" +} # Error class for all exceptions to do with Seq. class SeqError < StandardError; end class Seq - #include Patscan - include PatternMatcher + # Quality scores bases + SCORE_BASE = 33 + SCORE_MIN = 0 + SCORE_MAX = 40 + + include Digest + include Trim attr_accessor :seq_name, :seq, :type, :qual @@ -51,7 +88,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) @@ -62,9 +99,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 @@ -104,9 +141,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 @@ -114,6 +151,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. @@ -129,26 +167,48 @@ class Seq self.seq.scan(regex).size end + # Method to remove indels from seq and qual if qual. + def indels_remove + if self.qual.nil? + self.seq.delete!(Regexp.escape(INDELS.join(''))) + else + na_seq = NArray.to_na(self.seq, "byte") + na_qual = NArray.to_na(self.qual, "byte") + mask = NArray.byte(self.length) + + INDELS.each do |c| + mask += na_seq.eq(c.ord) + end + + mask = mask.eq(0) + + self.seq = na_seq[mask].to_s + self.qual = na_qual[mask].to_s + end + + self + end + # 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 @@ -156,11 +216,60 @@ 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 length must be a multiplum of 3 - was: #{self.length}" unless (self.length % 3) == 0 + + case trans_tab + when 11 + codon_start_hash = TRANS_TAB11_START + codon_hash = TRANS_TAB11 + else + raise SeqError, "Unknown translation table: #{trans_tab}" + end + + codon = self.seq[0 ... 3].upcase + + aa = codon_start_hash[codon] + + raise SeqError, "Unknown start codon: #{codon}" if aa.nil? + + protein = aa + + i = 3 + + while i < self.length + codon = self.seq[i ... i + 3].upcase + + aa = codon_hash[codon] + + raise SeqError, "Unknown codon: #{codon}" if aa.nil? + + protein << aa + + i += 3 + end + + self.seq = protein + self.qual = nil + self.type = :protein + + self + end + + alias :to_protein! :translate! + + def translate(trans_tab = 11) + self.dup.translate!(trans_tab) + end + + alias :to_protein :translate + # Method that given a Seq entry returns a Biopieces record (a hash). def to_bp raise SeqError, "Missing seq_name" if self.seq_name.nil? @@ -176,8 +285,8 @@ class Seq # Method that given a Seq entry returns a FASTA entry (a string). def to_fasta(wrap = nil) - raise SeqError, "Missing seq_name" if self.seq_name.nil? - raise SeqError, "Missing seq" if self.seq.nil? + 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 @@ -226,90 +335,166 @@ class Seq key end - # Method to reverse complement sequence. - def reverse_complement - self.reverse - self.complement + # 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 end # Method that complements sequence including ambiguity codes. 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? - self.seq.tr!( 'AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'TCGAAYRWSKMDHBVNtcgaayrwskmdhbvn' ) + entry.seq = self.seq.tr('AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'TCGAAYRWSKMDHBVNtcgaayrwskmdhbvn') elsif self.is_rna? - self.seq.tr!( 'AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'UCGAAYRWSKMDHBVNucgaayrwskmdhbvn' ) + 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? + self.seq.tr!('AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'UCGAAYRWSKMDHBVNucgaayrwskmdhbvn') + 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) + self.seq.upcase.hamming_distance(entry.seq.upcase) + 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(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 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 + + # 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 + + # Index method for Seq objects. + def [](*args) + entry = Seq.new + entry.seq_name = self.seq_name + entry.seq = self.seq[*args] + entry.type = self.type + entry.qual = self.qual[*args] unless self.qual.nil? + + 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 # 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} < 1" if length <= 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 - stop = start + length - 1 + if length == 0 + seq = "" + qual = "" unless self.qual.nil? + else + stop = start + length - 1 + + seq = self.seq[start .. stop] + qual = self.qual[start .. stop] unless self.qual.nil? + end - seq = self.seq[start .. stop] - qual = self.qual[start .. stop] unless self.qual.nil? + seq_name = self.seq_name.nil? ? nil : self.seq_name.dup - Seq.new(self.seq_name, seq, self.type, qual) + Seq.new(seq_name, seq, self.type, qual) 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} < 1" if length <= 0 - raise SeqError, "subsequence start + length > Seq.length: #{start} + #{length} > #{self.length}" if start + length > self.length + s = subseq(start, length) - stop = start + length - 1 + self.seq_name = s.seq_name + self.seq = s.seq + self.type = s.type + self.qual = s.qual - self.seq = self.seq[start .. stop] - self.qual = self.qual[start .. stop] unless self.qual.nil? + self end # Method that returns a subsequence of a given length @@ -337,23 +522,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 @@ -366,39 +534,149 @@ class Seq ((self.seq.scan(/[a-z]/).size.to_f / (self.len - self.indels).to_f) * 100).round(2) end - # Method to convert the quality scores from a specified base - # to another base. - def convert_phred2illumina! - self.qual.gsub!(/./) do |score| - score_phred = score.ord - SCORE_PHRED - raise SeqError, "Bad Phred score: #{score} (#{score_phred})" unless (0 .. 41).include? score_phred - score_illumina = score_phred + SCORE_ILLUMINA - score = score_illumina.chr + # Hard masks sequence residues where the corresponding quality score + # is below a given cutoff. + def mask_seq_hard!(cutoff) + raise SeqError, "seq is nil" if self.seq.nil? + raise SeqError, "qual is nil" if self.qual.nil? + raise SeqError, "cufoff value: #{cutoff} out of range #{SCORE_MIN} .. #{SCORE_MAX}" unless (SCORE_MIN .. SCORE_MAX).include? cutoff + + na_seq = NArray.to_na(self.seq, "byte") + na_qual = NArray.to_na(self.qual, "byte") + mask = (na_qual - SCORE_BASE) < cutoff + mask *= na_seq.ne("-".ord) + + na_seq[mask] = 'N'.ord + + self.seq = na_seq.to_s + + self + end + + # Soft masks sequence residues where the corresponding quality score + # is below a given cutoff. + def mask_seq_soft!(cutoff) + raise SeqError, "seq is nil" if self.seq.nil? + raise SeqError, "qual is nil" if self.qual.nil? + raise SeqError, "cufoff value: #{cutoff} out of range #{SCORE_MIN} .. #{SCORE_MAX}" unless (SCORE_MIN .. SCORE_MAX).include? cutoff + + na_seq = NArray.to_na(self.seq, "byte") + na_qual = NArray.to_na(self.qual, "byte") + mask = (na_qual - SCORE_BASE) < cutoff + mask *= na_seq.ne("-".ord) + + na_seq[mask] ^= ' '.ord + + self.seq = na_seq.to_s + + self + end + + # 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 + + # 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 + 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 the quality scores from Solexa odd/ratio to - # Illumina format. - def convert_solexa2illumina! - self.qual.gsub!(/./) do |score| - score = solexa_char2illumina_char(score) + # 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 + + self end - private + # Method to calculate and return the mean quality score. + def scores_mean + raise SeqError, "Missing qual in entry" if self.qual.nil? + + na_qual = NArray.to_na(self.qual, "byte") + na_qual -= SCORE_BASE - # Method to convert a Solexa score (odd ratio) to - # a phred (probability) integer score. - def solexa2phred(score) - (10.0 * Math.log(10.0 ** (score / 10.0) + 1.0, 10)).to_i + na_qual.mean end - # Method to convert a Solexa score encoded using base - # 64 ASCII to a Phred score encoded using base 64 ASCII. - def solexa_char2illumina_char(char) - score_solexa = char.ord - 64 - score_phred = solexa2phred(score_solexa) - (score_phred + 64).chr + # Method to find open reading frames (ORFs). + def each_orf(size_min, size_max, start_codons, stop_codons, pick_longest = false) + orfs = [] + pos_beg = 0 + + 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 + + if (length % 3) == 0 + if size_min <= length and length <= size_max + subseq = self.subseq(pos_beg, length) + + orfs << [subseq, pos_beg, pos_end + 3] + end + end + end + + pos_beg += 1 + end + end + + if pick_longest + orf_hash = {} + + orfs.each { |orf| orf_hash[orf.last] = orf unless orf_hash[orf.last] } + + orfs = orf_hash.values + end + + if block_given? + orfs.each { |orf| yield orf } + else + return orfs + end end end