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