diff --git a/CommandLines.h b/CommandLines.h index f8c7308..e7371c1 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.14-r310" +#define HA_VERSION "0.14-r311" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index 615812d..6c24438 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -11553,10 +11553,9 @@ kvec_asg_arc_t_warp* new_rtg_edges, int max_hang, int min_ovlp) uint8_t* primary_flag = (uint8_t*)calloc(sg->n_seq, sizeof(uint8_t)); int tmp_cov = asm_opt.hom_global_coverage; asm_opt.hom_global_coverage = -1; - fprintf(stderr, "xxx-0-xxx\n"); + purge_dups(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, 0, 0, 0, 1, NULL); - fprintf(stderr, "xxx-1-xxx\n"); dip_thre_max = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE)*0.70; asm_opt.hom_global_coverage = tmp_cov; ///fprintf(stderr, "dip_thre_max: %lu\n", dip_thre_max); @@ -25201,10 +25200,8 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, uint32_t chai memset(&bub, 0, sizeof(bubble_type)); copy_ug = copy_untig_graph(ug); - fprintf(stderr, "\nstart: sbsbsb-0-sbsbsb\n"); reset_bub(&bub, ug, sg, copy_ug, &link, ©_link, ruIndex, coverage_cut, sources, reverse_sources, max_hang, min_ovlp, &new_rtg_edges); - fprintf(stderr, "start: sbsbsb-0-sbsbsb\n"); beg_idx = bub.f_bub; occ = bub.b_bub + bub.b_end_bub + bub.tangle_bub; rescue_bubbles_by_contained_reads(ug, sg, sources, coverage_cut, ruIndex, max_hang, min_ovlp, chainLenThres, beg_idx, occ, &bub); ///output_unitig_graph(sg, coverage_cut, (char*)"debug_1.hic", sources, ruIndex, max_hang, min_ovlp); @@ -25212,20 +25209,16 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, uint32_t chai ma_ug_destroy(ug); ug = NULL; ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); - fprintf(stderr, "\nstart: sbsbsb-1-sbsbsb\n"); reset_bub(&bub, ug, sg, copy_ug, &link, ©_link, ruIndex, coverage_cut, sources, reverse_sources, max_hang, min_ovlp, &new_rtg_edges); - fprintf(stderr, "start: sbsbsb-1-sbsbsb\n"); beg_idx = bub.f_bub; occ = bub.b_bub + bub.b_end_bub + bub.tangle_bub; rescue_bubbles_by_missing_ovlp(ug, sg, sources, coverage_cut, ruIndex, max_hang, min_ovlp, chainLenThres, beg_idx, occ, &bub); ///output_unitig_graph(sg, coverage_cut, (char*)"debug_2.hic", sources, ruIndex, max_hang, min_ovlp); ma_ug_destroy(ug); ug = NULL; ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); - fprintf(stderr, "\nstart: sbsbsb-2-sbsbsb\n"); reset_bub(&bub, ug, sg, copy_ug, &link, ©_link, ruIndex, coverage_cut, sources, reverse_sources, max_hang, min_ovlp, &new_rtg_edges); - fprintf(stderr, "start: sbsbsb-2-sbsbsb\n"); beg_idx = bub.f_bub; occ = bub.b_bub + bub.b_end_bub + bub.tangle_bub; rescue_bubbles_by_missing_ovlp_backward(ug, sg, sources, coverage_cut, ruIndex, max_hang, min_ovlp, chainLenThres, beg_idx, occ, &bub); ///output_unitig_graph(sg, coverage_cut, (char*)"debug_3.hic", sources, ruIndex, max_hang, min_ovlp); diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index 8044a7b..b0e36a6 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -1500,18 +1500,6 @@ void quick_LIS(asg_arc_t_offset* x, uint32_t n, kvec_t_i32_warp* tailIndex, kvec tailIndex->a.a[0] = 0; prevIndex->a.a[0] = -1; - /*******************************for debug************************************/ - if(n == 3) - { - fprintf(stderr, "n=%u\n", n); - for(i = 0; i < (long long)n; i++) - { - fprintf(stderr, "i: %lld, Get_xOff: %lld, Get_yOff: %lld\n", - i, Get_xOff(x[i].Off), Get_yOff(x[i].Off)); - } - } - /*******************************for debug************************************/ - ///x has already sorted by x_pos for(i = 1; i < (long long)n; i++) { @@ -1532,13 +1520,6 @@ void quick_LIS(asg_arc_t_offset* x, uint32_t n, kvec_t_i32_warp* tailIndex, kvec // future subsequence // It will replace ceil value in tailIndices pos = GetCeilIndex(x, tailIndex, -1, len - 1, Get_yOff(x[i].Off)); - /*******************************for debug************************************/ - if(n == 3) - { - fprintf(stderr, "++++i: %lld, Get_xOff: %lld, Get_yOff: %lld, pos: %lld\n", - i, Get_xOff(x[i].Off), Get_yOff(x[i].Off), pos); - } - /*******************************for debug************************************/ prevIndex->a.a[i] = pos > 0? tailIndex->a.a[pos - 1] : -1; tailIndex->a.a[pos] = i; } @@ -1644,9 +1625,7 @@ uint32_t* xBeg, uint32_t* xEnd, uint32_t* yBeg, uint32_t* yEnd) u_buffer->a.n = m; ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "after"); - fprintf(stderr, "+sbsbsbsbsbsb-quick_LIS, u_buffer->a.n: %u\n", (uint32_t)u_buffer->a.n); quick_LIS(u_buffer->a.a, u_buffer->a.n, tailIndex, prevIndex); - fprintf(stderr, "-sbsbsbsbsbsb-quick_LIS, u_buffer->a.n: %u\n", (uint32_t)u_buffer->a.n); if(tailIndex->a.n == 0) return; @@ -1664,13 +1643,10 @@ kvec_t_i32_warp* prevIndex, long long* r_x_pos_beg, long long* r_x_pos_end, long { uint32_t x_pos_beg, y_pos_beg, x_pos_end, y_pos_end; /*************************x***************************/ - fprintf(stderr, "+sbsbsbsbsbsb-xUid: %u, yUid: %u, Get_x_beg(*hap_can): %u, Get_x_end(*hap_can): %u, Get_y_beg(*hap_can): %u, Get_y_end(*hap_can): %u, Get_rev(*hap_can): %u\n", - xUid, yUid, Get_x_beg(*hap_can), Get_x_end(*hap_can), Get_y_beg(*hap_can), Get_y_end(*hap_can), Get_rev(*hap_can)); get_base_boundary_advance(ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, min_ovlp, xReads, yReads, xUid, yUid, Get_x_beg(*hap_can), Get_x_end(*hap_can), Get_y_beg(*hap_can), Get_y_end(*hap_can), Get_rev(*hap_can), u_buffer, tailIndex, prevIndex, &x_pos_beg, &x_pos_end, &y_pos_beg, &y_pos_end); - fprintf(stderr, "-sbsbsbsbsbsb-xUid: %u, yUid: %u\n", xUid, yUid); /*************************x***************************/ if(x_pos_beg == (uint32_t)-1 || y_pos_beg == (uint32_t)-1 @@ -1948,12 +1924,10 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon hap_can->y_beg_pos = r_y_interval_beg; hap_can->y_end_pos = r_y_interval_end; - fprintf(stderr, "+inner-xUid: %u, yUid: %u\n", xUid, yUid); hap_can->index_end = determine_hap_overlap_type_advance(hap_can, xReads, yReads, ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, min_ovlp, xUid, yUid, u_buffer, tailIndex, prevIndex, r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end); - fprintf(stderr, "-inner-xUid: %u, yUid: %u\n", xUid, yUid); if(hap_can->index_end == XCY && yReads->len > (xReads->len*2)) return NON_PLOID; if(hap_can->index_end == YCX && xReads->len > (yReads->len*2)) return NON_PLOID; @@ -2825,7 +2799,6 @@ void hap_alignment_worker(void *_data, long eid, int tid) static void hap_alignment_advance_worker(void *_data, long eid, int tid) { - fprintf(stderr, "\neid: %ld, tid: %d\n", eid, tid); hap_alignment_struct_pip* hap_buf = (hap_alignment_struct_pip*)_data; ma_ug_t *ug = hap_buf->ug; asg_t *read_g = hap_buf->read_g; @@ -2997,20 +2970,14 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) is_update = 0; if(u_can->a.a[k].weight < Get_match(hap_can)*Hap_rate) continue; - fprintf(stderr, "+k: %u, u_can->a.n: %u, xUid: %u, yUid: %u\n", - k, (uint32_t)u_can->a.n, xUid, yUid); if(calculate_pair_hap_similarity_advance(&(u_can->a.a[k]), position_index, xUid, yUid, xReads, yReads, sources, reverse_sources, read_g, ruIndex, coverage_cut, Hap_rate, max_hang, min_ovlp, cov_threshold, u_buffer, score_vc, prevIndex_vec, &r_x_pos_beg, &r_x_pos_end, &r_y_pos_beg, &r_y_pos_end)!=PLOID) { - fprintf(stderr, "-0-k: %u, u_can->a.n: %u, xUid: %u, yUid: %u\n", - k, (uint32_t)u_can->a.n, xUid, yUid); continue; } - fprintf(stderr, "-1-k: %u, u_can->a.n: %u, xUid: %u, yUid: %u\n", - k, (uint32_t)u_can->a.n, xUid, yUid); if(Get_match(hap_can) < Get_match(u_can->a.a[k])) { @@ -4212,13 +4179,8 @@ uint32_t minLen, double purge_threshold) { asg_t* nsg = ug->g; uint64_t v, k, total_bases = 0, alter_bases = 0, primary_bases = 0, purge_bases = 0; - /*******************************for debug************************************/ - fprintf(stderr, "bbb-0-bbb\n"); - fprintf(stderr, "asm_opt.thread_num: %d\n", asm_opt.thread_num); - ///kt_for(asm_opt.thread_num, hap_alignment_advance_worker, hap_buf, nsg->n_seq); - kt_for(1, hap_alignment_advance_worker, hap_buf, nsg->n_seq); - fprintf(stderr, "bbb-1-bbb\n"); - /*******************************for debug************************************/ + kt_for(asm_opt.thread_num, hap_alignment_advance_worker, hap_buf, nsg->n_seq); + filter_hap_overlaps_by_length(all_ovlp, minLen); normalize_hap_overlaps_advance(all_ovlp, back_all_ovlp, ug, read_g, reverse_sources, ruIndex); @@ -4333,7 +4295,7 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link) init_hap_alignment_struct_pip(&hap_buf, asm_opt.thread_num, nsg->n_seq, ug, read_g, sources, reverse_sources, ruIndex, coverage_cut, position_index, density, max_hang, min_ovlp, 0.05, &all_ovlp); - fprintf(stderr, "aaa-0-aaa\n"); + if(hap_buf.cov_threshold < 0) { if(if_ploid_sample(ug, read_g, ruIndex, sources, reverse_sources, coverage_cut, @@ -4351,10 +4313,9 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link) if(asm_opt.hom_global_coverage == -1) asm_opt.hom_global_coverage = hap_buf.cov_threshold; fprintf(stderr, "[M::%s] purge duplication coverage threshold: %lld\n", __func__, hap_buf.cov_threshold); if(just_coverage) goto end_coverage; - fprintf(stderr, "aaa-1-aaa\n"); + ///kt_for(asm_opt.thread_num, hap_alignment_worker, &hap_buf, nsg->n_seq); kt_for(asm_opt.thread_num, hap_alignment_advance_worker, &hap_buf, nsg->n_seq); - fprintf(stderr, "aaa-2-aaa\n"); ///if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp);