diff --git a/Assembly.cpp b/Assembly.cpp index 7798f24..1b09dd9 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -645,6 +645,44 @@ static void worker_ovec(void *data, long i, int tid) } +static void worker_ovec_cal0(void *data, long i, int tid) +{ + ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid]; + int fully_cov, abnormal; + // if(i != 33) return; + // fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i); + // if (memcmp("m64012_190920_173625/88015004/ccs", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { + // fprintf(stderr, "[M::%s-beg] rid->%ld\n", __func__, i); + // } else { + // return; + // } + + ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, + 0.02, asm_opt.max_n_chain, 1, NULL/**&(b->k_flag)**/, &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region), NULL, &(b->sp)); + + clear_Cigar_record(&b->cigar1); + clear_Round2_alignment(&b->round2); + + correct_overlap(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon, + &b->cigar1, &b->hap, &b->round2, &b->r_buf, &(b->tmp_region.w_list), 0, 0, &fully_cov, &abnormal); + + b->num_read_base += b->self_read.length; + b->num_correct_base += b->correct.corrected_base; + b->num_recorrect_base += b->round2.dumy.corrected_base; + + + + R_INF.paf[i].is_fully_corrected = 0; + R_INF.paf[i].is_abnormal = abnormal; + R_INF.trio_flag[i] = AMBIGU; + + push_final_overlaps(&(R_INF.paf[i]), R_INF.reverse_paf, &b->olist, 1); + push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2); + + if(het_cnt) het_cnt[i] = get_het_cnt(&b->hap); +} + + static void worker_ovec_related_reads(void *data, long i, int tid) { ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid]; @@ -982,6 +1020,45 @@ void ha_overlap_and_correct(int round) ///debug_print_pob_regions(); } +void ha_overlap_cal(int round) +{ + int i, hom_cov, het_cov; + ha_ovec_buf_t **b; + ha_flt_tab_hp = ha_idx_hp = NULL; + + // overlap and correct reads + CALLOC(b, asm_opt.thread_num); + for (i = 0; i < asm_opt.thread_num; ++i) + b[i] = ha_ovec_init(0, (round == asm_opt.number_of_round - 1),0); + if(ha_idx) hom_cov = asm_opt.hom_cov; + if(ha_idx == NULL) ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, 0, &R_INF, &hom_cov, &het_cov); // build the index + ///debug_adapter(&asm_opt, &R_INF); + if (round == 0 && ha_flt_tab == 0) // then asm_opt.hom_cov hasn't been updated + ha_opt_update_cov(&asm_opt, hom_cov); + het_cnt = NULL; + if(round == asm_opt.number_of_round-1 && asm_opt.is_dbg_het_cnt) CALLOC(het_cnt, R_INF.total_reads); + // fprintf(stderr, "[M::%s-start]\n", __func__); + kt_for(asm_opt.thread_num, worker_ovec_cal0, b, R_INF.total_reads);///debug_for_fix + // fprintf(stderr, "[M::%s-end]\n", __func__); + + ha_pt_destroy(ha_idx); + ha_idx = NULL; + + if(het_cnt) { + print_het_cnt_log(het_cnt); free(het_cnt); het_cnt = NULL; + } + + // collect statistics + for (i = 0; i < asm_opt.thread_num; ++i) { + 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 += ha_ovec_mem(b[i], NULL); + ha_ovec_destroy(b[i]); + } + free(b); +} + void update_overlaps(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf, UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact) @@ -1518,6 +1595,52 @@ void Output_PAF() fprintf(stderr, "PAF has been written.\n"); } + +void Output_PAF0(ma_hit_t_alloc* sources, const char *prefix) +{ + fprintf(stderr, "Writing PAF to disk ...... \n"); + char* paf_name = (char*)malloc(strlen(asm_opt.output_file_name)+strlen(prefix)+50); + sprintf(paf_name, "%s.%s.ovlp.paf", asm_opt.output_file_name, prefix); + FILE* output_file = fopen(paf_name, "w"); + uint64_t i, j; + + for (i = 0; i < R_INF.total_reads; i++) + { + for (j = 0; j < sources[i].length; j++) + { + fwrite(Get_NAME(R_INF, Get_qn(sources[i].buffer[j])), 1, + Get_NAME_LENGTH(R_INF, Get_qn(sources[i].buffer[j])), output_file); + fwrite("\t", 1, 1, output_file); + fprintf(output_file, "%lu\t", (unsigned long)Get_READ_LENGTH(R_INF, Get_qn(sources[i].buffer[j]))); + fprintf(output_file, "%d\t", Get_qs(sources[i].buffer[j])); + fprintf(output_file, "%d\t", Get_qe(sources[i].buffer[j])); + if(sources[i].buffer[j].rev) + { + fprintf(output_file, "-\t"); + } + else + { + fprintf(output_file, "+\t"); + } + fwrite(Get_NAME(R_INF, Get_tn(sources[i].buffer[j])), 1, + Get_NAME_LENGTH(R_INF, Get_tn(sources[i].buffer[j])), output_file); + fwrite("\t", 1, 1, output_file); + fprintf(output_file, "%lu\t", (unsigned long)Get_READ_LENGTH(R_INF, Get_tn(sources[i].buffer[j]))); + fprintf(output_file, "%d\t", Get_ts(sources[i].buffer[j])); + fprintf(output_file, "%d\t", Get_te(sources[i].buffer[j])); + fprintf(output_file, "%d\t", sources[i].buffer[j].ml); + fprintf(output_file, "%d\t", sources[i].buffer[j].bl); + fprintf(output_file, "255\n"); + + } + } + + free(paf_name); + fclose(output_file); + + fprintf(stderr, "PAF has been written.\n"); +} + int check_cluster(uint64_t* list, long long listLen, ma_hit_t_alloc* paf, float threshold) { long long i, k; @@ -1747,6 +1870,37 @@ void ug_idx_build(ma_ug_t *ug, int hap_n) exit(1); } +int ha_assemble_ovec(void) +{ + extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o); + int r = 0, hom_cov = -1; + + ha_flt_tab = ha_idx = NULL; + + // construct hash table for high occurrence k-mers + if (!(asm_opt.flag & HA_F_NO_KMER_FLT) && ha_flt_tab == NULL) + { + ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF, &hom_cov, 0); + ha_opt_update_cov(&asm_opt, hom_cov); + } + // error correction + assert(asm_opt.number_of_round > 0); + ha_opt_reset_to_round(&asm_opt, r); // this update asm_opt.roundID and a few other fields + ha_overlap_cal(r); + fprintf(stderr, "[M::%s] size of buffer: %.3fGB\n", __func__, asm_opt.mem_buf / 1073741824.0); + + + ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads); + ha_ft_destroy(ha_flt_tab); + + Output_PAF0(R_INF.paf, "0"); + Output_PAF0(R_INF.reverse_paf, "1"); + if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF(); + + destory_All_reads(&R_INF); + return 0; +} + int ha_assemble(void) { // debug_mc_g_t(MC_NAME); diff --git a/Assembly.h b/Assembly.h index 68ab3ff..37ea0ec 100644 --- a/Assembly.h +++ b/Assembly.h @@ -48,7 +48,7 @@ void ug_idx_build(ma_ug_t *ug, int hap_n); ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov, int is_ug); ha_ovec_buf_t *ha_ovec_buf_init(void *km, int is_final, int save_ov, int is_ug); void ha_ovec_destroy(ha_ovec_buf_t *b); -int64_t ha_ovec_mem(const ha_ovec_buf_t *b, int64_t *mem_a) -; +int64_t ha_ovec_mem(const ha_ovec_buf_t *b, int64_t *mem_a); +int ha_assemble_ovec(void); #endif diff --git a/CommandLines.cpp b/CommandLines.cpp index 761a18b..8672ced 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -57,6 +57,7 @@ static ko_longopt_t long_options[] = { { "ul-round", ko_required_argument, 342}, { "prt-raw", ko_no_argument, 343}, { "integer-correct", ko_required_argument, 344}, + { "dbg-ovec", ko_no_argument, 345}, { 0, 0, 0 } }; @@ -269,6 +270,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->ul_clean_round = 1; asm_opt->prt_dbg_gfa = 0; asm_opt->integer_correct_round = 0; + asm_opt->dbg_ovec_cal = 0; } void destory_enzyme(enzyme* f) @@ -804,6 +806,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 342) asm_opt->ul_clean_round = atol(opt.arg); else if (c == 343) asm_opt->prt_dbg_gfa = 1; else if (c == 344) asm_opt->integer_correct_round = atol(opt.arg); + else if (c == 345) asm_opt->dbg_ovec_cal = 1; else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg); } diff --git a/CommandLines.h b/CommandLines.h index e737492..84de370 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.19.0-r554" +#define HA_VERSION "0.19.0-r555" #define VERBOSE 0 @@ -141,6 +141,7 @@ typedef struct { int32_t ul_clean_round; int32_t prt_dbg_gfa; int32_t integer_correct_round; + uint8_t dbg_ovec_cal; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/gfa_ut.cpp b/gfa_ut.cpp index f20c097..06c884c 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -16488,7 +16488,7 @@ ul2ul_idx_t *gen_ul2ul(ul_resolve_t *uidx, ug_opt_t *uopt, ulg_opt_t *ulopt, uin u2g_clean(uidx, ulopt, keep_raw_utg, is_bridg); // renew_ul2_utg(uidx); - output_integer_graph(uidx, z->i_ug, asm_opt.output_file_name, 0); + // output_integer_graph(uidx, z->i_ug, asm_opt.output_file_name, 0); return z; } diff --git a/main.cpp b/main.cpp index d58a8e1..c4a0816 100644 --- a/main.cpp +++ b/main.cpp @@ -61,8 +61,9 @@ int main(int argc, char *argv[]) // ed_band_cal_global((char*)"ACTTTTTT", 8, (char*)"AATTTT", 6, 3), // ed_band_cal_global_128bit((char*)"ACTTTTTT", 8, (char*)"AATTTT", 6, 3)); // exit(1); + if(!(asm_opt.dbg_ovec_cal)) ret = ha_assemble(); + else ret = ha_assemble_ovec(); - ret = ha_assemble(); destory_opt(&asm_opt); fprintf(stderr, "[M::%s] Version: %s\n", __func__, HA_VERSION); fprintf(stderr, "[M::%s] CMD:", __func__);