From 68179ed195c54c8340c207722360a0da2d790a69 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sat, 12 Oct 2024 22:29:27 -0400 Subject: [PATCH] r1215: scoring apparently works --- index.c | 8 +++----- ksw2_exts2_sse.c | 7 ++++--- minimap.h | 2 +- 3 files changed, 8 insertions(+), 9 deletions(-) diff --git a/index.c b/index.c index aa8db8f..2a85f46 100644 --- a/index.c +++ b/index.c @@ -800,6 +800,7 @@ int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc) fp = fn && strcmp(fn, "-") != 0? gzopen(fn, "rb") : gzdopen(0, "rb"); if (fp == 0) return -1; + if (idx->h == 0) mm_idx_index_name(idx); if (max_sc > 63) max_sc = 63; idx->spsc = Kcalloc(0, mm_idx_spsc_t, idx->n_seq * 2); ks = ks_init(fp); @@ -866,14 +867,11 @@ static int32_t mm_idx_find_intv(int32_t n, const uint64_t *a, int64_t x) assert(0); } -int64_t mm_idx_spsc_get(const mm_idx_t *db, int32_t cid, int64_t st0, int64_t en0, int32_t rev, uint8_t *sc) +int64_t mm_idx_spsc_get(const mm_idx_t *db, int32_t cid, int64_t st, int64_t en, int32_t rev, uint8_t *sc) { - int64_t st, en; const mm_idx_spsc_t *s; if (cid >= db->n_seq || cid < 0 || db->spsc == 0) return -1; - if (en0 < 0 || en0 > db->seq[cid].len) en0 = db->seq[cid].len; - if (!rev) st = st0, en = en0; - else st = db->seq[cid].len - en0, en = db->seq[cid].len - st0; + if (en < 0 || en > db->seq[cid].len) en = db->seq[cid].len; memset(sc, 0xff, en - st); s = &db->spsc[cid << 1 | (!!rev)]; if (s->n > 0) { diff --git a/ksw2_exts2_sse.c b/ksw2_exts2_sse.c index c16866d..7892317 100644 --- a/ksw2_exts2_sse.c +++ b/ksw2_exts2_sse.c @@ -192,11 +192,12 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin } if (junc && (flag & KSW_EZ_SPLICE_SCORE)) { // junc[] keeps the donor score - uint8_t donor_val = !(flag & KSW_EZ_REV_CIGAR)? 0 : 1; + uint8_t donor_val = !!(flag & KSW_EZ_SPLICE_FOR) == !(flag & KSW_EZ_REV_CIGAR)? 0 : 1; for (t = 0; t < tlen - 1; ++t) - ((int8_t*)donor)[t] += junc[t+1] == 0xff? -junc_pen : (junc[t+1]&1) == donor_val? (int8_t)(junc[t+1]>>1) - (int8_t)KSW_SPSC_OFFSET : 0; + ((int8_t*)donor)[t] += junc[t+1] == 0xff || (junc[t+1]&1) != donor_val? -junc_pen : (int8_t)(junc[t+1]>>1) - (int8_t)KSW_SPSC_OFFSET; for (t = 0; t < tlen - 1; ++t) - ((int8_t*)acceptor)[t] += junc[t+1] == 0xff? -junc_pen : (junc[t+1]&1) == !donor_val? (int8_t)(junc[t+1]>>1) - (int8_t)KSW_SPSC_OFFSET : 0; + ((int8_t*)acceptor)[t] += junc[t+1] == 0xff || (junc[t+1]&1) != !donor_val? -junc_pen : (int8_t)(junc[t+1]>>1) - (int8_t)KSW_SPSC_OFFSET; + //for (t = 0; t < tlen - 1; ++t) if (junc[t+1] != 0xff) fprintf(stderr, "Y2\t%d\t%d\t%c\t%d\n", ((int8_t*)donor)[t], ((int8_t*)acceptor)[t], "DA"[junc[t+1]&1], (int8_t)(junc[t+1]>>1) - (int8_t)KSW_SPSC_OFFSET); } else if (junc) { // junc[] keeps the splice sites if (!(flag & KSW_EZ_REV_CIGAR)) { for (t = 0; t < tlen - 1; ++t) diff --git a/minimap.h b/minimap.h index 8b4a7a7..6311379 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#define MM_VERSION "2.28-r1211-dirty" +#define MM_VERSION "2.28-r1215-dirty" #define MM_F_NO_DIAG (0x001LL) // no exact diagonal hit #define MM_F_NO_DUAL (0x002LL) // skip pairs where query name is lexicographically larger than target name