diff --git a/CommandLines.h b/CommandLines.h index be50045..baea653 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.19.3-r572" +#define HA_VERSION "0.19.4-r575" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index 0372f9d..2dc3f08 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -15208,11 +15208,10 @@ void purge_ovlp_cov_adv(uint32_t id, kv_u_trans_t *ta, asg64_v *b64, ug_rid_cov_ } } -void trans_sec_cut_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc* src) +void trans_sec_cut_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc* src, ug_rid_cov_t *in) { - uint64_t k; - asg64_v b64; kv_init(b64); - ug_rid_cov_t *cc = gen_ug_rid_cov_t(ug, sg, src); + uint64_t k; asg64_v b64; kv_init(b64); + ug_rid_cov_t *cc = ((in)?(in):(gen_ug_rid_cov_t(ug, sg, src))); fprintf(stderr, "+[M::%s]\thom_cov::%lu\thet_cov::%lu\thom_cut::%lu\n", __func__, cc->hom_cov, cc->het_cov, cc->hom_max); @@ -15221,7 +15220,7 @@ void trans_sec_cut_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *sg, ma purge_ovlp_cov_adv(k, ta, &b64, cc, cc->hom_max); } - destory_ug_rid_cov_t(cc); free(cc); kv_destroy(b64); + if(!in) {destory_ug_rid_cov_t(cc); free(cc);} kv_destroy(b64); } static void worker_for_trans_sec_simple_cut(void *data, long i, int tid) // callback for kt_for() @@ -15251,8 +15250,18 @@ static void worker_for_trans_sec_simple_cut(void *data, long i, int tid) // call } } +void gen_ug_rid_cov_t_by_ovlp(kv_u_trans_t *ta, ug_rid_cov_t *cc) +{ + uint64_t z, k, id, n; u_trans_t *a; + for (z = 0; z < cc->ug->g->n_seq; z++) { + id = z; a = u_trans_a(*ta, id); n = u_trans_n(*ta, id); + for (k = 0; k < n; k++) { + append_cov_line_ug_rid_cov_t(id, cc->cov.a+cc->idx[id], &(a[k]), cc, ((uint64_t)-1), -1); + } + } +} -void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* src) +void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* src, ug_rid_cov_t *in) { u_trans_clean_t sl; uint64_t k, i, l, st, occ; ha_mzl_t *tz; ha_mzl_v srt_a; kv_u_trans_t *bl; u_trans_t *z; @@ -15318,7 +15327,7 @@ void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t * free(sl.res); kv_destroy(srt_a); kt_u_trans_t_idx(ta, ug->g->n_seq); // dbg_prt_utg_trans(ta, ug, "after"); - trans_sec_cut_filter_mmhap_adv(ta, ug, read_g, src); + trans_sec_cut_filter_mmhap_adv(ta, ug, read_g, src, in); kt_for(sl.n_thread, worker_for_trans_sec_simple_cut, &sl, sl.ta->idx.n); // CALLOC(sl.srt, sl.n_thread); sl.sec_rate = 0.5; @@ -15338,6 +15347,7 @@ void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t * st = i; } } + // dbg_prt_utg_extra_trans(ta, ug, asm_opt.output_file_name); } @@ -15660,9 +15670,9 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, long long gap // if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_trans_chain(cov->t_ch, output_file_name); } - dbg_prt_utg_trans(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, "pre"); - clean_u_trans_t_idx_filter_mmhap_adv(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, sg, opt->sources); - dbg_prt_utg_trans(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, "after"); + // dbg_prt_utg_trans(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, "pre"); + clean_u_trans_t_idx_filter_mmhap_adv(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, sg, opt->sources, NULL); + // dbg_prt_utg_trans(&(cov?cov->t_ch->k_trans:t_ch->k_trans), ug, "after"); // refine_hic_trans_mmhap(opt, &(cov?cov->t_ch->k_trans:t_ch->k_trans), sg, ug); ///for debug @@ -15723,7 +15733,7 @@ float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, long long gap kv_destroy(d_edges.a); asg_cleanup(sg); - dbg_prt_trio_mmhap_label(ug, rh, output_file_name); + // dbg_prt_trio_mmhap_label(ug, rh, output_file_name); output_trio_mmhap(sg, coverage_cut, output_file_name, sources, reverse_sources, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, max_hang, min_ovlp, gap_fuzz, b_mask_t, opt, ug, rh, asm_opt.polyploidy); @@ -18596,46 +18606,62 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov) { asg_t *g = ug->g; uint32_t is_first = 1; - + // if(trio_flag == MOTHER) { + // print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len, 0, 1, 0); + // exit(1); + // } redo: - ///print_untig((ug), 61955, "i-0:", 0); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-0:", 0); // fprintf(stderr, "[M::%s] 0\n", __func__); asg_pop_bubble_primary_trio(ug, NULL, trio_flag, DROP, cov, NULL, 1); + + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-1:", 0); // fprintf(stderr, "[M::%s] 1\n", __func__); magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate); // fprintf(stderr, "[M::%s] 2\n", __func__); + + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-2:", 0); /**********debug**********/ if(just_bubble_pop == 0) { cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2); } // fprintf(stderr, "[M::%s] 3\n", __func__); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-3:", 0); /**********debug**********/ long long pre_cons = get_graph_statistic(g); long long cur_cons = 0; while(pre_cons != cur_cons) { // fprintf(stderr, "[M::%s] 4\n", __func__); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-4:", 0); pre_cons = get_graph_statistic(g); // fprintf(stderr, "[M::%s] 5\n", __func__); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-5:", 0); ///need consider tangles asg_pop_bubble_primary_trio(ug, NULL, trio_flag, DROP, cov, NULL, 1); // fprintf(stderr, "[M::%s] 6\n", __func__); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-6:", 0); /**********debug**********/ if(just_bubble_pop == 0) { ///need consider tangles asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, trio_flag, cov, NULL); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-7:", 0); // fprintf(stderr, "[M::%s] 7\n", __func__); // if(trio_flag == MOTHER) print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, cov, NULL); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-8:", 0); // fprintf(stderr, "[M::%s] 8\n", __func__); asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, cov, NULL, trio_flag); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-9:", 0); // fprintf(stderr, "[M::%s] 9\n", __func__); asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, cov, NULL, trio_flag); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-10:", 0); // fprintf(stderr, "[M::%s] 10\n", __func__); detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex, NULL, cov->is_r_het); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-11:", 0); // fprintf(stderr, "[M::%s] 11\n", __func__); ///need consider tangles ///note we need both the read graph and the untig graph @@ -18643,34 +18669,42 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov) /**********debug**********/ cur_cons = get_graph_statistic(g); // fprintf(stderr, "[M::%s] 12\n", __func__); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-12:", 0); } if(just_bubble_pop == 0) { // fprintf(stderr, "[M::%s] 13\n", __func__); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-13:", 0); cut_trio_tip_primary(g, ug, tipsLen, trio_flag, 0, read_g, reverse_sources, ruIndex, cov->is_r_het, 2); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-14:", 0); // fprintf(stderr, "[M::%s] 14\n", __func__); } // print_debug_gfa(read_g, ug, coverage_cut, "debug_dups", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); // fprintf(stderr, "[M::%s] 15\n", __func__); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-15:", 0); magic_trio_phasing(g, ug, read_g, coverage_cut, sources, reverse_sources, 2, ruIndex, trio_flag, trio_drop_rate); // fprintf(stderr, "[M::%s] 16\n", __func__); - + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-16:", 0); // print_debug_gfa(read_g, ug, coverage_cut, "resolve_tangles", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len, 0, 0, 0); // exit(1); ///bug here resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, cov->is_r_het, trio_flag, drop_ratio); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-17:", 0); // fprintf(stderr, "[M::%s] 17\n", __func__); drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex, cov->is_r_het); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-18:", 0); // fprintf(stderr, "[M::%s] 18\n", __func__); all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, cov->is_r_het, DOUBLE_CHECK_THRES, asm_opt.trio_flag_occ_thres); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-19:", 0); // fprintf(stderr, "[M::%s] 19\n", __func__); // if(trio_flag == MOTHER) print_untig_by_read(ug, "m54329U_190827_173812/30214441/ccs", (uint32_t)-1, NULL, NULL, "bf-16"); if(is_first) { is_first = 0; unitig_arc_del_short_diploid_by_length(ug->g, drop_ratio); + // if(trio_flag == MOTHER) print_untig((ug), 425, "i-20:", 0); // fprintf(stderr, "[M::%s] 20\n", __func__); goto redo; } @@ -19332,16 +19366,105 @@ void delete_useless_nodes(ma_ug_t **ug) asg_cleanup(nsg); } +inline uint32_t is_useful_node(uint32_t flag_occ, uint32_t non_flag_occ, uint32_t drop_occ, uint32_t tot_occ, +float flag_rate, float used_rate, uint32_t min_occ) +{ + if((flag_occ > 0) && (flag_occ >= min_occ) && (flag_occ >= ((non_flag_occ+flag_occ+drop_occ)*flag_rate)) + && (drop_occ <= (tot_occ*used_rate))) { + return 1; + } + return 0; +} +void recover_chain_nodes(buf_t *in, ma_ug_t *ug, uint32_t flag, float flag_rate, float used_rate, uint32_t min_occ) +{ + ma_utg_t *u = NULL; uint32_t flag_occ, non_flag_occ, drop_occ, rid; + uint32_t tot_flag_occ, tot_non_flag_occ, tot_drop_occ, tot_occ, i, k, z; + tot_flag_occ = tot_non_flag_occ = tot_drop_occ = tot_occ = 0; + + for (i = 0; i < in->b.n; i++) { + u = &(ug->u.a[in->b.a[i]>>1]); + flag_occ = non_flag_occ = drop_occ = 0; + for (k = 0; k < u->n; k++) { + rid = u->a[k]>>33; + if(R_INF.trio_flag[rid] == AMBIGU) continue; + else if(R_INF.trio_flag[rid] == DROP) drop_occ++; + else if(R_INF.trio_flag[rid] == flag) flag_occ++; + else if(R_INF.trio_flag[rid] != flag) non_flag_occ++; + } + + if(is_useful_node(flag_occ, non_flag_occ, drop_occ, u->n, flag_rate, used_rate, min_occ)) { + ug->g->seq[in->b.a[i]>>1].c = PRIMARY_LABLE; + } + tot_flag_occ += flag_occ; + tot_non_flag_occ += non_flag_occ; + tot_drop_occ += drop_occ; + tot_occ += u->n; + if(is_useful_node(tot_flag_occ, tot_non_flag_occ, tot_drop_occ, tot_occ, flag_rate, used_rate, min_occ)) { + for (z = 0; z <= i; z++) { + if(ug->g->seq[in->b.a[z]>>1].c == ALTER_LABLE) { + u = &(ug->u.a[in->b.a[z]>>1]); + flag_occ = non_flag_occ = drop_occ = 0; + for (k = 0; k < u->n; k++) { + rid = u->a[k]>>33; + if(R_INF.trio_flag[rid] == AMBIGU) continue; + else if(R_INF.trio_flag[rid] == DROP) drop_occ++; + else if(R_INF.trio_flag[rid] == flag) flag_occ++; + else if(R_INF.trio_flag[rid] != flag) non_flag_occ++; + } + if(is_useful_node(flag_occ, non_flag_occ, drop_occ, u->n, flag_rate, used_rate, 0/**min_occ**/)) { + ug->g->seq[in->b.a[z]>>1].c = PRIMARY_LABLE; + } + } + } + } + } +} + +void rescue_useless_trio_nodes(ma_ug_t *ug, uint32_t flag, float flag_rate, float used_rate, uint32_t min_occ) +{ + asg_t* nsg = ug->g; + uint32_t v, w, n_vtx = nsg->n_seq<<1, i, k, z; + long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen; + buf_t b; memset(&b, 0, sizeof(buf_t)); + + for (v = 0; v < n_vtx; v++) { + if(nsg->seq[v>>1].del) continue; + if(nsg->seq[v>>1].c != ALTER_LABLE) continue; + ///check if beg is the tig end + if(get_real_length(nsg, v^1, NULL) == 1) { + get_real_length(nsg, v^1, &w); + if(get_real_length(nsg, w^1, NULL) == 1) continue; + } + + b.b.n = 0; + get_unitig(nsg, ug, v, &w, &nodeLen, &baseLen, &max_stop_nodeLen, &max_stop_baseLen, 1, &b); + recover_chain_nodes(&b, ug, flag, flag_rate, used_rate, min_occ); + if(b.b.n > 1) { + k = b.b.n>>1; + for (i = 0; i < k; i++) { + z = b.b.a[i]; b.b.a[i] = b.b.a[b.b.n-i-1]; b.b.a[b.b.n-i-1] = z; + } + recover_chain_nodes(&b, ug, flag, flag_rate, used_rate, min_occ); + } + } + + free(b.b.a); + asg_cleanup(nsg); +} void delete_useless_trio_nodes(ma_ug_t **ug, asg_t* read_g, ma_sub_t* coverage_cut, -ma_hit_t_alloc* sources, R_to_U* ruIndex) +ma_hit_t_alloc* sources, R_to_U* ruIndex, uint32_t flag, float flag_rate, float used_rate, uint32_t min_occ) { asg_t* nsg = (*ug)->g; uint32_t v, n_vtx = nsg->n_seq; uint8_t* primary_flag = (uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t)); + if(flag_rate > 0 && used_rate > 0 && min_occ > 0) { + rescue_useless_trio_nodes(*ug, flag, flag_rate, used_rate, min_occ); + } + for (v = 0; v < n_vtx; ++v) { if(nsg->seq[v].del) continue; @@ -19630,7 +19753,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) // fprintf(stderr, "[M::%s] 4\n", __func__); ///delete_useless_nodes(ug); - delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex); + delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex, flag, 0.8, 0.15, 16); // fprintf(stderr, "[M::%s] 5\n", __func__); update_hap_label(*ug, read_g); @@ -19678,7 +19801,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) // fprintf(stderr, "[M::%s] 16\n", __func__); ///delete_useless_nodes(ug); - delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex); + delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex, flag, 0.8, 0.15, 16); // fprintf(stderr, "[M::%s] 17\n", __func__); @@ -19688,7 +19811,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, 1, 0, cov, 0, 0); ///delete_useless_nodes(ug); - delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex); + delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex, flag, 0.8, 0.15, 16); } // fprintf(stderr, "[M::%s] 18\n", __func__); diff --git a/Overlaps.h b/Overlaps.h index 65d116e..cb99f5e 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -1233,5 +1233,7 @@ void destory_ug_rid_cov_t(ug_rid_cov_t *p); uint32_t append_cov_line_ug_rid_cov_t(uint64_t uid, uint64_t *qcc, u_trans_t *p, ug_rid_cov_t *idx, uint64_t hom_cut, double cut_rate); uint64_t infer_mmhap_copy(ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc *src, uint8_t *ff, uint64_t uid, uint64_t het_cov, uint64_t n_hap); uint64_t trans_sec_cut0(kv_u_trans_t *ta, asg64_v *srt, uint32_t id, double sec_rate, uint64_t bd, ma_ug_t *ug); +void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* src, ug_rid_cov_t *in); +void gen_ug_rid_cov_t_by_ovlp(kv_u_trans_t *ta, ug_rid_cov_t *cc); #endif diff --git a/README.md b/README.md index c081c5a..b3babde 100644 --- a/README.md +++ b/README.md @@ -43,11 +43,8 @@ See [tutorial][tutorial] for more details. ## Introduction -Hifiasm is a fast haplotype-resolved de novo assembler for PacBio HiFi reads. -It can assemble a human genome in several hours and assemble a ~30Gb California -redwood genome in a few days. Hifiasm emits partially phased assemblies of -quality competitive with the best assemblers. Given parental short reads or -Hi-C data, it produces arguably the best haplotype-resolved assemblies so far. +Hifiasm is a fast haplotype-resolved de novo assembler initially designed for PacBio HiFi reads. +Its latest release could support the telomere-to-telomere assembly by utilizing ultralong Oxford Nanopore reads. Hifiasm produces arguably the best single-sample telomere-to-telomere assemblies combing HiFi, ultralong and Hi-C reads, and it is one of the best haplotype-resolved assemblers for the trio-binning assembly given parental short reads. For a human genome, hifiasm can produce the telomere-to-telomere assembly in one day. ## Why Hifiasm? @@ -146,11 +143,18 @@ The second command line will run much faster than the first. ### Ultra-long ONT integration -Hifiasm could integrate ultra-long ONT reads to improve the assembly quality: +Hifiasm could integrate ultra-long ONT reads to produce the telomere-to-telomere assembly: ```sh hifiasm -o NA12878.asm -t32 --ul ul.fq.gz HiFi-reads.fq.gz ``` -Please note that this mode is not stable right now. We have only tested with >=100kb UL reads. +For the single-sample telomere-to-telomere assembly with Hi-C reads: +```sh +hifiasm -o NA12878.asm -t32 --ul ul.fq.gz --h1 read1.fq.gz --h2 read2.fq.gz HiFi-reads.fq.gz +``` +For the trio-binning telomere-to-telomere assembly； +```sh +hifiasm -o NA12878.asm -t32 --ul ul.fq.gz -1 pat.yak -2 mat.yak HiFi-reads.fq.gz +``` ### Output files diff --git a/inter.cpp b/inter.cpp index 2cf8218..b093f6a 100644 --- a/inter.cpp +++ b/inter.cpp @@ -10434,6 +10434,202 @@ static void worker_for_trans_ovlp_mmhap_adv(void *data, long i, int tid) // call s->free_cnt[tid]++; } +uint32_t tranfor_ovlp(u_trans_t *qovlp, u_trans_t *tovlp, asg_t *g, ul_ov_t *res, uint32_t adjust_rev) +{ + int64_t os, oe, s_shift, e_shift, tt, qs, qe, ts, te; + os = MAX(qovlp->ts, tovlp->qs); + oe = MIN(qovlp->te, tovlp->qe); + if(oe <= os) return 0; + + ///[os, oe) -> qovlp->t* + s_shift = get_offset_adjust(os-qovlp->ts, qovlp->te-qovlp->ts, qovlp->qe-qovlp->qs); + e_shift = get_offset_adjust(qovlp->te-oe, qovlp->te-qovlp->ts, qovlp->qe-qovlp->qs); + if(qovlp->rev) { + tt = s_shift; s_shift = e_shift; e_shift = tt; + } + qs = qovlp->qs+s_shift; qe = ((int64_t)qovlp->qe)-e_shift; + if(qs >= qe) return 0; + + ///[os, oe) -> tovlp->q* + s_shift = get_offset_adjust(os-tovlp->qs, tovlp->qe-tovlp->qs, tovlp->te-tovlp->ts); + e_shift = get_offset_adjust(tovlp->qe-oe, tovlp->qe-tovlp->qs, tovlp->te-tovlp->ts); + if(tovlp->rev) { + tt = s_shift; s_shift = e_shift; e_shift = tt; + } + ts = tovlp->ts+s_shift; te = ((int64_t)tovlp->te)-e_shift; + if(ts >= te) return 0; + + memset(res, 0, sizeof(*res)); + res->qn = qovlp->qn; res->qs = qs; res->qe = qe; + res->tn = tovlp->tn; res->ts = ts; res->te = te; + res->rev = ((qovlp->rev == tovlp->rev)?0:1); + if(adjust_rev && res->rev) {///for linear chaining + res->ts = g->seq[res->tn].len - te; + res->te = g->seq[res->tn].len - ts; + } + return 1; +} + +uint32_t rescue_adject_ovlp(asg_t *g, uint32_t id, kv_u_trans_t *ta, kv_ul_ov_t *out, st_mt_t *buf) +{ + u_trans_t *a, *b; ul_ov_t rr; uint64_t a_n, b_n, k, l, i, z, m; + a = u_trans_a((*ta), id); a_n = u_trans_n((*ta), id); + for (i = out->n = buf->n = 0; i < a_n; i++) { + b = u_trans_a((*ta), a[i].tn); b_n = u_trans_n((*ta), a[i].tn); + z = a[i].tn; z <<= 32; kv_push(uint64_t, *buf, z); + for (k = 0; k < b_n; k++) { + if(b[k].tn == id) continue; + if(!tranfor_ovlp(&(a[i]), &(b[k]), g, &rr, 1)) continue; + z = rr.tn; z <<= 32; z |= out->n; z |= ((uint64_t)0x80000000); + rr.tn <<= 1; rr.tn |= rr.rev; kv_push(ul_ov_t, *out, rr); + } + } + if(out->n == 0) return 1; + + radix_sort_gfa64(buf->a, buf->a + buf->n); + for (k = 1, l = m = 0; k <= buf->n; k++) { + if(k == buf->n || (buf->a[l]>>32)!=(buf->a[k]>>32)) { + if((k - l > 1) && (!(buf->a[l]&((uint64_t)0x80000000)))) {///overlap within bck + for (z = l; z < k; z++) { + if(buf->a[z]&((uint64_t)0x80000000)) { + out->a[(uint32_t)(buf->a[z]-((uint64_t)0x80000000))].tn = (uint32_t)-1; + m++; + } + } + } + l = k; + } + } + + if(m) { + for (k = m = 0; k < out->n; k++) { + if(out->a[k].tn == (uint32_t)-1) continue; + out->a[m++] = out->a[k]; + } + out->n = m; + } + if(out->n == 0) return 1; + + radix_sort_ul_ov_srt_tn(out->a, out->a+out->n); + for (k = 0; k < out->n; k++) out->a[k].tn >>= 1; + + return 0; +} + +/** +uint64_t gen_trans_chain_mmhap(ug_trans_t *s, uint64_t rid, ha_ovec_buf_t *b, kv_ul_ov_t *bl, char *seq, uint64_t len, +double err, double bw) +{ + uint64_t cnt = ((s->idx_n.a[rid+1]-s->idx_n.a[rid])), ol_h = 0, pass_aln = 0; + uint32_t high_occ = asm_opt.polyploidy + 1; overlap_region *aux_o = NULL; + ///note: high_occ is different + ug_map_lchain(b->abl, rid, seq, len, s->w, s->k, &(s->udb), &b->olist, &b->clist, bw, bw, + s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3, + s->is_HPC, s->idx_a.a + s->idx_n.a[rid], cnt, s->srt_a.a, s->srt_a.n, s->mini_cut, s->chain_cut, NULL); + // if(rid == 57) { + // fprintf(stderr, "-1-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length); + // } + ///remove candidate chains that have been calculated + if(!fi) backward_dedup_ol(rid, bl, &(b->sp), &b->olist);///it is ok + // if(rid == 57) { + // fprintf(stderr, "-2-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length); + // } + filter_by_reliable_ovlp_mmhap_adv(rid, s->filter, &(b->sp), &b->olist, &(s->udb), s->sec_cutoff, 1, 1, s->ccov, &ol_h); + clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2); + if(!fi) ol_h = 0; + + // if(rid == 57) { + // fprintf(stderr, "-3-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length); + // } + + ol_h = split_ug_lalign(ol_h, &b->olist, err_high, err_low, + &b->clist, &(s->udb), s->uopt, seq, len, &b->self_read, &b->ovlp_read, + &b->correct, &b->exz, aux_o, rid, s->k, s->chain_cut, NULL); + + // if(rid == 57) { + // fprintf(stderr, "-4-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length); + // } + + aux_o = gen_aux_ovlp(&b->olist);///must be here + + // if(rid == 57) { + // fprintf(stderr, "-5-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length); + // } + + ol_h = split_ug_lalign(ol_h, &b->olist, err_high, err_low, + &b->clist, &(s->udb), s->uopt, seq, len, &b->self_read, &b->ovlp_read, + &b->correct, &b->exz, aux_o, rid, s->k, s->chain_cut, NULL); + + // if(rid == 57) { + // fprintf(stderr, "-6-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length); + // } + + if(fi) {///first round + pass_aln = test_het_aln_mmhap(rid, s->ccov, u_trans_a((*(s->filter)), rid), u_trans_n((*(s->filter)), rid), &b->olist, &(b->sp)); + push_ul_ov_t(&(s->udb), u_trans_a((*(s->filter)), rid), u_trans_n((*(s->filter)), rid), rid, &(b->sp), &b->olist, len, pass_aln, err_high, bl); + // fprintf(stderr, "-1-[M::%s] utg%.6lu%c, rid::%lu, pass_aln::%lu\n", + // __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, pass_aln); + } else {///second round + push_ul_ov_t(&(s->udb), NULL, 0, rid, &(b->sp), &b->olist, len, 0, err_high, bl); + remove_trans_ovlp_connect(s->udb.ug, rid, bl); + } + return pass_aln; +} +**/ + + + +static void worker_for_trans_chain_mmhap_adv(void *data, long i, int tid) // callback for kt_for() +{ + ug_trans_t *s = (ug_trans_t*)data; + ha_ovec_buf_t *b = s->hab[tid]; kv_ul_ov_t *bl = &(s->ll[tid].tk); + uint32_t high_occ = asm_opt.polyploidy + 1; uint64_t cnt; + char *seq = s->ug->u.a[i].s; int64_t len = s->ug->u.a[i].len; + if((!s->is_ovlp) && (s->is_cnt)) s->idx_n.a[i] = 0; + if(s->ug->g->seq[i].del) return; + if(is_mmhom_node(s->ccov->cov.a+s->ccov->idx[i], &(s->ug->u.a[i]), s->ccov->rg, s->ccov->hom_min, 0.9)) return; + // asprintf(&as, "\n[M::%s] rid::%ld, len::%lu, name::%.*s\n", __func__, s->id+i, s->len[i], (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a); + // push_vlog(&(overall_zdbg->a[s->id+i]), as); free(as); as = NULL; + // if(rescue_adject_ovlp(s->ug->g, i, s->filter, &(s->ll[tid].lo))) return; + + // gen_trans_chain_mmhap(s, i, b, bl, seq, len, 0.8, 0.8); + + if(!s->is_ovlp) { + if(s->is_cnt) { + s->idx_n.a[i] = ug_map_lchain(b->abl, i, seq, len, s->w, s->k, &(s->udb), NULL, NULL, s->bw_thres, s->bw_thres_double, + s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3, s->is_HPC, NULL, 0, NULL, 0, s->mini_cut, s->chain_cut, NULL); + } else { + cnt = ug_map_lchain(b->abl, i, seq, len, s->w, s->k, &(s->udb), NULL, NULL, s->bw_thres, s->bw_thres_double, + s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3, s->is_HPC, s->idx_a.a + s->idx_n.a[i], 0, NULL, 0, s->mini_cut, s->chain_cut, NULL); + assert(cnt == ((s->idx_n.a[i+1]-s->idx_n.a[i]))); + } + if(s->free_cnt[tid] >= FREE_BATCH) { + clear_count_buf(s, tid, 1); s->free_cnt[tid] = 0; + } + s->free_cnt[tid]++; + return; + } + + // if(i == 58) { + // fprintf(stderr, "\n-1-[M::%s] utg%.6u%c, rid::%ld, is_ovlp::%d, is_cnt::%d, len::%ld, str::%u\n", + // __func__, (uint32_t)i+1, "lc"[s->ug->u.a[i].circ], i, s->is_ovlp, s->is_cnt, len, (uint32_t)(!!seq)); + // } + + if(!gen_trans_adaptive_mmhap_aln(s, i, b, bl, seq, len, s->filter, s->diff_ec_ul, s->diff_ec_ul_double, s->bw_thres, s->bw_thres_double)) { + gen_trans_adaptive_mmhap_aln(s, i, b, bl, seq, len, NULL, s->diff_ec_ul_double, s->diff_ec_ul_double, s->bw_thres_double, s->bw_thres_double); + } + if(s->free_cnt[tid] >= FREE_BATCH) { + clear_count_buf(s, tid, 0); s->free_cnt[tid] = 0; + } + s->free_cnt[tid]++; +} + int64_t retrieve_cigar_err_dir(bit_extz_t *ez, int64_t s, int64_t e, int64_t *xk, int64_t *ck, int64_t is_back) { ///[ez->ts, ez->te]/[ez->qs, ez->qe]/[s, e) @@ -20178,6 +20374,145 @@ void gen_trans_base_count_comp(ug_trans_t *p, kv_u_trans_t *res) fprintf(stderr, "[M::%s::%.3f] ==> Qualification\n", __func__, yak_realtime()-index_time); } +void clean_trans_base_count_mmhap_comp_rmap(ug_trans_t *p, kv_u_trans_t *res) +{ + uint64_t i, k, l, occ, idx_n; ha_mzl_t *tz; u_trans_t *z; + kv_ul_ov_t *bl; double ww; ha_mzl_t *idx; + ///make results consistent + kv_resize(ha_mzl_t, p->srt_a, p->srt_a.n+p->ug->u.n); + idx = p->srt_a.a + p->srt_a.n; idx_n = p->ug->u.n; + for (i = 0; i < idx_n; i++) { + tz = &(idx[i]); + tz->x = (uint64_t)-1; tz->rev = 0; + tz->pos = tz->rid = tz->span = 0; + } + + for (i = 0, occ = res->n; (int64_t)i < p->n_thread; i++) { + bl = &(p->ll[i].tk); + if(!(bl->n)) continue; + for (k = 1, l = 0; k <= bl->n; k++) { + if(k == bl->n || bl->a[k].qn != bl->a[l].qn) { + if(k > l) { + tz = &(idx[bl->a[l].qn]); + tz->x = bl->a[l].qn; tz->x <<= 32; tz->x |= i; + tz->rid = l>>32; tz->pos = (uint32_t)l; tz->rev = 1; + occ += (k - l); + } + l = k; + } + } + } + + kv_resize(u_trans_t, *res, occ); + for (i = 0; i < idx_n; i++) { + tz = &(idx[i]); + if(!(tz->rev)) continue; + bl = &(p->ll[(uint32_t)(tz->x)].tk); + k = tz->rid; k <<= 32; k += tz->pos; + assert(bl->a[k].qn == (tz->x>>32)); + for (; (k < bl->n) && (bl->a[k].qn == (tz->x>>32)); k++) { + if(bl->a[k].qn == bl->a[k].tn) continue; + ww = cal_trans_ov_w(&(bl->a[k])); + if(ww <= 0) continue; + + kv_pushp(u_trans_t, *res, &z); + z->f = RC_3; z->rev = bl->a[k].rev; z->del = 0; + z->qn = bl->a[k].qn; z->qs = bl->a[k].qs; z->qe = bl->a[k].qe; + z->tn = bl->a[k].tn; z->ts = bl->a[k].ts; z->te = bl->a[k].te; + z->nw = ww; + } + } + destory_ug_rid_cov_t(p->ccov); free(p->ccov); + p->ccov = gen_ug_rid_cov_t(p->ug, p->rg, p->uopt->sources); + + clean_u_trans_t_idx_filter_mmhap_adv(res, p->ug, p->rg, p->uopt->sources, p->ccov); + gen_ug_rid_cov_t_by_ovlp(res, p->ccov); +} + +void gen_trans_base_count_mmhap_comp_rmap(ug_trans_t *p, kv_u_trans_t *res) +{ + double index_time = yak_realtime(); + uint64_t i, k, l, occ, m, cc; + p->ccov = gen_ug_rid_cov_t(p->ug, p->rg, p->uopt->sources); + clean_u_trans_t_idx_adv(res, p->ug, p->rg); p->filter = res; + + p->is_cnt = 1; p->is_ovlp = 0; + memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread); + kt_for(p->n_thread, worker_for_trans_ovlp_mmhap_adv, p, p->ug->u.n); + for (i = l = 0; i < p->ug->u.n; i++) { + occ = p->idx_n.a[i]; p->idx_n.a[i] = l; l += occ; + } + + p->idx_n.a[i] = l; + p->idx_a.n = p->idx_a.m = l; MALLOC(p->idx_a.a, p->idx_a.n); + p->is_cnt = 0; p->is_ovlp = 0; + memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread); + kt_for(p->n_thread, worker_for_trans_ovlp_mmhap_adv, p, p->ug->u.n); + p->srt_a.n = p->srt_a.m = p->idx_a.n; MALLOC(p->srt_a.a, p->srt_a.n); + + for (i = 0; i < p->srt_a.n; i++) { + p->srt_a.a[i] = p->idx_a.a[i]; + p->srt_a.a[i].pos = (uint32_t)i; + p->srt_a.a[i].rid = i>>32; + } + radix_sort_ha_mzl_t_srt(p->srt_a.a, p->srt_a.a + p->srt_a.n); + kvec_t(uint64_t) cut; kv_init(cut); + for (k = 1, l = 0; k <= p->srt_a.n; k++) { + if(k == p->srt_a.n || p->srt_a.a[l].x != p->srt_a.a[k].x) { + for (i = l; i < k; i++) { + m = p->srt_a.a[i].rid; m <<= 32; m |= p->srt_a.a[i].pos; + assert(p->srt_a.a[i].x == p->idx_a.a[m].x); + p->srt_a.a[i] = p->idx_a.a[m]; p->idx_a.a[m].x = i; + } + kv_push(uint64_t, cut, (k - l)); + l = k; + } + } + + if(cut.n > 0) { + radix_sort_gfa64(cut.a, cut.a + cut.n); + m = cut.n * 0.0002; cc = cut.a[cut.n-1] + 1; + if(m > 0 && m <= cut.n) cc = cut.a[cut.n-m] + 1; + if(cc < (uint64_t)p->mini_cut) p->mini_cut = cc; + } + kv_destroy(cut); + + p->is_cnt = 0; p->is_ovlp = 1; + memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread); + kt_for(p->n_thread, worker_for_trans_ovlp_mmhap_adv, p, p->ug->u.n); + + clean_trans_base_count_mmhap_comp_rmap(p, res); + + for (i = 0; (int64_t)i < p->n_thread; i++) p->ll[i].tk.n = 0; + + p->is_cnt = 0; p->is_ovlp = 1; + memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread); + kt_for(p->n_thread, worker_for_trans_chain_mmhap_adv, p, p->ug->u.n); + + + + + + for (i = 0; (int64_t)i < p->n_thread; i++) { + ha_ovec_destroy(p->hab[i]); + free(p->ll[i].lo.a); free(p->ll[i].srt.a.a); free(p->ll[i].tc.a); + } + free(p->idx_a.a); free(p->idx_n.a); free(p->hab); free(p->free_cnt); + destory_ug_rid_cov_t(p->ccov); free(p->ccov); + + + + + + + + + + + for (i = 0; (int64_t)i < p->n_thread; i++) free(p->ll[i].tk.a); + free(p->srt_a.a); free(p->ll); + fprintf(stderr, "[M::%s::%.3f] ==> Qualification\n", __func__, yak_realtime()-index_time); +} void gen_trans_base_count_mmhap_comp(ug_trans_t *p, kv_u_trans_t *res) {