-1
$\begingroup$

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


I could retrieve the features using:

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

and then scan them to find the relevant ones, but it feels like reinventing the wheel.


Is there a function in Biopython that does that?
Or in any other well-maintained package?

$\endgroup$
0
$\begingroup$

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)]
$\endgroup$
4
  • 2
    $\begingroup$ This is not a homework question. A list in python comes with the builtin function index. It felt to me a bit like implementing list.index, so I asked. This seems to me like an extremely common operation, so I thought it is probable that such a function already exists, and I just failed to find it. Thanks anyway for your answer! $\endgroup$ Dec 21 '19 at 13:41
  • $\begingroup$ Sorry, I am really bad at figuring out when it's a homework question and when it isn't... $\endgroup$ Dec 21 '19 at 16:22
  • $\begingroup$ No harm done. And clearly you did your best to help :) $\endgroup$ Dec 21 '19 at 16:28
  • $\begingroup$ I removed my statement about no homework and copypasted a chunk of code to convert it to a pandas dataframe which makes life easier. $\endgroup$ Sep 13 at 17:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.