This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Biopython SeqIO: remove gene annotation if not in a list from GenBank file

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
seqio biopython
• 3,858 views • link
updated by Ram Ram • written by Jon Jon

1 answer

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')​
• 0 views • link
written by Felix_Sim Felix_Sim

Log in to answer this question.

More posts like this