2
$\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$
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$
  • $\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 '15 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 '15 at 12:30
  • $\begingroup$ Actually, I turned that into a question in its own right. $\endgroup$ – TheChymera Jul 14 '15 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$ – peterjc Aug 24 '15 at 15:54
1
$\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$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.