Given a .gb file and a locus - how to get relevant annotations in Python?

Question

Given a .gb file and a specific locus in the genome - how can I parse the relevant annotations in Python (i.e., annotations that include that locus)?

I import the features using:

SeqIO.read('my_gb_file.gb', 'gb').features

I then wish scan/filter sequence meta data.

Is there a function or method in Biopython that does this?

Answer 1

Biopython is the main package for this. It's only a few lines, so it is not reinventing the wheel.

So you want to iterate across the features table (a list in Biopython) of a record and find the case where the qualifier['locus'][0] matches your query.

The things to watch out for are:

• filter by type (CDS?)
• entries where there is not key among the qualifiers, so add error catching (try: ... except Exception: pass
• the values of the qualifiers of a feature are a list. so add a [0]
• multithreading.Pool might help if speed is a worry (which it really shouldn't). Asyncio and thread are not the way (both single core).

A nice handy way to do operations on the features is to convert it to a pandas dataframe.

from Bio.SeqFeature import SeqFeature
from Bio import SeqIO
from Bio.Seq import Seq
import pandas as pd

        def convert(feat: SeqFeature, flatten:bool=True):
    if flatten:
        clean = lambda v: v[0] if isinstance(v, list) else v
        qualifiers = {k: clean(v) for k, v in feat.qualifiers.items()}
    else:
        qualifiers = feat.qualifiers
    return {'id': feat.id,
            'type': feat.type,
            'location_strand': feat.location.strand, 
            'location_start': feat.location.start,
            'location_end': feat.location.end,
            **qualifiers}

transcript = SeqIO.read('👾👾👾.gb', 'genbank')
feats = pd.DataFrame([convert(feat) for feat in transcript.features])

This means that in a jupyter notebook one can display the table to look at it (feats as last line in a cell or IPython.display.display(feats)), export it to a csv (feats.to_cvs('.csv')) or subset it:

subfeats = feats.loc[(feats['type'] == 'exon') & (feats['location_start']  > 100)& (feats['location_end']  < 500)]