From 84361ba6ddca2c9fe2adeb1703b7febf71f24601 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Tue, 24 Mar 2020 23:05:38 -0400 Subject: [PATCH] dense minimizer apparently works --- Assembly.h | 2 +- Trio.cpp | 2 +- main.cpp | 4 ++- sketch.cpp | 2 +- yak-count.cpp | 67 ++++++++++++++++++++++++++++++--------------------- yak.h | 2 +- 6 files changed, 46 insertions(+), 33 deletions(-) diff --git a/Assembly.h b/Assembly.h index 81e1104..44c7ef9 100644 --- a/Assembly.h +++ b/Assembly.h @@ -8,7 +8,7 @@ #define Get_Cigar_Type(RECORD) (RECORD&3) #define Get_Cigar_Length(RECORD) (RECORD>>2) -void ha_count_high(const hifiasm_opt_t *asm_opt); +void *ha_gen_flt_tab(const hifiasm_opt_t *asm_opt); void Counting_multiple_thr(); void Build_hash_table_multiple_thr(); void Overlap_calculate_multipe_thr(); diff --git a/Trio.cpp b/Trio.cpp index ff2f7eb..27b3dc8 100644 --- a/Trio.cpp +++ b/Trio.cpp @@ -98,7 +98,7 @@ static yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, .. yak_ht_t *h = ch->h[i].h; fread(t, 4, 2, fp); if (ch0 == 0) yak_ht_resize(h, t[0]); - for (j = 0; j < t[1]; ++j) { + for (j = 0; j < (int)t[1]; ++j) { uint64_t key; fread(&key, 8, 1, fp); if (mode == YAK_LOAD_ALL) { diff --git a/main.cpp b/main.cpp index c666d00..4ffd8a7 100644 --- a/main.cpp +++ b/main.cpp @@ -9,13 +9,15 @@ int main(int argc, char *argv[]) { int i; + void *flt_tab; init_opt(&asm_opt); if (!CommandLine_process(argc, argv, &asm_opt)) return 1; yak_reset_realtime(); - ha_count_high(&asm_opt); + flt_tab = ha_gen_flt_tab(&asm_opt); + ha_hf_destroy(flt_tab); if (0) { Correct_Reads(asm_opt.number_of_round); } diff --git a/sketch.cpp b/sketch.cpp index 05ddf80..d2270f8 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(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, void *hf) +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) { static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 }; uint64_t shift1 = k - 1, mask = (1ULL<bf_shift = 0; o->bf_n_hash = 4; o->k = 31; + o->w = 1; o->pre = YAK_COUNTER_BITS; o->n_thread = 4; o->chunk_size = 10000000; @@ -283,10 +284,10 @@ static void count_seq_buf_HPC(ch_buf_t *buf, int k, int p, int len, const char * typedef struct { // global data structure for kt_pipeline() const yak_copt_t *opt; - int create_new, is_mz, is_store, mz_win; + const void *flt_tab; + int create_new, is_store; kseq_t *ks; yak_ch_t *h; - void *hf; } pl_data_t; typedef struct { // data structure for each step in kt_pipeline() @@ -312,9 +313,9 @@ static void worker_for_mz(void *data, long i, int tid) st_data_t *s = (st_data_t*)data; ha_mz1_v *b = &s->mz_buf[tid]; s->mz_buf[tid].n = 0; - ha_sketch(s->seq[i], s->len[i], s->p->mz_win, s->p->opt->k, 0, s->p->opt->is_HPC, b, s->p->hf); + ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, 0, s->p->opt->is_HPC, b, s->p->flt_tab); s->mz[i].n = s->mz[i].m = b->n; - MALLOC(s->mz[i].a, s->mz[i].n); + MALLOC(s->mz[i].a, b->n); memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t)); } @@ -346,33 +347,42 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip else return s; } else if (step == 1) { // step 2: extract k-mers st_data_t *s = (st_data_t*)in; - int i, n = 1<opt->pre, m; - if (!p->is_mz) { // enumerate all k-mers - CALLOC(s->buf, n); - m = (int)(s->nk * 1.2 / n) + 1; - for (i = 0; i < n; ++i) { - s->buf[i].m = m; - MALLOC(s->buf[i].a, m); - } + int i, n_pre = 1<opt->pre, m; + // allocate the k-mer buffer + CALLOC(s->buf, n_pre); + m = (int)(s->nk * 1.2 / n_pre) + 1; + for (i = 0; i < n_pre; ++i) { + s->buf[i].m = m; + MALLOC(s->buf[i].a, m); + } + // fill the buffer + if (p->opt->w == 1) { // enumerate all k-mers for (i = 0; i < s->n_seq; ++i) { if (p->opt->is_HPC) count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]); else count_seq_buf(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]); - if (!p->is_store) - free(s->seq[i]); + if (!p->is_store) free(s->seq[i]); } } else { // minimizers only - CALLOC(s->mz_buf, p->opt->n_thread); + // compute minimizers CALLOC(s->mz, s->n_seq); + CALLOC(s->mz_buf, p->opt->n_thread); kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq); for (i = 0; i < p->opt->n_thread; ++i) free(s->mz_buf[i].a); free(s->mz_buf); - if (!p->is_store) { - for (i = 0; i < s->n_seq; ++i) - free(s->seq[i]); + // insert minimizers + for (i = 0; i < s->n_seq; ++i) { + uint32_t j; + for (j = 0; j < s->mz[i].n; ++j) + ch_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x); } + for (i = 0; i < s->n_seq; ++i) { + free(s->mz[i].a); + if (!p->is_store) free(s->seq[i]); + } + free(s->mz); } free(s->seq); free(s->len); s->seq = 0, s->len = 0; @@ -395,13 +405,14 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip return 0; } -static yak_ch_t *yak_count(const char *fn, const yak_copt_t *opt, yak_ch_t *h0) +static yak_ch_t *yak_count(const char *fn, const yak_copt_t *opt, yak_ch_t *h0, const void *flt_tab) { pl_data_t pl; gzFile fp; if ((fp = gzopen(fn, "r")) == 0) return 0; memset(&pl, 0, sizeof(pl_data_t)); pl.ks = kseq_init(fp); + pl.flt_tab = flt_tab; pl.opt = opt; if (h0) { pl.h = h0, pl.create_new = 0; @@ -416,12 +427,12 @@ static yak_ch_t *yak_count(const char *fn, const yak_copt_t *opt, yak_ch_t *h0) return pl.h; } -static yak_ch_t *yak_count_file(const yak_copt_t *opt, int n_fn, char **fn) +static yak_ch_t *yak_count_file(const yak_copt_t *opt, int n_fn, char **fn, const void *flt_tab) { int i; yak_ch_t *h = 0; for (i = 0; i < n_fn; ++i) - h = yak_count(fn[i], opt, h); + h = yak_count(fn[i], opt, h, flt_tab); if (opt->bf_shift > 0) yak_ch_destroy_bf(h); return h; @@ -447,9 +458,9 @@ static yak_hh_t *gen_hh(const yak_ch_t *h) return hh; } -void *ha_count_high(const hifiasm_opt_t *asm_opt) +void *ha_gen_flt_tab(const hifiasm_opt_t *asm_opt) { - yak_hh_t *high_ht; + yak_hh_t *flt_tab; int64_t cnt[YAK_N_COUNTS]; int peak_hom, peak_het, cutoff; yak_copt_t opt; @@ -459,18 +470,18 @@ void *ha_count_high(const hifiasm_opt_t *asm_opt) opt.k = asm_opt->k_mer_length; opt.n_thread = asm_opt->thread_num; opt.bf_shift = asm_opt->bf_shift; - h = yak_count_file(&opt, asm_opt->num_reads, asm_opt->read_file_names); + h = yak_count_file(&opt, asm_opt->num_reads, asm_opt->read_file_names, 0); yak_ch_hist(h, cnt, opt.n_thread); peak_hom = yak_analyze_count(YAK_N_COUNTS, cnt, &peak_het); if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het); cutoff = (int)(peak_hom * asm_opt->high_factor); if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1; yak_ch_shrink(h, cutoff, YAK_MAX_COUNT, opt.n_thread); - high_ht = gen_hh(h); + flt_tab = gen_hh(h); yak_ch_destroy(h); fprintf(stderr, "[M::%s] filtered out %ld k-mers occurring %d or more times\n", - __func__, (long)kh_size(high_ht), cutoff); - return (void*)high_ht; + __func__, (long)kh_size(flt_tab), cutoff); + return (void*)flt_tab; } int ha_hf_isflt(const void *hh, uint64_t y) diff --git a/yak.h b/yak.h index 1218dc0..e06987b 100644 --- a/yak.h +++ b/yak.h @@ -30,7 +30,7 @@ void yak_reset_realtime(void); double yak_realtime(void); long yak_peakrss(void); -void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, void *hf); +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 yak_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het); yak_bf_t *yak_bf_init(int n_shift, int n_hashes);