zzh
/
minimap2
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

r89: added minimal number of minimizer counts

This commit is contained in:
Heng Li 2017-06-27 18:43:15 -04:00
parent 8977f07269
commit fa80177e58
5 changed files with 19 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
chain.c
View File

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

13
main.c
View File

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

5
map.c
View File

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

2
minimap.h
View File

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

2
mmpriv.h
View File

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