Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-12-21T13:41:15.790Z
Answers.Text: 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!
Answers.UserDisplayName: Oren Milman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-12-21T16:22:42.737Z
Answers.Text: Sorry, I am really bad at figuring out when it's a homework question and when it isn't...
Answers.UserDisplayName: Matteo Ferla
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-12-21T16:28:38.027Z
Answers.Text: No harm done. And clearly you did your best to help :)
Answers.UserDisplayName: Oren Milman
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-13T17:14:50.760Z
Answers.Text: I removed my statement about no homework and copypasted a chunk of code to convert it to a pandas dataframe which makes life easier.
Answers.UserDisplayName: Matteo Ferla
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-12-21T13:20:58.357Z
Answers.LastActivityDate: 2021-09-13T17:13:20.563Z
Answers.LastEditDate: 2021-09-13T17:13:20.563Z
Answers.UserDisplayName: Matteo Ferla
Body: <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>
  
Tags: <sequence-annotation><python>
Title: Given a .gb file and a locus - how to get relevant annotations in Python?
ThreadText: # 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.  
  
  ---
  
  