对代码进行了重构，基本完成了fmt-index数据结构的建立和序列查找代码，保留着一些调试信息
This commit is contained in:
Gitea
parent
6161808f60
commit
8e1576a2cd
|
|
@ -1,3 +1,7 @@
|
||||||
|
# specific for bwa_perf
|
||||||
|
*.txt
|
||||||
|
fmtidx
|
||||||
|
|
||||||
# ---> C
|
# ---> C
|
||||||
# Prerequisites
|
# Prerequisites
|
||||||
*.d
|
*.d
|
||||||
|
|
|
||||||
5
Makefile
5
Makefile
|
|
@ -1,10 +1,11 @@
|
||||||
CC= g++
|
CC= g++
|
||||||
CFLAGS= -g -Wall -Wno-unused-function -O2
|
NOWARN= -Wno-unused-result -Wno-unused-function
|
||||||
|
CFLAGS= -g -Wall $(NOWARN) -O2
|
||||||
WRAP_MALLOC=-DUSE_MALLOC_WRAPPERS
|
WRAP_MALLOC=-DUSE_MALLOC_WRAPPERS
|
||||||
SHOW_PERF= -DSHOW_PERF
|
SHOW_PERF= -DSHOW_PERF
|
||||||
AR= ar
|
AR= ar
|
||||||
DFLAGS= -DHAVE_PTHREAD $(WRAP_MALLOC) $(SHOW_PERF)
|
DFLAGS= -DHAVE_PTHREAD $(WRAP_MALLOC) $(SHOW_PERF)
|
||||||
AOBJS= sa.o fmt_index.o
|
AOBJS= util.o sa.o fmt_index.o bwt.o
|
||||||
PROG= fmtidx
|
PROG= fmtidx
|
||||||
INCLUDES=
|
INCLUDES=
|
||||||
LIBS= -lm -lz -lpthread
|
LIBS= -lm -lz -lpthread
|
||||||
|
|
|
||||||
|
|
@ -0,0 +1,310 @@
|
||||||
|
#include <stdint.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <string.h>
|
||||||
|
#include <string>
|
||||||
|
#include <iostream>
|
||||||
|
|
||||||
|
#include "util.h"
|
||||||
|
#include "bwt.h"
|
||||||
|
|
||||||
|
using namespace std;
|
||||||
|
|
||||||
|
// 计算一个字节构成的T,G,C,A序列，对应的每个碱基的个数(按T,G,C,A顺序存储在32位整数中，每个占8位)
|
||||||
|
void bwt_gen_cnt_table(bwt_t *bwt)
|
||||||
|
{
|
||||||
|
int i, j;
|
||||||
|
for (i = 0; i != 256; ++i)
|
||||||
|
{
|
||||||
|
uint32_t x = 0;
|
||||||
|
for (j = 0; j != 4; ++j)
|
||||||
|
x |= (((i & 3) == j) + ((i >> 2 & 3) == j) + ((i >> 4 & 3) == j) + (i >> 6 == j)) << (j << 3);
|
||||||
|
bwt->cnt_table[i] = x;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// 解析两bit的bwt碱基序列
|
||||||
|
bwt_t *restore_bwt_str(const char *fn)
|
||||||
|
{
|
||||||
|
bwt_t *bwt;
|
||||||
|
bwt = (bwt_t *)calloc(1, sizeof(bwt_t));
|
||||||
|
FILE *fp = fopen(fn, "rb");
|
||||||
|
|
||||||
|
fseek(fp, 0, SEEK_END);
|
||||||
|
bwt->bwt_size = (ftell(fp) - sizeof(bwtint_t) * 5) >> 2; // 以32位word为单位计算的size
|
||||||
|
bwt->bwt = (uint32_t *)calloc(bwt->bwt_size, 4);
|
||||||
|
fseek(fp, 0, SEEK_SET);
|
||||||
|
fread(&bwt->primary, sizeof(bwtint_t), 1, fp);
|
||||||
|
fread(bwt->L2 + 1, sizeof(bwtint_t), 4, fp);
|
||||||
|
fread_fix(fp, bwt->bwt_size << 2, bwt->bwt);
|
||||||
|
bwt->seq_len = bwt->L2[4];
|
||||||
|
|
||||||
|
// char *buf = (char *)calloc(bwt->seq_len + 1, 1);
|
||||||
|
// for (bwtint_t i = 0; i < bwt->seq_len; ++i)
|
||||||
|
// {
|
||||||
|
// buf[i] = BASE[bwt->bwt[i >> 4] >> ((15 - (i & 15)) << 1) & 3];
|
||||||
|
// cout << buf[i];
|
||||||
|
// }
|
||||||
|
// cout << endl;
|
||||||
|
|
||||||
|
fclose(fp);
|
||||||
|
bwt_gen_cnt_table(bwt); // 字节所能表示的各种碱基组合中，各个碱基的累积数量
|
||||||
|
return bwt;
|
||||||
|
}
|
||||||
|
|
||||||
|
// 根据原始的字符串bwt创建interval-bwt
|
||||||
|
void create_interval_occ_bwt(bwt_t *bwt)
|
||||||
|
{
|
||||||
|
bwtint_t i, k, c[4], n_occ;
|
||||||
|
uint32_t *buf;
|
||||||
|
|
||||||
|
n_occ = (bwt->seq_len + OCC_INTERVAL - 1) / OCC_INTERVAL + 1;
|
||||||
|
bwt->bwt_size += n_occ * sizeof(bwtint_t); // the new size
|
||||||
|
buf = (uint32_t *)calloc(bwt->bwt_size, 4); // will be the new bwt
|
||||||
|
c[0] = c[1] = c[2] = c[3] = 0;
|
||||||
|
// 计算occ，生成naive bwt
|
||||||
|
for (i = k = 0; i < bwt->seq_len; ++i)
|
||||||
|
{
|
||||||
|
// cout << i << '\t';
|
||||||
|
// cout << c[0] << ' ' << c[1] << ' ' << c[2] << ' ' << c[3] << endl;
|
||||||
|
if (i % OCC_INTERVAL == 0)
|
||||||
|
{
|
||||||
|
memcpy(buf + k, c, sizeof(bwtint_t) * 4);
|
||||||
|
k += sizeof(bwtint_t); // in fact: sizeof(bwtint_t)=4*(sizeof(bwtint_t)/4) 每个c包含多少个32位
|
||||||
|
// cout << "i: " << i << "\tc: " << c[0] << '\t' << c[1] << '\t' << c[2] << '\t' << c[3] << endl;
|
||||||
|
}
|
||||||
|
if (i % 16 == 0)
|
||||||
|
buf[k++] = bwt->bwt[i / 16]; // 16 == sizeof(uint32_t)/2, 2个bit表示一个碱基
|
||||||
|
++c[bwt_B00(bwt, i)];
|
||||||
|
}
|
||||||
|
// the last element
|
||||||
|
// cout << c[0] << '\t' << c[1] << '\t' << c[2] << '\t' << c[3] << endl;
|
||||||
|
memcpy(buf + k, c, sizeof(bwtint_t) * 4);
|
||||||
|
xassert(k + sizeof(bwtint_t) == bwt->bwt_size, "inconsistent bwt_size");
|
||||||
|
// update bwt
|
||||||
|
free(bwt->bwt);
|
||||||
|
bwt->bwt = buf;
|
||||||
|
}
|
||||||
|
|
||||||
|
// 对64位整型数据y，计算碱基c的累积个数
|
||||||
|
static inline int __occ_aux(uint64_t y, int c)
|
||||||
|
{
|
||||||
|
// reduce nucleotide counting to bits counting
|
||||||
|
y = ((c & 2) ? y : ~y) >> 1 & ((c & 1) ? y : ~y) & 0x5555555555555555ull;
|
||||||
|
// count the number of 1s in y
|
||||||
|
y = (y & 0x3333333333333333ull) + (y >> 2 & 0x3333333333333333ull);
|
||||||
|
return ((y + (y >> 4)) & 0xf0f0f0f0f0f0f0full) * 0x101010101010101ull >> 56;
|
||||||
|
}
|
||||||
|
|
||||||
|
// k行(包含)之前，碱基c的累计总数, interval大于等于32才能正确计算
|
||||||
|
bwtint_t bwt_occ(const bwt_t *bwt, bwtint_t k, uint8_t c)
|
||||||
|
{
|
||||||
|
bwtint_t n;
|
||||||
|
uint32_t *p, *end;
|
||||||
|
|
||||||
|
if (k == bwt->seq_len)
|
||||||
|
return bwt->L2[c + 1] - bwt->L2[c];
|
||||||
|
if (k == (bwtint_t)(-1))
|
||||||
|
return 0;
|
||||||
|
k -= (k >= bwt->primary); // because $ is not in bwt
|
||||||
|
|
||||||
|
// retrieve Occ at k/OCC_INTERVAL
|
||||||
|
n = ((bwtint_t *)(p = bwt_occ_intv(bwt, k)))[c];
|
||||||
|
// cout << "bwt_occ - 1: " << (int)c << '\t' << k << '\t' << n << endl;
|
||||||
|
p += sizeof(bwtint_t); // jump to the start of the first BWT cell
|
||||||
|
|
||||||
|
// calculate Occ up to the last k/32
|
||||||
|
end = p + (((k >> 5) - ((k & ~OCC_INTV_MASK) >> 5)) << 1);
|
||||||
|
for (; p < end; p += 2)
|
||||||
|
n += __occ_aux((uint64_t)p[0] << 32 | p[1], c);
|
||||||
|
|
||||||
|
// cout << "bwt_occ - 2: " << (int)c << '\t' << k << '\t' << n << endl;
|
||||||
|
// calculate Occ
|
||||||
|
n += __occ_aux(((uint64_t)p[0] << 32 | p[1]) & ~((1ull << ((~k & 31) << 1)) - 1), c);
|
||||||
|
if (c == 0)
|
||||||
|
n -= ~k & 31; // corrected for the masked bits
|
||||||
|
// cout << "bwt_occ - 3: " << (int)c << '\t' << k << '\t' << n << endl;
|
||||||
|
return n;
|
||||||
|
}
|
||||||
|
|
||||||
|
// 统计k行（bwt mtx行）之前4种碱基累积数量，这里的k是bwt矩阵里的行，比bwt字符串多1
|
||||||
|
void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4])
|
||||||
|
{
|
||||||
|
bwtint_t x;
|
||||||
|
uint32_t *p, tmp, *end;
|
||||||
|
if (k == (bwtint_t)(-1))
|
||||||
|
{
|
||||||
|
memset(cnt, 0, 4 * sizeof(bwtint_t));
|
||||||
|
return;
|
||||||
|
}
|
||||||
|
k -= (k >= bwt->primary); // because $ is not in bwt
|
||||||
|
p = bwt_occ_intv(bwt, k);
|
||||||
|
memcpy(cnt, p, 4 * sizeof(bwtint_t));
|
||||||
|
p += sizeof(bwtint_t); // sizeof(bwtint_t) = 4*(sizeof(bwtint_t)/sizeof(uint32_t))
|
||||||
|
end = p + ((k >> 4) - ((k & ~OCC_INTV_MASK) >> 4)); // this is the end point of the following loop
|
||||||
|
for (x = 0; p < end; ++p)
|
||||||
|
x += __occ_aux4(bwt, *p);
|
||||||
|
tmp = *p & ~((1U << ((~k & 15) << 1)) - 1);
|
||||||
|
x += __occ_aux4(bwt, tmp) - (~k & 15); // 这里多算了A，要减去
|
||||||
|
cnt[0] += x & 0xff;
|
||||||
|
cnt[1] += x >> 8 & 0xff;
|
||||||
|
cnt[2] += x >> 16 & 0xff;
|
||||||
|
cnt[3] += x >> 24;
|
||||||
|
}
|
||||||
|
|
||||||
|
// an analogy to bwt_occ4() but more efficient, requiring k <= l
|
||||||
|
void bwt_2occ4(const bwt_t *bwt, bwtint_t k, bwtint_t l, bwtint_t cntk[4], bwtint_t cntl[4])
|
||||||
|
{
|
||||||
|
bwtint_t _k, _l;
|
||||||
|
_k = k - (k >= bwt->primary);
|
||||||
|
_l = l - (l >= bwt->primary);
|
||||||
|
if (_l >> OCC_INTV_SHIFT != _k >> OCC_INTV_SHIFT || k == (bwtint_t)(-1) || l == (bwtint_t)(-1))
|
||||||
|
{
|
||||||
|
bwt_occ4(bwt, k, cntk);
|
||||||
|
bwt_occ4(bwt, l, cntl);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
bwtint_t x, y;
|
||||||
|
uint32_t *p, tmp, *endk, *endl;
|
||||||
|
k -= (k >= bwt->primary); // because $ is not in bwt
|
||||||
|
l -= (l >= bwt->primary);
|
||||||
|
p = bwt_occ_intv(bwt, k);
|
||||||
|
memcpy(cntk, p, 4 * sizeof(bwtint_t));
|
||||||
|
p += sizeof(bwtint_t); // sizeof(bwtint_t) = 4*(sizeof(bwtint_t)/sizeof(uint32_t))
|
||||||
|
// prepare cntk[]
|
||||||
|
endk = p + ((k >> 4) - ((k & ~OCC_INTV_MASK) >> 4));
|
||||||
|
endl = p + ((l >> 4) - ((l & ~OCC_INTV_MASK) >> 4));
|
||||||
|
for (x = 0; p < endk; ++p)
|
||||||
|
x += __occ_aux4(bwt, *p);
|
||||||
|
y = x;
|
||||||
|
tmp = *p & ~((1U << ((~k & 15) << 1)) - 1);
|
||||||
|
x += __occ_aux4(bwt, tmp) - (~k & 15);
|
||||||
|
// calculate cntl[] and finalize cntk[]
|
||||||
|
for (; p < endl; ++p)
|
||||||
|
y += __occ_aux4(bwt, *p);
|
||||||
|
tmp = *p & ~((1U << ((~l & 15) << 1)) - 1);
|
||||||
|
y += __occ_aux4(bwt, tmp) - (~l & 15);
|
||||||
|
memcpy(cntl, cntk, 4 * sizeof(bwtint_t));
|
||||||
|
cntk[0] += x & 0xff;
|
||||||
|
cntk[1] += x >> 8 & 0xff;
|
||||||
|
cntk[2] += x >> 16 & 0xff;
|
||||||
|
cntk[3] += x >> 24;
|
||||||
|
cntl[0] += y & 0xff;
|
||||||
|
cntl[1] += y >> 8 & 0xff;
|
||||||
|
cntl[2] += y >> 16 & 0xff;
|
||||||
|
cntl[3] += y >> 24;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void bwt_extend(const bwt_t *bwt, const bwtintv_t *ik, bwtintv_t ok[4], int is_back)
|
||||||
|
{
|
||||||
|
bwtint_t tk[4], tl[4];
|
||||||
|
int i;
|
||||||
|
bwt_2occ4(bwt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl); // tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积
|
||||||
|
// 这里是反向扩展
|
||||||
|
for (i = 0; i != 4; ++i)
|
||||||
|
{
|
||||||
|
ok[i].x[!is_back] = bwt->L2[i] + 1 + tk[i]; // 起始行位置，互补链
|
||||||
|
ok[i].x[2] = tl[i] - tk[i]; // 间隔
|
||||||
|
}
|
||||||
|
// 因为计算的是互补碱基，所以3对应着0,2对应1，下边是正向扩展
|
||||||
|
ok[3].x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= bwt->primary && ik->x[!is_back] + ik->x[2] - 1 >= bwt->primary);
|
||||||
|
ok[2].x[is_back] = ok[3].x[is_back] + ok[3].x[2];
|
||||||
|
ok[1].x[is_back] = ok[2].x[is_back] + ok[2].x[2];
|
||||||
|
ok[0].x[is_back] = ok[1].x[is_back] + ok[1].x[2];
|
||||||
|
}
|
||||||
|
|
||||||
|
// 利用bwt搜索seed，完整搜索，只需要单向搜索
|
||||||
|
void bwt_search(bwt_t *bwt, const string &q)
|
||||||
|
{
|
||||||
|
bwtintv_t ik, ok[4];
|
||||||
|
int i, j, c, x = 0;
|
||||||
|
|
||||||
|
bwt_set_intv(bwt, bval(q[x]), ik);
|
||||||
|
ik.info = x + 1;
|
||||||
|
cout << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||||
|
|
||||||
|
for (i = x + 1; i < q.size(); ++i)
|
||||||
|
{
|
||||||
|
if (bval(q[i]) < 4)
|
||||||
|
{
|
||||||
|
c = cbval(q[i]);
|
||||||
|
bwt_extend(bwt, &ik, ok, 0);
|
||||||
|
ik = ok[c];
|
||||||
|
ik.info = i + 1;
|
||||||
|
cout << "bwt-1: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// 扩展两次
|
||||||
|
void bwt_extend2(const bwt_t *bwt, bwtintv_t *ik, bwtintv_t ok[4], int is_back, int c1)
|
||||||
|
{
|
||||||
|
bwtint_t tk[4], tl[4];
|
||||||
|
int i;
|
||||||
|
bwt_2occ4(bwt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl); // tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积
|
||||||
|
// 这里是反向扩展
|
||||||
|
for (i = 0; i != 4; ++i)
|
||||||
|
{
|
||||||
|
ok[i].x[!is_back] = bwt->L2[i] + 1 + tk[i]; // 起始行位置，互补链
|
||||||
|
ok[i].x[2] = tl[i] - tk[i]; // 间隔
|
||||||
|
}
|
||||||
|
// 因为计算的是互补碱基，所以3对应着0,2对应1，下边是正向扩展
|
||||||
|
ok[3].x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= bwt->primary && ik->x[!is_back] + ik->x[2] - 1 >= bwt->primary);
|
||||||
|
ok[2].x[is_back] = ok[3].x[is_back] + ok[3].x[2];
|
||||||
|
ok[1].x[is_back] = ok[2].x[is_back] + ok[2].x[2];
|
||||||
|
ok[0].x[is_back] = ok[1].x[is_back] + ok[1].x[2];
|
||||||
|
|
||||||
|
*ik = ok[c1];
|
||||||
|
bwt_2occ4(bwt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl);
|
||||||
|
for (i = 0; i != 4; ++i)
|
||||||
|
{
|
||||||
|
ok[i].x[!is_back] = bwt->L2[i] + 1 + tk[i]; // 起始行位置，互补链
|
||||||
|
ok[i].x[2] = tl[i] - tk[i]; // 间隔
|
||||||
|
}
|
||||||
|
// 因为计算的是互补碱基，所以3对应着0,2对应1，下边是正向扩展
|
||||||
|
ok[3].x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= bwt->primary && ik->x[!is_back] + ik->x[2] - 1 >= bwt->primary);
|
||||||
|
ok[2].x[is_back] = ok[3].x[is_back] + ok[3].x[2];
|
||||||
|
ok[1].x[is_back] = ok[2].x[is_back] + ok[2].x[2];
|
||||||
|
ok[0].x[is_back] = ok[1].x[is_back] + ok[1].x[2];
|
||||||
|
}
|
||||||
|
|
||||||
|
// 每次扩展两步
|
||||||
|
void bwt_search2(bwt_t *bwt, const string &q)
|
||||||
|
{
|
||||||
|
bwtintv_t ik, ok[4];
|
||||||
|
int i, j, c1, c2, x = 0;
|
||||||
|
|
||||||
|
bwt_set_intv(bwt, bval(q[x]), ik);
|
||||||
|
ik.info = x + 1;
|
||||||
|
cout << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||||
|
|
||||||
|
for (i = x + 1; i + 1 < q.size(); i += 2)
|
||||||
|
{
|
||||||
|
if (bval(q[i]) < 4 && bval(q[i + 1]) < 4)
|
||||||
|
{
|
||||||
|
|
||||||
|
c1 = cbval(q[i]);
|
||||||
|
c2 = cbval(q[i + 1]);
|
||||||
|
|
||||||
|
bwt_extend2(bwt, &ik, ok, 0, c1);
|
||||||
|
|
||||||
|
ik = ok[c2];
|
||||||
|
ik.info = i + 1;
|
||||||
|
cout << "bwt-2: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (i < q.size() && bval(q[i]) < 4)
|
||||||
|
{ // 最后一次扩展
|
||||||
|
c1 = cbval(q[i]);
|
||||||
|
bwt_extend(bwt, &ik, ok, 0);
|
||||||
|
ik = ok[c1];
|
||||||
|
ik.info = i + 1;
|
||||||
|
cout << "bwt-2: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
@ -0,0 +1,68 @@
|
||||||
|
#ifndef BWT_H_
|
||||||
|
#define BWT_H_
|
||||||
|
#include <stdint.h>
|
||||||
|
#include <string>
|
||||||
|
#include "util.h"
|
||||||
|
|
||||||
|
using std::string;
|
||||||
|
|
||||||
|
// occ间隔对应的右移位数，5表示间隔32行保存一次
|
||||||
|
#define OCC_INTV_SHIFT 5
|
||||||
|
#define OCC_INTERVAL (1LL << OCC_INTV_SHIFT)
|
||||||
|
#define OCC_INTV_MASK (OCC_INTERVAL - 1)
|
||||||
|
|
||||||
|
// 从最原始的bwt碱基串里获取k行对应的碱基，这里的bwt不包括occ check point数据
|
||||||
|
#define bwt_B00(b, k) ((b)->bwt[(k) >> 4] >> ((~(k) & 0xf) << 1) & 3)
|
||||||
|
|
||||||
|
/* retrieve a character from the $-removed BWT string. Note that
|
||||||
|
* bwt_t::bwt is not exactly the BWT string and therefore this macro is
|
||||||
|
* called bwt_B0 instead of bwt_B */
|
||||||
|
// 从构建完成的bwt（包含occ check point）获取k行（不含$，这里的k不输出bwt mtx的行，是bwt字符串的行）的碱基
|
||||||
|
#define bwt_B0(b, k) (bwt_bwt(b, k) >> ((~(k) & 0xf) << 1) & 3)
|
||||||
|
|
||||||
|
// 获取碱基c（待查找序列的首）和对应的互补碱基对应的行，以及间隔
|
||||||
|
#define bwt_set_intv(bwt, c, ik) ((ik).x[0] = (bwt)->L2[(int)(c)] + 1, (ik).x[2] = (bwt)->L2[(int)(c) + 1] - (bwt)->L2[(int)(c)], (ik).x[1] = (bwt)->L2[3 - (c)] + 1, (ik).info = 0)
|
||||||
|
|
||||||
|
// The following two lines are ONLY correct when OCC_INTERVAL==0x80
|
||||||
|
// #define bwt_bwt(b, k) ((b)->bwt[((k) >> 7 << 4) + sizeof(bwtint_t) + (((k) & 0x7f) >> 4)])
|
||||||
|
// #define bwt_occ_intv(b, k) ((b)->bwt + ((k) >> 7 << 4))
|
||||||
|
///* For general OCC_INTERVAL, the following is correct:
|
||||||
|
// k行对应的bwt str碱基，对应的check point bwt str数据起始地址
|
||||||
|
#define bwt_bwt(b, k) ((b)->bwt[(k) / OCC_INTERVAL * (OCC_INTERVAL / (sizeof(uint32_t) * 8 / 2) + sizeof(bwtint_t) / 4 * 4) + sizeof(bwtint_t) / 4 * 4 + (k) % OCC_INTERVAL / 16])
|
||||||
|
// k行（bwt str行（不包含$））对应的check point occ数据起始地址（小于k且是OCC_INTERVAL的整数倍）
|
||||||
|
#define bwt_occ_intv(b, k) ((b)->bwt + (k) / OCC_INTERVAL * (OCC_INTERVAL / (sizeof(uint32_t) * 8 / 2) + sizeof(bwtint_t) / 4 * 4))
|
||||||
|
// 字节b中包含的T G C A（按顺序保存在32位整数里（每个占8bit））的数量
|
||||||
|
#define __occ_aux4(bwt, b) \
|
||||||
|
((bwt)->cnt_table[(b) & 0xff] + (bwt)->cnt_table[(b) >> 8 & 0xff] + (bwt)->cnt_table[(b) >> 16 & 0xff] + (bwt)->cnt_table[(b) >> 24])
|
||||||
|
|
||||||
|
// 原始fm-index (bwt)结构
|
||||||
|
struct bwt_t
|
||||||
|
{
|
||||||
|
bwtint_t primary; // S^{-1}(0), or the primary index of BWT
|
||||||
|
bwtint_t L2[5]; // C(), cumulative count
|
||||||
|
bwtint_t seq_len; // sequence length
|
||||||
|
bwtint_t bwt_size; // size of bwt, about seq_len/4
|
||||||
|
uint32_t *bwt; // BWT
|
||||||
|
// occurance array, separated to two parts
|
||||||
|
uint32_t cnt_table[256];
|
||||||
|
// suffix array
|
||||||
|
int sa_intv;
|
||||||
|
bwtint_t n_sa;
|
||||||
|
uint8_t *sa;
|
||||||
|
};
|
||||||
|
|
||||||
|
// k行(包含，bwt mtx行)之前，碱基c的累计总数, interval大于等于32才能正确计算
|
||||||
|
bwtint_t bwt_occ(const bwt_t *bwt, bwtint_t k, uint8_t c);
|
||||||
|
// 统计k行（bwt mtx行，包含k行本身）之前4种碱基累积数量，这里的k是bwt矩阵里的行，比bwt字符串多1
|
||||||
|
void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4]);
|
||||||
|
// 解析两bit的bwt碱基序列
|
||||||
|
bwt_t *restore_bwt_str(const char *fn);
|
||||||
|
// 根据原始的字符串bwt创建interval-bwt
|
||||||
|
void create_interval_occ_bwt(bwt_t *bwt);
|
||||||
|
|
||||||
|
// 利用bwt搜索seed，完整搜索，只需要单向搜索
|
||||||
|
void bwt_search(bwt_t *bwt, const string &q);
|
||||||
|
// 每次扩展两步
|
||||||
|
void bwt_search2(bwt_t *bwt, const string &q);
|
||||||
|
|
||||||
|
#endif
|
||||||
34
common.h
34
common.h
|
|
@ -1,34 +0,0 @@
|
||||||
#ifndef COMMON_H_
|
|
||||||
#define COMMON_H_
|
|
||||||
|
|
||||||
#define OCC_INTV_SHIFT 5
|
|
||||||
#define OCC_INTERVAL (1LL << OCC_INTV_SHIFT)
|
|
||||||
#define OCC_INTV_MASK (OCC_INTERVAL - 1)
|
|
||||||
#define bwt_B00(b, k) ((b)->bwt[(k) >> 4] >> ((~(k) & 0xf) << 1) & 3)
|
|
||||||
|
|
||||||
/* retrieve a character from the $-removed BWT string. Note that
|
|
||||||
* bwt_t::bwt is not exactly the BWT string and therefore this macro is
|
|
||||||
* called bwt_B0 instead of bwt_B */
|
|
||||||
#define bwt_B0(b, k) (bwt_bwt(b, k) >> ((~(k) & 0xf) << 1) & 3)
|
|
||||||
|
|
||||||
#define bwt_set_intv(bwt, c, ik) ((ik).x[0] = (bwt)->L2[(int)(c)] + 1, (ik).x[2] = (bwt)->L2[(int)(c) + 1] - (bwt)->L2[(int)(c)], (ik).x[1] = (bwt)->L2[3 - (c)] + 1, (ik).info = 0)
|
|
||||||
|
|
||||||
|
|
||||||
#define SA_BYTES_33(n_sa) ((((33 * (n_sa) + 7) / 8) & (~7L)) + 8)
|
|
||||||
#define SA_BYTES_40(n_sa) ((((40 * (n_sa) + 7) / 8) & (~7L)) + 8)
|
|
||||||
|
|
||||||
#define xassert(cond, msg) \
|
|
||||||
if ((cond) == 0) \
|
|
||||||
_err_fatal_simple_core(__func__, msg)
|
|
||||||
|
|
||||||
typedef uint64_t bwtint_t;
|
|
||||||
|
|
||||||
double realtime(void);
|
|
||||||
|
|
||||||
void bwt_set_sa_33(uint8_t *sa_arr, bwtint_t k, bwtint_t val);
|
|
||||||
bwtint_t bwt_get_sa_33(uint8_t *sa_arr, bwtint_t k);
|
|
||||||
|
|
||||||
int main_sa(int argc, char **argv);
|
|
||||||
int main_fmtidx(int argc, char **argv);
|
|
||||||
|
|
||||||
#endif
|
|
||||||
412
fmt_index.cpp
412
fmt_index.cpp
|
|
@ -8,112 +8,12 @@
|
||||||
#include <algorithm>
|
#include <algorithm>
|
||||||
#include <string.h>
|
#include <string.h>
|
||||||
|
|
||||||
#include "common.h"
|
#include "util.h"
|
||||||
|
#include "fmt_index.h"
|
||||||
using namespace std;
|
using namespace std;
|
||||||
|
|
||||||
///* For general OCC_INTERVAL, the following is correct:
|
|
||||||
#define bwt_bwt(b, k) ((b)->bwt[(k) / OCC_INTERVAL * (OCC_INTERVAL / (sizeof(uint32_t) * 8 / 2) + sizeof(bwtint_t) / 4 * 4) + sizeof(bwtint_t) / 4 * 4 + (k) % OCC_INTERVAL / 16])
|
|
||||||
#define bwt_occ_intv(b, k) ((b)->bwt + (k) / OCC_INTERVAL * (OCC_INTERVAL / (sizeof(uint32_t) * 8 / 2) + sizeof(bwtint_t) / 4 * 4))
|
|
||||||
//*/
|
|
||||||
|
|
||||||
// The following two lines are ONLY correct when OCC_INTERVAL==0x80
|
|
||||||
// #define bwt_bwt(b, k) ((b)->bwt[((k) >> 7 << 4) + sizeof(bwtint_t) + (((k) & 0x7f) >> 4)])
|
|
||||||
// #define bwt_occ_intv(b, k) ((b)->bwt + ((k) >> 7 << 4))
|
|
||||||
|
|
||||||
/* retrieve a character from the $-removed BWT string. Note that
|
|
||||||
* bwt_t::bwt is not exactly the BWT string and therefore this macro is
|
|
||||||
* called bwt_B0 instead of bwt_B */
|
|
||||||
#define bwt_B0(b, k) (bwt_bwt(b, k) >> ((~(k) & 0xf) << 1) & 3)
|
|
||||||
|
|
||||||
#define bwt_set_intv(bwt, c, ik) ((ik).x[0] = (bwt)->L2[(int)(c)] + 1, (ik).x[2] = (bwt)->L2[(int)(c) + 1] - (bwt)->L2[(int)(c)], (ik).x[1] = (bwt)->L2[3 - (c)] + 1, (ik).info = 0)
|
|
||||||
|
|
||||||
#define __occ_aux4(bwt, b) \
|
|
||||||
((bwt)->cnt_table[(b) & 0xff] + (bwt)->cnt_table[(b) >> 8 & 0xff] + (bwt)->cnt_table[(b) >> 16 & 0xff] + (bwt)->cnt_table[(b) >> 24])
|
|
||||||
|
|
||||||
const char BASE[4] = {'A', 'C', 'G', 'T'};
|
const char BASE[4] = {'A', 'C', 'G', 'T'};
|
||||||
|
|
||||||
// base转成2bit值
|
|
||||||
inline int bval(char b)
|
|
||||||
{
|
|
||||||
if (b == 'A')
|
|
||||||
return 0;
|
|
||||||
if (b == 'C')
|
|
||||||
return 1;
|
|
||||||
if (b == 'G')
|
|
||||||
return 2;
|
|
||||||
if (b == 'T')
|
|
||||||
return 3;
|
|
||||||
return 4;
|
|
||||||
}
|
|
||||||
|
|
||||||
// 互补碱基值
|
|
||||||
inline int cbval(char b)
|
|
||||||
{
|
|
||||||
return 3 - bval(b);
|
|
||||||
}
|
|
||||||
|
|
||||||
struct bwtintv_t
|
|
||||||
{
|
|
||||||
bwtint_t x[3], info; // x[0]表示正链位置，x[1]表示互补链位置，x[2]表示间隔长度，info 表示read的起始结束位置
|
|
||||||
};
|
|
||||||
|
|
||||||
// 原始fm-index结构
|
|
||||||
struct bwt_t
|
|
||||||
{
|
|
||||||
bwtint_t primary; // S^{-1}(0), or the primary index of BWT
|
|
||||||
bwtint_t L2[5]; // C(), cumulative count
|
|
||||||
bwtint_t seq_len; // sequence length
|
|
||||||
bwtint_t bwt_size; // size of bwt, about seq_len/4
|
|
||||||
uint32_t *bwt; // BWT
|
|
||||||
// occurance array, separated to two parts
|
|
||||||
uint32_t cnt_table[256];
|
|
||||||
// suffix array
|
|
||||||
int sa_intv;
|
|
||||||
bwtint_t n_sa;
|
|
||||||
uint8_t *sa;
|
|
||||||
};
|
|
||||||
|
|
||||||
// fm-index, twice calc in one memory access
|
|
||||||
struct FMTIndex {
|
|
||||||
bwtint_t primary; // S^{-1}(0), or the primary index of BWT
|
|
||||||
bwtint_t sec_primary; // second primary line
|
|
||||||
bwtint_t L2[5]; // C(), cumulative count
|
|
||||||
bwtint_t seq_len; // sequence length
|
|
||||||
bwtint_t bwt_size; // size of bwt, about seq_len/4
|
|
||||||
uint32_t *bwt; // BWT
|
|
||||||
// occurance array, separated to two parts
|
|
||||||
uint32_t cnt_table[5][256]; // 4对应原来的cnt_table，0,1,2,3,分别对应该碱基的扩展
|
|
||||||
int sec_bcp; // base couple for sec primary line, AA=>0, AC=>1 ... TT=>15
|
|
||||||
int first_base; // 序列的第一个碱基2bit的int类型，0,1,2,3
|
|
||||||
int last_base; // dollar转换成的base
|
|
||||||
// suffix array
|
|
||||||
int sa_intv;
|
|
||||||
bwtint_t n_sa;
|
|
||||||
uint8_t *sa;
|
|
||||||
};
|
|
||||||
|
|
||||||
void _err_fatal_simple_core(const char *func, const char *msg)
|
|
||||||
{
|
|
||||||
fprintf(stderr, "[%s] %s Abort!\n", func, msg);
|
|
||||||
abort();
|
|
||||||
}
|
|
||||||
|
|
||||||
// 读取最原始的bwt，
|
|
||||||
static bwtint_t fread_fix(FILE *fp, bwtint_t size, void *a)
|
|
||||||
{ // Mac/Darwin has a bug when reading data longer than 2GB. This function fixes this issue by reading data in small chunks
|
|
||||||
const int bufsize = 0x1000000; // 16M block
|
|
||||||
bwtint_t offset = 0;
|
|
||||||
while (size)
|
|
||||||
{
|
|
||||||
int x = bufsize < size ? bufsize : size;
|
|
||||||
if ((x = fread((uint8_t *)a + offset, 1, x, fp)) == 0)
|
|
||||||
break;
|
|
||||||
size -= x;
|
|
||||||
offset += x;
|
|
||||||
}
|
|
||||||
return offset;
|
|
||||||
}
|
|
||||||
|
|
||||||
// 求反向互补序列
|
// 求反向互补序列
|
||||||
string calc_reverse_seq(string &seq)
|
string calc_reverse_seq(string &seq)
|
||||||
{
|
{
|
||||||
|
|
@ -133,69 +33,15 @@ string calc_reverse_seq(string &seq)
|
||||||
return rseq;
|
return rseq;
|
||||||
}
|
}
|
||||||
|
|
||||||
static inline int __occ_aux(uint64_t y, int c)
|
// 打印32位整型数据中包含的pre-bwt：bwt
|
||||||
|
void print_base_uint32(uint32_t p)
|
||||||
{
|
{
|
||||||
// reduce nucleotide counting to bits counting
|
for (int i = 30; i > 0; i -= 4)
|
||||||
y = ((c & 2) ? y : ~y) >> 1 & ((c & 1) ? y : ~y) & 0x5555555555555555ull;
|
|
||||||
// count the number of 1s in y
|
|
||||||
y = (y & 0x3333333333333333ull) + (y >> 2 & 0x3333333333333333ull);
|
|
||||||
return ((y + (y >> 4)) & 0xf0f0f0f0f0f0f0full) * 0x101010101010101ull >> 56;
|
|
||||||
}
|
|
||||||
|
|
||||||
// k行(包含)之前，碱基c的累计总数, interval大于等于32才能正确计算
|
|
||||||
bwtint_t bwt_occ(const bwt_t *bwt, bwtint_t k, uint8_t c)
|
|
||||||
{
|
|
||||||
bwtint_t n;
|
|
||||||
uint32_t *p, *end;
|
|
||||||
|
|
||||||
if (k == bwt->seq_len)
|
|
||||||
return bwt->L2[c + 1] - bwt->L2[c];
|
|
||||||
if (k == (bwtint_t)(-1))
|
|
||||||
return 0;
|
|
||||||
k -= (k >= bwt->primary); // because $ is not in bwt
|
|
||||||
|
|
||||||
// retrieve Occ at k/OCC_INTERVAL
|
|
||||||
n = ((bwtint_t *)(p = bwt_occ_intv(bwt, k)))[c];
|
|
||||||
// cout << "bwt_occ - 1: " << (int)c << '\t' << k << '\t' << n << endl;
|
|
||||||
p += sizeof(bwtint_t); // jump to the start of the first BWT cell
|
|
||||||
|
|
||||||
// calculate Occ up to the last k/32
|
|
||||||
end = p + (((k >> 5) - ((k & ~OCC_INTV_MASK) >> 5)) << 1);
|
|
||||||
for (; p < end; p += 2)
|
|
||||||
n += __occ_aux((uint64_t)p[0] << 32 | p[1], c);
|
|
||||||
|
|
||||||
//cout << "bwt_occ - 2: " << (int)c << '\t' << k << '\t' << n << endl;
|
|
||||||
// calculate Occ
|
|
||||||
n += __occ_aux(((uint64_t)p[0] << 32 | p[1]) & ~((1ull << ((~k & 31) << 1)) - 1), c);
|
|
||||||
if (c == 0)
|
|
||||||
n -= ~k & 31; // corrected for the masked bits
|
|
||||||
//cout << "bwt_occ - 3: " << (int)c << '\t' << k << '\t' << n << endl;
|
|
||||||
return n;
|
|
||||||
}
|
|
||||||
|
|
||||||
// 这里的k是bwt矩阵里的行，比bwt字符串多1
|
|
||||||
void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4])
|
|
||||||
{
|
|
||||||
bwtint_t x;
|
|
||||||
uint32_t *p, tmp, *end;
|
|
||||||
if (k == (bwtint_t)(-1))
|
|
||||||
{
|
{
|
||||||
memset(cnt, 0, 4 * sizeof(bwtint_t));
|
int b1 = p >> i & 3;
|
||||||
return;
|
int b2 = p >> (i - 2) & 3;
|
||||||
|
cout << BASE[b1] << BASE[b2] << endl;
|
||||||
}
|
}
|
||||||
k -= (k >= bwt->primary); // because $ is not in bwt
|
|
||||||
p = bwt_occ_intv(bwt, k);
|
|
||||||
memcpy(cnt, p, 4 * sizeof(bwtint_t));
|
|
||||||
p += sizeof(bwtint_t); // sizeof(bwtint_t) = 4*(sizeof(bwtint_t)/sizeof(uint32_t))
|
|
||||||
end = p + ((k >> 4) - ((k & ~OCC_INTV_MASK) >> 4)); // this is the end point of the following loop
|
|
||||||
for (x = 0; p < end; ++p)
|
|
||||||
x += __occ_aux4(bwt, *p);
|
|
||||||
tmp = *p & ~((1U << ((~k & 15) << 1)) - 1);
|
|
||||||
x += __occ_aux4(bwt, tmp) - (~k & 15); // 这里多算了A，要减去
|
|
||||||
cnt[0] += x & 0xff;
|
|
||||||
cnt[1] += x >> 8 & 0xff;
|
|
||||||
cnt[2] += x >> 16 & 0xff;
|
|
||||||
cnt[3] += x >> 24;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// 创建bwt矩阵
|
// 创建bwt矩阵
|
||||||
|
|
@ -208,7 +54,6 @@ void create_bwt_mtx(string &seq)
|
||||||
{
|
{
|
||||||
sarr[i] = sarr[0].substr(i) + sarr[0].substr(0, i);
|
sarr[i] = sarr[0].substr(i) + sarr[0].substr(0, i);
|
||||||
}
|
}
|
||||||
|
|
||||||
std::sort(sarr, sarr + seq.size() + 1);
|
std::sort(sarr, sarr + seq.size() + 1);
|
||||||
|
|
||||||
// bwt matrix
|
// bwt matrix
|
||||||
|
|
@ -224,7 +69,6 @@ void create_bwt_mtx(string &seq)
|
||||||
// cout << sarr[i].back();
|
// cout << sarr[i].back();
|
||||||
// }
|
// }
|
||||||
// cout << endl;
|
// cout << endl;
|
||||||
//
|
|
||||||
// cout << "pre bwt string" << endl;
|
// cout << "pre bwt string" << endl;
|
||||||
// for (int i = 0; i < sarr[0].size(); ++i)
|
// for (int i = 0; i < sarr[0].size(); ++i)
|
||||||
// {
|
// {
|
||||||
|
|
@ -233,20 +77,7 @@ void create_bwt_mtx(string &seq)
|
||||||
// cout << endl;
|
// cout << endl;
|
||||||
}
|
}
|
||||||
|
|
||||||
// 计算一个字节构成的A,T,C,G序列，对应的每个碱基的个数，因为最多有4个相同的碱基，所以每次左移3位就行
|
// fmt-index的count table，4对应着bwt碱基的累积量，0,1,2,3分别对应着bwt是A,C,G,T，pre-bwt的累积量
|
||||||
void bwt_gen_cnt_table(bwt_t *bwt)
|
|
||||||
{
|
|
||||||
int i, j;
|
|
||||||
for (i = 0; i != 256; ++i)
|
|
||||||
{
|
|
||||||
uint32_t x = 0;
|
|
||||||
for (j = 0; j != 4; ++j)
|
|
||||||
x |= (((i & 3) == j) + ((i >> 2 & 3) == j) + ((i >> 4 & 3) == j) + (i >> 6 == j)) << (j << 3);
|
|
||||||
bwt->cnt_table[i] = x;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// fmt-index的count table
|
|
||||||
void fmt_gen_cnt_table(FMTIndex *fmt)
|
void fmt_gen_cnt_table(FMTIndex *fmt)
|
||||||
{
|
{
|
||||||
int i, j, k;
|
int i, j, k;
|
||||||
|
|
@ -267,79 +98,7 @@ void fmt_gen_cnt_table(FMTIndex *fmt)
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
void print_base_uint32(uint32_t p)
|
|
||||||
{
|
|
||||||
for (int i = 30; i > 0; i -= 4)
|
|
||||||
{
|
|
||||||
int b1 = p >> i & 3;
|
|
||||||
int b2 = p >> (i - 2) & 3;
|
|
||||||
cout << BASE[b1] << BASE[b2] << endl;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// 解析两bit碱基序列
|
|
||||||
bwt_t *restore_bwt_str(const char *fn)
|
|
||||||
{
|
|
||||||
bwt_t *bwt;
|
|
||||||
bwt = (bwt_t *)calloc(1, sizeof(bwt_t));
|
|
||||||
FILE *fp = fopen(fn, "rb");
|
|
||||||
char *buf;
|
|
||||||
|
|
||||||
fseek(fp, 0, SEEK_END);
|
|
||||||
bwt->bwt_size = (ftell(fp) - sizeof(bwtint_t) * 5) >> 2; // 以32位word为单位计算的size
|
|
||||||
bwt->bwt = (uint32_t *)calloc(bwt->bwt_size, 4);
|
|
||||||
fseek(fp, 0, SEEK_SET);
|
|
||||||
fread(&bwt->primary, sizeof(bwtint_t), 1, fp);
|
|
||||||
fread(bwt->L2 + 1, sizeof(bwtint_t), 4, fp);
|
|
||||||
fread_fix(fp, bwt->bwt_size << 2, bwt->bwt);
|
|
||||||
bwt->seq_len = bwt->L2[4];
|
|
||||||
|
|
||||||
// buf = (char *)calloc(bwt->seq_len + 1, 1);
|
|
||||||
// for (bwtint_t i = 0; i < bwt->seq_len; ++i)
|
|
||||||
// {
|
|
||||||
// buf[i] = BASE[bwt->bwt[i >> 4] >> ((15 - (i & 15)) << 1) & 3];
|
|
||||||
// cout << buf[i];
|
|
||||||
// }
|
|
||||||
// cout << endl;
|
|
||||||
|
|
||||||
fclose(fp);
|
|
||||||
bwt_gen_cnt_table(bwt); // 字节所能表示的各种碱基组合中，各个碱基的累积数量
|
|
||||||
return bwt;
|
|
||||||
}
|
|
||||||
|
|
||||||
// 根据原始的字符串bwt创建interval-bwt
|
|
||||||
void create_interval_occ_bwt(bwt_t *bwt)
|
|
||||||
{
|
|
||||||
bwtint_t i, k, c[4], n_occ;
|
|
||||||
uint32_t *buf;
|
|
||||||
|
|
||||||
n_occ = (bwt->seq_len + OCC_INTERVAL - 1) / OCC_INTERVAL + 1;
|
|
||||||
bwt->bwt_size += n_occ * sizeof(bwtint_t); // the new size
|
|
||||||
buf = (uint32_t *)calloc(bwt->bwt_size, 4); // will be the new bwt
|
|
||||||
c[0] = c[1] = c[2] = c[3] = 0;
|
|
||||||
// 计算occ，生成naive bwt
|
|
||||||
for (i = k = 0; i < bwt->seq_len; ++i)
|
|
||||||
{
|
|
||||||
// cout << i << '\t';
|
|
||||||
// cout << c[0] << ' ' << c[1] << ' ' << c[2] << ' ' << c[3] << endl;
|
|
||||||
if (i % OCC_INTERVAL == 0)
|
|
||||||
{
|
|
||||||
memcpy(buf + k, c, sizeof(bwtint_t) * 4);
|
|
||||||
k += sizeof(bwtint_t); // in fact: sizeof(bwtint_t)=4*(sizeof(bwtint_t)/4) 每个c包含多少个32位
|
|
||||||
// cout << "i: " << i << "\tc: " << c[0] << '\t' << c[1] << '\t' << c[2] << '\t' << c[3] << endl;
|
|
||||||
}
|
|
||||||
if (i % 16 == 0)
|
|
||||||
buf[k++] = bwt->bwt[i / 16]; // 16 == sizeof(uint32_t)/2, 2个bit表示一个碱基
|
|
||||||
++c[bwt_B00(bwt, i)];
|
|
||||||
}
|
|
||||||
// the last element
|
|
||||||
// cout << c[0] << '\t' << c[1] << '\t' << c[2] << '\t' << c[3] << endl;
|
|
||||||
memcpy(buf + k, c, sizeof(bwtint_t) * 4);
|
|
||||||
xassert(k + sizeof(bwtint_t) == bwt->bwt_size, "inconsistent bwt_size");
|
|
||||||
// update bwt
|
|
||||||
free(bwt->bwt);
|
|
||||||
bwt->bwt = buf;
|
|
||||||
}
|
|
||||||
|
|
||||||
// 根据interval-bwt创建fmt-index
|
// 根据interval-bwt创建fmt-index
|
||||||
FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
|
FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
|
||||||
|
|
@ -523,160 +282,7 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
|
||||||
return fmt;
|
return fmt;
|
||||||
}
|
}
|
||||||
|
|
||||||
// an analogy to bwt_occ4() but more efficient, requiring k <= l
|
|
||||||
void bwt_2occ4(const bwt_t *bwt, bwtint_t k, bwtint_t l, bwtint_t cntk[4], bwtint_t cntl[4])
|
|
||||||
{
|
|
||||||
bwtint_t _k, _l;
|
|
||||||
_k = k - (k >= bwt->primary);
|
|
||||||
_l = l - (l >= bwt->primary);
|
|
||||||
if (_l >> OCC_INTV_SHIFT != _k >> OCC_INTV_SHIFT || k == (bwtint_t)(-1) || l == (bwtint_t)(-1))
|
|
||||||
{
|
|
||||||
bwt_occ4(bwt, k, cntk);
|
|
||||||
bwt_occ4(bwt, l, cntl);
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
bwtint_t x, y;
|
|
||||||
uint32_t *p, tmp, *endk, *endl;
|
|
||||||
k -= (k >= bwt->primary); // because $ is not in bwt
|
|
||||||
l -= (l >= bwt->primary);
|
|
||||||
p = bwt_occ_intv(bwt, k);
|
|
||||||
memcpy(cntk, p, 4 * sizeof(bwtint_t));
|
|
||||||
p += sizeof(bwtint_t); // sizeof(bwtint_t) = 4*(sizeof(bwtint_t)/sizeof(uint32_t))
|
|
||||||
// prepare cntk[]
|
|
||||||
endk = p + ((k >> 4) - ((k & ~OCC_INTV_MASK) >> 4));
|
|
||||||
endl = p + ((l >> 4) - ((l & ~OCC_INTV_MASK) >> 4));
|
|
||||||
for (x = 0; p < endk; ++p)
|
|
||||||
x += __occ_aux4(bwt, *p);
|
|
||||||
y = x;
|
|
||||||
tmp = *p & ~((1U << ((~k & 15) << 1)) - 1);
|
|
||||||
x += __occ_aux4(bwt, tmp) - (~k & 15);
|
|
||||||
// calculate cntl[] and finalize cntk[]
|
|
||||||
for (; p < endl; ++p)
|
|
||||||
y += __occ_aux4(bwt, *p);
|
|
||||||
tmp = *p & ~((1U << ((~l & 15) << 1)) - 1);
|
|
||||||
y += __occ_aux4(bwt, tmp) - (~l & 15);
|
|
||||||
memcpy(cntl, cntk, 4 * sizeof(bwtint_t));
|
|
||||||
cntk[0] += x & 0xff;
|
|
||||||
cntk[1] += x >> 8 & 0xff;
|
|
||||||
cntk[2] += x >> 16 & 0xff;
|
|
||||||
cntk[3] += x >> 24;
|
|
||||||
cntl[0] += y & 0xff;
|
|
||||||
cntl[1] += y >> 8 & 0xff;
|
|
||||||
cntl[2] += y >> 16 & 0xff;
|
|
||||||
cntl[3] += y >> 24;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
void bwt_extend(const bwt_t *bwt, const bwtintv_t *ik, bwtintv_t ok[4], int is_back)
|
|
||||||
{
|
|
||||||
bwtint_t tk[4], tl[4];
|
|
||||||
int i;
|
|
||||||
bwt_2occ4(bwt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl); // tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积
|
|
||||||
// 这里是反向扩展
|
|
||||||
for (i = 0; i != 4; ++i)
|
|
||||||
{
|
|
||||||
ok[i].x[!is_back] = bwt->L2[i] + 1 + tk[i]; // 起始行位置，互补链
|
|
||||||
ok[i].x[2] = tl[i] - tk[i]; // 间隔
|
|
||||||
}
|
|
||||||
// 因为计算的是互补碱基，所以3对应着0,2对应1，下边是正向扩展
|
|
||||||
ok[3].x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= bwt->primary && ik->x[!is_back] + ik->x[2] - 1 >= bwt->primary);
|
|
||||||
ok[2].x[is_back] = ok[3].x[is_back] + ok[3].x[2];
|
|
||||||
ok[1].x[is_back] = ok[2].x[is_back] + ok[2].x[2];
|
|
||||||
ok[0].x[is_back] = ok[1].x[is_back] + ok[1].x[2];
|
|
||||||
}
|
|
||||||
|
|
||||||
// 利用bwt搜索seed，完整搜索，只需要单向搜索
|
|
||||||
void bwt_search(bwt_t *bwt, const string &q)
|
|
||||||
{
|
|
||||||
bwtintv_t ik, ok[4];
|
|
||||||
int i, j, c, x = 0;
|
|
||||||
|
|
||||||
bwt_set_intv(bwt, bval(q[x]), ik);
|
|
||||||
ik.info = x + 1;
|
|
||||||
cout << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
|
||||||
|
|
||||||
for (i = x + 1; i < q.size(); ++i)
|
|
||||||
{
|
|
||||||
if (bval(q[i]) < 4)
|
|
||||||
{
|
|
||||||
c = cbval(q[i]);
|
|
||||||
bwt_extend(bwt, &ik, ok, 0);
|
|
||||||
ik = ok[c];
|
|
||||||
ik.info = i + 1;
|
|
||||||
cout << "bwt-1: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// 扩展两次
|
|
||||||
void bwt_extend2(const bwt_t *bwt, bwtintv_t *ik, bwtintv_t ok[4], int is_back, int c1)
|
|
||||||
{
|
|
||||||
bwtint_t tk[4], tl[4];
|
|
||||||
int i;
|
|
||||||
bwt_2occ4(bwt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl); // tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积
|
|
||||||
// 这里是反向扩展
|
|
||||||
for (i = 0; i != 4; ++i)
|
|
||||||
{
|
|
||||||
ok[i].x[!is_back] = bwt->L2[i] + 1 + tk[i]; // 起始行位置，互补链
|
|
||||||
ok[i].x[2] = tl[i] - tk[i]; // 间隔
|
|
||||||
}
|
|
||||||
// 因为计算的是互补碱基，所以3对应着0,2对应1，下边是正向扩展
|
|
||||||
ok[3].x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= bwt->primary && ik->x[!is_back] + ik->x[2] - 1 >= bwt->primary);
|
|
||||||
ok[2].x[is_back] = ok[3].x[is_back] + ok[3].x[2];
|
|
||||||
ok[1].x[is_back] = ok[2].x[is_back] + ok[2].x[2];
|
|
||||||
ok[0].x[is_back] = ok[1].x[is_back] + ok[1].x[2];
|
|
||||||
|
|
||||||
*ik = ok[c1];
|
|
||||||
bwt_2occ4(bwt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl);
|
|
||||||
for (i = 0; i != 4; ++i)
|
|
||||||
{
|
|
||||||
ok[i].x[!is_back] = bwt->L2[i] + 1 + tk[i]; // 起始行位置，互补链
|
|
||||||
ok[i].x[2] = tl[i] - tk[i]; // 间隔
|
|
||||||
}
|
|
||||||
// 因为计算的是互补碱基，所以3对应着0,2对应1，下边是正向扩展
|
|
||||||
ok[3].x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= bwt->primary && ik->x[!is_back] + ik->x[2] - 1 >= bwt->primary);
|
|
||||||
ok[2].x[is_back] = ok[3].x[is_back] + ok[3].x[2];
|
|
||||||
ok[1].x[is_back] = ok[2].x[is_back] + ok[2].x[2];
|
|
||||||
ok[0].x[is_back] = ok[1].x[is_back] + ok[1].x[2];
|
|
||||||
}
|
|
||||||
|
|
||||||
void bwt_search2(bwt_t *bwt, const string &q)
|
|
||||||
{
|
|
||||||
bwtintv_t ik, ok[4];
|
|
||||||
int i, j, c1, c2, x = 0;
|
|
||||||
|
|
||||||
bwt_set_intv(bwt, bval(q[x]), ik);
|
|
||||||
ik.info = x + 1;
|
|
||||||
cout << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
|
||||||
|
|
||||||
for (i = x + 1; i + 1 < q.size(); i += 2)
|
|
||||||
{
|
|
||||||
if (bval(q[i]) < 4 && bval(q[i+1]) < 4)
|
|
||||||
{
|
|
||||||
|
|
||||||
c1 = cbval(q[i]);
|
|
||||||
c2 = cbval(q[i + 1]);
|
|
||||||
|
|
||||||
bwt_extend2(bwt, &ik, ok, 0, c1);
|
|
||||||
|
|
||||||
ik = ok[c2];
|
|
||||||
ik.info = i + 1;
|
|
||||||
cout << "bwt-2: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
|
||||||
}
|
|
||||||
else
|
|
||||||
{
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if (i < q.size() && bval(q[i]) < 4) { // 最后一次扩展
|
|
||||||
c1 = cbval(q[i]);
|
|
||||||
bwt_extend(bwt, &ik, ok, 0);
|
|
||||||
ik = ok[c1];
|
|
||||||
ik.info = i + 1;
|
|
||||||
cout << "bwt-2: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
#define fmt_set_intv(fmt, c, ik) ((ik).x[0] = (fmt)->L2[(int)(c)] + 1, (ik).x[2] = (fmt)->L2[(int)(c) + 1] - (fmt)->L2[(int)(c)], (ik).x[1] = (fmt)->L2[3 - (c)] + 1, (ik).info = 0)
|
#define fmt_set_intv(fmt, c, ik) ((ik).x[0] = (fmt)->L2[(int)(c)] + 1, (ik).x[2] = (fmt)->L2[(int)(c) + 1] - (fmt)->L2[(int)(c)], (ik).x[1] = (fmt)->L2[3 - (c)] + 1, (ik).info = 0)
|
||||||
#define fmt_occ_intv(b, k) ((b)->bwt + (k) / OCC_INTERVAL * (OCC_INTERVAL / 8 + 20))
|
#define fmt_occ_intv(b, k) ((b)->bwt + (k) / OCC_INTERVAL * (OCC_INTERVAL / 8 + 20))
|
||||||
|
|
|
||||||
|
|
@ -0,0 +1,26 @@
|
||||||
|
#ifndef FMT_INDEX_H_
|
||||||
|
#define FMT_INDEX_H_
|
||||||
|
|
||||||
|
#include "bwt.h"
|
||||||
|
|
||||||
|
// fm-index, extend twice in one search step (one memory access)
|
||||||
|
struct FMTIndex
|
||||||
|
{
|
||||||
|
bwtint_t primary; // S^{-1}(0), or the primary index of BWT
|
||||||
|
bwtint_t sec_primary; // second primary line
|
||||||
|
bwtint_t L2[5]; // C(), cumulative count
|
||||||
|
bwtint_t seq_len; // sequence length
|
||||||
|
bwtint_t bwt_size; // size of bwt, about seq_len/4
|
||||||
|
uint32_t *bwt; // BWT
|
||||||
|
// occurance array, separated to two parts
|
||||||
|
uint32_t cnt_table[5][256]; // 4对应原来的cnt_table，0,1,2,3,分别对应该碱基的扩展
|
||||||
|
int sec_bcp; // base couple for sec primary line, AA=>0, AC=>1 ... TT=>15
|
||||||
|
int first_base; // 序列的第一个碱基2bit的int类型，0,1,2,3
|
||||||
|
int last_base; // dollar转换成的base
|
||||||
|
// suffix array
|
||||||
|
int sa_intv;
|
||||||
|
bwtint_t n_sa;
|
||||||
|
uint8_t *sa;
|
||||||
|
};
|
||||||
|
|
||||||
|
#endif
|
||||||
6
main.cpp
6
main.cpp
|
|
@ -4,9 +4,13 @@
|
||||||
#include <vector>
|
#include <vector>
|
||||||
#include <sys/time.h>
|
#include <sys/time.h>
|
||||||
|
|
||||||
#include "common.h"
|
|
||||||
using namespace std;
|
using namespace std;
|
||||||
|
|
||||||
|
// 各个分部的入口函数
|
||||||
|
int main_sa(int argc, char **argv);
|
||||||
|
int main_fmtidx(int argc, char **argv);
|
||||||
|
|
||||||
|
// 整个程序的入口
|
||||||
int main(int argc, char **argv)
|
int main(int argc, char **argv)
|
||||||
{
|
{
|
||||||
// main_sa(argc, argv);
|
// main_sa(argc, argv);
|
||||||
|
|
|
||||||
12
sa.cpp
12
sa.cpp
|
|
@ -2,20 +2,12 @@
|
||||||
#include <stdint.h>
|
#include <stdint.h>
|
||||||
#include <stdlib.h>
|
#include <stdlib.h>
|
||||||
#include <vector>
|
#include <vector>
|
||||||
#include <sys/time.h>
|
|
||||||
|
|
||||||
#include "common.h"
|
#include "util.h"
|
||||||
|
#include "sa.h"
|
||||||
|
|
||||||
using namespace std;
|
using namespace std;
|
||||||
|
|
||||||
double realtime(void)
|
|
||||||
{
|
|
||||||
struct timeval tp;
|
|
||||||
struct timezone tzp;
|
|
||||||
gettimeofday(&tp, &tzp);
|
|
||||||
return tp.tv_sec + tp.tv_usec * 1e-6;
|
|
||||||
}
|
|
||||||
|
|
||||||
inline void bwt_set_sa_33(uint8_t *sa_arr, bwtint_t k, bwtint_t val)
|
inline void bwt_set_sa_33(uint8_t *sa_arr, bwtint_t k, bwtint_t val)
|
||||||
{
|
{
|
||||||
const bwtint_t block_idx = (k >> 3) * 33; // 8个数为一组，共享33个字节
|
const bwtint_t block_idx = (k >> 3) * 33; // 8个数为一组，共享33个字节
|
||||||
|
|
|
||||||
|
|
@ -0,0 +1,14 @@
|
||||||
|
#ifndef SA_H_
|
||||||
|
#define SA_H_
|
||||||
|
|
||||||
|
#include "util.h"
|
||||||
|
|
||||||
|
void bwt_set_sa_33(uint8_t *sa_arr, bwtint_t k, bwtint_t val);
|
||||||
|
bwtint_t bwt_get_sa_33(uint8_t *sa_arr, bwtint_t k);
|
||||||
|
|
||||||
|
// 用33个bit来表示bwt行信息，所需的总字节数
|
||||||
|
#define SA_BYTES_33(n_sa) ((((33 * (n_sa) + 7) / 8) & (~7L)) + 8)
|
||||||
|
// 用40个bit来表示bwt行信息，所需的总字节数
|
||||||
|
#define SA_BYTES_40(n_sa) ((((40 * (n_sa) + 7) / 8) & (~7L)) + 8)
|
||||||
|
|
||||||
|
#endif
|
||||||
|
|
@ -0,0 +1,56 @@
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <sys/time.h>
|
||||||
|
|
||||||
|
#include "util.h"
|
||||||
|
|
||||||
|
// base转成2bit值
|
||||||
|
int bval(char b)
|
||||||
|
{
|
||||||
|
if (b == 'A')
|
||||||
|
return 0;
|
||||||
|
if (b == 'C')
|
||||||
|
return 1;
|
||||||
|
if (b == 'G')
|
||||||
|
return 2;
|
||||||
|
if (b == 'T')
|
||||||
|
return 3;
|
||||||
|
return 4;
|
||||||
|
}
|
||||||
|
|
||||||
|
// 互补碱基值
|
||||||
|
int cbval(char b)
|
||||||
|
{
|
||||||
|
return 3 - bval(b);
|
||||||
|
}
|
||||||
|
|
||||||
|
double realtime(void)
|
||||||
|
{
|
||||||
|
struct timeval tp;
|
||||||
|
struct timezone tzp;
|
||||||
|
gettimeofday(&tp, &tzp);
|
||||||
|
return tp.tv_sec + tp.tv_usec * 1e-6;
|
||||||
|
}
|
||||||
|
|
||||||
|
// 打印故障信息，并终止程序
|
||||||
|
void _err_fatal_simple_core(const char *func, const char *msg)
|
||||||
|
{
|
||||||
|
fprintf(stderr, "[%s] %s Abort!\n", func, msg);
|
||||||
|
abort();
|
||||||
|
}
|
||||||
|
|
||||||
|
// 读取数据
|
||||||
|
bwtint_t fread_fix(FILE *fp, bwtint_t size, void *a)
|
||||||
|
{ /* Mac/Darwin has a bug when reading data longer than 2GB. This function fixes this issue by reading data in small chunks */
|
||||||
|
const int bufsize = 0x1000000; // 16M block
|
||||||
|
bwtint_t offset = 0;
|
||||||
|
while (size)
|
||||||
|
{
|
||||||
|
int x = bufsize < size ? bufsize : size;
|
||||||
|
if ((x = fread((uint8_t *)a + offset, 1, x, fp)) == 0)
|
||||||
|
break;
|
||||||
|
size -= x;
|
||||||
|
offset += x;
|
||||||
|
}
|
||||||
|
return offset;
|
||||||
|
}
|
||||||
|
|
@ -0,0 +1,31 @@
|
||||||
|
#ifndef COMMON_H_
|
||||||
|
#define COMMON_H_
|
||||||
|
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <stdint.h>
|
||||||
|
|
||||||
|
typedef uint64_t bwtint_t;
|
||||||
|
|
||||||
|
#define xassert(cond, msg) \
|
||||||
|
if ((cond) == 0) \
|
||||||
|
_err_fatal_simple_core(__func__, msg)
|
||||||
|
|
||||||
|
double realtime(void);
|
||||||
|
|
||||||
|
// 在fm-indexv(或者bwt)查找过程中，记录结果
|
||||||
|
struct bwtintv_t
|
||||||
|
{
|
||||||
|
bwtint_t x[3], info; // x[0]表示正链位置，x[1]表示互补链位置，x[2]表示间隔长度，info 表示read的起始结束位置
|
||||||
|
};
|
||||||
|
|
||||||
|
// 读取二进制数据
|
||||||
|
bwtint_t fread_fix(FILE *fp, bwtint_t size, void *a);
|
||||||
|
// 给出问题信息并终止程序
|
||||||
|
void _err_fatal_simple_core(const char *func, const char *msg);
|
||||||
|
|
||||||
|
// base转成2bit值
|
||||||
|
int bval(char b);
|
||||||
|
// 互补碱基值
|
||||||
|
int cbval(char b);
|
||||||
|
|
||||||
|
#endif
|
||||||
Loading…
Reference in New Issue