diff --git a/.vscode/launch.json b/.vscode/launch.json index 279cc62..22ae4d4 100644 --- a/.vscode/launch.json +++ b/.vscode/launch.json @@ -10,7 +10,7 @@ "request": "launch", "program": "${workspaceRoot}/fmtidx", "args": [ - "/mnt/d/data/reference/human_g1k_v37_decoy.fasta" + "/home/zzh/reference/human_g1k_v37_decoy.fasta" ], "cwd": "${workspaceFolder}", // 当前工作路径:当前文件所在的工作空间 } diff --git a/bwt.cpp b/bwt.cpp index 5454d1e..5ed0cb9 100644 --- a/bwt.cpp +++ b/bwt.cpp @@ -354,7 +354,7 @@ bwtintv_t bwt_search2(bwt_t *bwt, const string &q) bwt_set_intv(bwt, bval(q[x]), ik); ik.info = x + 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; for (i = x + 1; i + 1 < (int)q.size(); i += 2) { diff --git a/fmt_index.cpp b/fmt_index.cpp index 3cb1def..30114ac 100644 --- a/fmt_index.cpp +++ b/fmt_index.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include "util.h" #include "bwt.h" @@ -203,12 +204,16 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt) fmt->bwt_size += n_occ * 20; // A,C,G,T和AA,AC.....TG,TT共20个 #ifdef FMT_MID_INTERVAL + uint32_t s1; bwtint_t mn_occ = (bwt->seq_len >> FMT_OCC_INTV_SHIFT) * (FMT_OCC_INTERVAL / FMT_MID_INTERVAL - 1); 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); + cout << (FMT_OCC_INTERVAL / FMT_MID_INTERVAL - 1) << ' ' << last_seq_len << ' ' << (last_seq_len + FMT_MID_INTERVAL - 1) / FMT_MID_INTERVAL - 1 << endl; + cout << "mn_occ: " << mn_occ << ' ' << mn_occ * 4 << endl; + i = 0; + //cout << ((i + 16 & FMT_MID_INTV_MASK) == 0) << endl; + // exit(0); #endif buf = (uint32_t *)calloc(fmt->bwt_size, 4); // 开辟计算fmt用到的缓存 @@ -304,47 +309,35 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt) 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; + +#ifdef FMT_MID_INTERVAL + s1 = pre_and_bwt_seq; + for (m = 0; m < 4; ++m) + mc[m] += cnt_table[m][s1 & 0xff] + cnt_table[m][s1 >> 8 & 0xff] + cnt_table[m][s1 >> 16 & 0xff] + cnt_table[m][s1 >> 24]; +#endif #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 - } -#if FMT_MID_INTERVAL == 16 - // 加入中间的check point - 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]; - 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]; - } - if ((i + 16) % FMT_OCC_INTERVAL != 0 && j > 8) // 序列最后的部分不用保存 - { for (m = 0; m < 4; ++m) buf[k++] = mc[m]; - } #endif + + buf[k++] = pre_and_bwt_seq_2; + +#ifdef FMT_MID_INTERVAL + s1 = pre_and_bwt_seq_2; + for (m = 0; m < 4; ++m) + mc[m] += cnt_table[m][s1 & 0xff] + cnt_table[m][s1 >> 8 & 0xff] + cnt_table[m][s1 >> 16 & 0xff] + cnt_table[m][s1 >> 24]; + + if ((i + 16) % FMT_OCC_INTERVAL != 0 && j == 16 && ((i + 16) & FMT_MID_INTV_MASK) == 0) + for (m = 0; m < 4; ++m) + buf[k++] = mc[m]; +#endif + } } } // the last element @@ -354,6 +347,7 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt) k += 16; // cout << "n occ: " << n_occ << endl; cout << "size: " << k << '\t' << fmt->bwt_size << endl; + // exit(0); xassert(k == fmt->bwt_size, "inconsistent bwt_size"); // update fmt fmt->bwt = buf; @@ -386,7 +380,7 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4] cnt[3] = q[b2]; // b2的occ p += 20; -#ifdef FMT_MID_INTERVAL +#ifdef FMT_MID_INTERVAL // 加入了middle checkpoint // 使用mid interval信息 int mk = k % FMT_OCC_INTERVAL; int n_mintv = mk >> FMT_MID_INTV_SHIFT; @@ -405,29 +399,34 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4] x = 0; p += 4; } -#if FMT_MID_INTERVAL == 16 +#if FMT_MID_INTERVAL == 16 // middle check point interval等于16时候,只需要判断一下是不是要计算两个uint32表示的碱基序列 if ((mk & FMT_MID_INTV_MASK) >> 3) { x += __fmt_occ_e2_aux2(fmt, ti, *p); ++p; } #elif FMT_MID_INTERVAL > 16 // 该地址是bwt和pre_bwt字符串数据的首地址 - uint32_t *end = p + ((k >> 3) - ((k & ~FMT_OCC_INTV_MASK) >> 3)); + uint32_t *end = p + ((k >> 3) - ((k & ~FMT_MID_INTV_MASK) >> 3)); for (; p < end; ++p) { x += __fmt_occ_e2_aux2(fmt, ti, *p); } - tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); - x += __fmt_occ_e2_aux2(fmt, ti, tmp); #endif -#else +#else // 没有加入middle check point interval +#if FMT_OCC_INTERVAL > 16 uint32_t *end = p + ((k >> 3) - ((k & ~FMT_OCC_INTV_MASK) >> 3)); + // p = end - (end - p) / 4; for (; p < end; ++p) { x += __fmt_occ_e2_aux2(fmt, ti, *p); } - tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); - x += __fmt_occ_e2_aux2(fmt, ti, tmp); +#else // FMT_OCC_INTERVAL等于16的时候,只需要判断一次就可以 + if ((k & FMT_OCC_INTV_MASK) >> 3) + { + x += __fmt_occ_e2_aux2(fmt, ti, *p); + ++p; + } +#endif #endif tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); x += __fmt_occ_e2_aux2(fmt, ti, tmp); @@ -452,119 +451,13 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4] cnt[3] += x >> 24 & 0xff; } -void calc_intv_occ(const FMTIndex *fmt, bwtint_t k, bwtint_t str_line, bwtint_t cp_line, - uint32_t *p, int b1, int b2, bwtint_t cnt[4]) -{ - - uint32_t x = 0; - uint32_t *q, tmp; - int i, ti = b1 << 2 | b2; -#ifdef FMT_MID_INTERVAL - // 使用mid interval信息 - int mk = k % FMT_OCC_INTERVAL; - int n_mintv = mk >> FMT_MID_INTV_SHIFT; - if (n_mintv > 0) // 至少超过了第一个mid interval - { - 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之和 - cnt[0] += __fmt_mid_sum(x); - x = *q; - cnt[1] += __fmt_mid_sum(x); // b1的occ - for (i = 3; i > b2; --i) - cnt[2] += x >> (i << 3) & 0xff; // 大于b2的occ之和 - cnt[3] += x >> (b2 << 3) & 0xff; // b2的occ - x = 0; - p += 4; - } -#if FMT_MID_INTERVAL == 16 - if ((mk & FMT_MID_INTV_MASK) >> 3) - { - x += __fmt_occ_e2_aux2(fmt, ti, *p); - ++p; - } -#elif FMT_MID_INTERVAL > 16 // 该地址是bwt和pre_bwt字符串数据的首地址 - uint32_t *end = p + ((k >> 3) - ((k & ~FMT_OCC_INTV_MASK) >> 3)); - for (; p < end; ++p) - { - x += __fmt_occ_e2_aux2(fmt, ti, *p); - } - tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); - x += __fmt_occ_e2_aux2(fmt, ti, tmp); -#endif -#else - uint32_t *end = p + ((k >> 3) - ((k & ~FMT_OCC_INTV_MASK) >> 3)); - for (; p < end; ++p) - { - x += __fmt_occ_e2_aux2(fmt, ti, *p); - } - tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); - x += __fmt_occ_e2_aux2(fmt, ti, tmp); -#endif - 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 && cp_line < fmt->sec_primary && str_line >= fmt->sec_primary) - { - if (b2 < fmt->last_base) - cnt[2] -= 1; - else if (b2 == fmt->last_base) - cnt[3] -= 1; - } - 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 tk[4], bwtint_t tl[4]) -{ - bwtint_t _k, _l; - _k = k - (k >= fmt->primary); - _l = l - (l >= fmt->primary); - if (_l >> FMT_OCC_INTV_SHIFT != _k >> FMT_OCC_INTV_SHIFT || k == (bwtint_t)(-1) || l == (bwtint_t)(-1)) - { - fmt_e2_occ(fmt, k, b1, b2, tk); - fmt_e2_occ(fmt, l, b1, b2, tl); - } - else - { - uint32_t *p, *q; - bwtint_t cp_line = _k >> FMT_OCC_INTV_SHIFT << FMT_OCC_INTV_SHIFT; - int i; - tk[0] = 0; - tk[2] = 0; - p = fmt_occ_intv(fmt, _k); - for (i = b1 + 1; i < 4; ++i) - tk[0] += p[i]; // 大于b1的碱基的occ之和 - tk[1] = p[b1]; // b1的occ - q = p + 4 + b1 * 4; - for (i = b2 + 1; i < 4; ++i) - tk[2] += q[i]; // 大于b2的occ之和 - tk[3] = q[b2]; // b2的occ - p += 20; - memcpy(tl, tk, 32); - calc_intv_occ(fmt, _k, k, cp_line, p, b1, b2, tk); - calc_intv_occ(fmt, _l, l, cp_line, p, b1, b2, tl); - } -} - // 扩展两个碱基 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_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_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); @@ -587,9 +480,8 @@ void fmt_extend1(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t *ok, int is_back, 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_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); // 这里是反向扩展 @@ -655,9 +547,9 @@ uint32_t str2bit(string &str) int main_fmtidx(int argc, char **argv) { // string seq("ACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTA"); - string seq("ACCCT"); - string rseq = calc_reverse_seq(seq); - seq = seq + rseq; + //string seq("ACCCT"); + //string rseq = calc_reverse_seq(seq); + //seq = seq + rseq; //create_bwt_mtx(seq); //cout << seq << endl; @@ -682,92 +574,93 @@ int main_fmtidx(int argc, char **argv) // } // return 0; + string prefix = argv[1]; + string bwt_str = prefix + ".bwt.str"; + string bwt_idx, fmt_idx; + bwt_t *bwt; + FMTIndex *fmt; + ostringstream oss_bwt, oss_fmt; + oss_bwt << '.' << OCC_INTERVAL; + oss_fmt << '.' << FMT_OCC_INTERVAL; +#ifdef FMT_MID_INTERVAL + oss_fmt << '.' << FMT_MID_INTERVAL; +#endif - // 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); + bwt_idx = prefix + oss_bwt.str() + ".bwt"; + fmt_idx = prefix + oss_fmt.str() + ".fmt"; - //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]) + ".64.fmt"; - // string fmt_idx = string(argv[1]) + ".32.fmt"; + cout << bwt_idx << endl; + cout << fmt_idx << endl; - t11 = realtime(); - bwt_t *bwt = restore_bwt(bwt_idx.c_str()); - t11 = realtime() - t11; - 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); - t10 = realtime() - t10; - cout << "[time gen fmt:] " << t10 << "s" << endl; +// 生成或读取bwt索引文件 + double time_read_bwt = realtime(); +#if 0 + bwt = restore_bwt(bwt_str.c_str()); + create_interval_occ_bwt(bwt); + dump_bwt(bwt_idx.c_str(), bwt); +#else + bwt = restore_bwt(bwt_idx.c_str()); +#endif + time_read_bwt = realtime() - time_read_bwt; + cout << "[time gen/read bwt:] " << time_read_bwt << "s" << endl; - vector seed_arr(10000000); - seed_arr[0] = "GCGATACTAAGA"; - seed_arr[0] = "GAGAGCTGTTCCCGTGTTTTCCATGGTTT"; +// 生成或读取fmt索引文件 + double time_read_fmt = realtime(); +#if 0 + fmt = create_fmt_from_bwt(bwt); + dump_fmt(fmt_idx.c_str(), fmt); +#else + fmt = restore_fmt(fmt_idx.c_str()); +#endif + time_read_fmt = realtime() - time_read_fmt; + cout << "[time gen/read fmt:] " << time_read_fmt << "s" << endl; + +// 生成随机序列进行测试 + int seq_size = 10000000; + // int seq_size = 1; + int seq_len = 19; + vector seq_arr(seq_size); + seq_arr[0] = "GCGATACTAAGA"; + // seq_arr[0] = "GAGAGCTGTTCCCGTGTTTTCCATGGTTT"; +#if 1 srand(time(NULL)); + double time_gen_seq = realtime(); + for (int i = 0; i < (int)seq_arr.size(); ++i) + seq_arr[i] = generate_rand_seq(seq_len); + time_gen_seq = realtime() - time_gen_seq; + cout << "[time gen seq:] " << time_gen_seq << "s" << endl; +#endif - t1 = realtime(); - for (int i = 0; i < (int)seed_arr.size(); ++i) - seed_arr[i] = generate_rand_seq(29); - t1 = realtime() - t1; - cout << "[time gen seed:] " << t1 << "s" << endl; +// 对比bwt和fmt搜索的时间 +#if 1 + double time_bwt_search = realtime(); + for (int i = 0; i < (int)seq_arr.size(); ++i) + bwt_search(bwt, seq_arr[i]); + time_bwt_search = realtime() - time_bwt_search; + cout << "[time bwt search:] " << time_bwt_search << "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; + double time_fmt_search = realtime(); + for (int i = 0; i < (int)seq_arr.size(); ++i) + fmt_search(fmt, seq_arr[i]); + time_fmt_search = realtime() - time_fmt_search; + cout << "[time fmt search:] " << time_fmt_search << "s" << endl; +#endif - //return 0; - - t4 = realtime(); - for (int i = 0; i < (int)seed_arr.size(); ++i) +// 对比bwt和fmt搜索的结果 +#if 1 + double time_cmp = realtime(); + for (int i = 0; i < (int)seq_arr.size(); ++i) { - bwtintv_t p1 = - bwt_search2(bwt, seed_arr[i]); - bwtintv_t p2 = fmt_search(fmt, seed_arr[i]); + bwtintv_t p1 = bwt_search2(bwt, seq_arr[i]); + bwtintv_t p2 = fmt_search(fmt, seq_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 + cout << seq_arr[i] << endl << p1.x[0] << ' ' << p1.x[1] << ' ' << p1.x[2] << endl - << p2.x[0] << ' ' << p2.x[1] << ' ' << p2.x[2] << endl; + << p2.x[0] << ' ' << p2.x[1] << ' ' << p2.x[2] << endl; } - t4 = realtime() - t4; - cout << "[time bwt search 2:] " << t4 << "s" << endl; + time_cmp = realtime() - time_cmp; + cout << "[time compare bwt & fmt:] " << time_cmp << "s" << endl; +#endif - //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) - // { - // cout << i << '\t' << bwt_B0(bwt, i) << endl; - // } - // TGGGAT - // FMTIndex *fmt = create_fmt_from_bwt(bwt); - // dump_fmt("ref.fmt", fmt); - // FMTIndex *fmt = restore_fmt("tiny.fmt"); - - //fmt_search(fmt, s); - // cout << bwt->bwt_size << endl; - // cout << bwt->seq_len << endl; - //cout << "sec_: " << fmt->sec_bcp << '\t' << fmt->sec_primary << endl; - - //uint8_t b8 = 2 << 4 | 2; - //cout << "AGAG: " << fmt->cnt_table[2][b8] << endl; - //cout << (((b8 >> 6) == 0 && (b8 >> 4 & 3) == 2) + ((b8 >> 2 & 3) == 0 && (b8 & 3) == 2)) << endl; return 0; } \ No newline at end of file diff --git a/fmt_index.h b/fmt_index.h index d4d58ca..aca2472 100644 --- a/fmt_index.h +++ b/fmt_index.h @@ -5,18 +5,26 @@ #define FMT_OCC_INTERVAL (1LL << FMT_OCC_INTV_SHIFT) #define FMT_OCC_INTV_MASK (FMT_OCC_INTERVAL - 1) -#define FMT_MID_INTV_SHIFT 4 +#define FMT_MID_INTV_SHIFT 7 #define FMT_MID_INTERVAL (1LL << FMT_MID_INTV_SHIFT) #define FMT_MID_INTV_MASK (FMT_MID_INTERVAL - 1) +// #undef FMT_MID_INTERVAL + // 获取碱基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的整数倍) -#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 + +#if FMT_MID_INTERVAL == 8 #define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 144)) +#elif FMT_MID_INTERVAL == 16 +#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 80)) +#elif FMT_MID_INTERVAL == 32 +#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 48)) +#elif FMT_MID_INTERVAL == 64 +#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 32)) +#elif FMT_MID_INTERVAL == 128 +#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 24)) #else #define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 20)) #endif