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

解决了创建fmt过程中的bug，把bwt中保存occ的类型从uint64改成了uint32

This commit is contained in:
zzh 2024-01-31 02:19:08 +08:00
parent 555e39df94
commit fe34be5d3a
6 changed files with 155 additions and 425 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.vscode/settings.json vendored
View File

@ -1,6 +1,7 @@
{ {
"files.associations": { "files.associations": {
"ostream": "cpp", "ostream": "cpp",
"iostream": "cpp" "iostream": "cpp",
"string": "cpp"
} }
} }

2
Makefile
View File

@ -1,7 +1,7 @@
CC= g++ CC= g++
NOWARN= -Wno-unused-result -Wno-unused-function NOWARN= -Wno-unused-result -Wno-unused-function
CFLAGS= #-g -Wall $(NOWARN) #-O2 CFLAGS= #-g -Wall $(NOWARN) #-O2
CPPFLAGS= -g -Wall $(NOWARN) -O3 CPPFLAGS= -g -Wall $(NOWARN) -O2
WRAP_MALLOC=-DUSE_MALLOC_WRAPPERS WRAP_MALLOC=-DUSE_MALLOC_WRAPPERS
SHOW_PERF= -DSHOW_PERF SHOW_PERF= -DSHOW_PERF
AR= ar AR= ar

95
bwt.cpp
View File

@ -10,6 +10,17 @@
using namespace std; 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) void dump_bwt(const char *fn, const bwt_t *bwt)
{ {
FILE *fp; FILE *fp;
@ -63,26 +74,29 @@ bwt_t *restore_bwt(const char *fn)
return bwt; return bwt;
} }
// 根据原始的字符串bwt创建interval-bwt // 根据原始的字符串bwt创建interval-bwt，用uint32_t来表示occ
void create_interval_occ_bwt(bwt_t *bwt) 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 *buf;
uint32_t c[4];
n_occ = (bwt->seq_len + OCC_INTERVAL - 1) / OCC_INTERVAL + 1; 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 buf = (uint32_t *)calloc(bwt->bwt_size, 4); // will be the new bwt
c[0] = c[1] = c[2] = c[3] = 0; c[0] = c[1] = c[2] = c[3] = 0;
// 计算occ，生成naive bwt // 计算occ，生成naive bwt
for (i = k = 0; i < bwt->seq_len; ++i) for (i = k = 0; i < bwt->seq_len; ++i)
{ {
// cout << i << '\t'; if (i % 1000000 == 0)
// cout << c[0] << ' ' << c[1] << ' ' << c[2] << ' ' << c[3] << endl; {
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) if (i % OCC_INTERVAL == 0)
{ {
memcpy(buf + k, c, sizeof(bwtint_t) * 4); memcpy(buf + k, c, sizeof(uint32_t) * 4);
k += sizeof(bwtint_t); // in fact: sizeof(bwtint_t)=4*(sizeof(bwtint_t)/4) 每个c包含多少个32位 k += 4; // 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;
} }
if (i % 16 == 0) if (i % 16 == 0)
buf[k++] = bwt->bwt[i / 16]; // 16 == sizeof(uint32_t)/2, 2个bit表示一个碱基 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 // the last element
// cout << c[0] << '\t' << c[1] << '\t' << c[2] << '\t' << c[3] << endl; // cout << c[0] << '\t' << c[1] << '\t' << c[2] << '\t' << c[3] << endl;
memcpy(buf + k, c, sizeof(bwtint_t) * 4); memcpy(buf + k, c, sizeof(uint32_t) * 4);
xassert(k + sizeof(bwtint_t) == bwt->bwt_size, "inconsistent bwt_size"); xassert(k + 4 == bwt->bwt_size, "inconsistent bwt_size");
// update bwt // update bwt
free(bwt->bwt); free(bwt->bwt);
bwt->bwt = buf; bwt->bwt = buf;
@ -143,6 +157,7 @@ void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4])
{ {
bwtint_t x; bwtint_t x;
uint32_t *p, tmp, *end; uint32_t *p, tmp, *end;
// bwtint_t bwt_k_base_line = k >> OCC_INTV_SHIFT << OCC_INTV_SHIFT;
if (k == (bwtint_t)(-1)) if (k == (bwtint_t)(-1))
{ {
memset(cnt, 0, 4 * sizeof(bwtint_t)); 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 k -= (k >= bwt->primary); // because $ is not in bwt
p = bwt_occ_intv(bwt, k); p = bwt_occ_intv(bwt, k);
// cout << "occ: " << p[0] << ' ' << p[1] << ' ' << p[2] << ' ' << p[3] << endl; //cout << "k: " << k << "; occ: " << p[0] << ' ' << p[1] << ' ' << p[2] << ' ' << p[3] << endl;
memcpy(cnt, p, 4 * sizeof(bwtint_t)); cnt[0] = p[0];
p += sizeof(bwtint_t); // sizeof(bwtint_t) = 4*(sizeof(bwtint_t)/sizeof(uint32_t)) 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 + ((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;
for (x = 0; p < end; ++p) 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); tmp = *p & ~((1U << ((~k & 15) << 1)) - 1);
//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;
@ -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) 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 tk[4], tl[4];
bwtint_t interval = 0; // oc1, oc2;
// int int1 = 0, int2 = 0;
int i; int i;
bwt_2occ4(bwt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl); // tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积 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) 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[!is_back] = bwt->L2[i] + 1 + tk[i]; // 起始行位置，互补链
ok[i].x[2] = tl[i] - 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，下边是正向扩展 // 因为计算的是互补碱基，所以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[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[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]; 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; //cout << (ik->x[!is_back] <= bwt->primary && ik->x[!is_back] + ik->x[2] - 1 >= bwt->primary) << endl;
*ik = ok[c1]; *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) for (i = 0; i != 4; ++i)
{ {
ok[i].x[!is_back] = bwt->L2[i] + 1 + tk[i]; // 起始行位置，互补链 ok[i].x[!is_back] = bwt->L2[i] + 1 + tk[i]; // 起始行位置，互补链
ok[i].x[2] = tl[i] - tk[i]; // 间隔 ok[i].x[2] = tl[i] - tk[i]; // 间隔
} }
// cout << "occ: " << ok[c2].x[2] << endl;
// 因为计算的是互补碱基，所以3对应着0,2对应1，下边是正向扩展 // 因为计算的是互补碱基，所以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[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[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[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]; 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); bwt_set_intv(bwt, bval(q[x]), ik);
ik.info = x + 1; 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) 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; // t10 += realtime() - tm_t;
ik = ok[c1]; ik = ok[c1];
ik.info = i + 1; 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; // cout << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
return ik; return ik;

8
bwt.h
View File

@ -7,7 +7,7 @@
using std::string; using std::string;
// occ间隔对应的右移位数，5表示间隔32行保存一次 // occ间隔对应的右移位数，5表示间隔32行保存一次
#define OCC_INTV_SHIFT 7 #define OCC_INTV_SHIFT 6
#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)
@ -28,9 +28,9 @@ using std::string;
// #define bwt_occ_intv(b, k) ((b)->bwt + ((k) >> 7 << 4)) // #define bwt_occ_intv(b, k) ((b)->bwt + ((k) >> 7 << 4))
///* For general OCC_INTERVAL, the following is correct: ///* For general OCC_INTERVAL, the following is correct:
// k行对应的bwt str碱基，对应的check point bwt str数据起始地址 // 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的整数倍） // 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））的数量 // 字节b中包含的T G C A（按顺序保存在32位整数里（每个占8bit））的数量
#define __occ_aux4(bwt, b) \ #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]) ((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_t *restore_bwt(const char *fn);
// 根据原始的字符串bwt创建interval-bwt // 根据原始的字符串bwt创建interval-bwt
void create_interval_occ_bwt(bwt_t *bwt); void create_interval_occ_bwt(bwt_t *bwt);
void dump_bwt(const char *fn, const bwt_t *bwt);
// 利用bwt搜索seed，完整搜索，只需要单向搜索 // 利用bwt搜索seed，完整搜索，只需要单向搜索
void bwt_search(bwt_t *bwt, const string &q); void bwt_search(bwt_t *bwt, const string &q);
// 每次扩展两步 // 每次扩展两步

460
fmt_index.cpp
View File

@ -10,6 +10,7 @@
#include <immintrin.h> #include <immintrin.h>
#include "util.h" #include "util.h"
#include "bwt.h"
#include "fmt_index.h" #include "fmt_index.h"
using namespace std; using namespace std;
@ -92,27 +93,6 @@ void create_bwt_mtx(string &seq)
// cout << endl; // 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 // 生成occ，每个字节对应一个pattern
void fmt_gen_cnt_occ(FMTIndex *fmt) 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); fread_fix(fp, fmt->bwt_size << 2, fmt->bwt);
fmt->seq_len = fmt->L2[4]; fmt->seq_len = fmt->L2[4];
fclose(fp); fclose(fp);
fmt_gen_cnt_table(fmt); // 字节所能表示的各种碱基组合中，各个碱基的累积数量 fmt_gen_cnt_occ(fmt); // 字节所能表示的各种碱基组合中，各个碱基的累积数量
fmt_gen_cnt_occ(fmt);
return fmt; return fmt;
} }
@ -184,7 +163,7 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
{ {
// FILE *fmt_out = fopen("fmt.txt", "w"); // FILE *fmt_out = fopen("fmt.txt", "w");
FMTIndex *fmt = (FMTIndex *)calloc(1, sizeof(FMTIndex)); 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]; 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*/ 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->L2[i] = bwt->L2[i]; // 每个碱基的总累积值
fmt->primary = bwt->primary; // $在末尾的行，在bwt matrix行中的排序位置 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 = (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个
buf = (uint32_t *)calloc(fmt->bwt_size, 4); // 开辟计算fmt用到的缓存 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) for (i = k = 0; i < bwt->seq_len; ++i)
{ {
// 记录occ // 记录occ
if (i % OCC_INTERVAL == 0) if (i % FMT_OCC_INTERVAL == 0)
{ {
memcpy(buf + k, c, sizeof(uint32_t) * 4); // bwt str中各个碱基的occ memcpy(buf + k, c, sizeof(uint32_t) * 4); // bwt str中各个碱基的occ
k += 4; k += 4;
@ -224,8 +203,8 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
if (i % 16 == 0) // 每个32位整数可以包含16个碱基，每次需要处理16个碱基，也就是间隔最小可以设置为16 if (i % 16 == 0) // 每个32位整数可以包含16个碱基，每次需要处理16个碱基，也就是间隔最小可以设置为16
{ {
uint32_t pre_bwt_16_seq = 0; // 16个pre-bwt碱基串 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信息 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 % OCC_INTERVAL) / 16; // 每OCC_INTERVAL个碱基共享同一个基准地址，每16个碱基共用一个uint32整型，因此需要偏移量来获取当前碱基串的首地址 int offset = (i % FMT_OCC_INTERVAL) / 16; // 每OCC_INTERVAL个碱基共享同一个基准地址，每16个碱基共用一个uint32整型，因此需要偏移量来获取当前碱基串的首地址
uint32_t bwt_16_seq = *(bwt_addr + offset); // 待处理的当前16个碱基串的首地址 uint32_t bwt_16_seq = *(bwt_addr + offset); // 待处理的当前16个碱基串的首地址
for (j = 0; j < 16; ++j) // 对于bwt碱基串，一个一个碱基分别处理 for (j = 0; j < 16; ++j) // 对于bwt碱基串，一个一个碱基分别处理
{ {
@ -309,300 +288,12 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
return fmt; 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 // 扩展两个个碱基，计算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;
uint32_t *p, *q, tmp, *end; 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; int i, ti;
cnt[0] = 0; cnt[0] = 0;
cnt[1] = 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)) if (k == (bwtint_t)(-1))
{ {
p = fmt->bwt + 4 + b1 * 4; p = fmt->bwt + 4 + b1 * 4;
for (i = 3; i > b2; --i) for (i = 3; i > b2; --i) cnt[2] += p[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); // _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 << "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 = 3; i > b1; --i) cnt[0] += p[i];
cnt[1] = p[b1]; cnt[1] = p[b1];
q = p + 4 + b1 * 4; q = p + 4 + b1 * 4;
for (i = 3; i > b2; --i) cnt[2] += q[i]; for (i = 3; i > b2; --i) cnt[2] += q[i];
cnt[3] = q[b2]; cnt[3] = q[b2];
p += 20; // 该地址是bwt和pre_bwt字符串数据的首地址 p += 20; // 该地址是bwt和pre_bwt字符串数据的首地址
end = p + ((k >> 3) - ((k & ~OCC_INTV_MASK) >> 3)); // this is the end point of the following loop end = p + ((k >> 3) - ((k & ~FMT_OCC_INTV_MASK) >> 3)); // this is the end point of the following loop
//if ((int)(end - p) > 8) // p = end - (end - p) / 4;
// cout << "interval: " << (int)(end - p) << endl; // cout << "k - kbase: " << k - bwt_k_base_line << endl;
// end = p + 2 + (end - p) / 6;
// end = p + 2;
for (x = 0; p < end; ++p) for (x = 0; p < end; ++p)
{ {
x += __fmt_occ_e2_aux2(fmt, ti, *p); x += __fmt_occ_e2_aux2(fmt, ti, *p);
// print_base_uint32(*p);
} }
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;
// 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)
{ {
x -= (~k & 7) << 8; x -= (~k & 7) << 8;
if (b2 == 0) if (b2 == 0)
x -= (~k & 7) << 24; x -= (~k & 7) << 24;
} }
// 如果跨过了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)
{ {
@ -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) else if (b2 < fmt->last_base)
x -= 1 << 16; x -= 1 << 16;
} }
cnt[0] += x & 0xff; cnt[0] += x & 0xff;
cnt[1] += x >> 8 & 0xff; cnt[1] += x >> 8 & 0xff;
cnt[2] += x >> 16 & 0xff; cnt[2] += x >> 16 & 0xff;
cnt[3] += x >> 24 & 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) 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; bwtint_t tk[4], tl[4], first_pos;
// tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积 // 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[!is_back] = fmt->L2[b2] + 1 + tk[3];
ok->x[2] = tl[3] - 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]; 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]; 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]; 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，完整搜索，只需要单向搜索 // 利用fmt搜索seed，完整搜索，只需要单向搜索
bwtintv_t fmt_search(FMTIndex *fmt, const string &q) bwtintv_t fmt_search(FMTIndex *fmt, const string &q)
{ {
bwtintv_t ik; bwtintv_t ik;
// bwtintv_t ok[4];
bwtintv_t ok; bwtintv_t ok;
int i, c1, c2, x = 0; int i, c1, c2, x = 0;
int qlen = (int)q.size(); int qlen = (int)q.size();
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 << 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)
{ {
@ -727,7 +416,6 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q)
c2 = cbval(q[i + 1]); c2 = cbval(q[i + 1]);
//double tm_t = realtime(); //double tm_t = realtime();
fmt_extend2(fmt, &ik, &ok, 0, c1, c2); fmt_extend2(fmt, &ik, &ok, 0, c1, c2);
//fmt_extend2_1(fmt, &ik, ok, 0, c1, c2);
ik = ok; ik = ok;
// t8 += realtime() - tm_t; // t8 += realtime() - tm_t;
ik.info = i + 1; ik.info = i + 1;
@ -741,11 +429,12 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q)
if (i < qlen && bval(q[i]) < 4) if (i < qlen && bval(q[i]) < 4)
{ // 最后一次扩展 { // 最后一次扩展
c1 = cbval(q[i]); c1 = cbval(q[i]);
//double tm_t = realtime(); // double tm_t = realtime();
// fmt_extend1(fmt, &ik, ok, 0, c1); fmt_extend1(fmt, &ik, &ok, 0, c1);
//t9 += realtime() - tm_t; ik = ok;
// t9 += realtime() - tm_t;
ik.info = i + 1; 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; // cout << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
return ik; return ik;
@ -760,14 +449,29 @@ int main_fmtidx(int argc, char **argv)
//create_bwt_mtx(seq); //create_bwt_mtx(seq);
//cout << seq << endl; //cout << seq << endl;
string bwt_idx = string(argv[1]) + ".bwt"; //string bwt_str = string(argv[1]) + ".bwt.str";
string fmt_idx = string(argv[1]) + ".fmt"; // 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()); bwt_t *bwt = restore_bwt(bwt_idx.c_str());
FMTIndex *fmt = restore_fmt(fmt_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<string> seed_arr(10000000); vector<string> seed_arr(10000000);
seed_arr[0] = "CCGTCATCATCCG"; seed_arr[0] = "GCGATACTAAGA";
srand(time(NULL)); srand(time(NULL));
t1 = realtime(); t1 = realtime();
@ -776,11 +480,16 @@ int main_fmtidx(int argc, char **argv)
t1 = realtime() - t1; t1 = realtime() - t1;
cout << "[time gen seed:] " << t1 << "s" << endl; cout << "[time gen seed:] " << t1 << "s" << endl;
// t2 = realtime(); t2 = realtime();
// for (int i = 0; i < seed_arr.size(); ++i) for (int i = 0; i < (int)seed_arr.size(); ++i)
// bwt_search(bwt, seed_arr[i]); bwt_search(bwt, seed_arr[i]);
// t2 = realtime() - t2; t2 = realtime() - t2;
// cout << "[time bwt search:] " << t2 << "s" << endl; 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(); t4 = realtime();
for (int i = 0; i < (int)seed_arr.size(); ++i) for (int i = 0; i < (int)seed_arr.size(); ++i)
@ -795,11 +504,6 @@ int main_fmtidx(int argc, char **argv)
} }
t4 = realtime() - t4; t4 = realtime() - t4;
cout << "[time bwt search 2:] " << t4 << "s" << endl; 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 << "bwt occ: " << t5 << "s; " << t7 << '\t' << t10 << endl;
//cout << "fmt occ: " << t6 << "s; " << t11 << '\t' << t8 << '\t' << t9 << endl; //cout << "fmt occ: " << t6 << "s; " << t11 << '\t' << t8 << '\t' << t9 << endl;

6
fmt_index.h
View File

@ -1,12 +1,14 @@
#ifndef FMT_INDEX_H_ #ifndef FMT_INDEX_H_
#define 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（待查找序列的首个碱基）和对应的互补碱基对应的行，以及间隔 // 获取碱基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的整数倍）
#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））的数量 // 字节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])