code_ruby/Maasha/lib/genbank.rb

   1 # This program is free software; you can redistribute it and/or
   2 # modify it under the terms of the GNU General Public License
   3 # as published by the Free Software Foundation; either version 2
   4 # of the License, or (at your option) any later version.
   5
   6 # This program is distributed in the hope that it will be useful,
   7 # but WITHOUT ANY WARRANTY; without even the implied warranty of
   8 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   9 # GNU General Public License for more details.
  10
  11 # You should have received a copy of the GNU General Public License
  12 # along with this program; if not, write to the Free Software
  13 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  14
  15 # http://www.gnu.org/copyleft/gpl.html
  16
  17 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  18
  19 # This software is part of the Biopieces framework (www.biopieces.org).
  20
  21 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  22
  23 require 'seq'
  24 require 'zlib'
  25
  26 # Error class for all exceptions to do with Genbank.
  27 class GenbankError < StandardError; end
  28
  29 class Genbank
  30   include Enumerable
  31
  32   # Class method allowing open to be used on (zipped) files.
  33   # See File.open.
  34   def self.open(*args)
  35     ios = self.zopen(*args)
  36
  37     if block_given?
  38       begin
  39         yield ios
  40       ensure
  41         ios.close
  42       end
  43     else
  44       return ios
  45     end
  46   end
  47
  48   def initialize(io)
  49     @io      = io
  50     @entry   = []
  51   end
  52
  53   # Method to close ios.
  54   def close
  55     @io.close
  56   end
  57
  58   # Iterator method for parsing Genbank entries.
  59   def each(hash_keys, hash_feats, hash_quals)
  60     while @entry = get_entry do
  61       keys = get_keys(hash_keys)
  62       seq  = get_seq
  63
  64       features = GenbankFeatures.new(@entry, hash_feats, hash_quals)
  65
  66       features.each do |record|
  67         keys.each_pair { |key,val| record[key] = val }
  68         loc = Locator.new(record[:LOCATOR], seq)
  69         record[:SEQ]    = loc.subseq.seq
  70         record[:STRAND] = loc.strand
  71         record[:S_BEG]  = loc.s_beg
  72         record[:S_END]  = loc.s_end
  73
  74         yield record
  75       end
  76     end
  77   end
  78
  79   private
  80
  81   # Helper method to return an ios to a file that may be zipped in which case
  82   # the ios is unzipped on the fly. See File.open.
  83   def self.zopen(*args)
  84     ios = File.open(*args)
  85
  86     begin
  87       ios = Zlib::GzipReader.new(ios)
  88     rescue
  89       ios.rewind
  90     end
  91
  92     self.new(ios)
  93   end
  94
  95   # Method to get the next Genbank entry form an ios and return this.
  96   def get_entry
  97     block = @io.gets("//" + $/)
  98     return nil if block.nil?
  99
 100     block.chomp!("//" + $/ )
 101
 102     entry = block.split $/
 103
 104     return nil if entry.empty?
 105
 106     entry
 107   end
 108
 109   # Method to get the DNA sequence from a Genbank entry and return
 110   # this as a Seq object.
 111   def get_seq
 112     seq = Seq.new
 113     i   = @entry.size
 114
 115     while @entry[i] !~ /^ORIGIN/
 116       i -= 1
 117     end
 118
 119     seq.seq  = @entry[i + 1 .. @entry.size].join.delete( " 0123456789")
 120     seq.type = "dna"
 121     seq
 122   end
 123
 124   # Method to get the base keys from Genbank entry and return these
 125   # in a hash.
 126   def get_keys(hash_keys)
 127     keys = {}
 128     i    = 0
 129     j    = 0
 130
 131     while @entry[i] and @entry[i] !~ /^FEATURES/
 132       if @entry[i] =~ /^\s{0,3}([A-Z]{2})/
 133         if want_key?(hash_keys, $1)
 134           j = i + 1
 135
 136           key, val = @entry[i].lstrip.split(/\s+/, 2)
 137
 138           while @entry[j] and @entry[j] !~ /^\s{0,3}[A-Z]/
 139             val << @entry[j].lstrip
 140             j += 1
 141           end
 142
 143           if keys[key.to_sym]
 144             keys[key.to_sym] << ";" + val
 145           else
 146             keys[key.to_sym] = val
 147           end
 148         end
 149       end
 150
 151       j > i ? i = j : i += 1
 152     end
 153
 154     keys
 155   end
 156
 157   def want_key?(hash_keys, key)
 158     if hash_keys
 159       if hash_keys[key.to_sym]
 160         return true
 161       else
 162         return false
 163       end
 164     else
 165       return true
 166     end
 167   end
 168 end
 169
 170 class GenbankFeatures
 171   @@i = 0
 172   @@j = 0
 173
 174   def initialize(entry, hash_feats, hash_quals)
 175     @entry      = entry
 176     @hash_feats = hash_feats
 177     @hash_quals = hash_quals
 178   end
 179
 180   def each
 181     while @entry[@@i] and @entry[@@i] !~ /^ORIGIN/
 182       if @entry[@@i] =~ /^\s{5}([A-Za-z_-]+)/
 183         if want_feat? $1
 184           record = {}
 185
 186           feat, loc = @entry[@@i].lstrip.split(/\s+/, 2)
 187
 188           @@j = @@i + 1
 189
 190           while @entry[@@j] and @entry[@@j] !~ /^(\s{21}\/|\s{5}[A-Za-z_-]|[A-Z])/
 191             loc << @entry[@@j].lstrip
 192             @@j += 1
 193           end
 194
 195           get_quals.each_pair { |k,v|
 196             record[k.upcase.to_sym] = v
 197           }
 198
 199           record[:FEATURE] = feat
 200           record[:LOCATOR] = loc
 201
 202           yield record
 203         end
 204       end
 205
 206       @@j > @@i ? @@i = @@j : @@i += 1
 207     end
 208   end
 209
 210   private
 211
 212   def get_quals
 213     quals = {}
 214     k     = 0
 215
 216     while @entry[@@j] and @entry[@@j] !~ /^\s{5}[A-Za-z_-]|^[A-Z]/
 217       if @entry[@@j] =~ /^\s{21}\/([^=]+)="([^"]+)/
 218         qual = $1
 219         val  = $2
 220
 221         if want_qual? qual
 222           k = @@j + 1
 223
 224           while @entry[k] and @entry[k] !~ /^(\s{21}\/|\s{5}[A-Za-z_-]|[A-Z])/
 225             val << @entry[k].lstrip.chomp('"')
 226             k += 1
 227           end
 228
 229           if quals[qual]
 230             quals[qual] << ";" + val
 231           else
 232             quals[qual] = val
 233           end
 234         end
 235       end
 236
 237       k > @@j ? @@j = k : @@j += 1
 238     end
 239
 240     quals
 241   end
 242
 243   def want_feat?(feat)
 244     if @hash_feats
 245       if @hash_feats[feat.upcase.to_sym]
 246         return true
 247       else
 248         return false
 249       end
 250     else
 251       return true
 252     end
 253   end
 254
 255   def want_qual?(qual)
 256     if @hash_quals
 257       if @hash_quals[qual.upcase.to_sym]
 258         return true
 259       else
 260         return false
 261       end
 262     else
 263       return true
 264     end
 265   end
 266 end
 267
 268
 269 # Error class for all exceptions to do with Genbank/EMBL/DDBJ feature table locators.
 270 class LocatorError < StandardError; end
 271
 272 class Locator
 273   attr_accessor :locator, :seq, :subseq
 274
 275   def initialize(locator, seq)
 276     @loc_orig = locator
 277     @locator  = locator
 278     @seq      = seq
 279     @subseq   = Seq.new(nil, "", "dna")
 280   end
 281
 282   def subseq
 283     parse_locator
 284   end
 285
 286   def strand
 287     if @loc_orig.match("complement")
 288       return "-"
 289     else
 290       return "+"
 291     end
 292   end
 293
 294   def s_beg
 295     if @loc_orig =~ /(\d+)/
 296       return $1.to_i - 1
 297     end
 298   end
 299
 300   def s_end
 301     if @loc_orig.reverse =~ /(\d+)/
 302       return $1.reverse.to_i - 1
 303     end
 304   end
 305
 306   private
 307
 308   # Method that uses recursion to parse a locator string from a feature
 309   # table and fetches the appropriate subsequence. the operators
 310   # join(), complement(), and order() are handled.
 311   # the locator string is broken into a comma separated lists, and
 312   # modified if the params donnot balance. otherwise the comma separated
 313   # list of ranges are stripped from operators, and the subsequence are
 314   # fetched and handled according to the operators.
 315   # SNP locators are also dealt with (single positions).
 316   def parse_locator(join = nil, comp = nil, order = nil)
 317     intervals = @locator.split(",")
 318
 319     unless balance_params?(intervals.first)   # locator includes a join/comp/order of several ranges
 320       case @locator
 321       when /^join\((.*)\)$/
 322         @locator = $1
 323         join     = true
 324       when /^complement\((.*)\)$/
 325         @locator = $1
 326         comp     = true
 327       when /^order\((.*)\)$/
 328         @locator = $1
 329         order    = true
 330       end
 331
 332       parse_locator(join, comp, order)
 333     else
 334       intervals.each do |interval|
 335         case interval
 336         when /^join\((.*)\)$/
 337           @locator = $1
 338           join     = true
 339           parse_locator(join, comp, order)
 340         when /^complement\((.*)\)$/
 341           @locator = $1
 342           comp     = true
 343           parse_locator(join, comp, order)
 344         when /^order\((.*)\)$/
 345           @locator = $1
 346           order    = true
 347           parse_locator(join, comp, order)
 348         when /^[<>]?(\d+)[^\d]+(\d+)$/
 349           int_beg = $1.to_i - 1
 350           int_end = $2.to_i - 1
 351
 352           newseq = Seq.new(nil, @seq.seq[int_beg...int_end], "dna")
 353           newseq.revcomp if comp
 354
 355           @subseq.seq << (order ? " " + newseq.seq : newseq.seq)
 356         when /^(\d+)$/
 357           pos = $1.to_i - 1
 358
 359           newseq = Seq.new(nil, @seq.seq[pos], "dna")
 360           newseq.revcomp if comp
 361
 362           @subseq.seq << (order ? " " + newseq.seq : newseq.seq)
 363         else
 364           $stderr.puts "WARNING: Could not match locator -> #{locator}";
 365           @subseq.seq << ""
 366         end
 367       end
 368     end
 369
 370     return @subseq
 371   end
 372
 373   def balance_params?(locator)
 374     parens = 0
 375
 376     locator.each_char do |char|
 377       case char
 378       when '(' then parens += 1
 379       when ')' then parens -= 1
 380       end
 381     end
 382
 383     if parens == 0
 384       return true
 385     else
 386       return false
 387     end
 388   end
 389 end
 390
 391
 392 __END__