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