From 5f181677036d3bd311df38fab250fcaebfcb697e Mon Sep 17 00:00:00 2001 From: zzh Date: Fri, 2 Feb 2024 03:59:35 +0800 Subject: [PATCH] =?UTF-8?q?mid=20interval=E6=94=AF=E6=8C=818=E5=92=8C16?= =?UTF-8?q?=EF=BC=8C=E4=BB=A3=E7=A0=81=E4=B8=AD=E5=90=AB=E6=9C=89=E6=B3=A8?= =?UTF-8?q?=E9=87=8A=E6=8E=89=E7=9A=84=E8=B0=83=E8=AF=95=E4=BF=A1=E6=81=AF?= =?UTF-8?q?=EF=BC=8C=E5=8F=91=E7=8E=B0=E5=B9=B6=E8=A7=A3=E5=86=B3=E4=BA=86?= =?UTF-8?q?=E4=B8=80=E4=B8=AA=E6=95=B0=E6=8D=AE=E6=BA=A2=E5=87=BA=E7=9A=84?= =?UTF-8?q?bug?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Makefile | 2 +- bwt.cpp | 24 ++++++++++----- fmt_index.cpp | 84 ++++++++++++++++++++++++++++++++++++++++----------- fmt_index.h | 4 ++- 4 files changed, 87 insertions(+), 27 deletions(-) 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 bd9f45e..b680eda 100644 --- a/bwt.cpp +++ b/bwt.cpp @@ -10,6 +10,8 @@ using namespace std; +const static char BASE[4] = {'A', 'C', 'G', 'T'}; + // 打印32位整型数据中包含的pre-bwt：bwt static void print_base_uint32(uint32_t p) { @@ -165,7 +167,7 @@ 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 << "k: " << k << "; occ: " << p[0] << ' ' << p[1] << ' ' << p[2] << ' ' << p[3] << endl; + //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]; @@ -176,16 +178,20 @@ void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4]) // cout << "k - kbase: " << k - bwt_k_base_line << endl; for (; p < end; ++p) { - x += __occ_aux4(bwt, *p); - //print_base_uint32(*p); + x += __occ_aux4(bwt, *p); + //if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550) + // print_base_uint32(*p); } tmp = *p & ~((1U << ((~k & 15) << 1)) - 1); - // print_base_uint32(tmp); + //if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550) + // print_base_uint32(tmp); 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; + //if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550) + // cout << "final-occ: " << cnt[0] << ' ' << cnt[1] << ' ' << cnt[2] << ' ' << cnt[3] << endl; } // an analogy to bwt_occ4() but more efficient, requiring k <= l @@ -210,8 +216,11 @@ void bwt_2occ4(const bwt_t *bwt, bwtint_t k, bwtint_t l, bwtint_t cntk[4], bwtin 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)) + cntk[0] = p[0]; + cntk[1] = p[1]; + cntk[2] = p[2]; + cntk[3] = p[3]; + p += 4; // 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)); @@ -356,12 +365,13 @@ bwtintv_t bwt_search2(bwt_t *bwt, const string &q) c2 = cbval(q[i + 1]); // double tm_t = realtime(); + //cout << BASE[c2] << BASE[c1] << endl; bwt_extend2(bwt, &ik, ok, 0, c1, c2); // t7 += realtime() - tm_t; ik = ok[c2]; 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; } else { diff --git a/fmt_index.cpp b/fmt_index.cpp index e994ef8..57031a2 100644 --- a/fmt_index.cpp +++ b/fmt_index.cpp @@ -19,7 +19,7 @@ double t1 = 0, t2 = 0, t3 = 0, t4 = 0, t5 = 0, t6 = 0, t7 = 0, t8 = 0, t9 = 0, t long f1 = 0, f2 = 0, f3 = 0, f4 = 0, f5 = 0; -const char BASE[4] = {'A', 'C', 'G', 'T'}; +const static char BASE[4] = {'A', 'C', 'G', 'T'}; // 求反向互补序列 string calc_reverse_seq(string &seq) @@ -207,6 +207,8 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt) bwtint_t last_seq_len = bwt->seq_len % FMT_OCC_INTERVAL; mn_occ += (last_seq_len + FMT_MID_INTERVAL - 1) / FMT_MID_INTERVAL - 1; fmt->bwt_size += mn_occ * 4; + //cout << "mn_occ: " << mn_occ << endl; + //exit(0); #endif buf = (uint32_t *)calloc(fmt->bwt_size, 4); // 开辟计算fmt用到的缓存 @@ -230,7 +232,9 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt) k += 4; memcpy(buf + k, c2, sizeof(uint32_t) * 16); // pre-bwt:bwt碱基对的occ k += 16; +#ifdef FMT_MID_INTERVAL mc[0] = mc[1] = mc[2] = mc[3] = 0; +#endif } // 每个32位整数保存8个倒数第二列碱基（pre-bwt）和8个倒数第一列(bwt)碱基 if (i % 16 == 0) // 每个32位整数可以包含16个碱基，每次需要处理16个碱基，也就是间隔最小可以设置为16 @@ -297,16 +301,36 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt) pre_and_bwt_seq |= (((pre_bwt_16_seq & (3 << lshift_bit)) >> (m * 2)) | ((bwt_16_seq & (3 << lshift_bit)) >> ((m * 2) + 2))); } buf[k++] = pre_and_bwt_seq; + if (j > 8) { + +#if FMT_MID_INTERVAL == 8 + for (m = 0; m < 4; ++m) + { + uint32_t s1 = pre_and_bwt_seq; + mc[m] += cnt_table[m][s1 & 0xff] + cnt_table[m][s1 >> 8 & 0xff] + cnt_table[m][s1 >> 16 & 0xff] + cnt_table[m][s1 >> 24]; + buf[k++] = mc[m]; + } +#endif for (m = 8; m > 0; --m) { int lshift_bit = 2 * m - 2; pre_and_bwt_seq_2 |= (((pre_bwt_16_seq & (3 << lshift_bit)) << (m * 2)) | ((bwt_16_seq & (3 << lshift_bit)) << (m * 2 - 2))); } buf[k++] = pre_and_bwt_seq_2; + +#if FMT_MID_INTERVAL == 8 + if ((i + 16) % FMT_OCC_INTERVAL != 0 && j == 16) + for (m = 0; m < 4; ++m) + { + uint32_t s1 = pre_and_bwt_seq_2; + mc[m] += cnt_table[m][s1 & 0xff] + cnt_table[m][s1 >> 8 & 0xff] + cnt_table[m][s1 >> 16 & 0xff] + cnt_table[m][s1 >> 24]; + buf[k++] = mc[m]; + } +#endif } -#ifdef FMT_MID_INTERVAL +#if FMT_MID_INTERVAL == 16 // 加入中间的check point for (m = 0; m < 4; ++m) { @@ -329,7 +353,7 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt) memcpy(buf + k, c2, sizeof(uint32_t) * 16); k += 16; // cout << "n occ: " << n_occ << endl; - // cout << "size: " << k << '\t' << fmt->bwt_size << endl; + cout << "size: " << k << '\t' << fmt->bwt_size << endl; xassert(k == fmt->bwt_size, "inconsistent bwt_size"); // update fmt fmt->bwt = buf; @@ -356,13 +380,17 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4] ti = b1 << 2 | b2; k -= (k >= fmt->primary); // k由bwt矩阵对应的行转换成bwt字符串对应的行（去掉了$，所以大于$的行，都减掉1） p = fmt_occ_intv(fmt, k); - // cout << "k-base: " << k << "; 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 = b1 + 1; i < 4; ++i) cnt[0] += p[i]; // 大于b1的碱基的occ之和 cnt[1] = p[b1]; // b1的occ q = p + 4 + b1 * 4; for (i = b2 + 1; i < 4 ; ++i) cnt[2] += q[i]; // 大于b2的occ之和 cnt[3] = q[b2]; // b2的occ p += 20; + + //if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550) { + // cout << "sec base-occ: " << q[0] << ' ' << q[1] << ' ' << q[2] << ' ' << q[3] << endl; + //} #ifdef FMT_MID_INTERVAL // 使用mid interval信息 int mk = k % FMT_OCC_INTERVAL; @@ -387,7 +415,7 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4] //} if (n_mintv > 0) // 至少超过了第一个mid interval { - p += n_mintv * 6 - 4; // 对应的mid interval check point的首地址，即A C G T的局部累积量 + p += n_mintv * (4 + (FMT_MID_INTERVAL >> 3)) - 4; // 对应的mid interval check point的首地址，即A C G T的局部累积量 q = p + b1; for (i = b1 + 1; i < 4; ++i) x += p[i]; // 大于b1的碱基的occ之和 @@ -406,18 +434,24 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4] x = 0; p += 4; } - + // cout << "mid-occ: " << cnt[0] << ' ' << cnt[1] << ' ' << cnt[2] << ' ' << cnt[3] << endl; // cout << "k: " << k << ' ' << cnt[1] << endl; - +#if FMT_MID_INTERVAL == 16 if ((mk & FMT_MID_INTV_MASK) >> 3) { x += __fmt_occ_e2_aux2(fmt, ti, *p); // print_base_uint32(*p); ++p; } +#endif tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); - // print_base_uint32(tmp); x += __fmt_occ_e2_aux2(fmt, ti, tmp); + //if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550) + //{ + // print_base_uint32(tmp); + // cout << "多的: " << (~k & 7) << endl; + // cout << (x >> 24 & 0xff) << ' ' << (x >> 16 & 0xff) << ' ' << (x >> 8 & 0xff) << ' ' << (x & 0xff) << endl; + //} #else // 该地址是bwt和pre_bwt字符串数据的首地址 uint32_t *end = p + ((k >> 3) - ((k & ~FMT_OCC_INTV_MASK) >> 3)); // this is the end point of the following loop // p = end - (end - p) / 8; @@ -437,19 +471,30 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4] x -= (~k & 7) << 8; if (b2 == 0) x -= (~k & 7) << 24; + // cout << "here-1" << endl; } // 如果跨过了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 << 16; + cnt[2] -= 1 ; else if (b2 == fmt->last_base) - x -= 1 << 24; + cnt[3] -= 1; + // cout << "here-2" << endl; + } +// if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550) +// { +// cout << BASE[b1] << BASE[b2] << endl; +// cout << "final-x: " << (x >> 24 & 0xff) << ' ' << (x >> 16 & 0xff) << ' ' << (x >> 8 & 0xff) << ' ' << (x & 0xff) << endl; +// } + cnt[0] += x & 0xff; cnt[1] += x >> 8 & 0xff; cnt[2] += x >> 16 & 0xff; cnt[3] += x >> 24 & 0xff; + // if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550) + // cout << "final-occ: " << cnt[0] << ' ' << cnt[1] << ' ' << cnt[2] << ' ' << cnt[3] << endl; } // 扩展两个碱基 @@ -501,7 +546,7 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q) fmt_set_intv(fmt, bval(q[x]), ik); ik.info = x + 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; for (i = x + 1; i + 1 < qlen; i += 2) { @@ -509,12 +554,13 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q) { c1 = cbval(q[i]); c2 = cbval(q[i + 1]); + //cout << BASE[c2] << BASE[c1] << endl; //double tm_t = realtime(); fmt_extend2(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 { @@ -576,7 +622,7 @@ int main_fmtidx(int argc, char **argv) // 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]) + ".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()); @@ -584,8 +630,9 @@ int main_fmtidx(int argc, char **argv) //dump_bwt(bwt_idx.c_str(), bwt); //string fmt_idx = string(argv[1]) + ".fmt"; + //string fmt_idx = string(argv[1]) + ".256.8.fmt"; string fmt_idx = string(argv[1]) + ".256.fmt"; - // string fmt_idx = string(argv[1]) + ".128.fmt"; + //string fmt_idx = string(argv[1]) + ".128.fmt"; // string fmt_idx = string(argv[1]) + ".64.fmt"; // string fmt_idx = string(argv[1]) + ".32.fmt"; @@ -595,18 +642,19 @@ int main_fmtidx(int argc, char **argv) cout << "[time read bwt:] " << t11 << "s" << endl; t10 = realtime(); FMTIndex *fmt = restore_fmt(fmt_idx.c_str()); - // FMTIndex *fmt = create_fmt_from_bwt(bwt); - // dump_fmt(fmt_idx.c_str(), fmt); + //FMTIndex *fmt = create_fmt_from_bwt(bwt); + //dump_fmt(fmt_idx.c_str(), fmt); t10 = realtime() - t10; cout << "[time gen fmt:] " << t10 << "s" << endl; vector seed_arr(10000000); seed_arr[0] = "GCGATACTAAGA"; + seed_arr[0] = "GAGAGCTGTTCCCGTGTTTTCCATGGTTT"; srand(time(NULL)); t1 = realtime(); for (int i = 0; i < (int)seed_arr.size(); ++i) - seed_arr[i] = generate_rand_seq(13); + seed_arr[i] = generate_rand_seq(7); t1 = realtime() - t1; cout << "[time gen seed:] " << t1 << "s" << endl; @@ -621,7 +669,7 @@ int main_fmtidx(int argc, char **argv) t3 = realtime() - t3; cout << "[time fmt search:] " << t3 << "s" << endl; - return 0; + //return 0; t4 = realtime(); for (int i = 0; i < (int)seed_arr.size(); ++i) diff --git a/fmt_index.h b/fmt_index.h index e825505..d4d58ca 100644 --- a/fmt_index.h +++ b/fmt_index.h @@ -12,9 +12,11 @@ // 获取碱基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的整数倍） -#ifdef FMT_MID_INTERVAL +#if FMT_MID_INTERVAL == 16 // 包含mid interval #define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 80)) +#elif FMT_MID_INTERVAL == 8 +#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 144)) #else #define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 20)) #endif