zzh
/
bwa_perf
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

解决了各种bug，现在结果是准确的，比原始的扩展快了30%，比64间隔的扩展快了20%

This commit is contained in:
zzh 2024-02-01 13:04:57 +08:00
parent a500d76eaa
commit e8ceb3ff58
4 changed files with 185 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
bwt.cpp
View File

@ -155,7 +155,7 @@ bwtint_t bwt_occ(const bwt_t *bwt, bwtint_t k, uint8_t c)
// 统计k行（bwt mtx行）之前4种碱基累积数量，这里的k是bwt矩阵里的行，比bwt字符串多1 // 统计k行（bwt mtx行）之前4种碱基累积数量，这里的k是bwt矩阵里的行，比bwt字符串多1
void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4]) void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4])
{ {
bwtint_t x; bwtint_t x = 0;
uint32_t *p, tmp, *end; uint32_t *p, tmp, *end;
// bwtint_t bwt_k_base_line = k >> OCC_INTV_SHIFT << OCC_INTV_SHIFT; // bwtint_t bwt_k_base_line = k >> OCC_INTV_SHIFT << OCC_INTV_SHIFT;
if (k == (bwtint_t)(-1)) if (k == (bwtint_t)(-1))
@ -165,7 +165,7 @@ void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4])
} }
k -= (k >= bwt->primary); // because $ is not in bwt k -= (k >= bwt->primary); // because $ is not in bwt
p = bwt_occ_intv(bwt, k); 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[0] = p[0];
cnt[1] = p[1]; cnt[1] = p[1];
cnt[2] = p[2]; cnt[2] = p[2];
@ -173,14 +173,14 @@ void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4])
p += 4; // sizeof(bwtint_t) = 4*(sizeof(bwtint_t)/sizeof(uint32_t)) 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 + ((k >> 4) - ((k & ~OCC_INTV_MASK) >> 4)); // this is the end point of the following loop
//end = p + 1; //end = p + 1;
//cout << "k - kbase: " << k - bwt_k_base_line << endl; // cout << "k - kbase: " << k - bwt_k_base_line << endl;
for (x = 0; p < end; ++p) for (; p < end; ++p)
{ {
x += __occ_aux4(bwt, *p); x += __occ_aux4(bwt, *p);
//print_base_uint32(*p); //print_base_uint32(*p);
} }
tmp = *p & ~((1U << ((~k & 15) << 1)) - 1); tmp = *p & ~((1U << ((~k & 15) << 1)) - 1);
//print_base_uint32(tmp); // print_base_uint32(tmp);
x += __occ_aux4(bwt, tmp) - (~k & 15); // 这里多算了A，要减去 x += __occ_aux4(bwt, tmp) - (~k & 15); // 这里多算了A，要减去
cnt[0] += x & 0xff; cnt[0] += x & 0xff;
cnt[1] += x >> 8 & 0xff; cnt[1] += x >> 8 & 0xff;

2
bwt.h
View File

@ -7,7 +7,7 @@
using std::string; using std::string;
// occ间隔对应的右移位数，5表示间隔32行保存一次 // occ间隔对应的右移位数，5表示间隔32行保存一次
#define OCC_INTV_SHIFT 6 #define OCC_INTV_SHIFT 7
#define OCC_INTERVAL (1LL << OCC_INTV_SHIFT) #define OCC_INTERVAL (1LL << OCC_INTV_SHIFT)
#define OCC_INTV_MASK (OCC_INTERVAL - 1) #define OCC_INTV_MASK (OCC_INTERVAL - 1)

203
fmt_index.cpp
View File

@ -120,6 +120,23 @@ void fmt_gen_cnt_occ(FMTIndex *fmt)
} }
} }
// fmt-index的count table，4对应着bwt碱基的累积量，0,1,2,3分别对应着bwt是A,C,G,T，pre-bwt的累积量
void fmt_gen_cnt_table(uint32_t cnt_table[4][256])
{
int i, j, k;
uint32_t x = 0;
for (i = 0; i != 256; ++i) // 遍历单个字节的各种情况
{
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);
cnt_table[k][i] = x;
}
}
}
void dump_fmt(const char *fn, const FMTIndex *fmt) void dump_fmt(const char *fn, const FMTIndex *fmt)
{ {
FILE *fp; FILE *fp;
@ -169,6 +186,13 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
uint32_t c[4], c2[16]; /*c用来保存原来的bwt碱基串的累积值，c2用来保存pre-bwt和bwt碱基对的累计值，如AA..TT*/ uint32_t c[4], c2[16]; /*c用来保存原来的bwt碱基串的累积值，c2用来保存pre-bwt和bwt碱基对的累计值，如AA..TT*/
uint32_t *buf; /* 计算之后变成fmt结构中bwt部分 */ uint32_t *buf; /* 计算之后变成fmt结构中bwt部分 */
#ifdef FMT_MID_INTERVAL
// 加入中间的check point
uint32_t mc[4] = {0};
uint32_t cnt_table[4][256]; // 4对应原来的cnt_table，0,1,2,3,分别对应该碱基的扩展
fmt_gen_cnt_table(cnt_table);
#endif
fmt->seq_len = bwt->seq_len; // bwt碱基序列的长度，不包含$字符，也就是该长度比bwt matrix长度少1 fmt->seq_len = bwt->seq_len; // bwt碱基序列的长度，不包含$字符，也就是该长度比bwt matrix长度少1
for (i = 0; i < 5; ++i) for (i = 0; i < 5; ++i)
fmt->L2[i] = bwt->L2[i]; // 每个碱基的总累积值 fmt->L2[i] = bwt->L2[i]; // 每个碱基的总累积值
@ -177,6 +201,14 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
n_occ = (bwt->seq_len + FMT_OCC_INTERVAL - 1) / FMT_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 = (fmt->seq_len * 2 + 15) >> 4; // 要保存最后两列碱基
fmt->bwt_size += n_occ * 20; // A,C,G,T和AA,AC.....TG,TT共20个 fmt->bwt_size += n_occ * 20; // A,C,G,T和AA,AC.....TG,TT共20个
#ifdef FMT_MID_INTERVAL
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;
#endif
buf = (uint32_t *)calloc(fmt->bwt_size, 4); // 开辟计算fmt用到的缓存 buf = (uint32_t *)calloc(fmt->bwt_size, 4); // 开辟计算fmt用到的缓存
c[0] = c[1] = c[2] = c[3] = 0; c[0] = c[1] = c[2] = c[3] = 0;
@ -198,6 +230,7 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
k += 4; k += 4;
memcpy(buf + k, c2, sizeof(uint32_t) * 16); // pre-bwt:bwt碱基对的occ memcpy(buf + k, c2, sizeof(uint32_t) * 16); // pre-bwt:bwt碱基对的occ
k += 16; k += 16;
mc[0] = mc[1] = mc[2] = mc[3] = 0;
} }
// 每个32位整数保存8个倒数第二列碱基（pre-bwt）和8个倒数第一列(bwt)碱基 // 每个32位整数保存8个倒数第二列碱基（pre-bwt）和8个倒数第一列(bwt)碱基
if (i % 16 == 0) // 每个32位整数可以包含16个碱基，每次需要处理16个碱基，也就是间隔最小可以设置为16 if (i % 16 == 0) // 每个32位整数可以包含16个碱基，每次需要处理16个碱基，也就是间隔最小可以设置为16
@ -257,6 +290,7 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
} }
// 保存bwt和pre_bwt // 保存bwt和pre_bwt
uint32_t pre_and_bwt_seq = 0; uint32_t pre_and_bwt_seq = 0;
uint32_t pre_and_bwt_seq_2 = 0;
for (m = 0; m < 8; ++m) for (m = 0; m < 8; ++m)
{ {
int lshift_bit = 30 - 2 * m; int lshift_bit = 30 - 2 * m;
@ -265,14 +299,28 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
buf[k++] = pre_and_bwt_seq; buf[k++] = pre_and_bwt_seq;
if (j > 8) if (j > 8)
{ {
pre_and_bwt_seq = 0;
for (m = 8; m > 0; --m) for (m = 8; m > 0; --m)
{ {
int lshift_bit = 2 * m - 2; int lshift_bit = 2 * m - 2;
pre_and_bwt_seq |= (((pre_bwt_16_seq & (3 << lshift_bit)) << (m * 2)) | ((bwt_16_seq & (3 << lshift_bit)) << (m * 2 - 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; buf[k++] = pre_and_bwt_seq_2;
} }
#ifdef FMT_MID_INTERVAL
// 加入中间的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
} }
} }
// the last element // the last element
@ -284,15 +332,15 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
// cout << "size: " << k << '\t' << fmt->bwt_size << endl; // cout << "size: " << k << '\t' << fmt->bwt_size << endl;
xassert(k == fmt->bwt_size, "inconsistent bwt_size"); xassert(k == fmt->bwt_size, "inconsistent bwt_size");
// update fmt // update fmt
fmt->bwt = buf; fmt->bwt = buf;
return fmt; return fmt;
} }
// 扩展两个个碱基，计算bwt base为b的pre-bwt str中各个碱基的occ // 扩展两个个碱基，计算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]) void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4])
{ {
uint32_t x; uint32_t x = 0;
uint32_t *p, *q, tmp, *end; uint32_t *p, *q, tmp;
bwtint_t bwt_k_line = k, bwt_k_base_line = k >> FMT_OCC_INTV_SHIFT << FMT_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; int i, ti;
cnt[0] = 0; cnt[0] = 0;
@ -301,25 +349,80 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4]
if (k == (bwtint_t)(-1)) if (k == (bwtint_t)(-1))
{ {
p = fmt->bwt + 4 + b1 * 4; p = fmt->bwt + 4 + b1 * 4;
for (i = 3; i > b2; --i) cnt[2] += p[i]; for (i = b2 + 1; i < 4; ++i) cnt[2] += p[i];
cnt[3] = p[b2]; cnt[3] = p[b2];
return; return;
} }
ti = b1 << 2 | b2; ti = b1 << 2 | b2;
// _mm_prefetch((const char *)(&fmt->cnt_occ[ti]), _MM_HINT_T0);
k -= (k >= fmt->primary); // k由bwt矩阵对应的行转换成bwt字符串对应的行（去掉了$，所以大于$的行，都减掉1） k -= (k >= fmt->primary); // k由bwt矩阵对应的行转换成bwt字符串对应的行（去掉了$，所以大于$的行，都减掉1）
p = fmt_occ_intv(fmt, k); 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 = 3; i > b1; --i) cnt[0] += p[i]; for (i = b1 + 1; i < 4; ++i) cnt[0] += p[i]; // 大于b1的碱基的occ之和
cnt[1] = p[b1]; cnt[1] = p[b1]; // b1的occ
q = p + 4 + b1 * 4; q = p + 4 + b1 * 4;
for (i = 3; i > b2; --i) cnt[2] += q[i]; for (i = b2 + 1; i < 4 ; ++i) cnt[2] += q[i]; // 大于b2的occ之和
cnt[3] = q[b2]; cnt[3] = q[b2]; // b2的occ
p += 20; // 该地址是bwt和pre_bwt字符串数据的首地址 p += 20;
end = p + ((k >> 3) - ((k & ~FMT_OCC_INTV_MASK) >> 3)); // this is the end point of the following loop #ifdef FMT_MID_INTERVAL
//p = end - (end - p) / 4; // 使用mid interval信息
int mk = k % FMT_OCC_INTERVAL;
int n_mintv = mk >> FMT_MID_INTV_SHIFT;
//if (k == 3137454504)
//{
// for (i = 0; i < n_mintv; ++i)
// {
// q = p + i * 6;
// print_base_uint32(*q);
// print_base_uint32(*(q + 1));
// x = *(q + 2);
// cout << ((x) >> 24 & 0xff) << ' ' << ((x) >> 16 & 0xff) << ' ' << ((x) >> 8 & 0xff) << ' ' << ((x) & 0xff) << ' ' << __fmt_mid_sum(x) << endl;
// x = *(q + 3);
// cout << ((x) >> 24 & 0xff) << ' ' << ((x) >> 16 & 0xff) << ' ' << ((x) >> 8 & 0xff) << ' ' << ((x) & 0xff) << ' ' << __fmt_mid_sum(x) << endl;
// x = *(q + 4);
// cout << ((x) >> 24 & 0xff) << ' ' << ((x) >> 16 & 0xff) << ' ' << ((x) >> 8 & 0xff) << ' ' << ((x) & 0xff) << ' ' << __fmt_mid_sum(x) << endl;
// x = *(q + 5);
// cout << ((x) >> 24 & 0xff) << ' ' << ((x) >> 16 & 0xff) << ' ' << ((x) >> 8 & 0xff) << ' ' << ((x) & 0xff) << ' ' << __fmt_mid_sum(x) << endl;
// }
// x = 0;
//}
if (n_mintv > 0) // 至少超过了第一个mid interval
{
p += n_mintv * 6 - 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;
// if (k == 3137454504)
// {
// cout << ((x) >> 24 & 0xff) << '\t' << ((x) >> 16 & 0xff)
// << '\t' << ((x) >> 8 & 0xff) << '\t' << ((x) & 0xff) << endl;
// cout << __fmt_mid_sum(x) << endl;
// }
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;
}
// cout << "k: " << k << ' ' << cnt[1] << endl;
if ((mk & FMT_MID_INTV_MASK) >> 3)
{
x += __fmt_occ_e2_aux2(fmt, ti, *p);
// print_base_uint32(*p);
++p;
}
tmp = *p & ~((1U << ((~k & 7) << 2)) - 1);
// print_base_uint32(tmp);
x += __fmt_occ_e2_aux2(fmt, ti, tmp);
#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;
// cout << "k - kbase: " << k - bwt_k_base_line << endl; // cout << "k - kbase: " << k - bwt_k_base_line << endl;
for (x = 0; p < end; ++p) for (; p < end; ++p)
{ {
x += __fmt_occ_e2_aux2(fmt, ti, *p); x += __fmt_occ_e2_aux2(fmt, ti, *p);
// print_base_uint32(*p); // print_base_uint32(*p);
@ -327,14 +430,7 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4]
tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); tmp = *p & ~((1U << ((~k & 7) << 2)) - 1);
// print_base_uint32(tmp); // print_base_uint32(tmp);
x += __fmt_occ_e2_aux2(fmt, ti, tmp); x += __fmt_occ_e2_aux2(fmt, ti, tmp);
// end = p + (end - p) / 4; #endif
// 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) if (b1 == 0)
{ {
@ -345,10 +441,10 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4]
// 如果跨过了second_primary,那么可能需要减掉一次累积值 // 如果跨过了second_primary,那么可能需要减掉一次累积值
if (b1 == fmt->first_base && bwt_k_base_line < fmt->sec_primary && bwt_k_line >= fmt->sec_primary) if (b1 == fmt->first_base && bwt_k_base_line < fmt->sec_primary && bwt_k_line >= fmt->sec_primary)
{ {
if (b2 == fmt->last_base) if (b2 < fmt->last_base)
x -= 1 << 24;
else if (b2 < fmt->last_base)
x -= 1 << 16; x -= 1 << 16;
else if (b2 == fmt->last_base)
x -= 1 << 24;
} }
cnt[0] += x & 0xff; cnt[0] += x & 0xff;
cnt[1] += x >> 8 & 0xff; cnt[1] += x >> 8 & 0xff;
@ -405,13 +501,12 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q)
fmt_set_intv(fmt, bval(q[x]), ik); fmt_set_intv(fmt, bval(q[x]), ik);
ik.info = x + 1; 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) for (i = x + 1; i + 1 < qlen; i += 2)
{ {
if (bval(q[i]) < 4 && bval(q[i + 1]) < 4) if (bval(q[i]) < 4 && bval(q[i + 1]) < 4)
{ {
c1 = cbval(q[i]); c1 = cbval(q[i]);
c2 = cbval(q[i + 1]); c2 = cbval(q[i + 1]);
//double tm_t = realtime(); //double tm_t = realtime();
@ -440,6 +535,14 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q)
return ik; return ik;
} }
uint32_t str2bit(string &str)
{
uint32_t pac = 0;
for (int i = 0; i < 16; ++i)
pac = (pac << 2) | bval(str[i]);
return pac;
}
int main_fmtidx(int argc, char **argv) int main_fmtidx(int argc, char **argv)
{ {
// string seq("ACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTA"); // string seq("ACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTA");
@ -449,9 +552,31 @@ int main_fmtidx(int argc, char **argv)
//create_bwt_mtx(seq); //create_bwt_mtx(seq);
//cout << seq << endl; //cout << seq << endl;
//string bwt_str = string(argv[1]) + ".bwt.str"; // uint32_t cnt_table[4][256];
// string bwt_idx = string(argv[1]) + ".128.bwt"; // fmt_gen_cnt_table(cnt_table);
string bwt_idx = string(argv[1]) + ".64.bwt"; // string ss1("ATATATCAATTCTCTT");
// string ss2("ATTGCCTCTGAATTAC");
// uint32_t bs1 = str2bit(ss1), bs2 = str2bit(ss2);
// uint32_t s1 = bs1;
// uint32_t x = 0;
// for (int m = 0; m < 4; ++m)
// {
// x = cnt_table[m][s1 & 0xff] + cnt_table[m][s1 >> 8 & 0xff] + cnt_table[m][s1 >> 16 & 0xff] + cnt_table[m][s1 >> 24];
// cout << ((x) >> 24 & 0xff) << ' ' << ((x) >> 16 & 0xff) << ' ' << ((x) >> 8 & 0xff) << ' ' << ((x) & 0xff) << ' ' << __fmt_mid_sum(x) << endl;
// }
// x = 0;
// s1 = bs2;
// for (int m = 0; m < 4; ++m)
// {
// x = cnt_table[m][s1 & 0xff] + cnt_table[m][s1 >> 8 & 0xff] + cnt_table[m][s1 >> 16 & 0xff] + cnt_table[m][s1 >> 24];
// cout << ((x) >> 24 & 0xff) << ' ' << ((x) >> 16 & 0xff) << ' ' << ((x) >> 8 & 0xff) << ' ' << ((x) & 0xff) << ' ' << __fmt_mid_sum(x) << endl;
// }
// return 0;
// 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]) + ".16.bwt";
// string bwt_idx = string(argv[1]) + ".bwt"; // string bwt_idx = string(argv[1]) + ".bwt";
//bwt_t *bwt = restore_bwt(bwt_str.c_str()); //bwt_t *bwt = restore_bwt(bwt_str.c_str());
@ -459,17 +584,21 @@ int main_fmtidx(int argc, char **argv)
//dump_bwt(bwt_idx.c_str(), bwt); //dump_bwt(bwt_idx.c_str(), bwt);
//string fmt_idx = string(argv[1]) + ".fmt"; //string fmt_idx = string(argv[1]) + ".fmt";
// string fmt_idx = string(argv[1]) + ".256.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]) + ".64.fmt";
// string fmt_idx = string(argv[1]) + ".32.fmt"; // string fmt_idx = string(argv[1]) + ".32.fmt";
t11 = realtime();
bwt_t *bwt = restore_bwt(bwt_idx.c_str()); 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 = restore_fmt(fmt_idx.c_str());
// FMTIndex *fmt = create_fmt_from_bwt(bwt); // FMTIndex *fmt = create_fmt_from_bwt(bwt);
// dump_fmt(fmt_idx.c_str(), fmt); // dump_fmt(fmt_idx.c_str(), fmt);
t10 = realtime() - t10;
cout << "[time gen fmt:] " << t10 << "s" << endl;
vector<string> seed_arr(10000000); vector<string> seed_arr(10000000);
seed_arr[0] = "GCGATACTAAGA"; seed_arr[0] = "GCGATACTAAGA";
@ -492,6 +621,8 @@ int main_fmtidx(int argc, char **argv)
t3 = realtime() - t3; t3 = realtime() - t3;
cout << "[time fmt search:] " << t3 << "s" << endl; cout << "[time fmt search:] " << t3 << "s" << endl;
return 0;
t4 = realtime(); t4 = realtime();
for (int i = 0; i < (int)seed_arr.size(); ++i) for (int i = 0; i < (int)seed_arr.size(); ++i)
{ {

15
fmt_index.h
View File

@ -1,16 +1,23 @@
#ifndef FMT_INDEX_H_ #ifndef FMT_INDEX_H_
#define FMT_INDEX_H_ #define FMT_INDEX_H_
#define FMT_OCC_INTV_SHIFT 6 #define FMT_OCC_INTV_SHIFT 8
#define FMT_OCC_INTERVAL (1LL << FMT_OCC_INTV_SHIFT) #define FMT_OCC_INTERVAL (1LL << FMT_OCC_INTV_SHIFT)
#define FMT_OCC_INTV_MASK (FMT_OCC_INTERVAL - 1) #define FMT_OCC_INTV_MASK (FMT_OCC_INTERVAL - 1)
#define FMT_MID_INTV 16 #define FMT_MID_INTV_SHIFT 4
#define FMT_MID_INTERVAL (1LL << FMT_MID_INTV_SHIFT)
#define FMT_MID_INTV_MASK (FMT_MID_INTERVAL - 1)
// 获取碱基c（待查找序列的首个碱基）和对应的互补碱基对应的行，以及间隔 // 获取碱基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) #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的整数倍） // k行（bwt str行（不包含$））对应的check point occ数据起始地址（小于k且是OCC_INTERVAL的整数倍）
#ifdef FMT_MID_INTERVAL
// 包含mid interval
#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 80))
#else
#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 20)) #define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 20))
#endif
// 字节val中包含bwt base为b的pre-bwt中T G C A（按顺序保存在32位整数里（每个占8bit））的数量 // 字节val中包含bwt base为b的pre-bwt中T G C A（按顺序保存在32位整数里（每个占8bit））的数量
#define __fmt_occ_e2_aux4(fmt, b, val) \ #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]) ((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])
@ -18,6 +25,9 @@
#define __fmt_occ_e2_aux2(fmt, b, val) \ #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]) ((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])
#define __fmt_mid_sum(x) \
((x) >> 24 & 0xff) + ((x) >> 16 & 0xff) + ((x) >> 8 & 0xff) + ((x) & 0xff)
// fm-index, extend twice in one search step (one memory access) // fm-index, extend twice in one search step (one memory access)
struct FMTIndex struct FMTIndex
{ {
@ -28,7 +38,6 @@ struct FMTIndex
bwtint_t bwt_size; // size of bwt, about seq_len/4 bwtint_t bwt_size; // size of bwt, about seq_len/4
uint32_t *bwt; // BWT uint32_t *bwt; // BWT
// occurance array, separated to two parts // 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格式 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 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 first_base; // 序列的第一个碱基2bit的int类型，0,1,2,3