Print gene locations from a .gbk file

Question

[kxd419]

Hello,

I am trying extract certain information from a gbk file I can extract the locus tag and the amino acid sequence however I am struggling to extract the gene location as it not in the same format in the file e.g.: /locus_tag="NCTC86_00002"

This is my script so far:

from Bio import GenBank
from Bio import SeqIO

gbk_filename = "HS.gb"
faa_filename = "HS_converted.faa"

input_handle = open(gbk_filename, "r")
output_handle = open(faa_filename, "w")

for seq_record in SeqIO.parse(input_handle, "genbank"):
    print "Dealing with GenBank record %s" % seq_record.id
    for seq_feature in seq_record.features :
        if seq_feature.type=="CDS" :
            assert len(seq_feature.qualifiers['translation'])==1
            output_handle.write(">%s from %s\n%s\n" % (
                   seq_feature.qualifiers['locus_tag'][0],
                   seq_record.name,
                   seq_feature.qualifiers['translation'][0]))

output_handle.close()
input_handle.close()
print "Done"

[kxd419]

So also if I wanted to print the gene annotation, not every CDS entry contains a /gene=""

Do I need to put in a if there is no /gene="" clause?

Answer 1

[mgalactus]

Hi,

the SeqFeature objects have a "location" attribute that contains the start/stop position of the feature.

from Bio import SeqIO

gbk_filename = "HS.gb"
faa_filename = "HS_converted.faa"

output_handle = open(faa_filename, "w")

for seq_record in SeqIO.parse(gbk_filename, "genbank") :
    print "Dealing with GenBank record %s" % seq_record.id
    for seq_feature in seq_record.features :
        if seq_feature.type=="CDS":
            assert len(seq_feature.qualifiers['translation'])==1
            output_handle.write(">%s from %s\n%s\n" % (
                   seq_feature.qualifiers['locus_tag'][0],
                   seq_record.name,
                   seq_feature.qualifiers['translation'][0]))
            print('Start: %d, Stop: %d, Strand: %d'%(int(seq_feature.location.start),
                                                     int(seq_feature.location.end),
                                                     seq_feature.strand))

output_handle.close()
print "Done"

Hope this helps

[kxd419]

Hi mgalactus

Thank you so much for your help. Im brand new to python and am trying to learn my best.

When I run that script I get this error:

Traceback (most recent call last):

  File "gbkaafasta.py", line 18, in <module>
    print('Start: %d, Stop: %d, Strand: %d'%(int(seq_feature.locations.start),
AttributeError: 'SeqFeature' object has no attribute 'locations'

Thank you so much for your help

[mgalactus]

My fault, it should have been "location" and not "locations" (I've updated the reply)

[kxd419]

You are fantastic! It works perfectly. Thank you so much

Answer 2

[nhaituan]

Thank you very much, here my modified codes:

from Bio import SeqIO
gbk_filename = "N16961_Ch1.gbk"
faa_filename = "N16961_Ch1.faa"
output_handle = open(faa_filename, "w")
for seq_record in SeqIO.parse(gbk_filename, "genbank") :
    print "Dealing with GenBank record %s" % seq_record.id
    for seq_feature in seq_record.features :
        if seq_feature.type=="CDS":
            assert len(seq_feature.qualifiers['translation'])==1
            output_handle.write(">%s Start: %s Stop: %s Strand: %s From %s\n" % (
                   seq_feature.qualifiers['locus_tag'][0],
                   seq_feature.location.start,
                   seq_feature.location.end,
                   seq_feature.strand,
                   seq_record.name))

output_handle.close()
print "Done"