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?

$\endgroup$
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

$\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$ Sep 8 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$ Sep 8 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))]
$\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.

$\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 at 11:14
  • 1
    $\begingroup$ I thought that the -w flag would prevent that? $\endgroup$ Sep 7 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 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 at 18:58
  • $\begingroup$ Great suggestion Ram RS; I have tried it but it gives me an empty file? $\endgroup$ Sep 8 at 8:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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