快速查询区间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;
}