Biopython SeqIO: remove gene annotation if not in a list from GenBank file

Question

[Jon]

I'm trying to "slice" a genbank record based on the annotation (record.features) and running into some problems. I have a genome genbank file, I do some filtering of gene models elsewhere which results in a list of genes that I would like to keep in the annotation, while those not in the list I'd like to remove (delete). I'm currently trying to do it by deleting features that are not in my "keep" list, however I'm getting some errors. Here's what I have thus far, and the corresponding error:

from Bio import SeqIO
input = 'test.gbk'
keep = ['gene1', 'gene2', 'gene3', 'gene4']
with open(input, 'rU') as input:
    for record in SeqIO.parse(input, 'genbank'):
        for feature in record.features:
            if "locus_tag" in feature.qualifiers:
                if not feature.qualifiers['locus_tag'][0] in keep:
                    del feature[:]
        SeqIO.write(record, 'output.gbk', 'genbank')​

However this seems to tell me I can't delete SeqRecord features, is there another way to do this?

Traceback (most recent call last):
File "<stdin>", line 6, in <module>
TypeError: 'SeqFeature' object does not support item deletion

Answer 1

[Felix_Sim]

I'm unable to test this implementation right now but maybe something along the following lines:

from Bio import SeqIO
input = 'test.gbk'
keep = ['gene1', 'gene2', 'gene3', 'gene4']
with open(input, 'rU') as input:
    for record in SeqIO.parse(input, 'genbank'):
        for index, feature in enumerate(record.features):
            if "locus_tag" in feature.qualifiers:
                if not feature.qualifiers['locus_tag'][0] in keep:
                    record.features.pop(index)
        SeqIO.write(record, 'output.gbk', 'genbank')​