lh3 / Cgranges

Introduction

cgranges is a small C library for genomic interval overlap queries: given a genomic region r and a set of regions R, finding all regions in R that overlaps r. Although this library is based on interval tree, a well known data structure, the core algorithm of cgranges is distinct from all existing implementations to the best of our knowledge. Specifically, the interval tree in cgranges is implicitly encoded as a plain sorted array (similar to binary heap but packed differently). Tree traversal is achieved by jumping between array indices. This treatment makes cgranges very efficient and compact in memory. The core algorithm can be implemented in ~50 lines of C++ code, much shorter than others as well. Please see the code comments in cpp/IITree.h for details.

Usage

Test with BED coverage

For testing purposes, this repo implements the bedtools coverage tool with cgranges. The source code is located in the test/ directory. You can compile and run the test with:

cd test && make
./bedcov-cr test1.bed test2.bed

The first BED file is loaded into RAM and indexed. The depth and the breadth of coverage of each region in the second file is computed by query against the index of the first file.

The test/ directory also contains a few other implementations based on IntervalTree.h in C++, quicksect in Cython and ncls in Cython. The table below shows timing and peak memory on two test BEDs available in the release page. The first BED contains GenCode annotations with ~1.2 million lines, mixing all types of features. The second contains ~10 million direct-RNA mappings. Time1a/Mem1a indexes the GenCode BED into memory. Time1b adds whole chromosome intervals to the GenCode BED when indexing. Time2/Mem2 indexes the RNA-mapping BED into memory. Numbers are averaged over 5 runs.

Here, IAITree = implicit augmented interval tree, used by cgranges; CITree = centered interval tree, used by Erik Garrison's IntervalTree; AIList = augmented interval list, by Feng et al; NCList = nested containment list, taken from ncls by Feng et al; AITree = augmented interval tree, from kerneltree. "Cov" indicates whether the program calculates breadth of coverage. Comments:

• AIList keeps start and end only. IAITree and CITree addtionally store a 4-byte "ID" field per interval to reference the source of interval. This is partly why AIList uses the least memory.

• IAITree is more sensitive to the worse case: the presence of an interval spanning the whole chromosome.

• IAITree uses an efficient radix sort. CITree uses std::sort from STL, which is ok. AIList and NCList use qsort from libc, which is slow. Faster sorting leads to faster indexing.

• IAITree in C++ uses identical core algorithm to the C version, but limited by its APIs, it wastes time on memory locality and management. CITree has a similar issue.

• Computing coverage is better done when the returned list of intervals are start sorted. IAITree returns sorted list. CITree doesn't. Not sure about others. Computing coverage takes a couple of seconds. Sorting will be slower.

• Printing intervals also takes a noticeable fraction of time. Custom printf equivalent would be faster.

• IAITree+Cython is a wrapper around the C version of cgranges. Cython adds significant overhead.

• Bedtools is designed for a variety of applications in addition to computing coverage. It may keep other information in its internal data structure. This micro-benchmark may be unfair to bedtools.

• In general, the performance is affected a lot by subtle implementation details. CITree, IAITree, NCList and AIList are all broadly comparable in performance. AITree is not recommended when indexed intervals are immutable.

Use cgranges as a C library

cgranges_t *cr = cr_init(); // initialize a cgranges_t object
cr_add(cr, "chr1", 20, 30, 0); // add a genomic interval
cr_add(cr, "chr2", 10, 30, 1);
cr_add(cr, "chr1", 10, 25, 2);
cr_index(cr); // index

int64_t i, n, *b = 0, max_b = 0;
n = cr_overlap(cr, "chr1", 15, 22, &b, &max_b); // overlap query; output array b[] can be reused
for (i = 0; i < n; ++i) // traverse overlapping intervals
	printf("%d\t%d\t%d\n", cr_start(cr, b[i]), cr_end(cr, b[i]), cr_label(cr, b[i]));
free(b); // b[] is allocated by malloc() inside cr_overlap(), so needs to be freed with free()

cr_destroy(cr);

Use IITree as a C++ library

IITree<int, int> tree;
tree.add(12, 34, 0); // add an interval
tree.add(0, 23, 1);
tree.add(34, 56, 2);
tree.index(); // index
std::vector<size_t> a;
tree.overlap(22, 25, a); // retrieve overlaps
for (size_t i = 0; i < a.size(); ++i)
	printf("%d\t%d\t%d\n", tree.start(a[i]), tree.end(a[i]), tree.data(a[i]));