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