8000 feat(map): expose --max-sub-occ and --min-sub-occ CLI ablation flags by nh13 · Pull Request #2 · nh13/minibwa · GitHub
[go: up one dir, main page]
Skip to content

feat(map): expose --max-sub-occ and --min-sub-occ CLI ablation flags #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
nh13 wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

feat(map): expose --max-sub-occ and --min-sub-occ CLI ablation flags #2

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions map-algo.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -685,7 +685,7 @@ mb_hit_t *mb_map(const mb_opt_t *opt, const mb_idx_t *idx, int32_t qlen, const c
seq[i] = kom_nt4_table[(uint8_t)seq0[i]];
if (mt != L2B_METH_NONE)
l2b_meth_convert(mt, qlen, seq);
mb_seed_intv(b->km, idx->bwt, qlen, seq, opt->min_len, opt->max_sub_occ, &u);
mb_seed_intv(b->km, idx->bwt, qlen, seq, opt->min_len, opt->max_sub_occ, opt->min_sub_occ, &u);
kfree(b->km, seq);
ret = mb_map_sai(&opt_adap, idx, qlen, seq0, mt, &u, n_hit_, b, qname);
if (b0 == 0) mb_tbuf_destroy(b);
Expand Down Expand Up @@ -730,7 +730,7 @@ mb_hit_t **mb_map_batch(const mb_opt_t *opt, const mb_idx_t *idx, int32_t n_seq,

// batch SMEM for sub-batch
memset(sai, 0, sb_n * sizeof(mb_sai_v));
mb_seed_intv_batch(km, idx->bwt, sb_n, &qlen[sb_st], seq4, opt->min_len, opt->max_sub_occ, sai);
mb_seed_intv_batch(km, idx->bwt, sb_n, &qlen[sb_st], seq4, opt->min_len, opt->max_sub_occ, opt->min_sub_occ, sai);
for (k = 0; k < sb_n; ++k) kfree(km, seq4[k]);

// map each sequence in sub-batch
Expand Down
14 changes: 13 additions & 1 deletion map-main.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@ static void worker_for_se_batch(void *data, long i, int tid)
}
}
assert(p == n);
mb_seed_intv_batch(km, idx->bwt, n, len, seq, opt->min_len, opt->max_sub_occ, sai);
mb_seed_intv_batch(km, idx->bwt, n, len, seq, opt->min_len, opt->max_sub_occ, opt->min_sub_occ, sai);
kfree(km, seq);
kfree(km, len);
kfree(km, buf);
Expand Down Expand Up @@ -346,6 +346,8 @@ static ko_longopt_t long_options[] = {
{ "mmap", ko_optional_argument, 313 },
{ "xa-ratio", ko_required_argument, 314 },
{ "outs", ko_required_argument, 315 },
{ "max-sub-occ", ko_required_argument, 320 }, // ablation: 0 disables Pass-2 sub-SMEM reseeding
{ "min-sub-occ", ko_required_argument, 321 }, // ablation: skip Pass-2 for SMEMs with SA-interval size < N (default 1)
{ "dbg-aln-seq", ko_no_argument, 601 },
{ "dbg-anchor", ko_no_argument, 602 },
{ "dbg-seed", ko_no_argument, 603 },
Expand Down Expand Up @@ -379,6 +381,8 @@ static int usage_map(FILE *fp, const mb_opt_t *opt)
fprintf(fp, " -p FLOAT min secondary-to-primary score ratio [%g]\n", opt->pri_ratio);
fprintf(fp, " -N INT retain at most INT secondary alignments [%d]\n", opt->best_n);
fprintf(fp, " --chain-only perform chaining only without base alignment; force -fP\n");
fprintf(fp, " --max-sub-occ=INT reseed Pass-2 sub-SMEMs only when SA-interval size <= INT (0 disables Pass-2) [%d]\n", opt->max_sub_occ);
fprintf(fp, " --min-sub-occ=INT skip Pass-2 reseeding when SA-interval size < INT [%d]\n", opt->min_sub_occ);
fprintf(fp, " -x STR preset (sr, lr or adap for mixed short/long reads) [adap]\n");
fprintf(fp, " Alignment:\n");
fprintf(fp, " -A INT matching score [%d]\n", opt->a);
Expand Down Expand Up @@ -512,6 +516,10 @@ int main_map(int argc, char *argv[])
if (o.arg != 0 && strcmp(o.arg, "lite") == 0) mmap_preload = 0;
} else if (c == 314 || c == 315) { // --outs or --xa-ratio
mo.out_s = atof(o.arg);
} else if (c == 320) { // --max-sub-occ
mo.max_sub_occ = atoi(o.arg);
} else if (c == 321) { // --min-sub-occ
mo.min_sub_occ = atoi(o.arg);
} else if (c == 601) { // --dbg-aln-seq
kom_dbg_flag |= MB_DBG_ALN_SEQ;
} else if (c == 602) { // --dbg-anchor
Expand Down Expand Up @@ -552,6 +560,10 @@ int main_map(int argc, char *argv[])
}
}
if (mo.flag & MB_F_NO_ALN) mo.flag |= MB_F_NO_PAIRING | MB_F_PAF;
if (mo.max_sub_occ > 0 && mo.min_sub_occ > mo.max_sub_occ) {
fprintf(stderr, "[ERROR] --min-sub-occ (%d) must not exceed --max-sub-occ (%d)\n", mo.min_sub_occ, mo.max_sub_occ);
return 1;
}
if (argc - o.ind < 2)
return usage_map(stderr, &mo);

Expand Down
4 changes: 2 additions & 2 deletions mbpriv.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,8 +53,8 @@ void mb_opt_adap(const mb_opt_t *opt0, int32_t len, mb_opt_t *opt);
void mb_bwtgen(const char *fn_pac, const char *fn_bwt, int block_size);

// defined in seed.c
void mb_seed_intv(void *km, const mb_bwt_t *bwt, int32_t len, const uint8_t *seq, int32_t min_len, int32_t max_sub_occ, mb_sai_v *v);
void mb_seed_intv_batch(void *km, const mb_bwt_t *bwt, int32_t n_seq, const int32_t *len, uint8_t *const* seq, int32_t min_len, int32_t max_sub_occ, mb_sai_v *v);
void mb_seed_intv(void *km, const mb_bwt_t *bwt, int32_t len, const uint8_t *seq, int32_t min_len, int32_t max_sub_occ, int32_t min_sub_occ, mb_sai_v *v);
void mb_seed_intv_batch(void *km, const mb_bwt_t *bwt, int32_t n_seq, const int32_t *len, uint8_t *const* seq, int32_t min_len, int32_t max_sub_occ, int32_t min_sub_occ, mb_sai_v *v);
double mb_anchor(void *km, const mb_idx_t *idx, mb_sai_v *u, int32_t min_len, int32_t qlen, const uint8_t *qseq, l2b_meth_t mt, int32_t max_occ, mb_anchor_v *v);
void mb_anchor_sort(const l2b_t *l2b, int64_t n_a, mb_anchor_t *a);

Expand Down
3 changes: 2 additions & 1 deletion minibwa.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,8 @@ typedef struct {
uint64_t flag;
// seeding options
int32_t min_len; // min seed length
int32_t max_sub_occ; // look for shorter seed if smem occ below this value
int32_t max_sub_occ; // run Pass-2 sub-SMEM reseeding only when SA-interval size <= this value (0 disables Pass-2)
int32_t min_sub_occ; // run Pass-2 sub-SMEM reseeding only when SA-interval size >= this value (default 1 = no lower gate)
int32_t max_occ; // max interval occurrence
// general algorithm options
int32_t bw, bw_long; // bandwidth
Expand Down
1 change: 1 addition & 0 deletions options.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ static void mb_opt_reset(mb_opt_t *opt)
opt->max_sr_len = 325;
// seeding options
opt->max_sub_occ = 10;
opt->min_sub_occ = 1;
opt->max_occ = 250;
// chaining options
opt->max_chain_skip = 25;
Expand Down
8 changes: 4 additions & 4 deletions seed.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ KRADIX_SORT_INIT(mb_anchor, mb_anchor_t, key_anchor, 8)
* Seeding *
***********/

void mb_seed_intv(void *km, const mb_bwt_t *bwt, int32_t len, const uint8_t *seq, int32_t min_len, int32_t max_sub_occ, mb_sai_v *v)
void mb_seed_intv(void *km, const mb_bwt_t *bwt, int32_t len, const uint8_t *seq, int32_t min_len, int32_t max_sub_occ, int32_t min_sub_occ, mb_sai_v *v)
{
int64_t x = 0, i, n_a0;
mb_sai_t p;
Expand All @@ -39,7 +39,7 @@ void mb_seed_intv(void *km, const mb_bwt_t *bwt, int32_t len, const uint8_t *seq
for (i = 0; i < n_a0; ++i) { // pass 2: sub-SMEMs
int32_t sub_min_len;
uint32_t st = v->a[i].info>>32, en = (uint32_t)v->a[i].info;
if (en - st < min_len * 2 || v->a[i].size > max_sub_occ)
if (en - st < min_len * 2 || v->a[i].size > max_sub_occ || v->a[i].size < min_sub_occ)
continue;
x = st;
sub_min_len = (en - st) / 2 > min_len? (en - st) / 2 : min_len;
Expand All @@ -53,7 +53,7 @@ void mb_seed_intv(void *km, const mb_bwt_t *bwt, int32_t len, const uint8_t *seq
}
}

void mb_seed_intv_batch(void *km, const mb_bwt_t *bwt, int32_t n_seq, const int32_t *len, uint8_t *const* seq, int32_t min_len, int32_t max_sub_occ, mb_sai_v *v)
void mb_seed_intv_batch(void *km, const mb_bwt_t *bwt, int32_t n_seq, const int32_t *len, uint8_t *const* seq, int32_t min_len, int32_t max_sub_occ, int32_t min_sub_occ, mb_sai_v *v)
{ // identical to mb_seed_intv() though the order of intervals is often different
const int max_batch_size = 50;
mb_smem_entry_t *s;
Expand Down Expand Up @@ -82,7 +82,7 @@ void mb_seed_intv_batch(void *km, const mb_bwt_t *bwt, int32_t n_seq, const int3
for (j = 0; j < nv[i]; ++j) {
mb_smem_entry_t *t;
uint32_t st = v[i].a[j].info>>32, en = (uint32_t)v[i].a[j].info;
if (en - st < min_len * 2 || v[i].a[j].size > max_sub_occ)
if (en - st < min_len * 2 || v[i].a[j].size > max_sub_occ || v[i].a[j].size < min_sub_occ)
continue;
t = &s[n_s++];
t->min_len = (en - st) / 2 > min_len? (en - st) / 2 : min_len;
Expand Down
0