Counting how many genes are annotated in gtf.gz file

Question

[sallyey2]

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!

[GenoMax]

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

[sallyey2]

Sure it looks like:

##description: evidence-based annotation of the human genome (GRCh38), version 33 (Ensembl 99)
##provider: GENCODE
##contact: [e-mail]
##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

[Ram]

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.

[sallyey2]

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

[sallyey2]

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

[wm]

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

[sallyey2]

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?

[Ram]

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)

[sallyey2]

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!

[Ram]

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.