From 98a04c168ac550200102bc0aa269c1926184934a Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Thu, 18 Mar 2021 05:26:33 -0400 Subject: [PATCH] update r313 --- CommandLines.cpp | 30 +++++--- CommandLines.h | 8 ++- Overlaps.cpp | 174 +++++++---------------------------------------- Overlaps.h | 5 +- Purge_Dups.cpp | 2 + hic.cpp | 10 +-- hifiasm.1 | 10 +-- 7 files changed, 62 insertions(+), 177 deletions(-) diff --git a/CommandLines.cpp b/CommandLines.cpp index f69f93e..5586fa8 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -23,7 +23,7 @@ static ko_longopt_t long_options[] = { { "ex-iter", ko_required_argument, 308 }, { "purge-cov", ko_required_argument, 309 }, { "pri-range", ko_required_argument, 310 }, - { "high-het", ko_no_argument, 311 }, + ///{ "high-het", ko_no_argument, 311 }, { "lowQ", ko_required_argument, 312 }, { "min-hist-cnt", ko_required_argument, 313 }, { "h1", ko_required_argument, 314 }, @@ -90,13 +90,13 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " Purge-dups:\n"); fprintf(stderr, " -l INT purge level. 0: no purging; 1: light; 2/3: aggressive [0 for trio; 2 for unzip]\n"); - fprintf(stderr, " -s FLOAT similarity threshold for duplicate haplotigs [%g]\n", - asm_opt->purge_simi_rate); + fprintf(stderr, " -s FLOAT similarity threshold for duplicate haplotigs [%g for -l1/-l2, %g for -l3]\n", + asm_opt->purge_simi_rate_l2, asm_opt->purge_simi_rate_l3); fprintf(stderr, " -O INT min number of overlapped reads for duplicate haplotigs [%d]\n", asm_opt->purge_overlap_len); fprintf(stderr, " --purge-cov INT\n"); fprintf(stderr, " coverage upper bound of Purge-dups [auto]\n"); - fprintf(stderr, " --high-het enable this mode for high heterozygosity sample [experimental, not stable]\n"); + ///fprintf(stderr, " --high-het enable this mode for high heterozygosity sample [experimental, not stable]\n"); fprintf(stderr, " Hi-C-partition [experimental, not stable]:\n"); fprintf(stderr, " --h1 FILEs file names of Hi-C R1 [r1_1.fq,r1_2.fq,...]\n"); @@ -149,10 +149,11 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->mid_cnt = 5; asm_opt->purge_level_primary = 2; asm_opt->purge_level_trio = 0; - asm_opt->purge_simi_rate = 0.75; - asm_opt->purge_simi_rate_hic = 0.85; + asm_opt->purge_simi_rate_l2 = 0.75; + asm_opt->purge_simi_rate_l3 = 0.55; + ///asm_opt->purge_simi_rate_hic = 0.85; asm_opt->purge_overlap_len = 1; - asm_opt->purge_overlap_len_hic = 50; + ///asm_opt->purge_overlap_len_hic = 50; asm_opt->recover_atg_cov_min = -1024; asm_opt->recover_atg_cov_max = INT_MAX; asm_opt->hom_global_coverage = -1; @@ -453,6 +454,12 @@ int check_option(hifiasm_opt_t* asm_opt) return 0; } + if(asm_opt->purge_simi_thres < 0) + { + fprintf(stderr, "[ERROR] [-s] must >= 0\n"); + return 0; + } + // fprintf(stderr, "input file num: %d\n", asm_opt->num_reads); // fprintf(stderr, "output file: %s\n", asm_opt->output_file_name); // fprintf(stderr, "number of threads: %d\n", asm_opt->thread_num); @@ -470,7 +477,7 @@ int check_option(hifiasm_opt_t* asm_opt) // fprintf(stderr, "mid_cnt: %d\n", asm_opt->mid_cnt); // fprintf(stderr, "purge_level_primary: %d\n", asm_opt->purge_level_primary); // fprintf(stderr, "purge_level_trio: %d\n", asm_opt->purge_level_trio); - // fprintf(stderr, "purge_simi_rate: %f\n", asm_opt->purge_simi_rate); + // fprintf(stderr, "purge_simi_thres: %f\n", asm_opt->purge_simi_thres); // fprintf(stderr, "purge_overlap_len: %d\n", asm_opt->purge_overlap_len); return 1; @@ -620,7 +627,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) asm_opt->recover_atg_cov_min = asm_opt->recover_atg_cov_max = -1; } } - else if (c == 311) asm_opt->flag |= HA_F_HIGH_HET; + ///else if (c == 311) asm_opt->flag |= HA_F_HIGH_HET; else if (c == 312) asm_opt->bed_inconsist_rate = atoi(opt.arg); else if (c == 313) asm_opt->min_hist_kmer_cnt = atoi(opt.arg); else if (c == 314) get_hic_enzymes(opt.arg, &(asm_opt->hic_reads[0]), 0); @@ -633,7 +640,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) { ///0: disable purge_dup; 1: purge containment; 2: purge overlap asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg); } - else if (c == 's') asm_opt->purge_simi_rate = atof(opt.arg); + else if (c == 's') asm_opt->purge_simi_rate_l2 = asm_opt->purge_simi_rate_l3 = atof(opt.arg); else if (c == 'O') asm_opt->purge_overlap_len = atoll(opt.arg); else if (c == ':') { @@ -647,6 +654,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) } } + if(asm_opt->purge_level_primary > 2) asm_opt->purge_simi_thres = asm_opt->purge_simi_rate_l3; + else asm_opt->purge_simi_thres = asm_opt->purge_simi_rate_l2; + if (argc == opt.ind) { diff --git a/CommandLines.h b/CommandLines.h index 253e6a8..4f263f5 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -70,7 +70,7 @@ typedef struct { int purge_level_primary; int purge_level_trio; int purge_overlap_len; - int purge_overlap_len_hic; + ///int purge_overlap_len_hic; int recover_atg_cov_min; int recover_atg_cov_max; int hom_global_coverage; @@ -80,8 +80,10 @@ typedef struct { float max_hang_rate; float min_drop_rate; float max_drop_rate; - float purge_simi_rate; - float purge_simi_rate_hic; + float purge_simi_rate_l2; + float purge_simi_rate_l3; + float purge_simi_thres; + ///float purge_simi_rate_hic; long long small_pop_bubble_size; long long large_pop_bubble_size; diff --git a/Overlaps.cpp b/Overlaps.cpp index 5be775b..2d9782c 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -11093,27 +11093,6 @@ uint32_t get_num_trio_flag(ma_ug_t *ug, uint32_t v, uint32_t flag) } -void set_pre_uid(buf_t* b, hc_links* link, ma_ug_t *ug) -{ - uint32_t k = 0, i = 0, m = 0, rId, pre = (uint32_t)-1; - ma_utg_t* u = NULL; - for (i = 0; i < b->b.n; i++) - { - u = &(ug->u.a[b->b.a[i]>>1]); - if(u->m == 0) continue; - for (k = 0; k < u->n; k++) - { - rId = u->a[k]>>33; - if(link->u_idx[rId] == (uint32_t)-1) continue; - if(pre == link->u_idx[rId]) continue; - pre = link->u_idx[rId]; - b->b.a[m] = link->u_idx[rId]; - m++; - } - } - b->b.n = m; -} - inline uint64_t get_utg_len(buf_t* b, ma_ug_t *ug, asg_t *read_sg, uint64_t ignore_end, uint64_t* len_thre, uint64_t* occ) { @@ -11217,77 +11196,6 @@ inline uint64_t get_utg_len(buf_t* b, ma_ug_t *ug, asg_t *read_sg, uint64_t igno return len; } -void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug, asg_t *read_sg) -{ - uint32_t b_0_i, b_0_k, b_1_i, b_1_k, pre_0, pre_1, rId_0, rId_1, ori_0, ori_1; - uint64_t d = RC_1, len_0, len_1, thre_0, thre_1; - ma_utg_t* u_b_0 = NULL; - ma_utg_t* u_b_1 = NULL; - if(b_0->b.n == 0 || b_1->b.n == 0) return; - - len_0 = get_utg_len(b_0, ug, read_sg, 1, NULL, NULL); - len_1 = get_utg_len(b_1, ug, read_sg, 1, NULL, NULL); - - len_0 = MIN(len_0, len_1); - get_utg_len(b_0, ug, read_sg, 0, &len_0, &thre_0); - get_utg_len(b_1, ug, read_sg, 0, &len_0, &thre_1); - - - for (b_0_i = len_0 = 0, pre_0 = (uint32_t)-1; b_0_i < b_0->b.n; b_0_i++) - { - ori_0 = b_0->b.a[b_0_i]&1; - u_b_0 = &(ug->u.a[b_0->b.a[b_0_i]>>1]); - if(u_b_0->n == 0) continue; - for (b_0_k = 0; b_0_k < u_b_0->n; b_0_k++) - { - len_0++; - if(len_0 > thre_0) return; - if(ori_0 == 1) - { - rId_0 = u_b_0->a[u_b_0->n - b_0_k - 1]>>33; - } - else - { - rId_0 = u_b_0->a[b_0_k]>>33; - } - - if(link->u_idx[rId_0] == (uint32_t)-1) continue; - if(pre_0 == link->u_idx[rId_0]) continue; - pre_0 = link->u_idx[rId_0]; - - for (b_1_i = len_1 = 0, pre_1 = (uint32_t)-1; b_1_i < b_1->b.n; b_1_i++) - { - ori_1 = b_1->b.a[b_1_i]&1; - u_b_1 = &(ug->u.a[b_1->b.a[b_1_i]>>1]); - if(u_b_1->n == 0) continue; - for (b_1_k = 0; b_1_k < u_b_1->n; b_1_k++) - { - len_1++; - if(len_1 > thre_1) goto b_1_i_end; - if(ori_1 == 1) - { - rId_1 = u_b_1->a[u_b_1->n - b_1_k - 1]>>33; - } - else - { - rId_1 = u_b_1->a[b_1_k]>>33; - } - if(link->u_idx[rId_1] == (uint32_t)-1) continue; - if(pre_1 == link->u_idx[rId_1]) continue; - pre_1 = link->u_idx[rId_1]; - - - push_hc_edge(&(link->a.a[pre_0]), pre_1, 1, 1, &d); - push_hc_edge(&(link->a.a[pre_1]), pre_0, 1, 1, &d); - } - } - - b_1_i_end:; - } - } -} - - uint32_t set_utg_offset(buf_t* b, ma_ug_t *ug, asg_t *read_sg, hap_cov_t *cov, uint32_t is_clear) { uint32_t ori, uid, v, nv, l, k; @@ -11810,7 +11718,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, int max_hang, int min_ovlp) asm_opt.hom_global_coverage = -1; 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, 1, cov); + asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, 0, 0, 1, cov); 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); @@ -11909,9 +11817,9 @@ void init_hc_links(hc_links* link, uint64_t ug_num, uint64_t r_num) kv_init(link->a.a[i].e); kv_init(link->a.a[i].f); } - MALLOC(link->u_idx, r_num); - memset(link->u_idx, -1, r_num*sizeof(uint32_t)); - link->r_num = r_num; + // MALLOC(link->u_idx, r_num); + // memset(link->u_idx, -1, r_num*sizeof(uint32_t)); + // link->r_num = r_num; kv_malloc(link->bed, ug_num); link->bed.n = ug_num; for (i = 0; i < link->bed.n; i++) { @@ -11928,7 +11836,7 @@ void destory_hc_links(hc_links* link) kv_destroy(link->a.a[i].f); } kv_destroy(link->a); - free(link->u_idx); + ///free(link->u_idx); for (i = 0; i < link->bed.n; i++) { kv_destroy(link->bed.a[i]); @@ -12061,8 +11969,8 @@ bub_label_t* b_mask_t) init_hc_links(&link, ug->g->n_seq, R_INF.total_reads); asg_t *copy_sg = copy_read_graph(sg); ma_ug_t *copy_ug = copy_untig_graph(ug); - asm_opt.purge_overlap_len = asm_opt.purge_overlap_len_hic; - asm_opt.purge_simi_rate = asm_opt.purge_simi_rate_hic; + ///asm_opt.purge_overlap_len = asm_opt.purge_overlap_len_hic; + ///asm_opt.purge_simi_thres = asm_opt.purge_simi_rate_hic; adjust_utg_by_primary(©_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, max_hang, min_ovlp, &new_rtg_edges, &link, b_mask_t); @@ -12363,7 +12271,7 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) uint32_t print_untig_by_read(ma_ug_t *g, const char* name, uint32_t in, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, const char* info) { - uint32_t i, k, rId = (uint32_t)-1; + uint32_t i, k, rId = (uint32_t)-1, flag = 0; if(in != (uint32_t)-1) { rId = in; @@ -12405,7 +12313,8 @@ ma_hit_t_alloc* reverse_sources, const char* info) { fprintf(stderr, "%s: %s is the %u-th read at %u-th unitig (label: %u, occ: %u)\n", info, name, k, i, g->g->seq[i].c, u->n); - return i; + flag = 1; + ///return i; } } } @@ -12415,7 +12324,7 @@ ma_hit_t_alloc* reverse_sources, const char* info) - fprintf(stderr, "%s: %s is not at any unitig\n", info, name); + if(flag == 0) fprintf(stderr, "%s: %s is not at any unitig\n", info, name); return (uint32_t)-1; } @@ -13810,10 +13719,9 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov) redo: ///print_untig((ug), 61955, "i-0:", 0); - asg_pop_bubble_primary_trio(ug, NULL, trio_flag, DROP, cov); untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, trio_flag, DROP, cov); - magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate); + magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate); ///drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex); /**********debug**********/ if(just_bubble_pop == 0) @@ -13836,7 +13744,6 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov) asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, cov); asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov); asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov); - ///print_debug_gfa(read_g, ug, coverage_cut, "debug_chimeric", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex); ///need consider tangles ///note we need both the read graph and the untig graph @@ -13845,16 +13752,16 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov) cur_cons = get_graph_statistic(g); } untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, trio_flag, DROP, cov); - if(just_bubble_pop == 0) { cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, 2); } + ///print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); + resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, trio_flag, drop_ratio); drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex); all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, DOUBLE_CHECK_THRES); - if(is_first) { is_first = 0; @@ -14591,7 +14498,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp, 0); purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, - asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, + asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, 1, 1, cov); if(asm_opt.recover_atg_cov_min == -1024) { @@ -14612,13 +14519,11 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) } - + adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex); ///primary_flag = get_utg_attributes(*ug, read_g, coverage_cut, sources, ruIndex); update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 0, DOUBLE_CHECK_THRES, flag, drop_rate); - adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex); - nsg = (*ug)->g; n_vtx = nsg->n_seq; for (v = 0; v < n_vtx; ++v) @@ -14671,7 +14576,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) if(asm_opt.purge_level_trio == 1) { purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, - asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, 1, 0, + asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, 1, 0, cov); ///delete_useless_nodes(ug); delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex); @@ -14718,7 +14623,6 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, int is_bench, kvec_asg_arc_t_warp new_rtg_edges; kv_init(new_rtg_edges.a); - ///print_untig_by_read(ug, "m64011_190830_220126/117834372/ccs", 865264, sources, reverse_sources, "beg"); adjust_utg_by_trio(&ug, sg, flag, TRIO_THRES, sources, reverse_sources, coverage_cut, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, max_hang, min_ovlp, &new_rtg_edges, b_mask_t); @@ -19879,7 +19783,7 @@ float drop_ratio) } } - adjust_utg_advance(read_g, src, reverse_sources, ruIndex); + ///adjust_utg_advance(read_g, src, reverse_sources, ruIndex); ///note: we must reset start for each unitig n_vtx = src->g->n_seq; for (v = 0; v < n_vtx; ++v) @@ -19974,7 +19878,6 @@ float drop_ratio) nsu->n = m; } - n_vtx = nsg->n_seq; for (i = 0; i < n_vtx; ++i) { @@ -20108,7 +20011,6 @@ float drop_ratio) free(b_0.b.a); free(b_1.b.a); - ///note: we must reset start for each unitig n_vtx = src->g->n_seq; for (v = 0; v < n_vtx; ++v) @@ -20124,7 +20026,7 @@ float drop_ratio) src->u.a[v].start = UINT32_MAX; } } - adjust_utg_advance(read_g, src, reverse_sources, ruIndex); + ///adjust_utg_advance(read_g, src, reverse_sources, ruIndex); ///note: we must reset start for each unitig n_vtx = src->g->n_seq; for (v = 0; v < n_vtx; ++v) @@ -20133,6 +20035,7 @@ float drop_ratio) if(src->u.a[v].m==0) continue; EvaluateLen(src->u, v) = src->u.a[v].n; } + ///print_untig_by_read(src, "m64076_200203_181219/82511682/ccs", 2429597, NULL, NULL, "end-1"); } @@ -21759,31 +21662,6 @@ R_to_U* ruIndex) } -void reset_reverse_unitigs(hc_links* link, ma_utg_t *u) -{ - uint32_t k = 0, i = 0, rId, pre = (uint32_t)-1; - hc_edge *e = NULL; - if(u->n == 0 || u->m == 0) return; - for (k = 0; k < u->n; k++) - { - rId = u->a[k]>>33; - if(link->u_idx[rId] == (uint32_t)-1) continue; - if(pre == link->u_idx[rId]) continue; - pre = link->u_idx[rId]; - if(link->a.a[pre].f.n == 0) continue; - - for (i = 0; i < link->a.a[pre].f.n; i++) - { - if(link->a.a[pre].f.a[i].del) continue; - link->a.a[pre].f.a[i].del = 1; - e = get_hc_edge(link, link->a.a[pre].f.a[i].uID, pre, 1); - if(e == NULL) continue; - e->del = 1; - } - } -} - - void reset_trans_chain(trans_chain* t_ch, ma_utg_t *u) { uint32_t k = 0, i = 0, p_uId = (uint32_t)-1, c_uId; @@ -21834,7 +21712,6 @@ void append_utg(ma_ug_t* ptg, ma_ug_t* atg, trans_chain* t_ch) for (v = 0; v < atg->g->n_seq; ++v) { if(atg->g->seq[v].del || atg->u.a[v].m == 0) continue; - ///if(link) reset_reverse_unitigs(link, &(atg->u.a[v])); if(t_ch) reset_trans_chain(t_ch, &(atg->u.a[v])); p = &(ptg->u.a[ptg->u.n]); @@ -22070,9 +21947,8 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link, bub_label_t* b_mask_t) ///print_utg_coverage(*ug, coverage_cut, 440, sources); ///exit(0); - drop_semi_circle((*ug), nsg, read_g, reverse_sources, ruIndex); - - asg_cleanup(nsg); + // drop_semi_circle((*ug), nsg, read_g, reverse_sources, ruIndex); + // asg_cleanup(nsg); adjust_utg_advance(read_g, (*ug), reverse_sources, ruIndex); nsg = (*ug)->g; @@ -22094,7 +21970,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link, bub_label_t* b_mask_t) just_contain = 0; if(asm_opt.purge_level_primary == 1) just_contain = 1; purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, - asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, + asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, just_contain, 0, cov); delete_useless_nodes(ug); renew_utg(ug, read_g, new_rtg_edges); @@ -22114,7 +21990,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link, bub_label_t* b_mask_t) just_contain = 0; if(asm_opt.purge_level_primary == 1) just_contain = 1; purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, - asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, + asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, just_contain, 0, cov); delete_useless_nodes(ug); renew_utg(ug, read_g, new_rtg_edges); @@ -22124,7 +22000,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link, bub_label_t* b_mask_t) if(asm_opt.purge_level_primary == 0) { purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, - asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, 0, + asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, 0, 1, cov); } diff --git a/Overlaps.h b/Overlaps.h index d49c203..606b518 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -1052,8 +1052,8 @@ typedef struct{ kvec_t(hc_linkeage) a; kvec_t(uint64_t) enzymes; kvec_t(bed_in) bed; - uint32_t* u_idx; - uint64_t r_num; + ///uint32_t* u_idx; + ///uint64_t r_num; } hc_links; @@ -1103,7 +1103,6 @@ ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link, bub_label_t* b_mask_t); -void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug, asg_t *read_sg); ma_ug_t* copy_untig_graph(ma_ug_t *src); ma_ug_t* output_trio_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, uint8_t flag, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index e6fe0a6..b8f5663 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -4197,6 +4197,7 @@ hap_cov_t *cov) } } +/** void collect_reverse_unitig_pair(hc_links* link, ma_ug_t *ug, hap_overlaps* t) { uint32_t i = 0, k = 0, rId_0, rId_1, pre_0, pre_1, b_0 = t->xUid, b_1 = t->yUid; @@ -4238,6 +4239,7 @@ void collect_reverse_unitigs_purge(buf_t* b_0, hc_links* link, ma_ug_t *ug, hap_ collect_reverse_unitig_pair(link, ug, &(all_ovlp->x[b_0->b.a[k]>>1].a.a[index])); } } +**/ void link_unitigs(asg_t *purge_g, ma_ug_t *ug, hap_overlaps_list* all_ovlp, diff --git a/hic.cpp b/hic.cpp index 61c8304..ad4e13e 100644 --- a/hic.cpp +++ b/hic.cpp @@ -3808,8 +3808,8 @@ void write_hc_links(hc_links* link, const char *fn) fwrite(&link->enzymes.n, sizeof(link->enzymes.n), 1, fp); fwrite(link->enzymes.a, sizeof(uint64_t), link->enzymes.n, fp); - fwrite(&link->r_num, sizeof(link->r_num), 1, fp); - fwrite(link->u_idx, sizeof(uint32_t), 1, fp); + // fwrite(&link->r_num, sizeof(link->r_num), 1, fp); + // fwrite(link->u_idx, sizeof(uint32_t), 1, fp); fwrite(&(link->bed.n), sizeof(link->bed.n), 1, fp); for (k = 0; k < link->bed.n; k++) @@ -3858,9 +3858,9 @@ int load_hc_links(hc_links* link, const char *fn) flag += fread(&link->enzymes.n, sizeof(link->enzymes.n), 1, fp); link->enzymes.m = link->enzymes.n; MALLOC(link->enzymes.a, link->enzymes.n); flag += fread(link->enzymes.a, sizeof(uint64_t), link->enzymes.n, fp); - fread(&link->r_num, sizeof(link->r_num), 1, fp); - MALLOC(link->u_idx, link->r_num); - fread(link->u_idx, sizeof(uint32_t), 1, fp); + // fread(&link->r_num, sizeof(link->r_num), 1, fp); + // MALLOC(link->u_idx, link->r_num); + // fread(link->u_idx, sizeof(uint32_t), 1, fp); diff --git a/hifiasm.1 b/hifiasm.1 index 245a7b1..c4bcb54 100644 --- a/hifiasm.1 +++ b/hifiasm.1 @@ -286,13 +286,14 @@ times in the other sample. .TP 10 .BI -l \ INT Level of purge-dup. 0 to disable purge-dup, 1 to only purge contained haplotigs, -2 to purge all types of haplotigs, 3 to purge all types of haplotigs in most aggressive way. +2 to purge all types of haplotigs, 3 to purge all types of haplotigs in most aggressive way +for high heterozygosity sample. In default, [2] for non-trio assembly, [0] for trio assembly. For trio assembly, only level 0 and level 1 are allowed. .TP .BI -s \ FLOAT -Similarity threshold for duplicate haplotigs that should be purged [0.75]. +Similarity threshold for duplicate haplotigs that should be purged [0.75 for -l1/-l2, 0.55 for -l3]. .TP .BI -O \ FLOAT @@ -303,11 +304,6 @@ Min number of overlapped reads for duplicate haplotigs that should be purged [1] Coverage upper bound of Purge-dups, which is inferred automatically in default. If the coverage of a contig is higher than this bound, don't apply Purge-dups. -.TP -.BI --high-het \ INT -Enable this mode for high heterozygosity sample, which will increase running time. -For ordinary samples, no need to enable this mode [experimental, not stable]. - .SS Debugging options