Answers.Body: <p>It turned out to be a lot easier if I used <code>grep</code> to do this:</p>
  <pre><code>grep -w -f genes.txt gencode.v19.annotation.gtf &gt; sub_set.gtf
  </code></pre>
  <p><code>genes.txt</code> contains gene symbols on new lines.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-07T11:14:18.713Z
Answers.Text: the only risk you are taking with this approach is that some of the gene names might be sub-strings of other gene names and as a consequence, those would get pulled as well.
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-07T11:26:07.943Z
Answers.Text: I thought that the `-w` flag would prevent that?
Answers.UserDisplayName: justinian482
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-07T13:00:53.447Z
Answers.Text: Ah, you are right, I overlooked `-w`. Should be fine, unless the genes names contain non-word constituent characters (like - or .).
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-07T18:58:47.037Z
Answers.Text: Maybe a better way would be: `grep -wf <(awk '{print "gene_name \"$0\""}' genes.txt) gencode.v19.annotation,gtf > subset.gtf`. This will ensure that the strings are compared only to the `gene_name` tag in each feature.
Answers.UserDisplayName: Ram RS
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-08T08:55:42.353Z
Answers.Text: Great suggestion Ram RS; I have tried it but it gives me an empty file?
Answers.UserDisplayName: justinian482
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-09T15:44:23.763Z
Answers.Text: I'm not sure what the problem could be. Try just the awk command on the first 1000 lines of the big GTF, then try the grep with one of the output lines from the awk command. Step by step debugging like this is the way to solve these problems.
Answers.UserDisplayName: Ram RS
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-11T12:21:18.763Z
Answers.Text: @justinian482 have you ever edited the `genes.txt` file in Windows? Could there be Windows line endings?
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-12T14:20:50.170Z
Answers.Text: @terdon no I am using Ubuntu 20.04.3
Answers.UserDisplayName: justinian482
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-07T09:21:45.087Z
Answers.LastActivityDate: 2021-09-07T09:21:45.087Z
Answers.UserDisplayName: justinian482
Answers.Body: <p>Following up on zorbax's answer, you could read in and filter the GTF file in this way, among others:</p>
  <pre><code>#!/usr/bin/env python
  
  import gtfparse as gp
  
  gtf_file = &quot;test.gtf&quot;
  test_list = [&quot;PCNA&quot;, &quot;USP21&quot;, &quot;USP1&quot;]
  
  df = gp.read_gtf(gtf_file)
  subset = df[df['gene_name'].str.contains('|'.join(test_list))]
  print(subset)
  </code></pre>
  <p>You could further cull this subset by other attributes or columns, e.g.:</p>
  <pre><code>genes_only = subset[subset[&quot;feature&quot;] == &quot;gene&quot;]
  </code></pre>
  <p>etc.</p>
  <p>The <code>gtfparse</code> library is available via PyPi: <a href="https://pypi.org/project/gtfparse/" rel="nofollow noreferrer">https://pypi.org/project/gtfparse/</a> or Bioconda: <a href="https://bioconda.github.io/recipes/gtfparse/README.html" rel="nofollow noreferrer">https://bioconda.github.io/recipes/gtfparse/README.html</a></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-08T09:28:41.480Z
Answers.Text: How would I save the df as a proper GTF after subsetting? The separators of the final columns are different compared to the first ones (; vs \t)?
Answers.UserDisplayName: justinian482
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-08T10:00:38.027Z
Answers.Text: You could add `import pandas as pd` and then try `df.to_csv(out_filename, sep='\t')` to write out a tab-delimited file from various data frame columns. You'll probably need to first reduce the attribute columns to a single column (https://stackoverflow.com/questions/27298178/concatenate-strings-from-several-rows-using-pandas-groupby), and then remove those additional columns. Or create a new data frame from the non-attribute and attribute columns and export that. Perhaps try it out and if you get stuck, ask another question showing what you've tried.
Answers.UserDisplayName: Alex Reynolds
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-07T18:26:26.147Z
Answers.LastActivityDate: 2021-09-07T18:26:26.147Z
Answers.UserDisplayName: Alex Reynolds
Answers.Body: <p>If you have your gtf file like a DataFrame you can use:</p>
  <pre><code> df[df['gene_name'].str.contains('|'.join(test_list))]
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-07T14:21:54.620Z
Answers.LastActivityDate: 2021-09-07T14:21:54.620Z
Answers.UserDisplayName: zorbax
Body: <p>I have a GTF file (in 'test.gtf' file):</p>
  <pre><code>chr1    HAVANA  gene    10069   18412   .   +   .   gene_id &quot;ENSG00000223972.4&quot;; transcript_id &quot;ENSG00000223972.4&quot;; gene_type &quot;pseudogene&quot;; gene_status &quot;KNOWN&quot;; gene_name &quot;DDX11L1&quot;; transcript_type &quot;pseudogene&quot;; transcript_status &quot;KNOWN&quot;; transcript_name &quot;DDX11L1&quot;; level 2; havana_gene &quot;OTTHUMG00000000961.2&quot;;
  chr2    HAVANA  gene    11869   14412   .   +   .   gene_id &quot;ENSG0000023172.1&quot;; transcript_id &quot;ENSG000023172.1&quot;; gene_type &quot;pseudogene&quot;; gene_status &quot;KNOWN&quot;; gene_name &quot;PCNA&quot;; transcript_type &quot;pseudogene&quot;; transcript_status &quot;KNOWN&quot;; transcript_name &quot;PCNA&quot;; level 2; havana_gene &quot;OTTHUMG00000001931.2&quot;;
  chr3    HAVANA  gene    15769   19412   .   +   .   gene_id &quot;ENSG00000433972.1&quot;; transcript_id &quot;ENSG000023172.1&quot;; gene_type &quot;pseudogene&quot;; gene_status &quot;KNOWN&quot;; gene_name &quot;USP21&quot;; transcript_type &quot;pseudogene&quot;; transcript_status &quot;KNOWN&quot;; transcript_name &quot;USP21&quot;; level 2; havana_gene &quot;OTTHUMG00000002361.2&quot;;
  </code></pre>
  <p>I have loaded it as a pyranges object and I want to subset it so that it only contains gene from <code>test_list</code>:</p>
  <pre><code>import pyranges as pr
  import pandas as pd
  
  gtf_file = &quot;test.gtf&quot;
  test_list = [&quot;PCNA&quot;,&quot;USP21&quot;, &quot;USP1&quot;]
  
  gtf = pr.read_gtf(gtf_file)
  
  sub_gtf = gtf
  
  sub_gtf.df = gtf.df[gtf.df.gene_name.isin(test_list)]
  
  </code></pre>
  <p>But the last line throws this error:</p>
  <pre><code>File &quot;&lt;ipython-input-24-af0b4b5c7576&gt;&quot;, line 1, in &lt;module&gt;
      irongtf.df = irongtf.df[irongtf.df.gene_name.isin(test_list)]
  
    File &quot;/home/niek/.local/lib/python3.8/site-packages/pyranges/pyranges.py&quot;, line 314, in __setattr__
      _setattr(self, column_name, column)
  
    File &quot;/home/niek/.local/lib/python3.8/site-packages/pyranges/methods/attr.py&quot;, line 48, in _setattr
      df.insert(pos, column_name, column)
  
    File &quot;/home/niek/.local/lib/python3.8/site-packages/pandas/core/frame.py&quot;, line 4413, in insert
      value = self._sanitize_column(value)
  
    File &quot;/home/niek/.local/lib/python3.8/site-packages/pandas/core/frame.py&quot;, line 4504, in _sanitize_column
      com.require_length_match(value, self.index)
  
    File &quot;/home/niek/.local/lib/python3.8/site-packages/pandas/core/common.py&quot;, line 531, in require_length_match
      raise ValueError(
  
  ValueError: Length of values (2) does not match length of index (1)
  </code></pre>
  <p>What would be the right way to subset the gtf file so that it only contains the rows where <code>gene_name</code> contains the gene symbols from <code>test_list</code>?
  Or is there another tool that I could use to do this?</p>
  
Tags: <python><gtf><pyranges>
Title: Subset GTF file for specific genes
ThreadText: # Subset GTF file for specific genes
  
  (Asked by: username_0 on 2021-09-06)
  
  I have a GTF file (in 'test.gtf' file):
  
  ```
  chr1    HAVANA  gene    10069   18412   .   +   .   gene_id "ENSG00000223972.4"; transcript_id "ENSG00000223972.4"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "DDX11L1"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "DDX11L1"; level 2; havana_gene "OTTHUMG00000000961.2";
  chr2    HAVANA  gene    11869   14412   .   +   .   gene_id "ENSG0000023172.1"; transcript_id "ENSG000023172.1"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "PCNA"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "PCNA"; level 2; havana_gene "OTTHUMG00000001931.2";
  chr3    HAVANA  gene    15769   19412   .   +   .   gene_id "ENSG00000433972.1"; transcript_id "ENSG000023172.1"; gene_type "pseudogene"; gene_status "KNOWN"; gene_name "USP21"; transcript_type "pseudogene"; transcript_status "KNOWN"; transcript_name "USP21"; level 2; havana_gene "OTTHUMG00000002361.2";
  
  ```
  
  I have loaded it as a pyranges object and I want to subset it so that it only contains gene from `test_list`:
  
  ```
  import pyranges as pr
  import pandas as pd
  
  gtf_file = "test.gtf"
  test_list = ["PCNA","USP21", "USP1"]
  
  gtf = pr.read_gtf(gtf_file)
  
  sub_gtf = gtf
  
  sub_gtf.df = gtf.df[gtf.df.gene_name.isin(test_list)]
  
  ```
  
  But the last line throws this error:
  
  ```
  File "<ipython-input-24-af0b4b5c7576>", line 1, in <module>
      irongtf.df = irongtf.df[irongtf.df.gene_name.isin(test_list)]
  
    File "/home/niek/.local/lib/python3.8/site-packages/pyranges/pyranges.py", line 314, in __setattr__
      _setattr(self, column_name, column)
  
    File "/home/niek/.local/lib/python3.8/site-packages/pyranges/methods/attr.py", line 48, in _setattr
      df.insert(pos, column_name, column)
  
    File "/home/niek/.local/lib/python3.8/site-packages/pandas/core/frame.py", line 4413, in insert
      value = self._sanitize_column(value)
  
    File "/home/niek/.local/lib/python3.8/site-packages/pandas/core/frame.py", line 4504, in _sanitize_column
      com.require_length_match(value, self.index)
  
    File "/home/niek/.local/lib/python3.8/site-packages/pandas/core/common.py", line 531, in require_length_match
      raise ValueError(
  
  ValueError: Length of values (2) does not match length of index (1)
  
  ```
  
  What would be the right way to subset the gtf file so that it only contains the rows where `gene_name` contains the gene symbols from `test_list`?
  Or is there another tool that I could use to do this?
  
  **Tags:** python, gtf, pyranges
  
  
  ---
  
  ## 3 answers
  
  ### Answer by username_0 on 2021-09-07. Score: 1 | Accepted answer
  
  It turned out to be a lot easier if I used `grep` to do this:
  
  ```
  grep -w -f genes.txt gencode.v19.annotation.gtf > sub_set.gtf
  
  ```
  
  `genes.txt` contains gene symbols on new lines.
  
  > **Comments:**
  > **username_1**: the only risk you are taking with this approach is that some of the gene names might be sub-strings of other gene names and as a consequence, those would get pulled as well.  
  > **username_0**: I thought that the `-w` flag would prevent that?  
  > **username_1**: Ah, you are right, I overlooked `-w`. Should be fine, unless the genes names contain non-word constituent characters (like - or .).  
  > **username_2**: Maybe a better way would be: `grep -wf <(awk '{print "gene_name \"$0\""}' genes.txt) gencode.v19.annotation,gtf > subset.gtf`. This will ensure that the strings are compared only to the `gene_name` tag in each feature.  
  > **username_0**: Great suggestion username_2; I have tried it but it gives me an empty file?  
  > **username_2**: I'm not sure what the problem could be. Try just the awk command on the first 1000 lines of the big GTF, then try the grep with one of the output lines from the awk command. Step by step debugging like this is the way to solve these problems.  
  > **username_3**: @username_0 have you ever edited the `genes.txt` file in Windows? Could there be Windows line endings?  
  > **username_0**: @username_3 no I am using Ubuntu 20.04.3  
  
  ---
  
  ### Answer by username_4 on 2021-09-07. Score: 3
  
  Following up on username_5's answer, you could read in and filter the GTF file in this way, among others:
  
  ```
  #!/usr/bin/env python
  
  import gtfparse as gp
  
  gtf_file = "test.gtf"
  test_list = ["PCNA", "USP21", "USP1"]
  
  df = gp.read_gtf(gtf_file)
  subset = df[df['gene_name'].str.contains('|'.join(test_list))]
  print(subset)
  
  ```
  
  You could further cull this subset by other attributes or columns, e.g.:
  
  ```
  genes_only = subset[subset["feature"] == "gene"]
  
  ```
  
  etc.
  
  The `gtfparse` library is available via PyPi: <https://pypi.org/project/gtfparse/> or Bioconda: <https://bioconda.github.io/recipes/gtfparse/README.html>
  
  > **Comments:**
  > **username_0**: How would I save the df as a proper GTF after subsetting? The separators of the final columns are different compared to the first ones (; vs \t)?  
  > **username_4**: You could add `import pandas as pd` and then try `df.to_csv(out_filename, sep='\t')` to write out a tab-delimited file from various data frame columns. You'll probably need to first reduce the attribute columns to a single column (https://stackoverflow.com/questions/27298178/concatenate-strings-from-several-rows-using-pandas-groupby), and then remove those additional columns. Or create a new data frame from the non-attribute and attribute columns and export that. Perhaps try it out and if you get stuck, ask another question showing what you've tried.  
  
  ---
  
  ### Answer by username_5 on 2021-09-07. Score: 2
  
  If you have your gtf file like a DataFrame you can use:
  
  ```
   df[df['gene_name'].str.contains('|'.join(test_list))]
  
  ```
  
  
  ---
  
  