Answers.Body: <p>I think this simple script that uses <a href="https://github.com/biocore-ntnu/pyranges" rel="nofollow noreferrer">pyranges</a> solves it. Half of the below is just setup to make a minimal reproducible example.</p>
  <pre><code># pip install pyranges
  # or
  # conda install -c bioconda pyranges
  
  import pandas as pd
  from io import StringIO 
  import pyranges as pr
  
  #########
  # setup #
  #########
  
  gff_name = &quot;gtfo.gtf&quot;
  contents = StringIO(&quot;&quot;&quot;Scaffold_1   JGI gene    22901   45904   .   +   .   ID=Genecv11000001m.g;Name=Genecv11000001m.g
  Scaffold_1  JGI mRNA    22901   45904   .   +   .   ID=PAC4GC:50510902;Name=Genecv11000001m;longest=1;Parent=Genecv11000001m.g
  Scaffold_1  JGI five_prime_UTR  22901   23284   .   +   .   ID=PAC4GC:50510902.five_prime_UTR.1;Parent=PAC4GC:50510902
  Scaffold_1  JGI CDS 23285   23423   .   +   0   ID=PAC4GC:50510902.CDS.1;Parent=PAC4GC:50510902
  Scaffold_1  JGI CDS 24031   24062   .   +   2   ID=PAC4GC:50510902.CDS.2;Parent=PAC4GC:50510902
  Scaffold_1  JGI gene    22901   45904   .   +   .   ID=Gene2.g;Name=Gene2.g
  Scaffold_1  JGI mRNA    22901   45904   .   +   .   ID=PAC4GC:WHATEVZ;Name=Gene2;longest=1;Parent=Gene2.g
  Scaffold_1  JGI five_prime_UTR  22901   23284   .   +   .   ID=PAC4GC:WHATEVZ.five_prime_UTR.1;Parent=PAC4GC:WHATEVZ
  Scaffold_1  JGI CDS 23285   23423   .   +   0   ID=PAC4GC:WHATEVZ.CDS.1;Parent=PAC4GC:WHATEVZ&quot;&quot;&quot;)
  
  pd.read_table(contents, sep=&quot;\s+&quot;, header=None).to_csv(gff_name, sep=&quot;\t&quot;, index=False, header=False)
  
  genes_to_keep = set([&quot;Genecv11000001m.g&quot;])
  
  ############
  # solution #
  ############
  
  df = pr.read_gff3(gff_name, as_df=True)
  
  def all_descendants(df, genes_to_keep):
      old_len, new_len = -1, len(genes_to_keep)
      while old_len != new_len:
          genes_to_keep.update(df[df.Parent.isin(genes_to_keep)].ID.drop_duplicates())
          old_len, new_len = new_len, len(genes_to_keep)
      return df[df.ID.isin(genes_to_keep)]
  
  desired_gff = all_descendants(df, genes_to_keep)
  print(desired_gff)
  #    Chromosome Source         Feature  Start    End Score Strand Frame                                ID               Name longest             Parent
  # 0  Scaffold_1    JGI            gene  22900  45904     .      +     .                 Genecv11000001m.g  Genecv11000001m.g     NaN                NaN
  # 1  Scaffold_1    JGI            mRNA  22900  45904     .      +     .                   PAC4GC:50510902    Genecv11000001m       1  Genecv11000001m.g
  # 2  Scaffold_1    JGI  five_prime_UTR  22900  23284     .      +     .  PAC4GC:50510902.five_prime_UTR.1                NaN     NaN    PAC4GC:50510902
  # 3  Scaffold_1    JGI             CDS  23284  23423     .      +     0             PAC4GC:50510902.CDS.1                NaN     NaN    PAC4GC:50510902
  # 4  Scaffold_1    JGI             CDS  24030  24062     .      +     2             PAC4GC:50510902.CDS.2                NaN     NaN    PAC4GC:50510902
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-03-03T21:04:13.463Z
Answers.Text: This is super nice, and I believe way simpler than the script I had come up with! Many thanks!
Answers.UserDisplayName: Felipe Almeida
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-03-03T21:05:05.500Z
Answers.Text: My scipt outsources the parsing of the gtf to a library, which is probably why it looks nicer. Much complexity is hidden behind that `read_gtf` call XD
Answers.UserDisplayName: The Unfun Cat
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-03-02T20:27:52.727Z
Answers.LastActivityDate: 2021-03-02T21:00:30.827Z
Answers.LastEditDate: 2021-03-02T21:00:30.827Z
Answers.UserDisplayName: The Unfun Cat
Answers.Body: <p>Given the nested structure of a GFF file, it may be easier to make a python script using a GFF file parser, like <a href="https://pypi.org/project/gff3/" rel="nofollow noreferrer">gff3</a>, <a href="https://daler.github.io/gffutils/" rel="nofollow noreferrer">gffutils</a>, or <a href="https://biopython.org/wiki/GFF_Parsing" rel="nofollow noreferrer">BioPython</a>.
  Dealing with nested data in <code>bash</code> or with some fancy <code>awk</code> script is probably going to trickier to write than it's worth.</p>
  <p>You can start by ensuring your list of target genes is sorted in the same order as the GFF file, then iterate through the file. When you encounter a like with <code>mRNA</code> in the 3rd column, store that ID value and extract all the subsequent lines with that same ID.</p>
  <p>If the line does not have that ID, you know you've collected all the genes, mRNAs, etc associated with that gene and can pop it from your queue.</p>
  <p>Repeat until you reach the end of the GFF file or your list of target genes is empty.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-21T00:33:46.813Z
Answers.Text: I have actually done it in python using the "flow" you suggested. Thank you!
Answers.UserDisplayName: Felipe Almeida
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-21T12:34:13.520Z
Answers.Text: Answer as a [python script](https://github.com/fmalmeida/pythonScripts/blob/master/filter_gff.py) for this task specifically
Answers.UserDisplayName: Felipe Almeida
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-21T15:04:52.800Z
Answers.Text: That script looks good! My only suggestion is that, depending on the GFF file you're loading, you may get slowed down when trying to load the entire file as a pandas `DataFrame`. If you are trying to parse a multi-gigabyte GFF file, going line-by-line keeps the memory usage low. Otherwise, what you've done looks good
Answers.UserDisplayName: James Hawley
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-21T18:37:52.733Z
Answers.Text: Precisely, I have already felt that pain :). I am actually working on a better solution for this. Thanks for the help!
Answers.UserDisplayName: Felipe Almeida
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-20T22:52:46.793Z
Answers.LastActivityDate: 2020-07-20T22:52:46.793Z
Answers.UserDisplayName: James Hawley
Body: <p>I would like to subset a GFF file (gene and nested features) from a gene list.</p>
  <p>The GFF file looks like this</p>
  <pre><code>##gff-version 3
  Scaffold_1      JGI     gene    22901   45904   .       +       .       ID=Genecv11000001m.g;Name=Genecv11000001m.g
  Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510902;Name=Genecv11000001m;longest=1;Parent=Genecv11000001m.g
  Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510902.five_prime_UTR.1;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510902.CDS.1;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510902.CDS.2;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     24192   24254   .       +       0       ID=PAC4GC:50510902.CDS.3;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510902.CDS.4;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510902.CDS.5;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510902.CDS.6;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510902.CDS.7;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510902.CDS.8;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510902.CDS.9;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510902.CDS.10;Parent=PAC4GC:50510902
  Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510902.three_prime_UTR.1;Parent=PAC4GC:50510902
  Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510903;Name=Genecv11000002m;longest=0;Parent=Genecv11000001m.g
  Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510903.five_prime_UTR.1;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510903.CDS.1;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510903.CDS.2;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     24198   24254   .       +       0       ID=PAC4GC:50510903.CDS.3;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510903.CDS.4;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510903.CDS.5;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510903.CDS.6;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510903.CDS.7;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510903.CDS.8;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510903.CDS.9;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510903.CDS.10;Parent=PAC4GC:50510903
  Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510903.three_prime_UTR.1;Parent=PAC4GC:50510903
  </code></pre>
  <p>And a have the target genes in a list such as</p>
  <pre><code>Genecv11033552m
  Genecv11003131m
  Genecv11036683m
  Genecv11012576m
  Genecv11003654m
  Genecv11012587m
  </code></pre>
  <p>I know that is possible to subset the gff using <code>grep -f gene_list.txt &lt;gff_file&gt;</code>. However this extracts only the gene and mRNA features, missing the CDS and UTR entries, while I would like to subset the gene together with all its children features (mRNA, five_prime_UTR, CDS, three_prime_UTR).</p>
  <p>This happens because the ID in CDS and UTR features are the same of the mRNA and not as in the gene feature.</p>
  <p>Any ideas?</p>
  
Tags: <python><bash>
Title: How to subset genes and its nested features from a GFF file using a gene list
ThreadText: # How to subset genes and its nested features from a GFF file using a gene list
  
  (Asked by: username_0 on 2020-07-20)
  
  I would like to subset a GFF file (gene and nested features) from a gene list.
  
  The GFF file looks like this
  
  ```
  ##gff-version 3
  Scaffold_1      JGI     gene    22901   45904   .       +       .       ID=Genecv11000001m.g;Name=Genecv11000001m.g
  Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510902;Name=Genecv11000001m;longest=1;Parent=Genecv11000001m.g
  Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510902.five_prime_UTR.1;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510902.CDS.1;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510902.CDS.2;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     24192   24254   .       +       0       ID=PAC4GC:50510902.CDS.3;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510902.CDS.4;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510902.CDS.5;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510902.CDS.6;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510902.CDS.7;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510902.CDS.8;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510902.CDS.9;Parent=PAC4GC:50510902
  Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510902.CDS.10;Parent=PAC4GC:50510902
  Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510902.three_prime_UTR.1;Parent=PAC4GC:50510902
  Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510903;Name=Genecv11000002m;longest=0;Parent=Genecv11000001m.g
  Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510903.five_prime_UTR.1;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510903.CDS.1;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510903.CDS.2;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     24198   24254   .       +       0       ID=PAC4GC:50510903.CDS.3;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510903.CDS.4;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510903.CDS.5;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510903.CDS.6;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510903.CDS.7;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510903.CDS.8;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510903.CDS.9;Parent=PAC4GC:50510903
  Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510903.CDS.10;Parent=PAC4GC:50510903
  Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510903.three_prime_UTR.1;Parent=PAC4GC:50510903
  
  ```
  
  And a have the target genes in a list such as
  
  ```
  Genecv11033552m
  Genecv11003131m
  Genecv11036683m
  Genecv11012576m
  Genecv11003654m
  Genecv11012587m
  
  ```
  
  I know that is possible to subset the gff using `grep -f gene_list.txt <gff_file>`. However this extracts only the gene and mRNA features, missing the CDS and UTR entries, while I would like to subset the gene together with all its children features (mRNA, five\_prime\_UTR, CDS, three\_prime\_UTR).
  
  This happens because the ID in CDS and UTR features are the same of the mRNA and not as in the gene feature.
  
  Any ideas?
  
  **Tags:** python, bash
  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_1 on 2021-03-02. Score: 2 | Accepted answer
  
  I think this simple script that uses [pyranges](https://github.com/biocore-ntnu/pyranges) solves it. Half of the below is just setup to make a minimal reproducible example.
  
  ```
  # pip install pyranges
  # or
  # conda install -c bioconda pyranges
  
  import pandas as pd
  from io import StringIO 
  import pyranges as pr
  
  #########
  # setup #
  #########
  
  gff_name = "gtfo.gtf"
  contents = StringIO("""Scaffold_1   JGI gene    22901   45904   .   +   .   ID=Genecv11000001m.g;Name=Genecv11000001m.g
  Scaffold_1  JGI mRNA    22901   45904   .   +   .   ID=PAC4GC:50510902;Name=Genecv11000001m;longest=1;Parent=Genecv11000001m.g
  Scaffold_1  JGI five_prime_UTR  22901   23284   .   +   .   ID=PAC4GC:50510902.five_prime_UTR.1;Parent=PAC4GC:50510902
  Scaffold_1  JGI CDS 23285   23423   .   +   0   ID=PAC4GC:50510902.CDS.1;Parent=PAC4GC:50510902
  Scaffold_1  JGI CDS 24031   24062   .   +   2   ID=PAC4GC:50510902.CDS.2;Parent=PAC4GC:50510902
  Scaffold_1  JGI gene    22901   45904   .   +   .   ID=Gene2.g;Name=Gene2.g
  Scaffold_1  JGI mRNA    22901   45904   .   +   .   ID=PAC4GC:WHATEVZ;Name=Gene2;longest=1;Parent=Gene2.g
  Scaffold_1  JGI five_prime_UTR  22901   23284   .   +   .   ID=PAC4GC:WHATEVZ.five_prime_UTR.1;Parent=PAC4GC:WHATEVZ
  Scaffold_1  JGI CDS 23285   23423   .   +   0   ID=PAC4GC:WHATEVZ.CDS.1;Parent=PAC4GC:WHATEVZ""")
  
  pd.read_table(contents, sep="\s+", header=None).to_csv(gff_name, sep="\t", index=False, header=False)
  
  genes_to_keep = set(["Genecv11000001m.g"])
  
  ############
  # solution #
  ############
  
  df = pr.read_gff3(gff_name, as_df=True)
  
  def all_descendants(df, genes_to_keep):
      old_len, new_len = -1, len(genes_to_keep)
      while old_len != new_len:
          genes_to_keep.update(df[df.Parent.isin(genes_to_keep)].ID.drop_duplicates())
          old_len, new_len = new_len, len(genes_to_keep)
      return df[df.ID.isin(genes_to_keep)]
  
  desired_gff = all_descendants(df, genes_to_keep)
  print(desired_gff)
  #    Chromosome Source         Feature  Start    End Score Strand Frame                                ID               Name longest             Parent
  # 0  Scaffold_1    JGI            gene  22900  45904     .      +     .                 Genecv11000001m.g  Genecv11000001m.g     NaN                NaN
  # 1  Scaffold_1    JGI            mRNA  22900  45904     .      +     .                   PAC4GC:50510902    Genecv11000001m       1  Genecv11000001m.g
  # 2  Scaffold_1    JGI  five_prime_UTR  22900  23284     .      +     .  PAC4GC:50510902.five_prime_UTR.1                NaN     NaN    PAC4GC:50510902
  # 3  Scaffold_1    JGI             CDS  23284  23423     .      +     0             PAC4GC:50510902.CDS.1                NaN     NaN    PAC4GC:50510902
  # 4  Scaffold_1    JGI             CDS  24030  24062     .      +     2             PAC4GC:50510902.CDS.2                NaN     NaN    PAC4GC:50510902
  
  ```
  
  > **Comments:**
  > **username_0**: This is super nice, and I believe way simpler than the script I had come up with! Many thanks!  
  > **username_1**: My scipt outsources the parsing of the gtf to a library, which is probably why it looks nicer. Much complexity is hidden behind that `read_gtf` call XD  
  
  ---
  
  ### Answer by username_2 on 2020-07-20. Score: 1
  
  Given the nested structure of a GFF file, it may be easier to make a python script using a GFF file parser, like [gff3](https://pypi.org/project/gff3/), [gffutils](https://daler.github.io/gffutils/), or [BioPython](https://biopython.org/wiki/GFF_Parsing).
  Dealing with nested data in `bash` or with some fancy `awk` script is probably going to trickier to write than it's worth.
  
  You can start by ensuring your list of target genes is sorted in the same order as the GFF file, then iterate through the file. When you encounter a like with `mRNA` in the 3rd column, store that ID value and extract all the subsequent lines with that same ID.
  
  If the line does not have that ID, you know you've collected all the genes, mRNAs, etc associated with that gene and can pop it from your queue.
  
  Repeat until you reach the end of the GFF file or your list of target genes is empty.
  
  > **Comments:**
  > **username_0**: I have actually done it in python using the "flow" you suggested. Thank you!  
  > **username_0**: Answer as a [python script](https://github.com/fmalmeida/pythonScripts/blob/master/filter_gff.py) for this task specifically  
  > **username_2**: That script looks good! My only suggestion is that, depending on the GFF file you're loading, you may get slowed down when trying to load the entire file as a pandas `DataFrame`. If you are trying to parse a multi-gigabyte GFF file, going line-by-line keeps the memory usage low. Otherwise, what you've done looks good  
  > **username_0**: Precisely, I have already felt that pain :). I am actually working on a better solution for this. Thanks for the help!  
  
  ---
  
  