3
$\begingroup$

I can create a new Genbank record in Biopython with the following code:

from Bio import SeqIO, SeqFeature
from Bio.SeqRecord import SeqRecord
from Bio.Seq import Seq
from Bio.Alphabet import generic_dna
from Bio.Alphabet.IUPAC import IUPACAmbiguousDNA

my_sequence = "ATAGGGACUCATAGATA"
my_record = SeqRecord(Seq(my_sequence, IUPACAmbiguousDNA()))
my_feature_location = SeqFeature.FeatureLocation(2,5, strand=1)
my_feature = SeqFeature.SeqFeature(my_feature_location, type="CDS", id="Cre")
my_feature.qualifiers["foo"] = "bar"
my_record.features.append(my_feature)

print my_record.format("gb")

The corresponding genbank file then looks like:

LOCUS       .                         17 bp    DNA              UNK 01-JAN-1980
DEFINITION  .
ACCESSION   <unknown id>
VERSION     <unknown id>
KEYWORDS    .
SOURCE      .
  ORGANISM  .
            .
FEATURES             Location/Qualifiers
     CDS             3..5
                     /foo="bar"
ORIGIN
        1 atagggacuc atagata
//

How can I change this bit UNK 01-JAN-1980?

$\endgroup$
0

2 Answers 2

5
$\begingroup$

So, the problem is that you are probably using wrong class for your record. Compare Bio.SeqRecord class with Bio.GenBank.Record class. Reason why in your GenBank format you have this reference date (01-JAN-1980) is because your record has intrinsically no date attribute, so SeqRecord.format set undefined fields to default (look through the source code to see where this happens).

What you want to do for creating GenBank record is following:

$ python
Python 2.7.10 |Anaconda 2.3.0 (x86_64)| (default, May 28 2015, 17:04:42) 
[GCC 4.2.1 (Apple Inc. build 5577)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
Anaconda is brought to you by Continuum Analytics.
>>> from Bio import GenBank
>>> from Bio.GenBank import Record
>>> my_rec = Record.Record()
>>> my_rec
<Bio.GenBank.Record.Record object at 0x1030b1890>
>>> my_rec.sequence = Seq(my_sequence)
>>> my_rec.sequence
Seq('ATAGGGACUCATAGATA', Alphabet())
>>> print(my_rec)
LOCUS                         bp                                         
DEFINITION  .
KEYWORDS    .
SOURCE      .
  ORGANISM  .
            .
ORIGIN      
        1 atagggacuc atagata
//
>>> my_rec.date = '13-JUL-2015'
>>> my_rec.features.append(my_feature)
>>> my_rec.data_file_division = "ROD"
>>> print(my_rec)
LOCUS                         bp                    ROD       13-JUL-2015
DEFINITION  .
KEYWORDS    .
SOURCE      .
  ORGANISM  .
            .
FEATURES             Location/Qualifiers
type: CDS
location: [2:5](+)
id: Cre
qualifiers:
    Key: foo, Value: bar
ORIGIN      
        1 atagggacuc atagata
//

Notice, how print(my_seq) outputs data in GenBank format by default (since that is a GenBank record after all). I am sorry I cannot really talk on differences between GenBank.Record and Bio.SeqRecord and why/when use one object or another.

Anyways, looking up reference documentation is always a good habit: http://biopython.org/DIST/docs/api/Bio.GenBank.Record.Record-class.html

$\endgroup$
4
  • $\begingroup$ have you tried the code you suggest? appending my_record.date = "13-JUL-2015" or my_record.data_file_division = "ROD" does not modify my resulting record. $\endgroup$
    – TheChymera
    Jul 13, 2015 at 9:09
  • $\begingroup$ Many thanks! As a follow-up question, do you know how I could actually save a Bio.GenBank.Record object? The SeqIO.write(main_record, destination+".gb", "genbank") I used before gives me TypeError: 'Record' object is not iterable. The obvious solution, SeqIO.write([main_record], destination+".gb", "genbank") gives me TypeError: object of type 'Record' has no len()... $\endgroup$
    – TheChymera
    Jul 14, 2015 at 12:30
  • $\begingroup$ Actually, I turned that into a question in its own right. $\endgroup$
    – TheChymera
    Jul 14, 2015 at 12:50
  • $\begingroup$ For Bio.SeqIO.write(...)``` you have to give it SeqRecord`` object(s), it does not understand Bio.GenBank.Record objects. $\endgroup$
    – Peter Cock
    Aug 24, 2015 at 15:54
2
$\begingroup$

Using the SeqRecord object approach with Bio.SeqIO, you would do something like this:

my_record.annotations["data_file_division"] = "PLN"
my_record.annotations["date"] = "24-DEC-2015"

Note that only valid NCBI GenBank divisions, or their EMBL equivalent, will be accepted.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .