From ef8c04755920df0dd7b2d06128b1926ea048861c Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 29 Oct 2020 11:28:43 -0400 Subject: [PATCH 01/10] removed unused ha_sketch --- Overlaps.cpp | 5 +- anchor.cpp | 9 ++-- htab.cpp | 20 +++++-- htab.h | 2 +- sketch.cpp | 148 +++++++-------------------------------------------- 5 files changed, 43 insertions(+), 141 deletions(-) diff --git a/Overlaps.cpp b/Overlaps.cpp index f86ac7c..c026767 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -9974,8 +9974,7 @@ void clean_weak_ma_hit_t(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_source -void debug_info_of_specfic_read(const char* name, ma_hit_t_alloc* sources, -ma_hit_t_alloc* reverse_sources, int id, const char* command) +void debug_info_of_specfic_read(const char* name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int id, const char* command) { long long i, j, Len; uint32_t tn; @@ -9999,7 +9998,7 @@ ma_hit_t_alloc* reverse_sources, int id, const char* command) fprintf(stderr, "\n\n\nafter %s\n", command); fprintf(stderr, "****************ma_hit_t (%lld)ref_read: %.*s, len: %lu****************\n", - i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), Get_READ_LENGTH(R_INF, i)); + i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i), (unsigned long)Get_READ_LENGTH(R_INF, i)); fprintf(stderr, "sources Len: %d, is_fully_corrected: %d\n", diff --git a/anchor.cpp b/anchor.cpp index 50a644f..ddf2589 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -57,8 +57,8 @@ int ha_ov_type(const overlap_region *r, uint32_t len) else return r->x_pos_s == 0? 0 : 1; } -void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, -kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct) +void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain, + kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct) { uint32_t i, rlen; uint64_t k, l; @@ -205,7 +205,8 @@ void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf) void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, -int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct) + int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar, + kvec_t_u64_warp* dbg_ct) { extern void *ha_flt_tab; extern ha_pt_t *ha_idx; @@ -302,4 +303,4 @@ int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list) { ks_introsort_or_xs(overlap_list->length, overlap_list->list); -} \ No newline at end of file +} diff --git a/htab.cpp b/htab.cpp index 7082582..2194d23 100644 --- a/htab.cpp +++ b/htab.cpp @@ -557,7 +557,7 @@ static void worker_for_mz(void *data, long i, int tid) ha_mz1_v *b = &s->mz_buf[tid]; s->mz_buf[tid].n = 0; ///s->p->opt->w = 51, s->p->opt->k - ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab); + ha_sketch_query(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, 0, 0); s->mz[i].n = s->mz[i].m = b->n; MALLOC(s->mz[i].a, b->n); memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t)); @@ -803,7 +803,7 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi * High count filter table * ***************************/ -KHASHL_SET_INIT(static klib_unused, yak_ft_t, yak_ft, uint64_t, kh_hash_dummy, kh_eq_generic) +KHASHL_MAP_INIT(static klib_unused, yak_ft_t, yak_ft, uint64_t, int16_t, kh_hash_dummy, kh_eq_generic) static yak_ft_t *gen_hh(const ha_ct_t *h) { @@ -813,12 +813,14 @@ static yak_ft_t *gen_hh(const ha_ct_t *h) yak_ft_resize(hh, h->tot * 2); for (i = 0; i < 1<pre; ++i) { yak_ct_t *ht = h->h[i].h; - khint_t k; + khint_t k, l; for (k = 0; k < kh_end(ht); ++k) { if (kh_exist(ht, k)) { uint64_t y = kh_key(ht, k) >> h->pre << YAK_COUNTER_BITS | i; int absent; - yak_ft_put(hh, y, &absent); + l = yak_ft_put(hh, y, &absent); + if (absent) + kh_val(hh, l) = kh_key(ht, k)&YAK_MAX_COUNT; } } } @@ -833,6 +835,14 @@ int ha_ft_isflt(const void *hh, uint64_t y) return k == kh_end(h)? 0 : 1; } +int ha_ft_cnt(const void *hh, uint64_t y) +{ + yak_ft_t *h = (yak_ft_t*)hh; + khint_t k; + k = yak_ft_get(h, y); + return k == kh_end(h)? 0 : kh_val(h, k); +} + void ha_ft_destroy(void *h) { if (h) yak_ft_destroy((yak_ft_t*)h); @@ -1092,7 +1102,7 @@ int load_pt_index(void **r_flt_tab, ha_pt_t **r_ha_idx, All_reads* r, hifiasm_op ha_pt_t *ha_idx = NULL; char mode = 0; - int f_flag, absent, i; + int f_flag = 0, absent, i; double index_time, index_s_time, pos_time, pos_s_time; diff --git a/htab.h b/htab.h index bb53bb6..c703cb8 100644 --- a/htab.h +++ b/htab.h @@ -33,7 +33,7 @@ extern void *ha_ct_table; void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode); -int ha_ft_isflt(const void *hh, uint64_t y); +int ha_ft_cnt(const void *hh, uint64_t y); void ha_ft_destroy(void *h); ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov); diff --git a/sketch.cpp b/sketch.cpp index 9bd31ec..2fe0b22 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -36,12 +36,13 @@ static inline int tq_shift(tiny_queue_t *q) * @param is_hpc homopolymer-compressed or not * @param p minimizers */ -void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf) +void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct) { ///in default, w = 51, k = 51, is_hpc = 1 /** uint64_t x; 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 }; uint64_t shift1 = k - 1, mask = (1ULL< 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63)); + if (dbg_ct != NULL) dbg_ct->a.n = 0; + if (k_flag != NULL) { + kv_resize(uint8_t, k_flag->a, (uint64_t)len); + k_flag->a.n = len; + memset(k_flag->a.a, 0, k_flag->a.n); + } + ///sizeof(ha_mz1_t) = 16 memset(buf, 0xff, w * 16); memset(&tq, 0, sizeof(tiny_queue_t)); @@ -70,12 +78,16 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, } tq_push(&tq, skip_len); kmer_span += skip_len; + ///how many bases that are covered by this HPC k-mer + ///kmer_span includes at most k HPC elements if (tq.count > k) kmer_span -= tq_shift(&tq); } else kmer_span = l + 1 < k? l + 1 : k; ///kmer_span should be used for HPC k-mer - ///so for non-HPC k-mer, kmer_span should be k in any case? + ///non-HPC k-mer, kmer_span should be k ///kmer_span is used to calculate anchor pos on reverse complementary strand + if (k_flag != NULL) k_flag->a.a[i] = 1;///lable all useful base, which are not ignored by HPC + kmer[0] = (kmer[0] << 1 | (c&1)) & mask; // forward k-mer kmer[1] = (kmer[1] << 1 | (c>>1)) & mask; kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; // reverse k-mer @@ -85,13 +97,17 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ++l; if (l >= k && kmer_span < 256) { uint64_t y; + int filtered = 0; y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]); - if (hf == 0 || ha_ft_isflt(hf, y) == 0) + if (hf != 0) filtered = (ha_ft_cnt(hf, y) > 0); + 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 (filtered == 0) info.x = y, info.rid = rid, 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]++; } } else l = 0, tq.count = tq.front = 0, kmer_span = 0; - //for non-HPC k-mer, l = i; but for HPC k-mer, l is always less than i //i is the real base iterator, while l is the HPC base iterator //only if l >= k, info is a useful minimizer (ha_mz1_t.x != UINT64_MAX) @@ -135,127 +151,3 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, if (min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); } - - - -void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, -kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct) -{ ///in default, w = 51, k = 51, is_hpc = 1 - /** - uint64_t x; - uint64_t rid:28, pos:27, rev:1, span:8; - **/ - extern void *ha_ct_table; - if(dbg_ct != NULL) dbg_ct->a.n = 0; - - static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 }; - uint64_t shift1 = k - 1, mask = (1ULL<a, (uint64_t)len); - k_flag->a.n = len; - memset(k_flag->a.a, 0, k_flag->a.n); - } - - - assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63)); - ///sizeof(ha_mz1_t) = 16 - memset(buf, 0xff, w * 16); - memset(&tq, 0, sizeof(tiny_queue_t)); - ///len/w is the evaluated minimizer numbers - kv_resize(ha_mz1_t, *p, p->n + len/w); - - for (i = l = buf_pos = min_pos = 0; i < len; ++i) { - int c = seq_nt4_table[(uint8_t)str[i]]; - ha_mz1_t info = dummy; - if (c < 4) { // not an ambiguous base - int z; - if (is_hpc) { - int skip_len = 1; - if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) { - for (skip_len = 2; i + skip_len < len; ++skip_len) - if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c) - break; - i += skip_len - 1; // put $i at the end of the current homopolymer run - } - tq_push(&tq, skip_len); - kmer_span += skip_len; - ///how many bases that are covered by this HPC k-mer - ///kmer_span includes at most k HPC elements - if (tq.count > k) kmer_span -= tq_shift(&tq); - } else kmer_span = l + 1 < k? l + 1 : k; - ///kmer_span should be used for HPC k-mer - ///non-HPC k-mer, kmer_span should be k - ///kmer_span is used to calculate anchor pos on reverse complementary strand - - if(k_flag != NULL) k_flag->a.a[i] = 1;///lable all useful base, which are not ignored by HPC - - kmer[0] = (kmer[0] << 1 | (c&1)) & mask; // forward k-mer - kmer[1] = (kmer[1] << 1 | (c>>1)) & mask; - kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; // reverse k-mer - kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1; - if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand - z = kmer[1] < kmer[3]? 0 : 1; // strand - ++l; - if (l >= k && kmer_span < 256) { - uint64_t y; - y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]); - - filtered = 0; - if(hf != 0) filtered = ha_ft_isflt(hf, y); - 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 (hf == 0 || ha_ft_isflt(hf, y) == 0) - if(filtered == 0) - info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span; - if(k_flag != NULL) k_flag->a.a[i]++; - if(k_flag != NULL && filtered == 1) k_flag->a.a[i]++; - } - } else l = 0, tq.count = tq.front = 0, kmer_span = 0; - - - //for non-HPC k-mer, l = i; but for HPC k-mer, l is always less than i - //i is the real base iterator, while l is the HPC base iterator - //only if l >= k, info is a useful minimizer (ha_mz1_t.x != UINT64_MAX) - //but even if l < k, infor is still stored into buf - 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]); - 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]); - } - /** - * There are three cases: - * 1. info.x <= min.x, means info is a new minimizer - * 2. info.x > min.x, info is not a new minimizer - * (1) buf_pos != min_pos, do nothing - * (2) buf_pos == min_pos, means current minimizer has moved outside the window - * **/ - ///three cases: 1. - if (info.x <= min.x) { // a new minimum; then write the old min - if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); - min = info, min_pos = buf_pos; - } else if (buf_pos == min_pos) { // old min has moved outside the window - if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); - ///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 - for (j = 0; j <= buf_pos; ++j) - if (min.x >= buf[j].x) 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]); - 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 (++buf_pos == w) buf_pos = 0; - } - if (min.x != UINT64_MAX) - kv_push(ha_mz1_t, *p, min); -} \ No newline at end of file From f37c4e29fa00dac964c7dbf010bda1e19561145e Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 29 Oct 2020 18:24:47 -0400 Subject: [PATCH 02/10] a bit refactoring for the next change --- CommandLines.cpp | 1 + CommandLines.h | 1 + anchor.cpp | 2 +- htab.cpp | 29 ++++++++++++----------------- htab.h | 6 ++---- sketch.cpp | 16 +++++++++------- 6 files changed, 26 insertions(+), 29 deletions(-) diff --git a/CommandLines.cpp b/CommandLines.cpp index 2d27297..bb828cd 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -102,6 +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->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 86b351b..cdff95e 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -35,6 +35,7 @@ typedef struct { int k_mer_length; int mz_win; int bf_shift; + int max_kmer_cnt; double high_factor; // coverage cutoff set to high_factor*hom_cov double max_ov_diff_ec; double max_ov_diff_final; diff --git a/anchor.cpp b/anchor.cpp index ddf2589..6888be0 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -73,7 +73,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg rlen = Get_READ_LENGTH(R_INF, rid); // read length // get the list of anchors - ha_sketch_query(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag, dbg_ct); + ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag, dbg_ct); // minimizer of queried read if (ab->mz.m > ab->old_mz_m) { ab->old_mz_m = ab->mz.m; diff --git a/htab.cpp b/htab.cpp index 2194d23..d92b097 100644 --- a/htab.cpp +++ b/htab.cpp @@ -1,4 +1,3 @@ -#include #include #include #include @@ -557,7 +556,7 @@ static void worker_for_mz(void *data, long i, int tid) ha_mz1_v *b = &s->mz_buf[tid]; s->mz_buf[tid].n = 0; ///s->p->opt->w = 51, s->p->opt->k - ha_sketch_query(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, 0, 0); + ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, 0, 0); s->mz[i].n = s->mz[i].m = b->n; MALLOC(s->mz[i].a, b->n); memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t)); @@ -803,12 +802,15 @@ ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const voi * High count filter table * ***************************/ +// Warning: the max count is 32767 KHASHL_MAP_INIT(static klib_unused, yak_ft_t, yak_ft, uint64_t, int16_t, kh_hash_dummy, kh_eq_generic) -static yak_ft_t *gen_hh(const ha_ct_t *h) +static yak_ft_t *gen_hh(const ha_ct_t *h, int max_cnt) { int i; yak_ft_t *hh; + if (max_cnt > YAK_MAX_COUNT - 1) max_cnt = YAK_MAX_COUNT - 1; + if (max_cnt > INT16_MAX - 1) max_cnt = INT16_MAX - 1; hh = yak_ft_init(); yak_ft_resize(hh, h->tot * 2); for (i = 0; i < 1<pre; ++i) { @@ -819,28 +821,22 @@ static yak_ft_t *gen_hh(const ha_ct_t *h) uint64_t y = kh_key(ht, k) >> h->pre << YAK_COUNTER_BITS | i; int absent; l = yak_ft_put(hh, y, &absent); - if (absent) - kh_val(hh, l) = kh_key(ht, k)&YAK_MAX_COUNT; + if (absent) { + int cnt = kh_key(ht, k) & YAK_MAX_COUNT; + kh_val(hh, l) = cnt > max_cnt? INT16_MAX : cnt; + } } } } return hh; } -int ha_ft_isflt(const void *hh, uint64_t y) +int32_t ha_ft_cnt(const void *hh, uint64_t y) { yak_ft_t *h = (yak_ft_t*)hh; khint_t k; k = yak_ft_get(h, y); - return k == kh_end(h)? 0 : 1; -} - -int ha_ft_cnt(const void *hh, uint64_t y) -{ - yak_ft_t *h = (yak_ft_t*)hh; - khint_t k; - k = yak_ft_get(h, y); - return k == kh_end(h)? 0 : kh_val(h, k); + return k == kh_end(h)? 0 : kh_val(h, k) == INT16_MAX? INT32_MAX : kh_val(h, k); } void ha_ft_destroy(void *h) @@ -906,7 +902,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int i if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1; } ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num); - flt_tab = gen_hh(h); + flt_tab = gen_hh(h, asm_opt->max_kmer_cnt); ha_ct_destroy(h); fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> filtered out %ld k-mers occurring %d or more times\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb(), (long)kh_size(flt_tab), cutoff); @@ -970,7 +966,6 @@ int query_ct_index(void* ct_idx, uint64_t hash) return kh_key(g->h, k)&YAK_MAX_COUNT; } - int write_ct_index(void *i_ct_idx, char* file_name) { char* gfa_name = (char*)malloc(strlen(file_name)+25); diff --git a/htab.h b/htab.h index c703cb8..4a8de45 100644 --- a/htab.h +++ b/htab.h @@ -31,9 +31,8 @@ extern void *ha_flt_tab_hp; extern ha_pt_t *ha_idx_hp; extern void *ha_ct_table; - void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode); -int ha_ft_cnt(const void *hh, uint64_t y); +int32_t ha_ft_cnt(const void *hh, uint64_t y); void ha_ft_destroy(void *h); ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, int is_hp_mode, All_reads *rs, int *hom_cov, int *het_cov); @@ -59,8 +58,7 @@ double yak_cpu_usage(void); void ha_triobin(const hifiasm_opt_t *opt); -void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf); -void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct); +void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct); int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het); static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function diff --git a/sketch.cpp b/sketch.cpp index 2fe0b22..b9ff3bd 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -36,7 +36,7 @@ static inline int tq_shift(tiny_queue_t *q) * @param is_hpc homopolymer-compressed or not * @param p minimizers */ -void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct) +void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct) { ///in default, w = 51, k = 51, is_hpc = 1 /** uint64_t x; @@ -57,8 +57,7 @@ void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int i memset(k_flag->a.a, 0, k_flag->a.n); } - ///sizeof(ha_mz1_t) = 16 - memset(buf, 0xff, w * 16); + memset(buf, 0xff, w * sizeof(ha_mz1_t)); memset(&tq, 0, sizeof(tiny_queue_t)); ///len/w is the evaluated minimizer numbers kv_resize(ha_mz1_t, *p, p->n + len/w); @@ -97,12 +96,13 @@ void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int i ++l; if (l >= k && kmer_span < 256) { uint64_t y; - int filtered = 0; + int32_t cnt = 0, filtered = 0; y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]); - if (hf != 0) filtered = (ha_ft_cnt(hf, y) > 0); + if (hf != 0) cnt = ha_ft_cnt(hf, y); + filtered = (cnt > 0); 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 (filtered == 0) - info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span; + if (cnt < 1<<28 && filtered == 0) + 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]++; } @@ -150,4 +150,6 @@ void ha_sketch_query(const char *str, int len, int w, int k, uint32_t rid, int i } if (min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); + for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts + p->a[i].rid = rid; } From f02eb66e6a0774c34b2d92b30177a1deda3c30c8 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 29 Oct 2020 22:36:54 -0400 Subject: [PATCH 03/10] r326: select high-occ k-mers (experimental) --- CommandLines.cpp | 2 +- CommandLines.h | 2 +- htab.cpp | 13 +++++----- sketch.cpp | 62 ++++++++++++++++++++++++++++++++++++++---------- 4 files changed, 59 insertions(+), 20 deletions(-) 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; } From 6fa4d44426959fd813c8dbbdc0f30230cc5d5a12 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sun, 1 Nov 2020 10:44:02 -0500 Subject: [PATCH 04/10] r327: choose up to 3 high-freq minimizers --- CommandLines.h | 2 +- sketch.cpp | 24 +++++++++++++++++------- 2 files changed, 18 insertions(+), 8 deletions(-) diff --git a/CommandLines.h b/CommandLines.h index 249f3eb..53b27cf 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.13-r326-dirty" +#define HA_VERSION "0.13-r327-dirty" #define VERBOSE 0 diff --git a/sketch.cpp b/sketch.cpp index 31240bf..1a7b7ba 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -44,14 +44,24 @@ static void select_mz(ha_mz1_v *p, int len) 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; + ha_mz1_t min1 = dummy, min2 = dummy, min3 = dummy; + int32_t min1_i = -1, min2_i = -1, min3_i = -1; + for (j = st; j < en; ++j) { // choose up to three minimum k-mers + if (mzcmp(&p->a[j], &min1) < 0) + min3 = min2, min3_i = min2_i, min2 = min1, min2_i = min1_i, min1 = p->a[j], min1_i = j; + else if (mzcmp(&p->a[j], &min2) < 0) + min3 = min2, min3_i = min2_i, min2 = p->a[j], min2_i = j; + else if (mzcmp(&p->a[j], &min3) < 0) + min3 = p->a[j], min3_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; + if (min3_i >= 0 && p->a[min3_i].rid < pe - ps) + p->a[min3_i].rid = 0; + } } - 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; } From dcb583bd697c5eed08472fe9d409ce95dc37905e Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sun, 1 Nov 2020 22:49:06 -0500 Subject: [PATCH 05/10] r328: use binomial heap for general cases --- CommandLines.h | 2 +- ksort.h | 19 +++++++++++++++++- sketch.cpp | 52 +++++++++++++++++++++++++++++--------------------- 3 files changed, 49 insertions(+), 24 deletions(-) diff --git a/CommandLines.h b/CommandLines.h index 53b27cf..d631ffb 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.13-r327-dirty" +#define HA_VERSION "0.13-r328-dirty" #define VERBOSE 0 diff --git a/ksort.h b/ksort.h index 0d3d82c..1736abc 100644 --- a/ksort.h +++ b/ksort.h @@ -40,7 +40,24 @@ typedef struct { #define KSORT_SWAP(type_t, a, b) { register type_t t=(a); (a)=(b); (b)=t; } -#define KSORT_INIT(name, type_t, __sort_lt) \ +#define KSORT_INIT(name, type_t, __sort_lt) \ + void ks_heapdown_##name(size_t i, size_t n, type_t l[]) \ + { \ + size_t k = i; \ + type_t tmp = l[i]; \ + while ((k = (k << 1) + 1) < n) { \ + if (k != n - 1 && __sort_lt(l[k], l[k+1])) ++k; \ + if (__sort_lt(l[k], tmp)) break; \ + l[i] = l[k]; i = k; \ + } \ + l[i] = tmp; \ + } \ + void ks_heapmake_##name(size_t lsize, type_t l[]) \ + { \ + size_t i; \ + for (i = (lsize >> 1) - 1; i != (size_t)(-1); --i) \ + ks_heapdown_##name(i, lsize, l); \ + } \ static inline void __ks_insertsort_##name(type_t *s, type_t *t) \ { \ type_t *i, *j, swap_tmp; \ diff --git a/sketch.cpp b/sketch.cpp index 1a7b7ba..a9d8594 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -4,6 +4,10 @@ #include #include "kvec.h" #include "htab.h" +#include "ksort.h" + +#define MAX_HIGH_OCC 3 // TODO: don't hard code if we need to tune this parameter +#define MAX_MAX_HIGH_OCC 16 typedef struct { // a simplified version of kdq int front, count; @@ -30,11 +34,18 @@ 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 }; +#define mz_lt(a, b) (mzcmp(&(a), &(b)) < 0) +KSORT_INIT(mz, ha_mz1_t, mz_lt) + +static void select_mz(ha_mz1_v *p, int len, int max_high_occ) +{ // for high-occ minimizers, choose up to max_high_occ in each high-occ streak int32_t i, last0 = -1, n = (int32_t)p->n, m = 0; + ha_mz1_t b[MAX_MAX_HIGH_OCC]; // this is to avoid a heap allocation + if (n == 0 || n == 1) return; + assert(n < 1<<27); // 27 is the number of bits for ha_mz1_t::pos; this should be safe as there are more bases than minimizers + if (max_high_occ > MAX_MAX_HIGH_OCC) + max_high_occ = MAX_MAX_HIGH_OCC; for (i = 0; i < n; ++i) if (p->a[i].rid != 0) ++m; if (m == 0) return; // no high-frequency k-mers; do nothing @@ -43,25 +54,22 @@ static void select_mz(ha_mz1_v *p, int len) 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, min3 = dummy; - int32_t min1_i = -1, min2_i = -1, min3_i = -1; - for (j = st; j < en; ++j) { // choose up to three minimum k-mers - if (mzcmp(&p->a[j], &min1) < 0) - min3 = min2, min3_i = min2_i, min2 = min1, min2_i = min1_i, min1 = p->a[j], min1_i = j; - else if (mzcmp(&p->a[j], &min2) < 0) - min3 = min2, min3_i = min2_i, min2 = p->a[j], min2_i = j; - else if (mzcmp(&p->a[j], &min3) < 0) - min3 = p->a[j], min3_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; - if (min3_i >= 0 && p->a[min3_i].rid < pe - ps) - p->a[min3_i].rid = 0; + int32_t j, k, st = last0 + 1, en = i; + for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k) + b[k] = p->a[j], b[k].pos = j; // b[].pos keeps the index in p->a[] + if (j < en) { // if there are more, choose top max_high_occ + assert(k == max_high_occ); + ks_heapmake_mz(max_high_occ, b); // initialize the binomial heap + for (; j < en; ++j) { + if (mz_lt(p->a[j], b[0])) { // then update the heap + b[0] = p->a[j], b[0].pos = j; + ks_heapdown_mz(0, max_high_occ, b); + } } } + for (j = 0; j < k; ++j) + if (b[j].rid < pe - ps) + p->a[b[j].pos].rid = 0; } last0 = i; } @@ -149,7 +157,7 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, 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 (!filtered) 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; // initially ha_mz1_t::rid keeps the k-mer count if (k_flag != NULL) k_flag->a.a[i]++; if (k_flag != NULL && filtered > 0) k_flag->a.a[i]++; } @@ -197,7 +205,7 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, } if (min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); - select_mz(p, len); + select_mz(p, len, MAX_HIGH_OCC); for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts p->a[i].rid = rid; } From 02f0b5dd22cea9f82e7cfa3afcfead57f346fe1a Mon Sep 17 00:00:00 2001 From: Heng Li Date: Tue, 3 Nov 2020 00:44:15 -0500 Subject: [PATCH 06/10] r329: a bit code cleanup --- CommandLines.h | 2 +- ksort.h | 8 ++++++++ sketch.cpp | 15 ++++++--------- 3 files changed, 15 insertions(+), 10 deletions(-) diff --git a/CommandLines.h b/CommandLines.h index d631ffb..454b07e 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.13-r328-dirty" +#define HA_VERSION "0.13-r329-dirty" #define VERBOSE 0 diff --git a/ksort.h b/ksort.h index 1736abc..67836fb 100644 --- a/ksort.h +++ b/ksort.h @@ -58,6 +58,14 @@ typedef struct { for (i = (lsize >> 1) - 1; i != (size_t)(-1); --i) \ ks_heapdown_##name(i, lsize, l); \ } \ + void ks_heapsort_##name(size_t lsize, type_t l[]) \ + { \ + size_t i; \ + for (i = lsize - 1; i > 0; --i) { \ + type_t tmp; \ + tmp = *l; *l = l[i]; l[i] = tmp; ks_heapdown_##name(0, i, l); \ + } \ + } \ static inline void __ks_insertsort_##name(type_t *s, type_t *t) \ { \ type_t *i, *j, swap_tmp; \ diff --git a/sketch.cpp b/sketch.cpp index a9d8594..5e40680 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -57,16 +57,14 @@ static void select_mz(ha_mz1_v *p, int len, int max_high_occ) int32_t j, k, st = last0 + 1, en = i; for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k) b[k] = p->a[j], b[k].pos = j; // b[].pos keeps the index in p->a[] - if (j < en) { // if there are more, choose top max_high_occ - assert(k == max_high_occ); - ks_heapmake_mz(max_high_occ, b); // initialize the binomial heap - for (; j < en; ++j) { - if (mz_lt(p->a[j], b[0])) { // then update the heap - b[0] = p->a[j], b[0].pos = j; - ks_heapdown_mz(0, max_high_occ, b); - } + ks_heapmake_mz(k, b); // initialize the binomial heap + for (; j < en; ++j) { // if there are more, choose top max_high_occ + if (mz_lt(p->a[j], b[0])) { // then update the heap + b[0] = p->a[j], b[0].pos = j; + ks_heapdown_mz(0, k, b); } } + //ks_heapsort_mz(k, b); // sorting is not needed for now for (j = 0; j < k; ++j) if (b[j].rid < pe - ps) p->a[b[j].pos].rid = 0; @@ -77,7 +75,6 @@ static void select_mz(ha_mz1_v *p, int len, int max_high_occ) 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; } From da672118e54a37567be973e1b795887e9ce29c30 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 21 Jan 2021 11:05:37 -0500 Subject: [PATCH 07/10] r330: unimap heuristic for mz sampling --- CommandLines.cpp | 1 + CommandLines.h | 3 ++- anchor.cpp | 2 +- htab.cpp | 2 +- htab.h | 2 +- sketch.cpp | 16 ++++++++-------- 6 files changed, 14 insertions(+), 12 deletions(-) diff --git a/CommandLines.cpp b/CommandLines.cpp index 5037ae6..3568aa7 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -101,6 +101,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->thread_num = 1; asm_opt->k_mer_length = 51; asm_opt->mz_win = 51; + asm_opt->sample_dist = 200; asm_opt->bf_shift = 37; asm_opt->max_kmer_cnt = 2000; asm_opt->high_factor = 5.0; diff --git a/CommandLines.h b/CommandLines.h index 454b07e..082e141 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.13-r329-dirty" +#define HA_VERSION "0.13-r330-dirty" #define VERBOSE 0 @@ -34,6 +34,7 @@ typedef struct { int thread_num; int k_mer_length; int mz_win; + int sample_dist; int bf_shift; int max_kmer_cnt; double high_factor; // coverage cutoff set to high_factor*hom_cov diff --git a/anchor.cpp b/anchor.cpp index 6888be0..9f0e83d 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -73,7 +73,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg rlen = Get_READ_LENGTH(R_INF, rid); // read length // get the list of anchors - ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, k_flag, dbg_ct); + ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.sample_dist, k_flag, dbg_ct); // minimizer of queried read if (ab->mz.m > ab->old_mz_m) { ab->old_mz_m = ab->mz.m; diff --git a/htab.cpp b/htab.cpp index 0a9cb6e..936446c 100644 --- a/htab.cpp +++ b/htab.cpp @@ -556,7 +556,7 @@ static void worker_for_mz(void *data, long i, int tid) ha_mz1_v *b = &s->mz_buf[tid]; s->mz_buf[tid].n = 0; ///s->p->opt->w = 51, s->p->opt->k - ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, 0, 0); + ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, 0, 0, 0); s->mz[i].n = s->mz[i].m = b->n; MALLOC(s->mz[i].a, b->n); memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t)); diff --git a/htab.h b/htab.h index 4a8de45..9a801a3 100644 --- a/htab.h +++ b/htab.h @@ -58,7 +58,7 @@ double yak_cpu_usage(void); void ha_triobin(const hifiasm_opt_t *opt); -void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct); +void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct); int ha_analyze_count(int n_cnt, int start_cnt, const int64_t *cnt, int *peak_het); static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function diff --git a/sketch.cpp b/sketch.cpp index 5e40680..f053073 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -6,7 +6,7 @@ #include "htab.h" #include "ksort.h" -#define MAX_HIGH_OCC 3 // TODO: don't hard code if we need to tune this parameter +#define MAX_HIGH_OCC 8 // TODO: don't hard code if we need to tune this parameter #define MAX_MAX_HIGH_OCC 16 typedef struct { // a simplified version of kdq @@ -37,15 +37,13 @@ static inline int mzcmp(const ha_mz1_t *a, const ha_mz1_t *b) #define mz_lt(a, b) (mzcmp(&(a), &(b)) < 0) KSORT_INIT(mz, ha_mz1_t, mz_lt) -static void select_mz(ha_mz1_v *p, int len, int max_high_occ) +static void select_mz(ha_mz1_v *p, int len, int sample_dist) { // for high-occ minimizers, choose up to max_high_occ in each high-occ streak int32_t i, last0 = -1, n = (int32_t)p->n, m = 0; ha_mz1_t b[MAX_MAX_HIGH_OCC]; // this is to avoid a heap allocation if (n == 0 || n == 1) return; assert(n < 1<<27); // 27 is the number of bits for ha_mz1_t::pos; this should be safe as there are more bases than minimizers - if (max_high_occ > MAX_MAX_HIGH_OCC) - max_high_occ = MAX_MAX_HIGH_OCC; for (i = 0; i < n; ++i) if (p->a[i].rid != 0) ++m; if (m == 0) return; // no high-frequency k-mers; do nothing @@ -55,6 +53,9 @@ static void select_mz(ha_mz1_v *p, int len, int max_high_occ) int32_t ps = last0 < 0? 0 : p->a[last0].pos; int32_t pe = i == n? len : p->a[i].pos; int32_t j, k, st = last0 + 1, en = i; + int32_t max_high_occ = (int32_t)((double)(pe - ps) / sample_dist + .499); + if (max_high_occ > MAX_MAX_HIGH_OCC) + max_high_occ = MAX_MAX_HIGH_OCC; for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k) b[k] = p->a[j], b[k].pos = j; // b[].pos keeps the index in p->a[] ks_heapmake_mz(k, b); // initialize the binomial heap @@ -66,8 +67,7 @@ static void select_mz(ha_mz1_v *p, int len, int max_high_occ) } //ks_heapsort_mz(k, b); // sorting is not needed for now for (j = 0; j < k; ++j) - if (b[j].rid < pe - ps) - p->a[b[j].pos].rid = 0; + p->a[b[j].pos].rid = 0; } last0 = i; } @@ -89,7 +89,7 @@ static void select_mz(ha_mz1_v *p, int len, int max_high_occ) * @param is_hpc homopolymer-compressed or not * @param p minimizers */ -void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct) +void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf, int sample_dist, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct) { ///in default, w = 51, k = 51, is_hpc = 1 /** uint64_t x; @@ -202,7 +202,7 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, } if (min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); - select_mz(p, len, MAX_HIGH_OCC); + if (sample_dist > w) select_mz(p, len, MAX_HIGH_OCC); for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts p->a[i].rid = rid; } From 3b8f73fbab121ab343efff887fd36a7918675304 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Fri, 22 Jan 2021 10:31:29 -0500 Subject: [PATCH 08/10] r332: increased mz sample dist from 200 to 500 --- CommandLines.cpp | 2 +- CommandLines.h | 4 ++-- anchor.cpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/CommandLines.cpp b/CommandLines.cpp index 3568aa7..85a37e6 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -101,7 +101,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->thread_num = 1; asm_opt->k_mer_length = 51; asm_opt->mz_win = 51; - asm_opt->sample_dist = 200; + asm_opt->mz_sample_dist = 500; asm_opt->bf_shift = 37; asm_opt->max_kmer_cnt = 2000; asm_opt->high_factor = 5.0; diff --git a/CommandLines.h b/CommandLines.h index 082e141..bd61e93 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.13-r330-dirty" +#define HA_VERSION "0.13-r332-dirty" #define VERBOSE 0 @@ -34,7 +34,7 @@ typedef struct { int thread_num; int k_mer_length; int mz_win; - int sample_dist; + int mz_sample_dist; int bf_shift; int max_kmer_cnt; double high_factor; // coverage cutoff set to high_factor*hom_cov diff --git a/anchor.cpp b/anchor.cpp index 9f0e83d..fc7fa56 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -73,7 +73,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg rlen = Get_READ_LENGTH(R_INF, rid); // read length // get the list of anchors - ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.sample_dist, k_flag, dbg_ct); + ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct); // minimizer of queried read if (ab->mz.m > ab->old_mz_m) { ab->old_mz_m = ab->mz.m; From 61befd357aad3bf318270e2675d874b6ad9199df Mon Sep 17 00:00:00 2001 From: Heng Li Date: Fri, 22 Jan 2021 23:17:25 -0500 Subject: [PATCH 09/10] r333: also subsample high-freq mz during indexing --- CommandLines.h | 2 +- htab.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/CommandLines.h b/CommandLines.h index bd61e93..56f23d5 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.13-r332-dirty" +#define HA_VERSION "0.13-r333-dirty" #define VERBOSE 0 diff --git a/htab.cpp b/htab.cpp index 936446c..37c476a 100644 --- a/htab.cpp +++ b/htab.cpp @@ -556,7 +556,7 @@ static void worker_for_mz(void *data, long i, int tid) ha_mz1_v *b = &s->mz_buf[tid]; s->mz_buf[tid].n = 0; ///s->p->opt->w = 51, s->p->opt->k - ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, 0, 0, 0); + ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab, asm_opt.mz_sample_dist, 0, 0); s->mz[i].n = s->mz[i].m = b->n; MALLOC(s->mz[i].a, b->n); memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t)); From 205eadf3467170d79563db324dbaf9349a615012 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sat, 23 Jan 2021 16:25:52 -0500 Subject: [PATCH 10/10] r334: ignore a high-occ seed in a small window --- CommandLines.h | 2 +- sketch.cpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/CommandLines.h b/CommandLines.h index 56f23d5..d7e753b 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.13-r333-dirty" +#define HA_VERSION "0.13-r334-dirty" #define VERBOSE 0 diff --git a/sketch.cpp b/sketch.cpp index f053073..72f9d39 100644 --- a/sketch.cpp +++ b/sketch.cpp @@ -67,7 +67,8 @@ static void select_mz(ha_mz1_v *p, int len, int sample_dist) } //ks_heapsort_mz(k, b); // sorting is not needed for now for (j = 0; j < k; ++j) - p->a[b[j].pos].rid = 0; + if (b[j].rid < pe - ps) + p->a[b[j].pos].rid = 0; } last0 = i; }