diff --git a/bwa.c b/bwa.c index ed5f771..3a37a32 100644 --- a/bwa.c +++ b/bwa.c @@ -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; diff --git a/bwamem.c b/bwamem.c index 448bc4d..e48f89b 100644 --- a/bwamem.c +++ b/bwamem.c @@ -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) diff --git a/fastmap.c b/fastmap.c index d31340b..b1c36da 100644 --- a/fastmap.c +++ b/fastmap.c @@ -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); diff --git a/fmt_idx.c b/fmt_idx.c index 202caf2..739ee82 100644 --- a/fmt_idx.c +++ b/fmt_idx.c @@ -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,7 +782,14 @@ backward_search: for (j = 0; j < curr->n; ++j) { bwtintv_t *p = &curr->a[j]; // 前向扩展的种子 - for (i = x - 1; i > 0; i -= 2) + //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) diff --git a/fmt_idx.h b/fmt_idx.h index 730b355..7459a88 100644 --- a/fmt_idx.h +++ b/fmt_idx.h @@ -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]); diff --git a/utils.h b/utils.h index 881f912..8a1e919 100644 --- a/utils.h +++ b/utils.h @@ -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,