diff --git a/CommandLines.cpp b/CommandLines.cpp index c9e6eda..5e519f5 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -49,6 +49,8 @@ static ko_longopt_t long_options[] = { { "kpt-rate", ko_required_argument, 335}, { "ul-rate", ko_required_argument, 336}, { "dbg-het-cnt", ko_no_argument, 337}, + { "ul-tip", ko_required_argument, 338}, + { "low-het", ko_no_argument, 339}, { 0, 0, 0 } }; @@ -146,8 +148,10 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " Ultra-Long-integration (beta):\n"); fprintf(stderr, " --ul FILEs file names of Ultra-Long reads [r1.fq,r2.fq,...]\n"); fprintf(stderr, " --ul-rate FLOAT\n"); - fprintf(stderr, " similarity threshold for UL-to-HiFi alignment [%.3g]\n", asm_opt->ul_error_rate); - + fprintf(stderr, " error rate of Ultra-Long reads [%.3g]\n", asm_opt->ul_error_rate); + fprintf(stderr, " --ul-tip INT\n"); + fprintf(stderr, " remove tip unitigs composed of <=INT reads for the UL assembly [%d]\n", asm_opt->max_short_ul_tip); + // fprintf(stderr, " --low-het enable it for genomes with very low het heterozygosity rate (<0.0001%%)\n"); fprintf(stderr, "Example: ./hifiasm -o NA12878.asm -t 32 NA12878.fq.gz\n"); fprintf(stderr, "See `https://hifiasm.readthedocs.io/en/latest/' or `man ./hifiasm.1' for complete documentation.\n"); @@ -204,6 +208,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->min_overlap_Len = 50; asm_opt->min_overlap_coverage = 0; asm_opt->max_short_tip = 3; + asm_opt->max_short_ul_tip = 6; asm_opt->min_cnt = 2; asm_opt->mid_cnt = 5; asm_opt->purge_level_primary = 3; @@ -246,6 +251,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->ul_error_rate_hpc = 0.2; asm_opt->ul_ec_round = 3; asm_opt->is_dbg_het_cnt = 0; + asm_opt->is_low_het_ul = 0; } void destory_enzyme(enzyme* f) @@ -763,6 +769,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 335) asm_opt->kpt_rate = atof(opt.arg); else if (c == 336) asm_opt->ul_error_rate = atof(opt.arg); else if (c == 337) asm_opt->is_dbg_het_cnt = 1; + else if (c == 338) asm_opt->max_short_ul_tip = atol(opt.arg); + else if (c == 339) asm_opt->is_low_het_ul = 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 eb062e8..013c611 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.18.1-r466" +#define HA_VERSION "0.18.2-r467" #define VERBOSE 0 @@ -79,6 +79,7 @@ typedef struct { int min_overlap_Len; int min_overlap_coverage; int max_short_tip; + int max_short_ul_tip; int min_cnt; int mid_cnt; int purge_level_primary; @@ -128,6 +129,7 @@ typedef struct { double ul_error_rate, ul_error_rate_low, ul_error_rate_hpc; int32_t ul_ec_round; uint8_t is_dbg_het_cnt; + uint8_t is_low_het_ul; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Overlaps.cpp b/Overlaps.cpp index f86158a..95c7eb6 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -3594,7 +3594,7 @@ uint32_t startNode, uint32_t endNode, int max_dist, buf_t* bub) { Ns_first[0] = av[0].v; Ns_first[1] = av[1].v; f1 = asg_bub_end_finder_with_del_advance(g, Ns_first, 2, max_dist, - bub, 0, (u_int32_t)-1, &convex1); + bub, 0, (uint32_t)-1, &convex1); l1 = min_thres + 10; } @@ -3602,7 +3602,7 @@ uint32_t startNode, uint32_t endNode, int max_dist, buf_t* bub) { Ns_second[0] = aw[0].v; Ns_second[1] = aw[1].v; f2 = asg_bub_end_finder_with_del_advance(g, Ns_second, 2, max_dist, - bub, 0, (u_int32_t)-1, &convex2); + bub, 0, (uint32_t)-1, &convex2); l2 = min_thres + 10; } @@ -4034,7 +4034,7 @@ uint32_t startNode, uint32_t endNode, int max_dist, buf_t* bub) { Ns_first[0] = Nodes1[0]^1; Ns_first[1] = Nodes2[1]^1; f1 = asg_bub_end_finder_with_del_advance(g, Ns_first, 2, max_dist, - bub, 0, (u_int32_t)-1, &convex1); + bub, 0, (uint32_t)-1, &convex1); l1 = min_thres + 10; } @@ -4042,7 +4042,7 @@ uint32_t startNode, uint32_t endNode, int max_dist, buf_t* bub) { Ns_second[0] = Nodes2[0]^1; Ns_second[1] = Nodes1[1]^1; f2 = asg_bub_end_finder_with_del_advance(g, Ns_second, 2, max_dist, - bub, 0, (u_int32_t)-1, &convex2); + bub, 0, (uint32_t)-1, &convex2); l2 = min_thres + 10; } @@ -6896,7 +6896,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex) kv_push(uint32_t, b_r, w); aw = asg_arc_a(g, w); - min_edge = (u_int32_t)-1; + min_edge = (uint32_t)-1; for (t = 0; t < nw; t++) { if(aw[t].del) continue; @@ -6996,7 +6996,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex) } //check the length if(to_del_l == 0 && asg_bub_end_finder_with_del_advance(g, - b_f.a, b_f.n, max_dist, &bub, 0, (u_int32_t)-1, &sink)==1) + b_f.a, b_f.n, max_dist, &bub, 0, (uint32_t)-1, &sink)==1) { to_del_l = 1; } @@ -7143,7 +7143,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex) kv_push(uint32_t, b_r, w); aw = asg_arc_a(g, w); - min_edge = (u_int32_t)-1; + min_edge = (uint32_t)-1; for (t = 0; t < nw; t++) { if(aw[t].del) continue; @@ -7187,7 +7187,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex) kv_push(uint32_t, b_r, w); aw = asg_arc_a(g, w); - min_edge = (u_int32_t)-1; + min_edge = (uint32_t)-1; for (t = 0; t < nw; t++) { if(aw[t].del) continue; @@ -7288,7 +7288,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex) } //check the length if(to_del_l == 0 && asg_bub_end_finder_with_del_advance(g, - b_f.a, b_f.n, max_dist, &bub, 0, (u_int32_t)-1, &sink)==1) + b_f.a, b_f.n, max_dist, &bub, 0, (uint32_t)-1, &sink)==1) { to_del_l = 1; } @@ -7431,7 +7431,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex) kv_push(uint32_t, b_r, w); aw = asg_arc_a(g, w); - min_edge = (u_int32_t)-1; + min_edge = (uint32_t)-1; for (t = 0; t < nw; t++) { if(aw[t].del) continue; @@ -7475,7 +7475,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex) kv_push(uint32_t, b_r, w); aw = asg_arc_a(g, w); - min_edge = (u_int32_t)-1; + min_edge = (uint32_t)-1; for (t = 0; t < nw; t++) { if(aw[t].del) continue; @@ -7575,7 +7575,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex) } //check the length if(to_del_l == 0 && asg_bub_end_finder_with_del_advance(g, - b_f.a, b_f.n, max_dist, &bub, 0, (u_int32_t)-1, &sink)==1) + b_f.a, b_f.n, max_dist, &bub, 0, (uint32_t)-1, &sink)==1) { to_del_l = 1; } @@ -31884,7 +31884,7 @@ int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file) nopt.cov = coverage_cut; nopt.b_mask_t = b_mask_t; nopt.max_hang = max_hang_length; nopt.mini_ovlp = mini_overlap_length; ul_realignment_gfa(uopt, *sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, - asm_opt.max_short_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file, &nopt);; + asm_opt.max_short_tip, asm_opt.max_short_ul_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file, &nopt);; // ma_ug_t *iug = ul_realignment_gfa(uopt, *sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, // asm_opt.max_short_tip, b_mask_t, ha_opt_triobin(&asm_opt), o_file); // asg_t *ng = gen_ng(iug, *sg, uopt, coverage_cut, ruIndex, 256); diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 9f997d8..7892c67 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -14385,9 +14385,12 @@ void u2g_hybrid_detan_iter(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, uint asg64_v b64, ub64; kv_init(b64); kv_init(ub64); asg64_v tx = {0,0,0}, tb = {0,0,0}, *ob = NULL, *ub = NULL; ob = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); ob->n = ub->n = 0; + // fprintf(stderr, "\n[M::%s::] asm_opt.is_low_het_ul::%u, max_ext::%u\n", + // __func__, asm_opt.is_low_het_ul, max_ext); // prt_usg_t(uidx, ng, "ng0"); for (k = 0; k < clean_round; k++) { ncut += ug_ext_strict(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64); + // if(asm_opt.is_low_het_ul) break; ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 48, 1); // prt_usg_t(uidx, ng, "ng_python"); ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 16, 1); @@ -14403,7 +14406,7 @@ void u2g_hybrid_detan_iter(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, uint // prt_usg_t(uidx, ng, "ng.db"); // ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 10); // ncut += ug_ext_strict(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64); - // prt_usg_t(uidx, ng, "ng1"); + // prt_usg_t(uidx, ng, "ng.db"); if(ncut) { usg_cleanup(ng); usg_arc_cut_tips(ng, max_ext, 1, ub); } @@ -16479,8 +16482,38 @@ void destroy_ul_resolve_t(ul_resolve_t *uidx) free(uidx); } +static void clear_ma_hit_t_alloc(void *data, long i, int tid) +{ + ma_hit_t_alloc *src = (ma_hit_t_alloc *)data; + ma_hit_t_alloc *z = &(src[i]); uint32_t k; + for (k = 0; k < z->length; k++) z->buffer[k].del = 0; +} + +static void reset_ma_sub_t(void *data, long i, int tid) +{ + sset_aux *s = (sset_aux *)data; + if(s->g && s->g->seq[i].del) s->cov[i].del = 1; + else s->cov[i].del = 0; +} + +void renew_R_to_U(asg_t *ng, ma_hit_t_alloc* src, ma_hit_t_alloc* r_src, int64_t n_read, ma_sub_t *coverage_cut, +R_to_U* ruIndex, int max_hang, int min_ovlp) +{ + sset_aux s; memset(&s, 0, sizeof(s)); + kt_for(asm_opt.thread_num, clear_ma_hit_t_alloc, src, n_read); + kt_for(asm_opt.thread_num, clear_ma_hit_t_alloc, r_src, n_read); + + s.g = NULL; s.cov = coverage_cut; + kt_for(asm_opt.thread_num, reset_ma_sub_t, &s, n_read); + + ma_hit_contained_advance(src, n_read, coverage_cut, ruIndex, max_hang, min_ovlp); + + // s.g = ng; s.cov = coverage_cut; + // kt_for(asm_opt.thread_num, reset_ma_sub_t, &s, n_read); +} + void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio, -double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, +double max_ovlp_drop_ratio, int64_t max_tip, int64_t max_ul_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, ul_renew_t *ropt) { uint64_t i; uint8_t *r_het = NULL; bubble_type *bub = NULL; ulg_opt_t uu; @@ -16502,7 +16535,7 @@ ul_renew_t *ropt) // print_raw_uls_seq(uidx, asm_opt.output_file_name); // print_raw_uls_aln(uidx, asm_opt.output_file_name); ul_re_correct(uidx, 3); - init_ulg_opt_t(&uu, uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.55, max_tip, max_tip<<1, b_mask_t, is_trio); + init_ulg_opt_t(&uu, uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.55, max_tip, max_ul_tip, b_mask_t, is_trio); // print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug0", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1); /**ul2ul_idx_t *u2o = **/gen_ul2ul(uidx, uopt, &uu, 0); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-3"); @@ -16529,7 +16562,7 @@ ul_renew_t *ropt) asg_destroy((*(ropt->sg))); (*(ropt->sg)) = ng; (*(ropt->src)) = R_INF.paf; (*(ropt->r_src)) = R_INF.reverse_paf; (*(ropt->n_read)) = R_INF.total_reads; (*(ropt->readLen)) = R_INF.read_length; - ma_hit_contained_advance((*(ropt->src)), (*(ropt->n_read)), (*(ropt->cov)), ropt->ruIndex, ropt->max_hang, ropt->mini_ovlp); + renew_R_to_U(ng, (*(ropt->src)), (*(ropt->r_src)), (*(ropt->n_read)), (*(ropt->cov)), ropt->ruIndex, ropt->max_hang, ropt->mini_ovlp); post_rescue(uopt, (*(ropt->sg)), (*(ropt->src)), (*(ropt->r_src)), ropt->ruIndex, ropt->b_mask_t, 0); // print_raw_uls_aln(uidx, asm_opt.output_file_name); // exit(0); diff --git a/gfa_ut.h b/gfa_ut.h index 8fc2274..5823e21 100644 --- a/gfa_ut.h +++ b/gfa_ut.h @@ -6,6 +6,7 @@ typedef struct { asg_t *g; ma_hit_t_alloc *src; + ma_sub_t *cov; R_to_U* ruIndex; int64_t max_hang; int64_t min_ovlp; @@ -25,7 +26,7 @@ void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float o uint32_t asg_cut_large_indel(asg_t *g, asg64_v *in, int32_t max_ext, float ou_rat, uint32_t is_ou); uint32_t asg_cut_semi_circ(asg_t *g, uint32_t lim_len, uint32_t is_clean); void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio, -double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, ul_renew_t *ropt); +double max_ovlp_drop_ratio, int64_t max_tip, int64_t max_ul_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file, ul_renew_t *ropt); void recover_contain_g(asg_t *g, ma_hit_t_alloc *src, R_to_U* ruIndex, int64_t max_hang, int64_t min_ovlp, int64_t ul_occ); void normalize_gou(asg_t *g); void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd); diff --git a/hic.h b/hic.h index 0f2ffda..a6bc80c 100644 --- a/hic.h +++ b/hic.h @@ -110,4 +110,6 @@ void set_utg_by_dis(uint32_t v, pdq* pq, asg_t *g, kvec_t_u32_warp *res, uint32_ void dedup_hits(kvec_pe_hit* hits, uint64_t is_dup); void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, uint32_t is_poy, kvec_pe_hit **rhits); spg_t *hic_pre_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, kvec_pe_hit **rhits); +void prt_bubble_gfa_adv(FILE *fp, bubble_type *bub, const char* utg_pre, const char* bub_pre, const char* chain_pre); + #endif