diff --git a/CommandLines.cpp b/CommandLines.cpp index bb828cd..5037ae6 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -102,7 +102,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->k_mer_length = 51; asm_opt->mz_win = 51; asm_opt->bf_shift = 37; - asm_opt->max_kmer_cnt = 4000; + asm_opt->max_kmer_cnt = 2000; asm_opt->high_factor = 5.0; asm_opt->max_ov_diff_ec = 0.04; asm_opt->max_ov_diff_final = 0.03; diff --git a/CommandLines.h b/CommandLines.h index cdff95e..249f3eb 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.13-r308" +#define HA_VERSION "0.13-r326-dirty" #define VERBOSE 0 diff --git a/htab.cpp b/htab.cpp index d92b097..0a9cb6e 100644 --- a/htab.cpp +++ b/htab.cpp @@ -189,7 +189,7 @@ static int ha_ct_insert_list(ha_ct_t *h, int create_new, int n, const uint64_t * ///so low 12 bits are not useful uint64_t x = a[j] >> h->pre; khint_t k; - if ((a[j]&mask) != (a[0]&mask)) continue; + assert((a[j]&mask) == (a[0]&mask)); if (create_new) { ///for 0-th counting, g->b = NULL if (g->b) @@ -365,9 +365,9 @@ int ha_pt_insert_list(ha_pt_t *h, int n, const ha_mz1_t *a) khint_t k; int n; ha_idxpos_t *p; - if ((a[j].x&mask) != (a[0].x&mask)) continue; + assert((a[j].x&mask) == (a[0].x&mask)); k = yak_pt_get(g->h, x<h)) continue; + if (k == kh_end(g->h)) continue; // TODO: understand why we sometimes come here n = kh_key(g->h, k) & YAK_MAX_COUNT; assert(n < YAK_MAX_COUNT); p = &g->a[kh_val(g->h, k) + n]; @@ -518,7 +518,7 @@ typedef struct { // global data structure for kt_pipeline() const yak_copt_t *opt; const void *flt_tab; int flag, create_new, is_store; - uint64_t n_seq; ///number of total reads + uint64_t n_mz, n_seq; ///number of total reads kseq_t *ks; UC_Read ucr; ha_ct_t *ct; @@ -573,8 +573,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip s->n_seq0 = p->n_seq; if (p->rs_in && (p->flag & HAF_RS_READ)) { while (p->n_seq < p->rs_in->total_reads) { - if((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU) - { + if ((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU) { ++p->n_seq; continue; } @@ -685,6 +684,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x); } for (i = 0; i < s->n_seq; ++i) { + p->n_mz += s->mz[i].n; free(s->mz[i].a); if (!p->is_store) free(s->seq[i]); } @@ -763,6 +763,7 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt gzclose(fp); } *n_seq = pl.n_seq; + if (pl.opt->w > 1) fprintf(stderr, "[M::%s] collected %ld minimizers\n", __func__, (long)pl.n_mz); return pl.ct; } diff --git a/sketch.cpp b/sketch.cpp index b9ff3bd..31240bf 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -25,6 +25,44 @@ static inline int tq_shift(tiny_queue_t *q) return x; } +static inline int mzcmp(const ha_mz1_t *a, const ha_mz1_t *b) +{ + return a->rid < b->rid? -1 : a->rid > b->rid? 1 : ((a->x > b->x) - (a->x < b->x)); +} + +static void select_mz(ha_mz1_v *p, int len) +{ + static const ha_mz1_t dummy = { UINT64_MAX, (1<<28) - 1, 0, 0 }; + int32_t i, last0 = -1, n = (int32_t)p->n, m = 0; + if (n == 0 || n == 1) return; + for (i = 0; i < n; ++i) + if (p->a[i].rid != 0) ++m; + if (m == 0) return; // no high-frequency k-mers; do nothing + for (i = 0; i <= n; ++i) { + if (i == n || p->a[i].rid == 0) { + if (i - last0 > 1) { + int32_t ps = last0 < 0? 0 : p->a[last0].pos; + int32_t pe = i == n? len : p->a[i].pos; + int32_t j, st = last0 + 1, en = i; + ha_mz1_t min1 = dummy, min2 = dummy; + int32_t min1_i = -1, min2_i = -1; + for (j = st; j < en; ++j) { // choose up to two minimum k-mers + if (mzcmp(&p->a[j], &min1) < 0) min2 = min1, min2_i = min1_i, min1 = p->a[j], min1_i = j; + else if (mzcmp(&p->a[j], &min2) < 0) min2 = p->a[j], min2_i = j; + } + if (min1_i >= 0 && p->a[min1_i].rid < pe - ps) p->a[min1_i].rid = 0; + if (min2_i >= 0 && p->a[min2_i].rid < pe - ps) p->a[min2_i].rid = 0; + } + last0 = i; + } + } + for (i = n = 0; i < (int32_t)p->n; ++i) // squeeze out filtered minimizers + if (p->a[i].rid == 0) + p->a[n++] = p->a[i]; +// fprintf(stderr, "X\tn0=%d,n1=%d,m=%d\n", p->n, n, m); + p->n = n; +} + /** * Find symmetric (w,k)-minimizers on a DNA sequence * @@ -43,7 +81,7 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, uint64_t rid:28, pos:27, rev:1, span:8; **/ extern void *ha_ct_table; - static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 }; + static const ha_mz1_t dummy = { UINT64_MAX, (1<<28) - 1, 0, 0 }; uint64_t shift1 = k - 1, mask = (1ULL<= k && kmer_span < 256) { uint64_t y; - int32_t cnt = 0, filtered = 0; + int32_t cnt, filtered; y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]); - if (hf != 0) cnt = ha_ft_cnt(hf, y); - filtered = (cnt > 0); + cnt = hf? ha_ft_cnt(hf, y) : 0; + filtered = (cnt >= 1<<28); if (dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i)); - if (cnt < 1<<28 && filtered == 0) - info.x = y, info.rid = cnt, info.pos = i, info.rev = z, info.span = kmer_span; + if (!filtered) info.x = y, info.rid = cnt, info.pos = i, info.rev = z, info.span = kmer_span; if (k_flag != NULL) k_flag->a.a[i]++; if (k_flag != NULL && filtered > 0) k_flag->a.a[i]++; } @@ -115,9 +152,9 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet for (j = buf_pos + 1; j < w; ++j) - if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]); + if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]); for (j = 0; j < buf_pos; ++j) - if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]); + if (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]); } /** * There are three cases: @@ -135,21 +172,22 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ///buf_pos == min_pos, means current minimizer has moved outside the window ///so for now we need to find a new minimizer at the current window (w k-mers) for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers - if (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer + if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer for (j = 0; j <= buf_pos; ++j) - if (min.x >= buf[j].x) min = buf[j], min_pos = j; + if (mzcmp(&min, &buf[j]) >= 0) min = buf[j], min_pos = j; if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted - if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]); + if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]); for (j = 0; j <= buf_pos; ++j) - if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]); + if (mzcmp(&min, &buf[j]) == 0 && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]); } } if (++buf_pos == w) buf_pos = 0; } if (min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); + select_mz(p, len); for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts p->a[i].rid = rid; }