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     # Handle remaining sequence after last cut.
 172     if @offset > 0 and @offset < @seq.seq.length
 173       seq = Seq.new("#{@seq.seq_name}[#{@offset + 1}-#{@seq.seq.length}]", @seq.seq[@offset .. @seq.seq.length], @seq.type)
 174
 175       yield seq
 176     end
 177
 178     self # conventionally
 179   end
 180
 181   private
 182
 183   # Method that returns a regexp object with a restriction
 184   # enzyme pattern with ambiguity codes substituted to the
 185   # appropriate regexp.
 186   def disambiguate(pattern)
 187     ambiguity = {
 188       'A' => "A",
 189       'T' => "T",
 190       'U' => "T",
 191       'C' => "C",
 192       'G' => "G",
 193       'M' => "[AC]",
 194       'R' => "[AG]",
 195       'W' => "[AT]",
 196       'S' => "[CG]",
 197       'Y' => "[CT]",
 198       'K' => "[GT]",
 199       'V' => "[ACG]",
 200       'H' => "[ACT]",
 201       'D' => "[AGT]",
 202       'B' => "[CGT]",
 203       'N' => "[GATC]"
 204     }
 205
 206     new_pattern = ""
 207
 208     pattern.upcase.each_char do |char|
 209       if ambiguity.has_key? char
 210         new_pattern << ambiguity[char]
 211       else
 212         raise DigestError, "Could not disambiguate residue: #{char}"
 213       end
 214     end
 215
 216     Regexp.new(new_pattern)
 217   end
 218 end
 219
 220
 221 __END__
 222
 223
 224
 225
 226
 227 # Class containing generic sequence methods and nucleic acid and amino acid subclasses.
 228 class Seq < String
 229   # Guess the sequence type by analyzing the first 100 residues allowing for ambiguity codes.
 230   def guess_type
 231     raise ArgumentError, "No sequence." if self.empty?
 232
 233     seq_beg = self[0, 100].upcase
 234
 235     if seq_beg.count( "FLPQIE" ) > 0
 236       Seq::AA.new(self)
 237     elsif seq_beg.count("U") > 0
 238       Seq::NA::RNA.new(self)
 239     else
 240       Seq::NA::DNA.new(self)
 241     end
 242   end
 243
 244   # Method to wrap a sequence to a given width using a given delimiter.
 245   def wrap(width = 80, delimit = $/)
 246     raise ArgumentError, "Cannot wrap sequence to negative width: #{width}." if width <= 0
 247
 248     self.delete!(" \t\n\r")
 249     self.gsub(/.{#{width}}(?!$)/, "\\0#{delimit}")
 250   end
 251
 252   # Method to wrap and replace a sequence to a given width using a given delimiter.
 253   def wrap!(width = 80, delimit = $/)
 254     self.replace(self.wrap(width, delimit))
 255   end
 256
 257   # Method that generates a random sequence of a given length.
 258   def generate(length)
 259     raise ArgumentError, "Cannot generate negative sequence length: #{length}." if length <= 0
 260
 261     alph = self.residues
 262     Array.new(length) { alph[rand(alph.size)] }.join("")
 263   end
 264
 265   # Method that replaces sequence with a randomly generated sequence of a given length.
 266   def generate!(length)
 267     self.replace(self.generate(length))
 268   end
 269
 270   # Class containing methods specific for amino acid (AA) sequences.
 271   class AA < Seq
 272     # Method that returns an array of amino acid residues.
 273     def residues
 274       %w{ F L S Y C W P H Q R I M T N K V A D E G }
 275     end
 276
 277     # Calculate the molecular weight of an amino acid seuqunce.
 278     # The caluculation is only approximate since there is no correction
 279     # for amino bond formation and the MW used are somewhat imprecise:
 280     # http://www.expasy.ch/tools/pscale/Molecularweight.html
 281     def mol_weight
 282       raise ArgumentError, "invalid residues found: #{self.delete("#{residues.join( "" )}")}" if self.upcase =~ /[^#{residues.join( "" )}]/
 283
 284       mol_weight_aa = {
 285         "A" => 89.000,    # Ala
 286         "R" => 174.000,   # Arg
 287         "N" => 132.000,   # Asn
 288         "D" => 133.000,   # Asp
 289         "C" => 121.000,   # Cys
 290         "Q" => 146.000,   # Gln
 291         "E" => 147.000,   # Glu
 292         "G" => 75.000,    # Gly
 293         "H" => 155.000,   # His
 294         "I" => 131.000,   # Ile
 295         "L" => 131.000,   # Leu
 296         "K" => 146.000,   # Lys
 297         "M" => 149.000,   # Met
 298         "F" => 165.000,   # Phe
 299         "P" => 115.000,   # Pro
 300         "S" => 105.000,   # Ser
 301         "T" => 119.000,   # Thr
 302         "W" => 204.000,   # Trp
 303         "Y" => 181.000,   # Tyr
 304         "V" => 117.000,   # Val
 305       }
 306
 307       mw = 0.0
 308
 309       self.upcase.each_char { |c| mw += mol_weight_aa[ c ] }
 310
 311       mw
 312     end
 313   end