This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Counting how many genes are annotated in gtf.gz file

Hi, I want to count how many genes are annotated in my gtf.gz file. I found a open code from to parse my gtf.gz file and I copied the code to a file called gtf.py and also created a new Python script init_gtf.py including the following code:

#!/usr/bin/env python
import gtf
df = gtf.dataframe("my gtf filename.gz")
df.head()

And I run the codes with a linux command: python init_gtf.py

But, after the init_gtf.py is complete with running (it took few minutes), I was not able to see any result. Could you advise me how I can see the result from the code?

My goal is to count how many genes are annotated in my gtf.gz file, and I found the open code for the reference.

Since I am not familiar with both linux and python, I am struggling a lot. Please advise me. I really do appreciate your time and efforts. Thank you so much!

rna-seq rna-seq gene python
• 5,146 views • link
written by sallyey2 • sallyey2

Can you show us what your GTF file looks like? Post output of head -10 your.gtf here.

• 0 views • link
written by GenoMax GenoMax

Sure it looks like:

##description: evidence-based annotation of the human genome (GRCh38), version 33 (Ensembl 99)
##provider: GENCODE
##contact: gencode-help@ebi.ac.uk
##format: gtf
##date: 2019-12-13 

chr1    HAVANA  gene    11869   14409   .       +       .       gene_id "ENSG00000223972.5"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; 
chr1    HAVANA  transcript      11869   14409   .       +       .       gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_un 
chr1    HAVANA  exon    11869   12227   .    +       .       gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocesse 
chr1    HAVANA  exon    12613   12721   .       +       .       gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocesse 
chr1    HAVANA  exon    13221   14409   .      +       .       gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocesse 
chr1    HAVANA  transcript      12010   13670   .       +       .       gene_id "ENSG00000223972.5"; transcript_id "ENST00000450305.2"; gene_type "transcribed_unprocesse

Sorry, I am new to this website and my output looks very messy and I coudlnt include all of the output from my gtf.gz data.

it's public gtf.gz file from ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_33/gencode.v33.annotation.gtf.gz

• 0 views • link
written by sallyey2 • sallyey2

You seem to be running into line ending problems. What OS are you on?

Ideally, you should be able to pick all lines where the third column is gene and that will give you a count of genes.

• 0 views • link
written by Ram Ram

I am using Windows 10 and Ubuntu 18.04 for analysis... Is the parsing code supposed to work perfectly?

• 0 views • link
written by sallyey2 • sallyey2

Can I count gene without parsing my gtf.gz file?

• 0 views • link
written by sallyey2 • sallyey2

As above mentioned. Count the gene feature line in GTF file.

$ zcat gencode.v33.annotation.gtf.gz | awk '{a[$3]++}END{for(k in a){print k,a[k]}}'
 4
exon 1377112
CDS 763824
UTR 311868
gene 60662
start_codon 87848
Selenocysteine 119
stop_codon 80187
annotation 1
transcript 227912
• 0 views • link
written by wm wm

Thank you so much for the visualization!

I don't need to store gene_id and gene_name into a dictionary though? I still wanna practice how to store each key and values into dictionary for further analysis.

Could you help me with the parsing Python code above?

• 0 views • link
written by sallyey2 • sallyey2

Can you try these lines in your init_gtf.py script instead of what you have right now? I assume the other file is called gtf.py (like you've imported) and not GTF.py (as in the gist).

import gtf
l = gtf.lines("git_file.gtf.gz")
print(l)
• 0 views • link
written by Ram Ram

Hi thank you so much for trying this for me. When I did 

import gtf
df = gtf.dataframe("my gtf filename.gtf.gz")
print(df)

I was finally able to see the truncated version of the dataframe. Before this, I also tried df.head() instead of print(df), but it also didnt print the output. I still dont know why though..

I really do appreciate your recommendations!

• 0 views • link
written by sallyey2 • sallyey2

Glad your problem is resolved.

Please do not close/delete your post if you find the question no longer needs answered or that you'd made a typo/mistake that looks silly in hindsight. It's still a lesson and will be useful to other people that find themselves in your situation.

• 0 views • link
written by Ram Ram

0 answers

No answers yet.

Log in to answer this question.

More posts like this