From fe34be5d3a5e8c6d7b65c843e8464ec43aef6edd Mon Sep 17 00:00:00 2001 From: zzh Date: Wed, 31 Jan 2024 02:19:08 +0800 Subject: [PATCH] =?UTF-8?q?=E8=A7=A3=E5=86=B3=E4=BA=86=E5=88=9B=E5=BB=BAfm?= =?UTF-8?q?t=E8=BF=87=E7=A8=8B=E4=B8=AD=E7=9A=84bug=EF=BC=8C=E6=8A=8Abwt?= =?UTF-8?q?=E4=B8=AD=E4=BF=9D=E5=AD=98occ=E7=9A=84=E7=B1=BB=E5=9E=8B?= =?UTF-8?q?=E4=BB=8Euint64=E6=94=B9=E6=88=90=E4=BA=86uint32?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .vscode/settings.json | 3 +- Makefile | 2 +- bwt.cpp | 97 +++++---- bwt.h | 8 +- fmt_index.cpp | 464 ++++++++---------------------------------- fmt_index.h | 6 +- 6 files changed, 155 insertions(+), 425 deletions(-) diff --git a/.vscode/settings.json b/.vscode/settings.json index 48bd760..e13a436 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,6 +1,7 @@ { "files.associations": { "ostream": "cpp", - "iostream": "cpp" + "iostream": "cpp", + "string": "cpp" } } \ No newline at end of file diff --git a/Makefile b/Makefile index ce158fc..27f3eef 100644 --- a/Makefile +++ b/Makefile @@ -1,7 +1,7 @@ CC= g++ NOWARN= -Wno-unused-result -Wno-unused-function CFLAGS= #-g -Wall $(NOWARN) #-O2 -CPPFLAGS= -g -Wall $(NOWARN) -O3 +CPPFLAGS= -g -Wall $(NOWARN) -O2 WRAP_MALLOC=-DUSE_MALLOC_WRAPPERS SHOW_PERF= -DSHOW_PERF AR= ar diff --git a/bwt.cpp b/bwt.cpp index e23d666..b2d2438 100644 --- a/bwt.cpp +++ b/bwt.cpp @@ -10,6 +10,17 @@ using namespace std; +// 打印32位整型数据中包含的pre-bwt：bwt +static void print_base_uint32(uint32_t p) +{ + const char BASE[4] = {'A', 'C', 'G', 'T'}; + for (int i = 30; i > 0; i -= 2) + { + int b1 = p >> i & 3; + cout << BASE[b1] << endl; + } +} + void dump_bwt(const char *fn, const bwt_t *bwt) { FILE *fp; @@ -63,26 +74,29 @@ bwt_t *restore_bwt(const char *fn) return bwt; } -// 根据原始的字符串bwt创建interval-bwt +// 根据原始的字符串bwt创建interval-bwt，用uint32_t来表示occ void create_interval_occ_bwt(bwt_t *bwt) { - bwtint_t i, k, c[4], n_occ; + bwtint_t i, k, n_occ; uint32_t *buf; + uint32_t c[4]; n_occ = (bwt->seq_len + OCC_INTERVAL - 1) / OCC_INTERVAL + 1; - bwt->bwt_size += n_occ * sizeof(bwtint_t); // the new size + bwt->bwt_size += n_occ * 4; // 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 % 1000000 == 0) + { + cout << "line: " << i << '\t' << ((uint64_t)c[0] + c[1] + c[2] + c[3]) << '\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; + memcpy(buf + k, c, sizeof(uint32_t) * 4); + k += 4; // in fact: sizeof(bwtint_t)=4*(sizeof(bwtint_t)/4) 每个c包含多少个32位 } if (i % 16 == 0) buf[k++] = bwt->bwt[i / 16]; // 16 == sizeof(uint32_t)/2, 2个bit表示一个碱基 @@ -90,8 +104,8 @@ void create_interval_occ_bwt(bwt_t *bwt) } // 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"); + memcpy(buf + k, c, sizeof(uint32_t) * 4); + xassert(k + 4 == bwt->bwt_size, "inconsistent bwt_size"); // update bwt free(bwt->bwt); bwt->bwt = buf; @@ -143,6 +157,7 @@ void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4]) { bwtint_t x; uint32_t *p, tmp, *end; + // bwtint_t bwt_k_base_line = k >> OCC_INTV_SHIFT << OCC_INTV_SHIFT; if (k == (bwtint_t)(-1)) { memset(cnt, 0, 4 * sizeof(bwtint_t)); @@ -150,14 +165,22 @@ 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)) + //cout << "k: " << k << "; occ: " << p[0] << ' ' << p[1] << ' ' << p[2] << ' ' << p[3] << endl; + cnt[0] = p[0]; + cnt[1] = p[1]; + cnt[2] = p[2]; + cnt[3] = p[3]; + p += 4; // 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; + //cout << "k - kbase: " << k - bwt_k_base_line << endl; for (x = 0; p < end; ++p) - x += __occ_aux4(bwt, *p); + { + x += __occ_aux4(bwt, *p); + //print_base_uint32(*p); + } tmp = *p & ~((1U << ((~k & 15) << 1)) - 1); + //print_base_uint32(tmp); x += __occ_aux4(bwt, tmp) - (~k & 15); // 这里多算了A，要减去 cnt[0] += x & 0xff; cnt[1] += x >> 8 & 0xff; @@ -263,14 +286,10 @@ 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, int c2) { bwtint_t tk[4], tl[4]; - bwtint_t interval = 0; + // oc1, oc2; + // int int1 = 0, int2 = 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) @@ -278,8 +297,7 @@ void bwt_extend2(const bwt_t *bwt, bwtintv_t *ik, bwtintv_t ok[4], int is_back, 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]; @@ -287,31 +305,36 @@ void bwt_extend2(const bwt_t *bwt, bwtintv_t *ik, bwtintv_t ok[4], int is_back, 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; + //oc1 = oc2 = 0; + //for (i = 3; i > c1; --i) + //{ + // int1 += tl[i] - tk[i]; + // oc1 += tk[i]; + // oc2 += tl[i]; + //} + //cout << "bwt-2 interval-1: " << int1 << ' ' << oc1 << ' ' << oc2 << endl; + + 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]; // 间隔 } - - // 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]; - + // oc1 = oc2 = 0; + // for (i = 3; i > c2; --i) + // { + // int2 += tl[i] - tk[i]; + // oc1 += tk[i]; + // oc2 += tl[i]; + // } + // cout << "bwt-2 interval-2: " << int2 << ' ' << oc1 << ' ' << oc2 << endl; } // 每次扩展两步 @@ -322,7 +345,7 @@ bwtintv_t bwt_search2(bwt_t *bwt, const string &q) 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; + // cout << "bwt-2: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl; for (i = x + 1; i + 1 < (int)q.size(); i += 2) { @@ -353,7 +376,7 @@ bwtintv_t bwt_search2(bwt_t *bwt, const string &q) // t10 += realtime() - tm_t; ik = ok[c1]; ik.info = i + 1; - // cout << "bwt-2: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl; + //cout << "bwt-2: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl; } // cout << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl; return ik; diff --git a/bwt.h b/bwt.h index 35caf27..e6c555a 100644 --- a/bwt.h +++ b/bwt.h @@ -7,7 +7,7 @@ using std::string; // occ间隔对应的右移位数，5表示间隔32行保存一次 -#define OCC_INTV_SHIFT 7 +#define OCC_INTV_SHIFT 6 #define OCC_INTERVAL (1LL << OCC_INTV_SHIFT) #define OCC_INTV_MASK (OCC_INTERVAL - 1) @@ -28,9 +28,9 @@ using std::string; // #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]) +#define bwt_bwt(b, k) ((b)->bwt[(k) / OCC_INTERVAL * (OCC_INTERVAL / 16 + 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)) +#define bwt_occ_intv(b, k) ((b)->bwt + (k) / OCC_INTERVAL * (OCC_INTERVAL / 16 + 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]) @@ -59,7 +59,7 @@ void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4]); bwt_t *restore_bwt(const char *fn); // 根据原始的字符串bwt创建interval-bwt void create_interval_occ_bwt(bwt_t *bwt); - +void dump_bwt(const char *fn, const bwt_t *bwt); // 利用bwt搜索seed，完整搜索，只需要单向搜索 void bwt_search(bwt_t *bwt, const string &q); // 每次扩展两步 diff --git a/fmt_index.cpp b/fmt_index.cpp index 147e5a1..da4913e 100644 --- a/fmt_index.cpp +++ b/fmt_index.cpp @@ -10,6 +10,7 @@ #include #include "util.h" +#include "bwt.h" #include "fmt_index.h" using namespace std; @@ -92,27 +93,6 @@ void create_bwt_mtx(string &seq) // cout << endl; } -// fmt-index的count table，4对应着bwt碱基的累积量，0,1,2,3分别对应着bwt是A,C,G,T，pre-bwt的累积量 -void fmt_gen_cnt_table(FMTIndex *fmt) -{ - int i, j, k; - for (i = 0; i != 256; ++i) // 遍历单个字节的各种情况 - { - uint32_t x = 0; - for (j = 0; j != 4; ++j) // 一个字节有8位，每个碱基是2位，所以一个字节包含4个碱基，从右向左数，第一个和第三个碱基数据bwt，第二个和第四个是对应的pre-bwt - x |= (((i & 3) == j) + ((i >> 4 & 3) == j)) << (j << 3); // 高位存pre-bwt，挨着存bwt，一一对应 - fmt->cnt_table[4][i] = x; // 保存单个字节中bwt碱基个数，每8位对应一个碱基的个数，从左到右依次是TGCA - - for (k = 0; k < 4;++k) // bwt碱基 - { - x = 0; // for [A,C,G,T][A,C,G,T] - for (j = 0; j != 4; ++j) // pre-bwt碱基 - x |= (((i >> 6 & 3) == j && (i >> 4 & 3) == k) + ((i >> 2 & 3) == j && (i & 3) == k)) << (j << 3); - fmt->cnt_table[k][i] = x; - } - } -} - // 生成occ，每个字节对应一个pattern void fmt_gen_cnt_occ(FMTIndex *fmt) { @@ -174,8 +154,7 @@ FMTIndex *restore_fmt(const char *fn) fread_fix(fp, fmt->bwt_size << 2, fmt->bwt); fmt->seq_len = fmt->L2[4]; fclose(fp); - fmt_gen_cnt_table(fmt); // 字节所能表示的各种碱基组合中，各个碱基的累积数量 - fmt_gen_cnt_occ(fmt); + fmt_gen_cnt_occ(fmt); // 字节所能表示的各种碱基组合中，各个碱基的累积数量 return fmt; } @@ -184,7 +163,7 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt) { // FILE *fmt_out = fopen("fmt.txt", "w"); FMTIndex *fmt = (FMTIndex *)calloc(1, sizeof(FMTIndex)); - fmt_gen_cnt_table(fmt); + fmt_gen_cnt_occ(fmt); bwtint_t i, j, k, m, n, n_occ, cnt[4], cnt2[4]; uint32_t c[4], c2[16]; /*c用来保存原来的bwt碱基串的累积值，c2用来保存pre-bwt和bwt碱基对的累计值，如AA..TT*/ @@ -195,7 +174,7 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt) fmt->L2[i] = bwt->L2[i]; // 每个碱基的总累积值 fmt->primary = bwt->primary; // $在末尾的行，在bwt matrix行中的排序位置 - n_occ = (bwt->seq_len + OCC_INTERVAL - 1) / OCC_INTERVAL + 1; // check point 个数 + n_occ = (bwt->seq_len + FMT_OCC_INTERVAL - 1) / FMT_OCC_INTERVAL + 1; // check point 个数 fmt->bwt_size = (fmt->seq_len * 2 + 15) >> 4; // 要保存最后两列碱基 fmt->bwt_size += n_occ * 20; // A,C,G,T和AA,AC.....TG,TT共20个 buf = (uint32_t *)calloc(fmt->bwt_size, 4); // 开辟计算fmt用到的缓存 @@ -213,7 +192,7 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt) for (i = k = 0; i < bwt->seq_len; ++i) { // 记录occ - if (i % OCC_INTERVAL == 0) + if (i % FMT_OCC_INTERVAL == 0) { memcpy(buf + k, c, sizeof(uint32_t) * 4); // bwt str中各个碱基的occ k += 4; @@ -224,8 +203,8 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt) if (i % 16 == 0) // 每个32位整数可以包含16个碱基，每次需要处理16个碱基，也就是间隔最小可以设置为16 { uint32_t pre_bwt_16_seq = 0; // 16个pre-bwt碱基串 - uint32_t *bwt_addr = bwt_occ_intv(bwt, i) + 8; // bwt字符串i对应的基准行，因为原始的bwt-cp（check point）包含由4个uint64_t(8个uint32_t)组成的occ信息 - int offset = (i % OCC_INTERVAL) / 16; // 每OCC_INTERVAL个碱基共享同一个基准地址，每16个碱基共用一个uint32整型，因此需要偏移量来获取当前碱基串的首地址 + uint32_t *bwt_addr = bwt_occ_intv(bwt, i) + 4; // 这里加4还是加8要看保存occ的是是uint32还是uint64，bwt字符串i对应的基准行，因为原始的bwt-cp（check point）包含由4个uint32_t(8个uint32_t)组成的occ信息 + int offset = (i % FMT_OCC_INTERVAL) / 16; // 每OCC_INTERVAL个碱基共享同一个基准地址，每16个碱基共用一个uint32整型，因此需要偏移量来获取当前碱基串的首地址 uint32_t bwt_16_seq = *(bwt_addr + offset); // 待处理的当前16个碱基串的首地址 for (j = 0; j < 16; ++j) // 对于bwt碱基串，一个一个碱基分别处理 { @@ -309,300 +288,12 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt) return fmt; } -// 扩展两个个碱基，计算bwt base为b的pre-bwt str中各个碱基的occ -void fmt_e2_occ4(const FMTIndex *fmt, bwtint_t k, int b, uint32_t cnt1[4], uint32_t cnt2[4]) -{ - uint32_t x1, x2; - uint32_t *p, tmp, *end; - bwtint_t bwt_k_line = k, bwt_k_base_line = k >> OCC_INTV_SHIFT << OCC_INTV_SHIFT; - if (k == (bwtint_t)(-1)) - { - p = fmt->bwt + 4 + b * 4; - memset(cnt1, 0, 4 * sizeof(uint32_t)); - memcpy(cnt2, p, 4 * sizeof(uint32_t)); - return; - } - k -= (k >= fmt->primary); // k由bwt矩阵对应的行转换成bwt字符串对应的行（去掉了$，所以大于$的行，都减掉1） - p = fmt_occ_intv(fmt, k); - memcpy(cnt1, p, 4 * sizeof(uint32_t)); - memcpy(cnt2, p + 4 + b * 4, 4 * sizeof(uint32_t)); - p += 20; // 该地址是bwt和pre_bwt字符串数据的首地址 - end = p + ((k >> 3) - ((k & ~OCC_INTV_MASK) >> 3)); // this is the end point of the following loop - - for (x1 = 0, x2 = 0; p < end; ++p) - { - x1 += __fmt_occ_e2_aux4(fmt, 4, *p); - x2 += __fmt_occ_e2_aux4(fmt, b, *p); - } - tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); - x1 += __fmt_occ_e2_aux4(fmt, 4, tmp) - (~k & 7); - x2 += __fmt_occ_e2_aux4(fmt, b, tmp); - if (b == 0) - x2 -= ~k & 7; - // 如果跨过了second_primary,那么可能需要减掉一次累积值 - if (b == fmt->first_base && bwt_k_base_line < fmt->sec_primary && bwt_k_line >= fmt->sec_primary) - { - x2 -= 1 << (fmt->last_base << 3); - } - - cnt1[0] += x1 & 0xff; - cnt1[1] += x1 >> 8 & 0xff; - cnt1[2] += x1 >> 16 & 0xff; - cnt1[3] += x1 >> 24; - cnt2[0] += x2 & 0xff; - cnt2[1] += x2 >> 8 & 0xff; - cnt2[2] += x2 >> 16 & 0xff; - cnt2[3] += x2 >> 24; - - // cout << "fmt-occ: " << k << '\t' << cnt1[0] << '\t' << cnt1[1] << '\t' << cnt1[2] << '\t' << cnt1[3] << endl; - // cout << "fmt-occ-2: " << k << '\t' << cnt2[0] << '\t' << cnt2[1] << '\t' << cnt2[2] << '\t' << cnt2[3] << endl; - // cout << "bwt_k_base_line: " << bwt_k_base_line << endl; - // cout << "bwt_k_line: " << bwt_k_line << endl; - // cout << "sec_primary: " << fmt->sec_primary << endl; -} - -// 对k行和l行同时计算occ，如果k和l落在同一个interval区间，可以减少一些计算量和访存 -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(); - 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_occ4(fmt, k, b, cntk1, cntk2); - fmt_e2_occ4(fmt, l, b, cntl1, cntl2); - } - else - { - uint32_t x1, x2, y1, y2; - uint32_t *p, tmp, *ek, *el; - bwtint_t bwt_k_line = k, bwt_l_line = l, bwt_base_line = k >> OCC_INTV_SHIFT << OCC_INTV_SHIFT; - - k -= (k >= fmt->primary); // because $ is not in bwt - l -= (l >= fmt->primary); - p = fmt_occ_intv(fmt, k); - memcpy(cntk1, p, 4 * sizeof(uint32_t)); - memcpy(cntk2, p + 4 + b * 4, 4 * sizeof(uint32_t)); - memcpy(cntl1, cntk1, 4 * sizeof(uint32_t)); - memcpy(cntl2, cntk2, 4 * sizeof(uint32_t)); - p += 20; - // prepare cntk[] - ek = p + ((k >> 3) - ((k & ~OCC_INTV_MASK) >> 3)); - el = p + ((l >> 3) - ((l & ~OCC_INTV_MASK) >> 3)); - for (x1 = 0, x2 = 0; p < ek; ++p) - { - x1 += __fmt_occ_e2_aux4(fmt, 4, *p); - x2 += __fmt_occ_e2_aux4(fmt, b, *p); - } - y1 = x1; - y2 = x2; - tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); - x1 += __fmt_occ_e2_aux4(fmt, 4, tmp) - (~k & 7); - x2 += __fmt_occ_e2_aux4(fmt, b, tmp); - if (b == 0) - x2 -= ~k & 7; - for (; p < el; ++p) - { - y1 += __fmt_occ_e2_aux4(fmt, 4, *p); - y2 += __fmt_occ_e2_aux4(fmt, b, *p); - } - tmp = *p & ~((1U << ((~l & 7) << 2)) - 1); - y1 += __fmt_occ_e2_aux4(fmt, 4, tmp) - (~l & 7); - y2 += __fmt_occ_e2_aux4(fmt, b, tmp); - if (b == 0) - y2 -= ~l & 7; - // 如果跨过了second_primary,那么可能需要减掉一次累积值 - if (b == fmt->first_base && bwt_base_line < fmt->sec_primary) - { - if (bwt_k_line >= fmt->sec_primary) - x2 -= 1 << (fmt->last_base << 3); - if (bwt_l_line >= fmt->sec_primary) - y2 -= 1 << (fmt->last_base << 3); - } - cntk1[0] += x1 & 0xff; - cntk1[1] += x1 >> 8 & 0xff; - cntk1[2] += x1 >> 16 & 0xff; - cntk1[3] += x1 >> 24; - cntk2[0] += x2 & 0xff; - cntk2[1] += x2 >> 8 & 0xff; - cntk2[2] += x2 >> 16 & 0xff; - cntk2[3] += x2 >> 24; - - cntl1[0] += y1 & 0xff; - cntl1[1] += y1 >> 8 & 0xff; - cntl1[2] += y1 >> 16 & 0xff; - cntl1[3] += y1 >> 24; - cntl2[0] += y2 & 0xff; - cntl2[1] += y2 >> 8 & 0xff; - cntl2[2] += y2 >> 16 & 0xff; - cntl2[3] += y2 >> 24; - - // cout << "fmt-occ: " << k << '\t' << cntk1[0] << '\t' << cntk1[1] << '\t' << cntk1[2] << '\t' << cntk1[3] << endl; - // cout << "fmt-occ-2: " << k << '\t' << cntk2[0] << '\t' << cntk2[1] << '\t' << cntk2[2] << '\t' << cntk2[3] << endl; - // 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; - // f2 += 1; -} - -// 扩展一个碱基，计算bwt str中各个碱基的occ -void fmt_e1_occ4(const FMTIndex *fmt, bwtint_t k, uint32_t cnt[4]) -{ - uint32_t x; - uint32_t *p, tmp, *end; - if (k == (bwtint_t)(-1)) - { - memset(cnt, 0, 4 * sizeof(uint32_t)); - return; - } - k -= (k >= fmt->primary); // k由bwt矩阵对应的行转换成bwt字符串对应的行（去掉了$，所以大于$的行，都减掉1） - p = fmt_occ_intv(fmt, k); - memcpy(cnt, p, 4 * sizeof(uint32_t)); - p += 20; // 该地址是bwt和pre_bwt字符串数据的首地址 - end = p + ((k >> 3) - ((k & ~OCC_INTV_MASK) >> 3)); // this is the end point of the following loop - - for (x = 0; p < end; ++p) - { - x += __fmt_occ_e2_aux4(fmt, 4, *p); - } - tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); - x += __fmt_occ_e2_aux4(fmt, 4, tmp) - (~k & 7); - cnt[0] += x & 0xff; - cnt[1] += x >> 8 & 0xff; - cnt[2] += x >> 16 & 0xff; - cnt[3] += x >> 24; -} - -// 对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(); - 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_e1_occ4(fmt, k, cntk); - fmt_e1_occ4(fmt, l, cntl); - } - else - { - uint32_t x1, y1; - uint32_t *p, tmp, *ek, *el; - - k -= (k >= fmt->primary); // because $ is not in bwt - l -= (l >= fmt->primary); - p = fmt_occ_intv(fmt, k); - memcpy(cntk, p, 4 * sizeof(uint32_t)); - memcpy(cntl, p, 4 * sizeof(uint32_t)); - p += 20; - // prepare cntk[] - ek = p + ((k >> 3) - ((k & ~OCC_INTV_MASK) >> 3)); - el = p + ((l >> 3) - ((l & ~OCC_INTV_MASK) >> 3)); - for (x1 = 0; p < ek; ++p) - { - x1 += __fmt_occ_e2_aux4(fmt, 4, *p); - } - y1 = x1; - tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); - x1 += __fmt_occ_e2_aux4(fmt, 4, tmp) - (~k & 7); - for (; p < el; ++p) - { - y1 += __fmt_occ_e2_aux4(fmt, 4, *p); - } - tmp = *p & ~((1U << ((~l & 7) << 2)) - 1); - y1 += __fmt_occ_e2_aux4(fmt, 4, tmp) - (~l & 7); - - cntk[0] += x1 & 0xff; - cntk[1] += x1 >> 8 & 0xff; - cntk[2] += x1 >> 16 & 0xff; - cntk[3] += x1 >> 24; - - cntl[0] += y1 & 0xff; - cntl[1] += y1 >> 8 & 0xff; - cntl[2] += y1 >> 16 & 0xff; - cntl[3] += y1 >> 24; - } - t11 += realtime() - tm_t; -} - -// 扩展一个碱基 -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; - - // fmt_e1_occ4(fmt, ik->x[!is_back] - 1, tk); - // fmt_e1_occ4(fmt, ik->x[!is_back] - 1 + ik->x[2], tl); - fmt_e1_2occ4(fmt, 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] = fmt->L2[i] + 1 + tk[i]; // 起始行位置，互补链 - ok[i].x[2] = tl[i] - tk[i]; // 间隔 - } - ok[3].x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= fmt->primary && ik->x[!is_back] + ik->x[2] - 1 >= fmt->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[b1]; -} - -// 扩展两个碱基 -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; - - // tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积 - fmt_e2_2occ4(fmt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], b1, tk1, tk2, tl1, tl2); - - // 这里是反向扩展 - for (i = 0; i != 4; ++i) - { - 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]; - - // 因为计算的是互补碱基，所以3对应着0,2对应1，下边是正向扩展 - - ok[3].x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= fmt->primary && ik->x[!is_back] + ik->x[2] - 1 >= fmt->primary); - ok[2].x[is_back] = ok[3].x[is_back] + tl1[3] - tk1[3]; - ok[1].x[is_back] = ok[2].x[is_back] + tl1[2] - tk1[2]; - ok[0].x[is_back] = ok[1].x[is_back] + tl1[1] - tk1[1]; - - // cout << "fmt-d: " << BASE[b1] << '\t' << ok[b1].x[is_back] << '\t' << ok[b1].x[2] << endl; - - ok[3].x[is_back] = ok[b1].x[is_back] + (ok[b1].x[!is_back] <= fmt->primary && ok[b1].x[!is_back] + ok[b1].x[2] - 1 >= fmt->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[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; + bwtint_t bwt_k_line = k, bwt_k_base_line = k >> FMT_OCC_INTV_SHIFT << FMT_OCC_INTV_SHIFT; int i, ti; cnt[0] = 0; cnt[1] = 0; @@ -610,42 +301,47 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4] if (k == (bwtint_t)(-1)) { p = fmt->bwt + 4 + b1 * 4; - for (i = 3; i > b2; --i) - cnt[2] += p[i]; + 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); + // _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; + // cout << "k-base: " << k << "; 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; + end = p + ((k >> 3) - ((k & ~FMT_OCC_INTV_MASK) >> 3)); // this is the end point of the following loop + // p = end - (end - p) / 4; + // cout << "k - kbase: " << k - bwt_k_base_line << endl; for (x = 0; p < end; ++p) { x += __fmt_occ_e2_aux2(fmt, ti, *p); + // print_base_uint32(*p); } tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); // print_base_uint32(tmp); x += __fmt_occ_e2_aux2(fmt, ti, tmp); + // end = p + (end - p) / 4; + // end = p + 2; + // p = end - (end - p) / 32; + // if (k % 2 == 0) + // x += __fmt_occ_e2_aux2(fmt, ti, *p); + // p += 1; + // tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); + // 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) { @@ -654,69 +350,62 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4] 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); - + fmt_e2_occ(fmt, ik->x[!is_back] - 1, b1, b2, tk); + fmt_e2_occ(fmt, ik->x[!is_back] - 1 + ik->x[2], b1, b2, tl); + //fmt_e2_occ_2way(fmt, ik->x[!is_back] - 1, b1, b2, tk); + //fmt_e2_occ_2way(fmt, ik->x[!is_back] - 1 + ik->x[2], b1, b2, tl); // 这里是反向扩展 + //cout << BASE[b1] << BASE[b2] << endl; + //cout << "fmt: interval-1: " << tl[0] - tk[0] << ' ' << tl[0] << ' ' << tk[0] << endl; + //cout << "fmt: interval-2: " << tl[2] - tk[2] << ' ' << tl[2] << ' ' << tk[2] << endl; 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]; } +// 扩展一个碱基 +void fmt_extend1(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t *ok, int is_back, int b1) +{ + bwtint_t tk[4], tl[4]; + int b2 = 3; + // tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积 + fmt_e2_occ(fmt, ik->x[!is_back] - 1, b1, b2, tk); + fmt_e2_occ(fmt, ik->x[!is_back] - 1 + ik->x[2], b1, b2, tl); + //fmt_e2_occ_2way(fmt, ik->x[!is_back] - 1, b1, b2, tk); + //fmt_e2_occ_2way(fmt, ik->x[!is_back] - 1 + ik->x[2], b1, b2, tl); + // 这里是反向扩展 + ok->x[!is_back] = fmt->L2[b1] + 1 + tk[1]; + ok->x[2] = tl[1] - tk[1]; + // 第一次正向扩展 + 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]; +} + // 利用fmt搜索seed，完整搜索，只需要单向搜索 bwtintv_t fmt_search(FMTIndex *fmt, const string &q) { bwtintv_t ik; - // bwtintv_t ok[4]; bwtintv_t ok; int i, c1, c2, x = 0; int qlen = (int)q.size(); fmt_set_intv(fmt, bval(q[x]), ik); ik.info = x + 1; - // cout << 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; for (i = x + 1; i + 1 < qlen; i += 2) { @@ -727,7 +416,6 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q) c2 = cbval(q[i + 1]); //double tm_t = realtime(); 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; @@ -741,11 +429,12 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q) if (i < qlen && bval(q[i]) < 4) { // 最后一次扩展 c1 = cbval(q[i]); - //double tm_t = realtime(); - // fmt_extend1(fmt, &ik, ok, 0, c1); - //t9 += realtime() - tm_t; + // double tm_t = realtime(); + fmt_extend1(fmt, &ik, &ok, 0, c1); + ik = ok; + // t9 += 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; } // cout << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl; return ik; @@ -760,14 +449,29 @@ int main_fmtidx(int argc, char **argv) //create_bwt_mtx(seq); //cout << seq << endl; - string bwt_idx = string(argv[1]) + ".bwt"; - string fmt_idx = string(argv[1]) + ".fmt"; + //string bwt_str = string(argv[1]) + ".bwt.str"; + // string bwt_idx = string(argv[1]) + ".128.bwt"; + string bwt_idx = string(argv[1]) + ".64.bwt"; + // string bwt_idx = string(argv[1]) + ".16.bwt"; + // string bwt_idx = string(argv[1]) + ".bwt"; + //bwt_t *bwt = restore_bwt(bwt_str.c_str()); + //create_interval_occ_bwt(bwt); + //dump_bwt(bwt_idx.c_str(), bwt); + + //string fmt_idx = string(argv[1]) + ".fmt"; + // string fmt_idx = string(argv[1]) + ".256.fmt"; + string fmt_idx = string(argv[1]) + ".64.fmt"; + // string fmt_idx = string(argv[1]) + ".32.fmt"; + + bwt_t *bwt = restore_bwt(bwt_idx.c_str()); FMTIndex *fmt = restore_fmt(fmt_idx.c_str()); + // FMTIndex *fmt = create_fmt_from_bwt(bwt); + // dump_fmt(fmt_idx.c_str(), fmt); vector seed_arr(10000000); - seed_arr[0] = "CCGTCATCATCCG"; + seed_arr[0] = "GCGATACTAAGA"; srand(time(NULL)); t1 = realtime(); @@ -776,16 +480,21 @@ int main_fmtidx(int argc, char **argv) t1 = realtime() - t1; cout << "[time gen seed:] " << t1 << "s" << endl; - // t2 = realtime(); - // for (int i = 0; i < seed_arr.size(); ++i) - // bwt_search(bwt, seed_arr[i]); - // t2 = realtime() - t2; - // cout << "[time bwt search:] " << t2 << "s" << endl; + t2 = realtime(); + for (int i = 0; i < (int)seed_arr.size(); ++i) + bwt_search(bwt, seed_arr[i]); + t2 = realtime() - t2; + cout << "[time bwt search:] " << t2 << "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; 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[0] != p2.x[0] || p1.x[1] != p2.x[1] || p1.x[2] != p2.x[2]) @@ -795,12 +504,7 @@ int main_fmtidx(int argc, char **argv) } 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; // bwt_search(bwt, s); diff --git a/fmt_index.h b/fmt_index.h index b24cab3..b5fedab 100644 --- a/fmt_index.h +++ b/fmt_index.h @@ -1,12 +1,14 @@ #ifndef FMT_INDEX_H_ #define FMT_INDEX_H_ -#include "bwt.h" +#define FMT_OCC_INTV_SHIFT 6 +#define FMT_OCC_INTERVAL (1LL << FMT_OCC_INTV_SHIFT) +#define FMT_OCC_INTV_MASK (FMT_OCC_INTERVAL - 1) // 获取碱基c（待查找序列的首个碱基）和对应的互补碱基对应的行，以及间隔 #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) // k行（bwt str行（不包含$））对应的check point occ数据起始地址（小于k且是OCC_INTERVAL的整数倍） -#define fmt_occ_intv(b, k) ((b)->bwt + (k) / OCC_INTERVAL * (OCC_INTERVAL / 8 + 20)) +#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 20)) // 字节val中包含bwt base为b的pre-bwt中T G C A（按顺序保存在32位整数里（每个占8bit））的数量 #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])