KeyError when getting features from a genbank file with biopython with some accessions but not others

Question

I'm very new to python but i've been using it to extract the sequence of a gene from a genbank file. The issue is is that sometimes i'll get the output i want (prints the sequence to a file) and sometimes it will return a key error. This depends on which accession i'm using. Does anyone know why it sometimes might give a key error? I thought it might be something to do with the genbank records themselves, but they look pretty similar and the gene is there (in the gene feature qualifier). EG works with HG738867.1 but not AP019703.1. Here's my code -

from Bio import Entrez, SeqIO

gi_genome = 'accession'
name = 'acrA'
Entrez.email = 'email'
handle = Entrez.efetch(db="nucleotide", id=gi_genome, rettype="gbwithparts", retmode="text")
record = SeqIO.read(handle, "gb")
handle.close()
element = 0
for feature in record.features:
    if feature.type == 'CDS' and name in feature.qualifiers["gene"]:
        report = 'record.features[%s]' % str(element)
        gene_sequence = feature.extract(record.seq)
        with open('output.fasta', 'a') as f:
            print('>' + gi_genome + ' ' + name, file=f)
            print(gene_sequence, file=f)
        break
    else:
        element = element + 1

Here's the traceback -

Traceback (most recent call last):
  File "/home/ubuntu/Documents/Git_Branches/Project_planning/Learning/In_progress/utils/data.py", line 11, in <module>
    if feature.type == 'CDS' and name in feature.qualifiers["gene"]:
KeyError: 'gene'

Process finished with exit code 1

Thanks in advance!

Answer 1

Good question.

feature.qualifiers is a dict, if the dict doesn't have that key it will throw a KeyError. The way your code works it expects feature.qualifiers to have 'gene' as a key every time feature.type == 'CDS'. However, it is possible to have feature.type == 'CDS' without a gene entry and in the case of AP019703.1 this comes up before you reach the gene acrA.

As simeple fix is to simply use get instead:

feature.qualifiers.get("gene",[])

Now returns a empty list if there is not key "gene", which also satisfies the 'name in' part of your if statement.

If you get stuck on stuff like this use more print statements, break up your if statements and look at the whole file (example below) to understand the structure.

gi_genome = 'AP019703.1'
name = 'acrA'
Entrez.email = 'email'
handle = Entrez.efetch(db="nucleotide", id=gi_genome, rettype="gbwithparts", retmode="text")
record = SeqIO.read(handle, "gb")
handle.close()
with open('AP019703.gb' ,'w') as fout:
    SeqIO.write(record, fout, 'genbank')