zzh
/
fast-bwa
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

kmer长度变为14，结果正确

This commit is contained in:
zzh 2024-02-12 20:54:57 +08:00
parent 463f7da138
commit d41b8da061
6 changed files with 245 additions and 121 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
bwa.c
View File

@ -303,7 +303,7 @@ FMTIndex *bwa_idx_load_fmt(const char *hint)
sprintf(suffix, ".256.%d.fmt", FMT_MID_INTERVAL);
strcpy(fmt_idx_fn, hint);
strcpy(fmt_idx_fn + l_hint, suffix);
sprintf(suffix, ".%d.kmer", KMER_LEN);
sprintf(suffix, ".14.xmer");
strcpy(kmer_idx_fn, hint);
strcpy(kmer_idx_fn + l_hint, suffix);
@ -313,12 +313,9 @@ FMTIndex *bwa_idx_load_fmt(const char *hint)
fprintf(stderr, "[E::%s] fail to locate the index files\n", __func__);
return 0;
}
fprintf(stderr, "zzh-fmt-1\n");
fmt = fmt_restore_fmt(fmt_idx_fn);
fprintf(stderr, "%s\n", kmer_idx_fn);
fmt->kmer_entry = fmt_restore_kmer_idx(kmer_idx_fn);
fprintf(stderr, "zzh-fmt-2\n");
fmt->kmer_hash = fmt_restore_kmer_idx(kmer_idx_fn);
strcpy(sa_fn, hint);
sprintf(suffix, ".33.%d.sa", SA_INTV);
@ -345,6 +342,7 @@ bwaidx_t *bwa_idx_load_from_disk(const char *hint, int which)
if (which & BWA_IDX_BWT) idx->bwt = bwa_idx_load_bwt(hint);
if (which & BWA_IDX_BWT) {
idx->fmt = bwa_idx_load_fmt(hint);
// 先和bwt共用sa
idx->fmt->sa = idx->bwt->sa;
idx->fmt->n_sa = idx->bwt->n_sa;
idx->fmt->sa_intv = idx->bwt->sa_intv;

10
bwamem.c
View File

@ -154,14 +154,15 @@ static void mem_collect_intv(const mem_opt_t *opt, const bwt_t *bwt, const FMTIn
// fprintf(stderr, "\n");
// first pass: find all SMEMs
//fprintf(fp1, "seq: %ld\n", dn++);
// dn ++;
fprintf(fp1, "seq: %ld\n", dn++);
//fprintf(stderr, "seq: %ld\n", dn++);
//dn ++;
while (x < len) {
if (seq[x] < 4) {
#ifdef SHOW_PERF
int64_t tmp_time = realtime_msec();
#endif
//x = bwt_smem1(bwt, len, seq, x, start_width, &a->mem1, a->tmpv);
// x = bwt_smem1(bwt, len, seq, x, start_width, &a->mem1, a->tmpv);
x = fmt_smem(fmt, len, seq, x, start_width, &a->mem1, a->tmpv);
#ifdef SHOW_PERF
tmp_time = realtime_msec() - tmp_time;
@ -169,7 +170,8 @@ static void mem_collect_intv(const mem_opt_t *opt, const bwt_t *bwt, const FMTIn
#endif
for (i = 0; i < a->mem1.n; ++i) {
bwtintv_t *p = &a->mem1.a[i];
//fprintf(fp1, "%ld %ld %ld %ld %d\n", p->x[0], p->x[1], p->x[2], p->info >> 32, (uint32_t)p->info);
fprintf(fp1, "%ld %ld %ld %ld %d\n", p->x[0], p->x[1], p->x[2], p->info >> 32, (uint32_t)p->info);
//fprintf(stderr, "%ld %ld %ld %ld %d\n", p->x[0], p->x[1], p->x[2], p->info >> 32, (uint32_t)p->info);
int slen = (uint32_t)p->info - (p->info >> 32); // seed length
max_seed_len = fmax(max_seed_len, slen);
if (slen >= opt->min_seed_len)

2
fastmap.c
View File

@ -48,6 +48,7 @@ int64_t time_ksw_extend2 = 0,
time_ksw_global2 = 0,
time_ksw_align2 = 0,
time_bwt_smem1a = 0,
time_fmt_smem_0 = 0,
time_bwt_extend = 0,
time_bwt_occ4 = 0,
time_bwt_sa = 0,
@ -431,6 +432,7 @@ int main_mem(int argc, char *argv[])
#ifdef SHOW_PERF
fprintf(stderr, "\n");
fprintf(stderr, "time_bwt_smem1a: %f s\n", time_bwt_smem1a / 1000.0 / opt->n_threads);
fprintf(stderr, "time_fmt_smem_0: %f s\n", time_fmt_smem_0 / 1000.0 / opt->n_threads);
fprintf(stderr, "time_bwt_extend: %f s\n", time_bwt_extend / 1000.0 / opt->n_threads);
fprintf(stderr, "time_bwt_sa: %f s\n", time_bwt_sa / 1000.0 / opt->n_threads);
fprintf(stderr, "time_ksw_extend2: %f s\n", time_ksw_extend2 / 1000.0 / opt->n_threads);

298
fmt_idx.c
View File

@ -117,24 +117,42 @@ FMTIndex *fmt_restore_fmt(const char *fn)
}
// 将kmer hash数据写入到文件
void dump_kmer_idx(const char *fn, const KmerEntry *kmer_entry)
void fmt_dump_kmer_idx(const char *fn, const KmerHash *kh)
{
FILE *fp;
fp = xopen(fn, "wb");
err_fwrite(kmer_entry, 1, KMER_ARR_SIZE * sizeof(KmerEntry), fp);
err_fwrite(kh->ke10, 1, (1 << (10 << 1)) * sizeof(KmerEntryArr), fp);
err_fwrite(kh->ke11, 1, (1 << (11 << 1)) * sizeof(KmerEntry), fp);
err_fwrite(kh->ke12, 1, (1 << (12 << 1)) * sizeof(KmerEntry), fp);
err_fwrite(kh->ke13, 1, (1 << (13 << 1)) * sizeof(KmerEntry), fp);
err_fwrite(kh->ke14, 1, (1 << (14 << 1)) * sizeof(KmerEntry), fp);
err_fflush(fp);
err_fclose(fp);
}
// 从文件中读取kmer hash信息
KmerEntry *fmt_restore_kmer_idx(const char *fn)
KmerHash fmt_restore_kmer_idx(const char *fn)
{
FILE *fp = xopen(fn, "rb");
uint32_t kmer_bytes = KMER_ARR_SIZE * sizeof(KmerEntry);
KmerEntry *kmer_entry = (KmerEntry *)malloc(kmer_bytes);
fread_fix(fp, kmer_bytes, kmer_entry);
KmerHash khash;
KmerHash *kh = &khash;
int len = 1 << (10 << 1);
kh->ke10 = (KmerEntryArr *)malloc(len * sizeof(KmerEntryArr));
fread_fix(fp, len * sizeof(KmerEntryArr), kh->ke10);
len = 1 << (11 << 1);
kh->ke11 = (KmerEntry *)malloc(len * sizeof(KmerEntry));
fread_fix(fp, len * sizeof(KmerEntry), kh->ke11);
len = 1 << (12 << 1);
kh->ke12 = (KmerEntry *)malloc(len * sizeof(KmerEntry));
fread_fix(fp, len * sizeof(KmerEntry), kh->ke12);
len = 1 << (13 << 1);
kh->ke13 = (KmerEntry *)malloc(len * sizeof(KmerEntry));
fread_fix(fp, len * sizeof(KmerEntry), kh->ke13);
len = 1 << (14 << 1);
kh->ke14 = (KmerEntry *)malloc(len * sizeof(KmerEntry));
fread_fix(fp, len * sizeof(KmerEntry), kh->ke14);
err_fclose(fp);
return kmer_entry;
return khash;
}
// 读取sa数据
@ -483,11 +501,11 @@ inline void fmt_extend1(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t *ok, int i
}
// 获取kmer的fmt匹配信息
inline void kmer_getval_at(KmerEntry *ke, bwtintv_t *ok, int pos)
inline void kmer_getval_at(uint8_t *mem_addr, bwtintv_t *ok, int pos)
{
bwtint_t x0, x1, x2;
int byte_idx = pos * 14;
uint8_t *arr = ke->intv_arr + byte_idx;
uint8_t *arr = mem_addr + byte_idx;
x0 = *arr;
x0 = (x0 << 32) | *((uint32_t *)(arr + 1));
arr += 5;
@ -501,10 +519,10 @@ inline void kmer_getval_at(KmerEntry *ke, bwtintv_t *ok, int pos)
}
// 设置kmer第pos个碱基对应的fmt匹配信息
inline void kmer_setval_at(KmerEntry *ke, bwtintv_t ik, int pos)
inline void kmer_setval_at(uint8_t *mem_addr, bwtintv_t ik, int pos)
{
int byte_idx = pos * 14;
uint8_t *arr = ke->intv_arr + byte_idx;
uint8_t *arr = mem_addr + byte_idx;
arr[0] = (uint8_t)(ik.x[0] >> 32);
*((uint32_t *)(arr + 1)) = (uint32_t)ik.x[0];
arr += 5;
@ -514,46 +532,29 @@ inline void kmer_setval_at(KmerEntry *ke, bwtintv_t ik, int pos)
*((uint32_t *)arr) = (uint32_t)ik.x[2];
}
// 查找并保存kmer中每扩展一个碱基对应的fmt位置信息
void fmt_search_store_kmer(FMTIndex *fmt, const char *q, int qlen, KmerEntry *ke)
// 获取kmer对应的fmt匹配信息, pos should be [0, 13]
inline void fmt_kmer_get(const FMTIndex *fmt, bwtintv_t *ok, uint32_t qbit, int pos)
{
bwtintv_t ik;
int i, c1, c2;
bwtint_t tk[4], tl[4];
fmt_set_intv(fmt, nst_nt4_table[(int)q[0]], ik);
kmer_setval_at(ke, ik, 0);
// 每次扩展两个碱基
for (i = 1; i + 1 < qlen; i += 2)
if (pos == 13)
{
// 默认kmer是由ACGT组成的，不含其他字符
c1 = 3 - nst_nt4_table[(int)q[i]];
c2 = 3 - nst_nt4_table[(int)q[i + 1]];
fmt_e2_occ(fmt, ik.x[1] - 1, c1, c2, tk);
fmt_e2_occ(fmt, ik.x[1] - 1 + ik.x[2], c1, c2, tl);
// 第一次扩展的结果
ik.x[0] = ik.x[0] + (ik.x[1] <= fmt->primary && ik.x[1] + ik.x[2] - 1 >= fmt->primary) + tl[0] - tk[0];
ik.x[1] = fmt->L2[c1] + 1 + tk[1];
ik.x[2] = tl[1] - tk[1];
kmer_setval_at(ke, ik, i);
// 第二次扩展的结果
ik.x[0] = ik.x[0] + (ik.x[1] <= fmt->primary && ik.x[1] + ik.x[2] - 1 >= fmt->primary) + tl[2] - tk[2];
ik.x[1] = fmt->L2[c2] + 1 + tk[3];
ik.x[2] = tl[3] - tk[3];
kmer_setval_at(ke, ik, i + 1);
kmer_getval_at(fmt->kmer_hash.ke14[qbit].intv_arr, ok, 0);
}
if (i < qlen)
{ // 最后一次扩展
c1 = 3 - nst_nt4_table[(int)q[i]];
c2 = 3;
fmt_e2_occ(fmt, ik.x[1] - 1, c1, c2, tk);
fmt_e2_occ(fmt, ik.x[1] - 1 + ik.x[2], c1, c2, tl);
// 第一次扩展的结果
ik.x[0] = ik.x[0] + (ik.x[1] <= fmt->primary && ik.x[1] + ik.x[2] - 1 >= fmt->primary) + tl[0] - tk[0];
ik.x[1] = fmt->L2[c1] + 1 + tk[1];
ik.x[2] = tl[1] - tk[1];
kmer_setval_at(ke, ik, i);
else if (pos == 12)
{
kmer_getval_at(fmt->kmer_hash.ke13[qbit >> 2].intv_arr, ok, 0);
}
else if (pos == 11)
{
kmer_getval_at(fmt->kmer_hash.ke12[qbit >> 4].intv_arr, ok, 0);
}
else if (pos == 10)
{
kmer_getval_at(fmt->kmer_hash.ke11[qbit >> 6].intv_arr, ok, 0);
}
else
{
kmer_getval_at(fmt->kmer_hash.ke10[qbit >> 8].intv_arr, ok, pos);
}
}
@ -585,67 +586,166 @@ static void fmt_reverse_intvs(bwtintv_v *p)
}
}
// 找smem（seed）
int fmt_smem(const FMTIndex *fmt, int len, const uint8_t *q, int x, int min_intv, bwtintv_v *mem, bwtintv_v *tmpvec[2])
// 创建正向的kmer
inline static uint32_t build_forward_kmer(const uint8_t *q, int qlen, int *base_consumed)
{
int i, j, ret, kmer_end = x + KMER_LEN;
bwtintv_t ik, ok1, ok2;
bwtintv_v a[2], *curr;
uint32_t qbit = 0;
//int only_forward = 0;
//if (x == 0 || q[x-1] > 3) only_forward = 1; // 只用向前扩展
mem->n = 0;
if (q[x] > 3) return x + 1;
if (min_intv < 1) min_intv = 1; // the interval size should be at least 1
kv_init(a[0]); kv_init(a[1]);
curr = tmpvec && tmpvec[1] ? tmpvec[1] : &a[1]; // use the temporary vector if provided
kmer_end = kmer_end > len ? len : kmer_end;
// 计算kmer hash key
for (i = x; i < kmer_end; ++i)
int i;
qlen = qlen < HASH_KMER_LEN ? qlen : HASH_KMER_LEN;
for (i = 0; i < qlen; ++i)
{
if (q[i] > 3) // 要考虑碱基是N
break;
qbit |= q[i] << ((KMER_LEN - 1 - (i - x)) << 1);
qbit |= q[i] << ((HASH_KMER_LEN - 1 - i) << 1);
}
KmerEntry *entry = &fmt->kmer_entry[qbit];
kmer_getval_at(entry, &ik, 0); // 初始碱基位置
*base_consumed = i;
return qbit;
}
// 创建f反向的kmer
inline static uint32_t build_backward_kmer(const uint8_t *q, int start_pos, int *base_consumed)
{
uint32_t qbit = 0;
int i, j, end_pos;
end_pos = start_pos - HASH_KMER_LEN;
end_pos = end_pos < 0 ? -1 : end_pos;
for (i = start_pos, j = 0; i > end_pos; --i, ++j)
{
if (q[i] > 3) // 要考虑碱基是N
break;
qbit |= q[i] << ((HASH_KMER_LEN - 1 - j) << 1);
}
*base_consumed = start_pos - i;
return qbit;
}
// 当x为0，或者q[x-1]为N时，只需要前向搜索即可
int fmt_smem_forward(const FMTIndex *fmt, int len, const uint8_t *q, int x, bwtintv_v *mem)
{
#ifdef SHOW_PERF
#if 1
int64_t tmp_time = realtime_msec();
#endif
#endif
int i, j = 1, ret, kmer_len;
const int min_intv = 1;
bwtintv_t ik, ok1, ok2;
uint32_t qbit = build_forward_kmer(&q[x], len - x, &kmer_len);
mem->n = 0;
fmt_kmer_get(fmt, &ik, qbit, kmer_len - j++);
while (ik.x[2] == 0) {
fmt_kmer_get(fmt, &ik, qbit, kmer_len - j++);
}
if (j != 2) { // kmer hash没有找到对应的匹配
ik.info = x + kmer_len - j + 2;
goto fmt_smem_forward_end;
}
ik.info = x + kmer_len;
// 继续向前扩展
for (i = x + kmer_len; i + 1 < len; i += 2) {
//if (ik.x[2] < 5)
// goto fmt_smem_forward_end;
if (q[i] < 4 && q[i + 1] < 4) // 两个都可以扩展
{
fmt_extend2(fmt, &ik, &ok1, &ok2, 0, 3 - q[i], 3 - q[i + 1]);
if (ok2.x[2] >= min_intv) { // 可以继续扩展
ik = ok2; ik.info = i + 2;
} else if (ok1.x[2] >= min_intv) { // 第二个间隔不够
ik = ok1; ik.info = i + 1;
goto fmt_smem_forward_end;
} else { // 两个间隔都不够
goto fmt_smem_forward_end;
}
}
else if (q[i] < 4) // q[i+1] >= 4，只能扩展一个
{
fmt_extend1(fmt, &ik, &ok1, 0, 3 - q[i]);
if (ok1.x[2] >= min_intv) {
ik = ok1; ik.info = i + 1;
}
goto fmt_smem_forward_end;
}
else { // q[i] >= 4
goto fmt_smem_forward_end;
}
}
if (i == len - 1) // 扩展到了最后一个碱基
{
if (q[i] < 4) {
fmt_extend1(fmt, &ik, &ok1, 0, 3 - q[i]);
if (ok1.x[2] >= min_intv) {
ik = ok1; ik.info = i + 1;
}
}
}
fmt_smem_forward_end:
ret = ik.info;
ik.info |= (uint64_t)x << 32;
kv_push(bwtintv_t, *mem, ik);
#ifdef SHOW_PERF
#if 1
tmp_time = realtime_msec() - tmp_time;
__sync_fetch_and_add(&time_fmt_smem_0, tmp_time);
#endif
#endif
return ret;
}
// 找smem（seed）
int fmt_smem(const FMTIndex *fmt, int len, const uint8_t *q, int x, int min_intv, bwtintv_v *mem, bwtintv_v *tmpvec[2])
{
int i, j, ret, kmer_len;
bwtintv_t ik, ok1, ok2;
// bwtintv_t tik, tok1, tok2;
bwtintv_v a[1], *curr;
uint32_t qbit = 0;
mem->n = 0;
if (q[x] > 3) return x + 1;
//if (x == 0 || q[x-1] > 3) return fmt_smem_forward(fmt, len, q, x, mem); // 只用向前扩展
if (min_intv < 1) min_intv = 1; // the interval size should be at least 1
kv_init(a[0]);
curr = tmpvec && tmpvec[0] ? tmpvec[0] : &a[0]; // use the temporary vector if provided
qbit = build_forward_kmer(&q[x], len - x, &kmer_len);
fmt_kmer_get(fmt, &ik, qbit, 0); // 初始碱基位置
ik.info = x + 1;
//fmt_set_intv(fmt, q[x], tik);
//tik.info = x + 1;
// check change of the interval size and whether the interval size is too small to be extended further
#define CHECK_INTV_CHANGE(iv, ov) \
if (ov.x[2] != iv.x[2]) { kv_push(bwtintv_t, *curr, iv); if (ov.x[2] < min_intv) break; }
#define CHECK_INTV_CHANGE(iv, ov, end_pos) \
if (ov.x[2] != iv.x[2]) { kv_push(bwtintv_t, *curr, iv); if (ov.x[2] < min_intv) break; } iv = ov; iv.info = end_pos
#define PUSH_VAL_AND_SKIP(iv) \
do { kv_push(bwtintv_t, *curr, iv); goto backward_search; } while(0)
// 处理kmer对应的匹配信息
for (j = 1, curr->n = 0; j < i - x; ++j)
{
kmer_getval_at(entry, &ok1, j);
CHECK_INTV_CHANGE(ik, ok1);
ik = ok1;
ik.info = x + j + 1;
for (j = 1, curr->n = 0; j < kmer_len; ++j) {
//fmt_extend1(fmt, &tik, &tok1, 0, 3 - q[x + j]);
//tik = tok1;
fmt_kmer_get(fmt, &ok1, qbit, j);
CHECK_INTV_CHANGE(ik, ok1, x + j + 1);
}
if (i != kmer_end) // 遇到了N
if (kmer_len != HASH_KMER_LEN) // 遇到了N或者到了序列最后
PUSH_VAL_AND_SKIP(ik);
// 扩展kmer之后的碱基
for (; i + 1 < len; i += 2)
for (i = (int)ik.info; i + 1 < len; i += 2)
{ // forward search
if (q[i] < 4 && q[i + 1] < 4)
{
fmt_extend2(fmt, &ik, &ok1, &ok2, 0, 3 - q[i], 3 - q[i + 1]);
CHECK_INTV_CHANGE(ik, ok1);
ok1.info = i + 1;
CHECK_INTV_CHANGE(ok1, ok2);
ik = ok2;
ik.info = i + 2;
CHECK_INTV_CHANGE(ik, ok1, i + 1);
CHECK_INTV_CHANGE(ik, ok2, i + 2);
} else if (q[i] < 4) // q[i+1] >= 4
{
fmt_extend1(fmt, &ik, &ok1, 0, 3 - q[i]);
CHECK_INTV_CHANGE(ik, ok1);
ik = ok1;
ik.info = i + 1;
CHECK_INTV_CHANGE(ik, ok1, i + 1);
PUSH_VAL_AND_SKIP(ik);
}
else // q[i] >= 4
@ -653,32 +753,23 @@ int fmt_smem(const FMTIndex *fmt, int len, const uint8_t *q, int x, int min_intv
PUSH_VAL_AND_SKIP(ik);
}
}
if (i == len - 1) // 扩展到了最后一个碱基
for (; i == len - 1; ++i) // 扩展到了最后一个碱基
{
if (q[i] < 4)
{
if (q[i] < 4) {
fmt_extend1(fmt, &ik, &ok1, 0, 3 - q[i]);
if (ok1.x[2] != ik.x[2]) {
kv_push(bwtintv_t, *curr, ik);
if (ok1.x[2] < min_intv)
goto backward_search;
}
ik = ok1;
ik.info = i + 1;
CHECK_INTV_CHANGE(ik, ok1, i + 1);
}
else
PUSH_VAL_AND_SKIP(ik);
++i;
}
if (i == len)
kv_push(bwtintv_t, *curr, ik); // push the last interval if we reach the end
backward_search:
fmt_reverse_intvs(curr); // s.t. smaller intervals (i.e. longer matches) visited first
ret = curr->a[0].info; // this will be the returned value，扩展到的最远的位置
// swap = curr;
// curr = prev;
// prev = swap;
// 按照种子进行遍历，反向扩展
#define CHECK_PUT_MEM(ok, pos, intv) \
if (ok.x[2] < min_intv) { \
@ -691,6 +782,13 @@ backward_search:
for (j = 0; j < curr->n; ++j)
{
bwtintv_t *p = &curr->a[j]; // 前向扩展的种子
//if (p->info - x < HASH_KMER_LEN) {
// // 创建反向kmer
// uint32_t qbit = build_backward_kmer(q, p->info - 1, &kmer_len);
// fmt_kmer_get(fmt, &ik, qbit, kmer_len - 1);
//}
// for (i = p->info - kmer_len; i > 0 i -= 2)
for (i = x - 1; i > 0; i -= 2)
{
if (q[i] < 4 && q[i - 1] < 4) // 两个都可以扩展
@ -704,7 +802,7 @@ backward_search:
}
else if (q[i] < 4) // 只能扩展一个
{
fmt_extend1(fmt, p, &ok1, 0, q[i]);
fmt_extend1(fmt, p, &ok1, 1, q[i]);
CHECK_PUT_MEM(ok1, i + 1, *p);
} else
{ // 不能扩展
@ -718,7 +816,7 @@ backward_search:
}
if (i == 0) { // 扩展到了第一个碱基
if (q[i] < 4) {
fmt_extend1(fmt, p, &ok1, 0, q[i]);
fmt_extend1(fmt, p, &ok1, 1, q[i]);
CHECK_PUT_MEM(ok1, i + 1, *p);
} else {
if (mem->n == 0 || (i + 1) < mem->a[mem->n - 1].info >> 32)

45
fmt_idx.h
View File

@ -51,16 +51,40 @@ Date : 2023/12/24
// sa存储的行间隔
#define SA_INTV 2
#define KMER_LEN 12
#define KMER_ARR_SIZE ((1 << (KMER_LEN << 1)))
#define HASH_KMER_LEN 14
#define BIT_KMER_LEN 17
// 用来保存kmer对应的fmt的位置信息
typedef struct
{
// 40+40+32 14个byte，这样好处理
uint8_t intv_arr[14 * KMER_LEN]; // 保存kmer中每扩展一个碱基对应的bwtintv_t数据
uint8_t intv_arr[14]; // 保存kmer中每扩展一个碱基对应的bwtintv_t数据
} KmerEntry;
typedef struct
{
uint8_t intv_arr[140]; // 保存长度为10的kmer，每个碱基对应的bwt匹配信息
} KmerEntryArr;
// 保存各个位置对应的bwt匹配信息
typedef struct
{
KmerEntryArr *ke10;
KmerEntry *ke11;
KmerEntry *ke12;
KmerEntry *ke13;
KmerEntry *ke14;
} KmerHash;
// 用来检测15,16,17这些长度的序列是否在bwt索引里有匹配
typedef struct
{
uint8_t *kb15;
uint8_t *kb16;
uint8_t *kb17;
} KmerBit;
// fm-index, extend twice in one search step (one memory access)
typedef struct
{
@ -76,7 +100,8 @@ typedef struct
uint8_t first_base; // 序列的第一个碱基2bit的int类型，0,1,2,3
uint8_t last_base; // dollar转换成的base
// 保存kmer对应的fmt位置信息
KmerEntry *kmer_entry;
KmerHash kmer_hash;
KmerBit kmer_bit; // 用来
// suffix array
int sa_intv;
bwtint_t n_sa;
@ -88,9 +113,9 @@ void dump_fmt(const char *fn, const FMTIndex *fmt);
// 从文件中读取fmt结构数据
FMTIndex *fmt_restore_fmt(const char *fn);
// 将kmer hash数据写入到文件
void dump_kmer_idx(const char *fn, const KmerEntry *kmer_entry);
void fmt_dump_kmer_idx(const char *fn, const KmerHash *kh);
// 从文件中读取kmer hash信息
KmerEntry *fmt_restore_kmer_idx(const char *fn);
KmerHash fmt_restore_kmer_idx(const char *fn);
// 读取sa数据
void fmt_restore_sa(const char *fn, FMTIndex *fmt);
// 根据interval-bwt创建fmt-index
@ -101,15 +126,13 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4]
void fmt_extend2(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t *ok1, bwtintv_t *ok2, int is_back, int b1, int b2);
// 扩展一个碱基
void fmt_extend1(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t *ok, int is_back, int b1);
// 查找并保存kmer中每扩展一个碱基对应的fmt位置信息
void fmt_search_store_kmer(FMTIndex *fmt, const char *q, int qlen, KmerEntry *ke);
// 生成所有KMER_LEN长度的序列，字符串表示
void gen_all_seq(char **seq_arr, int kmer_len);
// 设置kmer第pos个碱基对应的fmt匹配信息
void kmer_setval_at(KmerEntry *ke, bwtintv_t ik, int pos);
void kmer_setval_at(uint8_t *mem_addr, bwtintv_t ik, int pos);
// 获取kmer的fmt匹配信息
void kmer_getval_at(KmerEntry *ke, bwtintv_t *ok, int pos);
void kmer_getval_at(uint8_t *mem_addr, bwtintv_t *ok, int pos);
void fmt_kmer_get(const FMTIndex *fmt, bwtintv_t *ok, uint32_t qbit, int pos);
// 找smem（seed）
int fmt_smem(const FMTIndex *fmt, int len, const uint8_t *q, int x, int min_intv, bwtintv_v *mem, bwtintv_v *tmpvec[2]);

1
utils.h
View File

@ -37,6 +37,7 @@ extern int64_t time_ksw_extend2,
time_ksw_global2,
time_ksw_align2,
time_bwt_smem1a,
time_fmt_smem_0,
time_bwt_extend,
time_bwt_occ4,
time_bwt_sa,