3
$\begingroup$

Let's say I have a genbank file, e.g. the one downloaded from here, that contains entries as

 gene            190..255
                 /gene="thrL"
                 /locus_tag="b0001"
                 /gene_synonym="ECK0001; JW4367"
                 /db_xref="EcoGene:EG11277"
                 /db_xref="GeneID:944742"
 CDS             190..255
                 /gene="thrL"
                 /locus_tag="b0001"
                 /gene_synonym="ECK0001; JW4367"
                 /function="leader; Amino acid biosynthesis: Threonine"
                 /GO_process="GO:0009088 - threonine biosynthetic process"
                 /codon_start=1
                 /transl_table=11
                 /product="thr operon leader peptide"
                 /protein_id="NP_414542.1"
                 /db_xref="ASAP:ABE-0000006"
                 /db_xref="UniProtKB/Swiss-Prot:P0AD86"
                 /db_xref="EcoGene:EG11277"
                 /db_xref="GeneID:944742"
                 /translation="MKRISTTITTTITITTGNGAG"

and I want to extract the potein fasta from it, how would I do this?

I found a code sample here. When I then run:

from Bio import SeqIO

gbk_filename = 'NC_000913.3.gb'
faa_filename = 'NC_000913.3_warwick.fasta'
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()

it results in

     31     for seq_feature in seq_record.features :
     32         if seq_feature.type=="CDS" :
---> 33             assert len(seq_feature.qualifiers['translation'])==1
     34             output_handle.write(">%s from %s\n%s\n" % (
     35                    seq_feature.qualifiers['locus_tag'][0],

KeyError: 'translation'

How could this be fixed?

$\endgroup$
6
$\begingroup$

One can get it to work by using SeqIO.InsdcIO.GenBankCdsFeatureIterator:

from Bio import SeqIO

file_name = 'NC_000913.3.gb'

# stores all the CDS entries
all_entries = []

with open(file_name, 'r') as GBFile:

    GBcds = SeqIO.InsdcIO.GenBankCdsFeatureIterator(GBFile)

    for cds in GBcds:
        if cds.seq is not None:
            cds.id = cds.name
            cds.description = ''
            all_entries.append(cds)


# write file
SeqIO.write(all_entries, '{}.fasta'.format(file_name[:-3]), 'fasta')

This will then give the desired output (just the first two entries):

>b0001
MKRISTTITTTITITTGNGAG
>b0002
MRVLKFGGTSVANAERFLRVADILESNARQGQVATVLSAPAKITNHLVAMIEKTISGQDA
LPNISDAERIFAELLTGLAAAQPGFPLAQLKTFVDQEFAQIKHVLHGISLLGQCPDSINA
ALICRGEKMSIAIMAGVLEARGHNVTVIDPVEKLLAVGHYLESTVDIAESTRRIAASRIP
ADHMVLMAGFTAGNEKGELVVLGRNGSDYSAAVLAACLRADCCEIWTDVDGVYTCDPRQV
PDARLLKSMSYQEAMELSYFGAKVLHPRTITPIAQFQIPCLIKNTGNPQAPGTLIGASRD
EDELPVKGISNLNNMAMFSVSGPGMKGMVGMAARVFAAMSRARISVVLITQSSSEYSISF
CVPQSDCVRAERAMQEEFYLELKEGLLEPLAVTERLAIISVVGDGMRTLRGISAKFFAAL
ARANINIVAIAQGSSERSISVVVNNDDATTGVRVTHQMLFNTDQVIEVFVIGVGGVGGAL
LEQLKRQQSWLKNKHIDLRVCGVANSKALLTNVHGLNLENWQEELAQAKEPFNLGRLIRL
VKEYHLLNPVIVDCTSSQAVADQYADFLREGFHVVTPNKKANTSSMDYYHQLRYAAEKSR
RKFLYDTNVGAGLPVIENLQNLLNAGDELMKFSGILSGSLSYIFGKLDEGMSFSEATTLA
REMGYTEPDPRDDLSGMDVARKLLILARETGRELELADIEIEPVLPAEFNAEGDVAAFMA
NLSQLDDLFAARVAKARDEGKVLRYVGNIDEDGVCRVKIAEVDGNDPLFKVKNGENALAF
YSHYYQPLPLVLRGYGAGNDVTAAGVFADLLRTLSWKLGV
| improve this answer | |
$\endgroup$
  • $\begingroup$ Is there any way to include the proteins names annotated at the gbk file as the sequence description? $\endgroup$ – F.Lira Aug 26 at 8:00

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.