diff --git a/chain.c b/chain.c index a6357e2..038d1bc 100644 --- a/chain.c +++ b/chain.c @@ -19,7 +19,7 @@ static inline int ilog2_32(uint32_t v) return (t = v>>8) ? 8 + LogTable256[t] : LogTable256[v]; } -int mm_chain_dp(int max_dist, int bw, int max_skip, int min_sc, int n, mm128_t *a, uint64_t **_u, void *km) +int mm_chain_dp(int max_dist, int bw, int max_skip, int min_cnt, int min_sc, int n, mm128_t *a, uint64_t **_u, void *km) { int32_t st = 0, i, j, k, *f, *p, *t, *v, n_u, n_v; uint64_t *u; @@ -88,8 +88,10 @@ int mm_chain_dp(int max_dist, int bw, int max_skip, int min_sc, int n, mm128_t * t[j] = 1; j = p[j]; } while (j >= 0 && t[j] == 0); - if (j < 0) u[k++] = u[i]>>32<<32 | (n_v - n_v0); - else if ((int32_t)(u[i]>>32) - f[j] >= min_sc) u[k++] = ((u[i]>>32) - f[j]) << 32 | (n_v - n_v0); + if (j < 0 && n_v - n_v0 >= min_cnt) + u[k++] = u[i]>>32<<32 | (n_v - n_v0); + else if ((int32_t)(u[i]>>32) - f[j] >= min_sc && n_v - n_v0 >= min_cnt) + u[k++] = ((u[i]>>32) - f[j]) << 32 | (n_v - n_v0); else n_v = n_v0; } n_u = k, *_u = u; diff --git a/main.c b/main.c index 1d3bcbd..e5b4eaa 100644 --- a/main.c +++ b/main.c @@ -10,7 +10,7 @@ #include "minimap.h" #include "mmpriv.h" -#define MM_VERSION "2.0-r88-pre" +#define MM_VERSION "2.0-r89-pre" void liftrlimit() { @@ -46,6 +46,9 @@ static struct option long_options[] = { { "mb-size", required_argument, 0, 0 }, { "int-rname", no_argument, 0, 0 }, { "version", no_argument, 0, 'V' }, + { "min-count", required_argument, 0, 'n' }, + { "min-chain-score",required_argument, 0, 's' }, + { "mask-level", required_argument, 0, 'M' }, { 0, 0, 0, 0} }; @@ -63,7 +66,7 @@ int main(int argc, char *argv[]) mm_realtime0 = realtime(); mm_mapopt_init(&opt); - while ((c = getopt_long(argc, argv, "w:k:t:r:f:Vv:g:I:d:ST:s:x:Hcp:m:z:F:", long_options, &long_idx)) >= 0) { + while ((c = getopt_long(argc, argv, "w:k:t:r:f:Vv:g:I:d:ST:s:x:Hcp:M:n:z:F:", long_options, &long_idx)) >= 0) { if (c == 'w') w = atoi(optarg); else if (c == 'k') k = atoi(optarg); else if (c == 'H') is_hpc = 1; @@ -74,10 +77,11 @@ int main(int argc, char *argv[]) else if (c == 'v') mm_verbose = atoi(optarg); else if (c == 'g') opt.max_gap = atoi(optarg); else if (c == 'p') opt.pri_ratio = atof(optarg); - else if (c == 'm') opt.mask_level = atof(optarg); + else if (c == 'M') opt.mask_level = atof(optarg); else if (c == 'c') opt.flag |= MM_F_CIGAR; else if (c == 'S') opt.flag |= MM_F_AVA | MM_F_NO_SELF; else if (c == 'T') opt.sdust_thres = atoi(optarg); + else if (c == 'n') opt.min_cnt = atoi(optarg); else if (c == 's') opt.min_score = atoi(optarg); else if (c == 'z') opt.zdrop = atoi(optarg); else if (c == 0 && long_idx == 0) bucket_bits = atoi(optarg); // bucket-bits @@ -123,7 +127,8 @@ int main(int argc, char *argv[]) fprintf(stderr, " Mapping:\n"); fprintf(stderr, " -f FLOAT filter out top FLOAT fraction of repetitive minimizers [%.3f]\n", opt.mid_occ_frac); fprintf(stderr, " -r INT bandwidth [%d]\n", opt.bw); - fprintf(stderr, " -s INT min score [%d]\n", opt.min_score); + fprintf(stderr, " -n INT minimal number of minimizers [%d]\n", opt.min_cnt); + fprintf(stderr, " -s INT minimal chaining score (this is not SW score) [%d]\n", opt.min_score); fprintf(stderr, " -g INT split a mapping if there is a gap longer than INT [%d]\n", opt.max_gap); fprintf(stderr, " -T INT SDUST threshold; 0 to disable SDUST [%d]\n", opt.sdust_thres); fprintf(stderr, " -S skip self and dual mappings\n"); diff --git a/map.c b/map.c index c3cb1e3..98f1a4d 100644 --- a/map.c +++ b/map.c @@ -14,8 +14,9 @@ void mm_mapopt_init(mm_mapopt_t *opt) opt->mid_occ_frac = 2e-4f; opt->sdust_thres = 0; + opt->min_cnt = 3; opt->min_score = 40; - opt->bw = 1000; + opt->bw = 10000; opt->max_gap = 10000; opt->max_skip = 15; @@ -286,7 +287,7 @@ mm_reg1_t *mm_map_frag(const mm_mapopt_t *opt, const mm_idx_t *mi, mm_tbuf_t *b, kfree(b->km, m); //for (i = 0; i < n_a; ++i) printf("%c\t%s\t%d\t%d\n", "+-"[a[i].x>>63], mi->seq[a[i].x<<1>>33].name, (uint32_t)a[i].x, (uint32_t)a[i].y); - n_u = mm_chain_dp(opt->max_gap, opt->bw, opt->max_skip, opt->min_score, n_a, a, &u, b->km); + n_u = mm_chain_dp(opt->max_gap, opt->bw, opt->max_skip, opt->min_cnt, opt->min_score, n_a, a, &u, b->km); regs = mm_gen_reg(qlen, n_u, u, a); *n_regs = n_u; mm_select_sub(opt->mask_level, opt->pri_ratio, n_regs, regs, b->km); diff --git a/minimap.h b/minimap.h index 06fa758..08aa98a 100644 --- a/minimap.h +++ b/minimap.h @@ -71,7 +71,7 @@ typedef struct { int bw; // bandwidth int max_gap; // break a chain if there are no minimizers in a max_gap window int max_skip; - int min_score; + int min_cnt, min_score; float pri_ratio; float mask_level; int a, b, q, e; // matching score, mismatch, gap-open and gap-ext penalties diff --git a/mmpriv.h b/mmpriv.h index b74072b..cc56ad0 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -29,7 +29,7 @@ uint32_t ks_ksmall_uint32_t(size_t n, uint32_t arr[], size_t kk); void mm_write_paf(kstring_t *s, const mm_idx_t *mi, bseq1_t *t, int which, mm_reg1_t *r); void mm_write_sam(kstring_t *s, const mm_idx_t *mi, bseq1_t *t, int which, mm_reg1_t *r); -int mm_chain_dp(int max_dist, int bw, int max_skip, int min_sc, int n, mm128_t *a, uint64_t **_u, void *km); +int mm_chain_dp(int max_dist, int bw, int max_skip, int min_cnt, int min_sc, int n, mm128_t *a, uint64_t **_u, void *km); mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a); #ifdef __cplusplus