<p>Biopython is the main package for this. It's only a few lines, so it is not reinventing the wheel.</p>
<p>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.</p>
<p>The things to watch out for are:</p>
<ul>
<li>filter by <code>type</code> (CDS?)</li>
<li>entries where there is not key among the qualifiers, so add error catching (try: ... except Exception: pass</li>
<li>the values of the qualifiers of a feature are a list. so add a [0]</li>
<li>multithreading.Pool might help if speed is a worry (which it really shouldn't). Asyncio and thread are not the way (both single core).</li>
</ul>
<p>A nice handy way to do operations on the <code>features</code> is to convert it to a pandas dataframe.</p>
<pre><code>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])
</code></pre>
<p>This means that in a jupyter notebook one can display the table to look at it (<code>feats</code> as last line in a cell or <code>IPython.display.display(feats)</code>), export it to a csv (<code>feats.to_cvs('.csv')</code>) or subset it:</p>
<pre><code>subfeats = feats.loc[(feats['type'] == 'exon') &amp; (feats['location_start']  &gt; 100)&amp; (feats['location_end']  &lt; 500)]
</code></pre>

CC BY-SA 4.0
2019-12-21T13:41:15.790Z
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!
Oren Milman
CC BY-SA 4.0
2019-12-21T16:22:42.737Z
Sorry, I am really bad at figuring out when it's a homework question and when it isn't...
Matteo Ferla
CC BY-SA 4.0
2019-12-21T16:28:38.027Z
No harm done. And clearly you did your best to help :)
Oren Milman
CC BY-SA 4.0
2021-09-13T17:14:50.760Z
I removed my statement about no homework and copypasted a chunk of code to convert it to a pandas dataframe which makes life easier.
Matteo Ferla
CC BY-SA 4.0
2019-12-21T13:20:58.357Z
2021-09-13T17:13:20.563Z
2021-09-13T17:13:20.563Z
Matteo Ferla
<p>Given a <code>.gb</code> file and a specific locus in the genome - how can I parse the relevant annotations in Python (i.e., annotations that include that locus)?</p>
<br/>
<p>I import the features using:</p>
<pre><code>SeqIO.read('my_gb_file.gb', 'gb').features
</code></pre>
<p>I then wish scan/filter sequence meta data.</p>
<br/>
<p>Is there a function or method in <code>Biopython</code> that does this?<br/></p>

<sequence-annotation><python>
Given a .gb file and a locus - how to get relevant annotations in Python?
# Given a .gb file and a locus - how to get relevant annotations in Python?

(Asked by: username_0 on 2019-12-20)

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?

**Tags:** sequence-annotation, python

---

## 1 answers

### Answer by username_1 on 2019-12-21. Score: 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)]

> **Comments:**
> **username_0**: 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!  
> **username_1**: Sorry, I am really bad at figuring out when it's a homework question and when it isn't...  
> **username_0**: No harm done. And clearly you did your best to help :)  
> **username_1**: I removed my statement about no homework and copypasted a chunk of code to convert it to a pandas dataframe which makes life easier.  

---

