快速查询区间overlap的C/C++库：cgranges

前言

本库的核心数据结构是隐式区间树（implicit interval trees），本人之前写过一篇介绍，可看这里 说是数据结构，其实只是一种遍历策略 隐式区间树示例：（i=9的节点的Max打印错了，应该是450）

cgranges

源代码在这里（cgranges.h、cgranges.c） 适用于生物学上的区间查询（genomic interval overlap queries） 给定一系列区间R与查询区间r，可以返回R中与r有overlap的所有区间

一、理论基础

理论原文为Bedtk: finding interval overlap with implicit interval tree（原文很短小，理论核心可以说半页都不到，推荐大家都去看看）

二、数据结构主体 cgranges_t：主体结构，主要记录全局信息typedef struct { int64_t n_r, m_r; // number and max number of intervals cr_intv_t *r; // list of intervals (of size _n_r_) int32_t n_ctg, m_ctg; // number and max number of contigs cr_ctg_t *ctg; // list of contigs (of size _n_ctg_) void *hc; // dictionary for converting contig names to integers } cgranges_t; r：存储区间数组。整体上按contig（染色体）顺序排列，同一个contig下按区间起点排序(隐式区间树)ctg：存储contig数组。contig顺序按加入的顺序排列hc：存储contig名称与内部id的映射。实际是个用khash模块实例化的map<string, int> cr_intv_t：区间结构体，记录区间信息typedef struct { // an interval uint64_t x; // prior to cr_index(), x = ctg_id<<32|start_pos; after: x = start_pos<<32|end_pos uint32_t y:31, rev:1; int32_t label; // NOT used } cr_intv_t; x：区间坐标。建索引前后不一样：调用前为ctg_id<<32 | start_pos，调用后为start_pos<<32 | end_posy：建索引前后不一样：调用前为end_pos，调用后为额外信息当前子树的最大区间终点maxrev：标记区间是否反向（未使用） cr_ctg_t：contig（染色体）结构体，记录contig相关信息typedef struct { // a contig char *name; // name of the contig int32_t len; // max length seen in data int32_t root_k; int64_t n, off; // sum of lengths of previous contigs } cr_ctg_t; len：该contig已知的最大坐标（最大区间终点）root_k：该contig根节点的层级Kn：属于该contig的区间数off：属于该contig的区间在区间数组中的起始偏移 三、API操作 1）初始化/释放：cr_init、cr_destroy cgranges_t *cr_init(void); void cr_destroy(cgranges_t *cr); 2）添加区间：cr_add cr_add：主API，添加一个区间 cr：经初始化的主体结构ctg：contig名称（如“chr1”）st、en：区间起点、终点（从if (st > en) return 0;来看应该都是闭区间）label_int：label（未使用） cr_intv_t *cr_add(cgranges_t *cr, const char *ctg, int32_t st, int32_t en, int32_t label_int) { cr_intv_t *p; int32_t k; if (st > en) return 0; k = cr_add_ctg(cr, ctg, 0); if (cr->n_r == cr->m_r) EXPAND(cr->r, cr->m_r); p = &cr->r[cr->n_r++]; p->x = (uint64_t)k << 32 | st; p->y = en; p->label = label_int; if (cr->ctg[k].len < en) cr->ctg[k].len = en; return p; } cr_add_ctg：添加一个contig到hash表中（该hash表操作细节可看这里） 用kh_put尝试添加，如果没有key就先插入新的，然后更新contig已知最大坐标len // Add a contig and length. Call this for desired contig ordering. _len_ can be 0. int32_t cr_add_ctg(cgranges_t *cr, const char *ctg, int32_t len);

3）建索引：cr_index void cr_index(cgranges_t *cr) { int32_t i; cr_index_prepare(cr); for (i = 0; i < cr->n_ctg; ++i) cr->ctg[i].root_k = cr_index1(&cr->r[cr->ctg[i].off], cr->ctg[i].n); } 准备：先排序所有区间，然后遍历两遍所有区间：第一遍记录contig信息，第二遍把x更新为start_pos<<32 | end_pos，并把y置0（为了腾出来位置存额外信息max）void cr_index_prepare(cgranges_t *cr) { int64_t i, st; if (!cr_is_sorted(cr)) cr_sort(cr); for (st = 0, i = 1; i <= cr->n_r; ++i) { if (i == cr->n_r || cr->r[i].x>>32 != cr->r[st].x>>32) { ... } } for (i = 0; i < cr->n_r; ++i) { cr_intv_t *r = &cr->r[i]; r->x = r->x<<32 | r->y; r->y = 0; } } 构建：更新所有节点的额外信息子树包含的最大区间终点max（因为隐式区间树的这种节点间特殊的跳转关系，可看到每层节点间的跳转也非常有规律） 先把所有叶子节点for (i = 0; i < n; i += 2) ...的额外信息y进行更新（就等于区间终点） last_i、last：记录的是上一层最后节点的索引及其max，是为了处理右子节点为空的情况 然后一层一层往上进行更新 k：当前层级x：当前层节点与子节点间的距离i0：当前层的第一个节点step：当前层节点间的距离last_i、last：更新为当前层的最后节点的索引和max 最后返回当前contig的总层级K int32_t cr_index1(cr_intv_t *a, int64_t n) { ... for (i = 0; i < n; i += 2) last_i = i, last = a[i].y = (int32_t)a[i].x; for (k = 1; 1LL<<k <= n; ++k) { int64_t x = 1LL<<(k-1), i0 = (x<<1) - 1, step = x<<2; for (i = i0; i < n; i += step) { ... a[i].y = e; } last_i = last_i>>k&1? last_i - x : last_i + x; if (last_i < n && a[last_i].y > last) last = a[last_i].y; } return k - 1; }

4）排序：cr_sort

cgranges.c文件内部实现了一个“基数排序”的模块（与C通用库Klib的其他模块一样即插即用）

#define cr_intv_key(r) ((r).x) KRADIX_SORT_INIT(cr_intv, cr_intv_t, cr_intv_key, 8) void cr_sort(cgranges_t *cr) { if (cr->n_ctg == 0 || cr->n_r == 0) return; radix_sort_cr_intv(cr->r, cr->r + cr->n_r); } 基数排序：KRADIX_SORT_INIT #define RS_MIN_SIZE 64 #define KRADIX_SORT_INIT(name, rstype_t, rskey, sizeof_key) \ ... \ void radix_sort_##name(rstype_t *beg, rstype_t *end) { \ if (end - beg <= RS_MIN_SIZE) rs_insertsort_##name(beg, end); \ else rs_sort_##name(beg, end, RS_MAX_BITS, (sizeof_key - 1) * RS_MAX_BITS); \ }

虽说是基数排序，实际在要排序的元素个数少于等于RS_MIN_SIZE时会使用插入排序，个数较多时才用基数排序。插入排序我就不过多赘述了，这里重点是后者。

内部数据结构

主要的数据结构就是下面的结构体，可理解为一个“桶”，内部只有两个指针，一个指向桶包含的元素起始位置，另一个指向终点位置的后一个（经典的左闭右开的设计，两个指针相等时表示不包含元素）

#define KRADIX_SORT_INIT(name, rstype_t, rskey, sizeof_key) \ typedef struct { \ rstype_t *b, *e; \ } rsbucket_##name##_t; ... rstype_t：元素类型（比如区间类型cr_intv_t）rskey：获取元素key函数（比如取区间的起始位置宏cr_intv_key，注意此时key为ctg_id << 32 | start_pos）sizeof_key：元素key的字节数 算法流程

我今天才算是知道基数排序有两种模式：先排小位后排大位的LSD模式(Least Significant Digit 最低有效位)，和先大后小的MSD(Most Significant Digit 最高有效位)模式 LSD在菜鸟上有原理讲解、动图展示和实际代码，具体可看这里，下面讲解一下MSD原理

主体思路：（可以说跟快速排序很像，都属于分治法）每次排序时根据当前最高有效位的大小把所有元素分配到各个“桶”中，然后递归地排序每个“桶”（原地操作），这样“桶”的顺序是有序的，每个“桶”中的顺序也是有序的，最后就能保证整个数组有序参数说明：beg、end为要排序的区间起止位置，n_bits为排序时考虑的bit位数（就是MSD，关系到桶的大小），s为当前MSD的偏移位置关键变量：m为掩码，b为桶集合，be为桶的结束位置，k与l为两个桶指针，指向某个具体的桶算法流程： 第一个for循环初始化每个“桶”第二个for循环遍历一次数组，计算每个元素对应到“桶”的索引位置，把对应桶的计数进行更新（计数是利用桶的e与b之间的距离进行记录的）第三个for循环再次遍历桶，把所有的“桶”首尾“连接”起来（此时，可以预见的是第一个“桶”的b等于beg，最后一个桶的e等于end）此时相当于给每个元素都找好了对应的桶的位置，桶也准备好等元素放进来，就等着给元素换位置了第四个for循环就是对元素进行位置交换 k从第一个桶开始，往后遍历每个“桶”，当前“桶”还有元素需要交换时就继续处理l负责记录当前元素需要交换的桶位置tmp、swap与l->b合力完成元素交换，tmp存储下一个要处理的元素 经过上一步，每个“桶”的b == e，为了对每个桶进行分治法，需要把所有的“桶”重新“连接”起来，所以第五个for循环就是重置每个桶的b至此，元素都已分配到对应的“桶”中，减少偏移s，第六个for循环就是用下一个MSD对每个“桶”进行递归排序（元素个数较少时依旧使用插入排序） #define KRADIX_SORT_INIT(name, rstype_t, rskey, sizeof_key) \ ... \ void rs_sort_##name(rstype_t *beg, rstype_t *end, int n_bits, int s) { \ int size = 1<<n_bits, m = size - 1; \ rsbucket_##name##_t *k, b[1<<RS_MAX_BITS], *be = b + size; \ ... \ for (k = b; k != be; ++k) k->b = k->e = beg; \ for (rstype_t *i = beg; i != end; ++i) ++b[rskey(*i)>>s&m].e; \ for (k = b + 1; k != be; ++k) \ k->e += (k-1)->e - beg, k->b = (k-1)->e; \ for (k = b; k != be;) { \ if (k->b != k->e) { \ rsbucket_##name##_t *l; \ if ((l = b + (rskey(*k->b)>>s&m)) != k) { \ rstype_t tmp = *k->b, swap; \ do { \ swap = tmp; tmp = *l->b; *l->b++ = swap; \ l = b + (rskey(tmp)>>s&m); \ } while (l != k); \ *k->b++ = tmp; \ } else ++k->b; \ } else ++k; \ } \ for (b->b = beg, k = b + 1; k != be; ++k) k->b = (k-1)->e; \ if (s) { \ s = s > n_bits? s - n_bits : 0; \ for (k = b; k != be; ++k) \ if (k->e - k->b > RS_MIN_SIZE) rs_sort_##name(k->b, k->e, n_bits, s); \ else if (k->e - k->b > 1) rs_insertsort_##name(k->b, k->e); \ } \ } 5）查询overlap：cr_overlap

调用cr_get_ctg获取contig的ID后进行查询

int64_t cr_overlap(const cgranges_t *cr, const char *ctg, int32_t st, int32_t en, int64_t **b_, int64_t *m_b_) { return cr_overlap_int(cr, cr_get_ctg(cr, ctg), st, en, b_, m_b_); }

查询过程：

从该contig的根节点开始，用 栈+迭代 的方式遍历隐式区间树如果当前节点层数<= 3，即子树较小，则直接线性检查当前子树的所有节点如果当前节点层数较高，则按中序遍历顺序检查左子节点、当前节点、右子节点返回overlap结果个数，入参b_存储结果里每个overlap的索引，入参m_b_存储数组b_的容量

注1：因为b_指向的空间是API内部申请出来的，所以调用完API及处理后记得释放！ 注2：让我感觉有点怪的是：输出的条件是st < cr_en(&r[z.x]) && cr_st(&r[i]) < en，看起来查询区间是开区间，这与前面存储的闭区间不同，按理说如果存的是闭区间那么查询应该也是闭区间，那这样判定条件应该加上=才对，所以我感觉这里不太懂为啥这样设计

int64_t cr_overlap_int(const cgranges_t *cr, int32_t ctg_id, int32_t st, int32_t en, int64_t **b_, int64_t *m_b_) { ... c = &cr->ctg[ctg_id]; r = &cr->r[c->off]; p = &stack[t++]; p->k = c->root_k, p->x = (1LL<<p->k) - 1, p->w = 0; // push the root into the stack while (t) { // stack is not empyt istack_t z = stack[--t]; if (z.k <= 3) { // the subtree is no larger than (1<<(z.k+1))-1; do a linear scan int64_t i, i0 = z.x >> z.k << z.k, i1 = i0 + (1LL<<(z.k+1)) - 1; if (i1 >= c->n) i1 = c->n; for (i = i0; i < i1 && cr_st(&r[i]) < en; ++i) if (st < cr_en(&r[i])) { ... } } else if (z.w == 0) { // if left child not processed int64_t y = z.x - (1LL<<(z.k-1)); p = &stack[t++]; p->k = z.k, p->x = z.x, p->w = 1; // push current node back to the stack if (y >= c->n || r[y].y > st) { p = &stack[t++]; p->k = z.k - 1, p->x = y, p->w = 0; // push the left child to the stack } } else if (z.x < c->n && cr_st(&r[z.x]) < en) { if (st < cr_en(&r[z.x])) { // then z.x overlaps the query; write to the output array ... } p = &stack[t++]; p->k = z.k - 1, p->x = z.x + (1LL<<(z.k-1)), p->w = 0; // push the right child } } *b_ = b, *m_b_ = m_b; return n; } 6）查看指定索引区间：cr_start、cr_end

用例可看最后的代码示例

前两个用于知道区间指针时使用后两个用于知道区间索引位置i时使用 // retrieve start and end positions from a cr_intv_t object int32_t cr_st(const cr_intv_t *r) { return (int32_t)(r->x>>32); } int32_t cr_en(const cr_intv_t *r) { return (int32_t)r->x; } int32_t cr_start(const cgranges_t *cr, int64_t i) { return cr_st(&cr->r[i]); } int32_t cr_end(const cgranges_t *cr, int64_t i) { return cr_en(&cr->r[i]); } // NOT used int32_t cr_label(const cgranges_t *cr, int64_t i) { return cr->r[i].label; } 7）查询contig的ID：cr_get_ctg

直接用khash模块接口kh_get、kh_val获取指定contig的ID

int32_t cr_get_ctg(const cgranges_t *cr, const char *ctg) { khint_t k; strhash_t *h = (strhash_t*)cr->hc; k = kh_get(str, h, ctg); return k == kh_end(h)? -1 : kh_val(h, k); } 8）未测试的API cr_contain：获取完全被包含在区间[st, en]中的所有区间 先用二分法获取起始位置>= st的区间的索引然后往后直接顺序遍历，找出所有满足条件的区间当起始位置>= en时退出（后续不会有结果，因为是按起始位置递增排列的） int64_t cr_contain(const cgranges_t *cr, const char *ctg, int32_t st, int32_t en, int64_t **b_, int64_t *m_b_); 四、代码示例 #include <stdio.h> #include <stdlib.h> #include "cgranges.h" int main(void) { 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); return 0; }