code_ruby/lib/maasha/seq.rb

   1 # Copyright (C) 2007-2011 Martin A. Hansen.
   2
   3 # This program is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU General Public License
   5 # as published by the Free Software Foundation; either version 2
   6 # of the License, or (at your option) any later version.
   7
   8 # This program is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 # GNU General Public License for more details.
  12
  13 # You should have received a copy of the GNU General Public License
  14 # along with this program; if not, write to the Free Software
  15 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  16
  17 # http://www.gnu.org/copyleft/gpl.html
  18
  19 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  20
  21 # This software is part of the Biopieces framework (www.biopieces.org).
  22
  23 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  24
  25 require 'maasha/digest'
  26 require 'maasha/patternmatcher'
  27 require 'maasha/bits'
  28 #require 'maasha/patscan'
  29
  30 # Residue alphabets
  31 DNA     = %w[a t c g]
  32 RNA     = %w[a u c g]
  33 PROTEIN = %w[f l s y c w p h q r i m t n k v a d e g]
  34 INDELS  = %w[. - _ ~]
  35
  36 # Quality scores bases
  37 SCORE_PHRED    = 33
  38 SCORE_ILLUMINA = 64
  39
  40 # Error class for all exceptions to do with Seq.
  41 class SeqError < StandardError; end
  42
  43 class Seq
  44   #include Patscan
  45   include PatternMatcher
  46
  47   attr_accessor :seq_name, :seq, :type, :qual
  48
  49   # Class method to instantiate a new Sequence object given
  50   # a Biopiece record.
  51   def self.new_bp(record)
  52     seq_name = record[:SEQ_NAME]
  53     seq      = record[:SEQ]
  54     type     = record[:SEQ_TYPE]
  55     qual     = record[:SCORES]
  56
  57     self.new(seq_name, seq, type, qual)
  58   end
  59
  60   # Class method that generates all possible oligos of a specifed length and type.
  61   def self.generate_oligos(length, type)
  62     raise SeqError, "Cannot generate negative oligo length: #{length}" if length <= 0
  63
  64     case type.downcase
  65     when /dna/     then alph = DNA
  66     when /rna/     then alph = RNA
  67     when /protein/ then alph = PROTEIN
  68     else
  69       raise SeqError, "Unknown sequence type: #{type}"
  70     end
  71
  72     oligos = [""]
  73
  74     (1 .. length).each do
  75       list = []
  76
  77       oligos.each do |oligo|
  78         alph.each do |char|
  79           list << oligo + char
  80         end
  81       end
  82
  83       oligos = list
  84     end
  85
  86     oligos
  87   end
  88
  89   # Initialize a sequence object with the following arguments:
  90   # - seq_name: Name of the sequence.
  91   # - seq: The sequence.
  92   # - type: The sequence type - DNA, RNA, or protein
  93   # - qual: An Illumina type quality scores string.
  94   def initialize(seq_name = nil, seq = nil, type = nil, qual = nil)
  95     @seq_name = seq_name
  96     @seq      = seq
  97     @type     = type
  98     @qual     = qual
  99   end
 100
 101   # Method that guesses and returns the sequence type
 102   # by inspecting the first 100 residues.
 103   def type_guess
 104     raise SeqError, "Guess failed: sequence is nil" if self.seq.nil?
 105
 106     case self.seq[0 ... 100].downcase
 107     when /[flpqie]/ then return "protein"
 108     when /[u]/      then return "rna"
 109     else                 return "dna"
 110     end
 111   end
 112
 113   # Method that guesses and sets the sequence type
 114   # by inspecting the first 100 residues.
 115   def type_guess!
 116     self.type = self.type_guess
 117   end
 118
 119   # Returns the length of a sequence.
 120   def length
 121     self.seq.nil? ? 0 : self.seq.length
 122   end
 123
 124   alias :len :length
 125
 126   # Return the number indels in a sequence.
 127   def indels
 128     regex = Regexp.new(/[#{Regexp.escape(INDELS.join(""))}]/)
 129     self.seq.scan(regex).size
 130   end
 131
 132   # Method that returns true is a given sequence type is DNA.
 133   def is_dna?
 134     self.type == 'dna'
 135   end
 136
 137   # Method that returns true is a given sequence type is RNA.
 138   def is_rna?
 139     self.type == 'rna'
 140   end
 141
 142   # Method that returns true is a given sequence type is protein.
 143   def is_protein?
 144     self.type == 'protein'
 145   end
 146
 147   # Method to transcribe DNA to RNA.
 148   def to_rna
 149     raise SeqError, "Cannot transcribe 0 length sequence" if self.length == 0
 150     raise SeqError, "Cannot transcribe sequence type: #{self.type}" unless self.is_dna?
 151     self.type = 'rna'
 152     self.seq.tr!('Tt','Uu')
 153   end
 154
 155   # Method to reverse-transcribe RNA to DNA.
 156   def to_dna
 157     raise SeqError, "Cannot reverse-transcribe 0 length sequence" if self.length == 0
 158     raise SeqError, "Cannot reverse-transcribe sequence type: #{self.type}" unless self.is_rna?
 159
 160     self.type = 'dna'
 161     self.seq.tr!('Uu','Tt')
 162   end
 163
 164   # Method that given a Seq entry returns a Biopieces record (a hash).
 165   def to_bp
 166     raise SeqError, "Missing seq_name" if self.seq_name.nil?
 167     raise SeqError, "Missing seq"      if self.seq.nil?
 168
 169     record             = {}
 170     record[:SEQ_NAME] = self.seq_name
 171     record[:SEQ]      = self.seq
 172     record[:SEQ_LEN]  = self.length
 173     record[:SCORES]   = self.qual if self.qual
 174     record
 175   end
 176
 177   # Method that given a Seq entry returns a FASTA entry (a string).
 178   def to_fasta(wrap = nil)
 179     raise SeqError, "Missing seq_name" if self.seq_name.nil?
 180     raise SeqError, "Missing seq"      if self.seq.nil?
 181
 182     seq_name = self.seq_name.to_s
 183     seq      = self.seq.to_s
 184
 185     unless wrap.nil?
 186       seq.gsub!(/(.{#{wrap}})/) do |match|
 187         match << $/
 188       end
 189
 190       seq.chomp!
 191     end
 192
 193     ">" + seq_name + $/ + seq + $/
 194   end
 195
 196   # Method that given a Seq entry returns a FASTQ entry (a string).
 197   def to_fastq
 198     raise SeqError, "Missing seq_name" if self.seq_name.nil?
 199     raise SeqError, "Missing seq"      if self.seq.nil?
 200     raise SeqError, "Missing qual"     if self.qual.nil?
 201
 202     seq_name = self.seq_name.to_s
 203     seq      = self.seq.to_s
 204     qual     = self.qual.to_s
 205
 206     "@" + seq_name + $/ + seq + $/ + "+" + $/ + qual + $/
 207   end
 208
 209   # Method that generates a unique key for a
 210   # DNA sequence and return this key as a Fixnum.
 211   def to_key
 212     key = 0
 213
 214     self.seq.upcase.each_char do |char|
 215       key <<= 2
 216
 217       case char
 218       when 'A' then key |= 0
 219       when 'C' then key |= 1
 220       when 'G' then key |= 2
 221       when 'T' then key |= 3
 222       else raise SeqError, "Bad residue: #{char}"
 223       end
 224     end
 225
 226     key
 227   end
 228
 229   # Method to reverse complement sequence.
 230   def reverse_complement
 231     self.reverse
 232     self.complement
 233     self
 234   end
 235
 236   alias :revcomp :reverse_complement
 237
 238   # Method to reverse the sequence.
 239   def reverse
 240     self.seq.reverse!
 241     self.qual.reverse! if self.qual
 242     self
 243   end
 244
 245   # Method that complements sequence including ambiguity codes.
 246   def complement
 247     raise SeqError, "Cannot complement 0 length sequence" if self.length == 0
 248
 249     if self.is_dna?
 250       self.seq.tr!( 'AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'TCGAAYRWSKMDHBVNtcgaayrwskmdhbvn' )
 251     elsif self.is_rna?
 252       self.seq.tr!( 'AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'UCGAAYRWSKMDHBVNucgaayrwskmdhbvn' )
 253     else
 254       raise SeqError, "Cannot complement sequence type: #{self.type}"
 255     end
 256   end
 257
 258   # Method to determine the Hamming Distance between
 259   # two Sequence objects (case insensitive).
 260   def hamming_distance(seq)
 261     self.seq.upcase.hamming_distance(seq.seq.upcase)
 262   end
 263
 264   # Method that generates a random sequence of a given length and type.
 265   def generate(length, type)
 266     raise SeqError, "Cannot generate sequence length < 1: #{length}" if length <= 0
 267
 268     case type.downcase
 269     when "dna"
 270       alph = DNA
 271     when "rna"
 272       alph = RNA
 273     when "protein"
 274       alph = PROTEIN
 275     else
 276       raise SeqError, "Unknown sequence type: #{type}"
 277     end
 278
 279     seq_new   = Array.new(length) { alph[rand(alph.size)] }.join("")
 280     self.seq  = seq_new
 281     self.type = type.downcase
 282     seq_new
 283   end
 284
 285   # Method to shuffle a sequence readomly inline.
 286   def shuffle!
 287     self.seq = self.seq.split('').shuffle!.join
 288     self
 289   end
 290
 291   # Method that returns a subsequence of from a given start position
 292   # and of a given length.
 293   def subseq(start, length = self.length - start)
 294     raise SeqError, "subsequence start: #{start} < 0"                                                if start  < 0
 295     raise SeqError, "subsequence length: #{length} < 1"                                              if length <= 0
 296     raise SeqError, "subsequence start + length > Seq.length: #{start} + #{length} > #{self.length}" if start + length > self.length
 297
 298     stop = start + length - 1
 299
 300     seq  = self.seq[start .. stop]
 301     qual = self.qual[start .. stop] unless self.qual.nil?
 302
 303     Seq.new(self.seq_name, seq, self.type, qual)
 304   end
 305
 306   # Method that replaces a sequence with a subsequence from a given start position
 307   # and of a given length.
 308   def subseq!(start, length = self.length - start)
 309     raise SeqError, "subsequence start: #{start} < 0"                                                if start  < 0
 310     raise SeqError, "subsequence length: #{length} < 1"                                              if length <= 0
 311     raise SeqError, "subsequence start + length > Seq.length: #{start} + #{length} > #{self.length}" if start + length > self.length
 312
 313     stop = start + length - 1
 314
 315     self.seq  = self.seq[start .. stop]
 316     self.qual = self.qual[start .. stop] unless self.qual.nil?
 317   end
 318
 319   # Method that returns a subsequence of a given length
 320   # beginning at a random position.
 321   def subseq_rand(length)
 322     if self.length - length + 1 == 0
 323       start = 0
 324     else
 325       start = rand(self.length - length + 1)
 326     end
 327
 328     self.subseq(start, length)
 329   end
 330
 331   def quality_trim(min)
 332   end
 333
 334   def quality_trim_right(min)
 335   end
 336
 337   def quality_trim_left(min)
 338   end
 339
 340   def quality_trim!(min)
 341   end
 342
 343   def quality_trim_right!(min)
 344   end
 345
 346   def quality_trim_left!(min)
 347   end
 348
 349   # Method that returns the residue compositions of a sequence in
 350   # a hash where the key is the residue and the value is the residue
 351   # count.
 352   def composition
 353     comp = Hash.new(0);
 354
 355     self.seq.upcase.each_char do |char|
 356       comp[char] += 1
 357     end
 358
 359     comp
 360   end
 361
 362   # Method that returns the length of the longest homopolymeric stretch
 363   # found in a sequence.
 364   def homopol_max(min = 1)
 365     return 0 if self.seq.nil? or self.seq.empty?
 366
 367     found = false
 368
 369     self.seq.upcase.scan(/A{#{min},}|T{#{min},}|G{#{min},}|C{#{min},}|N{#{min},}/) do |match|
 370       found = true
 371       min   = match.size > min ? match.size : min
 372     end
 373
 374     return 0 unless found
 375
 376     min
 377   end
 378
 379   # Method that returns the percentage of hard masked residues
 380   # or N's in a sequence.
 381   def hard_mask
 382     ((self.seq.upcase.scan("N").size.to_f / (self.len - self.indels).to_f) * 100).round(2)
 383   end
 384
 385   # Method that returns the percentage of soft masked residues
 386   # or lower cased residues in a sequence.
 387   def soft_mask
 388     ((self.seq.scan(/[a-z]/).size.to_f / (self.len - self.indels).to_f) * 100).round(2)
 389   end
 390
 391   # Method to convert the quality scores from a specified base
 392   # to another base.
 393   def convert_phred2illumina!
 394     self.qual.gsub!(/./) do |score|
 395       score_phred  = score.ord - SCORE_PHRED
 396       raise SeqError, "Bad Phred score: #{score} (#{score_phred})" unless (0 .. 41).include? score_phred
 397       score_illumina = score_phred + SCORE_ILLUMINA
 398       score          = score_illumina.chr
 399     end
 400   end
 401
 402   # Method to convert the quality scores from Solexa odd/ratio to
 403   # Illumina format.
 404   def convert_solexa2illumina!
 405     self.qual.gsub!(/./) do |score|
 406       score = solexa_char2illumina_char(score)
 407     end
 408   end
 409
 410   private
 411
 412   # Method to convert a Solexa score (odd ratio) to
 413   # a phred (probability) integer score.
 414   def solexa2phred(score)
 415     (10.0 * Math.log(10.0 ** (score / 10.0) + 1.0, 10)).to_i
 416   end
 417
 418   # Method to convert a Solexa score encoded using base
 419   # 64 ASCII to a Phred score encoded using base 64 ASCII.
 420   def solexa_char2illumina_char(char)
 421     score_solexa = char.ord - 64
 422     score_phred  = solexa2phred(score_solexa)
 423     (score_phred + 64).chr
 424   end
 425 end
 426
 427 __END__