r1215: scoring apparently works

This commit is contained in:
Heng Li
2024-10-12 22:29:27 -04:00
parent d1f4c8d232
commit 68179ed195
3 changed files with 8 additions and 9 deletions

View File

@@ -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) {

View File

@@ -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)

View File

@@ -5,7 +5,7 @@
#include <stdio.h>
#include <sys/types.h>
#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