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

	#!/usr/bin/env python
	"""
	GTF.py
	Kamil Slowikowski
	December 24, 2013
	Read GFF/GTF files. Works with gzip compressed files and pandas.
	http://useast.ensembl.org/info/website/upload/gff.html
	LICENSE
	This is free and unencumbered software released into the public domain.
	Anyone is free to copy, modify, publish, use, compile, sell, or
	distribute this software, either in source code form or as a compiled
	binary, for any purpose, commercial or non-commercial, and by any
	means.
	In jurisdictions that recognize copyright laws, the author or authors
	of this software dedicate any and all copyright interest in the
	software to the public domain. We make this dedication for the benefit
	of the public at large and to the detriment of our heirs and
	successors. We intend this dedication to be an overt act of
	relinquishment in perpetuity of all present and future rights to this
	software under copyright law.
	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
	EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
	MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
	IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
	OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
	ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
	OTHER DEALINGS IN THE SOFTWARE.
	For more information, please refer to <http://unlicense.org/>
	"""
	from collections import defaultdict
	import gzip
	import pandas as pd
	import re
	GTF_HEADER = ['seqname', 'source', 'feature', 'start', 'end', 'score',
	'strand', 'frame']
	R_SEMICOLON = re.compile(r'\s*;\s*')
	R_COMMA = re.compile(r'\s*,\s*')
	R_KEYVALUE = re.compile(r'(\s+|\s*=\s*)')
	def dataframe(filename):
	"""Open an optionally gzipped GTF file and return a pandas.DataFrame.
	"""
	# Each column is a list stored as a value in this dict.
	result = defaultdict(list)
	for i, line in enumerate(lines(filename)):
	for key in line.keys():
	# This key has not been seen yet, so set it to None for all
	# previous lines.
	if key not in result:
	result[key] = [None] * i
	# Ensure this row has some value for each column.
	for key in result.keys():
	result[key].append(line.get(key, None))
	return pd.DataFrame(result)
	def lines(filename):
	"""Open an optionally gzipped GTF file and generate a dict for each line.
	"""
	fn_open = gzip.open if filename.endswith('.gz') else open
	with fn_open(filename) as fh:
	for line in fh:
	if line.startswith('#'):
	continue
	else:
	yield parse(line)
	def parse(line):
	"""Parse a single GTF line and return a dict.
	"""
	result = {}
	fields = line.rstrip().split('\t')
	for i, col in enumerate(GTF_HEADER):
	result[col] = _get_value(fields[i])
	# INFO field consists of "key1=value;key2=value;...".
	infos = [x for x in re.split(R_SEMICOLON, fields[8]) if x.strip()]
	for i, info in enumerate(infos, 1):
	# It should be key="value".
	try:
	key, _, value = re.split(R_KEYVALUE, info, 1)
	# But sometimes it is just "value".
	except ValueError:
	key = 'INFO{}'.format(i)
	value = info
	# Ignore the field if there is no value.
	if value:
	result[key] = _get_value(value)
	return result
	def _get_value(value):
	if not value:
	return None
	# Strip double and single quotes.
	value = value.strip('"\'')
	# Return a list if the value has a comma.
	if ',' in value:
	value = re.split(R_COMMA, value)
	# These values are equivalent to None.
	elif value in ['', '.', 'NA']:
	return None
	return value

view raw GTF.py hosted with ❤ by GitHub

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.