[Biopython-dev] GenBank parser -- first go

Andrew Dalke dalke at acm.org
Wed Dec 6 03:10:36 EST 2000

>I spent this weekend getting together a GenBank parser,
>which I hope is something that we could include in Biopython in the

Wow!  I'm glad that people other than me can use it.  I've been
working on it for so long now that I don't have a good idea of what it
means to come into it from scratch.

>we all definately have to give
>Andrew another big pat on the back for his awesome tool :-).

Thank you.

>5. Can the code be speeded up/improved in any ways? Suggestions to
>help me code better are always very welcome!

I look over it and it seems quite good.  I do have some comments,
which I've included here.  (Most of the points are relevant to
Python and Martel programming, so I've sent it to the list instead
of just you directly.)

You use

cur_record = iterator.next()
while cur_record:
  cur_record = iterator.next()

The standard Python idiom is

while 1:
  cur_record = iterator.next()
  if not cur_record:

Instead of
  indent_space = Martel.RepN(Martel.Str(" "), 2)

it is better to do (note: two spaces)
  indent_space = Martel.Str("  ")

They give the same result, but "  " gets checked with a single
test while two " "s is done with two tests.  This is even more
appropriate with
  qualifier_space = Martel.RepN(Martel.Str(" "), 21)

(Actually, there is a Martel.optimize module which contains the
function 'merge_strings'.  It should merge " "," " into " " but it
isn't automatically used.)

As an aside, I can see I've been focused on regexps compared to
you.  You have
  blank_space = Martel.Rep1(Martel.Str(" "))
where I would use
  blank_space = Martel.Re(" *")
There's no difference in implementation - the decision on which one to
use is based on readability/usability.  I (sadly?) know a lot about
regexps so I probably find them more usable :)

Case in point, I probably would have written the LOCUS definition with
more regexps.

def choice(tag, words):
  exp_list = map(Martel.Str, words)
  return Martel.Group(tag, exp_list)

residue_types = choice("residue_type",
                       ["DNA", "RNA", "mRNA", "PROTEIN"])

data_file_division = choice("data_file_division",
              ["PRI", "ROD", "MAM", "VRT", "INV", "PLN", "BCT", "RNA",
               "VRL", "PHG", "SYN", "UNA", "EST", "PAT", "STS", "GSS",

date = Martel.Group("date",

locus_line = Martel.Re("LOCUS +(?P<locus>\w+) +(?P<size>\d+) bp +") + \
             residue_types + blank_space + data_file_division + \
             blank_space + date + Martel.AnyEol()

Interestingly, I didn't use a regexp for residue_types like you do.
That's because I worry about people looking at a list of strings and
thinking they can add arbitrary text - forgetting about escapeing
regex characters.

Skipping ahead, the 'valid_f_keys' might someday include characters
like '+' and '.', so you really should use a function like the one I
gave above.

Hmmm.  I've been pickier than you about ignoring the leading
whitespace.  For example, you have

# definition line
# DEFINITION  Genomic sequence for Arabidopsis thaliana BAC T25K16 from
#             chromosome I, complete sequence.

definition = Martel.Group("definition",
                          Martel.Rep1(blank_space +

definition_line = Martel.Group("definition_line",
                               Martel.Str("DEFINITION") +

By comparison, I would enforced that the text be folded to start on
column 13 using

definition_line = Martel.Group("definition_line",
             Martel.Str("DEFINITION ") + Martel.ToEol("definition") + \
  Martel.Rep(Martel.Str("           ") + Martel.ToEol("definition")))

Here's a justification for this.  It's already common practice with
GenBank files to have subitems indented under the major item.  For

SOURCE      thale cress.
  ORGANISM  Arabidopsis thaliana

Suppose some day the powers that be add a subitem to the DEFINITION
DEFINITION  Genomic sequence for Arabidopsis thaliana BAC T25K16 from
            chromosome I, complete sequence.
  BLAH      Abcd Ef Ghijkl.

I consider it a good thing for the parser to break at this point,
rather than include the "BLAH      Abcd Ef Ghijkl." text as part of
the definition.

It seems the form of
  "LABEL"    text starts on column 13
             and may fold over multiple

is pretty common.  If that's the case, you can make things simpler by
using a functions to make the definition.

def make_line(label, line_name, data_name):
    assert len(label) < INDENT, "label text too long"
    first_line = Martel.Str(label + " " * (INDENT - len(label))) + \
    other_lines = Martel.Str(" " * INDENT) + \
    return Martel.Group(line_name, first_line + Rep(other_lines))

definition_line = make_line("DEFINITION", "definition_line",

accession_line = make_line("ACCESSION", "accession_line", "accession")

(If you were really trusting you could use just the line label and
string.lower it to get the data name, and add "_line" to that to get
the line name.  I'm usually more explicit than that.)

\w includes \d so you don't need to do [\w\d] for the nid

For that matter, [\d]+ is the same as \d+ (as in the gi)

I'm still undecided about the AtBeginning and AtEnd commands.  (You
use the former in the version_line definition.)  I don't know if they
should test for beginning/end of line or beginning/end of input text.
In fact, I would rather not have them at all, since they don't work
well with the RecordReader idea, where the text is broken up into
parts.  There's no real need for AtBeginning here so you can remove

With the function definition above, you can also replace

  keywords_line = make_line("KEYWORDS", "keywords_line", "keywords")
  source_line = make_line("SOURCE", "source_line", "source")

Out of curiosity, why do arab1.gb and cor6_6.gb have different

SOURCE      thale cress.
  ORGANISM  Arabidopsis thaliana
            Eukaryota; Viridiplantae; Streptophyta; Embryophyta;
            euphyllophytes; Spermatophyta; Magnoliophyta; eudicotyledons;
            Rosidae; Capparales; Brassicaceae; Arabidopsis.

SOURCE      thale cress.
  ORGANISM  Arabidopsis thaliana
            Eukaryota; Viridiplantae; Embryophyta; Tracheophyta;
            Magnoliophyta; eudicotyledons; core eudicots; Rosidae; eurosids
            Brassicales; Brassicaceae; Arabidopsis.

The first has "Streptophyta", "euphyllophytes" and "Capparales".

The second has "core eudicots", "eurosids II" and "Brassicales".

Martel includes a helper function called 'Integer' which can simplify
some of your definitions, as with

reference_num = Martel.Integer("reference_num")
pubmed_id = Martel.Integer("pubmed_id")

Here's another preference of mine.  You have

  sequence = Martel.Group("sequence",
                          Martel.Re("[\w ]+"))

where I would do
  sequence = Martel.ToEol("sequence")

The difference is that you only accept a-zA-Z0-9_ and the space
character.  It doesn't accept "-" or "*", which I can see possibly
getting into the data.  Since there's no need to validate the
characters, might as well just consume anything you find.

(Note to self - need to use \R for the ToEol definition.)

Looking at GenBank.py, you use the record definition thusly:

        parser = genbank_format.record.make_parser(debug_level = debug)


You really should cache the created parser.  It can take quite some
time to generate.  During my PIR tests, about 98% of time of some of
the tests were spent doing generation rather than parsing.

I think you do it because you want to allow different debug levels.
You can still support that in a couple of ways.  Here's one:

class _Scanner:
    def __init__(self):
        self._cached_parsers = {}
    def feed(self, handler, consumer, debug = 0)
        parser = self._cached_parsers.get(debug)
        if parser is None:
            parser = self._cached_parsers[debug] = \
                genbank_format.record.make_parser(debug_level = debug)


Here's another which is tuned for smaller memory use and the
assumption that you almost never change debug levels.

class _Scanner:
    def __init__(self):
        self._cached_parser = None
        self._cached_debug = None
    def feed(self, handler, consumer, debug = 0)
        if self._cached_debug == debug:
            parser = self._cached_parser
            parser = self._cached_parser = \
                genbank_format.record.make_parser(debug_level = debug)
            self._cached_debug = debug


You have a list of tags you are interested in receiving.  Martel has a
function to create a new expression tree based but only sending back
the events you are interested in receiving.

   expression = genbank_format.record
   expression = Martel.select_names(expression, interest_tags)
   parser = expression.make_parser()

BTW, this was implemented because Python's function call overhead is
pretty large so I wanted a way to reduce the number of calls if I knew
an event wasn't needed.

         fun_to_call = eval("self._consumer" + "." + name)
         fun_to_call = getattr("self._consumer",  name)

You should also do
instead of
         apply(fun_to_call, (info_to_pass,))

Here's a cute implementation for _PrintConsumer (untested)

# Let's you define new labels if you want them
event_name_converter = {
  "start_feature_table": "Starting feature table",
  "record_end": "End of Record!",

class _PrintWrapper:
  def __init__(self, name):
    self.name = name
  def __call__(self, content):
    print "%s: %s" % (self.name, content)

class _PrintConsumer:
  def __init__(self):
    self.data = 'blah'
  def __getattr__(self, name):
    if name[:2] == "__":
      raise AttributeError, name
    return _PrintWrapper(event_name_converter.get(name, name))

                    dalke at acm.org

More information about the Biopython-dev mailing list