验证了fmt-index的正确性，用两种方式进行了实现，第二种记录b1和b2的occ以及大于他们的碱基的occ，更快一些
This commit is contained in:
parent
22b75d335b
commit
555e39df94
|
|
@ -10,7 +10,7 @@
|
|||
"request": "launch",
|
||||
"program": "${workspaceRoot}/fmtidx",
|
||||
"args": [
|
||||
"/home/zzh/data/reference/human_g1k_v37_decoy.fasta.bwt"
|
||||
"/mnt/d/data/reference/human_g1k_v37_decoy.fasta"
|
||||
],
|
||||
"cwd": "${workspaceFolder}", // 当前工作路径：当前文件所在的工作空间
|
||||
}
|
||||
|
|
|
|||
2
Makefile
2
Makefile
|
|
@ -1,7 +1,7 @@
|
|||
CC= g++
|
||||
NOWARN= -Wno-unused-result -Wno-unused-function
|
||||
CFLAGS= #-g -Wall $(NOWARN) #-O2
|
||||
CPPFLAGS= -g -Wall $(NOWARN) -O2
|
||||
CPPFLAGS= -g -Wall $(NOWARN) -O3
|
||||
WRAP_MALLOC=-DUSE_MALLOC_WRAPPERS
|
||||
SHOW_PERF= -DSHOW_PERF
|
||||
AR= ar
|
||||
|
|
|
|||
66
bwt.cpp
66
bwt.cpp
|
|
@ -150,9 +150,11 @@ void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4])
|
|||
}
|
||||
k -= (k >= bwt->primary); // because $ is not in bwt
|
||||
p = bwt_occ_intv(bwt, k);
|
||||
// cout << "occ: " << p[0] << ' ' << p[1] << ' ' << p[2] << ' ' << p[3] << endl;
|
||||
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
|
||||
//end = p + 1;
|
||||
for (x = 0; p < end; ++p)
|
||||
x += __occ_aux4(bwt, *p);
|
||||
tmp = *p & ~((1U << ((~k & 15) << 1)) - 1);
|
||||
|
|
@ -167,6 +169,9 @@ void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4])
|
|||
void bwt_2occ4(const bwt_t *bwt, bwtint_t k, bwtint_t l, bwtint_t cntk[4], bwtint_t cntl[4])
|
||||
{
|
||||
// double tm_t = realtime();
|
||||
bwt_occ4(bwt, k, cntk);
|
||||
bwt_occ4(bwt, l, cntl);
|
||||
return;
|
||||
bwtint_t _k, _l;
|
||||
_k = k - (k >= bwt->primary);
|
||||
_l = l - (l >= bwt->primary);
|
||||
|
|
@ -235,13 +240,13 @@ void bwt_extend(const bwt_t *bwt, const bwtintv_t *ik, bwtintv_t ok[4], int is_b
|
|||
void bwt_search(bwt_t *bwt, const string &q)
|
||||
{
|
||||
bwtintv_t ik, ok[4];
|
||||
int i, j, c, x = 0;
|
||||
int i, 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)
|
||||
for (i = x + 1; i < (int)q.size(); ++i)
|
||||
{
|
||||
if (bval(q[i]) < 4)
|
||||
{
|
||||
|
|
@ -255,48 +260,71 @@ void bwt_search(bwt_t *bwt, const string &q)
|
|||
}
|
||||
|
||||
// 扩展两次
|
||||
void bwt_extend2(const bwt_t *bwt, bwtintv_t *ik, bwtintv_t ok[4], int is_back, int c1)
|
||||
void bwt_extend2(const bwt_t *bwt, bwtintv_t *ik, bwtintv_t ok[4], int is_back, int c1, int c2)
|
||||
{
|
||||
bwtint_t tk[4], tl[4];
|
||||
bwtint_t interval = 0;
|
||||
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 = 3; i > c1; --i)
|
||||
//{
|
||||
// interval += tl[i] - tk[i];
|
||||
//}
|
||||
//cout << "interval-1: " << interval << endl;
|
||||
|
||||
// 这里是反向扩展
|
||||
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]; // 间隔
|
||||
}
|
||||
// ok[c1].x[!is_back] = bwt->L2[c1] + 1 + tk[c1]; // 起始行位置，互补链
|
||||
// ok[c1].x[2] = tl[c1] - tk[c1]; // 间隔
|
||||
// 因为计算的是互补碱基，所以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];
|
||||
//cout << (ik->x[!is_back] <= bwt->primary && ik->x[!is_back] + ik->x[2] - 1 >= bwt->primary) << endl;
|
||||
*ik = ok[c1];
|
||||
bwt_2occ4(bwt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl);
|
||||
//cout << (ik->x[!is_back] <= bwt->primary && ik->x[!is_back] + ik->x[2] - 1 >= bwt->primary) << ' ' << ok[c1].x[is_back] << endl;
|
||||
//interval = 0;
|
||||
//for (i = 3; i > c2; --i)
|
||||
//{
|
||||
// interval += tl[i] - tk[i];
|
||||
//}
|
||||
//cout << "interval-2: " << interval << endl;
|
||||
|
||||
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]; // 间隔
|
||||
}
|
||||
|
||||
// cout << "occ: " << ok[c2].x[2] << endl;
|
||||
|
||||
|
||||
// 因为计算的是互补碱基，所以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];
|
||||
|
||||
}
|
||||
|
||||
// 每次扩展两步
|
||||
bwtintv_t bwt_search2(bwt_t *bwt, const string &q)
|
||||
{
|
||||
bwtintv_t ik, ok[4];
|
||||
int i, j, c1, c2, x = 0;
|
||||
int i, 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)
|
||||
for (i = x + 1; i + 1 < (int)q.size(); i += 2)
|
||||
{
|
||||
if (bval(q[i]) < 4 && bval(q[i + 1]) < 4)
|
||||
{
|
||||
|
|
@ -305,7 +333,7 @@ bwtintv_t bwt_search2(bwt_t *bwt, const string &q)
|
|||
c2 = cbval(q[i + 1]);
|
||||
|
||||
// double tm_t = realtime();
|
||||
bwt_extend2(bwt, &ik, ok, 0, c1);
|
||||
bwt_extend2(bwt, &ik, ok, 0, c1, c2);
|
||||
// t7 += realtime() - tm_t;
|
||||
|
||||
ik = ok[c2];
|
||||
|
|
@ -317,7 +345,7 @@ bwtintv_t bwt_search2(bwt_t *bwt, const string &q)
|
|||
break;
|
||||
}
|
||||
}
|
||||
if (i < q.size() && bval(q[i]) < 4)
|
||||
if (i < (int)q.size() && bval(q[i]) < 4)
|
||||
{ // 最后一次扩展
|
||||
c1 = cbval(q[i]);
|
||||
// double tm_t = realtime();
|
||||
|
|
|
|||
201
fmt_index.cpp
201
fmt_index.cpp
|
|
@ -7,6 +7,7 @@
|
|||
#include <stdio.h>
|
||||
#include <algorithm>
|
||||
#include <string.h>
|
||||
#include <immintrin.h>
|
||||
|
||||
#include "util.h"
|
||||
#include "fmt_index.h"
|
||||
|
|
@ -112,6 +113,33 @@ void fmt_gen_cnt_table(FMTIndex *fmt)
|
|||
}
|
||||
}
|
||||
|
||||
// 生成occ，每个字节对应一个pattern
|
||||
void fmt_gen_cnt_occ(FMTIndex *fmt)
|
||||
{
|
||||
// 0-8：大于a的occ，8-16：大于b的occ，16-24：b的occ
|
||||
int i, a, b, ti;
|
||||
uint32_t oa, ooa, ob, oob;
|
||||
for (i = 0; i != 256; ++i) // 遍历单个字节的各种情况
|
||||
{
|
||||
for (a = 0; a < 4; ++a) // ba格式
|
||||
{
|
||||
oa = 0;
|
||||
ooa = 0;
|
||||
oa += ((i >> 4 & 3) == a) + ((i & 3) == a);
|
||||
ooa += ((i >> 4 & 3) > a) + ((i & 3) > a);
|
||||
for (b = 0; b < 4; ++b)
|
||||
{
|
||||
oob = ob = 0;
|
||||
oob += ((i >> 6 & 3) > b && (i >> 4 & 3) == a) + ((i >> 2 & 3) > b && (i & 3) == a);
|
||||
ob += ((i >> 6 & 3) == b && (i >> 4 & 3) == a) + ((i >> 2 & 3) == b && (i & 3) == a);
|
||||
|
||||
ti = a << 2 | b;
|
||||
fmt->cnt_occ[ti][i]= ob << 24 | oob << 16 | oa << 8 | ooa ;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void dump_fmt(const char *fn, const FMTIndex *fmt)
|
||||
{
|
||||
FILE *fp;
|
||||
|
|
@ -147,6 +175,7 @@ FMTIndex *restore_fmt(const char *fn)
|
|||
fmt->seq_len = fmt->L2[4];
|
||||
fclose(fp);
|
||||
fmt_gen_cnt_table(fmt); // 字节所能表示的各种碱基组合中，各个碱基的累积数量
|
||||
fmt_gen_cnt_occ(fmt);
|
||||
return fmt;
|
||||
}
|
||||
|
||||
|
|
@ -336,10 +365,10 @@ void fmt_e2_occ4(const FMTIndex *fmt, bwtint_t k, int b, uint32_t cnt1[4], uint3
|
|||
void fmt_e2_2occ4(const FMTIndex *fmt, bwtint_t k, bwtint_t l, int b,
|
||||
uint32_t cntk1[4], uint32_t cntk2[4], uint32_t cntl1[4], uint32_t cntl2[4])
|
||||
{
|
||||
// fmt_e2_occ4(fmt, k, b, cntk1, cntk2);
|
||||
// fmt_e2_occ4(fmt, l, b, cntl1, cntl2);
|
||||
// return;
|
||||
//double tm_t = realtime();
|
||||
fmt_e2_occ4(fmt, k, b, cntk1, cntk2);
|
||||
fmt_e2_occ4(fmt, l, b, cntl1, cntl2);
|
||||
return;
|
||||
// double tm_t = realtime();
|
||||
bwtint_t _k, _l;
|
||||
_k = k - (k >= fmt->primary); // 换算成了seq的行
|
||||
_l = l - (l >= fmt->primary);
|
||||
|
|
@ -418,7 +447,7 @@ void fmt_e2_2occ4(const FMTIndex *fmt, bwtint_t k, bwtint_t l, int b,
|
|||
// cout << "fmt-occ: " << l << '\t' << cntl1[0] << '\t' << cntl1[1] << '\t' << cntl1[2] << '\t' << cntl1[3] << endl;
|
||||
// cout << "fmt-occ-2: " << l << '\t' << cntl2[0] << '\t' << cntl2[1] << '\t' << cntl2[2] << '\t' << cntl2[3] << endl;
|
||||
}
|
||||
//t6 += realtime() - tm_t;
|
||||
// t6 += realtime() - tm_t;
|
||||
// f2 += 1;
|
||||
}
|
||||
|
||||
|
|
@ -453,7 +482,7 @@ void fmt_e1_occ4(const FMTIndex *fmt, bwtint_t k, uint32_t cnt[4])
|
|||
// 对k行和l行同时计算bwt str的occ，如果k和l落在同一个interval区间，可以减少一些计算量和访存
|
||||
void fmt_e1_2occ4(const FMTIndex *fmt, bwtint_t k, bwtint_t l, uint32_t cntk[4], uint32_t cntl[4])
|
||||
{
|
||||
//double tm_t = realtime();
|
||||
double tm_t = realtime();
|
||||
bwtint_t _k, _l;
|
||||
_k = k - (k >= fmt->primary); // 换算成了seq的行
|
||||
_l = l - (l >= fmt->primary);
|
||||
|
|
@ -500,11 +529,11 @@ void fmt_e1_2occ4(const FMTIndex *fmt, bwtint_t k, bwtint_t l, uint32_t cntk[4],
|
|||
cntl[2] += y1 >> 16 & 0xff;
|
||||
cntl[3] += y1 >> 24;
|
||||
}
|
||||
//t11 += realtime() - tm_t;
|
||||
t11 += realtime() - tm_t;
|
||||
}
|
||||
|
||||
// 扩展一个碱基
|
||||
void fmt_extend1(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t ok[4], int is_back, int b1)
|
||||
void fmt_extend1_1(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t ok[4], int is_back, int b1)
|
||||
{
|
||||
uint32_t tk[4], tl[4];
|
||||
int i;
|
||||
|
|
@ -525,7 +554,7 @@ void fmt_extend1(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t ok[4], int is_bac
|
|||
}
|
||||
|
||||
// 扩展两个碱基
|
||||
void fmt_extend2(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t ok[4], int is_back, int b1, int b2)
|
||||
void fmt_extend2_1(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t ok[4], int is_back, int b1, int b2)
|
||||
{
|
||||
uint32_t tk1[4], tl1[4], tk2[4], tl2[4];
|
||||
int i;
|
||||
|
|
@ -539,8 +568,8 @@ void fmt_extend2(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t ok[4], int is_bac
|
|||
ok[i].x[!is_back] = fmt->L2[i] + 1 + tk2[i]; // 起始行位置，互补链
|
||||
ok[i].x[2] = tl2[i] - tk2[i]; // 间隔
|
||||
}
|
||||
ok[b2].x[!is_back] = fmt->L2[b2] + 1 + tk2[b2];
|
||||
ok[b2].x[2] = tl2[b2] - tk2[b2];
|
||||
// ok[b2].x[!is_back] = fmt->L2[b2] + 1 + tk2[b2];
|
||||
// ok[b2].x[2] = tl2[b2] - tk2[b2];
|
||||
|
||||
// 因为计算的是互补碱基，所以3对应着0,2对应1，下边是正向扩展
|
||||
|
||||
|
|
@ -558,10 +587,130 @@ void fmt_extend2(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t ok[4], int is_bac
|
|||
|
||||
*ik = ok[b2];
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
// 扩展两个个碱基，计算bwt base为b的pre-bwt str中各个碱基的occ
|
||||
void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4])
|
||||
{
|
||||
uint32_t x;
|
||||
uint32_t *p, *q, tmp, *end;
|
||||
bwtint_t bwt_k_line = k, bwt_k_base_line = k >> OCC_INTV_SHIFT << OCC_INTV_SHIFT;
|
||||
int i, ti;
|
||||
cnt[0] = 0;
|
||||
cnt[1] = 0;
|
||||
cnt[2] = 0;
|
||||
if (k == (bwtint_t)(-1))
|
||||
{
|
||||
p = fmt->bwt + 4 + b1 * 4;
|
||||
for (i = 3; i > b2; --i)
|
||||
cnt[2] += p[i];
|
||||
cnt[3] = p[b2];
|
||||
return;
|
||||
}
|
||||
ti = b1 << 2 | b2;
|
||||
_mm_prefetch((const char *)(&fmt->cnt_occ[ti]), _MM_HINT_T0);
|
||||
k -= (k >= fmt->primary); // k由bwt矩阵对应的行转换成bwt字符串对应的行（去掉了$，所以大于$的行，都减掉1）
|
||||
p = fmt_occ_intv(fmt, k);
|
||||
// cout << "occ: " << p[0] << ' ' << p[1] << ' ' << p[2] << ' ' << p[3] << endl;
|
||||
for (i = 3; i > b1; --i) cnt[0] += p[i];
|
||||
cnt[1] = p[b1];
|
||||
q = p + 4 + b1 * 4;
|
||||
for (i = 3; i > b2; --i) cnt[2] += q[i];
|
||||
cnt[3] = q[b2];
|
||||
p += 20; // 该地址是bwt和pre_bwt字符串数据的首地址
|
||||
end = p + ((k >> 3) - ((k & ~OCC_INTV_MASK) >> 3)); // this is the end point of the following loop
|
||||
//if ((int)(end - p) > 8)
|
||||
// cout << "interval: " << (int)(end - p) << endl;
|
||||
|
||||
// end = p + 2 + (end - p) / 6;
|
||||
// end = p + 2;
|
||||
for (x = 0; p < end; ++p)
|
||||
{
|
||||
x += __fmt_occ_e2_aux2(fmt, ti, *p);
|
||||
}
|
||||
tmp = *p & ~((1U << ((~k & 7) << 2)) - 1);
|
||||
// print_base_uint32(tmp);
|
||||
x += __fmt_occ_e2_aux2(fmt, ti, tmp);
|
||||
if (b1 == 0)
|
||||
{
|
||||
x -= (~k & 7) << 8;
|
||||
if (b2 == 0)
|
||||
x -= (~k & 7) << 24;
|
||||
}
|
||||
|
||||
// 如果跨过了second_primary,那么可能需要减掉一次累积值
|
||||
if (b1 == fmt->first_base && bwt_k_base_line < fmt->sec_primary && bwt_k_line >= fmt->sec_primary)
|
||||
{
|
||||
if (b2 == fmt->last_base)
|
||||
x -= 1 << 24;
|
||||
else if (b2 < fmt->last_base)
|
||||
x -= 1 << 16;
|
||||
}
|
||||
|
||||
cnt[0] += x & 0xff;
|
||||
cnt[1] += x >> 8 & 0xff;
|
||||
cnt[2] += x >> 16 & 0xff;
|
||||
cnt[3] += x >> 24 & 0xff;
|
||||
}
|
||||
|
||||
void fmt_e2_2occ(const FMTIndex *fmt, bwtint_t k, bwtint_t l,
|
||||
int b1, int b2, bwtint_t cntk[4], bwtint_t cntl[4])
|
||||
{
|
||||
// double tm_t = realtime();
|
||||
fmt_e2_occ(fmt, k, b1, b2, cntk);
|
||||
fmt_e2_occ(fmt, l, b1, b2, cntl);
|
||||
return;
|
||||
bwtint_t _k, _l;
|
||||
_k = k - (k >= fmt->primary); // 换算成了seq的行
|
||||
_l = l - (l >= fmt->primary);
|
||||
if (_l >> OCC_INTV_SHIFT != _k >> OCC_INTV_SHIFT || k == (bwtint_t)(-1) || l == (bwtint_t)(-1))
|
||||
{
|
||||
fmt_e2_occ(fmt, k, b1, b2, cntk);
|
||||
fmt_e2_occ(fmt, l, b1, b2, cntl);
|
||||
}
|
||||
else
|
||||
{
|
||||
}
|
||||
}
|
||||
|
||||
// 扩展两个碱基
|
||||
void fmt_extend2(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t *ok, int is_back, int b1, int b2)
|
||||
{
|
||||
bwtint_t tk[4], tl[4], first_pos;
|
||||
// tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积
|
||||
fmt_e2_2occ(fmt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], b1, b2, tk, tl);
|
||||
|
||||
// 这里是反向扩展
|
||||
ok->x[!is_back] = fmt->L2[b2] + 1 + tk[3];
|
||||
ok->x[2] = tl[3] - tk[3];
|
||||
|
||||
// 第一次正向扩展
|
||||
ok->x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= fmt->primary && ik->x[!is_back] + ik->x[2] - 1 >= fmt->primary) + tl[0] - tk[0];
|
||||
// 第二次正向扩展
|
||||
// cout << "inteval-1: " << tl[0] - tk[0] << endl;
|
||||
// cout << (ik->x[!is_back] <= fmt->primary && ik->x[!is_back] + ik->x[2] - 1 >= fmt->primary) << endl;
|
||||
// cout << (ok->x[!is_back] <= fmt->primary && ok->x[!is_back] + ok->x[2] - 1 >= fmt->primary) << ' ' << ok->x[is_back] << endl;
|
||||
// cout << "inteval-2: " << tl[2] - tk[2] << endl;
|
||||
// cout << "occ: " << tl[3] - tk[3] << endl;
|
||||
|
||||
first_pos = fmt->L2[b1] + 1 + tk[1];
|
||||
ok->x[is_back] = ok->x[is_back] + (first_pos <= fmt->primary && first_pos + tl[1] - tk[1] - 1 >= fmt->primary) + tl[2] - tk[2];
|
||||
}
|
||||
|
||||
// 利用fmt搜索seed，完整搜索，只需要单向搜索
|
||||
bwtintv_t fmt_search(FMTIndex *fmt, const string &q)
|
||||
{
|
||||
bwtintv_t ik, ok[4];
|
||||
bwtintv_t ik;
|
||||
// bwtintv_t ok[4];
|
||||
bwtintv_t ok;
|
||||
int i, c1, c2, x = 0;
|
||||
int qlen = (int)q.size();
|
||||
|
||||
|
|
@ -577,10 +726,12 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q)
|
|||
c1 = cbval(q[i]);
|
||||
c2 = cbval(q[i + 1]);
|
||||
//double tm_t = realtime();
|
||||
fmt_extend2(fmt, &ik, ok, 0, c1, c2);
|
||||
//t8 += realtime() - tm_t;
|
||||
fmt_extend2(fmt, &ik, &ok, 0, c1, c2);
|
||||
//fmt_extend2_1(fmt, &ik, ok, 0, c1, c2);
|
||||
ik = ok;
|
||||
// t8 += realtime() - tm_t;
|
||||
ik.info = i + 1;
|
||||
//cout << "fmt : " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||
// cout << "fmt : " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
|
|
@ -591,7 +742,7 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q)
|
|||
{ // 最后一次扩展
|
||||
c1 = cbval(q[i]);
|
||||
//double tm_t = realtime();
|
||||
fmt_extend1(fmt, &ik, ok, 0, c1);
|
||||
// fmt_extend1(fmt, &ik, ok, 0, c1);
|
||||
//t9 += realtime() - tm_t;
|
||||
ik.info = i + 1;
|
||||
// cout << "fmt : " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||
|
|
@ -616,6 +767,7 @@ int main_fmtidx(int argc, char **argv)
|
|||
FMTIndex *fmt = restore_fmt(fmt_idx.c_str());
|
||||
|
||||
vector<string> seed_arr(10000000);
|
||||
seed_arr[0] = "CCGTCATCATCCG";
|
||||
srand(time(NULL));
|
||||
|
||||
t1 = realtime();
|
||||
|
|
@ -633,25 +785,24 @@ int main_fmtidx(int argc, char **argv)
|
|||
t4 = realtime();
|
||||
for (int i = 0; i < (int)seed_arr.size(); ++i)
|
||||
{
|
||||
// bwtintv_t p1 =
|
||||
bwtintv_t p1 =
|
||||
bwt_search2(bwt, seed_arr[i]);
|
||||
// bwtintv_t p2 = fmt_search(fmt, seed_arr[i]);
|
||||
// if (p1.x[1] != p2.x[1] || p1.x[2] != p2.x[2] || p1.x[3] != p2.x[3])
|
||||
// cout << seed_arr[i] << endl
|
||||
// << p1.x[1] << ' ' << p1.x[2] << ' ' << p1.x[3] << endl
|
||||
// << p2.x[1] << ' ' << p2.x[2] << ' ' << p2.x[3] << endl;
|
||||
bwtintv_t p2 = fmt_search(fmt, seed_arr[i]);
|
||||
if (p1.x[0] != p2.x[0] || p1.x[1] != p2.x[1] || p1.x[2] != p2.x[2])
|
||||
cout << seed_arr[i] << endl
|
||||
<< p1.x[0] << ' ' << p1.x[1] << ' ' << p1.x[2] << endl
|
||||
<< p2.x[0] << ' ' << p2.x[1] << ' ' << p2.x[2] << endl;
|
||||
}
|
||||
t4 = realtime() - t4;
|
||||
cout << "[time bwt search 2:] " << t4 << "s" << endl;
|
||||
|
||||
t3 = realtime();
|
||||
for (int i = 0; i < (int)seed_arr.size(); ++i)
|
||||
fmt_search(fmt, seed_arr[i]);
|
||||
t3 = realtime() - t3;
|
||||
cout << "[time fmt search:] " << t3 << "s" << endl;
|
||||
|
||||
cout << "bwt occ: " << t5 << "s; " << t7 << '\t' << t10 << endl;
|
||||
cout << "fmt occ: " << t6 << "s; " << t11 << '\t' << t8 << '\t' << t9 << endl;
|
||||
//cout << "bwt occ: " << t5 << "s; " << t7 << '\t' << t10 << endl;
|
||||
//cout << "fmt occ: " << t6 << "s; " << t11 << '\t' << t8 << '\t' << t9 << endl;
|
||||
// bwt_search(bwt, s);
|
||||
// bwt_search2(bwt, s);
|
||||
// for (int i = 0; i < 120; ++i)
|
||||
|
|
|
|||
|
|
@ -11,6 +11,9 @@
|
|||
#define __fmt_occ_e2_aux4(fmt, b, val) \
|
||||
((fmt)->cnt_table[(b)][(val) & 0xff] + (fmt)->cnt_table[b][(val) >> 8 & 0xff] + (fmt)->cnt_table[b][(val) >> 16 & 0xff] + (fmt)->cnt_table[b][(val) >> 24])
|
||||
|
||||
#define __fmt_occ_e2_aux2(fmt, b, val) \
|
||||
((fmt)->cnt_occ[(b)][(val) & 0xff] + (fmt)->cnt_occ[b][(val) >> 8 & 0xff] + (fmt)->cnt_occ[b][(val) >> 16 & 0xff] + (fmt)->cnt_occ[b][(val) >> 24])
|
||||
|
||||
// fm-index, extend twice in one search step (one memory access)
|
||||
struct FMTIndex
|
||||
{
|
||||
|
|
@ -22,6 +25,7 @@ struct FMTIndex
|
|||
uint32_t *bwt; // BWT
|
||||
// occurance array, separated to two parts
|
||||
uint32_t cnt_table[5][256]; // 4对应原来的cnt_table，0,1,2,3,分别对应该碱基的扩展
|
||||
uint32_t cnt_occ[16][256]; // 前16-24位表示b（碱基）的occ，8-16位表示大于b的occ，0-8表示大于a的occ，ba格式
|
||||
uint8_t sec_bcp; // base couple for sec primary line, AA=>0, AC=>1 ... TT=>15
|
||||
uint8_t first_base; // 序列的第一个碱基2bit的int类型，0,1,2,3
|
||||
uint8_t last_base; // dollar转换成的base
|
||||
|
|
|
|||
Loading…
Reference in New Issue