Skip to content

Instantly share code, notes, and snippets.

@slowkow
Last active March 6, 2024 02:05
Show Gist options
  • Star 34 You must be signed in to star a gist
  • Fork 11 You must be signed in to fork a gist
  • Save slowkow/8101481 to your computer and use it in GitHub Desktop.
Save slowkow/8101481 to your computer and use it in GitHub Desktop.
GTF.py is a simple module for reading GTF and GFF files
#!/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
@hindlem
Copy link

hindlem commented Mar 13, 2015

Nice code two suggestions.
Line 70 causes problems when the key:value pair sequence ends in a ; leaving a blank string at the end of infos
infos = [x for x in re.split(R_SEMICOLON, fields[8]) if len(x.strip()) > 0]
Line 75 fails when there is a \s in a value...this is easily fixed by restricting the split to the first instance
key, _, value = re.split(R_KEYVALUE, info, 1)

@slowkow
Copy link
Author

slowkow commented Aug 4, 2015

Hi @hindlem, thanks for the suggestions -- they're now in my code.

(I didn't notice your comment until now because github doesn't notify the gist author about comments.)

@vmesel
Copy link

vmesel commented Aug 4, 2016

Pretty awesome code! Going to put it into my lib!

@dyl4nm4rsh4ll
Copy link

dyl4nm4rsh4ll commented Jan 25, 2018

Might I suggest pandas instead?

GFF3 = pd.read_csv(
    filepath_or_buffer=gff3_dir + 'Saccharomyces_cerevisiae.R64-1-1.91.gff3', 
    sep='\t', 
    header=None,
    names=['seqid', 'source', 'type', 'start', 'end', 'score', 'strand', 'phase', 'attributes'],
    skiprows=[i for i in xrange(25)])

GFF3 = GFF3[GFF3['source'].notnull()]

@slowkow
Copy link
Author

slowkow commented Jan 25, 2018

Dylan, my code should be able to handle lots of different types of GTF or GFF files and automatically parse the key=value items. The pandas read_csv function will not do this.

@endrebak
Copy link

pyranges parses and outputs gtf/gff.

Example:

Install instructions:

# pip install pyranges
# or 
# conda install -c bioconda pyranges

Example file:

# !head ensembl.gtf
# #!genome-build GRCh38.p10
# #!genome-version GRCh38
# #!genome-date 2013-12
# #!genome-build-accession NCBI:GCA_000001405.25
# #!genebuild-last-updated 2017-06
# 1       havana  gene    11869   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene";
# 1       havana  transcript      11869   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; tag "basic"; transcript_support_level "1";
# 1       havana  exon    11869   12227   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00002234944"; exon_version "1"; tag "basic"; transcript_support_level "1";
# 1       havana  exon    12613   12721   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00003582793"; exon_version "1"; tag "basic"; transcript_support_level "1";
# 1       havana  exon    13221   14409   .       +       .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "3"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00002312635"; exon_version "1"; tag "basic"; transcript_support_level "1";

Using pyranges:

import pyranges as pr

# as PyRanges-object
gr = pr.read_gtf("ensembl.gtf")
# +--------------+------------+--------------+-----------+-----------+------------+--------------+------------+-----------------+----------------+-------------+----------------+------------------------------------+-----------------+----------------------+-------+
# | Chromosome   | Source     | Feature      | Start     | End       | Score      | Strand       | Frame      | gene_id         | gene_version   | gene_name   | gene_source    | gene_biotype                       | transcript_id   | transcript_version   | +13   |
# | (category)   | (object)   | (category)   | (int32)   | (int32)   | (object)   | (category)   | (object)   | (object)        | (object)       | (object)    | (object)       | (object)                           | (object)        | (object)             | ...   |
# |--------------+------------+--------------+-----------+-----------+------------+--------------+------------+-----------------+----------------+-------------+----------------+------------------------------------+-----------------+----------------------+-------|
# | 1            | havana     | gene         | 11869     | 14409     | .          | +            | .          | ENSG00000223972 | 5              | DDX11L1     | havana         | transcribed_unprocessed_pseudogene | nan             | nan                  | ...   |
# | 1            | havana     | transcript   | 11869     | 14409     | .          | +            | .          | ENSG00000223972 | 5              | DDX11L1     | havana         | transcribed_unprocessed_pseudogene | ENST00000456328 | 2                    | ...   |
# | 1            | havana     | exon         | 11869     | 12227     | .          | +            | .          | ENSG00000223972 | 5              | DDX11L1     | havana         | transcribed_unprocessed_pseudogene | ENST00000456328 | 2                    | ...   |
# | 1            | havana     | exon         | 12613     | 12721     | .          | +            | .          | ENSG00000223972 | 5              | DDX11L1     | havana         | transcribed_unprocessed_pseudogene | ENST00000456328 | 2                    | ...   |
# | ...          | ...        | ...          | ...       | ...       | ...        | ...          | ...        | ...             | ...            | ...         | ...            | ...                                | ...             | ...                  | ...   |
# | 1            | ensembl    | transcript   | 120725    | 133723    | .          | -            | .          | ENSG00000238009 | 6              | AL627309.1  | ensembl_havana | lincRNA                            | ENST00000610542 | 1                    | ...   |
# | 1            | ensembl    | exon         | 133374    | 133723    | .          | -            | .          | ENSG00000238009 | 6              | AL627309.1  | ensembl_havana | lincRNA                            | ENST00000610542 | 1                    | ...   |
# | 1            | ensembl    | exon         | 129055    | 129223    | .          | -            | .          | ENSG00000238009 | 6              | AL627309.1  | ensembl_havana | lincRNA                            | ENST00000610542 | 1                    | ...   |
# | 1            | ensembl    | exon         | 120874    | 120932    | .          | -            | .          | ENSG00000238009 | 6              | AL627309.1  | ensembl_havana | lincRNA                            | ENST00000610542 | 1                    | ...   |
# +--------------+------------+--------------+-----------+-----------+------------+--------------+------------+-----------------+----------------+-------------+----------------+------------------------------------+-----------------+----------------------+-------+
# Stranded PyRanges object has 95 rows and 28 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.
# 13 hidden columns: transcript_name, transcript_source, transcript_biotype, tag, transcript_support_level, exon_number, exon_id, exon_version, (assigned, previous, ccds_id, protein_id, protein_version


# as DataFrame
df = gr.df

# Chromosome   Source     Feature   Start     End Score Strand Frame          gene_id gene_version   gene_name  ...    transcript_biotype    tag transcript_support_level exon_number          exon_id exon_version (assigned previous ccds_id protein_id protein_version
# 0           1   havana        gene   11869   14409     .      +     .  ENSG00000223972            5     DDX11L1  ...                   NaN    NaN                      NaN         NaN              NaN          NaN       NaN      NaN     NaN        NaN             NaN
# 1           1   havana  transcript   11869   14409     .      +     .  ENSG00000223972            5     DDX11L1  ...  processed_transcript  basic                        1         NaN              NaN          NaN       NaN      NaN     NaN        NaN             NaN
# 2           1   havana        exon   11869   12227     .      +     .  ENSG00000223972            5     DDX11L1  ...  processed_transcript  basic                        1           1  ENSE00002234944            1       NaN      NaN     NaN        NaN             NaN
# 3           1   havana        exon   12613   12721     .      +     .  ENSG00000223972            5     DDX11L1  ...  processed_transcript  basic                        1           2  ENSE00003582793            1       NaN      NaN     NaN        NaN             NaN
# 4           1   havana        exon   13221   14409     .      +     .  ENSG00000223972            5     DDX11L1  ...  processed_transcript  basic                        1           3  ENSE00002312635            1       NaN      NaN     NaN        NaN             NaN
# ..        ...      ...         ...     ...     ...   ...    ...   ...              ...          ...         ...  ...                   ...    ...                      ...         ...              ...          ...       ...      ...     ...        ...             ...
# 90          1   havana        exon  110953  111357     .      -     .  ENSG00000238009            6  AL627309.1  ...               lincRNA    NaN                        5           3  ENSE00001879696            1       NaN      NaN     NaN        NaN             NaN
# 91          1  ensembl  transcript  120725  133723     .      -     .  ENSG00000238009            6  AL627309.1  ...               lincRNA  basic                        5         NaN              NaN          NaN       NaN      NaN     NaN        NaN             NaN
# 92          1  ensembl        exon  133374  133723     .      -     .  ENSG00000238009            6  AL627309.1  ...               lincRNA  basic                        5           1  ENSE00003748456            1       NaN      NaN     NaN        NaN             NaN
# 93          1  ensembl        exon  129055  129223     .      -     .  ENSG00000238009            6  AL627309.1  ...               lincRNA  basic                        5           2  ENSE00003734824            1       NaN      NaN     NaN        NaN             NaN
# 94          1  ensembl        exon  120874  120932     .      -     .  ENSG00000238009            6  AL627309.1  ...               lincRNA  basic                        5           3  ENSE00003740919            1       NaN      NaN     NaN        NaN             NaN
# 
# [95 rows x 28 columns]

@sallyey
Copy link

sallyey commented Mar 17, 2020

Thank you so much. I am very impressed by your code.
I am new to Linux/Unix though. Could you let me know a Linux terminal command to execute your code interacting with the gtf.gz file?
It will be greatly appreciated! Thank you.

@slowkow
Copy link
Author

slowkow commented Mar 18, 2020

Hi @sallyey

To use this code, you'll need Python 2.

Check which version you have installed by running:

python --version

Copy the code into a file called GTF.py.

In a new Python script in the same folder as GTF.py, you can write:

import GTF
df = GTF.dataframe("filename.gtf")

@sallyey
Copy link

sallyey commented Mar 18, 2020

Hi, thank you so much for your prompt response.
I did save your code as GTF.py and I also created a new file, init_gtf.py, including the last two commands:
import GTF
df = GTF.dataframe("myfilename.gtf.gz")

And I run the codes with a linux command: python init_gtf.py
But, after the init_gtf.py is complete with running, I was not able to see any result.
Could you advise me how I can see the result from your code?

My goal is to count how many genes are annotated in my gtf.gz file, and I found your 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!

@slowkow
Copy link
Author

slowkow commented Mar 18, 2020

@sallyey You might consider asking for help on https://www.biostars.org

There, many people ask for help and provide answers every day. Hundreds of people will see your question, and some of them might offer suggestions.

Make sure you describe your problem, your goals, what you have tried, and what went wrong. The more information you provide, the more likely it is that other people will understand what you are hoping to do.

Good luck!

@sallyey
Copy link

sallyey commented Mar 18, 2020

Thanks so much!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment