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