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
  28 # Error class for all exceptions to do with Genbank.
  29 class GenbankError < StandardError; end
  30
  31 class Genbank < Filesys
  32   def initialize(io)
  33     @io      = io
  34     @entry   = []
  35   end
  36
  37   # Iterator method for parsing Genbank entries.
  38   def each(hash_keys = nil, hash_feats = nil, hash_quals = nil)
  39     while @entry = get_entry do
  40       keys = get_keys(hash_keys)
  41       seq  = get_seq
  42
  43       features = GenbankFeatures.new(@entry, hash_feats, hash_quals)
  44
  45       features.each do |record|
  46         keys.each_pair { |key,val| record[key] = val }
  47
  48                                 loc = Locator.new(record[:LOCATOR], seq)
  49                                 record[:SEQ]     = loc.subseq.seq
  50         record[:SEQ_LEN] = loc.subseq.length
  51                                 record[:STRAND]  = loc.strand
  52                                 record[:S_BEG]   = loc.s_beg
  53                                 record[:S_END]   = loc.s_end
  54
  55         yield record
  56       end
  57     end
  58   end
  59
  60   private
  61
  62   # Method to get the next Genbank entry form an ios and return this.
  63   def get_entry
  64     block = @io.gets("//" + $/)
  65     return nil if block.nil?
  66
  67     block.chomp!("//" + $/ )
  68
  69     entry = block.split $/
  70
  71     return nil if entry.empty?
  72
  73     entry
  74   end
  75
  76   # Method to get the DNA sequence from a Genbank entry and return
  77   # this as a Seq object.
  78   def get_seq
  79     i = @entry.size
  80     j = i - 1
  81
  82     while @entry[j] and @entry[j] !~ /^[A-Z]/
  83       j -= 1
  84     end
  85
  86     seq = @entry[j + 1 .. i].join.delete(" 0123456789")
  87
  88     Seq.new(nil, seq, "dna") if seq
  89   end
  90
  91   # Method to get the base keys from Genbank entry and return these
  92   # in a hash.
  93   def get_keys(hash_keys)
  94     keys = {}
  95     i    = 0
  96     j    = 0
  97
  98     while @entry[i] and @entry[i] !~ /^FEATURES/
  99       if @entry[i] =~ /^\s{0,3}([A-Z]{2})/
 100         if want_key?(hash_keys, $1)
 101           j = i + 1
 102
 103           key, val = @entry[i].lstrip.split(/\s+/, 2)
 104
 105           while @entry[j] and @entry[j] !~ /^\s{0,3}[A-Z]/
 106             val << @entry[j].lstrip
 107             j += 1
 108           end
 109
 110           if keys[key.to_sym]
 111             keys[key.to_sym] << ";" + val
 112           else
 113             keys[key.to_sym] = val
 114           end
 115         end
 116       end
 117
 118       j > i ? i = j : i += 1
 119     end
 120
 121     keys
 122   end
 123
 124   def want_key?(hash_keys, key)
 125     if hash_keys
 126       if hash_keys[key.to_sym]
 127         return true
 128       else
 129         return false
 130       end
 131     else
 132       return true
 133     end
 134   end
 135 end
 136
 137 class GenbankFeatures
 138   @@i = 0
 139   @@j = 0
 140
 141   def initialize(entry, hash_feats, hash_quals)
 142     @entry      = entry
 143     @hash_feats = hash_feats
 144     @hash_quals = hash_quals
 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__