3
$\begingroup$

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?

Improve this question
$\endgroup$

3 Answers 3

Reset to default
3
$\begingroup$

Following up on zorbax'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

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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)? $\endgroup$
    – justinian482
    Sep 8, 2021 at 9:28
  • 1
    $\begingroup$ 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 (stackoverflow.com/questions/27298178/…), 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. $\endgroup$
    – Alex Reynolds
    Sep 8, 2021 at 10:00
2
$\begingroup$

If you have your gtf file like a DataFrame you can use:

 df[df['gene_name'].str.contains('|'.join(test_list))]
Improve this answer
$\endgroup$
1
$\begingroup$

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.

Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ 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. $\endgroup$
    – Kamil S Jaron
    Sep 7, 2021 at 11:14
  • 1
    $\begingroup$ I thought that the -w flag would prevent that? $\endgroup$
    – justinian482
    Sep 7, 2021 at 11:26
  • 1
    $\begingroup$ Ah, you are right, I overlooked -w. Should be fine, unless the genes names contain non-word constituent characters (like - or .). $\endgroup$
    – Kamil S Jaron
    Sep 7, 2021 at 13:00
  • 1
    $\begingroup$ 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. $\endgroup$
    – Ram RS
    Sep 7, 2021 at 18:58
  • $\begingroup$ Great suggestion Ram RS; I have tried it but it gives me an empty file? $\endgroup$
    – justinian482
    Sep 8, 2021 at 8:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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