From f56edd07dd2f0b15609694b4814bb65b3a517c6d Mon Sep 17 00:00:00 2001 From: Heng Li Date: Tue, 25 Oct 2011 12:31:36 -0400 Subject: [PATCH] forward-backward search seems working --- Makefile | 1 + bwt.c | 13 ++++++------- bwt.h | 2 +- fastmap.c | 6 +++++- 4 files changed, 13 insertions(+), 9 deletions(-) diff --git a/Makefile b/Makefile index 20350506..e4b073fe 100644 --- a/Makefile +++ b/Makefile @@ -34,6 +34,7 @@ bwt1away.o:bwt.h bwtaln.h bwt2fmv.o:bwt.h bntseq.o:bntseq.h bwtgap.o:bwtgap.h bwtaln.h bwt.h +fastmap:bwt.h bwtsw2_core.o:bwtsw2.h bwt.h bwt_lite.h stdaln.h bwtsw2_aux.o:bwtsw2.h bwt.h bwt_lite.h stdaln.h diff --git a/bwt.c b/bwt.c index 6f5eb673..5b139d27 100644 --- a/bwt.c +++ b/bwt.c @@ -249,13 +249,13 @@ void bwt_extend(const bwt_t *bwt, const bwtintv_t *ik, bwtintv_t ok[4], int is_b int i; bwt_2occ4(bwt, 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] = bwt->L2[i] + tk[i] + 1; - ok[i].x[2] = (tl[i] -= tk[i]); + ok[i].x[!is_back] = bwt->L2[i] + 1 + tk[i]; + ok[i].x[2] = tl[i] - tk[i]; } - ok[3].x[is_back] = ik->x[is_back]; - ok[2].x[is_back] = ok[3].x[is_back] + tl[3]; - ok[1].x[is_back] = ok[2].x[is_back] + tl[2]; - ok[0].x[is_back] = ok[1].x[is_back] + tl[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[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]; } static void bwt_reverse_intvs(bwtintv_v *p) @@ -298,7 +298,6 @@ int bwt_smem1(const bwt_t *bwt, int len, const uint8_t *q, int x, bwtintv_v *mem ret = curr->a[0].info; // this will be the returned value swap = curr; curr = prev; prev = swap; - if (x == 40) printf("[%lld,%lld,%lld]\n", prev->a[0].x[0], prev->a[0].x[1], prev->a[0].x[2]); for (i = x - 1; i >= -1; --i) { // backward search for MEMs if (q[i] > 3) break; c = i < 0? 0 : q[i]; diff --git a/bwt.h b/bwt.h index e635b13c..75dd21c6 100644 --- a/bwt.h +++ b/bwt.h @@ -81,7 +81,7 @@ typedef struct { size_t n, m; bwtintv_t *a; } bwtintv_v; (bwt)->L2[bwt_B0(bwt, k)] + bwt_occ(bwt, k, bwt_B0(bwt, k)) \ : (bwt)->L2[bwt_B0(bwt, (k)-1)] + bwt_occ(bwt, k, bwt_B0(bwt, (k)-1))) -#define bwt_set_intv(bwt, c, ik) ((ik).x[0] = (bwt)->L2[(int)(c)], (ik).x[2] = (bwt)->L2[(int)(c)+1] - (bwt)->L2[(int)(c)], (ik).x[1] = (bwt)->L2[3-(c)], (ik).info = 0) +#define bwt_set_intv(bwt, c, ik) ((ik).x[0] = (bwt)->L2[(int)(c)]+1, (ik).x[2] = (bwt)->L2[(int)(c)+1]-(bwt)->L2[(int)(c)], (ik).x[1] = (bwt)->L2[3-(c)]+1, (ik).info = 0) #ifdef __cplusplus extern "C" { diff --git a/fastmap.c b/fastmap.c index cf952b77..de5bc0d9 100644 --- a/fastmap.c +++ b/fastmap.c @@ -40,7 +40,9 @@ int main_fastmap(int argc, char *argv[]) for (i = 0; i < seq->seq.l; ++i) seq->seq.s[i] = nst_nt4_table[(int)seq->seq.s[i]]; { - int beg = 98; + int beg = 97; + uint8_t str[3]; + bwtint_t b, e; bwtintv_t ik, ok[4]; bwt_set_intv(bwt, seq->seq.s[seq->seq.l - 1], ik); for (i = seq->seq.l - 2; i >= beg; --i) { @@ -50,6 +52,7 @@ int main_fastmap(int argc, char *argv[]) if (ik.x[2] == 0) break; } printf("[%lld,%lld,%lld] @ %d\n", ik.x[0], ik.x[1], ik.x[2], i+1); + //str[0] = '\1'; str[1] = '\3'; bwt_match_exact(bwt, 2, str, &b, &e); printf("%lld,%lld,%lld\n", b, e, e-b+1); printf("======================== %lld, [%lld,%lld,%lld,%lld]\n", bwt->primary, bwt->L2[1], bwt->L2[2]-bwt->L2[1], bwt->L2[3]-bwt->L2[2], bwt->L2[4]-bwt->L2[3]); bwt_set_intv(bwt, seq->seq.s[beg], ik); for (i = beg + 1; i < seq->seq.l; ++i) { @@ -59,6 +62,7 @@ int main_fastmap(int argc, char *argv[]) if (ik.x[2] == 0) break; } printf("[%lld,%lld,%lld] @ %d\n", ik.x[0], ik.x[1], ik.x[2], i-1); + //str[0] = '\0'; str[1] = '\2'; bwt_match_exact(bwt, 2, str, &b, &e); printf("%lld,%lld,%lld\n", b, e, e-b+1); } bwt_smem(bwt, seq->seq.l, (uint8_t*)seq->seq.s, &mem, tvec); printf(">%s\t%ld\n", seq->name.s, mem.n);