From a8094ad859a37b8dada546e625a49f51ca8d868f Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sun, 6 Apr 2025 10:33:57 -0400 Subject: [PATCH] r1256: working for one example, but still buggy --- jump.c | 30 ++++++++++++++++++++---------- main.c | 9 ++++++++- map.c | 4 ++++ minimap.h | 2 +- mmpriv.h | 3 +++ options.c | 2 +- 6 files changed, 37 insertions(+), 13 deletions(-) diff --git a/jump.c b/jump.c index f4afdef..4d5f27f 100644 --- a/jump.c +++ b/jump.c @@ -1,3 +1,4 @@ +#include #include "mmpriv.h" #include "kalloc.h" @@ -22,17 +23,26 @@ static int32_t mm_jump_check(void *km, const mm_idx_t *mi, int32_t qlen, const u static uint8_t *mm_jump_get_qseq_seq(void *km, int32_t qlen, const uint8_t *qseq0, const mm_reg1_t *r, int32_t is_left, int32_t ql0, uint8_t *qseq) { - int32_t i, k; + extern unsigned char seq_nt4_table[256]; + int32_t i, k = 0; if (!r->rev) { - if (is_left) memcpy(qseq, qseq0, ql0); - else memcpy(qseq, &qseq0[qlen - ql0], ql0); + if (is_left) + for (i = 0; i < ql0; ++i) + qseq[k++] = seq_nt4_table[(uint8_t)qseq0[i]]; + else + for (i = qlen - ql0; i < qlen; ++i) + qseq[k++] = seq_nt4_table[(uint8_t)qseq0[i]]; } else { if (is_left) - for (i = qlen - 1, k = 0; i >= qlen - ql0; --i) - qseq[k++] = qseq0[i] >= 4? qseq0[i] : 3 - qseq0[i]; + for (i = qlen - 1; i >= qlen - ql0; --i) { + uint8_t c = seq_nt4_table[(uint8_t)qseq0[i]]; + qseq[k++] = c >= 4? c : 3 - c; + } else - for (i = ql0 - 1, k = 0; i >= 0; --i) - qseq[k++] = qseq0[i] >= 4? qseq0[i] : 3 - qseq0[i]; + for (i = ql0 - 1; i >= 0; --i) { + uint8_t c = seq_nt4_table[(uint8_t)qseq0[i]]; + qseq[k++] = c >= 4? c : 3 - c; + } } return qseq; } @@ -73,7 +83,7 @@ static void mm_jump_split_left(void *km, const mm_idx_t *mi, const mm_mapopt_t * for (mm2 = 0; j < clip + ext; ++j) if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3) ++mm2; - if (mm1 == 0 && mm2 == 1) + if (mm1 == 0 && mm2 <= 1) i0 = i, ++m; // i0 points to the rightmost i } kfree(km, tseq); @@ -113,7 +123,7 @@ static void mm_jump_split_right(void *km, const mm_idx_t *mi, const mm_mapopt_t for (i = 0; i < n; ++i) { // traverse possible jumps const mm_idx_jjump1_t *ai = &a[i]; int32_t tlen, tl1, j, mm1, mm2; - assert(ai->off >= r->rs - extt && ai->off < r->rs + ext); + assert(ai->off >= r->re - ext && ai->off < r->re + extt); if (ts_strand * ai->strand < 0) continue; // wrong strand if (ai->off2 <= ai->off) continue; // wrong direction if (ai->off2 + clip + ext > mi->seq[r->rid].len) continue; // not long enough @@ -133,7 +143,7 @@ static void mm_jump_split_right(void *km, const mm_idx_t *mi, const mm_mapopt_t for (mm1 = 0; j < clip + ext; ++j) if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3) ++mm1; - if (mm1 == 0 && mm2 == 1) + if (mm1 == 0 && mm2 <= 1) i0 = i0 >= 0? i0 : i, ++m; // i0 points to the leftmost i } kfree(km, tseq); diff --git a/main.c b/main.c index b4d3ce4..0e44d3c 100644 --- a/main.c +++ b/main.c @@ -82,6 +82,7 @@ static ko_longopt_t long_options[] = { { "spsc", ko_required_argument, 357 }, { "junc-pen", ko_required_argument, 358 }, { "pe-ind-chain", ko_no_argument, 359 }, + { "jump-bed", ko_required_argument, 360 }, { "dbg-seed-occ", ko_no_argument, 501 }, { "help", ko_no_argument, 'h' }, { "max-intron-len", ko_required_argument, 'G' }, @@ -131,7 +132,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, *fn_spsc = 0, *s, *alt_list = 0; + char *fnw = 0, *rg = 0, *junc_bed = 0, *jump_bed = 0, *fn_spsc = 0, *s, *alt_list = 0; FILE *fp_help = stderr; mm_idx_reader_t *idx_rdr; mm_idx_t *mi; @@ -254,6 +255,7 @@ int main(int argc, char *argv[]) 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 == 359) opt.flag |= MM_F_PE_IND; // --pe-ind-chain + else if (c == 360) jump_bed = o.arg; // --jump-bed 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"); @@ -440,6 +442,11 @@ int main(int argc, char *argv[]) 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 (mi->J == 0 && mm_verbose >= 2) + fprintf(stderr, "[WARNING] failed to load the 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/map.c b/map.c index 385e1dd..2546c80 100644 --- a/map.c +++ b/map.c @@ -359,6 +359,10 @@ void mm_map_frag_core(const mm_idx_t *mi, int n_segs, const int *qlens, const ch kfree(b->km, u); kfree(b->km, mini_pos); + if (mi->J && n_segs == 1 && is_splice) + for (i = 0; i < n_regs0; ++i) + mm_jump_split(b->km, mi, opt, qlens[0], (const uint8_t*)seqs[0], ®s0[i], 0); + if (b->km) { km_stat(b->km, &kmst); if (mm_dbg_flag & MM_DBG_PRINT_QNAME) diff --git a/minimap.h b/minimap.h index 59d30cc..95496a4 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#define MM_VERSION "2.28-r1255-dirty" +#define MM_VERSION "2.28-r1256-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 bc71d7b..eec4882 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -84,6 +84,7 @@ 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); 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); 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, @@ -115,6 +116,8 @@ mm_seg_t *mm_seg_gen(void *km, uint32_t hash, int n_segs, const int *qlens, int void mm_seg_free(void *km, int n_segs, mm_seg_t *segs); void mm_pair(void *km, int max_gap_ref, int dp_bonus, int sub_diff, int match_sc, const int *qlens, int *n_regs, mm_reg1_t **regs); +void mm_jump_split(void *km, const mm_idx_t *mi, const mm_mapopt_t *opt, int32_t qlen, const uint8_t *qseq, mm_reg1_t *r, int32_t ts_strand); + FILE *mm_split_init(const char *prefix, const mm_idx_t *mi); mm_idx_t *mm_split_merge_prep(const char *prefix, int n_splits, FILE **fp, uint32_t *n_seq_part); int mm_split_merge(int n_segs, const char **fn, const mm_mapopt_t *opt, int n_split_idx); diff --git a/options.c b/options.c index 22927db..072eaa3 100644 --- a/options.c +++ b/options.c @@ -63,7 +63,7 @@ void mm_mapopt_init(mm_mapopt_t *opt) opt->pe_ori = 0; // FF opt->pe_bonus = 33; - opt->jump_min_alen = 10; + opt->jump_min_alen = 5; } void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi)