解决了各种bug，目前mid interval和普通的occ interval合并在一起处理了，通过编译指导语句，进行选择和调试
This commit is contained in:
zzh
parent
5ff04ec853
commit
45318fc95b
|
|
@ -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}", // 当前工作路径：当前文件所在的工作空间
|
||||
}
|
||||
|
|
|
|||
2
bwt.cpp
2
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)
|
||||
{
|
||||
|
|
|
|||
351
fmt_index.cpp
351
fmt_index.cpp
|
|
@ -8,6 +8,7 @@
|
|||
#include <algorithm>
|
||||
#include <string.h>
|
||||
#include <immintrin.h>
|
||||
#include <sstream>
|
||||
|
||||
#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<string> 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<string> 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;
|
||||
}
|
||||
18
fmt_index.h
18
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
|
||||
|
|
|
|||
Loading…
Reference in New Issue