fork/minimap2
1
0
Fork 0
mirror of https://github.com/lh3/minimap2.git synced 2026-10-07 21:08:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

r1274: a little better pass-1

Browse Source 
Also fixed a bug in jump
This commit is contained in:
Heng Li
2025-04-13 13:55:15 -04:00
parent e88110463a
commit 819b3bf017
3 changed files with 31 additions and 13 deletions
Download Patch File Download Diff File Expand all files Collapse all files
index.c
+7 -4
View File
@@ -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) 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_intv_t *I;
mm_idx_jjump_t *J; mm_idx_jjump_t *J;
if (mi->h == 0) mm_idx_index_name(mi); 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); free(mi->J); free(J);
mi->J = J2; mi->J = J2;
} else mi->J = J; } else mi->J = J;
for (i = 0; i < mi->n_seq; ++i) for (i = 0; i < mi->n_seq; ++i) {
n += mi->J[i].n; 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) 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; return 0;
} }
jump.c
+23 -8
View File
@@ -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; uint32_t cigar;
if (!r->p || r->p->n_cigar <= 0) return -1; // only working with CIGAR if (!r->p || r->p->n_cigar <= 0) return -1; // only working with CIGAR
clip = e == 0? r->qs : qlen - r->qe; 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; clen = (cigar&0xf) == MM_CIGAR_MATCH? cigar>>4 : 0;
if (clen <= ext) return -1; if (clen <= ext) return -1;
if (is_left) { 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) 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; 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 ext = 1 + (opt->b + opt->a - 1) / opt->a + 1;
int32_t clip = !r->rev? r->qs : qlen - r->qe; int32_t clip = !r->rev? r->qs : qlen - r->qe;
int32_t extt = clip < ext? clip : ext; 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) for (mm2 = 0; j < clip + ext; ++j)
if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3) if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3)
++mm2; ++mm2;
if (mm1 == 0 && mm2 <= 1) if (mm1 == 0 && mm2 <= 1) {
i0 = i, mm0 = mm1 + mm2, ++m; // i0 points to the rightmost i 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); kfree(km, tseq);
l = m > 0? a[i0].off - r->rs : 0; // may be negative 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) 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; 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 ext = 1 + (opt->b + opt->a - 1) / opt->a + 1;
int32_t clip = !r->rev? qlen - r->qe : r->qs; int32_t clip = !r->rev? qlen - r->qe : r->qs;
int32_t extt = clip < ext? clip : ext; int32_t extt = clip < ext? clip : ext;
const mm_idx_jjump1_t *a; 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); a = mm_idx_jump_get(mi, r->rid, r->re - ext, r->re + extt, &n);
if (n == 0) return; 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) if (qseq[j] != tseq[j] || qseq[j] > 3 || tseq[j] > 3)
++mm1; ++mm1;
if (mm1 == 0 && mm2 <= 1) { if (mm1 == 0 && mm2 <= 1) {
if (i0 < 0) i0 = i, mm0 = mm1 + mm2; if (ai->flag & MM_JUNC_ANNO) {
++m; 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); kfree(km, tseq);
l = m > 0? r->re - a[i0].off : 0; // may be negative l = m > 0? r->re - a[i0].off : 0; // may be negative
minimap.h
+1 -1
View File
@@ -5,7 +5,7 @@
#include <stdio.h> #include <stdio.h>
#include <sys/types.h> #include <sys/types.h>
#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_DIAG (0x001LL) // no exact diagonal hit
#define MM_F_NO_DUAL (0x002LL) // skip pairs where query name is lexicographically larger than target name #define MM_F_NO_DUAL (0x002LL) // skip pairs where query name is lexicographically larger than target name