GFFutils

Overview and motivation

This module is used for doing things with GFF files that are too complicated for a simple awk or grep command line call.

For example, to get a BED file of genes from a GFF file, you can use something simple like:

grep "gene" chr2.gff | awk '{print $1,$4,$5}' > out.bed

But how would you use commandline tools to get all the exons of a gene? Or exons that are found in all isoforms of a gene? Or a BED file of 3' exons from genes longer than 5 kb? How would you get the average number of isoforms for genes on the plus strand? These more complex questions are actually quite easy to answer using GFFutils.

A by-product of structuring a GFF file in this way is that it is easy to convert to something like a refFlat format -- use the GFFDB.refFlat() method for this.

See the Examples below to jump right in, or follow along for a step-by-step introduction.

Installation

Unzip the source code, and from the source directory run (with root priveliges):

python setup.py install

Now you're ready to create a GFF database and interact with it from a Python shell like IPython.

New to Python? Start here:

http://wiki.python.org/moin/BeginnersGuide

and here:

http://showmedo.com/videotutorials/python

Preparing to use GFFutils

For each GFF file you would like to use you need to create a GFF database. This database is stored in a file, and is simply a sqlite3 database. You only have to do this once for each GFF file. As long as you don't delete the database from your hard drive, you don't have to do this time-consuming step again.

The database will take roughly twice as much hard drive space as the original text file. This is the cost of interactivity and fast lookups.

You will need a GFF file to work with. If you don't already have one, here's one you can use (Drosophila melanogaster; chromosome 2L only; 42 MB):

ftp://ftp.flybase.net/genomes/Drosophila_melanogaster/dmel_r5.29_FB2010/gff/dmel-2L-r5.29.gff.gz

If you're using a FlyBase GFF file, you might want to take a look at the function GFFutils.clean_gff() in order to filter out features that may not be of interest. Assuming you have a cleaned-up GFF file, here's how to create a GFF database from either a Python script or a Python shell. To do this you simply specify the path to the GFF file and the path to the new database you'd like to create:

>>> import GFFutils
>>> GFFutils.create_gffdb('/data/annotations/dm3.gff', '/data/dm3.db')

This may take some time to run, but you only have to do this step once for every GFF file you use. Now dm3.db is the sqlite3 database that can be used in all sorts of weird and wonderful ways, outlined below.

You can find more information on exactly what's going on in the Strategy section below.

For power users, you can of course work on the sqlite3 database directly. Here's the schema:

CREATE TABLE features (
                            id text,
                            chrom text,
                            start int,
                            stop int,
                            strand text,
                            featuretype text,
                            value float,
                            source text,
                            phase text,
                            attributes text,
                            primary key (id)
                          );
CREATE TABLE relations (parent text, child text, level int, primary key(parent,child,level) );

CREATE INDEX childindex on relations (child);
CREATE INDEX featuretypes on features(featuretype);
CREATE INDEX ids ON features (id);
CREATE INDEX parentindex on relations (parent);
CREATE INDEX starts on features(start);
CREATE INDEX startstrand on features(start,strand);
CREATE INDEX stops on features(stop);
CREATE INDEX stopstrand on features(stop,strand);

Here's how to to connect to the new database from Python. First, wrap your new database in a GFFDB object:

>>> G = GFFutils.GFFDB('dm3.db')

From now on we'll be accessing the database using this new object, G, which is a GFFutils.GFFDB object.

The next couple of sections will take the form of a tutorial. If you're itching to get your hands dirty, all the methods should be documented so you can explore the object interactively. You might want to peek at the examples below, too.

Identifying what's in the database

For this section, I'm assuming that you've created a GFF database and have connected to it as described above. I'm also assuming you named the GFFDB object G.

I'm also not assuming much Python knowledge. If this sounds overly pedantic to you, feel free to jump right to the Examples!

As an introduction to using the database, let's start with answering a simple question: "What sorts of features are in the GFF file?" To do this, we'll use the features() method of the GFFDB object. The GFFDB.features() method returns a generator of the featuretypes that were in the GFF file (and which are now in the featuretype field of the sqlite3 database, which this method accesses).

Most methods in a GFFDB object return generators for performance.

Since this is the first example of using the generators returned by a GFFDB object, here are a few different ways to get the results from the generator.

Method 0: Convert iterator to a list. This is the most memory-intensive:

>>> featuretype_iterator = G.features()
>>> featuretypes = list(featuretype_iterator)

Method 1: Use iterator in a for-loop (preferred):

>>> featuretype_iterator = G.features()
>>> for featuretype in featuretype_iterator:
...     print featuretype

Method 2: Call next() incrementally on the iterator. This is the most awkward, but may sometimes be useful:

>>> featuretype_iterator = G.features()
>>> featuretype_1 = featuretype_iterator.next()
>>> featuretype_2 = featuretype_iterator.next()
>>> featuretype_3 = featuretype_iterator.next()
>>> featuretype_4 = featuretype_iterator.next()
...
...

>>> featuretypes = [featuretype1, featuretype2, ...]

It's mostly a matter of preference which method you use. However, using the for-loop approach is most memory-efficient, since only a single featuretype is in memory at one time. This is not too important for iterating through featuretypes (of which there are usually <50; typically 3-10). But when you want to iterate through 15,000 genes it can be useful.

In any case, we get something like the following. What you see on your screen depends entirely on the GFF file that you created your database from:

>>> print featuretypes

['BAC_cloned_genomic_insert',
 'CDS',
 'DNA_motif',
 'breakpoint',
 'chromosome_arm',
 'chromosome_band',
 'complex_substitution',
 'deletion',
 'enhancer',
 'exon',
 'five_prime_UTR',
 'gene',
 'insertion_site',
 'intron',
 ...
 ...
  'tRNA',
 'tandem_repeat',
 'three_prime_UTR',
 'transposable_element',
 'transposable_element_insertion_site',
 'uncharacterized_change_in_nucleotide_sequence']

Retrieving specific feature types

To retrieve just genes, just exons, or any other feature type that was in the GFF file, use the GFFDB.features_of_type() method. This will return an iterator of GFFFeature objects. GFFFeature objects are described in more detail in another section below.

'gene' was in the list of featuretypes above. Let's find out how many genes there were. In this method, we're not bringing ALL the genes into a giant list -- we'll just increment a counter. Only a single GFFFeature object is in memory at a time, which is the advantage of iterators . . .

>>> gene_count = 0
>>> for gene in G.features_of_type('gene'):
...     gene_count += 1
>>> print gene_count

This is something I found myself doing quite often, so there's a shortcut method that just does a count() in the SQL directly. Use it like this:

>>> gene_count = G.count_features_of_type('gene')

Feature types not found in the db will not return an error (maybe they should, eventually?); they just don't return anything:

>>> ncabbages = G.count_features_of_type('cabbage')
>>> print ncabbages  # zero cabbages.

Already know the ID of a feature? Get the GFFFeature object for that gene directly like this:

>>> my_favorite_feature = G['FBgn0002121']

The ID of a feature can be hard to remember. The name of a gene is often much easier to search by. However, GFF files are not consistent in how they store the name of a gene (for example, FlyBase GFF files have one name stored in the Name attribute, while other names may be stored in the Alias attribute). Nevertheless, there's a way to get named genes if the name is somewhere in the attributes field:

>>> candidates = G.attribute_search('Rm62')
>>> assert len(candidates) == 1
>>> my_favorite_gene = candidates[0]

This searches the attributes of all features of genes for the text 'Rm62'; the search is case-insensitive. Note that you get a list as a return value; that's because there may be more than one gene with that text in the attributes; it's up to you to figure out if the search returned the results you expected.

I found myself getting a gene to play around with by doing this:

>>> g = G.features_of_type('gene').next()

However, this always returns the same gene. For better testing, there's a random_feature() method that chooses a random feature out of the database. You can specify a featuretype if you'd like; otherwise you have a chance of getting any feature that was in the GFF file:

>>> g = G.random_feature('gene')

GFFFeatures in more detail

This section discusses GFFFeatures which are the things you get back when you query the database for a feature.

Just to make sure we're on the same page, here's the setup for this section, assuming you've created a GFF database called dm3.db:

>>> import GFFutils
>>> G = GFFutils.GFFDB('dm3.db')

Let's get a single GFFFeature to work with:

>>> gene = G.random_feature('gene')

GFFFeature objects, when printed, show useful information:

GFFFeature gene 'FBgn0031208': chr2L:7529-9484 (+)
#           ^          ^              ^         ^
#           |          |              |         |
# featuretype      accession   genomic coords   strand

GFFFeature objects have an attribute, id, which contains the accession in the attributes field of the original GFF file:

>>> print gene.id
'FBgn0031208'

If there was no unique ID in the original GFF file, then the ID will be the feature type plus the genomic coords (for example, "gene:chr2L:125-1340").

GFFFeature objects have many other properties:

>>> gene.start
7529

>>> gene.stop
9484

>>> gene.chrom
'chr2L'

>>> gene.featuretype
'gene'

>>> gene.strand
'+'

You can get the length of a feature with:

>>> gene_len = gene.stop - gene.start

or you can use the perhaps-more-convenient:

>>> gene_len = len(gene)

In a GFFFeature object, the GFFFeature.attributes attribute holds all the info that was in the attributes column of your GFF file. This will vary based on what was in your original GFF file. You can get a list of attribute names for a feature with:

>>> print gene.attributes._attrs

and you can access any of the attributes with a dot, then the attribute name. For example, in the GFF file I used, since the above code returned the following available attributes:

['ID', 'Name', 'Ontology_term', 'Dbxref', 'derived_computed_cyto', 'gbunit']

then we could get the ontology terms for this gene with:

>>> gene.attributes.Ontology_term
['SO:0000010', 'SO:0000087', 'GO:0008234', 'GO:0006508']

Or the DBxref (database cross-reference) for the gene with:

>>> gene.attributes.Dbxref
['FlyBase:FBan0011023',
 'FlyBase_Annotation_IDs:CG11023',
 'GB_protein:ACZ94128',
 'GB_protein:AAO41164',
 'GB:AI944728',
 'GB:AJ564667',
 'GB_protein:CAD92822',
 'GB:BF495604',
 'UniProt/TrEMBL:Q6KEV3',
 'UniProt/TrEMBL:Q86BM6',
 'INTERPRO:IPR003653',
 'EntrezGene:33155',
 'BIOGRID:59420',
 'FlyAtlas:CG11023-RA',
 'GenomeRNAi_gene:33155']

You now know enough to be able to generate a line for a BED-format file (note subtracting 1 from the start to convert to BED format's zero-based start):

>>>line = '%s\t%s\t%s\t%s\t%s\t%s\n' % (gene.chrom,
...                                     gene.start-1,
...                                     gene.stop,
...                                     gene.id,
...                                     gene.value,
...                                     gene.strand)
>>> print line

But GFFFeature objects have a convenience function, to_bed(), which also accepts a number from 3 to 6 so you can tell it how many BED fields you want returned (3 fields is the default).

So you could write a BED file of all the genes longer than 5 kb like so:

>>> fout = open('genes.bed','w')  # open a file for writing
>>> for gene in G.features_of_type('gene'):
...     if len(gene) > 5000:
...         fout.write(gene.to_bed())
>>> fout.close()

Other useful things in GFFFeature objects:

Reconstruct the GFF line for this feature, and automatically add a newline:

>>> feature.tostring()

Get the transcription start site of the feature. Note that all features have a TSS property, not just genes. It is simply the feature's start position if it's on the "+" strand or the feature's stop position if it's on the "-" strand:

>>> feature.TSS

Get the midpoint of the feature:

>>> feature.midpoint

See the Examples below for more info on this.

Navigating the hierarchy of features

Here's how to find the transcripts belonging to a gene. The GFFDB.children() and GFFDB.parents() methods take either a feature ID as an argument or a GFFFeature object. The return value is a generator of features that are children of the feature:

>>> for i in G.children(gene.id):
...     print i

Here's how to find the exons belonging to a gene. By default, level=1, which means a 'hierarchy distance' of 1 (direct parent/children, for example genes and transcripts). level=2 is analagous to grandparent/grandchild, which is used for the relationship between genes/exons. level=3 is not currently implemented:

>>> for i in G.children(gene_name, level=2):
...     print i

Note that, depending on your GFF file, you may have more than just exons as the children of genes (e.g., 3' UTRs, introns, 5' UTRs). If you just want the exons, then you can filter by feature type by specifying the featuretype keyword argument to children():

>>> for exon in G.children(gene.id, level=2, featuretype='exon'):
...     print exon

Similarly, you can get the parents with GFFDB.parents(). Here's how to get what gene an exon belongs to:

>>> exon = G.random_feature('exon')
>>> for grandparent_gene in G.parents(exon, level=2, featuretype='gene'):
...     print grandparent_gene

File format conversions

Converting features to BED files was described above; briefly:

>>> fout = open('genes.bed','w')
>>> for gene in G.features_of_type('gene'):
...     fout.write(gene.to_bed())
>>> fout.close()

Exporting a refFlat entry for one gene:

>>> print G.refFlat(gene_name)

Now create a new file, writing a refFlat entry for each gene. Note that the refFlat() method is set up such that it will return None if there were no CDSs for a particular gene. We don't want to write these to file, but do want to keep track of them.

This will take a few seconds to run:

>>> missing_cds = []
>>> fout = open('mydatabase.refFlat','w')
>>> for gene in G.features_of_type('gene'):
...     rflt = G.refFlat(gene.id)
...     if rflt is not None:
...         fout.write(rflt)
...     else:
...         missing_cds.append(gene)
>>> fout.close()

So, what were those genes that didn't have CDSs? Check the first 25:

>>> for g in missing_cds[:25]:
...     print g.attributes.Name[0]

A bunch of snoRNAs, tRNAs, etc.

GFFFeatures have a GFFFeature.tostring() method which prints back the GFF file entry as a string (with the newline included). This makes it very easy to write new GFF files containing a subset of the features in the original GFF file:

# new GFF file with genes > 5kb
>>> fout = open('big-genes.gff','w')
>>> for gene in G.features_of_type('gene'):
...     if len(gene) < 5000:
...         fout.write(gene.tostring())
>>> fout.close()

Examples

In each case, assume the following setup:

import GFFutils
GFFutils.create_gffdb('dm3.gff','dm3.db')
G = GFFutils.GFFDB('dm3.db')

Inspecting the database

print G.chromosomes()

print G.strands()

print list(G.features())

Gene count

G.count_features_of_type('gene')

Average gene length

gene_lengths = 0
gene_count = 0
for gene in G.features_of_type('gene'):
    gene_lengths += len(gene)
    gene_count += 1
mean_gene_length = float(gene_lengths) / gene_count

BED file of genes

An awk one-liner would be easier if all you needed were start/stop; adding the ID as the BED feature "name" would be a little more annoying. This example shows the automatic use the ID of the gene for the BED "name" field, making BED file creation pretty straightforward.

fout = open('genes.bed','w')
for gene in G.features_of_type('gene'):
    fout.write(gene.to_bed(6))
fout.close()

Longest gene

maxlen = 0
for gene in G.features_of_type('gene'):
    gene_len = len(gene)
    if gene_len > maxlen:
        maxlen = gene_len
        maxgene = gene
print maxlen
print maxgene

Longest gene, raw SQL

This version runs faster because it only ever looks at the start and stop columns as opposed to the above version, which returns a full GFFFeature object for each gene:

c = G.conn.cursor()
c.execute('''
    SELECT (stop-start) as LEN, *
    FROM features
    WHERE featuretype="gene"
    ORDER BY LEN DESC
''')
results = c.fetchone()
maxlen = results[0]

Average exon count

This takes several seconds to run, but as far as I know it's not something that can be done easily using grep or awk:

exon_count = 0
gene_count = 0
for gene in G.features_of_type('gene'):
    gene_exon_count = 0

    # get all grandchildren, only counting the exons
    for child in G.children(gene.id,2):
        if child.featuretype == 'exon':
            gene_exon_count += 1

    exon_count += gene_exon_count
    gene_count += 1
mean_exon_count = float(exon_count) / gene_count
print mean_exon_count

Histogram of exon lengths

(Assumes you have matplotlib installed)

from matplotlib import pyplot as p
lengths = [len(i) for i in G.features_of_type('exon')]
p.hist(lengths,bins=50)
p.xlabel('Length of exon')
p.ylabel('Frequency')
p.show()

Average number of isoforms for genes on plus strand

This assumes that transcripts are labeled as "mRNA" instead of "transcript" or something:

isoform_count = 0
gene_count = 0
for gene in G.features_of_type('gene'):
    if gene.strand == '-':
        continue
    isoforms = [i for i in G.children(gene.id) if i.featuretype=='mRNA']
    isoform_count += len(isoforms)
    gene_count += 1
mean_isoform_count = float(isoform_count) / gene_count

BED file of all exonic bases on chr2L

exons = G.features_of_type('exon', chrom='chr2L')
merged_exons = G.merge_features(exons,ignore_strand=True)
fout = open('out.bed','w')
for i in merged_exons:
    fout.write(i.to_bed())
fout.close()

Closest features

Get the closest gene (ignoring the gene you supply) and how far away it is:

g = G.random_feature('gene')
distance, closest_id = G.closest_feature(g.chr,
                                         g.start,
                                         featuretype='gene',
                                         ignore=g.id)

Get the closest upstream exon that belongs to a different gene from the one you supply:

g = G.random_feature('gene')
child_exons = G.children(g.id, level=2, featuretype='exon')
ignore = [exon.id for exon in child_exons]
distance, closest_exon = G.closest_feature(g.chr,
                                           g.start,
                                           featuretype='exon',
                                           ignore=ignore,
                                           strand=g.strand,
                                           direction='upstream')

Overlapping features

Get the exons in the first MB of chr2L that are on the plus strand:

exons_of_interest = G.overlapping_features(chrom='chr2L',
                                           start=1,
                                           stop=1e6,
                                           featuretype='exon',
                                           strand='+',
                                           completely_within=True)

Merging features

This is useful if you want to get a "meta-exon" feature that is all exons together. For example, say you have a gene with two isoforms, and you want to merge the exons together to get merged exons to indicate the presence of an exon in any isoform. Graphically:

isoform 1: [[[[[[[[[-----[[[[[[[[------------[[[
             exon1         exon2             exon3
isoform 2:     [[[[[[------------------------[[[[[[
                exon4                         exon5

merge    : [[[[[[[[[[----[[[[[[[[------------[[[[[[
            merged1         merged2           merged3

Code:

g = G.random_feature('gene')
exons = G.children(g.id, level=2, featuretype='exon')
merged_exons = G.merge_features(exons)

# If you want to create a new GFF file...
fout = open('new.gff','w')
for merged_exon in merged_exons:
    fout.write(merged_exon.tostring())
fout.close()

Imputing introns

Sometimes a GFF file doesn't explicitly include introns as features. You can construct them using the interfeatures() method. This is a pretty barebones method, so you'll have to add your own IDs and featuretypes after you have the introns created.

g = G.random_feature('gene')
exons = G.children(g.id, level=2, featuretype='exon')
introns = list(G.interfeatures(exons))

for i,intron in enumerate(introns):
    intron.featuretype='intron'
    intron.add_attribute('ID', '%s_intron:%s' % (g.id,i))

Promoter regions

Promoter regions, 1kb upstream and downstream of a gene's TSS:

g = G.random_feature('gene')
promoter = G.promoter(g.id)
g.TSS - promoter.start
promoter.stop - g.TSS

Promoter region defined as 2kb upstream:

g = G.random_feature('gene')
promoter = G.promoter(g.id, dist=2000, bidirectional=False)
g.TSS - promoter.start
promoter.stop - g.TSS

Coding genes

Useful for excluding tRNAs, rRNAs, etc . . . this returns a generator of all genes that have a CDS annotated as a child of level 2:

we_make_proteins = G.coding_genes()

Isoform counts

Useful for getting constitutive exons (i.e. exons found in all isoforms of a gene):

g = G.random_feature('gene')
n_gene_isos = G.n_gene_isoforms(g.id)
for exon in G.children(g.id,level=2,featuretype='exon'):
    if G.n_exon_isoforms(exon.id) == n_gene_isos:
        print exon.id, 'is found in all isoforms of', g.id

Arbitrary SQL commands

Note that this places a lot of trust in the user to not mess up the database!

Things at the beginning of chromosomes:

c = G.conn.cursor()
results = c.execute("""
SELECT id FROM features WHERE start BETWEEN 1 AND 100
""")
results = list(results)

Manually creating relationships:

c.execute("""
INSERT INTO relations VALUES ('fake_parent', 'fake_child', 100)
""")

Manually removing relationships:

c.execute("""
DELETE FROM relations WHERE parent='fake_parent'
""")

Strategy

The following is my reasoning for the design of this package. I'd be interested to hear any thoughts on this or ways to improve it.

A GFF database is built in several passes.

During the first pass, the lines from the GFF file are split up into fields and imported into the features table. If a "Parent" attribute is defined for the feature, then we know its first-order parent and we can enter this into the relations table.

For example, say we have the following GFF line:

chr2L FlyBase exon 8668 9276 .  + 0 ID=exon_1;Parent=mRNA_1

It will be entered into the features table like this:

ID     chrom source  type start stop  value strand phase attributes
------ ----- ------- ---- ----- ----- ----- ------ ----- -----------------------
exon_1 chr2L FlyBase CDS  8668  9276  .     +      0     ID=exon_1;mRNA_1

Since this CDS has an annotated parent, this relationship is entered into the relations table:

parent  child   level
------- ------- -----
mRNA_1  exon_1  1

Note that we can't assign any second-order parents. On this first pass, we can only add first-order parents because that's the only information that's available on a single line in the GFF file.

At some point in the GFF file though, the parent transcript is found. Here it is:

chr2L FlyBase mRNA 7529 9484 . + . ID=mRNA_1;Parent=gene_1

...and we import it into the features table, just as the exon feature was added:

ID     chrom source  type start stop  value strand phase attributes
------ ----- ------- ---- ----- ----- ----- ------ ----- -----------------------
exon_1 chr2L FlyBase CDS  8668  9276  .     +      0     ID=exon_1;mRNA_1
mRNA_1 chr2L FlyBase mRNA 7529  9484  .     +      .     ID=mRNA_1;Parent=gene_1

as well as the relations table, again just as the exon feature was added. Note however that the mRNA_1 is now in the child column. This will become important later

parent  child   level
------- ------- -----
mRNA_1  exon_1  1
gene_1  mRNA_1  1

The features table and the relations table continue to grow as the GFF file is parsed. Still, only first-order children/parents are added. When this first pass is done, indexes are created to speed up searching in the second pass.

The second pass looks at the relations table. Note that the current implementation only goes 2 levels deep; I still need to write a more general recursive form of this to support hierarchies of arbitrary depth.

In the second pass, we go through each ID in the features column, matching up IDs that are in the child column with the same ID in the parent column. In the example above, we find "exon_1" in the child column. Then we get its parent ("mRNA_1"). Then we take that parent and get it's parent by looking for "mRNA_1" in the child column and then grabbing its parent ("gene_1").

Now we know that gene_1 is the "grandparent" of exon_1, and we can enter it into the relations table as a parent of level 2:

parent  child   level
------- ------- -----
mRNA_1  exon_1  1
gene_1  mRNA_1  1
gene_1  exon_1  2

In practice, the results of the "parent search" are written to a temporary text file and then imported into the relations table as a batch in the end. This is to avoid recalculating the index each time a new row is added, something that would be extraordinarily time consuming.

Once the second pass is complete, indexes are built and the database is ready for use.

For a 130MB GFF file with 800,000+ features, the entire process takes a little under 10 mins to run. Luckily, you only need to make this time investment when you have a new GFF file; if you already have a database built then using GFFutils is quite fast.