--- /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