diff --git a/index.c b/index.c index cd8b40e..aa8db8f 100644 --- a/index.c +++ b/index.c @@ -12,6 +12,7 @@ #include "bseq.h" #include "minimap.h" #include "mmpriv.h" +#include "ksw2.h" #include "kvec.h" #include "khash.h" @@ -65,6 +66,7 @@ void mm_idx_destroy(mm_idx_t *mi) kh_destroy(idx, (idxhash_t*)mi->B[i].h); } } + if (mi->spsc) free(mi->spsc); if (mi->I) { for (i = 0; i < mi->n_seq; ++i) free(mi->I[i].a); @@ -657,6 +659,10 @@ int mm_idx_alt_read(mm_idx_t *mi, const char *fn) return n_alt; } +/******************* + * Known junctions * + *******************/ + #define sort_key_bed(a) ((a).st) KRADIX_SORT_INIT(bed, mm_idx_intv1_t, sort_key_bed, 4) @@ -774,3 +780,113 @@ int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uin } return left; } + +/**************** + * splice score * + ****************/ + +typedef struct mm_idx_spsc_s { + uint32_t n, m; + uint64_t *a; // pos<<56 | score<<1 | acceptor +} mm_idx_spsc_t; + +int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc) +{ + gzFile fp; + kstring_t str = {0,0,0}; + kstream_t *ks; + int32_t dret, j; + int64_t n_read = 0; + + fp = fn && strcmp(fn, "-") != 0? gzopen(fn, "rb") : gzdopen(0, "rb"); + if (fp == 0) return -1; + if (max_sc > 63) max_sc = 63; + idx->spsc = Kcalloc(0, mm_idx_spsc_t, idx->n_seq * 2); + ks = ks_init(fp); + while (ks_getuntil(ks, KS_SEP_LINE, &str, &dret) >= 0) { + mm_idx_spsc_t *s; + char *p, *q, *name = 0; + int32_t i, type = -1, strand = 0, cid = -1, score = -1; + int64_t pos = -1; + for (i = 0, p = q = str.s;; ++p) { + if (*p == '\t' || *p == 0) { + int c = *p; + *p = 0; + if (i == 0) { + name = q; + } else if (i == 1) { + pos = atol(q); + } else if (i == 2) { + strand = *q == '+'? 1 : '-'? -1 : 0; + } else if (i == 3) { + type = *q == 'D'? 0 : *q == 'A'? 1 : -1; + } else if (i == 4) { + score = atoi(q); + break; + } + if (c == 0) break; + q = p + 1, ++i; + } + } + if (i < 4) continue; // not enough fields + if (score > max_sc) score = max_sc; + if (score < -max_sc) score = -max_sc; + cid = mm_idx_name2id(idx, name); + if (cid < 0 || type < 0 || strand == 0 || pos < 0) continue; // FIXME: give a warning! + s = &idx->spsc[cid << 1 | (strand > 0? 0 : 1)]; + Kgrow(0, uint64_t, s->a, s->n, s->m); + if (pos > 0 && pos < idx->seq[cid].len) { // ignore scores at the ends + s->a[s->n++] = (uint64_t)pos << 8 | (score + KSW_SPSC_OFFSET) << 1 | type; + ++n_read; + } + } + ks_destroy(ks); + gzclose(fp); + for (j = 0; j < idx->n_seq * 2; ++j) { + mm_idx_spsc_t *s = &idx->spsc[j]; + if (s->n > 0) + radix_sort_64(s->a, s->a + s->n); + } + if (mm_verbose >= 3) + fprintf(stderr, "[M::%s] read %ld splice scores\n", __func__, (long)n_read); + return 0; +} + +static int32_t mm_idx_find_intv(int32_t n, const uint64_t *a, int64_t x) +{ + int32_t s = 0, e = n; + if (n == 0) return -1; + if (x < a[0]>>8) return -1; + while (s < e) { + int32_t mid = s + (e - s) / 2; + if (x >= a[mid]>>8 && (mid + 1 >= n || x < a[mid+1]>>8)) return mid; + else if (x < a[mid]>>8) e = mid; + else s = mid + 1; + } + 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 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; + memset(sc, 0xff, en - st); + s = &db->spsc[cid << 1 | (!!rev)]; + if (s->n > 0) { + int32_t j, l, r; + l = mm_idx_find_intv(s->n, s->a, st); + r = mm_idx_find_intv(s->n, s->a, en); + for (j = l + 1; j <= r; ++j) { + int64_t x = (s->a[j]>>8) - st; + uint8_t score = s->a[j] & 0xff; + assert(x <= en - st); + if (x == en - st) continue; + if (sc[x] == 0xff || sc[x] < score) sc[x] = score; + } + } + return en - st; +} diff --git a/kalloc.h b/kalloc.h index 4378672..87bdf28 100644 --- a/kalloc.h +++ b/kalloc.h @@ -31,6 +31,14 @@ void km_stat_print(const void *km); #define Kcalloc(km, type, cnt) ((type*)kcalloc((km), (cnt), sizeof(type))) #define Krealloc(km, type, ptr, cnt) ((type*)krealloc((km), (ptr), (cnt) * sizeof(type))) +#define Kgrow(km, type, ptr, __i, __m) do { \ + if ((__i) >= (__m)) { \ + (__m) = (__i) + 1; \ + (__m) += ((__m)>>1) + 16; \ + (ptr) = Krealloc(km, type, ptr, (__m)); \ + } \ + } while (0) + #define Kexpand(km, type, a, m) do { \ (m) = (m) >= 4? (m) + ((m)>>1) : 16; \ (a) = Krealloc(km, type, (a), (m)); \ diff --git a/ksw2.h b/ksw2.h index 1f94c6f..296030b 100644 --- a/ksw2.h +++ b/ksw2.h @@ -24,6 +24,8 @@ #define KSW_CIGAR_DEL 2 #define KSW_CIGAR_N_SKIP 3 +#define KSW_SPSC_OFFSET 64 + #ifdef __cplusplus extern "C" { #endif diff --git a/main.c b/main.c index 2fbfce1..004d0ba 100644 --- a/main.c +++ b/main.c @@ -79,6 +79,7 @@ static ko_longopt_t long_options[] = { { "secondary-seq", ko_no_argument, 354 }, { "ds", ko_no_argument, 355 }, { "rmq-inner", ko_required_argument, 356 }, + { "spsc", ko_required_argument, 357 }, { "dbg-seed-occ", ko_no_argument, 501 }, { "help", ko_no_argument, 'h' }, { "max-intron-len", ko_required_argument, 'G' }, @@ -128,7 +129,7 @@ int main(int argc, char *argv[]) mm_mapopt_t opt; mm_idxopt_t ipt; int i, c, n_threads = 3, n_parts, old_best_n = -1; - char *fnw = 0, *rg = 0, *junc_bed = 0, *s, *alt_list = 0; + char *fnw = 0, *rg = 0, *junc_bed = 0, *fn_spsc = 0, *s, *alt_list = 0; FILE *fp_help = stderr; mm_idx_reader_t *idx_rdr; mm_idx_t *mi; @@ -248,6 +249,7 @@ int main(int argc, char *argv[]) else if (c == 354) opt.flag |= MM_F_SECONDARY_SEQ; // --secondary-seq else if (c == 355) opt.flag |= MM_F_OUT_DS; // --ds else if (c == 356) opt.rmq_inner_dist = mm_parse_num(o.arg); // --rmq-inner + else if (c == 357) fn_spsc = o.arg; // --spsc else if (c == 501) mm_dbg_flag |= MM_DBG_SEED_FREQ; // --dbg-seed-occ else if (c == 330) { fprintf(stderr, "[WARNING] \033[1;31m --lj-min-ratio has been deprecated.\033[0m\n"); @@ -432,6 +434,7 @@ int main(int argc, char *argv[]) if (argc != o.ind + 1) mm_mapopt_update(&opt, mi); if (mm_verbose >= 3) mm_idx_stat(mi); if (junc_bed) mm_idx_bed_read(mi, junc_bed, 1); + if (fn_spsc) mm_idx_spsc_read(mi, fn_spsc, mm_max_spsc_bonus(&opt)); if (alt_list) mm_idx_alt_read(mi, alt_list); if (argc - (o.ind + 1) == 0) { mm_idx_destroy(mi); diff --git a/minimap.h b/minimap.h index 015f4a8..aec09e6 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#define MM_VERSION "2.28-r1209" +#define MM_VERSION "2.28-r1211-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 @@ -91,6 +91,7 @@ typedef struct { uint32_t *S; // 4-bit packed sequence struct mm_idx_bucket_s *B; // index (hidden) struct mm_idx_intv_s *I; // intervals (hidden) + struct mm_idx_spsc_s *spsc;// splice score (hidden) void *km, *h; } mm_idx_t; @@ -411,6 +412,10 @@ int mm_idx_alt_read(mm_idx_t *mi, const char *fn); int mm_idx_bed_read(mm_idx_t *mi, const char *fn, int read_junc); int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uint8_t *s); +int mm_max_spsc_bonus(const mm_mapopt_t *mo); +int32_t mm_idx_spsc_read(mm_idx_t *idx, const char *fn, int32_t max_sc); +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); + // deprecated APIs for backward compatibility void mm_mapopt_init(mm_mapopt_t *opt); mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads); diff --git a/options.c b/options.c index c0b67e6..1689471 100644 --- a/options.c +++ b/options.c @@ -179,6 +179,13 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo) return 0; } +int mm_max_spsc_bonus(const mm_mapopt_t *mo) +{ + int max_sc = (mo->q2 + 1) / 2 - 1; + max_sc = max_sc > mo->q2 - mo->q? max_sc : mo->q2 - mo->q; + return max_sc; +} + int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo) { if (mo->bw > mo->bw_long) {