code_ruby/lib/maasha/seq.rb

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