code_ruby/Maasha/lib/seq.rb

   1 # Residue alphabets
   2 DNA     = %w[a t c g]
   3 RNA     = %w[a u c g]
   4 PROTEIN = %w[f l s y c w p h q r i m t n k v a d e g]
   5
   6 # Quality scores bases
   7 SCORE_PHRED  = 33
   8 SCORE_SOLEXA = 64
   9
  10 # Error class for all exceptions to do with Seq.
  11 class SeqError < StandardError; end
  12
  13 class Seq
  14   attr_accessor :seq_name, :seq, :type, :qual
  15
  16   # Initialize a sequence object with the following arguments:
  17   # - seq_name: Name of the sequence.
  18   # - seq: The sequence.
  19   # - type: The sequence type - DNA, RNA, or protein
  20   # - qual: A Solexa type quality scores string.
  21   def initialize(seq_name = nil, seq = nil, type = nil, qual = nil)
  22     @seq_name = seq_name
  23     @seq      = seq
  24     @type     = type
  25     @qual     = qual
  26   end
  27
  28   # Returns the length of a sequence.
  29   def length
  30     self.seq.nil? ? 0 : self.seq.length
  31   end
  32
  33   alias len length
  34
  35   # Method that returns true is a given sequence type is DNA.
  36   def is_dna?
  37     self.type == 'dna'
  38   end
  39
  40   # Method that returns true is a given sequence type is RNA.
  41   def is_rna?
  42     self.type == 'rna'
  43   end
  44
  45   # Method that returns true is a given sequence type is protein.
  46   def is_protein?
  47     self.type == 'protein'
  48   end
  49
  50   # Method to transcribe DNA to RNA.
  51   def to_rna
  52     raise SeqError, "Cannot transcribe 0 length sequence" if self.length == 0
  53     raise SeqError, "Cannot transcribe sequence type: #{self.type}" unless self.is_dna?
  54     self.type = 'rna'
  55     self.seq.tr!('Tt','Uu')
  56   end
  57
  58   # Method to reverse-transcribe RNA to DNA.
  59   def to_dna
  60     raise SeqError, "Cannot reverse-transcribe 0 length sequence" if self.length == 0
  61     raise SeqError, "Cannot reverse-transcribe sequence type: #{self.type}" unless self.is_rna?
  62     self.type = 'dna'
  63     self.seq.tr!('Uu','Tt')
  64   end
  65
  66   # Method that given a Seq entry returns a Biopieces record (a hash).
  67   def to_bp
  68     raise SeqError, "Missing seq_name" if self.seq_name.nil?
  69     raise SeqError, "Missing seq"      if self.seq.nil?
  70     record             = {}
  71     record[:SEQ_NAME] = self.seq_name
  72     record[:SEQ]      = self.seq
  73     record[:SEQ_LEN]  = self.length
  74     record[:SCORES]   = self.qual if self.qual
  75     record
  76   end
  77
  78   # Method that given a Seq entry returns a FASTA entry (a string).
  79   def to_fasta(wrap = nil)
  80     raise SeqError, "Missing seq_name" if self.seq_name.nil?
  81     raise SeqError, "Missing seq"      if self.seq.nil?
  82
  83     seq_name = self.seq_name
  84     seq      = self.seq
  85
  86     unless wrap.nil?
  87       seq.gsub! /(.{#{wrap}})/ do |match|
  88         match << "\n"
  89       end
  90
  91       seq.chomp!
  92     end
  93
  94     ">#{seq_name}\n#{seq}\n"
  95   end
  96
  97   # Method to reverse complement sequence.
  98   def reverse_complement
  99     self.reverse
 100     self.complement
 101   end
 102
 103   alias revcomp reverse_complement
 104
 105   # Method to reverse the sequence.
 106   def reverse
 107     self.seq.reverse!
 108   end
 109
 110   # Method that complements sequence including ambiguity codes.
 111   def complement
 112     raise SeqError, "Cannot complement 0 length sequence" if self.length == 0
 113
 114     if self.is_dna?
 115       self.seq.tr!( 'AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'TCGAAYRWSKMDHBVNtcgaayrwskmdhbvn' )
 116     elsif self.is_rna?
 117       self.seq.tr!( 'AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'UCGAAYRWSKMDHBVNucgaayrwskmdhbvn' )
 118     else
 119       raise SeqError, "Cannot complement sequence type: #{self.type}"
 120     end
 121   end
 122
 123   # Method that generates a random sequence of a given length.
 124   def generate(length,type)
 125     raise SeqError, "Cannot generate negative sequence length: #{length}" if length <= 0
 126
 127     case type.downcase
 128     when "dna"
 129       alph = DNA
 130     when "rna"
 131       alph = RNA
 132     when "protein"
 133       alph = PROTEIN
 134     else
 135       raise SeqError, "Unknown sequence type: #{type}"
 136     end
 137
 138     seq_new   = Array.new(length) { alph[rand(alph.size)] }.join("")
 139     self.seq  = seq_new
 140     self.type = type.downcase
 141     seq_new
 142   end
 143
 144   # Method to convert the quality scores from a specified base
 145   # to another base.
 146   def convert!
 147     self.qual.gsub! /./ do |score|
 148       score_phred  = score.ord - SCORE_PHRED
 149       raise SeqError, "Bad Phred score: #{score} (#{score_phred})" unless (0 .. 40).include? score_phred
 150       score_solexa = score_phred + SCORE_SOLEXA
 151       score        = score_solexa.chr
 152     end
 153   end
 154 end
 155
 156 # Error class for all exceptions to do with Digest.
 157 class DigestError < StandardError; end
 158
 159 class Digest
 160   include Enumerable
 161
 162   # Initialize a digest object with the following arguments:
 163   # - seq: A sequence object.
 164   # - pattern: A restriction enzyme recognition pattern.
 165   # - cut_pos: Offset from match position where the enzyme cuts.
 166   def initialize(seq, pattern, cut_pos)
 167     @seq     = seq
 168     @pattern = disambiguate(pattern)
 169     @cut_pos = cut_pos
 170     @offset  = 0
 171   end
 172
 173   # Method to get the next digestion product from a sequence.
 174   def each
 175     @seq.seq.scan @pattern do
 176       pos = $`.length + @cut_pos - 1
 177
 178       if pos >= 0 and pos < @seq.seq.length - 2
 179         seq = Seq.new("#{@seq.seq_name}[#{@offset + 1}-#{pos + 1}]", @seq.seq[@offset .. pos], @seq.type)
 180
 181         yield seq
 182       end
 183
 184       @offset = pos + 1
 185     end
 186
 187     if @offset < 0
 188       @offset = 0
 189     elsif @offset > @seq.seq.length
 190       @offset = 0
 191     end
 192
 193     seq = Seq.new("#{@seq.seq_name}[#{@offset + 1}-#{@seq.seq.length}]", @seq.seq[@offset .. @seq.seq.length], @seq.type)
 194
 195     yield seq
 196
 197     self # conventionally
 198   end
 199
 200   private
 201
 202   # Method that returns a regexp object with a restriction
 203   # enzyme pattern with ambiguity codes substituted to the
 204   # appropriate regexp.
 205   def disambiguate(pattern)
 206     ambiguity = {
 207       'A' => "A",
 208       'T' => "T",
 209       'U' => "T",
 210       'C' => "C",
 211       'G' => "G",
 212       'M' => "[AC]",
 213       'R' => "[AG]",
 214       'W' => "[AT]",
 215       'S' => "[CG]",
 216       'Y' => "[CT]",
 217       'K' => "[GT]",
 218       'V' => "[ACG]",
 219       'H' => "[ACT]",
 220       'D' => "[AGT]",
 221       'B' => "[CGT]",
 222       'N' => "[GATC]"
 223     }
 224
 225     new_pattern = ""
 226
 227     pattern.upcase.each_char do |char|
 228       if ambiguity.has_key? char
 229         new_pattern << ambiguity[char]
 230       else
 231         raise DigestError, "Could not disambiguate residue: #{char}"
 232       end
 233     end
 234
 235     Regexp.new(new_pattern)
 236   end
 237 end
 238
 239
 240 __END__
 241
 242
 243
 244
 245
 246 # Class containing generic sequence methods and nucleic acid and amino acid subclasses.
 247 class Seq < String
 248   # Guess the sequence type by analyzing the first 100 residues allowing for ambiguity codes.
 249   def guess_type
 250     raise ArgumentError, "No sequence." if self.empty?
 251
 252     seq_beg = self[0, 100].upcase
 253
 254     if seq_beg.count( "FLPQIE" ) > 0
 255       Seq::AA.new(self)
 256     elsif seq_beg.count("U") > 0
 257       Seq::NA::RNA.new(self)
 258     else
 259       Seq::NA::DNA.new(self)
 260     end
 261   end
 262
 263   # Method to wrap a sequence to a given width using a given delimiter.
 264   def wrap(width = 80, delimit = $/)
 265     raise ArgumentError, "Cannot wrap sequence to negative width: #{width}." if width <= 0
 266
 267     self.delete!(" \t\n\r")
 268     self.gsub(/.{#{width}}(?!$)/, "\\0#{delimit}")
 269   end
 270
 271   # Method to wrap and replace a sequence to a given width using a given delimiter.
 272   def wrap!(width = 80, delimit = $/)
 273     self.replace(self.wrap(width, delimit))
 274   end
 275
 276   # Method that generates a random sequence of a given length.
 277   def generate(length)
 278     raise ArgumentError, "Cannot generate negative sequence length: #{length}." if length <= 0
 279
 280     alph = self.residues
 281     Array.new(length) { alph[rand(alph.size)] }.join("")
 282   end
 283
 284   # Method that replaces sequence with a randomly generated sequence of a given length.
 285   def generate!(length)
 286     self.replace(self.generate(length))
 287   end
 288
 289   # Class containing methods specific for amino acid (AA) sequences.
 290   class AA < Seq
 291     # Method that returns an array of amino acid residues.
 292     def residues
 293       %w{ F L S Y C W P H Q R I M T N K V A D E G }
 294     end
 295
 296     # Calculate the molecular weight of an amino acid seuqunce.
 297     # The caluculation is only approximate since there is no correction
 298     # for amino bond formation and the MW used are somewhat imprecise:
 299     # http://www.expasy.ch/tools/pscale/Molecularweight.html
 300     def mol_weight
 301       raise ArgumentError, "invalid residues found: #{self.delete("#{residues.join( "" )}")}" if self.upcase =~ /[^#{residues.join( "" )}]/
 302
 303       mol_weight_aa = {
 304         "A" => 89.000,    # Ala
 305         "R" => 174.000,   # Arg
 306         "N" => 132.000,   # Asn
 307         "D" => 133.000,   # Asp
 308         "C" => 121.000,   # Cys
 309         "Q" => 146.000,   # Gln
 310         "E" => 147.000,   # Glu
 311         "G" => 75.000,    # Gly
 312         "H" => 155.000,   # His
 313         "I" => 131.000,   # Ile
 314         "L" => 131.000,   # Leu
 315         "K" => 146.000,   # Lys
 316         "M" => 149.000,   # Met
 317         "F" => 165.000,   # Phe
 318         "P" => 115.000,   # Pro
 319         "S" => 105.000,   # Ser
 320         "T" => 119.000,   # Thr
 321         "W" => 204.000,   # Trp
 322         "Y" => 181.000,   # Tyr
 323         "V" => 117.000,   # Val
 324       }
 325
 326       mw = 0.0
 327
 328       self.upcase.each_char { |c| mw += mol_weight_aa[ c ] }
 329
 330       mw
 331     end
 332   end