Queries on Interval Data

TL;DR

When working with genomic data, you often need to perform queries on intervals.

For example, you might want to find peaks or SNVs that overlap with gene regions, or identify histone marks within open chromatin regions.

Naively performing such queries results in $O(N^2)$ complexity. In practice, this computational cost is far from ideal.

Solutions

CLI

CLI tools like bedtools can handle these queries reasonably efficiently. They generally achieve around $O(Nlog(N))$ complexity. If performance is not a concern, using these well-established tools is a good approach.

• bedtools
• bedops
• bedtk

Since these tools are well-tested, they provide peace of mind for routine operations. Another advantage is that you don't have to write your own parsers for formats like bed. Among them, bedtools is the most popular and feature-rich, while bedtk appears to be the fastest.

Building Your Own

When CLI tools cannot solve the problem, you'll need to write your own code. Data structures suitable for this purpose have naturally been studied for a long time. While the computational complexity in terms of order hasn't changed much, benchmarks show significant differences, suggesting that optimizations have been steadily improving.

The fundamental ideas are based on data structures such as Interval Tree and R-tree.

For implementations of these basic data structures, there are various options including the rust-bio data structures and a Python implementation. Since these are implemented as AVL trees, using them is much less effort than implementing your own.

More recently developed data structures that are faster and more memory-efficient include the following. cgranges is the data structure used by bedtk.

Based on the benchmarks published on GitHub, COITree is the fastest. Personally, since I have very little experience with C, I would probably just use COITree.