diff --git a/index.c b/index.c index 3bd4d1f..f99ba9d 100644 --- a/index.c +++ b/index.c @@ -669,9 +669,9 @@ int mm_idx_alt_read(mm_idx_t *mi, const char *fn) return n_alt; } -/******************* - * Known junctions * - *******************/ +/*************** + * BED reading * + ***************/ #define sort_key_bed(a) ((a).st) KRADIX_SORT_INIT(bed, mm_idx_intv1_t, sort_key_bed, 4) @@ -679,10 +679,7 @@ KRADIX_SORT_INIT(bed, mm_idx_intv1_t, sort_key_bed, 4) #define sort_key_end(a) ((a).en) KRADIX_SORT_INIT(end, mm_idx_intv1_t, sort_key_end, 4) -#define sort_key_jj(a) ((a).off) -KRADIX_SORT_INIT(jj, mm_idx_jjump1_t, sort_key_jj, 4) - -static mm_idx_intv_t *mm_idx_bed_read_core(const mm_idx_t *mi, const char *fn, int read_junc, int is_pass1) +static mm_idx_intv_t *mm_idx_bed_read_core(const mm_idx_t *mi, const char *fn, int read_junc, int min_sc) { gzFile fp; kstream_t *ks; @@ -729,7 +726,7 @@ static mm_idx_intv_t *mm_idx_bed_read_core(const mm_idx_t *mi, const char *fn, i } } if (id < 0 || t.st < 0 || t.st >= t.en) continue; // contig ID not found, or other problems - if (is_pass1 && t.score < 5) continue; // for pass-1 BED, ignore junctions with weak signals; NB: paired with pass-1! + if (min_sc > 0 && t.score < min_sc) continue; r = &I[id]; if (i >= 11 && read_junc) { // BED12 int32_t st, sz, en; @@ -762,12 +759,12 @@ static mm_idx_intv_t *mm_idx_bed_read_core(const mm_idx_t *mi, const char *fn, i return I; } -static mm_idx_intv_t *mm_idx_bed_read_merged(const mm_idx_t *mi, const char *fn, int read_junc, int is_pass1) +static mm_idx_intv_t *mm_idx_bed_read_merge(const mm_idx_t *mi, const char *fn, int read_junc, int min_sc) { long n = 0, n0 = 0; int32_t i; mm_idx_intv_t *I; - I = mm_idx_bed_read_core(mi, fn, read_junc, is_pass1); + I = mm_idx_bed_read_core(mi, fn, read_junc, min_sc); if (I == 0) return 0; for (i = 0; i < mi->n_seq; ++i) { int32_t j, j0, k; @@ -796,46 +793,11 @@ static mm_idx_intv_t *mm_idx_bed_read_merged(const mm_idx_t *mi, const char *fn, return I; } -static mm_idx_jjump_t *mm_idx_bed2jjump(const mm_idx_t *mi, const mm_idx_intv_t *I) -{ - int32_t i; - mm_idx_jjump_t *J; - J = CALLOC(mm_idx_jjump_t, mi->n_seq); - for (i = 0; i < mi->n_seq; ++i) { - int32_t j, k; - const mm_idx_intv_t *intv = &I[i]; - mm_idx_jjump_t *jj = &J[i]; - jj->n = intv->n * 2; - jj->a = CALLOC(mm_idx_jjump1_t, jj->n); - for (j = k = 0; j < intv->n; ++j) { - jj->a[k].off = intv->a[j].st, jj->a[k].off2 = intv->a[j].en, jj->a[k].cnt = intv->a[j].cnt, jj->a[k].strand = intv->a[j].strand, ++k; - jj->a[k].off = intv->a[j].en, jj->a[k].off2 = intv->a[j].st, jj->a[k].cnt = intv->a[j].cnt, jj->a[k].strand = intv->a[j].strand, ++k; - } - radix_sort_jj(jj->a, jj->a + jj->n); - } - return J; -} - -int mm_idx_bed_read2(mm_idx_t *mi, const char *fn, int read_junc, int for_score, int for_jump) -{ - int32_t i; - mm_idx_intv_t *I; - if (mi->h == 0) mm_idx_index_name(mi); - I = mm_idx_bed_read_merged(mi, fn, read_junc, 0); - if (I == 0) return 0; - if (for_jump) - mi->J = mm_idx_bed2jjump(mi, I); - if (!for_score) { - for (i = 0; i < mi->n_seq; ++i) - free(I[i].a); - free(I); - } else mi->I = I; - return 0; -} - int mm_idx_bed_read(mm_idx_t *mi, const char *fn, int read_junc) { - return mm_idx_bed_read2(mi, fn, read_junc, 1, 0); + if (mi->h == 0) mm_idx_index_name(mi); + mi->I = mm_idx_bed_read_merge(mi, fn, read_junc, -1); + return 0; } int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uint8_t *s) @@ -863,6 +825,96 @@ int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uin return left; } +/********************************* + * Reading junctions for jumping * + *********************************/ + +#define sort_key_jj(a) ((a).off) +KRADIX_SORT_INIT(jj, mm_idx_jjump1_t, sort_key_jj, 4) + +#define sort_key_jj2(a) ((a).off2) +KRADIX_SORT_INIT(jj2, mm_idx_jjump1_t, sort_key_jj2, 4) + +static mm_idx_jjump_t *mm_idx_bed2jjump(const mm_idx_t *mi, const mm_idx_intv_t *I, uint16_t flag) +{ + int32_t i; + mm_idx_jjump_t *J; + J = CALLOC(mm_idx_jjump_t, mi->n_seq); + for (i = 0; i < mi->n_seq; ++i) { + int32_t j, k; + const mm_idx_intv_t *intv = &I[i]; + mm_idx_jjump_t *jj = &J[i]; + jj->n = intv->n * 2; + jj->a = CALLOC(mm_idx_jjump1_t, jj->n); + for (j = k = 0; j < intv->n; ++j) { + jj->a[k].off = intv->a[j].st, jj->a[k].off2 = intv->a[j].en, jj->a[k].cnt = intv->a[j].cnt, jj->a[k].strand = intv->a[j].strand, jj->a[k++].flag = flag; + jj->a[k].off = intv->a[j].en, jj->a[k].off2 = intv->a[j].st, jj->a[k].cnt = intv->a[j].cnt, jj->a[k].strand = intv->a[j].strand, jj->a[k++].flag = flag; + } + radix_sort_jj(jj->a, jj->a + jj->n); + } + return J; +} + +static mm_idx_jjump_t *mm_idx_jjump_merge(const mm_idx_t *mi, const mm_idx_jjump_t *J0, const mm_idx_jjump_t *J1) +{ + int32_t i; + mm_idx_jjump_t *J2; + J2 = CALLOC(mm_idx_jjump_t, mi->n_seq); + for (i = 0; i < mi->n_seq; ++i) { + int32_t j, j0, k; + const mm_idx_jjump_t *jj0 = &J0[i], *jj1 = &J1[i]; + mm_idx_jjump_t *jj2 = &J2[i]; + // merge jj0 and jj1 into jj2; faster with sorted merge but the performance difference should be negligible + jj2->n = jj0->n + jj1->n; + jj2->a = CALLOC(mm_idx_jjump1_t, jj2->n); + for (j = k = 0; j < jj0->n; ++j) jj2->a[k++] = jj0->a[j]; + for (j = k = 0; j < jj1->n; ++j) jj2->a[k++] = jj1->a[j]; + radix_sort_jj(jj2->a, jj2->a + jj2->n); // sort by a[].off + // sort by a[].off and then by a[].off2 such that they can be merged later + for (j0 = 0, j = 1; j <= jj2->n; ++j) { + if (j == jj2->n || jj2->a[j0].off != jj2->a[j].off) { + radix_sort_jj2(jj2->a + j0, jj2->a + j); + j0 = j; + } + } + // the actual merge + for (j0 = 0, j = 1, k = 0; j <= jj2->n; ++j) { + if (j == jj2->n || jj2->a[j0].off != jj2->a[j].off || jj2->a[j0].off2 != jj2->a[j].off2) { + int32_t t, cnt = 0; + uint16_t flag = 0; + for (t = j0; t < j; ++t) cnt += jj2->a[t].cnt, flag |= jj2->a[t].flag; + jj2->a[k] = jj2->a[j0]; + jj2->a[k].cnt = cnt; + jj2->a[k++].flag = flag; + j0 = j; + } + } + } + return J2; +} + +int mm_idx_jjump_read(mm_idx_t *mi, const char *fn, int flag, int min_sc) +{ + int32_t i; + mm_idx_intv_t *I; + mm_idx_jjump_t *J; + if (mi->h == 0) mm_idx_index_name(mi); + I = mm_idx_bed_read_merge(mi, fn, 1, min_sc); + J = mm_idx_bed2jjump(mi, I, flag); + for (i = 0; i < mi->n_seq; ++i) free(I[i].a); + free(I); + if (mi->J) { + mm_idx_jjump_t *J2; + J2 = mm_idx_jjump_merge(mi, mi->J, J2); + for (i = 0; i < mi->n_seq; ++i) { + free(mi->J[i].a); free(J[i].a); + } + free(mi->J); free(J); + mi->J = J2; + } else mi->J = J; + return 0; +} + static int32_t mm_idx_jump_get_core(int32_t n, const mm_idx_jjump1_t *a, int32_t x) // similar to mm_idx_find_intv() { int32_t s = 0, e = n; diff --git a/main.c b/main.c index 2ca8e69..0d725bd 100644 --- a/main.c +++ b/main.c @@ -84,6 +84,7 @@ static ko_longopt_t long_options[] = { { "pairing", ko_required_argument, 359 }, { "jump-min-match", ko_required_argument, 360 }, { "write-junc", ko_no_argument, 361 }, + { "jump-pass1", ko_required_argument, 362 }, { "dbg-seed-occ", ko_no_argument, 501 }, { "help", ko_no_argument, 'h' }, { "max-intron-len", ko_required_argument, 'G' }, @@ -133,7 +134,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, *jump_bed = 0, *fn_spsc = 0, *s, *alt_list = 0; + char *fnw = 0, *rg = 0, *fn_bed_junc = 0, *fn_bed_jump = 0, *fn_bed_pass1 = 0, *fn_spsc = 0, *s, *alt_list = 0; FILE *fp_help = stderr; mm_idx_reader_t *idx_rdr; mm_idx_t *mi; @@ -195,7 +196,7 @@ int main(int argc, char *argv[]) else if (c == 'R') rg = o.arg; else if (c == 'h') fp_help = stdout; else if (c == '2') opt.flag |= MM_F_2_IO_THREADS; - else if (c == 'j') jump_bed = o.arg; + else if (c == 'j') fn_bed_jump = o.arg; else if (c == 'J') { int t; t = atoi(o.arg); @@ -237,7 +238,7 @@ int main(int argc, char *argv[]) else if (c == 336) opt.flag |= MM_F_HARD_MLEVEL; // --hard-mask-level else if (c == 337) opt.max_sw_mat = mm_parse_num(o.arg); // --cap-sw-mat else if (c == 338) opt.max_qlen = mm_parse_num(o.arg); // --max-qlen - else if (c == 340) junc_bed = o.arg; // --junc-bed + else if (c == 340) fn_bed_junc = o.arg; // --junc-bed else if (c == 341) opt.junc_bonus = atoi(o.arg); // --junc-bonus else if (c == 358) opt.junc_pen = atoi(o.arg); // --junc-pen else if (c == 342) opt.flag |= MM_F_SAM_HIT_ONLY; // --sam-hit-only @@ -257,6 +258,7 @@ int main(int argc, char *argv[]) else if (c == 357) fn_spsc = o.arg; // --spsc else if (c == 360) opt.jump_min_match = mm_parse_num(o.arg); // --jump-min-match else if (c == 361) opt.flag |= MM_F_OUT_JUNC | MM_F_CIGAR; // --write-junc + else if (c == 362) fn_bed_pass1 = o.arg; // --jump-pass1 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"); @@ -458,16 +460,21 @@ int main(int argc, char *argv[]) __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), mi->n_seq); 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_bed_junc) { + mm_idx_bed_read(mi, fn_bed_junc, 1); if (mi->I == 0 && mm_verbose >= 2) fprintf(stderr, "[WARNING] failed to load the junction BED file\n"); } - if (jump_bed) { - mm_idx_bed_read2(mi, jump_bed, 1, 0, 1); + if (fn_bed_jump) { + mm_idx_jjump_read(mi, fn_bed_jump, MM_JUNC_ANNO, -1); if (mi->J == 0 && mm_verbose >= 2) fprintf(stderr, "[WARNING] failed to load the jump BED file\n"); } + if (fn_bed_pass1) { + mm_idx_jjump_read(mi, fn_bed_pass1, MM_JUNC_MISC, 5); + if (mi->J == 0 && mm_verbose >= 2) + fprintf(stderr, "[WARNING] failed to load the pass-1 jump BED file\n"); + } if (fn_spsc) { mm_idx_spsc_read(mi, fn_spsc, mm_max_spsc_bonus(&opt)); if (mi->spsc == 0 && mm_verbose >= 2) diff --git a/minimap.h b/minimap.h index f4aeda6..61829ca 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#define MM_VERSION "2.28-r1271-dirty" +#define MM_VERSION "2.28-r1272-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 diff --git a/mmpriv.h b/mmpriv.h index 4331e3a..e00a463 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -24,6 +24,9 @@ #define MM_SEED_SEG_SHIFT 48 #define MM_SEED_SEG_MASK (0xffULL<<(MM_SEED_SEG_SHIFT)) +#define MM_JUNC_ANNO 0x1 +#define MM_JUNC_MISC 0x2 + #ifndef kroundup32 #define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x)) #endif @@ -54,8 +57,9 @@ typedef struct { } mm_seg_t; typedef struct { - int32_t off, off2; - int32_t cnt, strand; + int32_t off, off2, cnt; + int16_t strand; + uint16_t flag; } mm_idx_jjump1_t; double cputime(void); @@ -81,14 +85,17 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int64_t opt_flag, int rep_len); void mm_write_junc(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r); +// indexing related in index.c void mm_idxopt_init(mm_idxopt_t *opt); const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n); int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f); int mm_idx_getseq2(const mm_idx_t *mi, int is_rev, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq); mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a, int is_qstrand); -int mm_idx_bed_read2(mm_idx_t *mi, const char *fn, int read_junc, int for_score, int for_jump); +int mm_idx_bed_read(mm_idx_t *mi, const char *fn, int read_junc); +int mm_idx_jjump_read(mm_idx_t *mi, const char *fn, int flag, int min_sc); const mm_idx_jjump1_t *mm_idx_jump_get(const mm_idx_t *db, int32_t cid, int32_t st, int32_t en, int32_t *n); +// chaining in lchain.c mm128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km); mm128_t *mg_lchain_rmq(int max_dist, int max_dist_inner, int bw, int max_chn_skip, int cap_rmq_size, int min_cnt, int min_sc, float chn_pen_gap, float chn_pen_skip,