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

r1272: support pass-1 junction processing

Browse Source 
NOT tested yet
This commit is contained in:
Heng Li
2025-04-12 23:15:04 -04:00
parent d930ea94ad
commit fb81e150f2
4 changed files with 125 additions and 59 deletions
Download Patch File Download Diff File Expand all files Collapse all files

148
index.c
View File

@@ -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;

21
main.c
View File

@@ -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)

2
minimap.h
View File

@@ -5,7 +5,7 @@
#include <stdio.h>
#include <sys/types.h>
#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

13
mmpriv.h
View File

@@ -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,