diff --git a/Assembly.cpp b/Assembly.cpp index 3e10111..89e2a4f 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -382,7 +382,7 @@ typedef struct { overlap_region_alloc olist; ha_abuf_t *ab; // error correction related buffers - int64_t num_read_base, num_correct_base, num_recorrect_base, mem_buf; + int64_t num_read_base, num_correct_base, num_recorrect_base; Cigar_record cigar1; Graph POA_Graph; Graph DAGCon; @@ -626,7 +626,7 @@ void Output_corrected_reads() void ha_overlap_and_correct(int round) { - int i; + int i, hom_cov; ha_ovec_buf_t **b; ha_ecsave_buf_t *e; @@ -634,7 +634,9 @@ void ha_overlap_and_correct(int round) CALLOC(b, asm_opt.thread_num); for (i = 0; i < asm_opt.thread_num; ++i) b[i] = ha_ovec_init(0); - ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF); // build the index + ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF, &hom_cov); // build the index + if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated + ha_opt_update_cov(&asm_opt, hom_cov); if (asm_opt.required_read_name) kt_for(asm_opt.thread_num, worker_ovec_related_reads, b, R_INF.total_reads); else @@ -647,7 +649,7 @@ void ha_overlap_and_correct(int round) asm_opt.num_bases += b[i]->num_read_base; asm_opt.num_corrected_bases += b[i]->num_correct_base; asm_opt.num_recorrected_bases += b[i]->num_recorrect_base; - asm_opt.mem_buf += b[i]->mem_buf; + asm_opt.mem_buf += ha_ovec_mem(b[i]); ha_ovec_destroy(b[i]); } free(b); @@ -1135,12 +1137,12 @@ long long readNum, long long rescue_threshold, float cluster_threshold) void ha_overlap_final(void) { - int i; + int i, hom_cov; ha_ovec_buf_t **b; CALLOC(b, asm_opt.thread_num); for (i = 0; i < asm_opt.thread_num; ++i) b[i] = ha_ovec_init(1); - ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF); // build the index + ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF, &hom_cov); // build the index kt_for(asm_opt.thread_num, worker_ov_final, b, R_INF.total_reads); ha_pt_destroy(ha_idx); ha_idx = 0; @@ -1151,19 +1153,21 @@ void ha_overlap_final(void) int ha_assemble(void) { - int r, ovlp_loaded = 0; + int r, hom_cov = -1, ovlp_loaded = 0; if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) { ovlp_loaded = 1; fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage()); } if (!ovlp_loaded) { // construct hash table for high occurrence k-mers - if (!asm_opt.no_kmer_flt) - 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, &hom_cov); + ha_opt_update_cov(&asm_opt, hom_cov); + } // error correction assert(asm_opt.number_of_round > 0); for (r = 0; r < asm_opt.number_of_round; ++r) { - clear_opt(&asm_opt, r); // this update asm_opt.roundID and a few other fields + ha_opt_reset_to_round(&asm_opt, r); // this update asm_opt.roundID and a few other fields ha_overlap_and_correct(r); fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> corrected reads for round %d\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb(), r + 1); @@ -1174,7 +1178,7 @@ int ha_assemble(void) //Output_corrected_reads(); fprintf(stderr, "[M::%s::%.3f*%.2f] ==> written corrected reads to disk\n", __func__, yak_realtime(), yak_cpu_usage()); // overlap between corrected reads - clear_opt(&asm_opt, asm_opt.number_of_round); + ha_opt_reset_to_round(&asm_opt, asm_opt.number_of_round); ha_overlap_final(); fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb()); diff --git a/CommandLines.cpp b/CommandLines.cpp index 4ba6ba2..ec55c95 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -35,15 +35,11 @@ void Print_H(hifiasm_opt_t* asm_opt) 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 >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor); - fprintf(stderr, " -N INT consider up to INT overlaps for each oriented read [%d]\n", asm_opt->max_n_chain); - ///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, " -N INT consider up to max(-D*coverage,-N) overlaps for each oriented read [%d]\n", asm_opt->max_n_chain); fprintf(stderr, " -i ignore saved overlaps in *.ovlp* files\n"); fprintf(stderr, " -z INT length of adapters that should be removed [%d]\n", asm_opt->adapterLen); - fprintf(stderr, " -m INT size of popped large bubbles for contig graph [%lld]\n", - asm_opt->large_pop_bubble_size); - fprintf(stderr, " -p INT size of popped small bubbles for haplotype-resolved unitig graph [%lld]\n", - asm_opt->small_pop_bubble_size); + fprintf(stderr, " -m INT size of popped large bubbles for contig graph [%lld]\n", asm_opt->large_pop_bubble_size); + fprintf(stderr, " -p INT size of popped small bubbles for haplotype-resolved unitig graph [%lld]\n", asm_opt->small_pop_bubble_size); fprintf(stderr, " -n INT small removed unitig threshold [%d]\n", asm_opt->max_short_tip); fprintf(stderr, " -x FLOAT max overlap drop ratio [%.2g]\n", asm_opt->max_drop_rate); fprintf(stderr, " -y FLOAT min overlap drop ratio [%.2g]\n", asm_opt->min_drop_rate); @@ -51,8 +47,8 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " -h show help information\n"); fprintf(stderr, " Trio-partition:\n"); - fprintf(stderr, " -P FILE paternal trio index generated by \"yak count\" [NULL]\n"); - fprintf(stderr, " -M FILE maternal trio index generated by \"yak count\" [NULL]\n"); + fprintf(stderr, " -P FILE paternal trio index generated by \"yak count\" []\n"); + fprintf(stderr, " -M FILE maternal trio index generated by \"yak count\" []\n"); fprintf(stderr, " -c INT lower bound of the binned k-mer's frequency [%d]\n", asm_opt->min_cnt); fprintf(stderr, " -d INT upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt); @@ -73,10 +69,10 @@ 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->high_factor = 5.0f; asm_opt->no_HPC = 0; asm_opt->no_kmer_flt = 0; - asm_opt->max_n_chain = 400; + asm_opt->high_factor = 5.0f; + asm_opt->max_n_chain = 100; asm_opt->k_mer_min_freq = 3; asm_opt->k_mer_max_freq = 66; asm_opt->load_index_from_disk = 1; @@ -107,7 +103,7 @@ void destory_opt(hifiasm_opt_t* asm_opt) } } -void clear_opt(hifiasm_opt_t* asm_opt, int round) +void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round) { asm_opt->num_bases = 0; asm_opt->num_corrected_bases = 0; @@ -116,6 +112,14 @@ void clear_opt(hifiasm_opt_t* asm_opt, int round) asm_opt->roundID = round; } +void ha_opt_update_cov(hifiasm_opt_t *opt, int hom_cov) +{ + int max_n_chain = (int)(hom_cov * opt->high_factor + .499); + if (opt->max_n_chain < max_n_chain) + opt->max_n_chain = max_n_chain; + fprintf(stderr, "[M::%s] updated max_n_chain to %d\n", __func__, opt->max_n_chain); +} + int check_file(char* name, const char* opt) { if(!name) @@ -308,7 +312,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:FN:", long_options)) >= 0) { + while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:w:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:", long_options)) >= 0) { if (c == 'h') { Print_H(asm_opt); @@ -325,7 +329,6 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 'r') asm_opt->number_of_round = atoi(opt.arg); else if (c == 'k') asm_opt->k_mer_length = atoi(opt.arg); else if (c == 'i') asm_opt->load_index_from_disk = 0; - 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; diff --git a/CommandLines.h b/CommandLines.h index 3039737..7df0950 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.3.0-dirty-r200" +#define HA_VERSION "0.3.0-dirty-r201" #define VERBOSE 0 @@ -18,10 +18,10 @@ typedef struct { int k_mer_length; int mz_win; int bf_shift; - float high_factor; int no_HPC; int no_kmer_flt; - int max_n_chain; + float high_factor; // coverage cutoff set to high_factor*hom_cov + int max_n_chain; // fall-back max number of chains to consider int k_mer_min_freq; int k_mer_max_freq; int load_index_from_disk; @@ -56,7 +56,8 @@ extern hifiasm_opt_t asm_opt; void init_opt(hifiasm_opt_t* asm_opt); void destory_opt(hifiasm_opt_t* asm_opt); -void clear_opt(hifiasm_opt_t* asm_opt, int round); +void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round); +void ha_opt_update_cov(hifiasm_opt_t *opt, int hom_cov); int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt); double Get_T(void); diff --git a/htab.cpp b/htab.cpp index d96b21d..190ebed 100644 --- a/htab.cpp +++ b/htab.cpp @@ -771,7 +771,7 @@ void ha_ft_destroy(void *h) * High-level interfaces * *************************/ -void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs) +void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov) { yak_ft_t *flt_tab; int64_t cnt[YAK_N_COUNTS]; @@ -780,6 +780,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs) h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN, NULL, NULL, rs); ha_ct_hist(h, cnt, asm_opt->thread_num); peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het); + if (hom_cov) *hom_cov = peak_hom; 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; @@ -791,7 +792,7 @@ void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs) return (void*)flt_tab; } -ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, 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, int *hom_cov) { int64_t cnt[YAK_N_COUNTS], tot_cnt; int peak_hom, peak_het, i, extra_flag1, extra_flag2; @@ -812,6 +813,7 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f 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 (hom_cov) *hom_cov = peak_hom; if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het); if (flt_tab == 0) { int cutoff = (int)(peak_hom * asm_opt->high_factor); diff --git a/htab.h b/htab.h index 6298923..0355745 100644 --- a/htab.h +++ b/htab.h @@ -26,11 +26,11 @@ extern const unsigned char seq_nt4_table[256]; extern void *ha_flt_tab; extern ha_pt_t *ha_idx; -void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs); +void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov); int ha_ft_isflt(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, 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, int *hom_cov); void ha_pt_destroy(ha_pt_t *h); const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);