code_ruby/Maasha/lib/seq.rb

   1 # Class containing generic sequence methods and nucleic acid and amino acid subclasses.
   2 class Seq < String
   3         attr_accessor :seq, :seq_type
   4
   5         # Method to initialize a new sequence.
   6         def initialize( seq = "", seq_type = nil )
   7                 @seq      = seq
   8                 @seq_type = seq_type
   9         end
  10
  11         # Guess the sequence type by analyzing the first 100 residues allowing for ambiguity codes.
  12         def seq_type_guess
  13                 seq_beg = @seq[ 0, 100 ].upcase
  14
  15                 if seq_beg.count( "FLPQIE" ) > 0
  16                         :AA
  17                 elsif seq_beg.count( "U" ) > 0
  18                         :RNA
  19                 else
  20                         :DNA
  21                 end
  22         end
  23
  24         # Guess and replace the sequence type by analyzing the first 100 residues allowing for ambiguity codes.
  25         def seq_type_guess!
  26                 @seq_type = seq_type_guess
  27         end
  28
  29         # Method that return an array of the residue alphabet for a given sequence type.
  30         def seq_alph( seq_type )
  31                 hash = {
  32                         :DNA => %w{ A T C G },
  33                         :RNA => %w{ A U C G },
  34                         :AA  => %w{ F L S Y C W P H Q R I M T N K V A D E G },
  35                 }
  36
  37                 raise "ERROR: Sequence type '#{ seq_type }' not recognized." unless hash.include?( seq_type )
  38                 return hash[ seq_type ]
  39         end
  40
  41         # Method to wrap a sequence to a given width using a given delimiter.
  42         def wrap( width = 80, delimit = "\n" )
  43                 raise "ERROR: Wrap width must be an integer." unless width.is_a? Fixnum
  44                 raise "ERROR: Cannot wrap sequence to negative width: #{ width }." if width <= 0
  45                 @seq.tr!( " \t\n\r", '' )
  46                 @seq.gsub( /.{#{ width }}/, "\\0#{ delimit }" ).sub( /#{ delimit }$/, "" )
  47         end
  48
  49         # Method to wrap and replace a sequence to a given width using a given delimiter.
  50         def wrap!( width = 80, delimit = "\n" )
  51                 @seq = wrap( width, delimit )
  52         end
  53
  54         # Method that generates a random sequence of a given length.
  55         def generate( length )
  56                 raise "ERROR: Length must be an integer." unless length.is_a? Fixnum
  57                 raise "ERROR: Cannot generate negative sequence length: #{ length }." if length <= 0
  58
  59                 alph = seq_alph( @seq_type )
  60                 seq  = Array.new( length ) { alph[ rand( alph.size ) ] }.join
  61         end
  62
  63         # Method that replaces sequence with a randomly generated sequence of a given length.
  64         def generate!( length )
  65                 @seq = generate( length )
  66         end
  67
  68         # Class containing methods specific for amino acid (AA) sequences.
  69         class AA < Seq
  70                 # Method to initialize a new amino acid sequence.
  71                 def initialize( seq = "" )
  72                         @seq      = seq
  73                         @seq_type = :AA
  74                 end
  75
  76                 # Calculate the molecular weight of an amino acid seuqunce.
  77                 # The caluculation is only approximate since there is no correction
  78                 # for amino bond formation and the MW used are somewhat imprecise:
  79                 # http://www.expasy.ch/tools/pscale/Molecularweight.html
  80                 def mol_weight
  81                         mol_weight_aa = {
  82                                 "A" => 89.000,    # Ala
  83                                 "R" => 174.000,   # Arg
  84                                 "N" => 132.000,   # Asn
  85                                 "D" => 133.000,   # Asp
  86                                 "C" => 121.000,   # Cys
  87                                 "Q" => 146.000,   # Gln
  88                                 "E" => 147.000,   # Glu
  89                                 "G" => 75.000,    # Gly
  90                                 "H" => 155.000,   # His
  91                                 "I" => 131.000,   # Ile
  92                                 "L" => 131.000,   # Leu
  93                                 "K" => 146.000,   # Lys
  94                                 "M" => 149.000,   # Met
  95                                 "F" => 165.000,   # Phe
  96                                 "P" => 115.000,   # Pro
  97                                 "S" => 105.000,   # Ser
  98                                 "T" => 119.000,   # Thr
  99                                 "W" => 204.000,   # Trp
 100                                 "Y" => 181.000,   # Tyr
 101                                 "V" => 117.000,   # Val
 102                         }
 103
 104                         mw = 0.0
 105
 106                         @seq.upcase.each_char do |c|
 107                                 raise "ERROR: Unknown amino acid: #{ c }" unless mol_weight_aa.include?( c )
 108                                 mw += mol_weight_aa[ c ]
 109                         end
 110
 111                         mw
 112                 end
 113         end
 114
 115         # Class containing methods specific for nucleic acid (NA) sequences.
 116         class NA < Seq
 117                 # Class containing methods specific for DNA sequences.
 118                 class DNA < NA
 119                         # Method that complements DNA sequence including ambiguity codes.
 120                         def complement
 121                                 @seq.tr!( 'AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'TCGAAYRWSKMDHBVNtcgaayrwskmdhbvn' )
 122                         end
 123                 end
 124
 125                 # Class containing methods specific for RNA sequences.
 126                 class RNA < NA
 127                         # Method that complements RNA sequence including ambiguity codes.
 128                         def complement
 129                                 @seq.tr!( 'AGCUTRYWSMKHDVBNagcutrywsmkhdvbn', 'UCGAAYRWSKMDHBVNucgaayrwskmdhbvn' )
 130                         end
 131                 end
 132         end
 133 end