From 819b3bf017e51f908e7d0023367556ae89220eec Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sun, 13 Apr 2025 13:55:15 -0400 Subject: [PATCH] r1274: a little better pass-1 Also fixed a bug in jump --- index.c | 11 +++++++---- jump.c | 31 +++++++++++++++++++++++-------- minimap.h | 2 +- 3 files changed, 31 insertions(+), 13 deletions(-) diff --git a/index.c b/index.c index 0b74f3d..988ac6a 100644 --- a/index.c +++ b/index.c @@ -902,7 +902,7 @@ static mm_idx_jjump_t *mm_idx_jjump_merge(const mm_idx_t *mi, const mm_idx_jjump int mm_idx_jjump_read(mm_idx_t *mi, const char *fn, int flag, int min_sc) { - int32_t i, n = 0; + int32_t i, j, n_anno = 0, n_misc = 0; mm_idx_intv_t *I; mm_idx_jjump_t *J; if (mi->h == 0) mm_idx_index_name(mi); @@ -919,10 +919,13 @@ int mm_idx_jjump_read(mm_idx_t *mi, const char *fn, int flag, int min_sc) free(mi->J); free(J); mi->J = J2; } else mi->J = J; - for (i = 0; i < mi->n_seq; ++i) - n += mi->J[i].n; + for (i = 0; i < mi->n_seq; ++i) { + for (j = 0; j < mi->J[i].n; ++j) + if (mi->J[i].a[j].flag & MM_JUNC_ANNO) ++n_anno; + else ++n_misc; + } if (mm_verbose >= 3) - fprintf(stderr, "[%s] there are %d splice positions in the index\n", __func__, n); + fprintf(stderr, "[%s] there are %d annotated and %d other splice positions in the index\n", __func__, n_anno, n_misc); return 0; } diff --git a/jump.c b/jump.c index d7a26c1..83e23da 100644 --- a/jump.c +++ b/jump.c @@ -10,7 +10,7 @@ static int32_t mm_jump_check(void *km, const mm_idx_t *mi, int32_t qlen, const u uint32_t cigar; if (!r->p || r->p->n_cigar <= 0) return -1; // only working with CIGAR clip = e == 0? r->qs : qlen - r->qe; - cigar = r->p->cigar[e == 0? 0 : r->p->n_cigar - 1]; + cigar = r->p->cigar[is_left? 0 : r->p->n_cigar - 1]; clen = (cigar&0xf) == MM_CIGAR_MATCH? cigar>>4 : 0; if (clen <= ext) return -1; if (is_left) { @@ -50,7 +50,8 @@ static uint8_t *mm_jump_get_qseq_seq(void *km, int32_t qlen, const uint8_t *qseq static void mm_jump_split_left(void *km, const mm_idx_t *mi, const mm_mapopt_t *opt, int32_t qlen, const uint8_t *qseq0, mm_reg1_t *r, int32_t ts_strand) { uint8_t *tseq = 0, *qseq = 0; - int32_t i, n, i0 = -1, m = 0, mm0 = 0, l; + int32_t i, n, l, i0, m, mm0; + int32_t i0_anno = -1, n_anno = 0, mm0_anno = 0, i0_misc = -1, n_misc = 0, mm0_misc = 0; int32_t ext = 1 + (opt->b + opt->a - 1) / opt->a + 1; int32_t clip = !r->rev? r->qs : qlen - r->qe; int32_t extt = clip < ext? clip : ext; @@ -84,9 +85,15 @@ 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) - i0 = i, mm0 = mm1 + mm2, ++m; // i0 points to the rightmost i + if (mm1 == 0 && mm2 <= 1) { + if (ai->flag & MM_JUNC_ANNO) + i0_anno = i, mm0_anno = mm1 + mm2, ++n_anno; // i0 points to the rightmost i + else + i0_misc = i, mm0_misc = mm1 + mm2, ++n_misc; + } } + if (n_anno > 0) m = n_anno, i0 = i0_anno, mm0 = mm0_anno; + else m = n_misc, i0 = i0_misc, mm0 = mm0_misc; kfree(km, tseq); l = m > 0? a[i0].off - r->rs : 0; // may be negative @@ -115,13 +122,14 @@ static void mm_jump_split_left(void *km, const mm_idx_t *mi, const mm_mapopt_t * static void mm_jump_split_right(void *km, const mm_idx_t *mi, const mm_mapopt_t *opt, int32_t qlen, const uint8_t *qseq0, mm_reg1_t *r, int32_t ts_strand) { uint8_t *tseq = 0, *qseq = 0; - int32_t i, n, i0 = -1, m = 0, mm0 = 0, l; + int32_t i, n, l, i0, m, mm0; + int32_t i0_anno = -1, n_anno = 0, mm0_anno = 0, i0_misc = -1, n_misc = 0, mm0_misc = 0; int32_t ext = 1 + (opt->b + opt->a - 1) / opt->a + 1; int32_t clip = !r->rev? qlen - r->qe : r->qs; int32_t extt = clip < ext? clip : ext; const mm_idx_jjump1_t *a; - if (mm_jump_check(km, mi, qlen, qseq0, r, ext + MM_MIN_EXON_LEN, 1) < 0) return; + if (mm_jump_check(km, mi, qlen, qseq0, r, ext + MM_MIN_EXON_LEN, 0) < 0) return; a = mm_idx_jump_get(mi, r->rid, r->re - ext, r->re + extt, &n); if (n == 0) return; @@ -150,10 +158,17 @@ static void mm_jump_split_right(void *km, const mm_idx_t *mi, const mm_mapopt_t if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3) ++mm1; if (mm1 == 0 && mm2 <= 1) { - if (i0 < 0) i0 = i, mm0 = mm1 + mm2; - ++m; + if (ai->flag & MM_JUNC_ANNO) { + if (i0_anno < 0) i0_anno = i, mm0_anno = mm1 + mm2; + ++n_anno; + } else { + if (i0_misc < 0) i0_misc = i, mm0_misc = mm1 + mm2; + ++n_misc; + } } } + if (n_anno > 0) m = n_anno, i0 = i0_anno, mm0 = mm0_anno; + else m = n_misc, i0 = i0_misc, mm0 = mm0_misc; kfree(km, tseq); l = m > 0? r->re - a[i0].off : 0; // may be negative diff --git a/minimap.h b/minimap.h index addc818..ab70d75 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#define MM_VERSION "2.28-r1273-dirty" +#define MM_VERSION "2.28-r1274-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