daler/gffutils — pull request #16 by yarden (2013-05-11 17:12:49+00:00)
Title: cli-dev merge with refactor branch

These are the main changes I made:

## Library features
gffwriter.py: helper for serializing GFF database to a GFF text file.
Uses canonical output ordering to make uniform looking GFF files with
predictable order.

- Supports writing to a stream or filename

- Allows one to write to a GFF file "in place"; i.e. to overwrite
the current GFF file with an altered version of itself. Uses
named temporary files under the hood like pybedtools. This is useful
if you want to standardize/sanitize/clean a bunch of GFF files without
creating more GFF files.

interface.py:
- iter_by_parent_childs: method for iteration through db in
'parent-centric' way. See code for sanitize_gff_file to see
this function in action

db.py:
- added `__len__()` feature

helpers.py:
- sanitize_gff: sanitize a GFF file
- get_db_fname: returns a on-the-fly database for a given GFF filename if one does not
already exist for it. Probably not going to be very commonly used in light of ':memory:'
databases, but handy for serialized/named dbs.
- is_gff_db: helper function to determine if a path is a GFF db or not (goes by
'.db' extension)

Example of the kind of features added. sanitize_gff_db sanitizes GFF files
by ensuring that start < stop (like "clean_gff") but additionally adds a
gene id-like attribute that is propagated through all children nodes of each
gene. That makes the resulting GFF file grep-able without interfering with
the standard GFF format:

```
def sanitize_gff_db(db, gid_field="gid"):
def sanitized_iterator():
# Iterate through the database by each gene's records
for gene_recs in db.iter_by_parent_childs():
# The gene's ID
gene_id = gene_recs[0].id
for rec in gene_recs:
# Fixup coordinates if necessary
if rec.start > rec.stop:
rec.start, rec.stop = rec.stop, rec.start
# Add a gene id field to each gene's records
rec.attributes[gid_field] = [gene_id]
yield rec
# Return sanitized GFF database
sanitized_db = \
gffutils.create_db(sanitized_iterator(), ":memory:",
verbose=False)
return sanitized_db
```

Note that 'iter_by_parent_childs' returns an iterator of parent-child record lists
(where the first entry is the parent, the remaining entries are children.)
By default it uses "gene" entries as parents but takes a 'featuretype' argument
to use other features.

Coupled with GFFWriter, one can do this:

```
# Write to 'mygff.gff3' in place: the sanitized
# gff file will overwrite 'mygff.gff3'
gff_out = GFFWriter("mygff.gff3", in_place=True)
# Sanitize the database and write it out
sanitized_db = helpers.sanitize_gff_db()
# Write the GFF entries in canonical order
for gene_rec in sanitized_db.all_features(featuretype="gene"):
gff_out.write_gene_recs(sanitized_db, gene_rec.id)
gff_out.close()
```

This sanitizes the file and writes it out using canonical GFFWriter order.
Now 'mygff.gff3' is a sanitized version of itself. See below for detailed
example of 'sanitize'.

## Command line interface features

1. You can now sanitize a file using gffutils-cli (by default it gets printed into stdout.)
Suppose you want to pull out the Sox2 gene ENSMUSG00000074637:

```
$ grep ENSMUSG00000074637 small_ensembl_test.gff3
3	protein_coding	gene	34548927	34551383	.	+	.	ID=ENSMUSG00000074637;Name=Sox2;Name=ENSMUSG00000074637
3	protein_coding	mRNA	34548977	34551381	.	+	.	ID=ENSMUST00000169607;Name=Sox2-201;Parent=ENSMUSG00000074637
3	protein_coding	mRNA	34548927	34551383	.	+	.	ID=ENSMUST00000099151;Name=Sox2-001;Parent=ENSMUSG00000074637
3	protein_coding	mRNA	34549058	34551380	.	+	.	ID=ENSMUST00000171472;Name=Sox2-202;Parent=ENSMUSG00000074637
```

As we know, grepping does not work, since you'll only pull features that have the gene ID as a parent.
Let's sanitize the file, which will put it in a canonical order, fix potential coordinate issues,
and add a grep-able attribute.

```
$ gffutils-cli sanitize small_ensembl_test.gff3 > sanitized.gff3
Sanitizing GFF small_ensembl_test.gff3
- Loading GFF in memory
```

The sanitized file has a signature of gffutils:

```
$ head -n 2 sanitized.gff3
#GFF3 file (created by gffutils on 2013-05-11 12:55:28)
3	protein_coding	gene	34548927	34551383	.	+	.	ID=ENSMUSG00000074637;Name=Sox2,ENSMUSG00000074637;gid=ENSMUSG00000074637
```

It added the 'gid' optional attribute to the chilren of ENSMUSG00000074637. Now the file
is grep-able:

```
$ grep ENSMUSG00000074637 sanitized.gff3
3	protein_coding	gene	34548927	34551383	.	+	.	ID=ENSMUSG00000074637;Name=Sox2,ENSMUSG00000074637;gid=ENSMUSG00000074637
3	protein_coding	mRNA	34548927	34551383	.	+	.	ID=ENSMUST00000099151;Name=Sox2-001;Parent=ENSMUSG00000074637;gid=ENSMUSG00000074637
3	protein_coding	exon	34548927	34551383	.	+	.	ID=exon:ENSMUST00000099151:1;Parent=ENSMUST00000099151;gid=ENSMUSG00000074637
3	protein_coding	CDS	34549338	34550297	.	+	0	ID=CDS:ENSMUST00000099151:1;Parent=ENSMUST00000099151;gid=ENSMUSG00000074637
3	protein_coding	five_prime_UTR	34548927	34549337	.	+	.	ID=five_prime_UTR:ENSMUST00000099151:1;Parent=ENSMUST00000099151;gid=ENSMUSG00000074637
3	protein_coding	start_codon	34549338	34549340	.	+	0	ID=start_codon:ENSMUST00000099151:1;Parent=ENSMUST00000099151;gid=ENSMUSG00000074637
3	protein_coding	stop_codon	34550295	34550297	.	+	0	ID=stop_codon:ENSMUST00000099151:1;Parent=ENSMUST00000099151;gid=ENSMUSG00000074637
3	protein_coding	mRNA	34549058	34551380	.	+	.	ID=ENSMUST00000171472;Name=Sox2-202;Parent=ENSMUSG00000074637;gid=ENSMUSG00000074637
3	protein_coding	exon	34549058	34550408	.	+	.	ID=exon:ENSMUST00000171472:1;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
3	protein_coding	exon	34550440	34551380	.	+	.	ID=exon:ENSMUST00000171472:2;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
3	protein_coding	CDS	34549338	34550297	.	+	0	ID=CDS:ENSMUST00000171472:1;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
3	protein_coding	five_prime_UTR	34549058	34549337	.	+	.	ID=five_prime_UTR:ENSMUST00000171472:1;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
3	protein_coding	start_codon	34549338	34549340	.	+	0	ID=start_codon:ENSMUST00000171472:1;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
3	protein_coding	stop_codon	34550295	34550297	.	+	0	ID=stop_codon:ENSMUST00000171472:1;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
3	protein_coding	three_prime_UTR	34550440	34551380	.	+	.	ID=three_prime_UTR:ENSMUST00000171472:1;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
3	protein_coding	mRNA	34548977	34551381	.	+	.	ID=ENSMUST00000169607;Name=Sox2-201;Parent=ENSMUSG00000074637;gid=ENSMUSG00000074637
3	protein_coding	exon	34548977	34550405	.	+	.	ID=exon:ENSMUST00000169607:1;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	exon	34550440	34550508	.	+	.	ID=exon:ENSMUST00000169607:2;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	exon	34550596	34551381	.	+	.	ID=exon:ENSMUST00000169607:3;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	CDS	34549338	34550297	.	+	0	ID=CDS:ENSMUST00000169607:1;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	five_prime_UTR	34548977	34549337	.	+	.	ID=five_prime_UTR:ENSMUST00000169607:1;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	start_codon	34549338	34549340	.	+	0	ID=start_codon:ENSMUST00000169607:1;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	stop_codon	34550295	34550297	.	+	0	ID=stop_codon:ENSMUST00000169607:1;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	three_prime_UTR	34550440	34550508	.	+	.	ID=three_prime_UTR:ENSMUST00000169607:1;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	three_prime_UTR	34550596	34551381	.	+	.	ID=three_prime_UTR:ENSMUST00000169607:2;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	CDS	34549338	34550297	.	+	0	ID=CDS:ENSMUST00000099151:1;Parent=ENSMUST00000099151;gid=ENSMUSG00000074637
3	protein_coding	CDS	34549338	34550297	.	+	0	ID=CDS:ENSMUST00000169607:1;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	CDS	34549338	34550297	.	+	0	ID=CDS:ENSMUST00000171472:1;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
3	protein_coding	exon	34548927	34551383	.	+	.	ID=exon:ENSMUST00000099151:1;Parent=ENSMUST00000099151;gid=ENSMUSG00000074637
3	protein_coding	exon	34548977	34550405	.	+	.	ID=exon:ENSMUST00000169607:1;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	exon	34550440	34550508	.	+	.	ID=exon:ENSMUST00000169607:2;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	exon	34550596	34551381	.	+	.	ID=exon:ENSMUST00000169607:3;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	exon	34549058	34550408	.	+	.	ID=exon:ENSMUST00000171472:1;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
3	protein_coding	exon	34550440	34551380	.	+	.	ID=exon:ENSMUST00000171472:2;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
3	protein_coding	five_prime_UTR	34548927	34549337	.	+	.	ID=five_prime_UTR:ENSMUST00000099151:1;Parent=ENSMUST00000099151;gid=ENSMUSG00000074637
3	protein_coding	five_prime_UTR	34548977	34549337	.	+	.	ID=five_prime_UTR:ENSMUST00000169607:1;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	five_prime_UTR	34549058	34549337	.	+	.	ID=five_prime_UTR:ENSMUST00000171472:1;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
3	protein_coding	start_codon	34549338	34549340	.	+	0	ID=start_codon:ENSMUST00000099151:1;Parent=ENSMUST00000099151;gid=ENSMUSG00000074637
3	protein_coding	start_codon	34549338	34549340	.	+	0	ID=start_codon:ENSMUST00000169607:1;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	start_codon	34549338	34549340	.	+	0	ID=start_codon:ENSMUST00000171472:1;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
3	protein_coding	stop_codon	34550295	34550297	.	+	0	ID=stop_codon:ENSMUST00000099151:1;Parent=ENSMUST00000099151;gid=ENSMUSG00000074637
3	protein_coding	stop_codon	34550295	34550297	.	+	0	ID=stop_codon:ENSMUST00000169607:1;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	stop_codon	34550295	34550297	.	+	0	ID=stop_codon:ENSMUST00000171472:1;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
3	protein_coding	three_prime_UTR	34550440	34550508	.	+	.	ID=three_prime_UTR:ENSMUST00000169607:1;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	three_prime_UTR	34550596	34551381	.	+	.	ID=three_prime_UTR:ENSMUST00000169607:2;Parent=ENSMUST00000169607;gid=ENSMUSG00000074637
3	protein_coding	three_prime_UTR	34550440	34551380	.	+	.	ID=three_prime_UTR:ENSMUST00000171472:1;Parent=ENSMUST00000171472;gid=ENSMUSG00000074637
```

Note that GFFWriter writes the entries in a standard order: the gene, the first mRNA,
the mRNA's exons, the remaining mRNA's children, the second mRNA, etc.

You can optionally sanitize the file in place, which is very helpful if you want
to standardize a bunch of GFF files without creating more files.

```
$ gffutils-cli sanitize small_ensembl_test.gff3 --in-place
Sanitizing GFF small_ensembl_test.gff3
- Loading GFF in memory
- Sanitizing file in place (overwriting current file)
```

```
$ head small_ensembl_test.gff3
#GFF3 file (created by gffutils on 2013-05-11 12:57:58)
3	protein_coding	gene	34548927	34551383	.	+	.	ID=ENSMUSG00000074637;Name=Sox2,ENSMUSG00000074637;gid=ENSMUSG00000074637
3	protein_coding	mRNA	34548927	34551383	.	+	.	ID=ENSMUST00000099151;Name=Sox2-001;Parent=ENSMUSG00000074637;gid=ENSMUSG00000074637
3	protein_coding	exon	34548927	34551383	.	+	.	ID=exon:ENSMUST00000099151:1;Parent=ENSMUST00000099151;gid=ENSMUSG00000074637
3	protein_coding	CDS	34549338	34550297	.	+	0	ID=CDS:ENSMUST00000099151:1;Parent=ENSMUST00000099151;gid=ENSMUSG00000074637
```

Sanitize is inherently "gene-centric". It does not really make that much
sense if you use GFF files to store a bunch of features that don't have
gene records are parents -- e.g. if your GFF is just a list of CDS regions,
or a bunch of mRNAs. This is in part why I separated it from "clean_gff".

2. Revised fetch() in gffutils-cli: completed the fetch() function
to fetch IDs from a GFF file. It skips entries that are not found
instead of aborting.

3. Placeholders for:

- convert: command-line function to convert GTF <-> GFF

- annotate: command-line function to annotate a GFF file with entries
from another GFF file. Not sure if this is globally useful so I didn't put in
the implementation, but this is the kind of operation I use a lot: you have some
GFF file, and you want to annotate (e.g. by overlap) its features with features
from another GFF file. For example, you have a set of exons that meet some criteria
in one GFF, and a standard gene models annotation in another GFF. You want to then
output a GFF containing the exons along with the gene that they fall in based on
coordinate overlap in the second GFF. This can be done elegantly thanks to pybedtools.
The question is whether this is generic enough to be included, and whether you want to
add pybedtools as a dependency.

I have only tested my features with GFF files, not GTF, but it should theoretically
work identically for either format.

Some questions and misc. notes:

- GFF db creation is quite fast when using the ":memory:" features. However, when
I try to create dbs on disk with large-ish GFF files, it still takes a while. For example,
making a db for Ensembl gene models in the mouse genome (NCBIM37.65 release) takes about 15 minutes"

```
$ time gffutils-cli create Mus_musculus.NCBIM37.65.gff
Populating features table and first-order relations: 0 featuresatures

real	15m4.841s
user	6m42.829s
sys	0m32.678s
```

The file has ~1.5 million lines. Unfortunately this is going to be a very common
use case; it's just a gene annotation for the mouse genome, and the file is large only
because the GFF format is extremely verbose, making it potentially too large to load into memory.

Is it possible to speed this up? Any thoughts on alternative ways of dealing with these files?
I uploaded the file here:

http://genes.mit.edu/yarden/Mus_musculus.NCBIM37.65.gff

- Most functions that take a GFF input have to handle a plain GFF file (and create a database for it,
optionally in memory) or a GFF database. Handling all the cases can get a bit clunky and repetitive
if you have to do for each function. For example, in sanitize_gff_file(), I wrote this to
handle the argument "gff_fname" (which would be a database filename too):

```
if is_gff_db(gff_fname):
# It's a database filename, so load it
db = gffutils.FeatureDB(gff_fname)
else:
# Need to create a database for file
if in_memory:
db = gffutils.create_db(gff_fname, ":memory:",
verbose=False)
else:
db = get_gff_db(gff_fname)
```

The logic is: If the file is a DB, load it. If it isn't, create a db for it on the fly,
optionally in memory. It'd be helpful to come up with a standardized function/class that more uniformly
handles this.

- What is the right way to detect "orphan" features? By orphans I mean features that either
have no parent, or that have a parent whose ID is not present in the current GFF.
There are several functions (like "sanitize") where these features, if present, could be
dropped because they have no parent. It'd be very useful to be able to fetch these features
so that they can be specially handled. For example, you might sanitize all the features
that are within gene units, and then want to still output the orphan features along with
the sanitized file so that they are not lost.

- I tried to follow Ryan's coding conventions as much as I possible. There were
a few cases where I wasn't sure what the convention should be, like

- Since in-memory dbs are fast, is it easy to modify the argparser for gffutils-cli to optionally
read GFF files from stdin and not from a filename? This might be a bit tricky
since it only makes sense when the input and outputs are plain text GFFs and not
binary .db files.

- Finally, `"test.py:test_update()"` fails on my end, not sure if that's expected.