# Copyright 2000 by Jeffrey Chang. All rights reserved.
# This code is part of the Biopython distribution and governed by its
# license. Please see the LICENSE file that should have been included
# as part of this package.
"""Code to work with GenBank -- http://www.ncbi.nlm.nih.gov/
Classes:
Iterator Iterate through a file of GenBank entries
Dictionary Access a GenBank file using a dictionary interface.
ErrorFeatureParser Catch errors caused during parsing.
FeatureParser Parse GenBank data in Seq and SeqFeature objects.
RecordParser Parse GenBank data into a Record object.
NCBIDictionary Access GenBank using a dictionary interface.
_BaseGenBankConsumer A base class for GenBank consumer that implements
some helpful functions that are in common between
consumers.
_FeatureConsumer Create SeqFeature objects from info generated by
the Scanner
_RecordConsumer Create a GenBank record object from Scanner info.
_PrintingConsumer A debugging consumer.
_Scanner Set up a Martel based GenBank parser to parse a record.
ParserFailureError Exception indicating a failure in the parser (ie.
scanner or consumer)
LocationParserError Exception indiciating a problem with the spark based
location parser.
Functions:
index_file Get a GenBank file ready to be used as a Dictionary.
search_for Do a query against GenBank.
download_many Download many GenBank records.
"""
__all__ = [
'LocationParser',
'Record',
'genbank_format',
]
# standard library
import string
import os
import re
import sgmllib
import urlparse
# XML from python 2.0
from xml.sax import handler
# Martel
import Martel
from Martel import RecordReader
# other Biopython stuff
from Bio.SeqRecord import SeqRecord
from Bio import Alphabet
from Bio.Alphabet import IUPAC
from Bio.Seq import Seq
from Bio import File
from Bio import Index
from Bio.ParserSupport import AbstractConsumer
from Bio.WWW import NCBI
from Bio.WWW import RequestLimiter
from Bio.ParserSupport import EventGenerator
import genbank_format
import Record
from Bio.SeqFeature import Reference
from Bio import SeqFeature
from Bio.GenBank import LocationParser
class Dictionary:
"""Allow a GenBank file to be accessed using a dictionary interface.
"""
__filename_key = '__filename'
def __init__(self, index_file, parser = None):
"""Initialize and open up a GenBank dictionary.
Arguments:
o index_file - The name of the file that servers as the dictionary
index. You need to use the index_file function in this module
to creat this index.
o parser - An optional argument specifying a parser object that
the records should be run through before returning the output. If
parser is None then the unprocessed contents of the file will be
returned.
"""
self._index = Index.Index(index_file)
self._handle = open(self._index[Dictionary.__filename_key])
self._parser = parser
def __len__(self):
return len(self._index)
def __getitem__(self, key):
"""Retrieve an item from the dictionary.
"""
print "keys:", self._index.keys()
# get the location of the record of interest in the file
start, len = self._index[key]
print "start:", start, "len:", len
# read through and get the data from the file
self._handle.seek(start)
data = self._handle.read(len)
print "data:", data
# run the data through the parser if one is specified
if self._parser is not None:
return self._parser.parse(File.StringHandle(data))
return data
def __getattr__(self, name):
return getattr(self._index, name)
def keys(self):
"""Provide valid keys from the index.
If keys is just called on the index (using getattr) then
the index will return internal values such as '__filename' which
we just don't want to see. This just strips these values out
before returning.
"""
all_keys = self._index.keys()
keys_to_remove = []
for key in all_keys:
# if the key is internal (has a __) then add it to the remove list
if key[:2] == '__':
keys_to_remove.append(key)
for key in keys_to_remove:
all_keys.remove(key)
return all_keys
class Iterator:
"""Iterator interface to move over a file of GenBank entries one at a time.
"""
def __init__(self, handle, parser = None, has_header = 0):
"""Initialize the iterator.
Arguments:
o handle - A handle with GenBank entries to iterate through.
o parser - An optional parser to pass the entries through before
returning them. If None, then the raw entry will be returned.
o has_header - Whether or not the file to iterate over has one of
those GenBank headers (ie. if you downloaded it directly from
GenBank). If so, we'll iterate over the header to get past it, and
then the iterator will be set up to return the first record in
the file.
"""
if has_header:
first_line = handle.readline()
assert first_line.find("Genetic Sequence Data Bank") >= 0, \
"Doesn't seem to have a GenBank header."
while 1:
cur_line = handle.readline()
if cur_line.find("reported sequences") >= 0:
break
# read off two more lines and we are ready to go
handle.readline()
handle.readline()
self._reader = RecordReader.StartsWith(handle, "LOCUS")
self._parser = parser
def next(self):
"""Return the next GenBank record from the handle.
Will return None if we ran out of records.
"""
data = self._reader.next()
if self._parser is not None:
if data:
return self._parser.parse(File.StringHandle(data))
return data
class ParserFailureError(Exception):
"""Failure caused by some kind of problem in the parser.
"""
pass
class LocationParserError(Exception):
"""Could not Properly parse out a location from a GenBank file.
"""
pass
class ErrorParser:
"""Parse GenBank files and attempt to catch errors.
This is just a small wrapper class around a passed parser which
catches errors that occur in the parser.
When errors occur, we'll raise a ParserFailureError, which returns
(hopefully helpful) information about where the error occured
"""
def __init__(self, parser, bad_file_handle = None):
"""Initialize an ErrorFeatureParser.
Arguments:
o parser - The actual parser to use in parsing the records.
o bad_file_handle - A handle to write problem GenBank files to
If None, the files will not be saved.
"""
self._bad_file_handle = bad_file_handle
self._parser = parser
def parse(self, handle):
"""Parse the specified handle.
"""
record = handle.read()
try:
return self._parser.parse(File.StringHandle(record))
# deal with any problems as Parser Errors
# XXX Watch out for changes in exceptions, as noted in Andrew's
# comments in the Martel code.
except (Martel.Parser.ParserException,
Martel.Parser.ParserPositionException):
if self._bad_file_handle:
self._bad_file_handle.write(record)
raise ParserFailureError("Could not parse record %s" %
self._parser._consumer.data.id)
# SystemExits often come from location parsing with spark
except LocationParserError, msg:
raise ParserFailureError("Could not parse location %s in record %s"
% (msg, self._parser._consumer.data.id))
class FeatureParser:
"""Parse GenBank files into Seq + Feature objects.
"""
def __init__(self, debug_level = 0, use_fuzziness = 1):
"""Initialize a GenBank parser and Feature consumer.
Arguments:
o debug_level - An optional argument that species the amount of
debugging information Martel should spit out. By default we have
no debugging info (the fastest way to do things), but if you want
you can set this as high as two and see exactly where a parse fails.
o use_fuzziness - Specify whether or not to use fuzzy representations.
The default is 1 (use fuzziness).
"""
self._scanner = _Scanner(debug_level)
self.use_fuzziness = use_fuzziness
def parse(self, handle):
"""Parse the specified handle.
"""
self._consumer = _FeatureConsumer(self.use_fuzziness)
self._scanner.feed(handle, self._consumer)
return self._consumer.data
class RecordParser:
"""Parse GenBank files into Record objects
"""
def __init__(self, debug_level = 0):
"""Initialize the parser.
Arguments:
o debug_level - An optional argument that species the amount of
debugging information Martel should spit out. By default we have
no debugging info (the fastest way to do things), but if you want
you can set this as high as two and see exactly where a parse fails.
"""
self._scanner = _Scanner(debug_level)
def parse(self, handle):
"""Parse the specified handle into a GenBank record.
"""
self._consumer = _RecordConsumer()
self._scanner.feed(handle, self._consumer)
return self._consumer.data
class _BaseGenBankConsumer(AbstractConsumer):
"""Abstract GenBank consumer providing useful general functions.
This just helps to eliminate some duplication in things that most
GenBank consumers want to do.
"""
# Special keys in GenBank records that we should remove spaces from
# For instance, \translation keys have values which are proteins and
# should have spaces and newlines removed from them. This class
# attribute gives us more control over specific formatting problems.
remove_space_keys = ["translation"]
def __init__(self):
pass
def _split_keywords(self, keyword_string):
"""Split a string of keywords into a nice clean list.
"""
# process the keywords into a python list
if keyword_string[-1] == '.':
keywords = keyword_string[:-1]
else:
keywords = keyword_string
keyword_list = string.split(keywords, ';')
clean_keyword_list = map(string.strip, keyword_list)
return clean_keyword_list
def _split_accessions(self, accession_string):
"""Split a string of accession numbers into a list.
"""
# first replace all line feeds with spaces
accession = string.replace(accession_string, "\n", " ")
accession_list = string.split(accession, ' ')
clean_accession_list = map(string.strip, accession_list)
return clean_accession_list
def _split_taxonomy(self, taxonomy_string):
"""Split a string with taxonomy info into a list.
"""
if taxonomy_string[-1] == '.':
tax_info = taxonomy_string[:-1]
else:
tax_info = taxonomy_string
tax_list = string.split(tax_info, ';')
new_tax_list = []
for tax_item in tax_list:
new_items = tax_item.split("\n")
new_tax_list.extend(new_items)
while '' in new_tax_list:
new_tax_list.remove('')
clean_tax_list = map(string.strip, new_tax_list)
return clean_tax_list
def _clean_location(self, location_string):
"""Clean whitespace out of a location string.
The location parser isn't a fan of whitespace, so we clean it out
before feeding it into the parser.
"""
location_line = location_string
for ws in string.whitespace:
location_line = string.replace(location_line, ws, '')
return location_line
def _remove_newlines(self, text):
"""Remove any newlines in the passed text, returning the new string.
"""
# get rid of newlines in the qualifier value
newlines = ["\n", "\r"]
for ws in newlines:
text = text.replace(ws, "")
return text
def _normalize_spaces(self, text):
"""Replace multiple spaces in the passed text with single spaces.
"""
# get rid of excessive spaces
text_parts = text.split(" ")
while '' in text_parts:
text_parts.remove('')
return string.join(text_parts)
def _remove_spaces(self, text):
"""Remove all spaces from the passed text.
"""
return text.replace(" ", "")
class _FeatureConsumer(_BaseGenBankConsumer):
"""Create a SeqRecord object with Features to return.
Attributes:
o use_fuzziness - specify whether or not to parse with fuzziness in
feature locations.
"""
def __init__(self, use_fuzziness):
_BaseGenBankConsumer.__init__(self)
self.data = SeqRecord(None, id = None)
self._use_fuzziness = use_fuzziness
self._seq_type = ''
self._seq_data = ''
self._current_ref = None
self._cur_feature = None
self._cur_qualifier_key = None
self._cur_qualifier_value = None
def locus(self, locus_name):
"""Set the locus name is set as the name of the Sequence.
"""
self.data.name = locus_name
def size(self, content):
pass
def residue_type(self, type):
"""Record the sequence type so we can choose an appropriate alphabet.
"""
self._seq_type = type
def data_file_division(self, division):
self.data.annotations['data_file_division'] = division
def date(self, submit_date):
self.data.annotations['date'] = submit_date
def definition(self, definition):
"""Set the definition as the description of the sequence.
"""
self.data.description = definition
def accession(self, acc_num):
"""Set the accession number as the id of the sequence.
If we have multiple accession numbers, the first one passed is
used.
"""
new_acc_nums = self._split_accessions(acc_num)
# if we haven't set the id information yet, add the first acc num
if self.data.id is None:
if len(new_acc_nums) > 0:
self.data.id = new_acc_nums[0]
def nid(self, content):
self.data.annotations['nid'] = content
def version(self, version_id):
"""Set the version to overwrite the id.
Since the verison provides the same information as the accession
number, plus some extra info, we set this as the id if we have
a version.
"""
self.data.id = version_id
def gi(self, content):
self.data.annotations['gi'] = content
def keywords(self, content):
self.data.annotations['keywords'] = self._split_keywords(content)
def segment(self, content):
self.data.annotations['segment'] = content
def source(self, content):
if content[-1] == '.':
source_info = content[:-1]
else:
source_info = content
self.data.annotations['source'] = source_info
def organism(self, content):
self.data.annotations['organism'] = content
def taxonomy(self, content):
self.data.annotations['taxonomy'] = self._split_taxonomy(content)
def reference_num(self, content):
"""Signal the beginning of a new reference object.
"""
# if we have a current reference that hasn't been added to
# the list of references, add it.
if self._current_ref is not None:
self.data.annotations['references'].append(self._current_ref)
else:
self.data.annotations['references'] = []
self._current_ref = Reference()
def reference_bases(self, content):
"""Attempt to determine the sequence region the reference entails.
Possible types of information we may have to deal with:
(bases 1 to 86436)
(sites)
(bases 1 to 105654; 110423 to 111122)
"""
# first remove the parentheses
ref_base_info = content[1:-1]
all_locations = []
# only attempt to get out information if we find the words
# 'bases' and 'to'
if (string.find(ref_base_info, 'bases') != -1 and
string.find(ref_base_info, 'to') != -1):
# get rid of the beginning 'bases'
ref_base_info = ref_base_info[5:]
# split possibly multiple locations using the ';'
all_base_info = string.split(ref_base_info, ';')
for base_info in all_base_info:
start, end = string.split(base_info, 'to')
this_location = \
SeqFeature.FeatureLocation(int(string.strip(start)),
int(string.strip(end)))
all_locations.append(this_location)
# make sure if we are not finding information then we have
# the string 'sites' or the string 'bases'
elif (ref_base_info == 'sites' or
string.strip(ref_base_info) == 'bases'):
pass
# otherwise raise an error
else:
raise ValueError("Could not parse base info %s in record %s" %
(ref_base_info, self.data.id))
self._current_ref.location = all_locations
def authors(self, content):
self._current_ref.authors = content
def title(self, content):
self._current_ref.title = content
def journal(self, content):
self._current_ref.journal = content
def medline_id(self, content):
self._current_ref.medline_id = content
def pubmed_id(self, content):
self._current_ref.pubmed_id = content
def remark(self, content):
self._current_ref.comment = content
def comment(self, content):
self.data.annotations['comment'] = string.join(content, "\n")
def features_line(self, content):
"""Get ready for the feature table when we reach the FEATURE line.
"""
self.start_feature_table()
def start_feature_table(self):
"""Indicate we've got to the start of the feature table.
"""
# make sure we've added on our last reference object
if self._current_ref is not None:
self.data.annotations['references'].append(self._current_ref)
self._current_ref = None
def _add_feature(self):
"""Utility function to add a feature to the SeqRecord.
This does all of the appropriate checking to make sure we haven't
left any info behind, and that we are only adding info if it
exists.
"""
if self._cur_feature:
# if we have a left over qualifier, add it to the qualifiers
# on the current feature
self._add_qualifier()
self._cur_qualifier_key = ''
self._cur_qualifier_value = ''
self.data.features.append(self._cur_feature)
def feature_key(self, content):
# if we already have a feature, add it on
self._add_feature()
# start a new feature
self._cur_feature = SeqFeature.SeqFeature()
self._cur_feature.type = content
# assume positive strand to start with if we have DNA. The
# complement in the location will change this later.
if self._seq_type == "DNA":
self._cur_feature.strand = 1
def location(self, content):
"""Parse out location information from the location string.
This uses Andrew's nice spark based parser to do the parsing
for us, and translates the results of the parse into appropriate
Location objects.
"""
# --- first preprocess the location for the spark parser
# we need to clean up newlines and other whitespace inside
# the location before feeding it to the parser.
# locations should have no whitespace whatsoever based on the
# grammer
location_line = self._clean_location(content)
# Older records have junk like replace(266,"c") in the
# location line. Newer records just replace this with
# the number 266 and have the information in a more reasonable
# place. So we'll just grab out the number and feed this to the
# parser. We shouldn't really be losing any info this way.
if string.find(location_line, 'replace') != -1:
comma_pos = string.find(location_line, ',')
location_line = location_line[8:comma_pos]
# feed everything into the scanner and parser
try:
parse_info = \
LocationParser.parse(LocationParser.scan(location_line))
# spark raises SystemExit errors when parsing fails
except SystemExit:
raise LocationParserError(location_line)
# print "parse_info:", repr(parse_info)
# now add the information we get out of the parser to the feature
self._set_location_info(parse_info, self._cur_feature)
def _set_function(self, function, cur_feature):
"""Set the location information based on a function.
This handles all of the location functions like 'join', 'complement'
and 'order'.
Arguments:
o function - A LocationParser.Function object specifying the
function we are acting on.
o cur_feature - The feature to add information to.
"""
assert isinstance(function, LocationParser.Function), \
"Expected a Function object, got %s" % function
if function.name == "complement":
# mark the current feature as being on the opposite strand
cur_feature.strand = -1
# recursively deal with whatever is left inside the complement
for inner_info in function.args:
self._set_location_info(inner_info, cur_feature)
# deal with functions that have multipe internal segments that
# are connected somehow.
# join and order are current documented functions.
# one-of is something I ran across in old files. Treating it
# as a sub sequence feature seems appropriate to me.
elif (function.name == "join" or function.name == "order" or
function.name == "one-of"):
self._set_ordering_info(function, cur_feature)
else:
raise ValueError("Unexpected function name: %s" % function.name)
def _set_ordering_info(self, function, cur_feature):
"""Parse a join or order and all of the information in it.
This deals with functions that order a bunch of locations,
specifically 'join' and 'order'. The inner locations are
added as subfeatures of the top level feature
"""
# for each inner element, create a sub SeqFeature within the
# current feature, then get the information for this feature
for inner_element in function.args:
new_sub_feature = SeqFeature.SeqFeature()
# add _join or _order to the name to make the type clear
new_sub_feature.type = cur_feature.type + '_' + function.name
# inherit references and strand from the parent feature
new_sub_feature.ref = cur_feature.ref
new_sub_feature.ref_db = cur_feature.ref_db
new_sub_feature.strand = cur_feature.strand
# set the information for the inner element
self._set_location_info(inner_element, new_sub_feature)
# now add the feature to the sub_features
cur_feature.sub_features.append(new_sub_feature)
# set the location of the top -- this should be a combination of
# the start position of the first sub_feature and the end position
# of the last sub_feature
feature_start = cur_feature.sub_features[0].location.start
feature_end = cur_feature.sub_features[-1].location.end
cur_feature.location = SeqFeature.FeatureLocation(feature_start,
feature_end)
def _set_location_info(self, parse_info, cur_feature):
"""Set the location information for a feature from the parse info.
Arguments:
o parse_info - The classes generated by the LocationParser.
o cur_feature - The feature to add the information to.
"""
# base case -- we are out of information
if parse_info is None:
return
# parse a location -- this is another base_case -- we assume
# we have no information after a single location
elif isinstance(parse_info, LocationParser.AbsoluteLocation):
self._set_location(parse_info, cur_feature)
return
# parse any of the functions (join, complement, etc)
elif isinstance(parse_info, LocationParser.Function):
self._set_function(parse_info, cur_feature)
# otherwise we are stuck and should raise an error
else:
raise ValueError("Could not parse location info: %s"
% parse_info)
def _set_location(self, location, cur_feature):
"""Set the location information for a feature.
Arguments:
o location - An AbsoluteLocation object specifying the info
about the location.
o cur_feature - The feature to add the information to.
"""
# check to see if we have a cross reference to another accession
# ie. U05344.1:514..741
if location.path is not None:
cur_feature.ref = location.path.accession
cur_feature.ref_db = location.path.database
# now get the actual location information
cur_feature.location = self._get_location(location.local_location)
def _get_location(self, range_info):
"""Return a (possibly fuzzy) location from a Range object.
Arguments:
o range_info - A location range (ie. something like 67..100). This
may also be a single position (ie 27).
This returns a FeatureLocation object.
If parser.use_fuzziness is set at one, the positions for the
end points will possibly be fuzzy.
"""
# check if we just have a single base
if not(isinstance(range_info, LocationParser.Range)):
pos = self._get_position(range_info)
return SeqFeature.FeatureLocation(pos, pos)
# otherwise we need to get both sides of the range
else:
# get *Position objects for the start and end
start_pos = self._get_position(range_info.low)
end_pos = self._get_position(range_info.high)
return SeqFeature.FeatureLocation(start_pos, end_pos)
def _get_position(self, position):
"""Return a (possibly fuzzy) position for a single coordinate.
Arguments:
o position - This is a LocationParser.* object that specifies
a single coordinate. We will examine the object to determine
the fuzziness of the position.
This is used with _get_location to parse out a location of any
end_point of arbitrary fuzziness.
"""
# case 1 -- just a normal number
if (isinstance(position, LocationParser.Integer)):
final_pos = SeqFeature.ExactPosition(position.val)
# case 2 -- we've got a > sign
elif isinstance(position, LocationParser.LowBound):
final_pos = SeqFeature.AfterPosition(position.base.val)
# case 3 -- we've got a < sign
elif isinstance(position, LocationParser.HighBound):
final_pos = SeqFeature.BeforePosition(position.base.val)
# case 4 -- we've got 100^101
elif isinstance(position, LocationParser.Between):
final_pos = SeqFeature.BetweenPosition(position.low.val,
position.high.val)
# case 5 -- we've got (100.101)
elif isinstance(position, LocationParser.TwoBound):
final_pos = SeqFeature.WithinPosition(position.low.val,
position.high.val)
# if it is none of these cases we've got a problem!
else:
raise ValueError("Unexpected LocationParser object %r" %
position)
# if we are using fuzziness return what we've got
if self._use_fuzziness:
return final_pos
# otherwise return an ExactPosition equivalent
else:
return SeqFeature.ExactPosition(final_pos.location)
def _add_qualifier(self):
"""Add a qualifier to the current feature without loss of info.
If there are multiple qualifier keys with the same name we
would lose some info in the dictionary, so we append a unique
number to the end of the name in case of conflicts.
"""
# if we've got a key from before, add it to the dictionary of
# qualifiers
if self._cur_qualifier_key:
# get a unique name
unique_name = self._cur_qualifier_key
counter = 1
while self._cur_feature.qualifiers.has_key(unique_name):
unique_name = self._cur_qualifier_key + str(counter)
counter = counter + 1
self._cur_feature.qualifiers[unique_name] = \
self._cur_qualifier_value
def qualifier_key(self, content):
"""When we get a qualifier key, use it as a dictionary key.
"""
# add a qualifier if we've got one
self._add_qualifier()
# remove the / and = from the qualifier
qual_key = string.replace(content, '/', '')
qual_key = string.replace(qual_key, '=', '')
qual_key = string.strip(qual_key)
self._cur_qualifier_key = qual_key
self._cur_qualifier_value = ''
def qualifier_value(self, content):
# get rid of the quotes surrounding the qualifier if we've got 'em
qual_value = string.replace(content, '"', '')
self._cur_qualifier_value = qual_value
def origin_name(self, content):
pass
def base_count(self, content):
pass
def base_number(self, content):
pass
def sequence(self, content):
"""Add up sequence information as we get it.
"""
new_seq = string.replace(content, ' ', '')
new_seq = string.upper(new_seq)
self._seq_data += new_seq
def record_end(self, content):
"""Clean up when we've finished the record.
"""
# add the last feature in the table which hasn't been added yet
self._add_feature()
# add the sequence information
# first, determine the alphabet
# we default to an generic alphabet if we don't have a
# seq type or have strange sequence information.
seq_alphabet = Alphabet.generic_alphabet
if self._seq_type:
if string.find(self._seq_type, 'DNA') != -1:
seq_alphabet = Alphabet.generic_dna
elif string.find(self._seq_type, 'RNA') != -1:
seq_alphabet = Alphabet.generic_rna
elif self._seq_type == "PROTEIN":
seq_alphabet = Alphabet.generic_protein
# we have a bug if we get here
else:
raise ValueError("Could not determine alphabet for seq_type %s"
% self._seq_type)
# now set the sequence
self.data.seq = Seq(self._seq_data, seq_alphabet)
class _RecordConsumer(_BaseGenBankConsumer):
"""Create a GenBank Record object from scanner generated information.
"""
def __init__(self):
_BaseGenBankConsumer.__init__(self)
self.data = Record.Record()
self._cur_reference = None
self._cur_feature = None
self._cur_qualifier = None
def locus(self, content):
self.data.locus = content
def size(self, content):
self.data.size = content
def residue_type(self, content):
self.data.residue_type = content
def data_file_division(self, content):
self.data.data_file_division = content
def date(self, content):
self.data.date = content
def definition(self, content):
self.data.definition = content
def accession(self, content):
new_accessions = self._split_accessions(content)
self.data.accession.extend(new_accessions)
def nid(self, content):
self.data.nid = content
def version(self, content):
self.data.version = content
def gi(self, content):
self.data.gi = content
def keywords(self, content):
self.data.keywords = self._split_keywords(content)
def segment(self, content):
self.data.segment = content
def source(self, content):
self.data.source = content
def organism(self, content):
self.data.organism = content
def taxonomy(self, content):
self.data.taxonomy = self._split_taxonomy(content)
def reference_num(self, content):
"""Grab the reference number and signal the start of a new reference.
"""
# check if we have a reference to add
if self._cur_reference is not None:
self.data.references.append(self._cur_reference)
self._cur_reference = Record.Reference()
self._cur_reference.number = content
def reference_bases(self, content):
self._cur_reference.bases = content
def authors(self, content):
self._cur_reference.authors = content
def title(self, content):
self._cur_reference.title = content
def journal(self, content):
self._cur_reference.journal = content
def medline_id(self, content):
self._cur_reference.medline_id = content
def pubmed_id(self, content):
self._cur_reference.pubmed_id = content
def remark(self, content):
self._cur_reference.remark = content
def comment(self, content):
self.data.comment = string.join(content, "\n")
def features_line(self, content):
"""Get ready for the feature table when we reach the FEATURE line.
"""
self.start_feature_table()
def start_feature_table(self):
"""Signal the start of the feature table.
"""
# we need to add on the last reference
if self._cur_reference is not None:
self.data.references.append(self._cur_reference)
def feature_key(self, content):
"""Grab the key of the feature and signal the start of a new feature.
"""
# first add on feature information if we've got any
self._add_feature()
self._cur_feature = Record.Feature()
self._cur_feature.key = content
def _add_feature(self):
"""Utility function to add a feature to the Record.
This does all of the appropriate checking to make sure we haven't
left any info behind, and that we are only adding info if it
exists.
"""
if self._cur_feature is not None:
# if we have a left over qualifier, add it to the qualifiers
# on the current feature
if self._cur_qualifier is not None:
self._cur_feature.qualifiers.append(self._cur_qualifier)
self._cur_qualifier = None
self.data.features.append(self._cur_feature)
def location(self, content):
self._cur_feature.location = self._clean_location(content)
def qualifier_key(self, content):
# add on a qualifier if we've got one
if self._cur_qualifier is not None:
self._cur_feature.qualifiers.append(self._cur_qualifier)
self._cur_qualifier = Record.Qualifier()
self._cur_qualifier.key = content
def qualifier_value(self, content):
cur_content = self._remove_newlines(content)
# remove all spaces from the value if it is a type where spaces
# are not important
for remove_space_key in self.__class__.remove_space_keys:
if self._cur_qualifier.key.find(remove_space_key) >= 0:
cur_content = self._remove_spaces(cur_content)
self._cur_qualifier.value = self._normalize_spaces(cur_content)
def base_count(self, content):
self.data.base_counts = content
def origin_name(self, content):
self.data.origin = content
def sequence(self, content):
new_seq = string.replace(content, ' ', '')
self.data.sequence += string.upper(new_seq)
def record_end(self, content):
"""Signal the end of the record and do any necessary clean-up.
"""
# add on the last feature
self._add_feature()
def _strip_and_combine(line_list):
"""Combine multiple lines of content separated by spaces.
This function is used by the EventGenerator callback function to
combine multiple lines of information. The lines are first
stripped to remove whitepsace, and then combined so they are separated
by a space. This is a simple minded way to combine lines, but should
work for most cases.
"""
# first strip out extra whitespace
stripped_line_list = map(string.strip, line_list)
# now combine everything with spaces
return string.join(stripped_line_list, ' ')
class _Scanner:
"""Start up Martel to do the scanning of the file.
This initialzes the Martel based parser and connects it to a handler
that will generate events for a Feature Consumer.
"""
def __init__(self, debug = 0):
"""Initialize the scanner by setting up our caches.
Creating the parser takes a long time, so we want to cache it
to reduce parsing time.
Arguments:
o debug - The level of debugging that the parser should
display. Level 0 is no debugging, Level 2 displays the most
debugging info (but is much slower). See Martel documentation
for more info on this.
"""
# a listing of all tags we are interested in scanning for
# in the MartelParser
self.interest_tags = ["locus", "size", "residue_type",
"data_file_division", "date",
"definition", "accession", "nid", "version",
"gi", "keywords", "segment",
"source", "organism",
"taxonomy", "reference_num",
"reference_bases", "authors", "title",
"journal", "medline_id", "pubmed_id",
"remark", "comment",
"features_line", "feature_key",
"location", "qualifier_key",
"qualifier_value", "origin_name",
"base_count", "base_number",
"sequence", "record_end"]
# a listing of all tags which should be left alone with respect
# to whitespace handles
self.exempt_tags = ["comment"]
# make a parser that returns only the tags we are interested in
expression = Martel.select_names(genbank_format.record,
self.interest_tags)
self._parser = expression.make_parser(debug_level = debug)
def feed(self, handle, consumer):
"""Feeed a set of data into the scanner.
Arguments:
o handle - A handle with the information to parse.
o consumer - The consumer that should be informed of events.
"""
self._parser.setContentHandler(EventGenerator(consumer,
self.interest_tags,
_strip_and_combine,
self.exempt_tags))
self._parser.setErrorHandler(handler.ErrorHandler())
self._parser.parseFile(handle)
def index_file_db(genbank_file, db_name, db_directory,
identifier = "locus", aliases = ["accession"],
keywords = [], always_index = 0):
"""Index a GenBank file into a database for quick loading.
WARNING: This is very experimental and subject to change.
It requires the use of Andrew Dalke's mindy.
This is very similar to index_file, but uses a database instead
of a flat file to store the information about the genbank_file.
Arguments:
o genbank_file - The GenBank formatted file that we want to index.
o db_name - The name of the database to create. This name will allow you
to retrieve the file later.
o db_directory - The directory where the database information should be
stored.
o identifier - The primary identifier used to store records in the file
under. This will be used for retrieving them later.
o aliases - Secondary identifiers that point to the record. These can
be used for searching if a primary identifier is not found. This is
useful for GenBank since we'll index by a single identifier (the LOCUS
identifier by default) but might want to search by some other
identifier.
o keywords - More advanced Mindy features that I'm not positive
how to make full use of right now.
o always_index - A flag indicating whether or not to index a file even
if the file appears not to have changed. By default, the function will
try to skip indexing if it thinks the file hasn't changed.
"""
try:
from mindy import mindy_index, mindy_search
except ImportError:
raise SystemExit("You must have mindy installed:\n" +
"http://www.biopython.org/~dalke/mindy-0.1.tar.gz")
# try to skip the indexing if everything seems up to date
if not(always_index):
if os.path.exists(os.path.join(db_directory, db_name)):
# load up the database and see if the file size is the same
search_db = mindy_search.mindy_open(db_directory, db_name)
if search_db.mindy_data.has_key("file_sizes"):
file_size = search_db.mindy_data["file_sizes"][genbank_file]
if file_size == os.path.getsize(genbank_file):
print "File already indexed. Skipping...."
return
if not(os.path.exists(db_directory)):
os.makedirs(db_directory)
mindy_db = mindy_index.create(db_directory, db_name)
mindy_db.use_filename(genbank_file)
indexer = mindy_index.SimpleIndexer(mindy_db, "genbank_record", identifier,
aliases, keywords)
gb_format = genbank_format.record_format
# -- don't use this optimization right now, it causes genbank_format
# to be reset so that future calls to it don't return all of the item
# of interest
#if hasattr(indexer, "_wanted_elements"):
# gb_format = Martel.select_names(genbank_format.record_format,
# indexer._wanted_elements)
parser = gb_format.make_parser()
parser.setContentHandler(indexer)
mindy_db.use_filename(genbank_file)
parser.parseFile(open(genbank_file, "rb"))
class MindyDictionary:
"""Access a GenBank file using a dictionary interface, though a Mindy DB.
WARNING: This is very experimental and subject to change.
It requires the use of Andrew Dalke's mindy.
This is the Dictionary interface to use after you create an index
database using the function index_file_db.
"""
def __init__(self, db_name, db_directory, parser = None):
"""Initialize and open up a GenBank dictionary.
Arguments:
o db_name - The name of the database we should retrieve information
from.
o db_directory - The location of the database specified in db_name.
o parser - An optional argument specifying a parser object that
the records should be run through before returning the output. If
parser is None then the unprocessed contents of the file will be
returned.
"""
try:
from mindy import mindy_search
except ImportError:
raise SystemExit("You must have mindy installed:\n" +
"http://www.biopython.org/~dalke/mindy-0.1.tar.gz")
self._search = mindy_search.mindy_open(db_directory, db_name)
self._parser = parser
def __len__(self):
return len(self._search.identifiers)
def __getitem__(self, key):
"""Retrieve an item from the indexed file.
The key can be either a primary identifier or an alias. The lookup
will first try to get the file via the primary identifier, and if
it can't do this, will subsequently try to get it through the
aliases to these keys. If the aliases are ambigous, an error will
be raised.
Most of the time I find it easiest to search by aliases (the GenBank
accession numbers), but YMMV.
"""
try:
data = self._search[key]
except KeyError:
ids = self._search.aliases.get(key, [])
if len(ids) == 0:
raise KeyError("No records found for key %s" % key)
elif len(ids) != 1:
raise KeyError("Multiple records found for key %s" % key)
else:
data = self._search[ids[0]]
# run the data through the parser if one is specified
if self._parser is not None:
return self._parser.parse(File.StringHandle(data))
return File.StringHandle(data)
def __getattr__(self, name):
return getattr(self._index, name)
def keys(self):
"""Provide all identifiers for the current database.
"""
return self._search.identifiers.keys()
def aliases(self):
"""Provide all aliases in the current database.
"""
return self._search.aliases.keys()
def index_file(genbank_file, index_file, rec_to_key = None):
"""Index a GenBank file to prepare it for use as a dictionary.
Arguments:
o genbank_file - The name of the GenBank file to be index.
o index_name - The name of the index file which will be created.
o rec_to_key - A function object which, when called with a GenBank
record object, will return a key to be used for the record. If no
function is specified, then the accession numbers will be used as
the keys.
"""
if not os.path.exists(genbank_file):
raise ValueError("%s does not exist" % genbank_file)
index = Index.Index(index_file, truncate = 1)
index[Dictionary._Dictionary__filename_key] = genbank_file
gb_iter = Iterator(open(genbank_file), parser = RecordParser())
while 1:
start = gb_iter._reader.positions[gb_iter._reader.index]
rec = gb_iter.next()
length = gb_iter._reader.positions[gb_iter._reader.index] - start
# when we run out of records, stop looping
if rec is None:
break
# if we have a function to get a key, use it, otherwise use the default
if rec_to_key is not None:
key = rec_to_key(rec)
else:
# if we have accession numbers in the record use them
if len(rec.accession) >= 1:
key = rec.accession[0]
else:
key = None
# check for problems with the key
if key is None:
raise KeyError("Empty sequence key produced")
elif index.has_key(key):
raise KeyError("Duplicate key %s found" % key)
# finally apply the information to the key
index[key] = start, length
class NCBIDictionary:
"""Access GenBank using a read-only dictionary interface.
Methods:
"""
def __init__(self, database='Nucleotide', delay=5.0, parser=None):
"""NCBIDictionary([database][, delay][, parser])
Create a new Dictionary to access GenBank. database should be
either 'Nucleotide' or 'Protein'. delay is the number of
seconds to wait between each query (5 default). parser is an
optional parser object to change the results into another
form. If unspecified, then the raw contents of the file will
be returned.
"""
self.parser = parser
self.limiter = RequestLimiter(delay)
if database == 'Nucleotide':
self.format = 'GenBank'
elif database == 'Protein':
self.format = 'GenPept'
else:
raise ValueError, "database should be 'Nucleotide' or 'Protein'."
self.database = database
def __len__(self):
raise NotImplementedError, "GenBank contains lots of entries"
def clear(self):
raise NotImplementedError, "This is a read-only dictionary"
def __setitem__(self, key, item):
raise NotImplementedError, "This is a read-only dictionary"
def update(self):
raise NotImplementedError, "This is a read-only dictionary"
def copy(self):
raise NotImplementedError, "You don't need to do this..."
def keys(self):
raise NotImplementedError, "You don't really want to do this..."
def items(self):
raise NotImplementedError, "You don't really want to do this..."
def values(self):
raise NotImplementedError, "You don't really want to do this..."
def has_key(self, id):
"""S.has_key(id) -> bool"""
try:
self[id]
except KeyError:
return 0
return 1
def get(self, id, failobj=None):
try:
return self[id]
except KeyError:
return failobj
raise "How did I get here?"
def __getitem__(self, id):
"""S.__getitem__(id) -> object
Return the GenBank entry. id is the GenBank ID (gi) of the
entry. Raises a KeyError if there's an error.
"""
# First, check to see if enough time has passed since my
# last query.
self.limiter.wait()
try:
handle = NCBI.query(
'Text', self.database, dopt=self.format, uid=id)
except IOError, x:
# raise a KeyError instead of an IOError
# XXX I really should distinguish between a real IOError and
# if the id is not in the database.
raise KeyError, x
# If the id is not in the database, I get a message like:
# 'GenPept does not exist for GI "433174"\012\012'
line = handle.peekline()
if line.find('does not exist') >= 0:
raise KeyError, line
# if there is a problem with the server, we'll get back a line like:
# Please try again later. Server error for GI "7212005"
elif line.find("Please try again later.") >= 0:
raise KeyError, line
# if a sequuence has been withdrawn, we'll get a line like:
# The sequence has been intentionally withdrawn : GI "9993999"
elif line.find("intentionally withdrawn") >= 0:
raise KeyError, line
elif line.lower().find('html') >= 0:
raise KeyError, "I unexpectedly got back html-formatted data."
# Parse the record if a parser was passed in.
if self.parser is not None:
return self.parser.parse(handle)
return handle.read()
def search_for(search, database='Nucleotide', max_ids=500):
"""search_for(search[, database][, max_ids])
Search GenBank and return a list of GenBank identifiers (gi's).
search is the search string used to search the database. database
should be either 'Nucleotide' or 'Protein'. max_ids is the maximum
number of ids to retrieve (default 500).
"""
if database not in ['Nucleotide', 'Protein']:
raise ValueError, "database must be 'Nucleotide' or 'Protein'"
class ResultParser(sgmllib.SGMLParser):
# GenBank returns an HTML formatted page with lots of pretty stuff.
# I want to rip out all the genbank id's from the page. I'm going
# to do this by looking for links that retrieve records using
# the query.fcgi script.
# NM_006092
def __init__(self):
sgmllib.SGMLParser.__init__(self)
self.ids = []
def start_a(self, attrs):
# If I see a href back to the Nucleotide database from
# Entrez, then keep it.
href = None
for name, value in attrs:
if name == 'href':
href = value
if not href:
return
scheme, netloc, path, params, query, frag = urlparse.urlparse(href)
if path[-10:] != 'query.fcgi': # only want links to query.fcgi
return
# Valid queries look like (truncated):
# cmd=Retrieve&db=Nucleotide&list_uids=4503354&dopt=Gen
# I have to first split on '&' and then on '='.
params = query.split('&')
params = [x.split('=') for x in params]
list_uids = None
db = None
for name, value in params:
if name == 'list_uids':
list_uids = value
if list_uids is not None:
self.ids.append(list_uids)
parser = ResultParser()
handle = NCBI.query("Search", database, term=search, doptcmdl='DocSum',
dispmax=max_ids)
parser.feed(handle.read())
return parser.ids
def download_many(gis, callback_fn, broken_fn=None, db='Nucleotide',
delay=127.0, batchsize=500, parser=None):
"""download_many(gis, callback_fn[, delay][, batchsize])
Download many records from GenBank. gis is a list of Genbank
Gi's. Each time a record is downloaded, callback_fn is called
with the text of the record. delay is the number of seconds to
wait between requests. Waits 127 seconds by default. abatchsize
is the number of records to request each time. Default is 500
records, which is the maximum NCBI can handle.
This does not check to make sure all gi's are returned. The
client must make sure that the gi's are valid. This may be
implemented in the future.
"""
class _RecordExtractor(sgmllib.SGMLParser):
def __init__(self):
sgmllib.SGMLParser.__init__(self)
self.records = []
self._in_record = 0
def start_pre(self, attributes):
self._in_record = 1
self._current_record = []
def end_pre(self):
self.records.append(''.join(self._current_record))
self._in_record = 0
def handle_data(self, data):
if self._in_record:
self._current_record.append(data)
# parser is an undocumented parameter that allows people to
# specify an optional parser to handle each record. This is
# dangerous because the results may be malformed, and exceptions
# in the parser may disrupt the whole download process.
limiter = RequestLimiter(delay)
# Loop until all the gis are processed.
while gis:
gi_str = ','.join(gis[:batchsize])
# Make sure enough time has passed before I do another query.
limiter.wait()
# Query GenBank. This results in a HTML page that contains
# GenBank formatted records. If an ID is broken, this will
# return no data.
handle = NCBI.query('Retrieve', db, list_uids=gi_str,
dopt='GenBank', txt='on', dispmax=batchsize)
results = handle.read()
# Iterate through the results and pass the records to the
# callback. The GenBank records are between the
# tags.
extractor = _RecordExtractor()
extractor.feed(results)
for rec in extractor.records:
callback_fn(rec)
gis = gis[batchsize:]