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