diff --git a/CommandLines.cpp b/CommandLines.cpp index 7400188..ef8f9c9 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -29,7 +29,7 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " -k INT k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length); fprintf(stderr, " -w INT minimizer window size [%d]\n", asm_opt->mz_win); fprintf(stderr, " -f INT number of bits for bloom filter [%d]\n", asm_opt->bf_shift); - fprintf(stderr, " -D FLOAT drop k-mers occuring more than FLOAT*coverage times [%.1f]\n", asm_opt->high_factor); + fprintf(stderr, " -D FLOAT drop k-mers occuring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor); ///fprintf(stderr, " -l load all overlaps from disk, can avoid overlap calculation [%d]\n", asm_opt->load_index_from_disk); ///fprintf(stderr, " -i ignore saved overlaps in *.ovlp*.bin files\n"); fprintf(stderr, " -i ignore saved overlaps in *.ovlp* files\n"); @@ -69,6 +69,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->bf_shift = 37; asm_opt->high_factor = 7.0f; asm_opt->no_HPC = 0; + asm_opt->no_kmer_flt = 0; asm_opt->k_mer_min_freq = 3; asm_opt->k_mer_max_freq = 66; asm_opt->load_index_from_disk = 1; @@ -300,7 +301,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) int c; - while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:lw:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:", 0)) >= 0) { + while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:lw:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:F", 0)) >= 0) { if (c == 'h') { Print_H(asm_opt); @@ -320,6 +321,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 'l') asm_opt->load_index_from_disk = 1; else if (c == 'w') asm_opt->mz_win = atoi(opt.arg); else if (c == 'D') asm_opt->high_factor = atof(opt.arg); + else if (c == 'F') asm_opt->no_kmer_flt = 1; else if (c == 'a') asm_opt->clean_round = atoi(opt.arg); else if (c == 'z') asm_opt->adapterLen = atoi(opt.arg); else if (c == 'b') asm_opt->required_read_name = opt.arg; diff --git a/CommandLines.h b/CommandLines.h index 5812255..65978e0 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -19,6 +19,7 @@ typedef struct { int bf_shift; float high_factor; int no_HPC; + int no_kmer_flt; int k_mer_min_freq; int k_mer_max_freq; int load_index_from_disk; diff --git a/htab.cpp b/htab.cpp index 76cab23..1995a16 100644 --- a/htab.cpp +++ b/htab.cpp @@ -764,7 +764,7 @@ int ha_ft_isflt(const void *hh, uint64_t y) void ha_ft_destroy(void *h) { - yak_ft_destroy((yak_ft_t*)h); + if (h) yak_ft_destroy((yak_ft_t*)h); } /************************* @@ -794,20 +794,36 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs) ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs) { int64_t cnt[YAK_N_COUNTS], tot_cnt; - int peak_hom, peak_het, i, extra_flag = read_from_store? HAF_RS_READ : HAF_RS_WRITE_SEQ; + int peak_hom, peak_het, i, extra_flag1, extra_flag2; ha_ct_t *ct; ha_pt_t *pt; - ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag, NULL, flt_tab, rs); + if (read_from_store) { + extra_flag1 = extra_flag2 = HAF_RS_READ; + } else if (rs->total_reads == 0) { + extra_flag1 = HAF_RS_WRITE_LEN; + extra_flag2 = HAF_RS_WRITE_SEQ; + } else { + extra_flag1 = HAF_RS_WRITE_SEQ; + extra_flag2 = HAF_RS_READ; + } + ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs); fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__, yak_realtime(), yak_cputime() / yak_realtime(), (long)ct->tot); ha_ct_hist(ct, cnt, asm_opt->thread_num); fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]); peak_hom = ha_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); - ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num); - for (i = 2, tot_cnt = 0; i <= YAK_MAX_COUNT - 1; ++i) tot_cnt += cnt[i] * i; + if (flt_tab == 0) { + int cutoff = (int)(peak_hom * asm_opt->high_factor); + if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1; + ha_ct_shrink(ct, 2, cutoff, asm_opt->thread_num); + for (i = 2, tot_cnt = 0; i <= cutoff; ++i) tot_cnt += cnt[i] * i; + } else { + ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num); + for (i = 2, tot_cnt = 0; i <= YAK_MAX_COUNT - 1; ++i) tot_cnt += cnt[i] * i; + } pt = ha_pt_gen(ct, asm_opt->thread_num); - ha_count(asm_opt, HAF_COUNT_EXACT|HAF_RS_READ, pt, flt_tab, rs); + ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, pt, flt_tab, rs); assert((uint64_t)tot_cnt == pt->tot_pos); //ha_pt_sort(pt, asm_opt->thread_num); fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__, diff --git a/main.cpp b/main.cpp index dab9c41..a8219b0 100644 --- a/main.cpp +++ b/main.cpp @@ -12,7 +12,8 @@ int main(int argc, char *argv[]) init_opt(&asm_opt); if (!CommandLine_process(argc, argv, &asm_opt)) return 1; yak_reset_realtime(); - ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF); + if (!asm_opt.no_kmer_flt) + ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF); Correct_Reads(asm_opt.number_of_round); ha_ft_destroy(ha_flt_tab); destory_All_reads(&R_INF);