added --idx-min-occ and --idx-max-occ

This commit is contained in:
Heng Li
2019-12-23 22:58:44 -05:00
parent d90583b83c
commit 43b0399991
4 changed files with 47 additions and 21 deletions

61
index.c
View File

@@ -1,4 +1,5 @@
#include <stdlib.h> #include <stdlib.h>
#include <limits.h>
#include <assert.h> #include <assert.h>
#if defined(WIN32) || defined(_WIN32) #if defined(WIN32) || defined(_WIN32)
#include <io.h> // for open(2) #include <io.h> // for open(2)
@@ -188,12 +189,18 @@ int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f)
* Sort and generate hash tables * * Sort and generate hash tables *
*********************************/ *********************************/
typedef struct {
mm_idx_t *mi;
int min_occ, max_occ;
} idx_post_t;
static void worker_post(void *g, long i, int tid) static void worker_post(void *g, long i, int tid)
{ {
int n, n_keys; int n, n_keys;
size_t j, start_a, start_p; size_t j, start_a, start_p;
idxhash_t *h; idxhash_t *h;
mm_idx_t *mi = (mm_idx_t*)g; idx_post_t *o = (idx_post_t*)g;
mm_idx_t *mi = o->mi;
mm_idx_bucket_t *b = &mi->B[i]; mm_idx_bucket_t *b = &mi->B[i];
if (b->a.n == 0) return; if (b->a.n == 0) return;
@@ -203,8 +210,10 @@ static void worker_post(void *g, long i, int tid)
// count and preallocate // count and preallocate
for (j = 1, n = 1, n_keys = 0, b->n = 0; j <= b->a.n; ++j) { for (j = 1, n = 1, n_keys = 0, b->n = 0; j <= b->a.n; ++j) {
if (j == b->a.n || b->a.a[j].x>>8 != b->a.a[j-1].x>>8) { if (j == b->a.n || b->a.a[j].x>>8 != b->a.a[j-1].x>>8) {
++n_keys; if (n >= o->min_occ && n <= o->max_occ) {
if (n > 1) b->n += n; ++n_keys;
if (n > 1) b->n += n;
}
n = 1; n = 1;
} else ++n; } else ++n;
} }
@@ -218,18 +227,20 @@ static void worker_post(void *g, long i, int tid)
khint_t itr; khint_t itr;
int absent; int absent;
mm128_t *p = &b->a.a[j-1]; mm128_t *p = &b->a.a[j-1];
itr = kh_put(idx, h, p->x>>8>>mi->b<<1, &absent); if (n >= o->min_occ && n <= o->max_occ) {
assert(absent && j == start_a + n); itr = kh_put(idx, h, p->x>>8>>mi->b<<1, &absent);
if (n == 1) { assert(absent && j == start_a + n);
kh_key(h, itr) |= 1; if (n == 1) {
kh_val(h, itr) = p->y; kh_key(h, itr) |= 1;
} else { kh_val(h, itr) = p->y;
int k; } else {
for (k = 0; k < n; ++k) int k;
b->p[start_p + k] = b->a.a[start_a + k].y; for (k = 0; k < n; ++k)
radix_sort_64(&b->p[start_p], &b->p[start_p + n]); // sort by position; needed as in-place radix_sort_128x() is not stable b->p[start_p + k] = b->a.a[start_a + k].y;
kh_val(h, itr) = (uint64_t)start_p<<32 | n; radix_sort_64(&b->p[start_p], &b->p[start_p + n]); // sort by position; needed as in-place radix_sort_128x() is not stable
start_p += n; kh_val(h, itr) = (uint64_t)start_p<<32 | n;
start_p += n;
}
} }
start_a = j, n = 1; start_a = j, n = 1;
} else ++n; } else ++n;
@@ -242,9 +253,12 @@ static void worker_post(void *g, long i, int tid)
b->a.n = b->a.m = 0, b->a.a = 0; b->a.n = b->a.m = 0, b->a.a = 0;
} }
static void mm_idx_post(mm_idx_t *mi, int n_threads) static void mm_idx_post(mm_idx_t *mi, int n_threads, int min_occ, int max_occ)
{ {
kt_for(n_threads, worker_post, mi, 1<<mi->b); idx_post_t t;
if (max_occ <= 0 || max_occ < min_occ) max_occ = INT_MAX;
t.mi = mi, t.min_occ = min_occ, t.max_occ = max_occ;
kt_for(n_threads, worker_post, &t, 1<<mi->b);
} }
/****************** /******************
@@ -350,7 +364,7 @@ static void *worker_pipeline(void *shared, int step, void *in)
return 0; return 0;
} }
mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int flag, int mini_batch_size, int n_threads, uint64_t batch_size) mm_idx_t *mm_idx_gen2(mm_bseq_file_t *fp, int w, int k, int b, int flag, int mini_batch_size, int n_threads, uint64_t batch_size, int min_occ, int max_occ)
{ {
pipeline_t pl; pipeline_t pl;
if (fp == 0 || mm_bseq_eof(fp)) return 0; if (fp == 0 || mm_bseq_eof(fp)) return 0;
@@ -364,13 +378,18 @@ mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int flag, int mini
if (mm_verbose >= 3) if (mm_verbose >= 3)
fprintf(stderr, "[M::%s::%.3f*%.2f] collected minimizers\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0)); fprintf(stderr, "[M::%s::%.3f*%.2f] collected minimizers\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0));
mm_idx_post(pl.mi, n_threads); mm_idx_post(pl.mi, n_threads, min_occ, max_occ);
if (mm_verbose >= 3) if (mm_verbose >= 3)
fprintf(stderr, "[M::%s::%.3f*%.2f] sorted minimizers\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0)); fprintf(stderr, "[M::%s::%.3f*%.2f] sorted minimizers\n", __func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0));
return pl.mi; return pl.mi;
} }
mm_idx_t *mm_idx_gen(mm_bseq_file_t *fp, int w, int k, int b, int flag, int mini_batch_size, int n_threads, uint64_t batch_size)
{
return mm_idx_gen2(fp, w, k, b, flag, mini_batch_size, n_threads, batch_size, 0, INT_MAX);
}
mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads) // a simpler interface; deprecated mm_idx_t *mm_idx_build(const char *fn, int w, int k, int flag, int n_threads) // a simpler interface; deprecated
{ {
mm_bseq_file_t *fp; mm_bseq_file_t *fp;
@@ -427,7 +446,7 @@ mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const cha
} }
} }
free(a.a); free(a.a);
mm_idx_post(mi, 1); mm_idx_post(mi, 1, 0, 0);
return mi; return mi;
} }
@@ -588,7 +607,7 @@ mm_idx_t *mm_idx_reader_read(mm_idx_reader_t *r, int n_threads)
if (mi && mm_verbose >= 2 && (mi->k != r->opt.k || mi->w != r->opt.w || (mi->flag&MM_I_HPC) != (r->opt.flag&MM_I_HPC))) if (mi && mm_verbose >= 2 && (mi->k != r->opt.k || mi->w != r->opt.w || (mi->flag&MM_I_HPC) != (r->opt.flag&MM_I_HPC)))
fprintf(stderr, "[WARNING]\033[1;31m Indexing parameters (-k, -w or -H) overridden by parameters used in the prebuilt index.\033[0m\n"); fprintf(stderr, "[WARNING]\033[1;31m Indexing parameters (-k, -w or -H) overridden by parameters used in the prebuilt index.\033[0m\n");
} else } else
mi = mm_idx_gen(r->fp.seq, r->opt.w, r->opt.k, r->opt.bucket_bits, r->opt.flag, r->opt.mini_batch_size, n_threads, r->opt.batch_size); mi = mm_idx_gen2(r->fp.seq, r->opt.w, r->opt.k, r->opt.bucket_bits, r->opt.flag, r->opt.mini_batch_size, n_threads, r->opt.batch_size, r->opt.min_occ, r->opt.max_occ);
if (mi) { if (mi) {
if (r->fp_out) mm_idx_dump(r->fp_out, mi); if (r->fp_out) mm_idx_dump(r->fp_out, mi);
mi->index = r->n_parts++; mi->index = r->n_parts++;

4
main.c
View File

@@ -67,6 +67,8 @@ static ko_longopt_t long_options[] = {
{ "junc-bed", ko_required_argument, 340 }, { "junc-bed", ko_required_argument, 340 },
{ "junc-bonus", ko_required_argument, 341 }, { "junc-bonus", ko_required_argument, 341 },
{ "sam-hit-only", ko_no_argument, 342 }, { "sam-hit-only", ko_no_argument, 342 },
{ "idx-min-occ", ko_required_argument, 343 },
{ "idx-max-occ", ko_required_argument, 344 },
{ "help", ko_no_argument, 'h' }, { "help", ko_no_argument, 'h' },
{ "max-intron-len", ko_required_argument, 'G' }, { "max-intron-len", ko_required_argument, 'G' },
{ "version", ko_no_argument, 'V' }, { "version", ko_no_argument, 'V' },
@@ -210,6 +212,8 @@ int main(int argc, char *argv[])
else if (c == 338) opt.max_qlen = mm_parse_num(o.arg); // --max-qlen 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) junc_bed = o.arg; // --junc-bed
else if (c == 342) opt.flag |= MM_F_SAM_HIT_ONLY; // --sam-hit-only else if (c == 342) opt.flag |= MM_F_SAM_HIT_ONLY; // --sam-hit-only
else if (c == 343) ipt.min_occ = mm_parse_num(o.arg); // --idx-min-occ
else if (c == 344) ipt.max_occ = mm_parse_num(o.arg); // --idx-max-occ
else if (c == 314) { // --frag else if (c == 314) { // --frag
yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1); yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1);
} else if (c == 315) { // --secondary } else if (c == 315) { // --secondary

View File

@@ -101,6 +101,7 @@ typedef struct {
typedef struct { typedef struct {
short k, w, flag, bucket_bits; short k, w, flag, bucket_bits;
int mini_batch_size; int mini_batch_size;
int min_occ, max_occ;
uint64_t batch_size; uint64_t batch_size;
} mm_idxopt_t; } mm_idxopt_t;

View File

@@ -1,4 +1,5 @@
#include <stdio.h> #include <stdio.h>
#include <limits.h>
#include "mmpriv.h" #include "mmpriv.h"
void mm_idxopt_init(mm_idxopt_t *opt) void mm_idxopt_init(mm_idxopt_t *opt)
@@ -8,6 +9,7 @@ void mm_idxopt_init(mm_idxopt_t *opt)
opt->bucket_bits = 14; opt->bucket_bits = 14;
opt->mini_batch_size = 50000000; opt->mini_batch_size = 50000000;
opt->batch_size = 4000000000ULL; opt->batch_size = 4000000000ULL;
opt->min_occ = 0, opt->max_occ = INT_MAX;
} }
void mm_mapopt_init(mm_mapopt_t *opt) void mm_mapopt_init(mm_mapopt_t *opt)