code_ruby/lib/maasha/genbank.rb

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