From ef8c04755920df0dd7b2d06128b1926ea048861c Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 29 Oct 2020 11:28:43 -0400 Subject: [PATCH] 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