yarden
commented
11 years ago I've been thinking about two related GFF functionalities, which I believe reflect very common pattern usages for other users too, but unfortunately no tool that I know of has these features. I'm interested in: (1) an efficient library for parsing GFF files under the hood, from Python, which ideally is dependency-free (see more on this below). (2) a command line utility for quickly "grep"-ing, slicing and viewing GFF files.
A bit of background on what I'm trying to achieve and my current issues with GFF3: I generally think that GFF3 is not a good format, and so tools like #1 and #2 are needed to make it workable, whether you're parsing GFF files in your code or just using it to store / browse annotations. A primary problem with GFF3 is that even though it's a hierarchical format, there isn't a single identifier running through each gene tree, and this combined with the fact that the entries are unordered, makes it a pain to work with. Rather than doing multiple passes through the file to collect all the relevant records, I've usually defaulted on an inelegant and inefficient approach of loading all the genes into memory, which is suboptimal and obviously does not scale. (As an example, I include a script of mine below that uses the current library to add introns to a GFF. With the right GFF parsing library, this could code be much improved.)
Related to this, GFF is missing a random access indexing format, and so I wrote my own extremely primitive code (https://github.com/yarden/MISO/blob/fastmiso/misopy/index_gff.py) to at least make it so you can easily fetch an entire gene "unit" (all of the gene's records) by the ID= of the gene. Since the gene entries of GFF are independent, and since the problem I was working on is embarrassingly parallelizable (each gene can be processed independently on a cluster system), I needed a way to instantly retrieve a parsed gene unit in a way that wouldn't cause deadlock for multiple threads. I'm generally unhappy with the code and would like to find a library that replaces my current clunky solution. A library that uses C (or Cython) but has a Python wrapper sounds like the best solution to me. I'm hoping that this library could be used in my project (https://github.com/yarden/MISO) instead of the GFF parser.
Related to the library, I'd like to find a tool that allows slicing/viewing/retrieval of GFF3 files. I really liked the way you framed the core tasks on your gffutils README, and I was imagining the same set of features, except in a form also accessible from the command line. If I combine your requirements and mine, I think the core features would be something like:
Feature 1: A "grep" like feature for gene units:
# Return all child records of 'mygene' (including 'mygene' "gene" entry). The gene entry itself,
# the mRNAs, children of the mRNAs and children of the exons if available
gffutils --fetch-gene mygene mygff.gff3
# Return all child records of 'mygeneA', 'mygeneB', 'mygeneC'
gffutils --fetch-gene mygeneA,mygeneB,mygeneC mygff.gff3 Ideally, there would be a 'limit' feature (related to BCBio parser) that allows you to "clip" the hierarchy somewhere. For example:
# Return all mRNA child records of 'mygene'; don't return anything "below" mRNA in the hierarchy
gffutils --fetch-gene mygene mygff.gff3 --limit mRNA
# Return all child records of 'mygene' including mRNAs and exons, but nothing below
# (e.g. don't return CDSS or UTR annotations if those are given as children of exons)
gffutils --fetch-gene mygene mygff.gff3 --limit exonHowever, this can also be achieved by a feature analogous to grep -V that excludes out entries:
# Get all children records of a gene except CDSS and 3p_UTR, 5p_UTR
gffutils --fetch-gene mygene | gff-utils --exclude-types cdss,3p_utr,5p_utrFeature 2: Retrieve all exons common to each gene's mRNAs
gffutils --get-common-exons mygff.gff3(as you described in your README.)
Feature 3: This one is less essential in my view, but it'd be nice to retrieve entries by region coordinates:
# Get all entries on chr5:500-50000
gffutils --fetch-region chr5:500-50000 mygff.gff3In any case, combining features #1 and #2, you can easily now do things per gene or set of genes, which I think would be very helpful:
# Retrieve all exons common exons only for geneA and geneB
gffutils --fetch-gene mygeneA,mygeneB | gffutils --get-common-exons - Feature 4: Trivial but useful, a kind of "sanitize" or clean up (like your current "clean_gff") utility that checks the GFF for incomplete entries and fixes things like start > end. It would allow convenient piping like this:
# Sanitize a malformed GFF (that has start > end) and then retrieve entries for "mygene" from it
gffutils --clean bad.gff3 | gffutils --fetch-gene mygene - In summary, the ideal library would allow one to programmatically do the features required to achieve Features #1-#4. The only thing that (for my use case at least) would be indispensable in a library, but unrelated to the command line, is the random access index for GFF. I'd love to be able to run distributed code using the library and get random access to very large GFF files this way. The way I do it now is again pretty inelegant, using index_gff.py. It cannot result in deadlock on a distributed system as mentioned, but the disadvantage is that it creates a proliferation of files and probably additional I/O overhead. However, the indexing cost is minimal (takes ~10 mins per large GFF) and I think any user would be happy to pay this minimal startup cost to build an index if it later on allows random access by gene and/or coordinate regions.
Also, I definitely agree that your idea that universal GFF <-> GTF <-> ... etc conversions would be very useful. One approach to do this is to pick the format that is internally best (I think that might be GTF, since it might be the most "explicit") and just internally convert all inputs to that format. For example, if a GFF file is passed, it can be internally converted to GTF and represented that way and then optionally serialized as GFF or GTF.
Best, --Yarden
P.S. @chapmanb 's BCBio GFF parser (http://bcbio.wordpress.com/2009/03/08/initial-gff-parser-for-biopython/) looks really nice from what I can tell, but the dependency on BioPython makes it unusable for me as a library, since I'd like to package this with MISO and users are already complaining about the various Python dependencies (in fact, most people have trouble installing scipy and I'm working to remove this requirement - unfortunate since scipy is valuable.) I also don't use BioPython for anything else at the moment, so wouldn't want to tack it on as a dependency just for GFF parsing abilities.
##
## Add introns to GFF file
##
import misopy
import misopy.gff_utils as gff_utils
import misopy.Gene as gene_utils
def add_introns_to_gff(gff_filename, output_dir):
"""
Add 'intron' entries to GFF.
"""
output_basename = \
utils.trim_gff_ext(os.path.basename(gff_filename))
output_filename = \
os.path.join(output_dir,
"%s.with_introns.gff" %(output_basename))
print "Adding introns to GFF..."
print " - Input: %s" %(gff_filename)
print " - Output: %s" %(output_filename)
gff_out = gff_utils.Writer(open(output_filename, "w"))
gff_db = gff_utils.GFFDatabase(from_filename=gff_filename,
reverse_recs=True)
t1 = time.time()
genes = gene_utils.load_genes_from_gff(gff_filename)
for gene_id in genes:
gene_info = genes[gene_id]
gene_tree = gene_info["hierarchy"]
gene_obj = gene_info["gene_object"]
gene_rec = gene_tree[gene_id]["gene"]
# Write the GFF record
gff_out.write(gene_rec)
# Write out the mRNAs, their exons, and then
# input the introns
for mRNA_id in gene_tree[gene_id]["mRNAs"]:
curr_mRNA = gene_tree[gene_id]["mRNAs"][mRNA_id]
gff_out.write(curr_mRNA["record"])
# Write out the exons
curr_exons = gene_tree[gene_id]["mRNAs"][mRNA_id]["exons"]
for exon in curr_exons:
gff_out.write(curr_exons[exon]["record"])
# Now output the introns
for isoform in gene_obj.isoforms:
intron_coords = []
for first_exon, second_exon in zip(isoform.parts,
isoform.parts[1::1]):
# Intron start coordinate is the coordinate right after
# the end of the first exon, intron end coordinate is the
# coordinate just before the beginning of the second exon
intron_start = first_exon.end + 1
intron_end = second_exon.start - 1
if intron_start >= intron_end:
continue
intron_coords.append((intron_start, intron_end))
# Create record for this intron
intron_id = "%s:%d-%d:%s" %(gene_obj.chrom,
intron_start,
intron_end,
gene_obj.strand)
intron_rec = \
gff_utils.GFF(gene_obj.chrom, gene_rec.source, "intron",
intron_start, intron_end,
attributes={"ID": [gene_obj.label],
"Parent": [isoform.label]})
gff_out.write(intron_rec)
t2 = time.time()
print "Addition took %.2f minutes." %((t2 - t1)/60.)
chapmanb
daler
peterjc
What sort of functionality should be exposed to the CLI?
For example a universal GTF<->GFF<->refFlat conversion script would be useful.