From 18778182391a4913dc43cb797cdece757b7a4926 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sun, 6 Apr 2025 09:33:51 -0400 Subject: [PATCH] r1255: moved jump code to a separate file --- Makefile | 7 ++- hit.c | 167 ------------------------------------------------------ jump.c | 165 +++++++++++++++++++++++++++++++++++++++++++++++++++++ minimap.h | 2 +- 4 files changed, 170 insertions(+), 171 deletions(-) create mode 100644 jump.c diff --git a/Makefile b/Makefile index 17b13b6..d2d30e1 100644 --- a/Makefile +++ b/Makefile @@ -2,7 +2,7 @@ CFLAGS= -g -Wall -O2 -Wc++-compat #-Wextra CPPFLAGS= -DHAVE_KALLOC INCLUDES= OBJS= kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o \ - lchain.o align.o hit.o seed.o map.o format.o pe.o esterr.o splitidx.o \ + lchain.o align.o hit.o seed.o jump.o map.o format.o pe.o esterr.o splitidx.o \ ksw2_ll_sse.o PROG= minimap2 PROG_EXTRA= sdust minimap2-lite @@ -115,8 +115,9 @@ esterr.o: mmpriv.h minimap.h bseq.h kseq.h example.o: minimap.h kseq.h format.o: kalloc.h mmpriv.h minimap.h bseq.h kseq.h hit.o: mmpriv.h minimap.h bseq.h kseq.h kalloc.h khash.h -index.o: kthread.h bseq.h minimap.h mmpriv.h kseq.h kvec.h kalloc.h khash.h -index.o: ksort.h +index.o: kthread.h bseq.h minimap.h mmpriv.h kseq.h ksw2.h kalloc.h kvec.h +index.o: khash.h ksort.h +jump.o: mmpriv.h minimap.h bseq.h kseq.h kalloc.o: kalloc.h ksw2_extd2_sse.o: ksw2.h kalloc.h ksw2_exts2_sse.o: ksw2.h kalloc.h diff --git a/hit.c b/hit.c index ecfae91..f5a50a2 100644 --- a/hit.c +++ b/hit.c @@ -464,170 +464,3 @@ void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int ma } mm_set_inv_mapq(km, n_regs, regs); } - -/****************************** - * Instructed split alignment * - ******************************/ - -#define MM_MIN_EXON_LEN 20 - -static int32_t mm_jump_check(void *km, const mm_idx_t *mi, int32_t qlen, const uint8_t *qseq0, const mm_reg1_t *r, int32_t ext, int32_t is_left) // TODO: check close N -{ - int32_t clip, clen, e = !r->rev ^ !is_left; // 0 for left of the alignment; 1 for right - 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]; - clen = (cigar&0xf) == MM_CIGAR_MATCH? cigar>>4 : 0; - if (clen <= ext) return -1; - if (is_left) { - if (clip >= r->rs) return -1; // no space to jump - } else { - if (clip >= mi->seq[r->rid].len - r->re) return -1; // no space to jump - } - return 0; -} - -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; - if (!r->rev) { - if (is_left) memcpy(qseq, qseq0, ql0); - else memcpy(qseq, &qseq0[qlen - ql0], ql0); - } else { - if (is_left) - for (i = qlen - 1, k = 0; i >= qlen - ql0; --i) - qseq[k++] = qseq0[i] >= 4? qseq0[i] : 3 - qseq0[i]; - else - for (i = ql0 - 1, k = 0; i >= 0; --i) - qseq[k++] = qseq0[i] >= 4? qseq0[i] : 3 - qseq0[i]; - } - return 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, l; - 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; - const mm_idx_jjump1_t *a; - - if (mm_jump_check(km, mi, qlen, qseq0, r, ext + MM_MIN_EXON_LEN, 1) < 0) return; - a = mm_idx_jump_get(mi, r->rid, r->rs - extt, r->rs + ext, &n); - if (n == 0) return; - - 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); - if (ts_strand * ai->strand < 0) continue; // wrong strand - if (ai->off2 >= ai->off) continue; // wrong direction - if (ai->off2 < clip + ext) continue; // not long enough - if (tseq == 0) { - tseq = Kcalloc(km, uint8_t, (clip + ext) * 2); // tseq and qseq are allocated together - qseq = tseq + clip + ext; - mm_jump_get_qseq_seq(km, qlen, qseq0, r, 1, clip + ext, qseq); - } - tl1 = clip + (ai->off - r->rs); - tlen = mm_idx_getseq2(mi, 0, r->rid, ai->off, r->rs + ext, &tseq[tl1]); - assert(tlen == r->rs + ext - ai->off); - tlen = mm_idx_getseq2(mi, 0, r->rid, ai->off2 - tl1, ai->off2, tseq); - assert(tlen == tl1); - for (j = 0, mm1 = 0; j < tl1; ++j) - if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3) - ++mm1; - 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, ++m; // i0 points to the rightmost i - } - kfree(km, tseq); - - l = m > 0? a[i0].off - r->rs : 0; // may be negative - if (m == 1 && clip + l >= opt->jump_min_alen) { // add one more exon - mm_enlarge_cigar(r, 2); - memmove(r->p->cigar + 2, r->p->cigar, r->p->n_cigar * 4); - r->p->cigar[0] = (clip + l) << 4 | MM_CIGAR_MATCH; - r->p->cigar[1] = (a[i0].off - a[i0].off2) << 4 | MM_CIGAR_N_SKIP; - r->p->cigar[2] = ((r->p->cigar[2]>>4) - l) << 4 | MM_CIGAR_MATCH; - r->p->n_cigar += 2; - r->rs = a[i0].off2 - (clip + l); - if (!r->rev) r->qs = 0; - else r->qe = qlen; - } else if (m > 0 && a[i0].off > r->rs) { // trim by l; l is always positive - r->p->cigar[0] -= l << 4 | MM_CIGAR_MATCH; - r->rs += l; - if (!r->rev) r->qs += l; - else r->qe -= l; - } -} - -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, l; - 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; - a = mm_idx_jump_get(mi, r->rid, r->re - ext, r->re + extt, &n); - if (n == 0) return; - - 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); - 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 - if (tseq == 0) { - tseq = Kcalloc(km, uint8_t, (clip + ext) * 2); // tseq and qseq are allocated together - qseq = tseq + clip + ext; - mm_jump_get_qseq_seq(km, qlen, qseq0, r, 0, clip + ext, qseq); - } - tl1 = clip + (r->re - ai->off); - tlen = mm_idx_getseq2(mi, 0, r->rid, r->re - ext, ai->off, tseq); - assert(tlen == ai->off - (r->re - ext)); - tlen = mm_idx_getseq2(mi, 0, r->rid, ai->off2, ai->off2 + tl1, &tseq[clip + ext - tl1]); - assert(tlen == tl1); - for (j = 0, mm2 = 0; j < clip + ext - tl1; ++j) - if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3) - ++mm2; - for (mm1 = 0; j < clip + ext; ++j) - if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3) - ++mm1; - if (mm1 == 0 && mm2 == 1) - i0 = i0 >= 0? i0 : i, ++m; // i0 points to the leftmost i - } - kfree(km, tseq); - - l = m > 0? r->re - a[i0].off : 0; // may be negative - if (m == 1 && clip + l >= opt->jump_min_alen) { // add one more exon - mm_enlarge_cigar(r, 2); - memmove(r->p->cigar + 2, r->p->cigar, r->p->n_cigar * 4); - r->p->cigar[r->p->n_cigar - 1] = ((r->p->cigar[r->p->n_cigar - 1]>>4) - l) << 4 | MM_CIGAR_MATCH; - r->p->cigar[r->p->n_cigar] = (a[i0].off2 - a[i0].off) << 4 | MM_CIGAR_N_SKIP; - r->p->cigar[r->p->n_cigar + 1] = (clip + l) << 4 | MM_CIGAR_MATCH; - r->p->n_cigar += 2; - r->re = a[i0].off2 + (clip + l); - if (!r->rev) r->qe = qlen; - else r->qs = 0; - } else if (m > 0 && r->re > a[i0].off) { // trim by l; l is always positive - r->p->cigar[r->p->n_cigar - 1] -= l << 4 | MM_CIGAR_MATCH; - r->re -= l; - if (!r->rev) r->qe -= l; - else r->qs += l; - } -} - -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) -{ - assert((opt->flag & MM_F_EQX) == 0); - mm_jump_split_left(km, mi, opt, qlen, qseq, r, ts_strand); - mm_jump_split_right(km, mi, opt, qlen, qseq, r, ts_strand); -} diff --git a/jump.c b/jump.c new file mode 100644 index 0000000..f4afdef --- /dev/null +++ b/jump.c @@ -0,0 +1,165 @@ +#include "mmpriv.h" +#include "kalloc.h" + +#define MM_MIN_EXON_LEN 20 + +static int32_t mm_jump_check(void *km, const mm_idx_t *mi, int32_t qlen, const uint8_t *qseq0, const mm_reg1_t *r, int32_t ext, int32_t is_left) // TODO: check close N +{ + int32_t clip, clen, e = !r->rev ^ !is_left; // 0 for left of the alignment; 1 for right + 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]; + clen = (cigar&0xf) == MM_CIGAR_MATCH? cigar>>4 : 0; + if (clen <= ext) return -1; + if (is_left) { + if (clip >= r->rs) return -1; // no space to jump + } else { + if (clip >= mi->seq[r->rid].len - r->re) return -1; // no space to jump + } + return 0; +} + +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; + if (!r->rev) { + if (is_left) memcpy(qseq, qseq0, ql0); + else memcpy(qseq, &qseq0[qlen - ql0], ql0); + } else { + if (is_left) + for (i = qlen - 1, k = 0; i >= qlen - ql0; --i) + qseq[k++] = qseq0[i] >= 4? qseq0[i] : 3 - qseq0[i]; + else + for (i = ql0 - 1, k = 0; i >= 0; --i) + qseq[k++] = qseq0[i] >= 4? qseq0[i] : 3 - qseq0[i]; + } + return 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, l; + 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; + const mm_idx_jjump1_t *a; + + if (mm_jump_check(km, mi, qlen, qseq0, r, ext + MM_MIN_EXON_LEN, 1) < 0) return; + a = mm_idx_jump_get(mi, r->rid, r->rs - extt, r->rs + ext, &n); + if (n == 0) return; + + 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); + if (ts_strand * ai->strand < 0) continue; // wrong strand + if (ai->off2 >= ai->off) continue; // wrong direction + if (ai->off2 < clip + ext) continue; // not long enough + if (tseq == 0) { + tseq = Kcalloc(km, uint8_t, (clip + ext) * 2); // tseq and qseq are allocated together + qseq = tseq + clip + ext; + mm_jump_get_qseq_seq(km, qlen, qseq0, r, 1, clip + ext, qseq); + } + tl1 = clip + (ai->off - r->rs); + tlen = mm_idx_getseq2(mi, 0, r->rid, ai->off, r->rs + ext, &tseq[tl1]); + assert(tlen == r->rs + ext - ai->off); + tlen = mm_idx_getseq2(mi, 0, r->rid, ai->off2 - tl1, ai->off2, tseq); + assert(tlen == tl1); + for (j = 0, mm1 = 0; j < tl1; ++j) + if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3) + ++mm1; + 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, ++m; // i0 points to the rightmost i + } + kfree(km, tseq); + + l = m > 0? a[i0].off - r->rs : 0; // may be negative + if (m == 1 && clip + l >= opt->jump_min_alen) { // add one more exon + mm_enlarge_cigar(r, 2); + memmove(r->p->cigar + 2, r->p->cigar, r->p->n_cigar * 4); + r->p->cigar[0] = (clip + l) << 4 | MM_CIGAR_MATCH; + r->p->cigar[1] = (a[i0].off - a[i0].off2) << 4 | MM_CIGAR_N_SKIP; + r->p->cigar[2] = ((r->p->cigar[2]>>4) - l) << 4 | MM_CIGAR_MATCH; + r->p->n_cigar += 2; + r->rs = a[i0].off2 - (clip + l); + if (!r->rev) r->qs = 0; + else r->qe = qlen; + } else if (m > 0 && a[i0].off > r->rs) { // trim by l; l is always positive + r->p->cigar[0] -= l << 4 | MM_CIGAR_MATCH; + r->rs += l; + if (!r->rev) r->qs += l; + else r->qe -= l; + } +} + +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, l; + 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; + a = mm_idx_jump_get(mi, r->rid, r->re - ext, r->re + extt, &n); + if (n == 0) return; + + 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); + 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 + if (tseq == 0) { + tseq = Kcalloc(km, uint8_t, (clip + ext) * 2); // tseq and qseq are allocated together + qseq = tseq + clip + ext; + mm_jump_get_qseq_seq(km, qlen, qseq0, r, 0, clip + ext, qseq); + } + tl1 = clip + (r->re - ai->off); + tlen = mm_idx_getseq2(mi, 0, r->rid, r->re - ext, ai->off, tseq); + assert(tlen == ai->off - (r->re - ext)); + tlen = mm_idx_getseq2(mi, 0, r->rid, ai->off2, ai->off2 + tl1, &tseq[clip + ext - tl1]); + assert(tlen == tl1); + for (j = 0, mm2 = 0; j < clip + ext - tl1; ++j) + if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3) + ++mm2; + for (mm1 = 0; j < clip + ext; ++j) + if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3) + ++mm1; + if (mm1 == 0 && mm2 == 1) + i0 = i0 >= 0? i0 : i, ++m; // i0 points to the leftmost i + } + kfree(km, tseq); + + l = m > 0? r->re - a[i0].off : 0; // may be negative + if (m == 1 && clip + l >= opt->jump_min_alen) { // add one more exon + mm_enlarge_cigar(r, 2); + memmove(r->p->cigar + 2, r->p->cigar, r->p->n_cigar * 4); + r->p->cigar[r->p->n_cigar - 1] = ((r->p->cigar[r->p->n_cigar - 1]>>4) - l) << 4 | MM_CIGAR_MATCH; + r->p->cigar[r->p->n_cigar] = (a[i0].off2 - a[i0].off) << 4 | MM_CIGAR_N_SKIP; + r->p->cigar[r->p->n_cigar + 1] = (clip + l) << 4 | MM_CIGAR_MATCH; + r->p->n_cigar += 2; + r->re = a[i0].off2 + (clip + l); + if (!r->rev) r->qe = qlen; + else r->qs = 0; + } else if (m > 0 && r->re > a[i0].off) { // trim by l; l is always positive + r->p->cigar[r->p->n_cigar - 1] -= l << 4 | MM_CIGAR_MATCH; + r->re -= l; + if (!r->rev) r->qe -= l; + else r->qs += l; + } +} + +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) +{ + assert((opt->flag & MM_F_EQX) == 0); + mm_jump_split_left(km, mi, opt, qlen, qseq, r, ts_strand); + mm_jump_split_right(km, mi, opt, qlen, qseq, r, ts_strand); +} diff --git a/minimap.h b/minimap.h index cb2ff32..59d30cc 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#define MM_VERSION "2.28-r1251-dirty" +#define MM_VERSION "2.28-r1255-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