Stdpopsim: annotations: what filtering and interval operations should we support?

[grahamgower]

Assuming we have annotations as a set of intervals, what sorts of set operations do we really need to do on these intervals? See #560 for background. Can we get away with not having another dependency (pybedtools doesn't support windows, pyranges loads files very slowly)? I was curious whether naiive interval operations are sufficient for our use cases, so I threw together some example code that extracts ensembl exons for chr1, finds their intersection, and takes the complement (all intervals that aren't exons):
https://gist.github.com/grahamgower/c193443f88d0bdd4f3e52a1185393679

From the whole-genome gff.

$ time ./gff.py Homo_sapiens.GRCh37.87.gff3.gz
loading gff took 3.4992711640661582 seconds
getting chromosome extent took 0.005866469000466168 seconds
getting exons took 0.4508029669523239 seconds
taking complement took 0.5587347559630871 seconds

real    0m4,714s
user    0m4,540s
sys     0m0,170s

From chr1-only gff.

$ time ./gff.py Homo_sapiens.GRCh37.87.chromosome.1.gff3.gz
loading gff took 0.47997945395763963 seconds
getting chromosome extent took 0.0058168950490653515 seconds
getting exons took 0.4441574700176716 seconds
taking complement took 0.5427950490266085 seconds

real    0m1,672s
user    0m1,531s
sys     0m0,140s

[jeromekelleher]

An intermediate approach @grahamgower would be to use one of the generic Python interval libraries to implement the interval ops. Have you considered any of these?

[grahamgower]

I only took a quick look, but I'm not sure how appropriate they are for integer intervals. I played with ncls a little, which is used by pyranges. Its very fast indeed and has no external dependencies, but it would require a bit of digging as there's no documentation (I can't even guess what the functions do from their names alone).

[mufernando]

This looks great, @grahamgower!

I don't think we would need more than intersect and complement (and subtract, which is just these two operations nested together). This would suffice if we only allowed two complementary annotation sets (exon vs. non-exon). We should discuss this more in-depth in one of the calls, but it was my impression we wanted more than two annotations, e.g. exons, introns, UTR and intergenic. If this is the case, then we would need to apply these basic operations you implemented over intervals (it is my understanding as it is they only work one interval vs array of intervals, and not array of intervals vs array of intervals).

[mufernando]

I think we should use pyranges to do the interval operations. It's true that their utilities for reading files were slow, but given that we bypassed this with #560, I see no issue in using it for the interval operations. I don't see why we should reinvent the wheel here.

[petrelharp]

From pyrange's website:

PyRanges is in a beta state until any possible issues with the move to pandas 1.0 have been ironed out.

I don't know, this makes me kinda nervous? (both "beta" and "pandas")

[grahamgower]

I'm in two minds about this. If we use an external library (which I agree is preferred), we want something that is mature, performant, and available on all platforms. pyranges looks very nice, but is far from mature. Perhaps we can use it in the interim, and see how it goes?

[andrewkern]

looks like pyranges will let us create an object from a dict, so we can very easily read in the zarr from #560 as a PyRange object.

currently i'm shoving that info into a pandas.DataFrame