+++ /dev/null
-# Copyright 2007-2009 by Peter Cock. 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..
-
-"""Bio.SeqIO support for the "genbank" and "embl" file formats.
-
-You are expected to use this module via the Bio.SeqIO functions.
-Note that internally this module calls Bio.GenBank to do the actual
-parsing of both GenBank and EMBL files.
-
-See also:
-
-International Nucleotide Sequence Database Collaboration
-http://www.insdc.org/
-
-GenBank
-http://www.ncbi.nlm.nih.gov/Genbank/
-
-EMBL Nucleotide Sequence Database
-http://www.ebi.ac.uk/embl/
-
-DDBJ (DNA Data Bank of Japan)
-http://www.ddbj.nig.ac.jp/
-"""
-
-from Bio.Seq import UnknownSeq
-from Bio.GenBank.Scanner import GenBankScanner, EmblScanner
-from Bio import Alphabet
-from Interfaces import SequentialSequenceWriter
-
-# NOTE
-# ====
-# The "brains" for parsing GenBank and EMBL files (and any
-# other flat file variants from the INSDC in future) is in
-# Bio.GenBank.Scanner (plus the _FeatureConsumer in Bio.GenBank)
-
-def GenBankIterator(handle) :
- """Breaks up a Genbank file into SeqRecord objects.
-
- Every section from the LOCUS line to the terminating // becomes
- a single SeqRecord with associated annotation and features.
-
- Note that for genomes or chromosomes, there is typically only
- one record."""
- #This calls a generator function:
- return GenBankScanner(debug=0).parse_records(handle)
-
-def EmblIterator(handle) :
- """Breaks up an EMBL file into SeqRecord objects.
-
- Every section from the LOCUS line to the terminating // becomes
- a single SeqRecord with associated annotation and features.
-
- Note that for genomes or chromosomes, there is typically only
- one record."""
- #This calls a generator function:
- return EmblScanner(debug=0).parse_records(handle)
-
-def GenBankCdsFeatureIterator(handle, alphabet=Alphabet.generic_protein) :
- """Breaks up a Genbank file into SeqRecord objects for each CDS feature.
-
- Every section from the LOCUS line to the terminating // can contain
- many CDS features. These are returned as with the stated amino acid
- translation sequence (if given).
- """
- #This calls a generator function:
- return GenBankScanner(debug=0).parse_cds_features(handle, alphabet)
-
-def EmblCdsFeatureIterator(handle, alphabet=Alphabet.generic_protein) :
- """Breaks up a EMBL file into SeqRecord objects for each CDS feature.
-
- Every section from the LOCUS line to the terminating // can contain
- many CDS features. These are returned as with the stated amino acid
- translation sequence (if given).
- """
- #This calls a generator function:
- return EmblScanner(debug=0).parse_cds_features(handle, alphabet)
-
-class GenBankWriter(SequentialSequenceWriter) :
- HEADER_WIDTH = 12
- MAX_WIDTH = 80
-
- def _write_single_line(self, tag, text) :
- "Used in the the 'header' of each GenBank record."""
- assert len(tag) < self.HEADER_WIDTH
- assert len(text) < self.MAX_WIDTH - self.HEADER_WIDTH, \
- "Annotation %s too long for %s line" % (repr(text), tag)
- self.handle.write("%s%s\n" % (tag.ljust(self.HEADER_WIDTH),
- text.replace("\n"," ")))
-
- def _write_multi_line(self, tag, text) :
- "Used in the the 'header' of each GenBank record."""
- #TODO - Do the line spliting while preserving white space?
- max_len = self.MAX_WIDTH - self.HEADER_WIDTH
- assert len(tag) < self.HEADER_WIDTH
- text = text.strip()
- if len(text) < max_len :
- self._write_single_line(tag, text)
- return
-
- words = text.split()
- assert max([len(w) for w in words]) < max_len, \
- "Your description cannot be broken into nice lines!"
- text = ""
- while words and len(text) + 1 + len(words[0]) < max_len :
- text += " " + words.pop(0)
- text = text.strip()
- assert len(text) < max_len
- self._write_single_line(tag, text)
- while words :
- text = ""
- while words and len(text) + 1 + len(words[0]) < max_len :
- text += " " + words.pop(0)
- text = text.strip()
- assert len(text) < max_len
- self._write_single_line("", text)
- assert not words
-
- def _write_the_first_line(self, record) :
- """Write the LOCUS line."""
-
- locus = record.name
- if not locus or locus == "<unknown name>" :
- locus = record.id
- if not locus or locus == "<unknown id>" :
- locus = self._get_annotation_str(record, "accession", just_first=True)
- if len(locus) > 16 :
- raise ValueError("Locus identifier %s is too long" % repr(locus))
-
- if len(record) > 99999999999 :
- #Currently GenBank only officially support up to 350000, but
- #the length field can take eleven digits
- raise ValueError("Sequence too long!")
-
- #Get the base alphabet (underneath any Gapped or StopCodon encoding)
- a = Alphabet._get_base_alphabet(record.seq.alphabet)
- if not isinstance(a, Alphabet.Alphabet) :
- raise TypeError("Invalid alphabet")
- elif isinstance(a, Alphabet.ProteinAlphabet) :
- units = "bp"
- elif isinstance(a, Alphabet.NucleotideAlphabet) :
- units = "aa"
- else :
- #Must be something like NucleotideAlphabet or
- #just the generic Alphabet (default for fasta files)
- raise ValueError("Need a Nucleotide or Protein alphabet")
-
- #Get the molecule type
- #TODO - record this explicitly in the parser?
- if isinstance(a, Alphabet.ProteinAlphabet) :
- mol_type = ""
- elif isinstance(a, Alphabet.DNAAlphabet) :
- mol_type = "DNA"
- elif isinstance(a, Alphabet.RNAAlphabet) :
- mol_type = "RNA"
- else :
- #Must be something like NucleotideAlphabet or
- #just the generic Alphabet (default for fasta files)
- raise ValueError("Need a DNA, RNA or Protein alphabet")
-
- try :
- division = record.annotations["data_file_division"]
- except KeyError :
- division = "UNK"
- if division not in ["PRI","ROD","MAM","VRT","INV","PLN","BCT",
- "VRL","PHG","SYN","UNA","EST","PAT","STS",
- "GSS","HTG","HTC","ENV"] :
- division = "UNK"
-
- assert len(units) == 2
- assert len(division) == 3
- #TODO - date
- #TODO - mol_type
- line = "LOCUS %s %s %s %s %s 01-JAN-1980\n" \
- % (locus.ljust(16),
- str(len(record)).rjust(11),
- units,
- mol_type.ljust(6),
- division)
- assert len(line) == 79+1, repr(line) #plus one for new line
-
- assert line[12:28].rstrip() == locus, \
- 'LOCUS line does not contain the locus at the expected position:\n' + line
- assert line[28:29] == " "
- assert line[29:40].lstrip() == str(len(record)), \
- 'LOCUS line does not contain the length at the expected position:\n' + line
-
- #Tests copied from Bio.GenBank.Scanner
- assert line[40:44] in [' bp ', ' aa '] , \
- 'LOCUS line does not contain size units at expected position:\n' + line
- assert line[44:47] in [' ', 'ss-', 'ds-', 'ms-'], \
- 'LOCUS line does not have valid strand type (Single stranded, ...):\n' + line
- assert line[47:54].strip() == "" \
- or line[47:54].strip().find('DNA') != -1 \
- or line[47:54].strip().find('RNA') != -1, \
- 'LOCUS line does not contain valid sequence type (DNA, RNA, ...):\n' + line
- assert line[54:55] == ' ', \
- 'LOCUS line does not contain space at position 55:\n' + line
- assert line[55:63].strip() in ['','linear','circular'], \
- 'LOCUS line does not contain valid entry (linear, circular, ...):\n' + line
- assert line[63:64] == ' ', \
- 'LOCUS line does not contain space at position 64:\n' + line
- assert line[67:68] == ' ', \
- 'LOCUS line does not contain space at position 68:\n' + line
- assert line[70:71] == '-', \
- 'LOCUS line does not contain - at position 71 in date:\n' + line
- assert line[74:75] == '-', \
- 'LOCUS line does not contain - at position 75 in date:\n' + line
-
- self.handle.write(line)
-
- def _get_annotation_str(self, record, key, default=".", just_first=False) :
- """Get an annotation dictionary entry (as a string).
-
- Some entries are lists, in which case if just_first=True the first entry
- is returned. If just_first=False (default) this verifies there is only
- one entry before returning it."""
- try :
- answer = record.annotations[key]
- except KeyError :
- return default
- if isinstance(answer, list) :
- if not just_first : assert len(answer) == 1
- return str(answer[0])
- else :
- return str(answer)
-
- def _write_sequence(self, record):
- #Loosely based on code from Howard Salis
- #TODO - Force lower case?
- LETTERS_PER_LINE = 60
- SEQUENCE_INDENT = 9
-
- if isinstance(record.seq, UnknownSeq) :
- #We have already recorded the length, and there is no need
- #to record a long sequence of NNNNNNN...NNN or whatever.
- return
-
- data = self._get_seq_string(record) #Catches sequence being None
- seq_len = len(data)
- for line_number in range(0,seq_len,LETTERS_PER_LINE):
- self.handle.write(str(line_number+1).rjust(SEQUENCE_INDENT))
- for words in range(line_number,min(line_number+LETTERS_PER_LINE,seq_len),10):
- self.handle.write(" %s" % data[words:words+10])
- self.handle.write("\n")
-
- def write_record(self, record):
- """Write a single record to the output file."""
- handle = self.handle
- self._write_the_first_line(record)
-
- accession = self._get_annotation_str(record, "accession",
- record.id.split(".",1)[0],
- just_first=True)
- acc_with_version = accession
- if record.id.startswith(accession+".") :
- try :
- acc_with_version = "%s.%i" \
- % (accession, int(record.id.split(".",1)[1]))
- except ValueError :
- pass
- gi = self._get_annotation_str(record, "gi", just_first=True)
-
- descr = record.description
- if descr == "<unknown description>" : descr = "."
- self._write_multi_line("DEFINITION", descr)
-
- self._write_single_line("ACCESSION", accession)
- if gi != "." :
- self._write_single_line("VERSION", "%s GI:%s" % (acc_with_version,gi))
- else :
- self._write_single_line("VERSION", "%s" % (acc_with_version))
-
- try :
- #List of strings
- keywords = "; ".join(record.annotations["keywords"])
- except KeyError :
- keywords = "."
- self._write_multi_line("KEYWORDS", keywords)
-
- self._write_multi_line("SOURCE", \
- self._get_annotation_str(record, "source"))
- #The ORGANISM line MUST be a single line, as any continuation is the taxonomy
- org = self._get_annotation_str(record, "organism")
- if len(org) > self.MAX_WIDTH - self.HEADER_WIDTH :
- org = org[:self.MAX_WIDTH - self.HEADER_WIDTH-4]+"..."
- self._write_single_line(" ORGANISM", org)
- try :
- #List of strings
- taxonomy = "; ".join(record.annotations["taxonomy"])
- except KeyError :
- taxonomy = "."
- self._write_multi_line("", taxonomy)
-
- #TODO - References...
- handle.write("FEATURES Location/Qualifiers\n")
- for feature in record.features :
- self._write_feature(feature)
- handle.write("ORIGIN\n")
- self._write_sequence(record)
- handle.write("//\n")
-
- def _write_feature(self, feature):
- """Write a single SeqFeature object to features table.
-
- Not implemented yet, but this stub exists in the short term to
- facilitate working on writing GenBank files with a sub-class."""
- #TODO - Features...
- pass
-
-if __name__ == "__main__" :
- print "Quick self test"
- import os
- from StringIO import StringIO
-
- def check_genbank_writer(records) :
- handle = StringIO()
- GenBankWriter(handle).write_file(records)
- handle.seek(0)
-
- records2 = list(GenBankIterator(handle))
-
- assert len(records) == len(records2)
- for r1, r2 in zip(records, records2) :
- #The SwissProt parser may leave \n in the description...
- assert r1.description.replace("\n", " ") == r2.description
- assert r1.id == r2.id
- assert r1.name == r2.name
- assert str(r1.seq) == str(r2.seq)
- for key in ["gi", "keywords", "source", "taxonomy"] :
- if key in r1.annotations :
- assert r1.annotations[key] == r2.annotations[key], key
- for key in ["organism"] :
- if key in r1.annotations :
- v1 = r1.annotations[key]
- v2 = r2.annotations[key]
- assert isinstance(v1, str) and isinstance(v2, str)
- #SwissProt organism can be too long to record in GenBank format
- assert v1 == v2 or \
- (v2.endswith("...") and v1.startswith(v2[:-3])), key
-
- for filename in os.listdir("../../Tests/GenBank") :
- if not filename.endswith(".gbk") and not filename.endswith(".gb") :
- continue
- print filename
-
- handle = open("../../Tests/GenBank/%s" % filename)
- records = list(GenBankIterator(handle))
- handle.close()
-
- check_genbank_writer(records)
-
- for filename in os.listdir("../../Tests/EMBL") :
- if not filename.endswith(".embl") :
- continue
- print filename
-
- handle = open("../../Tests/EMBL/%s" % filename)
- records = list(EmblIterator(handle))
- handle.close()
-
- check_genbank_writer(records)
-
- from Bio import SeqIO
- for filename in os.listdir("../../Tests/SwissProt") :
- if not filename.startswith("sp") :
- continue
- print filename
-
- handle = open("../../Tests/SwissProt/%s" % filename)
- records = list(SeqIO.parse(handle,"swiss"))
- handle.close()
-
- check_genbank_writer(records)
-
git://source.jalview.org / jabaws.git / blobdiff