From da672118e54a37567be973e1b795887e9ce29c30 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 21 Jan 2021 11:05:37 -0500 Subject: [PATCH] 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; }