Skip to content

Instantly share code, notes, and snippets.

Created November 30, 2015 23:23
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
What would you like to do?
import tempfile
import gffutils
import os
from argparse import ArgumentParser
def get_gtf_db(gtf, in_memory=False):
create a gffutils DB
db_file = gtf + ".db"
if os.path.exists(db_file):
return gffutils.FeatureDB(db_file)
db_file = ":memory:" if in_memory else db_file
if in_memory or not os.path.exists(db_file):
infer_extent = guess_infer_extent(gtf)
db = gffutils.create_db(gtf, dbfn=db_file,
if in_memory:
return db
return gffutils.FeatureDB(db_file)
def guess_infer_extent(gtf_file):
guess if we need to use the gene extent option when making a gffutils
database by making a tiny database of 1000 lines from the original
GTF and looking for all of the features
_, ext = os.path.splitext(gtf_file)
tmp_out = tempfile.NamedTemporaryFile(suffix=".gtf", delete=False).name
with open(tmp_out, "w") as out_handle:
count = 0
in_handle = open(gtf_file) if ext != ".gz" else
for line in in_handle:
if count > 1000:
count += 1
db = gffutils.create_db(tmp_out, dbfn=":memory:", infer_gene_extent=False)
features = [x for x in db.featuretypes()]
if "gene" in features and "transcript" in features:
return False
return True
if __name__ == "__main__":
parser = ArgumentParser()
args = parser.parse_args()
db = get_gtf_db(args.gtf_file)
for x in db.featuretypes():
print x
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment