How to extract the protein fasta file from a genbank file?

Question

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?

Cleb · Accepted Answer · 2018-05-26T10:16:56.373

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

Is there any way to include the proteins names annotated at the gbk file as the sequence description? — F.Lira, Aug 26 '20 at 08:00

How to extract the protein fasta file from a genbank file?

1 Answers1