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