From 8aa87fdce8380991bca534f438f5f82218a13177 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Mon, 8 Mar 2021 20:52:11 -0500 Subject: [PATCH] purge_dups for high het --- CommandLines.cpp | 6 +- CommandLines.h | 2 +- Overlaps.cpp | 724 +++++-------- Overlaps.h | 43 +- Purge_Dups.cpp | 2534 +++++++++++++++++++++++++++------------------- Purge_Dups.h | 8 +- hic.cpp | 1211 +++++++++++++++++----- hifiasm.1 | 3 +- 8 files changed, 2766 insertions(+), 1765 deletions(-) diff --git a/CommandLines.cpp b/CommandLines.cpp index 4685b0d..caea7e4 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -89,7 +89,7 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " -4 FILE list of hap2/maternal read names []\n"); fprintf(stderr, " Purge-dups:\n"); - fprintf(stderr, " -l INT purge level. 0: no purging; 1: light; 2: aggressive [0 for trio; 2 for unzip]\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, " -O INT min number of overlapped reads for duplicate haplotigs [%d]\n", @@ -360,9 +360,9 @@ int check_option(hifiasm_opt_t* asm_opt) return 0; } - if(asm_opt->purge_level_primary < 0 || asm_opt->purge_level_primary > 2) + if(asm_opt->purge_level_primary < 0 || asm_opt->purge_level_primary > 3) { - fprintf(stderr, "[ERROR] the level of purge-dup should be [0, 2] (-l)\n"); + fprintf(stderr, "[ERROR] the level of purge-dup should be [0, 3] (-l)\n"); return 0; } diff --git a/CommandLines.h b/CommandLines.h index 309fd8d..253e6a8 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.14-r312" +#define HA_VERSION "0.14-r313" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index 3e320d8..6f3ba54 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -11062,7 +11062,6 @@ inline uint64_t get_utg_len(buf_t* b, ma_ug_t *ug, asg_t *read_sg, uint64_t igno for (r_i = 0; r_i < u->n; r_i++, idx++) { v = ((uint64_t)((u->a[u->n - r_i - 1])^(uint64_t)(0x100000000)))>>32; - ///w = ((uint64_t)((u->a[u->n - x->r_i - 2])^(uint64_t)(0x100000000)))>>32; if(p_v == (uint32_t)-1) { p_v = v; @@ -11085,7 +11084,7 @@ inline uint64_t get_utg_len(buf_t* b, ma_ug_t *ug, asg_t *read_sg, uint64_t igno len += l; if(len_thre && occ && len >= (*len_thre)) { - (*occ) = idx - 1; + (*occ) = idx; return len; } p_v = v; @@ -11120,7 +11119,7 @@ inline uint64_t get_utg_len(buf_t* b, ma_ug_t *ug, asg_t *read_sg, uint64_t igno len += l; if(len_thre && occ && len >= (*len_thre)) { - (*occ) = idx - 1; + (*occ) = idx; return len; } p_v = v; @@ -11133,7 +11132,7 @@ inline uint64_t get_utg_len(buf_t* b, ma_ug_t *ug, asg_t *read_sg, uint64_t igno len += read_sg->seq[p_v>>1].len; if(len_thre && occ && len >= (*len_thre)) { - (*occ) = idx - 1; + (*occ) = idx; return len; } } @@ -11208,22 +11207,6 @@ void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug 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); - - // if(pre_0 == 5 || pre_1 == 5) - // { - // fprintf(stderr, "\npre_0: utg%.6ul, len_0: %lu, thre_0: %lu, pre_1: utg%.6ul, len_1: %lu, thre_1: %lu\n", - // (int)(pre_0+1), len_0, thre_0, (int)(pre_1+1), len_1, thre_1); - // uint32_t xxx_i; - // for (xxx_i = 0; xxx_i < b_0->b.n; xxx_i++) - // { - // fprintf(stderr,"+: utg%.6ul\n", (int)((b_0->b.a[xxx_i]>>1)+1)); - // } - - // for (xxx_i = 0; xxx_i < b_1->b.n; xxx_i++) - // { - // fprintf(stderr,"-: utg%.6ul\n", (int)((b_1->b.a[xxx_i]>>1)+1)); - // } - // } } } @@ -11233,9 +11216,119 @@ void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug } +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; + uint32_t *a = b->b.a, a_n = b->b.n; + uint32_t u_i, r_i, len, p_v; + asg_arc_t *av = NULL; + ma_utg_t* u = NULL; + for (u_i = r_i = len = 0, p_v = (uint32_t)-1; u_i < a_n; u_i++) + { + uid = a[u_i] >> 1; + ori = a[u_i] & 1; + u = &(ug->u.a[uid]); + if(u->n == 0) continue; + + for (r_i = 0; r_i < u->n; r_i++) + { + l = 0; + v = (ori == 1?((uint64_t)((u->a[u->n - r_i - 1])^(uint64_t)(0x100000000)))>>32:((uint64_t)(u->a[r_i]))>>32); + + if(p_v != (uint32_t)-1 && is_clear == 0) + { + av = asg_arc_a(read_sg, p_v); + nv = asg_arc_n(read_sg, p_v); + + for (k = 0; k < nv; k++) + { + if(av[k].del) continue; + if(av[k].v == v) + { + l = asg_arc_len(av[k]); + break; + } + } + if(k == nv) fprintf(stderr, "ERROR\n"); + } + + p_v = v; len += l; + + if(is_clear == 1) + { + cov->pos_idx[v>>1] = (uint64_t)-1; + } + else + { + cov->pos_idx[v>>1] = len; + cov->pos_idx[v>>1] <<= 32; + cov->pos_idx[v>>1] |= (uint64_t)v; + } + } + } + + if(p_v != (uint32_t)-1) len += read_sg->seq[p_v>>1].len; + + return len; +} + + +void collect_trans_cov(buf_t* pri, buf_t* aux, ma_ug_t *ug, asg_t *read_sg, hap_cov_t *cov) +{ + uint32_t i, k, rid, occ, thre_pri; + uint64_t len_aux, uLen, uCov; + ma_utg_t* u = NULL; + if(pri->b.n == 0 || aux->b.n == 0) return; + + len_aux = set_utg_offset(aux, ug, read_sg, cov, 0); + chain_trans_ovlp(cov, ug, read_sg, pri, len_aux, &thre_pri); + if(thre_pri > 0) + { + for (i = uCov = 0; i < aux->b.n; i++) + { + u = &(ug->u.a[aux->b.a[i]>>1]); + if(u->n == 0) continue; + for (k = 0; k < u->n; k++) + { + rid = u->a[k]>>33; + uCov += cov->cov[rid]; + } + } + + + for (i = uLen = occ = 0; i < pri->b.n; i++) + { + u = &(ug->u.a[pri->b.a[i]>>1]); + if(u->n == 0) continue; + for (k = 0; k < u->n; k++, occ++) + { + if(occ >= thre_pri) break; + rid = u->a[k]>>33; + uLen += read_sg->seq[rid].len; + } + if(occ >= thre_pri) break; + } + + uCov = (uLen == 0? 0 : uCov / uLen); + + for (i = occ = 0; i < pri->b.n; i++) + { + u = &(ug->u.a[pri->b.a[i]>>1]); + if(u->n == 0) continue; + for (k = 0; k < u->n; k++, occ++) + { + if(occ >= thre_pri) break; + rid = u->a[k]>>33; + cov->cov[rid] += (uCov * read_sg->seq[rid].len); + } + if(occ >= thre_pri) break; + } + } + set_utg_offset(aux, ug, read_sg, cov, 1); +} int untig_asg_arc_simple_large_bubbles_trio(ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, -long long miniedgeLen, R_to_U* ruIndex, uint32_t positive_flag, uint32_t negative_flag, hc_links* link) +long long miniedgeLen, R_to_U* ruIndex, uint32_t positive_flag, uint32_t negative_flag, hc_links* link, hap_cov_t *cov) { asg_t *g = ug->g; double startTime = Get_T(); @@ -11360,6 +11453,7 @@ long long miniedgeLen, R_to_U* ruIndex, uint32_t positive_flag, uint32_t negativ } if(link) collect_reverse_unitigs(&b_0, &b_1, link, ug, read_sg); + if(cov) collect_trans_cov(&b_0, &b_1, ug, read_sg, cov); is_hap++; } @@ -11551,11 +11645,13 @@ kvec_asg_arc_t_warp* new_rtg_edges, int max_hang, int min_ovlp) { uint64_t i, dip_thre_max, dip_thres, n_utg; uint8_t* primary_flag = (uint8_t*)calloc(sg->n_seq, sizeof(uint8_t)); + hap_cov_t *cov = init_hap_cov_t(ug, sg, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp); + int tmp_cov = asm_opt.hom_global_coverage; 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, 0, 1, NULL); + asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, 0, 0, 0, 1, NULL, 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); @@ -11579,7 +11675,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, int max_hang, int min_ovlp) } } free(primary_flag); - + destory_hap_cov_t(&cov); ///fprintf(stderr, "[M::%s] diploid coverage threshold: %lu\n", __func__, dip_thres); } @@ -11622,92 +11718,6 @@ void destory_hc_links(hc_links* link) kv_destroy(link->enzymes); } -void pop_small_bub(ma_ug_t *ug) -{ - bubble_type bub; - uint64_t n_vtx = ug->g->n_seq*2, tLen, pathLen, nodeLen; - uint32_t i, k, v, mode = (((uint32_t)-1)<<2); - asg_cleanup(ug->g); if (!ug->g->is_symm) asg_symm(ug->g); - memset(&bub, 0, sizeof(bubble_type)); - CALLOC(bub.index, n_vtx); - kv_init(bub.list); kv_init(bub.num); kv_init(bub.pathLen); - buf_t b; memset(&b, 0, sizeof(buf_t)); b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); - for (i = 0, tLen = 1; i < ug->u.n; i++) tLen += ug->u.a[i].len; - for (v = 0; v < n_vtx; ++v) - { - if(ug->g->seq[v>>1].del) continue; - if(asg_arc_n(ug->g, v) < 2) continue; - if((bub.index[v]&(uint32_t)3) != 0) continue; - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL)) - { - //beg is v, end is b.S.a[0] - //note b.b include end, does not include beg - for (i = 0; i < b.b.n; i++) - { - if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue; - bub.index[b.b.a[i]] &= mode; bub.index[b.b.a[i]] += 1; - bub.index[b.b.a[i]^1] &= mode; bub.index[b.b.a[i]^1] += 1; - } - bub.index[v] &= mode; bub.index[v] += 2; - bub.index[b.S.a[0]^1] &= mode; bub.index[b.S.a[0]^1] += 3; - } - } - - for (v = 0; v < n_vtx; ++v) - { - if((bub.index[v]&(uint32_t)3) !=2) continue; - kv_push(uint32_t, bub.num, bub.list.n); - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL)) - { - kv_push(uint64_t, bub.pathLen, pathLen); - //beg is v, end is b.S.a[0] - kv_push(uint32_t, bub.list, v); - kv_push(uint32_t, bub.list, b.S.a[0]^1); - - //note b.b include end, does not include beg - for (i = 0; i < b.b.n; i++) - { - if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue; - kv_push(uint32_t, bub.list, b.b.a[i]); - } - } - } - - kv_push(uint32_t, bub.num, bub.list.n); - ///free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); - bub.f_bub = bub.num.n - 1; bub.b_bub = 0; - memset(bub.index, 0, n_vtx*sizeof(uint32_t)); - - uint32_t beg, sink, n, *a, total_nodes; - for (i = 0; i < bub.f_bub; i++) - { - get_bubbles(&bub, i, &beg, &sink, &a, &n, &pathLen); - for (k = total_nodes = 0; k < n; k++) - { - total_nodes += ug->u.a[a[k]>>1].n; - } - - for (k = 0; k < n; k++) - { - v = a[k]; - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, &nodeLen)) - { - - } - - - v = a[k]^1; - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, &nodeLen)) - { - - } - - } - - } - -} - void hic_clean(asg_t* read_g) { uint32_t n_vtx, v, u; @@ -11734,7 +11744,7 @@ void hic_clean(asg_t* read_g) if(ug->g->seq[v>>1].del) continue; if(asg_arc_n(ug->g, v) < 2) continue; if(bs_flag[v] != 0) continue; - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL)) + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL)) { //beg is v, end is b.S.a[0] //note b.b include end, does not include beg @@ -11751,7 +11761,7 @@ void hic_clean(asg_t* read_g) for (v = 0; v < n_vtx; ++v) { if(bs_flag[v] !=2) continue; - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL)) + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL)) { //note b.b include end, does not include beg for (i = v_occ = ax.n = 0; i < b.b.n; i++) @@ -11767,7 +11777,7 @@ void hic_clean(asg_t* read_g) { u = (ax.a[i]<<1) + k; if(asg_arc_n(ug->g, u) < 2) continue; - if(asg_bub_pop1_primary_trio(ug->g, NULL, u, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL)) + if(asg_bub_pop1_primary_trio(ug->g, NULL, u, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL)) { for (k_i = u_occ = utg_occ = 0; k_i < b.b.n; k_i++) { @@ -11779,7 +11789,7 @@ void hic_clean(asg_t* read_g) if(u_occ >= v_occ*bub_rate) continue; if(u_occ > 3) continue; if(utg_occ > 2) continue; - asg_bub_pop1_primary_trio(ug->g, NULL, u, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL); + asg_bub_pop1_primary_trio(ug->g, NULL, u, tLen, &b, (uint32_t)-1, (uint32_t)-1, 1, NULL, NULL, NULL); } } } @@ -11824,7 +11834,7 @@ R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ov hc_links link; - + if(load_hc_links(&link, output_file_name) == 0) { init_hc_links(&link, ug->g->n_seq, R_INF.total_reads); @@ -12513,7 +12523,7 @@ asg_t *read_sg, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint32_t min_e int asg_arc_cut_trio_long_tip_primary(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, -R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hc_links* link) +R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hc_links* link, hap_cov_t *cov) { double startTime = Get_T(); ///the reason is that each read has two direction (query->target, target->query) @@ -12613,6 +12623,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hc_links* link) } if(link && operation != CUT) collect_reverse_unitigs(&b_0, &b_1, link, ug, read_sg); + if(cov && operation != CUT) collect_trans_cov(&b_0, &b_1, ug, read_sg, cov); } } } @@ -12635,7 +12646,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, hc_links* link) } int asg_arc_cut_trio_long_tip_primary_complex(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, -R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_threshold, hc_links* link) +R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_threshold, hc_links* link, hap_cov_t *cov) { double startTime = Get_T(); uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, in, flag, operation; @@ -12711,6 +12722,7 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_thre } if(link && operation != CUT) collect_reverse_unitigs(&b_0, &b_1, link, ug, read_sg); + if(cov && operation != CUT) collect_trans_cov(&b_0, &b_1, ug, read_sg, cov); break; } @@ -12810,7 +12822,7 @@ long long* base_maxLen, long long* base_maxLen_i, uint32_t stops_threshold, buf_ int asg_arc_cut_trio_long_equal_tips_assembly(asg_t *g, ma_ug_t *ug, asg_t *read_sg, ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_t trio_flag, -hc_links* link) +hc_links* link, hap_cov_t *cov) { double startTime = Get_T(); uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, flag, is_hap, n_tips, return_flag, k; @@ -12898,6 +12910,7 @@ hc_links* link) } if(link) collect_reverse_unitigs(&b_0, &b_1, link, ug, read_sg); + if(cov) collect_trans_cov(&b_0, &b_1, ug, read_sg, cov); is_hap++; } @@ -13130,7 +13143,7 @@ R_to_U* ruIndex, uint32_t positive_flag, float drop_rate) } int asg_arc_cut_trio_long_equal_tips_assembly_complex(asg_t *g, ma_ug_t *ug, asg_t *read_sg, -ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_t stops_threshold, hc_links* link) +ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_t stops_threshold, hc_links* link, hap_cov_t *cov) { double startTime = Get_T(); uint32_t v, n_vtx = g->n_seq * 2, n_reduced = 0, convex, in, flag; @@ -13199,6 +13212,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_ } if(link) collect_reverse_unitigs(&b_0, &b_1, link, ug, read_sg); + if(cov) collect_trans_cov(&b_0, &b_1, ug, read_sg, cov); ///lable the primary one b_0.b.n = 0; @@ -13600,8 +13614,8 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate) redo: ///print_untig((ug), 61955, "i-0:", 0); - asg_pop_bubble_primary_trio(ug, bubble_dist, trio_flag, DROP); - untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, trio_flag, DROP, NULL); + asg_pop_bubble_primary_trio(ug, bubble_dist, trio_flag, DROP, NULL); + untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, trio_flag, DROP, NULL, NULL); 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**********/ @@ -13617,16 +13631,15 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate) { pre_cons = get_graph_statistic(g); ///need consider tangles - ///asg_pop_bubble_primary(g, bubble_dist); - asg_pop_bubble_primary_trio(ug, bubble_dist, trio_flag, DROP); + asg_pop_bubble_primary_trio(ug, bubble_dist, trio_flag, DROP, NULL); /**********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, NULL); - asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, NULL); - asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, NULL); - asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, NULL); + asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, NULL, NULL); + asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, trio_flag, NULL, NULL); + asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, NULL, NULL); + asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, NULL, NULL); ///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); @@ -13636,7 +13649,7 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate) /**********debug**********/ cur_cons = get_graph_statistic(g); } - untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, trio_flag, DROP, NULL); + untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, trio_flag, DROP, NULL, NULL); if(just_bubble_pop == 0) { @@ -13661,7 +13674,7 @@ void clean_primary_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* rever long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen, uint32_t miniBiGraph, float chimeric_rate, int is_final_clean, int just_bubble_pop, -float drop_ratio, hc_links* link) +float drop_ratio, hc_links* link, hap_cov_t *cov) { #define T_ROUND 2 asg_t *g = ug->g; @@ -13669,9 +13682,8 @@ float drop_ratio, hc_links* link) redo: - asg_pop_bubble_primary_trio(ug, bubble_dist, (uint32_t)-1, DROP); - untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, DROP, link); - + asg_pop_bubble_primary_trio(ug, bubble_dist, (uint32_t)-1, DROP, cov); + untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, DROP, link, cov); if(just_bubble_pop == 0) { cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, @@ -13683,17 +13695,16 @@ float drop_ratio, hc_links* link) while(pre_cons != cur_cons) { pre_cons = get_graph_statistic(g); - asg_pop_bubble_primary_trio(ug, bubble_dist, (uint32_t)-1, DROP); + asg_pop_bubble_primary_trio(ug, bubble_dist, (uint32_t)-1, DROP, cov); 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, link); - asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, link); - asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, link); - asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, link); + asg_arc_cut_trio_long_tip_primary(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, link, cov); + asg_arc_cut_trio_long_equal_tips_assembly(g, ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, link, cov); + asg_arc_cut_trio_long_tip_primary_complex(g, ug, read_g, reverse_sources, ruIndex, 2, tip_drop_ratio, stops_threshold, link, cov); + asg_arc_cut_trio_long_equal_tips_assembly_complex(g, ug, read_g, reverse_sources, 2, ruIndex, stops_threshold, link, cov); detect_chimeric_by_topo(g, ug, read_g, reverse_sources, 2, stops_threshold, chimeric_rate, ruIndex); - if(round != T_ROUND) { unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, @@ -13702,14 +13713,12 @@ float drop_ratio, hc_links* link) } cur_cons = get_graph_statistic(g); } - - untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, DROP, link); + untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, DROP, link, cov); if(just_bubble_pop == 0) { cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 2); } - resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, (uint32_t)-1, drop_ratio); drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex); unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, reverse_sources, 0, 1); @@ -14366,14 +14375,15 @@ kvec_asg_arc_t_warp* new_rtg_edges) { asg_t* nsg = (*ug)->g; uint32_t v, n_vtx = nsg->n_seq; - + hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp); + 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, bubble_dist, - drop_ratio, 1, 1, NULL); + drop_ratio, 1, 1, NULL, cov); if(asm_opt.recover_atg_cov_min == -1024) { asm_opt.recover_atg_cov_max = asm_opt.hom_global_coverage/HOM_PEAK_RATE; - asm_opt.recover_atg_cov_min = asm_opt.recover_atg_cov_max * 0.8; + asm_opt.recover_atg_cov_min = asm_opt.recover_atg_cov_max * 0.85; ///asm_opt.recover_atg_cov_max = asm_opt.recover_atg_cov_max * 1.2; asm_opt.recover_atg_cov_max = INT32_MAX; } @@ -14450,13 +14460,14 @@ kvec_asg_arc_t_warp* new_rtg_edges) { 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, bubble_dist, - drop_ratio, 1, 0, NULL); + drop_ratio, 1, 0, NULL, cov); ///delete_useless_nodes(ug); delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex); } set_drop_trio_flag(*ug); + destory_hap_cov_t(&cov); } @@ -14726,10 +14737,73 @@ static inline int count_out(const asg_t *g, uint32_t v) } +// in a resolved bubble, mark unused vertices and arcs as "reduced" +static void asg_bub_backtrack_primary_cov(ma_ug_t *ug, uint32_t v0, buf_t *b, hap_cov_t *cov) +{ + uint32_t i, k, v, u, uLen = 0, uCov = 0, uId, rId; + ma_utg_t* p = NULL; + ///b->S.a[0] is the sink of this bubble + + ///assert(b->S.n == 1); + ///first remove all nodes in this bubble + for (i = 0; i < b->b.n; ++i) + { + uId = b->b.a[i]>>1; + if(uId == (b->S.a[0]>>1)) continue; + p = &(ug->u.a[uId]); + if(p->n == 0) continue; + for (k = 0; k < p->n; k++) + { + rId = p->a[k]>>33; + uCov += cov->cov[rId]; + } + } + + ///v is the sink of this bubble + v = b->S.a[0]; + ///recover node + do { + u = b->a[v].p; // u->v + if(v != b->S.a[0]) + { + uId = v>>1; + p = &(ug->u.a[uId]); + if(p->n == 0) continue; + for (k = 0; k < p->n; k++) + { + rId = p->a[k]>>33; + uCov -= cov->cov[rId]; + uLen += cov->read_g->seq[rId].len; + } + } + v = u; + } while (v != v0); + + uCov = (uLen == 0? 0 : uCov / uLen); + + ///v is the sink of this bubble + v = b->S.a[0]; + ///recover node + do { + u = b->a[v].p; // u->v + if(v != b->S.a[0]) + { + uId = v>>1; + p = &(ug->u.a[uId]); + if(p->n == 0) continue; + for (k = 0; k < p->n; k++) + { + rId = p->a[k]>>33; + cov->cov[rId] += (uCov * cov->read_g->seq[rId].len); + } + } + v = u; + } while (v != v0); +} // in a resolved bubble, mark unused vertices and arcs as "reduced" -static void asg_bub_backtrack_primary(asg_t *g, uint32_t v0, buf_t *b) +void asg_bub_backtrack_primary(asg_t *g, uint32_t v0, buf_t *b) { uint32_t i, v, qn, tn; ///b->S.a[0] is the sink of this bubble @@ -15108,7 +15182,8 @@ pop_reset: uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, int max_dist, buf_t *b, -uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes) +uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes, +hap_cov_t *cov) { uint32_t i, n_pending = 0, is_first = 1, cur_m, cur_c, cur_np, cur_nc, to_replace, n_tips, tip_end; uint64_t n_pop = 0; @@ -15330,9 +15405,9 @@ uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop, uint64_t* path_ if (i < nv || b->S.n == 0) goto pop_reset; } while (b->S.n > 1 || n_pending); + if(cov && utg) asg_bub_backtrack_primary_cov(utg, v0, b, cov); if(is_pop) asg_bub_backtrack_primary(g, v0, b); if(path_base_len || path_nodes) asg_bub_backtrack_primary_length(g, utg, v0, b, path_base_len, path_nodes); - n_pop = 1; pop_reset: @@ -15345,7 +15420,7 @@ pop_reset: // pop bubbles -int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_flag, uint32_t negative_flag) +int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov) { asg_t *g = ug->g; uint32_t v, n_vtx = g->n_seq * 2; @@ -15365,7 +15440,7 @@ int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_fla for (i = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs if (!av[i].del) ++n_arc; if (n_arc > 1) - n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL); + n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL, cov); } free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); if (n_pop) asg_cleanup(g); @@ -15378,303 +15453,6 @@ int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_fla } - -// pop bubbles from vertex v0; the graph MJUST BE symmetric: if u->v present, v'->u' must be present as well -uint64_t asg_bub_pop1_primary_trio_debug(ma_ug_t *ug, uint32_t v0, int max_dist, buf_t *b, -uint32_t positive_flag, uint32_t negative_flag) -{ - asg_t *g = ug->g; - uint32_t i, n_pending = 0, is_first = 1, cur_m, cur_c, to_replace, n_tips, tip_end; - uint64_t n_pop = 0; - ///if this node has been deleted - if (g->seq[v0>>1].del || g->seq[v0>>1].c == ALTER_LABLE) return 0; // already deleted - ///asg_arc_n(n0) - if ((uint32_t)g->idx[v0] < 2) return 0; // no bubbles - ///S saves nodes with all incoming edges visited - b->S.n = b->T.n = b->b.n = b->e.n = 0; - ///for each node, b->a saves all related information - b->a[v0].c = b->a[v0].d = b->a[v0].m = 0; - ///b->S is the nodes with all incoming edges visited - kv_push(uint32_t, b->S, v0); - n_tips = 0; - tip_end = (uint32_t)-1; - - do { - ///v is a node that all incoming edges have been visited - ///d is the distance from v0 to v - uint32_t v = kv_pop(b->S), d = b->a[v].d, c = b->a[v].c, m = b->a[v].m; - uint32_t nv = asg_arc_n(g, v); - asg_arc_t *av = asg_arc_a(g, v); - ///why we have this assert? - ///assert(nv > 0); - ///all out-edges of v - for (i = 0; i < nv; ++i) { // loop through v's neighbors - /** - p->ul: |____________31__________|__________1___________|______________32_____________| - qn direction of overlap length of this node (not overlap length) - (in the view of query) - p->v : |___________31___________|__________1___________| - tn reverse direction of overlap - (in the view of target) - p->ol: overlap length - **/ - - uint32_t w = av[i].v, l = (uint32_t)av[i].ul; // v->w with length l - binfo_t *t = &b->a[w]; - ///that means there is a circle, directly terminate the whole bubble poping - ///if (w == v0) goto pop_reset; - if ((w>>1) == (v0>>1)) goto pop_reset; - /****************************may have bugs********************************/ - ///important when poping at long untig graph - if(is_first) l = 0; - /****************************may have bugs********************************/ - - ///if this edge has been deleted - if (av[i].del) continue; - - ///push the edge - ///high 32-bit of g->idx[v] is the start point of v's edges - //so here is the point of this specfic edge - kv_push(uint32_t, b->e, (g->idx[v]>>32) + i); - ///find a too far path? directly terminate the whole bubble poping - if (d + l > (uint32_t)max_dist) break; // too far - - ///if this node - if (t->s == 0) { // this vertex has never been visited - kv_push(uint32_t, b->b, w); // save it for revert - ///t->p is the parent node of - ///t->s = 1 means w has been visited - ///d is len(v0->v), l is len(v->w), so t->d is len(v0->w) - t->p = v, t->s = 1, t->d = d + l; - /****************************may have bugs********************************/ - cur_c = get_num_trio_flag(ug, w>>1, positive_flag); - cur_m = get_num_trio_flag(ug, w>>1, negative_flag); - t->c = c + cur_c; - t->m = m + cur_m; - /****************************may have bugs********************************/ - ///incoming edges of w - t->r = count_out(g, w^1); - ++n_pending; - } else { // visited before - /****************************may have bugs********************************/ - cur_c = get_num_trio_flag(ug, w>>1, positive_flag); - cur_m = get_num_trio_flag(ug, w>>1, negative_flag); - ///select the way with less negative_flag, more positive_flag, more distance - to_replace = 0; - if(m + cur_m < t->m) to_replace = 1; - if(to_replace == 0 && m + cur_m == t->m && c + cur_c > t->c) to_replace = 1; - if(to_replace == 0 && m + cur_m == t->m && c + cur_c == t->c && d + l > t->d) to_replace = 1; - if(to_replace) - { - t->p = v; - t->m = m + cur_m; - t->c = c + cur_c; - } - ///c is the weight (is very likely the number of node in this edge) of the parent node - ///select the longest edge (longest meams most reads/longest edge) - // if (c + 1 > t->c || (c + 1 == t->c && d + l > t->d)) t->p = v; - // if (c + 1 > t->c) t->c = c + 1; - /****************************may have bugs********************************/ - ///update len(v0->w) - ///node: t->d is not the length from this node's parent - ///it is the shortest edge - if (d + l < t->d) t->d = d + l; // update dist - } - ///assert(t->r > 0); - //if all incoming edges of w have visited - //push it to b->S - if (--(t->r) == 0) { - uint32_t x = asg_arc_n(g, w); - /****************************may have bugs for bubble********************************/ - /** - if (x) kv_push(uint32_t, b->S, w); - ///else kv_push(uint32_t, b->T, w); // a tip - else goto pop_reset; - **/ - if(x > 0) - { - kv_push(uint32_t, b->S, w); - } - else - { - ///at most one tip - if(n_tips != 0) goto pop_reset; - n_tips++; - tip_end = w; - } - /****************************may have bugs for bubble********************************/ - --n_pending; - } - } - is_first = 0; - //if found a tip - /****************************may have bugs for bubble********************************/ - if(n_tips == 1) - { - if(tip_end != (uint32_t)-1 && n_pending == 0 && b->S.n == 0) - { - kv_push(uint32_t, b->S, tip_end); - break; - } - else - { - goto pop_reset; - } - } - /****************************may have bugs for bubble********************************/ - ///if i < nv, that means (d + l > max_dist) - if (i < nv || b->S.n == 0) goto pop_reset; - } while (b->S.n > 1 || n_pending); - - n_pop = 1; -pop_reset: - for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices - binfo_t *t = &b->a[b->b.a[i]]; - t->s = t->c = t->d = t->m = 0; - } - return n_pop; -} - - - -// pop bubbles from vertex v0; the graph MJUST BE symmetric: if u->v present, v'->u' must be present as well -static uint64_t asg_bub_pop1_primary(asg_t *g, uint32_t v0, int max_dist, buf_t *b) -{ - uint32_t i, n_pending = 0, is_first = 1; - uint64_t n_pop = 0; - ///if this node has been deleted - if (g->seq[v0>>1].del || g->seq[v0>>1].c == ALTER_LABLE) return 0; // already deleted - ///asg_arc_n(n0) - if ((uint32_t)g->idx[v0] < 2) return 0; // no bubbles - ///S saves nodes with all incoming edges visited - b->S.n = b->T.n = b->b.n = b->e.n = 0; - ///for each node, b->a saves all related information - b->a[v0].c = b->a[v0].d = 0; - ///b->S is the nodes with all incoming edges visited - kv_push(uint32_t, b->S, v0); - - do { - ///v is a node that all incoming edges have been visited - ///d is the distance from v0 to v - uint32_t v = kv_pop(b->S), d = b->a[v].d, c = b->a[v].c; - uint32_t nv = asg_arc_n(g, v); - asg_arc_t *av = asg_arc_a(g, v); - ///why we have this assert? - ///assert(nv > 0); - ///all out-edges of v - for (i = 0; i < nv; ++i) { // loop through v's neighbors - /** - p->ul: |____________31__________|__________1___________|______________32_____________| - qn direction of overlap length of this node (not overlap length) - (in the view of query) - p->v : |___________31___________|__________1___________| - tn reverse direction of overlap - (in the view of target) - p->ol: overlap length - **/ - - uint32_t w = av[i].v, l = (uint32_t)av[i].ul; // v->w with length l - binfo_t *t = &b->a[w]; - ///that means there is a circle, directly terminate the whole bubble poping - ///if (w == v0) goto pop_reset; - if ((w>>1) == (v0>>1)) goto pop_reset; - /****************************may have bugs********************************/ - ///important when poping at long untig graph - if(is_first) l = 0; - /****************************may have bugs********************************/ - - ///if this edge has been deleted - if (av[i].del) continue; - - ///push the edge - ///high 32-bit of g->idx[v] is the start point of v's edges - //so here is the point of this specfic edge - kv_push(uint32_t, b->e, (g->idx[v]>>32) + i); - ///find a too far path? directly terminate the whole bubble poping - if (d + l > (uint32_t)max_dist) break; // too far - - ///if this node - if (t->s == 0) { // this vertex has never been visited - kv_push(uint32_t, b->b, w); // save it for revert - ///t->p is the parent node of - ///t->s = 1 means w has been visited - ///d is len(v0->v), l is len(v->w), so t->d is len(v0->w) - t->p = v, t->s = 1, t->d = d + l, t->c = c + 1; - ///incoming edges of w - t->r = count_out(g, w^1); - ++n_pending; - } else { // visited before - ///c is the weight (is very likely the number of node in this edge) of the parent node - ///select the longest edge (longest meams most reads/longest edge) - if (c + 1 > t->c || (c + 1 == t->c && d + l > t->d)) t->p = v; - if (c + 1 > t->c) t->c = c + 1; - ///update len(v0->w) - ///node: t->d is not the length from this node's parent - ///it is the shortest edge - if (d + l < t->d) t->d = d + l; // update dist - } - ///assert(t->r > 0); - //if all incoming edges of w have visited - //push it to b->S - if (--(t->r) == 0) { - uint32_t x = asg_arc_n(g, w); - if (x) kv_push(uint32_t, b->S, w); - ///else kv_push(uint32_t, b->T, w); // a tip - else goto pop_reset; - --n_pending; - } - } - is_first = 0; - ///if i < nv, that means (d + l > max_dist) - if (i < nv || b->S.n == 0) goto pop_reset; - } while (b->S.n > 1 || n_pending); - asg_bub_backtrack_primary(g, v0, b); - n_pop = 1; -pop_reset: - for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices - binfo_t *t = &b->a[b->b.a[i]]; - t->s = t->c = t->d = 0; - } - return n_pop; -} - - -// pop bubbles -int asg_pop_bubble_primary(asg_t *g, int max_dist) -{ - uint32_t v, n_vtx = g->n_seq * 2; - uint64_t n_pop = 0; - buf_t b; - if (!g->is_symm) asg_symm(g); - memset(&b, 0, sizeof(buf_t)); - ///set information for each node - b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); - //traverse all node with two directions - for (v = 0; v < n_vtx; ++v) { - uint32_t i, n_arc = 0, nv = asg_arc_n(g, v); - asg_arc_t *av = asg_arc_a(g, v); - ///some node could be deleted - if (nv < 2 || g->seq[v>>1].del || g->seq[v>>1].c == ALTER_LABLE) continue; - ///some edges could be deleted - for (i = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs - if (!av[i].del) ++n_arc; - if (n_arc > 1) - n_pop += asg_bub_pop1_primary(g, v, max_dist, &b); - } - free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); - if (n_pop) asg_cleanup(g); - - if(VERBOSE >= 1) - { - fprintf(stderr, "[M::%s] popped %lu bubbles\n", __func__, (unsigned long)n_pop); - } - return n_pop; -} - - - - - int test_triangular_directly(asg_t *g, uint32_t v, long long min_edge_length, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex) { @@ -18767,13 +18545,13 @@ uint32_t positive_flag, uint32_t negative_flag) v = beg; if((!g->seq[v>>1].del)&&(g->seq[v>>1].c!=ALTER_LABLE)&&get_real_length(g, v, NULL)>=2) { - n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL); + n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL, NULL); } v = end^1; if((!g->seq[v>>1].del)&&(g->seq[v>>1].c!=ALTER_LABLE)&&get_real_length(g, v, NULL)>=2) { - n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL); + n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL, NULL); } @@ -18788,7 +18566,7 @@ uint32_t positive_flag, uint32_t negative_flag) { v = v|k; if(get_real_length(g, v, NULL)<=1) continue; - n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL); + n_pop += asg_bub_pop1_primary_trio(ug->g, ug, v, max_dist, &b, positive_flag, negative_flag, 1, NULL, NULL, NULL); } } @@ -21323,6 +21101,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link) asg_t* nsg = (*ug)->g; uint32_t v, n_vtx = nsg->n_seq, k, rId, just_contain; ma_utg_t* u = NULL; + hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp); ///print_utg_coverage(*ug, coverage_cut, 440, sources); ///exit(0); @@ -21344,12 +21123,10 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link) } } } - 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; n_vtx = nsg->n_seq; for (v = 0; v < n_vtx; ++v) @@ -21358,22 +21135,17 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link) nsg->seq[v].c = PRIMARY_LABLE; EvaluateLen((*ug)->u, v) = (*ug)->u.a[v].n; } - - clean_primary_untig_graph(*ug, read_g, reverse_sources, bubble_dist, tipsLen, - tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, - chimeric_rate, 0, 0, drop_ratio, link); - + clean_primary_untig_graph(*ug, read_g, reverse_sources, bubble_dist, tipsLen, tip_drop_ratio, + stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, chimeric_rate, 0, 0, drop_ratio, link, cov); delete_useless_nodes(ug); renew_utg(ug, read_g, new_rtg_edges); - if(asm_opt.purge_level_primary > 0) { 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, bubble_dist, - drop_ratio, just_contain, 0, link); + drop_ratio, just_contain, 0, link, cov); delete_useless_nodes(ug); renew_utg(ug, read_g, new_rtg_edges); } @@ -21391,11 +21163,9 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link) { 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, bubble_dist, - drop_ratio, just_contain, 0, link); - + drop_ratio, just_contain, 0, link, cov); delete_useless_nodes(ug); renew_utg(ug, read_g, new_rtg_edges); } @@ -21405,7 +21175,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link) { 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, bubble_dist, - drop_ratio, 0, 1, link); + drop_ratio, 0, 1, link, cov); } n_vtx = read_g->n_seq; @@ -21442,7 +21212,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link) if(asm_opt.recover_atg_cov_min == -1024) { asm_opt.recover_atg_cov_max = asm_opt.hom_global_coverage/HOM_PEAK_RATE; - asm_opt.recover_atg_cov_min = asm_opt.recover_atg_cov_max * 0.8; + asm_opt.recover_atg_cov_min = asm_opt.recover_atg_cov_max * 0.85; ///asm_opt.recover_atg_cov_max = asm_opt.recover_atg_cov_max * 1.2; asm_opt.recover_atg_cov_max = INT32_MAX; } @@ -21474,6 +21244,8 @@ kvec_asg_arc_t_warp* new_rtg_edges, hc_links* link) link->a.a[v].f.n = m; } } + + destory_hap_cov_t(&cov); } @@ -21495,9 +21267,9 @@ long long tipsLen, R_to_U* ruIndex, int max_hang, int min_ovlp) nsg->seq[v].c = PRIMARY_LABLE; EvaluateLen(ug->u, v) = ug->u.a[v].n; } - asg_pop_bubble_primary_trio(ug, bubble_dist, (uint32_t)-1, DROP); + asg_pop_bubble_primary_trio(ug, bubble_dist, (uint32_t)-1, DROP, NULL); cut_trio_tip_primary(ug->g, ug, tipsLen, (uint32_t)-1, 0, sg, reverse_sources, ruIndex, 2); - asg_pop_bubble_primary_trio(ug, bubble_dist, (uint32_t)-1, DROP); + asg_pop_bubble_primary_trio(ug, bubble_dist, (uint32_t)-1, DROP, NULL); cut_trio_tip_primary(ug->g, ug, tipsLen, (uint32_t)-1, 0, sg, reverse_sources, ruIndex, 2); delete_useless_nodes(&ug); renew_utg(&ug, sg, &new_rtg_edges); @@ -22724,7 +22496,7 @@ void lable_all_bubbles(asg_t *r_g, long long bubble_dist) ///if this is a bubble ///if(asg_bub_finder_with_del_advance(r_g, v, bubble_dist, &b) == 1) - if(asg_bub_pop1_primary_trio(r_g, NULL, v, bubble_dist, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL)) + if(asg_bub_pop1_primary_trio(r_g, NULL, v, bubble_dist, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL)) { //beg is v, end is b.S.a[0] //note b.b include end, does not include beg @@ -27358,6 +27130,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g) // rescue_no_coverage_aggressive(sg, sources, reverse_sources, &coverage_cut, ruIndex, max_hang_length, // mini_overlap_length, bubble_dist, 10); + rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, bubble_dist, + (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz); if (asm_opt.flag & HA_F_VERBOSE_GFA) { @@ -27370,8 +27144,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g) if (ha_opt_triobin(&asm_opt) && ha_opt_hic(&asm_opt)) { - rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, bubble_dist, - (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz); + // rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, bubble_dist, + // (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz); char *buf = (char*)calloc(strlen(output_file_name) + 25, 1); sprintf(buf, "%s.hic.bench", output_file_name); @@ -27381,8 +27155,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g) } else if (ha_opt_triobin(&asm_opt)) { - rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, bubble_dist, - (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz); + // rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, bubble_dist, + // (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz); char *buf = (char*)calloc(strlen(output_file_name) + 25, 1); sprintf(buf, "%s.dip", output_file_name); @@ -27398,8 +27172,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g) } else if(ha_opt_hic(&asm_opt)) { - rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, bubble_dist, - (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz); + // rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, bubble_dist, + // (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz); char *buf = (char*)calloc(strlen(output_file_name) + 25, 1); sprintf(buf, "%s.hic", output_file_name); @@ -27419,8 +27193,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g) output_contig_graph_primary_pre(sg, coverage_cut, output_file_name, sources, reverse_sources, asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length); - rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, bubble_dist, - (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz); + // rescue_bubble_by_chain(sg, coverage_cut, sources, reverse_sources, bubble_dist, + // (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, 10, gap_fuzz); output_contig_graph_primary(sg, coverage_cut, output_file_name, sources, reverse_sources, bubble_dist, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, diff --git a/Overlaps.h b/Overlaps.h index ca18270..3e3b0c2 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -478,7 +478,6 @@ void set_R_to_U(R_to_U* x, uint32_t rID, uint32_t uID, uint32_t is_Unitig, uint8 void get_R_to_U(R_to_U* x, uint32_t rID, uint32_t* uID, uint32_t* is_Unitig); void transfor_R_to_U(R_to_U* x); void debug_utg_graph(ma_ug_t *ug, asg_t* read_g, kvec_asg_arc_t_warp* edge, int require_equal_nv, int test_tangle); -int asg_pop_bubble_primary(asg_t *g, int max_dist); long long asg_arc_del_simple_circle_untig(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *g, long long circleLen, int is_drop); typedef struct { @@ -493,9 +492,35 @@ typedef struct { uint32_t new_edges_i; } Edge_iter; +typedef struct { + asg_arc_t x; + uint64_t Off; + uint64_t weight; +}asg_arc_t_offset; + +typedef struct { + kvec_t(asg_arc_t_offset) a; + uint64_t i; +}kvec_asg_arc_t_offset; + +typedef struct { + uint32_t n; + uint32_t* cov; + uint64_t* pos_idx; + ma_hit_t_alloc* reverse_sources; + ma_sub_t *coverage_cut; + R_to_U* ruIndex; + asg_t *read_g; + int max_hang; + int min_ovlp; + kvec_asg_arc_t_offset u_buffer; + kvec_t_i32_warp tailIndex; + kvec_t_i32_warp prevIndex; +}hap_cov_t; + void init_Edge_iter(asg_t* g, uint32_t v, asg_arc_t* new_edges, uint32_t new_edges_n, Edge_iter* x); int get_arc_t(Edge_iter* x, asg_arc_t* get); -int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_flag, uint32_t negative_flag); +int asg_pop_bubble_primary_trio(ma_ug_t *ug, int max_dist, uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov); inline int get_real_length(asg_t *g, uint32_t v, uint32_t* v_s) @@ -1040,10 +1065,10 @@ R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, uint32_t bac uint32_t is_bubble_check, uint32_t is_primary_check); uint32_t get_edge_from_source(ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t query, uint32_t target, asg_arc_t* t); -uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, int max_dist, buf_t *b, -uint32_t positive_flag, uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes); +uint64_t asg_bub_pop1_primary_trio(asg_t *g, ma_ug_t *utg, uint32_t v0, int max_dist, buf_t *b, uint32_t positive_flag, +uint32_t negative_flag, uint32_t is_pop, uint64_t* path_base_len, uint64_t* path_nodes, hap_cov_t *cov); int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio); - +void asg_bub_backtrack_primary(asg_t *g, uint32_t v0, buf_t *b); typedef struct{ double weight; @@ -1084,11 +1109,6 @@ typedef struct{ void init_hc_links(hc_links* link, uint64_t ug_num, uint64_t r_num); void destory_hc_links(hc_links* link); -void clean_primary_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, -long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, -R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, uint32_t miniHapLen, -uint32_t miniBiGraph, float chimeric_rate, int is_final_clean, int just_bubble_pop, -float drop_ratio, hc_links* link); void adjust_utg_by_primary(ma_ug_t **ug, asg_t* read_g, float drop_rate, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, long long bubble_dist, long long tipsLen, float tip_drop_ratio, long long stops_threshold, @@ -1112,6 +1132,9 @@ inline int inter_interval(int a_s, int a_e, int b_s, int b_e, int* i_s, int* i_e return 1; } + + + #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index 7474a14..d4920be 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -25,18 +25,6 @@ KDQ_INIT(uint64_t) uint8_t debug_enable = 0; -typedef struct { - asg_arc_t x; - uint64_t Off; - uint64_t weight; -}asg_arc_t_offset; - -typedef struct { - kvec_t(asg_arc_t_offset) a; - uint64_t i; -}kvec_asg_arc_t_offset; - - typedef struct { uint64_t weight; uint32_t x_beg_pos; @@ -45,6 +33,7 @@ typedef struct { uint32_t y_end_pos; uint32_t index_beg; uint32_t index_end; + long long score; uint8_t rev; asg_arc_t t; }hap_candidates; @@ -73,6 +62,7 @@ typedef struct { uint32_t xUid; uint32_t yUid; uint32_t weight; + long long score; }hap_overlaps; typedef struct { @@ -115,6 +105,7 @@ typedef struct { float chain_rate; hap_overlaps_list* all_ovlp; long long cov_threshold; + hap_cov_t *cov; }hap_alignment_struct_pip; @@ -492,7 +483,7 @@ void destory_hap_alignment_struct(hap_alignment_struct* x) void init_hap_alignment_struct_pip(hap_alignment_struct_pip* x, uint32_t num_threads, uint32_t n_seq, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, ma_sub_t *coverage_cut, -uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chain_rate, hap_overlaps_list* all_ovlp) +uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chain_rate, hap_overlaps_list* all_ovlp, hap_cov_t *cov) { uint32_t i; x->num_threads = num_threads; @@ -514,6 +505,7 @@ uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chai x->min_ovlp = min_ovlp; x->chain_rate = chain_rate; x->all_ovlp = all_ovlp; + x->cov = cov; } @@ -857,10 +849,127 @@ uint64_t get_pair_hap_coverage(uint64_t* readIDs, uint32_t Len, ma_hit_t_alloc* return C_bases/R_bases; } + +uint64_t get_pair_purge_coverage(ma_utg_t *xReads, long long xPosBeg, long long xPosEnd, +ma_utg_t *yReads, long long yPosBeg, long long yPosEnd, uint32_t rev, asg_t *read_g, hap_cov_t *cov) +{ + long long offset, r_beg, r_end, i_beg, i_end, ovlp, IdxBeg, IdxEnd; + uint64_t i, rId, uCov, uLen; + ma_utg_t *x = NULL; + uCov = uLen = 0; + if(rev) + { + yPosBeg = yReads->len - yPosBeg - 1; + yPosEnd = yReads->len - yPosEnd - 1; + offset = yPosBeg; yPosBeg = yPosEnd; yPosEnd = offset; + } + + + IdxBeg = IdxEnd = -1; + x = xReads; i_beg = xPosBeg; i_end = xPosEnd; + for (i = 0, offset = 0; i < x->n; i++) + { + rId = x->a[i]>>33; + r_beg = offset; r_end = offset + (long long)(read_g->seq[rId].len) - 1; + offset += (uint32_t)x->a[i]; + ovlp = (long long)(MIN(r_end, i_end)) - (long long)(MAX(r_beg, i_beg)) + 1; + if(ovlp <= 0 || ovlp < read_g->seq[rId].len * 0.8) + { + if(IdxBeg != -1 && IdxEnd != -1) break; + continue; + } + + if(IdxBeg == -1) IdxBeg = i; + IdxEnd = i; + } + if(IdxBeg != -1 && IdxEnd != -1) + { + for (i = IdxBeg; (long long)i <= IdxEnd; i++) + { + rId = x->a[i]>>33; + uCov += cov->cov[rId]; + uLen += cov->read_g->seq[rId].len; + } + } + + + + IdxBeg = IdxEnd = -1; + x = yReads; i_beg = yPosBeg; i_end = yPosEnd; + for (i = 0, offset = 0; i < x->n; i++) + { + rId = x->a[i]>>33; + r_beg = offset; r_end = offset + (long long)(read_g->seq[rId].len) - 1; + offset += (uint32_t)x->a[i]; + ovlp = (long long)(MIN(r_end, i_end)) - (long long)(MAX(r_beg, i_beg)) + 1; + if(ovlp <= 0 || ovlp < read_g->seq[rId].len * 0.8) + { + if(IdxBeg != -1 && IdxEnd != -1) break; + continue; + } + + if(IdxBeg == -1) IdxBeg = i; + IdxEnd = i; + } + if(IdxBeg != -1 && IdxEnd != -1) + { + for (i = IdxBeg; (long long)i <= IdxEnd; i++) + { + rId = x->a[i]>>33; + uCov += cov->cov[rId]; + uLen += cov->read_g->seq[rId].len; + } + } + + return (uLen == 0? 0 : uCov / uLen); +} + + +void get_pair_hap_similarity_by_base(ma_utg_t *xReads, asg_t *read_g, uint32_t target_uId, +ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, long long xBegPos, long long xEndPos, +double* Match, double* Total) +{ + uint32_t i, j, qn, tn, is_Unitig, uId, min_count = 0, max_count = 0; + long long offset, r_beg, r_end, ovlp; + + for (i = 0, offset = 0; i < xReads->n; i++) + { + qn = xReads->a[i]>>33; + r_beg = offset; r_end = offset + (long long)(read_g->seq[qn].len) - 1; + offset += (uint32_t)xReads->a[i]; + + ovlp = (long long)(MIN(r_end, xEndPos)) - (long long)(MAX(r_beg, xBegPos)) + 1; + if(ovlp <= 0) continue; + + if(reverse_sources[qn].length > 0) min_count++; + if(reverse_sources[qn].length == 0) continue; + for (j = 0; j < reverse_sources[qn].length; j++) + { + tn = Get_tn(reverse_sources[qn].buffer[j]); + if(read_g->seq[tn].del == 1) + { + get_R_to_U(ruIndex, tn, &tn, &is_Unitig); + if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue; + } + + + get_R_to_U(ruIndex, tn, &uId, &is_Unitig); + if(uId!=(uint32_t)-1 && is_Unitig == 1 && uId == target_uId) + { + max_count++; + break; + } + } + } + + (*Match) = max_count; + (*Total) = min_count; +} + void get_pair_hap_similarity(uint64_t* readIDs, uint32_t Len, uint32_t target_uId, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match, double* Total) { - #define CUTOFF_THRES 100 + #define CUTOFF_THRES 1000 uint32_t i, j, qn, tn, is_Unitig, uId, min_count = 0, max_count = 0, cutoff = 0;; for (i = 0; i < Len; i++) { @@ -872,6 +981,7 @@ ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match, } qn = readIDs[i]>>33; if(reverse_sources[qn].length > 0) min_count++; + if(reverse_sources[qn].length == 0) continue; for (j = 0; j < reverse_sources[qn].length; j++) { tn = Get_tn(reverse_sources[qn].buffer[j]); @@ -1028,7 +1138,7 @@ uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U void determin_hap_alignment_boundary_single_side(uint64_t* readIDs, long long queryLen, long long targetBeg, long long targetEnd, long long targetID, long long eMatch, long long eTotal, long long dir, -float Hap_rate, uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, +float H_rate, int is_local, uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_count) { if(queryLen == 0) @@ -1037,19 +1147,27 @@ R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_co return; } long long i, maxId, min_count = eTotal, max_count = eMatch, matchLen = 0; + long long rLen, score = 0, max_score = 0; uint32_t is_found, is_match; if(dir == 0) { for (i = 0, maxId = 0; i < queryLen; i++) { + check_hap_match(readIDs[i]>>33, targetBeg, targetEnd, targetID, position_index, reverse_sources, read_g, ruIndex, &is_found, &is_match); min_count += is_found; max_count += is_match; - if(max_count > min_count*Hap_rate) maxId = i; - } + if(max_count > min_count*H_rate) maxId = i; + if(is_local && is_found) + { + rLen = read_g->seq[readIDs[i]>>33].len; + score += (is_match? rLen : (rLen*(-1))); + if(score >= max_score) max_score = score, maxId = i; + } + } for (i = maxId; i >= 0; i--) { @@ -1064,7 +1182,7 @@ R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_co min_count -= is_found; max_count -= is_match; } - + matchLen = i+1; } else @@ -1076,7 +1194,14 @@ R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_co min_count += is_found; max_count += is_match; - if(max_count > min_count*Hap_rate) maxId = i; + if(max_count > min_count*H_rate) maxId = i; + + if(is_local && is_found) + { + rLen = read_g->seq[readIDs[i]>>33].len; + score += (is_match? rLen : (rLen*(-1))); + if(score >= max_score) max_score = score, maxId = i; + } } for (i = maxId; i < queryLen; i++) @@ -1123,7 +1248,7 @@ long long* target_beg, long long* target_end) void bi_direction_hap_alignment_extention(ma_utg_t* xReads, uint32_t xLeftBeg, uint32_t xLeftLen, uint32_t xRightBeg, uint32_t xRightLen, uint32_t targetUid, uint32_t target_beg, uint32_t target_end, -float Hap_rate, uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, +float Hap_rate, int is_local, uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, uint32_t rev, long long* x_interval_beg, long long* x_interval_end) { if(rev) @@ -1136,27 +1261,27 @@ uint32_t rev, long long* x_interval_beg, long long* x_interval_end) uint32_t n_matchLenRight, x_max_countRight, x_min_countRight; n_matchLenLeft = x_max_countLeft = x_min_countLeft = 0; determin_hap_alignment_boundary_single_side(xReads->a+xLeftBeg, xLeftLen, - target_beg, target_end, targetUid, x_max_countLeft, x_min_countLeft, 1, Hap_rate, + target_beg, target_end, targetUid, x_max_countLeft, x_min_countLeft, 1, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, &n_matchLenLeft, &x_max_countLeft, &x_min_countLeft); n_matchLenRight = x_max_countRight = x_min_countRight = 0; determin_hap_alignment_boundary_single_side(xReads->a+xRightBeg, xRightLen, - target_beg, target_end, targetUid, x_max_countRight, x_min_countRight, 0, Hap_rate, + target_beg, target_end, targetUid, x_max_countRight, x_min_countRight, 0, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, &n_matchLenRight, &x_max_countRight, &x_min_countRight); if(x_max_countLeft >= x_max_countRight) { determin_hap_alignment_boundary_single_side(xReads->a+xRightBeg, xRightLen, - target_beg, target_end, targetUid, x_max_countLeft, x_min_countLeft, 0, Hap_rate, + target_beg, target_end, targetUid, x_max_countLeft, x_min_countLeft, 0, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, &n_matchLenRight, &x_max_countRight, &x_min_countRight); } else { determin_hap_alignment_boundary_single_side(xReads->a+xLeftBeg, xLeftLen, - target_beg, target_end, targetUid, x_max_countRight, x_min_countRight, 1, Hap_rate, + target_beg, target_end, targetUid, x_max_countRight, x_min_countRight, 1, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, &n_matchLenLeft, &x_max_countLeft, &x_min_countLeft); } @@ -1170,7 +1295,7 @@ uint32_t xLeftMatch, uint32_t xLeftTotal, uint32_t yLeftMatch, uint32_t yLeftTot uint32_t xRightMatch, uint32_t xRightTotal, uint32_t yRightMatch, uint32_t yRightTotal, uint32_t xLeftBeg, uint32_t xLeftLen, uint32_t yLeftBeg, uint32_t yLeftLen, uint32_t xRightBeg, uint32_t xRightLen, uint32_t yRightBeg, uint32_t yRightLen, -uint32_t xUid, uint32_t yUid, float Hap_rate, uint64_t* position_index, +uint32_t xUid, uint32_t yUid, float Hap_rate, int is_local, uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, uint32_t rev, long long* r_x_interval_beg, long long* r_x_interval_end, long long* r_y_interval_beg, long long* r_y_interval_end) @@ -1189,7 +1314,7 @@ long long* r_y_interval_beg, long long* r_y_interval_end) modify_target_interval(yLeftBeg, yLeftBeg+yLeftLen-1, yReads->n, &target_beg, &target_end); determin_hap_alignment_boundary_single_side(xReads->a+xLeftBeg, xLeftLen, /**yLeftBeg, yLeftBeg+yLeftLen-1,**/ target_beg, target_end, yUid, - x_max_count, x_min_count, 1, Hap_rate, position_index, reverse_sources, + x_max_count, x_min_count, 1, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, &n_matchLen, &x_max_count, &x_min_count); x_interval_beg = xLeftBeg + xLeftLen; x_interval_beg -= n_matchLen; @@ -1202,7 +1327,7 @@ long long* r_y_interval_beg, long long* r_y_interval_end) modify_target_interval(xRightBeg, xRightBeg+xRightLen-1, xReads->n, &target_beg, &target_end); determin_hap_alignment_boundary_single_side(yReads->a+yRightBeg, yRightLen, /**xRightBeg, xRightBeg+xRightLen-1,**/ target_beg, target_end, xUid, - y_max_count, y_min_count, rev, Hap_rate, position_index, reverse_sources, + y_max_count, y_min_count, rev, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, &n_matchLen, &y_max_count, &y_min_count); if(rev == 0) { @@ -1225,7 +1350,7 @@ long long* r_y_interval_beg, long long* r_y_interval_end) modify_target_interval(yRightBeg, yRightBeg+yRightLen-1, yReads->n, &target_beg, &target_end); determin_hap_alignment_boundary_single_side(xReads->a+xRightBeg, xRightLen, /**yRightBeg, yRightBeg+yRightLen-1,**/ target_beg, target_end, yUid, - x_max_count, x_min_count, 0, Hap_rate, position_index, reverse_sources, + x_max_count, x_min_count, 0, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, &n_matchLen, &x_max_count, &x_min_count); x_interval_beg = xLeftBeg; @@ -1238,7 +1363,7 @@ long long* r_y_interval_beg, long long* r_y_interval_end) modify_target_interval(xLeftBeg, xLeftBeg+xLeftLen-1, xReads->n, &target_beg, &target_end); determin_hap_alignment_boundary_single_side(yReads->a+yLeftBeg, yLeftLen, /**xLeftBeg, xLeftBeg+xLeftLen-1,**/ target_beg, target_end, xUid, - y_max_count, y_min_count, 1-rev, Hap_rate, position_index, reverse_sources, + y_max_count, y_min_count, 1-rev, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, &n_matchLen, &y_max_count, &y_min_count); if(rev == 0) { @@ -1256,7 +1381,7 @@ long long* r_y_interval_beg, long long* r_y_interval_end) { /********************x*********************/ bi_direction_hap_alignment_extention(xReads, xLeftBeg, xLeftLen, xRightBeg, xRightLen, - yUid, 0, yReads->n - 1, Hap_rate, position_index, reverse_sources, read_g, ruIndex, 0, + yUid, 0, yReads->n - 1, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, 0, &x_interval_beg, &x_interval_end); /********************x*********************/ @@ -1274,7 +1399,7 @@ long long* r_y_interval_beg, long long* r_y_interval_end) /********************y*********************/ bi_direction_hap_alignment_extention(yReads, yLeftBeg, yLeftLen, yRightBeg, yRightLen, - xUid, 0, xReads->n - 1, Hap_rate, position_index, reverse_sources, read_g, ruIndex, rev, + xUid, 0, xReads->n - 1, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, rev, &y_interval_beg, &y_interval_end); /********************y*********************/ } else abort(); @@ -1506,7 +1631,7 @@ void quick_LIS(asg_arc_t_offset* x, uint32_t n, kvec_t_i32_warp* tailIndex, kvec if(Get_yOff(x[i].Off) < Get_yOff(x[tailIndex->a.a[0]].Off)) { // new smallest value - tailIndex->a.a[0] = i; + tailIndex->a.a[0] = i; ///doesn't matter too much } else if(Get_yOff(x[i].Off) > Get_yOff(x[tailIndex->a.a[len - 1]].Off)) { @@ -1534,7 +1659,309 @@ void quick_LIS(asg_arc_t_offset* x, uint32_t n, kvec_t_i32_warp* tailIndex, kvec tailIndex->a.n = len; } -void get_base_boundary_advance(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, +inline uint64_t get_xy_pos_by_pos(asg_t *read_g, asg_arc_t* t, uint32_t v_in_unitig, uint32_t w_in_unitig, +uint32_t v_in_pos, uint32_t w_in_pos, uint32_t xUnitigLen, uint32_t yUnitigLen, uint8_t* rev) +{ + uint32_t x_pos, y_pos, x_dir = 0, y_dir = 0; + uint64_t tmp; + x_pos = y_pos = (uint32_t)-1; + if((t->ul>>32)==v_in_unitig)///end pos + { + x_pos = v_in_pos + read_g->seq[v_in_unitig>>1].len - 1; + x_dir = 0; + } + else if((t->ul>>32)==(v_in_unitig^1))///start pos + { + x_pos = v_in_pos; + x_dir = 1; + } + else + { + fprintf(stderr, "ERROR\n"); + } + + if(t->v == w_in_unitig) + { + y_pos = w_in_pos + t->ol - 1; + y_dir = 0; + } + else if(t->v == (w_in_unitig^1)) + { + y_pos = w_in_pos + read_g->seq[w_in_unitig>>1].len - t->ol; + y_dir = 1; + } + else + { + fprintf(stderr, "ERROR\n"); + } + + (*rev) = x_dir^y_dir; + if((*rev)) + { + if(yUnitigLen <= y_pos) + { + y_pos = (uint32_t)-1; + } + else + { + y_pos = yUnitigLen - y_pos - 1; + } + } + + if(x_pos>=xUnitigLen) x_pos = (uint32_t)-1; + if(y_pos>=yUnitigLen) y_pos = (uint32_t)-1; + + tmp = x_pos; tmp = tmp << 32; tmp = tmp | y_pos; + return tmp; +} + +void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd) +{ + ma_hit_t_alloc* reverse_sources = cov->reverse_sources; + ma_sub_t *coverage_cut = cov->coverage_cut; + int max_hang = cov->max_hang; + int min_ovlp = cov->min_ovlp; + kvec_asg_arc_t_offset* u_buffer = &(cov->u_buffer); + kvec_t_i32_warp* tailIndex = &(cov->tailIndex); + kvec_t_i32_warp* prevIndex = &(cov->prevIndex); + ma_hit_t_alloc *xR = NULL; + ma_hit_t *h = NULL; + ma_sub_t *sq = NULL, *st = NULL; + int32_t r; + asg_arc_t t; + uint32_t rId, v, w; + uint64_t tmp; + asg_arc_t_offset t_offset; + u_buffer->a.n = 0; + (*xEnd) = (uint32_t)-1; + uint32_t u_i, r_i, k, j, m, len, p_v, *a = xReads->b.a, uid, ori, l, aOcc, nv, xOcc = (uint32_t)-1; + ma_utg_t* u = NULL; + asg_arc_t *av = NULL; + + + for (u_i = r_i = len = aOcc = 0, xOcc = (uint32_t)-1, p_v = (uint32_t)-1; u_i < xReads->b.n; u_i++) + { + uid = a[u_i] >> 1; + ori = a[u_i] & 1; + u = &(ug->u.a[uid]); + if(u->n == 0) continue; + + for (r_i = 0; r_i < u->n; r_i++, aOcc++) + { + l = 0; + v = (ori == 1?((uint64_t)((u->a[u->n - r_i - 1])^(uint64_t)(0x100000000)))>>32:((uint64_t)(u->a[r_i]))>>32); + + if(p_v != (uint32_t)-1) + { + av = asg_arc_a(read_sg, p_v); + nv = asg_arc_n(read_sg, p_v); + + for (k = 0; k < nv; k++) + { + if(av[k].del) continue; + if(av[k].v == v) + { + l = asg_arc_len(av[k]); + break; + } + } + if(k == nv) fprintf(stderr, "ERROR\n"); + } + + p_v = v; len += l; + if(len >= targetBaseLen) + { + xOcc = aOcc; + break; + } + } + + if(xOcc != (uint32_t)-1) break; + } + if(xOcc == (uint32_t)-1) xOcc = aOcc; + if(xOcc == 0) xOcc = 1; + + + + for (u_i = r_i = len = aOcc = 0, p_v = (uint32_t)-1; u_i < xReads->b.n; u_i++) + { + uid = a[u_i] >> 1; + ori = a[u_i] & 1; + u = &(ug->u.a[uid]); + if(u->n == 0) continue; + + for (r_i = 0; r_i < u->n; r_i++, aOcc++) + { + if(aOcc >= xOcc) break; + l = 0; + v = (ori == 1?((uint64_t)((u->a[u->n - r_i - 1])^(uint64_t)(0x100000000)))>>32:((uint64_t)(u->a[r_i]))>>32); + + if(p_v != (uint32_t)-1) + { + av = asg_arc_a(read_sg, p_v); + nv = asg_arc_n(read_sg, p_v); + + for (k = 0; k < nv; k++) + { + if(av[k].del) continue; + if(av[k].v == v) + { + l = asg_arc_len(av[k]); + break; + } + } + if(k == nv) fprintf(stderr, "ERROR\n"); + } + + p_v = v; len += l; + + xR = &(reverse_sources[v>>1]); + for (j = 0; j < xR->length; j++) + { + h = &(xR->buffer[j]); + sq = &(coverage_cut[Get_qn(*h)]); + st = &(coverage_cut[Get_tn(*h)]); + if(st->del || read_sg->seq[Get_tn(*h)].del) continue; + r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, + asm_opt.max_hang_rate, min_ovlp, &t); + ///if it is a contained overlap, skip + if(r < 0) continue; + + rId = t.v>>1; + if(read_sg->seq[rId].del == 1) continue; + if(cov->pos_idx[rId] == (uint64_t)-1) continue; + w = (uint32_t)(cov->pos_idx[rId]); + if(rId != (w>>1)) continue; + + tmp = get_xy_pos_by_pos(read_sg, &t, v, w, len, cov->pos_idx[w>>1]>>32, + (uint32_t)-1, targetBaseLen, &(t.el)); + if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue; + if(t.el) continue; ///must + + t_offset.Off = tmp; + t_offset.x = t; + t_offset.weight = 1; + kv_push(asg_arc_t_offset, u_buffer->a, t_offset); + } + } + + if(aOcc >= xOcc) break; + } + + if(u_buffer->a.n == 0) return; + + qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment_chaining); + + ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "before"); + + + for (k = 1, l = 0, m = 0; k <= u_buffer->a.n; ++k) + { + if (k == u_buffer->a.n || u_buffer->a.a[k].x.el != u_buffer->a.a[l].x.el || + u_buffer->a.a[k].Off != u_buffer->a.a[l].Off) + { + u_buffer->a.a[m] = u_buffer->a.a[l]; + for (l += 1; l < k; l++) + { + u_buffer->a.a[m].weight += u_buffer->a.a[l].weight; + if(u_buffer->a.a[l].x.ol > u_buffer->a.a[m].x.ol) + { + u_buffer->a.a[m].x = u_buffer->a.a[l].x; + } + } + l = k; + m++; + } + } + u_buffer->a.n = m; + + ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "after"); + quick_LIS(u_buffer->a.a, u_buffer->a.n, tailIndex, prevIndex); + if(tailIndex->a.n == 0) return; + + + uint32_t xLen_thres = (uint32_t)-1; + asg_arc_t_offset* best = &(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]]); + for (u_i = r_i = len = aOcc = 0, p_v = (uint32_t)-1; u_i < xReads->b.n; u_i++) + { + uid = a[u_i] >> 1; + ori = a[u_i] & 1; + u = &(ug->u.a[uid]); + if(u->n == 0) continue; + + for (r_i = 0; r_i < u->n; r_i++, aOcc++) + { + l = 0; + v = (ori == 1?((uint64_t)((u->a[u->n - r_i - 1])^(uint64_t)(0x100000000)))>>32:((uint64_t)(u->a[r_i]))>>32); + + if(p_v != (uint32_t)-1) + { + av = asg_arc_a(read_sg, p_v); + nv = asg_arc_n(read_sg, p_v); + + for (k = 0; k < nv; k++) + { + if(av[k].del) continue; + if(av[k].v == v) + { + l = asg_arc_len(av[k]); + break; + } + } + if(k == nv) fprintf(stderr, "ERROR\n"); + } + + p_v = v; len += l; + + if((v>>1) == (best->x.ul>>33) && xLen_thres == (uint32_t)-1) + { + ///cov->pos_idx[v>>1] = len; + xR = &(reverse_sources[v>>1]); + for (j = 0; j < xR->length; j++) + { + h = &(xR->buffer[j]); + sq = &(coverage_cut[Get_qn(*h)]); + st = &(coverage_cut[Get_tn(*h)]); + if(st->del || read_sg->seq[Get_tn(*h)].del) continue; + r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, + asm_opt.max_hang_rate, min_ovlp, &t); + ///if it is a contained overlap, skip + if(r < 0) continue; + + rId = t.v>>1; + if(read_sg->seq[rId].del == 1) continue; + if(cov->pos_idx[rId] == (uint64_t)-1) continue; + w = (uint32_t)(cov->pos_idx[rId]); + if(rId != (w>>1)) continue; + + tmp = get_xy_pos_by_pos(read_sg, &t, v, w, len, cov->pos_idx[w>>1]>>32, + (uint32_t)-1, targetBaseLen, &(t.el)); + if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue; + if(t.el) continue; ///must + + t_offset.Off = tmp; + t_offset.x = t; + t_offset.weight = 1; + if(t_offset.Off == best->Off && t_offset.x.v == best->x.v && t_offset.x.ul == best->x.ul) + { + xLen_thres = Get_xOff(best->Off) + targetBaseLen - Get_yOff(best->Off); + } + } + } + + if(len >= xLen_thres) + { + (*xEnd) = aOcc; + return; + } + } + } + + (*xEnd) = aOcc; +} + + +void get_base_boundary_advance_back(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads, uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, long long yEndIndex, uint32_t rev, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, @@ -1635,15 +2062,14 @@ uint32_t* xBeg, uint32_t* xEnd, uint32_t* yBeg, uint32_t* yEnd) (*yEnd) = Get_yOff(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].Off); } - -uint32_t determine_hap_overlap_type_advance(hap_candidates* hap_can, ma_utg_t *xReads, ma_utg_t *yReads, +uint32_t determine_hap_overlap_type_advance_back(hap_candidates* hap_can, ma_utg_t *xReads, ma_utg_t *yReads, R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, uint32_t xUid, uint32_t yUid, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end) { uint32_t x_pos_beg, y_pos_beg, x_pos_end, y_pos_end; /*************************x***************************/ - get_base_boundary_advance(ruIndex, reverse_sources, coverage_cut, read_g, position_index, + get_base_boundary_advance_back(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); @@ -1666,6 +2092,321 @@ kvec_t_i32_warp* prevIndex, long long* r_x_pos_beg, long long* r_x_pos_end, long r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end); } +void get_idx_by_base(ma_utg_t *x, asg_t *read_g, long long beg_base, long long end_base, + long long* beg_idx, long long* end_idx) +{ + long long offset, r_beg, r_end; + uint64_t i, rId; + (*beg_idx) = (*end_idx) = -1; + for (i = 0, offset = 0; i < x->n; i++) + { + rId = x->a[i]>>33; + r_beg = offset; r_end = offset + (long long)(read_g->seq[rId].len) - 1; + offset += (uint32_t)x->a[i]; + if(beg_base > r_end || r_beg > end_base) + { + if((*beg_idx) != -1 && (*end_idx) != -1) break; + continue; + } + if((*beg_idx) == -1) (*beg_idx) = i; + (*end_idx) = i; + } +} + +int get_base_boundary_chain(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, +asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads, +uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, long long yEndIndex, +uint32_t rev, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex) +{ + long long k, j, l, offset, m; + ma_hit_t_alloc *xR = NULL; + ma_hit_t *h = NULL; + ma_sub_t *sq = NULL, *st = NULL; + int32_t r; + asg_arc_t t; + uint32_t rId, Hap_uId, is_Unitig, v, w, v_dir, w_dir; + uint64_t tmp; + asg_arc_t_offset t_offset; + u_buffer->a.n = 0; + for (k = xBegIndex; k <= xEndIndex; k++) + { + xR = &(reverse_sources[xReads->a[k]>>33]); + for (j = 0; j < xR->length; j++) + { + h = &(xR->buffer[j]); + sq = &(coverage_cut[Get_qn(*h)]); + st = &(coverage_cut[Get_tn(*h)]); + if(st->del || read_g->seq[Get_tn(*h)].del) continue; + + r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, + asm_opt.max_hang_rate, min_ovlp, &t); + ///if it is a contained overlap, skip + if(r < 0) continue; + + rId = t.v>>1; + if(read_g->seq[rId].del == 1) continue; + ///there are two cases: + ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id + ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1 + get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig); + if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; + if(Hap_uId != yUid) continue; + + v = xReads->a[k]>>32; + get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig); + if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; + if(Hap_uId != xUid) continue; + if((uint32_t)(position_index[v>>1]) != k) continue; + + w = (yReads->a[(uint32_t)(position_index[rId])])>>32; + + v_dir = ((t.ul>>32)==v)?1:0; + w_dir = (t.v == w)?1:0; + if(rev == 0 && v_dir != w_dir) continue; + if(rev == 1 && v_dir == w_dir) continue; + + /****************************may have bugs********************************/ + offset = (uint32_t)(position_index[rId]); + if(offset < yBegIndex || offset > yEndIndex) continue; + /****************************may have bugs********************************/ + + tmp = get_xy_pos(read_g, &t, v, w, xReads->len, yReads->len, position_index, &(t.el)); + if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue; + + t_offset.Off = tmp; + t_offset.x = t; + t_offset.weight = 1; + kv_push(asg_arc_t_offset, u_buffer->a, t_offset); + } + } + if(u_buffer->a.n == 0) return 0; + + qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment_chaining); + + ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "before"); + for (k = 1, l = 0, m = 0; k <= (long long)u_buffer->a.n; ++k) + { + if (k == (long long)u_buffer->a.n || u_buffer->a.a[k].Off != u_buffer->a.a[l].Off) + { + u_buffer->a.a[m] = u_buffer->a.a[l]; + for (l += 1; l < k; l++) + { + u_buffer->a.a[m].weight += u_buffer->a.a[l].weight; + if(u_buffer->a.a[l].x.ol > u_buffer->a.a[m].x.ol) + { + u_buffer->a.a[m].x = u_buffer->a.a[l].x; + } + } + l = k; + m++; + } + } + u_buffer->a.n = m; + + ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "after"); + quick_LIS(u_buffer->a.a, u_buffer->a.n, tailIndex, prevIndex); + + if(tailIndex->a.n == 0) return 0; + return 1; +} + + +void get_base_boundary_advance(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, +asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads, +uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, long long yEndIndex, +uint32_t rev, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, +uint32_t* xBeg, uint32_t* xEnd, uint32_t* yBeg, uint32_t* yEnd) +{ + long long offset; + long long new_xBeg, new_yBeg, new_xEnd, new_yEnd; + long long new_xIdxBeg, new_yIdxBeg, new_xIdxEnd, new_yIdxEnd; + + (*xBeg) = (*xEnd) = (*yBeg) = (*yEnd) = (uint32_t)-1; + if(!get_base_boundary_chain(ruIndex, reverse_sources, coverage_cut, read_g, position_index, + max_hang, min_ovlp, xReads, yReads, xUid, yUid, xBegIndex, xEndIndex, yBegIndex, yEndIndex, + rev, u_buffer, tailIndex, prevIndex)) + { + return; + } + ///base + new_xBeg = Get_xOff(u_buffer->a.a[tailIndex->a.a[0]].Off); + new_yBeg = Get_yOff(u_buffer->a.a[tailIndex->a.a[0]].Off); + new_xEnd = Get_xOff(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].Off); + new_yEnd = Get_yOff(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].Off); + + if(new_xBeg > new_xEnd || new_yBeg > new_yEnd) return; + + classify_hap_overlap(new_xBeg, new_xEnd, xReads->len, new_yBeg, new_yEnd, yReads->len, + &new_xBeg, &new_xEnd, &new_yBeg, &new_yEnd); + + if(rev) + { + new_yBeg = yReads->len - new_yBeg - 1; + new_yEnd = yReads->len - new_yEnd - 1; + offset = new_yBeg; new_yBeg = new_yEnd; new_yEnd = offset; + } + ///idx + get_idx_by_base(xReads, read_g, new_xBeg, new_xEnd, &new_xIdxBeg, &new_xIdxEnd); + get_idx_by_base(yReads, read_g, new_yBeg, new_yEnd, &new_yIdxBeg, &new_yIdxEnd); + if(new_xIdxBeg == -1 || new_xIdxEnd == -1 || new_yIdxBeg == -1 || new_yIdxEnd == -1) return; + + if(!get_base_boundary_chain(ruIndex, reverse_sources, coverage_cut, read_g, position_index, + max_hang, min_ovlp, xReads, yReads, xUid, yUid, new_xIdxBeg, new_xIdxEnd, new_yIdxBeg, + new_yIdxEnd, rev, u_buffer, tailIndex, prevIndex)) + { + return; + } + + (*xBeg) = Get_xOff(u_buffer->a.a[tailIndex->a.a[0]].Off); + (*yBeg) = Get_yOff(u_buffer->a.a[tailIndex->a.a[0]].Off); + (*xEnd) = Get_xOff(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].Off); + (*yEnd) = Get_yOff(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].Off); +} + +#define generic_key(x) (x) +KRADIX_SORT_INIT(i32, int32_t, generic_key, sizeof(int32_t)) + +long long get_chain_score(ma_utg_t *xReads, asg_t *read_g, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* idx, +ma_hit_t_alloc* reverse_sources, long long xBegPos, long long xEndPos) +{ + long long offset, r_beg, r_end, inp_beg, inp_end, hap_beg, hap_end, inp_match, hap_match, ovlp; + uint64_t i, k, rId; + idx->a.n = 0; + + for (i = k = 0; i < tailIndex->a.n; i++) + { + rId = u_buffer->a.a[tailIndex->a.a[i]].x.ul>>33; + + + for (; k < xReads->n; k++) + { + if(rId == (xReads->a[k]>>33)) break; + } + + if(k >= xReads->n) + { + for (k = 0; k < xReads->n; k++) + { + if(rId == (xReads->a[k]>>33)) break; + } + } + + if(k < xReads->n) kv_push(int32_t, idx->a, k); + else + { + fprintf(stderr, "\nERROR-get_chain_score: tailIndex->a.n: %lu, xReads->n: %lu\n", (uint64_t)tailIndex->a.n, (uint64_t)xReads->n); + } + } + + + radix_sort_i32(idx->a.a, idx->a.a + idx->a.n); + inp_beg = -1; inp_end = -2; + hap_beg = -1; hap_end = -2; + for (i = k = 0, offset = 0, inp_match = hap_match = 0; i < xReads->n; i++) + { + rId = xReads->a[i]>>33; + r_beg = offset; r_end = offset + (long long)(read_g->seq[rId].len) - 1; + offset += (uint32_t)xReads->a[i]; + + if(reverse_sources[rId].length > 0) + { + if(r_beg <= hap_end) + { + hap_end = MAX(hap_end, r_end); + } + else + { + ///match += (hap_end - hap_beg + 1); + ovlp = (long long)(MIN(hap_end, xEndPos)) - (long long)(MAX(hap_beg, xBegPos)) + 1; + hap_match += (ovlp >= 0? ovlp : 0); + hap_beg = r_beg; hap_end = r_end; + } + } + + for (; k < idx->a.n; k++) + { + if(i <= (uint64_t)idx->a.a[k]) break; + } + + if(k >= idx->a.n) continue; + + if(i == (uint64_t)idx->a.a[k]) + { + if(r_beg <= inp_end) + { + inp_end = MAX(inp_end, r_end); + } + else + { + ovlp = (long long)(MIN(inp_end, xEndPos)) - (long long)(MAX(inp_beg, xBegPos)) + 1; + inp_match += (ovlp >= 0? ovlp : 0); + inp_beg = r_beg; inp_end = r_end; + } + } + } + + ovlp = (long long)(MIN(inp_end, xEndPos)) - (long long)(MAX(inp_beg, xBegPos)) + 1; + inp_match += (ovlp >= 0? ovlp : 0); + + ovlp = (long long)(MIN(hap_end, xEndPos)) - (long long)(MAX(hap_beg, xBegPos)) + 1; + hap_match += (ovlp >= 0? ovlp : 0); + + // if(inp_match > (xEndPos - xBegPos + 1)) fprintf(stderr, "ERROR1\n"); + // if(hap_match > (xEndPos - xBegPos + 1)) fprintf(stderr, "ERRO2\n"); + // if(inp_match > hap_match) fprintf(stderr, "ERROR3\n"); + // fprintf(stderr, "tailIndex->a.n: %u, xReads->n: %u, total_match: %lld, hap_match: %lld, inp_match: %lld\n", + // tailIndex->a.n, xReads->n, (xEndPos - xBegPos + 1), hap_match, inp_match); + + return (inp_match*3) - ((hap_match-inp_match)*1); +} + + +uint32_t determine_hap_overlap_type_advance(hap_candidates* hap_can, ma_utg_t *xReads, ma_utg_t *yReads, +R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, uint64_t* position_index, +int max_hang, int min_ovlp, uint32_t xUid, uint32_t yUid, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, +kvec_t_i32_warp* prevIndex, long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end) +{ + uint32_t x_pos_beg, y_pos_beg, x_pos_end, y_pos_end; + /*************************x***************************/ + 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); + /*************************x***************************/ + if(x_pos_beg == (uint32_t)-1 || y_pos_beg == (uint32_t)-1 + || x_pos_end == (uint32_t)-1 || y_pos_end == (uint32_t)-1) + { + return (uint32_t)-1; + } + if(x_pos_beg > x_pos_end || y_pos_beg > y_pos_end) return (uint32_t)-1; + + /** + #define X2Y 0 + #define Y2X 1 + #define XCY 2 + #define YCX 3 + **/ + hap_can->index_end = classify_hap_overlap(x_pos_beg, x_pos_end, xReads->len, y_pos_beg, y_pos_end, yReads->len, + r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end); + hap_can->x_beg_pos = MIN((uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[0]].x.ul>>33]), + (uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].x.ul>>33])); + hap_can->x_end_pos = MAX((uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[0]].x.ul>>33]), + (uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].x.ul>>33])); + hap_can->y_beg_pos = MIN((uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[0]].x.v>>1]), + (uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].x.v>>1])); + hap_can->y_end_pos = MAX((uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[0]].x.v>>1]), + (uint32_t)(position_index[u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].x.v>>1])); + double xLeftMatch, xLeftTotal; + get_pair_hap_similarity(xReads->a + hap_can->x_beg_pos, hap_can->x_end_pos + 1 - hap_can->x_beg_pos, + yUid, reverse_sources, read_g, ruIndex, &xLeftMatch, &xLeftTotal); + if(xLeftMatch == 0 || xLeftTotal == 0) return (uint32_t)-1; + hap_can->weight = xLeftMatch; + hap_can->index_beg = xLeftTotal; + hap_can->score = get_chain_score(xReads, read_g, u_buffer, tailIndex, prevIndex, reverse_sources, + (*r_x_pos_beg), (*r_x_pos_end)); + return hap_can->index_end; +} + uint32_t determine_hap_overlap_type(hap_candidates* hap_can, ma_utg_t *xReads, ma_utg_t *yReads, @@ -1702,133 +2443,12 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon } -uint32_t calculate_pair_hap_similarity(kvec_asg_arc_t_offset* u_buffer, hap_candidates* hap_can, -uint64_t* position_index, uint32_t xUid, uint32_t yUid, ma_utg_t* xReads, ma_utg_t* yReads, -ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, ma_sub_t *coverage_cut, -float Hap_rate, int max_hang, int min_ovlp, long long* r_x_pos_beg, long long* r_x_pos_end, -long long* r_y_pos_beg, long long* r_y_pos_end) -{ - uint32_t max_count = 0, min_count = 0, i, flag; - uint32_t xLen = xReads->n, xIndex/**, xBasePos**/; - uint32_t yLen = yReads->n, yIndex/**, yBasePos**/; - uint32_t xLeftBeg, xLeftLen, yLeftBeg, yLeftLen; - uint32_t xRightBeg, xRightLen, yRightBeg, yRightLen; - uint64_t totalWeigth; - double xLeftMatch = 0, xLeftTotal = 0, yLeftMatch = 0, yLeftTotal = 0; - double xRightMatch = 0, xRightTotal = 0, yRightMatch = 0, yRightTotal = 0; - asg_arc_t_offset* arch = NULL; - - for (i = hap_can->index_beg, totalWeigth = 0; i <= hap_can->index_end; i++) - { - totalWeigth += u_buffer->a.a[i].weight; - if(totalWeigth >= (hap_can->weight/2)) break; - } - if(i > hap_can->index_end) i = hap_can->index_end; - - arch = &(u_buffer->a.a[i]); - xIndex = (uint32_t)(position_index[arch->x.ul>>33]); - yIndex = (uint32_t)(position_index[arch->x.v>>1]); - ///xBasePos = (uint32_t)(arch->Off>>32); - ///yBasePos = (uint32_t)(arch->Off); - - if(hap_can->rev == 0) - { - xLeftBeg = 0; xLeftLen = xIndex; xRightBeg = xIndex; xRightLen = xLen - xRightBeg; - yLeftBeg = 0; yLeftLen = yIndex; yRightBeg = yIndex; yRightLen = yLen - yRightBeg; - } - else - { - xLeftBeg = 0; xLeftLen = xIndex; xRightBeg = xIndex; xRightLen = xLen - xRightBeg; - - yLeftBeg = yIndex + 1; yLeftLen = yLen - yLeftBeg; - yRightBeg = 0; yRightLen = yIndex + 1; - } - - ///flag = classify_hap_overlap(xBasePos, xBasePos, xReads->len, yBasePos, yBasePos, yReads->len); - flag = vote_overlap_type(u_buffer, hap_can, position_index, xReads, yReads); - - - if(flag == XCY) - { - get_pair_hap_similarity(yReads->a, yLen, xUid, reverse_sources, read_g, ruIndex, - &yLeftMatch, &yLeftTotal); - max_count = yLeftMatch; - min_count = yLeftTotal; - } - else if(flag == YCX) - { - get_pair_hap_similarity(xReads->a, xLen, yUid, reverse_sources, read_g, ruIndex, - &xLeftMatch, &xLeftTotal); - max_count = xLeftMatch; - min_count = xLeftTotal; - } - else if(flag == X2Y) - { - get_pair_hap_similarity(yReads->a+yLeftBeg, yLeftLen, xUid, reverse_sources, read_g, ruIndex, - &yLeftMatch, &yLeftTotal); - get_pair_hap_similarity(xReads->a+xRightBeg, xRightLen, yUid, reverse_sources, read_g, ruIndex, - &xRightMatch, &xRightTotal); - max_count = yLeftMatch + xRightMatch; - min_count = yLeftTotal + xRightTotal; - } - else if(flag == Y2X) - { - get_pair_hap_similarity(xReads->a+xLeftBeg, xLeftLen, yUid, reverse_sources, read_g, ruIndex, - &xLeftMatch, &xLeftTotal); - get_pair_hap_similarity(yReads->a+yRightBeg, yRightLen, xUid, reverse_sources, read_g, ruIndex, - &yRightMatch, &yRightTotal); - max_count = xLeftMatch + yRightMatch; - min_count = xLeftTotal + yRightTotal; - } else abort(); - - hap_can->weight = hap_can->index_beg = 0; - if(min_count == 0) return NON_PLOID; - if(max_count > min_count*Hap_rate) - { - long long r_x_interval_beg, r_x_interval_end, r_y_interval_beg, r_y_interval_end; - - ///for containment, don't need to do anything - get_hap_alignment_boundary(xReads, yReads, flag, xLeftMatch, xLeftTotal, - yLeftMatch, yLeftTotal, xRightMatch, xRightTotal, yRightMatch, yRightTotal, - xLeftBeg, xLeftLen, yLeftBeg, yLeftLen, xRightBeg, xRightLen, yRightBeg, yRightLen, - xUid, yUid, Hap_rate, position_index, reverse_sources, read_g, ruIndex, hap_can->rev, - &r_x_interval_beg, &r_x_interval_end, &r_y_interval_beg, &r_y_interval_end); - - if(r_x_interval_beg < 0 || r_x_interval_end < 0 || r_y_interval_beg < 0 || r_y_interval_end < 0) - { - return NON_PLOID; - } - - get_pair_hap_similarity(xReads->a + r_x_interval_beg, r_x_interval_end + 1 - r_x_interval_beg, - yUid, reverse_sources, read_g, ruIndex, &xLeftMatch, &xLeftTotal); - if(xLeftMatch == 0 || xLeftTotal == 0) return NON_PLOID; - - hap_can->weight = xLeftMatch; - hap_can->index_beg = xLeftTotal; - hap_can->index_end = flag; - hap_can->x_beg_pos = r_x_interval_beg; - hap_can->x_end_pos = r_x_interval_end; - hap_can->y_beg_pos = r_y_interval_beg; - hap_can->y_end_pos = r_y_interval_end; - - hap_can->index_end = determine_hap_overlap_type(hap_can, xReads, yReads, ruIndex, - reverse_sources, coverage_cut, read_g, position_index, max_hang, min_ovlp, xUid, - yUid, r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end); - 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; - if(hap_can->index_end == (uint32_t)-1) return NON_PLOID; - - return PLOID; - } - return NON_PLOID; -} - - uint32_t calculate_pair_hap_similarity_advance(hap_candidates* hap_can, uint64_t* position_index, uint32_t xUid, uint32_t yUid, ma_utg_t* xReads, ma_utg_t* yReads, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, ma_sub_t *coverage_cut, -float Hap_rate, int max_hang, int min_ovlp, uint64_t cov_threshold, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, -long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end) +float Hap_rate, int is_local, int max_hang, int min_ovlp, uint64_t cov_threshold, kvec_asg_arc_t_offset* u_buffer, +kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, hap_cov_t *cov, long long* r_x_pos_beg, long long* r_x_pos_end, +long long* r_y_pos_beg, long long* r_y_pos_end) { uint32_t max_count = 0, min_count = 0, flag; uint32_t xLen = xReads->n, xIndex; @@ -1895,7 +2515,7 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon hap_can->weight = hap_can->index_beg = 0; if(min_count == 0) return NON_PLOID; - if(max_count > min_count*Hap_rate) + if((max_count > min_count*Hap_rate) || is_local) { long long r_x_interval_beg, r_x_interval_end, r_y_interval_beg, r_y_interval_end; uint64_t ploid_coverage = 0; @@ -1904,8 +2524,8 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon get_hap_alignment_boundary(xReads, yReads, flag, xLeftMatch, xLeftTotal, yLeftMatch, yLeftTotal, xRightMatch, xRightTotal, yRightMatch, yRightTotal, xLeftBeg, xLeftLen, yLeftBeg, yLeftLen, xRightBeg, xRightLen, yRightBeg, yRightLen, - xUid, yUid, Hap_rate, position_index, reverse_sources, read_g, ruIndex, hap_can->rev, - &r_x_interval_beg, &r_x_interval_end, &r_y_interval_beg, &r_y_interval_end); + xUid, yUid, Hap_rate, is_local, position_index, reverse_sources, read_g, ruIndex, + hap_can->rev, &r_x_interval_beg, &r_x_interval_end, &r_y_interval_beg, &r_y_interval_end); if(r_x_interval_beg < 0 || r_x_interval_end < 0 || r_y_interval_beg < 0 || r_y_interval_end < 0) { @@ -1914,7 +2534,11 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon get_pair_hap_similarity(xReads->a + r_x_interval_beg, r_x_interval_end + 1 - r_x_interval_beg, yUid, reverse_sources, read_g, ruIndex, &xLeftMatch, &xLeftTotal); - if(xLeftMatch == 0 || xLeftTotal == 0) return NON_PLOID; + if(xLeftMatch == 0 || xLeftTotal == 0 || (is_local == 0 && xLeftMatch <= xLeftTotal*Hap_rate)) + { + return NON_PLOID; + } + hap_can->weight = xLeftMatch; hap_can->index_beg = xLeftTotal; @@ -1923,25 +2547,26 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon hap_can->x_end_pos = r_x_interval_end; hap_can->y_beg_pos = r_y_interval_beg; hap_can->y_end_pos = r_y_interval_end; - 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); - 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; if(hap_can->index_end == (uint32_t)-1) return NON_PLOID; - ploid_coverage = 0; - ploid_coverage += get_pair_hap_coverage(xReads->a+r_x_interval_beg, r_x_interval_end+1-r_x_interval_beg, - sources, coverage_cut); - ploid_coverage += get_pair_hap_coverage(yReads->a+r_y_interval_beg, r_y_interval_end+1-r_y_interval_beg, - sources, coverage_cut); - ///fprintf(stderr, "ploid_coverage: %lu, cov_threshold: %lu\n", ploid_coverage, cov_threshold); + ploid_coverage = get_pair_purge_coverage(xReads, *r_x_pos_beg, *r_x_pos_end, + yReads, *r_y_pos_beg, *r_y_pos_end, hap_can->rev, read_g, cov); if(cov_threshold > 0 && ploid_coverage >= cov_threshold) return NON_PLOID; + get_pair_hap_similarity_by_base(xReads, read_g, yUid, reverse_sources, ruIndex, + *r_x_pos_beg, *r_x_pos_end, &xLeftMatch, &xLeftTotal); + if(xLeftMatch == 0 || xLeftTotal == 0 || xLeftMatch <= xLeftTotal*Hap_rate) + { + return NON_PLOID; + } + return PLOID; } return NON_PLOID; @@ -1957,277 +2582,6 @@ void print_hap_paf(ma_ug_t *ug, hap_overlaps* ovlp) ovlp->y_beg_pos, ovlp->y_beg_id, ovlp->y_end_pos, ovlp->y_end_id, ovlp->type, (uint32_t)ovlp->weight); } -void hap_alignment(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, -R_to_U* ruIndex, ma_sub_t *coverage_cut, uint64_t* position_index, uint64_t* vote_counting, -uint8_t* visit, kvec_t_u64_warp* u_vecs, kvec_asg_arc_t_offset* u_buffer, kvec_hap_candidates* u_can, -uint32_t Input_uId, float Hap_rate, int max_hang, int min_ovlp, float chain_rate, -hap_overlaps_list* all_ovlp) -{ - ma_utg_t *xReads = NULL, *yReads = NULL; - ma_hit_t_alloc *xR = NULL; - ma_hit_t *h = NULL; - ma_sub_t *sq = NULL, *st = NULL; - asg_t* nsg = ug->g; - uint32_t i, j, v, rId, k, is_Unitig, Hap_uId, xUid, yUid, seedOcc, xPos, yPos, is_update; - uint64_t tmp; - long long cur_offset, new_offset, interval_len; - long long r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end; - int32_t r; - asg_arc_t t; - asg_arc_t_offset t_offset; - hap_candidates hap_can; - memset(&hap_can, 0, sizeof(hap_candidates)); - hap_overlaps hap_align; - xUid = Input_uId; - if(nsg->seq[xUid].del || nsg->seq[xUid].c == ALTER_LABLE) return; - memset(vote_counting, 0, sizeof(uint64_t)*nsg->n_seq); - memset(visit, 0, nsg->n_seq); - u_vecs->a.n = 0; - u_can->a.n = 0; - - xReads = &(ug->u.a[xUid]); - for (i = 0; i < xReads->n; i++) - { - xR = &(reverse_sources[xReads->a[i]>>33]); - - for (k = 0; k < xR->length; k++) - { - rId = Get_tn(xR->buffer[k]); - - if(read_g->seq[rId].del == 1) - { - ///get the id of read that contains it - get_R_to_U(ruIndex, rId, &rId, &is_Unitig); - if(rId == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[rId].del == 1) continue; - } - - ///there are two cases: - ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id - ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1 - get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig); - if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; - ///here rId is the id of the read coming from the different haplotype - ///Hap_cId is the id of the corresponding contig (note here is the contig, instead of untig) - if(visit[Hap_uId]!=0) continue; - visit[Hap_uId] = 1; - if(vote_counting[Hap_uId] < UINT64_MAX) vote_counting[Hap_uId]++; - } - - clean_visit_flag(visit, read_g, ruIndex, nsg->n_seq, xR); - } - - - - u_vecs->a.n = 0; - for (i = 0; i < nsg->n_seq; i++) - { - if(i == xUid) continue; - if(vote_counting[i] == 0) continue; - tmp = vote_counting[i]; tmp = tmp << 32; tmp = tmp | (uint64_t)i; - kv_push(uint64_t, u_vecs->a, tmp); - } - - if(u_vecs->a.n == 0) return; - sort_kvec_t_u64_warp(u_vecs, 1); - - - ///scan each candidate unitig - for (i = 0; i < u_vecs->a.n; i++) - { - yUid = (uint32_t)u_vecs->a.a[i]; - seedOcc = u_vecs->a.a[i]>>32; - xReads = &(ug->u.a[xUid]); - yReads = &(ug->u.a[yUid]); - u_buffer->a.n = 0; - - for (k = 0; k < xReads->n; k++) - { - xR = &(reverse_sources[xReads->a[k]>>33]); - for (j = 0; j < xR->length; j++) - { - h = &(xR->buffer[j]); - sq = &(coverage_cut[Get_qn(*h)]); - st = &(coverage_cut[Get_tn(*h)]); - if(st->del || read_g->seq[Get_tn(*h)].del) continue; - - r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, - asm_opt.max_hang_rate, min_ovlp, &t); - ///if it is a contained overlap, skip - if(r < 0) continue; - - rId = t.v>>1; - if(read_g->seq[rId].del == 1) continue; - ///there are two cases: - ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id - ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1 - get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig); - if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; - if(Hap_uId != yUid) continue; - - v = xReads->a[k]>>32; - get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig); - if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; - if(Hap_uId != xUid) continue; - if((uint32_t)(position_index[v>>1]) != k) continue; - - - if((prefilter((uint32_t)(position_index[v>>1]), (uint32_t)(position_index[rId]), - xReads->n, yReads->n, 0, Hap_rate, seedOcc)==NON_PLOID) && - (prefilter((uint32_t)(position_index[v>>1]), (uint32_t)(position_index[rId]), - xReads->n, yReads->n, 1, Hap_rate, seedOcc)==NON_PLOID)) - { - continue; - } - - t_offset.Off = get_xy_pos(read_g, &t, v, (yReads->a[(uint32_t)(position_index[rId])])>>32, - xReads->len, yReads->len, position_index, &(t.el)); - if(((t_offset.Off>>32) == (uint32_t)-1) || (((uint32_t)t_offset.Off) == (uint32_t)-1)) continue; - - t_offset.x = t; - t_offset.weight = 1; - - kv_push(asg_arc_t_offset, u_buffer->a, t_offset); - } - - deduplicate_edge(u_buffer); - } - - if(u_buffer->a.n == 0) continue; - - qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment); - k = 0; - u_can->a.n = 0; - while (k < u_buffer->a.n) - { - hap_can.rev = u_buffer->a.a[k].x.el; - hap_can.index_beg = k; - hap_can.index_end = k; - hap_can.weight = u_buffer->a.a[k].weight; - hap_can.x_beg_pos = hap_can.x_end_pos = (uint32_t)(u_buffer->a.a[k].Off>>32); - hap_can.y_beg_pos = hap_can.y_end_pos = (uint32_t)(u_buffer->a.a[k].Off); - cur_offset = Cal_Off(u_buffer->a.a[k].Off); - interval_len = get_hap_overlapLen(hap_can.x_beg_pos, hap_can.x_end_pos, xReads->len, - hap_can.y_beg_pos, hap_can.y_end_pos, yReads->len, NULL, NULL, NULL, NULL); - - - k++; - while (k < u_buffer->a.n) - { - new_offset = Cal_Off(u_buffer->a.a[k].Off); - if(u_buffer->a.a[k].x.el != hap_can.rev) break; - if((new_offset - cur_offset)>(interval_len*chain_rate)) break; - - - hap_can.index_end = k; - hap_can.weight += u_buffer->a.a[k].weight; - - is_update = 0; - xPos = (uint32_t)(u_buffer->a.a[k].Off>>32); - yPos = (uint32_t)(u_buffer->a.a[k].Off); - if(xPos < hap_can.x_beg_pos) - { - hap_can.x_beg_pos = xPos; - is_update = 1; - } - - if(xPos > hap_can.x_end_pos) - { - hap_can.x_end_pos = xPos; - is_update = 1; - } - - if(yPos < hap_can.y_beg_pos) - { - hap_can.y_beg_pos = yPos; - is_update = 1; - } - - if(yPos > hap_can.y_end_pos) - { - hap_can.y_end_pos = yPos; - is_update = 1; - } - - if(new_offset == cur_offset) is_update = 0; - - if(is_update) - { - interval_len = get_hap_overlapLen(hap_can.x_beg_pos, hap_can.x_end_pos, xReads->len, - hap_can.y_beg_pos, hap_can.y_end_pos, yReads->len, NULL, NULL, NULL, NULL); - } - - k++; - } - - kv_push(hap_candidates, u_can->a, hap_can); - } - - if(u_can->a.n == 0) continue; - - qsort(u_can->a.a, u_can->a.n, sizeof(hap_candidates), cmp_hap_candidates); - - Get_match(hap_can) = Get_total(hap_can) = 0; - memset(&hap_align, 0, sizeof(hap_overlaps)); - - for (k = 0; k < u_can->a.n; k++) - { - is_update = 0; - if(u_can->a.a[k].weight < Get_match(hap_can)*Hap_rate) continue; - - if(calculate_pair_hap_similarity(u_buffer, &(u_can->a.a[k]), position_index, xUid, yUid, - xReads, yReads, reverse_sources, read_g, ruIndex, coverage_cut, Hap_rate, max_hang, - min_ovlp, &r_x_pos_beg, &r_x_pos_end, &r_y_pos_beg, &r_y_pos_end)!=PLOID) - { - continue; - } - - if(Get_match(hap_can) < Get_match(u_can->a.a[k])) - { - is_update = 1; - } - else if(Get_match(hap_can) == Get_match(u_can->a.a[k]) && - Get_total(hap_can) > Get_total(u_can->a.a[k])) - { - is_update = 1; - } - - if(is_update) - { - hap_can = u_can->a.a[k]; - hap_align.rev = Get_rev(hap_can); - hap_align.type = Get_type(hap_can); - hap_align.x_beg_id = Get_x_beg(hap_can); - hap_align.x_end_id = Get_x_end(hap_can) + 1; - hap_align.y_beg_id = Get_y_beg(hap_can); - hap_align.y_end_id = Get_y_end(hap_can) + 1; - hap_align.weight = Get_match(hap_can); - hap_align.x_beg_pos = r_x_pos_beg; - hap_align.x_end_pos = r_x_pos_end + 1; - if(hap_align.rev == 0) - { - hap_align.y_beg_pos = r_y_pos_beg; - hap_align.y_end_pos = r_y_pos_end + 1; - } - else - { - hap_align.y_beg_pos = yReads->len - r_y_pos_end - 1; - hap_align.y_end_pos = yReads->len - r_y_pos_beg - 1 + 1; - } - hap_align.xUid = xUid; - hap_align.yUid = yUid; - hap_align.status = SELF_EXIST; - } - } - - if(Get_match(hap_can) == 0 || Get_total(hap_can) == 0) continue; - - kv_push(hap_overlaps, all_ovlp->x[hap_align.xUid].a, hap_align); - } - -} - - - inline long long get_max_index(asg_arc_t_offset* x, int32_t* Scores, uint8_t* Flag, long long n, long long x_readLen, long long y_readLen) { @@ -2474,7 +2828,7 @@ long long y_readLen) if(is_merge == 0 && (Get_xOff(u_buffer->a.a[m-1].Off)==Get_xOff(u_buffer->a.a[i].Off))) { if((Get_yOff(u_buffer->a.a[i].Off)-(Get_yOff(u_buffer->a.a[m-1].Off))) == - (i-anchor_i)) + (i-anchor_i))///not sure why, does it use for tolerate indels in overlaps? { is_merge = 1; } @@ -2503,300 +2857,15 @@ long long y_readLen) begIndex_vec, flag_vec, band_width_threshold, max_skip, x_readLen, y_readLen, u_can); } - -///static void hap_alignment_worker(void *_data, long eid, int tid) -void hap_alignment_worker(void *_data, long eid, int tid) +int filter_secondary_chain(long long max_score, long long cur_score, double rate) { - 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; - ma_hit_t_alloc* reverse_sources = hap_buf->reverse_sources; - R_to_U* ruIndex = hap_buf->ruIndex; - ma_sub_t *coverage_cut = hap_buf->coverage_cut; - uint64_t* position_index = hap_buf->position_index; - float Hap_rate = hap_buf->Hap_rate; - int max_hang = hap_buf->max_hang; - int min_ovlp = hap_buf->min_ovlp; - float chain_rate = hap_buf->chain_rate; - hap_overlaps_list* all_ovlp = hap_buf->all_ovlp; - uint32_t Input_uId = eid; - uint64_t* vote_counting = hap_buf->buf[tid].vote_counting; - uint8_t* visit = hap_buf->buf[tid].visit; - kvec_t_u64_warp* u_vecs = &(hap_buf->buf[tid].u_vecs); - kvec_asg_arc_t_offset* u_buffer = &(hap_buf->buf[tid].u_buffer); - kvec_hap_candidates* u_can = &(hap_buf->buf[tid].u_can); - - - ma_utg_t *xReads = NULL, *yReads = NULL; - ma_hit_t_alloc *xR = NULL; - ma_hit_t *h = NULL; - ma_sub_t *sq = NULL, *st = NULL; - asg_t* nsg = ug->g; - uint32_t i, j, v, rId, k, is_Unitig, Hap_uId, xUid, yUid, seedOcc, xPos, yPos, is_update; - uint64_t tmp; - long long cur_offset, new_offset, interval_len; - long long r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end; - int32_t r; - asg_arc_t t; - asg_arc_t_offset t_offset; - hap_candidates hap_can; - memset(&hap_can, 0, sizeof(hap_candidates)); - hap_overlaps hap_align; - xUid = Input_uId; - if(nsg->seq[xUid].del || nsg->seq[xUid].c == ALTER_LABLE) return; - memset(vote_counting, 0, sizeof(uint64_t)*nsg->n_seq); - memset(visit, 0, nsg->n_seq); - u_vecs->a.n = 0; - u_can->a.n = 0; - - xReads = &(ug->u.a[xUid]); - for (i = 0; i < xReads->n; i++) - { - xR = &(reverse_sources[xReads->a[i]>>33]); - - for (k = 0; k < xR->length; k++) - { - rId = Get_tn(xR->buffer[k]); - - if(read_g->seq[rId].del == 1) - { - ///get the id of read that contains it - get_R_to_U(ruIndex, rId, &rId, &is_Unitig); - if(rId == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[rId].del == 1) continue; - } - - ///there are two cases: - ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id - ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1 - get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig); - if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; - ///here rId is the id of the read coming from the different haplotype - ///Hap_cId is the id of the corresponding contig (note here is the contig, instead of untig) - if(visit[Hap_uId]!=0) continue; - visit[Hap_uId] = 1; - if(vote_counting[Hap_uId] < UINT64_MAX) vote_counting[Hap_uId]++; - } - - clean_visit_flag(visit, read_g, ruIndex, nsg->n_seq, xR); - } - - - - u_vecs->a.n = 0; - for (i = 0; i < nsg->n_seq; i++) - { - if(i == xUid) continue; - if(vote_counting[i] == 0) continue; - tmp = vote_counting[i]; tmp = tmp << 32; tmp = tmp | (uint64_t)i; - kv_push(uint64_t, u_vecs->a, tmp); - } - - if(u_vecs->a.n == 0) return; - sort_kvec_t_u64_warp(u_vecs, 1); - - - ///scan each candidate unitig - for (i = 0; i < u_vecs->a.n; i++) - { - yUid = (uint32_t)u_vecs->a.a[i]; - seedOcc = u_vecs->a.a[i]>>32; - xReads = &(ug->u.a[xUid]); - yReads = &(ug->u.a[yUid]); - u_buffer->a.n = 0; - - for (k = 0; k < xReads->n; k++) - { - xR = &(reverse_sources[xReads->a[k]>>33]); - for (j = 0; j < xR->length; j++) - { - h = &(xR->buffer[j]); - sq = &(coverage_cut[Get_qn(*h)]); - st = &(coverage_cut[Get_tn(*h)]); - if(st->del || read_g->seq[Get_tn(*h)].del) continue; - - r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, - asm_opt.max_hang_rate, min_ovlp, &t); - ///if it is a contained overlap, skip - if(r < 0) continue; - - rId = t.v>>1; - if(read_g->seq[rId].del == 1) continue; - ///there are two cases: - ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id - ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1 - get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig); - if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; - if(Hap_uId != yUid) continue; - - v = xReads->a[k]>>32; - get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig); - if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; - if(Hap_uId != xUid) continue; - if((uint32_t)(position_index[v>>1]) != k) continue; - - - if((prefilter((uint32_t)(position_index[v>>1]), (uint32_t)(position_index[rId]), - xReads->n, yReads->n, 0, Hap_rate, seedOcc)==NON_PLOID) && - (prefilter((uint32_t)(position_index[v>>1]), (uint32_t)(position_index[rId]), - xReads->n, yReads->n, 1, Hap_rate, seedOcc)==NON_PLOID)) - { - continue; - } - - t_offset.Off = get_xy_pos(read_g, &t, v, (yReads->a[(uint32_t)(position_index[rId])])>>32, - xReads->len, yReads->len, position_index, &(t.el)); - if(((t_offset.Off>>32) == (uint32_t)-1) || (((uint32_t)t_offset.Off) == (uint32_t)-1)) continue; - - t_offset.x = t; - t_offset.weight = 1; - - kv_push(asg_arc_t_offset, u_buffer->a, t_offset); - } - - deduplicate_edge(u_buffer); - } - - if(u_buffer->a.n == 0) continue; - - // if(debug_enable) - // { - // print_debug_unitig(xReads, position_index, "xReads"); - // print_debug_unitig(yReads, position_index, "yReads"); - // } - - qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment); - k = 0; - u_can->a.n = 0; - while (k < u_buffer->a.n) - { - hap_can.rev = u_buffer->a.a[k].x.el; - hap_can.index_beg = k; - hap_can.index_end = k; - hap_can.weight = u_buffer->a.a[k].weight; - hap_can.x_beg_pos = hap_can.x_end_pos = (uint32_t)(u_buffer->a.a[k].Off>>32); - hap_can.y_beg_pos = hap_can.y_end_pos = (uint32_t)(u_buffer->a.a[k].Off); - cur_offset = Cal_Off(u_buffer->a.a[k].Off); - interval_len = get_hap_overlapLen(hap_can.x_beg_pos, hap_can.x_end_pos, xReads->len, - hap_can.y_beg_pos, hap_can.y_end_pos, yReads->len, NULL, NULL, NULL, NULL); - - - k++; - while (k < u_buffer->a.n) - { - new_offset = Cal_Off(u_buffer->a.a[k].Off); - if(u_buffer->a.a[k].x.el != hap_can.rev) break; - if((new_offset - cur_offset)>(interval_len*chain_rate)) break; - - - hap_can.index_end = k; - hap_can.weight += u_buffer->a.a[k].weight; - - is_update = 0; - xPos = (uint32_t)(u_buffer->a.a[k].Off>>32); - yPos = (uint32_t)(u_buffer->a.a[k].Off); - if(xPos < hap_can.x_beg_pos) - { - hap_can.x_beg_pos = xPos; - is_update = 1; - } - - if(xPos > hap_can.x_end_pos) - { - hap_can.x_end_pos = xPos; - is_update = 1; - } - - if(yPos < hap_can.y_beg_pos) - { - hap_can.y_beg_pos = yPos; - is_update = 1; - } - - if(yPos > hap_can.y_end_pos) - { - hap_can.y_end_pos = yPos; - is_update = 1; - } - - if(new_offset == cur_offset) is_update = 0; - - if(is_update) - { - interval_len = get_hap_overlapLen(hap_can.x_beg_pos, hap_can.x_end_pos, xReads->len, - hap_can.y_beg_pos, hap_can.y_end_pos, yReads->len, NULL, NULL, NULL, NULL); - } - - k++; - } - - kv_push(hap_candidates, u_can->a, hap_can); - } - - if(u_can->a.n == 0) continue; - - qsort(u_can->a.a, u_can->a.n, sizeof(hap_candidates), cmp_hap_candidates); - - Get_match(hap_can) = Get_total(hap_can) = 0; - memset(&hap_align, 0, sizeof(hap_overlaps)); - - for (k = 0; k < u_can->a.n; k++) - { - is_update = 0; - if(u_can->a.a[k].weight < Get_match(hap_can)*Hap_rate) continue; - - if(calculate_pair_hap_similarity(u_buffer, &(u_can->a.a[k]), position_index, xUid, yUid, - xReads, yReads, reverse_sources, read_g, ruIndex, coverage_cut, Hap_rate, max_hang, - min_ovlp, &r_x_pos_beg, &r_x_pos_end, &r_y_pos_beg, &r_y_pos_end)!=PLOID) - { - continue; - } - - if(Get_match(hap_can) < Get_match(u_can->a.a[k])) - { - is_update = 1; - } - else if(Get_match(hap_can) == Get_match(u_can->a.a[k]) && - Get_total(hap_can) > Get_total(u_can->a.a[k])) - { - is_update = 1; - } - - if(is_update) - { - hap_can = u_can->a.a[k]; - hap_align.rev = Get_rev(hap_can); - hap_align.type = Get_type(hap_can); - hap_align.x_beg_id = Get_x_beg(hap_can); - hap_align.x_end_id = Get_x_end(hap_can) + 1; - hap_align.y_beg_id = Get_y_beg(hap_can); - hap_align.y_end_id = Get_y_end(hap_can) + 1; - hap_align.weight = Get_match(hap_can); - hap_align.x_beg_pos = r_x_pos_beg; - hap_align.x_end_pos = r_x_pos_end + 1; - if(hap_align.rev == 0) - { - hap_align.y_beg_pos = r_y_pos_beg; - hap_align.y_end_pos = r_y_pos_end + 1; - } - else - { - hap_align.y_beg_pos = yReads->len - r_y_pos_end - 1; - hap_align.y_end_pos = yReads->len - r_y_pos_beg - 1 + 1; - } - hap_align.xUid = xUid; - hap_align.yUid = yUid; - hap_align.status = SELF_EXIST; - } - } - - if(Get_match(hap_can) == 0 || Get_total(hap_can) == 0) continue; - - kv_push(hap_overlaps, all_ovlp->x[hap_align.xUid].a, hap_align); - } - + if(cur_score >= max_score) return 1; + long long diff = max_score - cur_score; + if(max_score < 0) max_score *= -1; + if(diff >= max_score*(1-rate)) return 0; + return 1; } - static void hap_alignment_advance_worker(void *_data, long eid, int tid) { hap_alignment_struct_pip* hap_buf = (hap_alignment_struct_pip*)_data; @@ -2823,21 +2892,21 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) kvec_t_i32_warp* begIndex_vec = &(hap_buf->buf[tid].u_buffer_beg); kvec_t_u8_warp* flag_vec = &(hap_buf->buf[tid].u_buffer_flag); uint64_t cov_threshold = hap_buf->cov_threshold; + hap_cov_t *cov = hap_buf->cov; if(hap_buf->cov_threshold < 0) cov_threshold = (uint64_t)-1; - ma_utg_t *xReads = NULL, *yReads = NULL; ma_hit_t_alloc *xR = NULL; ma_hit_t *h = NULL; ma_sub_t *sq = NULL, *st = NULL; asg_t* nsg = ug->g; - uint32_t i, j, v, rId, k, is_Unitig, Hap_uId, xUid, yUid, seedOcc, is_update; - uint64_t tmp; - long long r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end; + uint32_t i, j, v, rId, k, is_Unitig, Hap_uId, xUid, yUid, seedOcc; + uint64_t tmp, max_weight, m; + long long r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end, max_score; int32_t r; asg_arc_t t; asg_arc_t_offset t_offset; - hap_candidates hap_can; hap_overlaps hap_align; + hap_overlaps *hap_align_x = NULL; xUid = Input_uId; if(nsg->seq[xUid].del || nsg->seq[xUid].c == ALTER_LABLE) return; memset(vote_counting, 0, sizeof(uint64_t)*nsg->n_seq); @@ -2868,7 +2937,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; ///here rId is the id of the read coming from the different haplotype ///Hap_cId is the id of the corresponding contig (note here is the contig, instead of untig) - if(visit[Hap_uId]!=0) continue; + if(visit[Hap_uId]!=0) continue; ///one read only has one vote for one hap unitig visit[Hap_uId] = 1; if(vote_counting[Hap_uId] < UINT64_MAX) vote_counting[Hap_uId]++; } @@ -2903,6 +2972,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) for (k = 0; k < xReads->n; k++) { xR = &(reverse_sources[xReads->a[k]>>33]); + for (j = 0; j < xR->length; j++) { h = &(xR->buffer[j]); @@ -2930,8 +3000,8 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) if(Hap_uId != xUid) continue; if((uint32_t)(position_index[v>>1]) != k) continue; - - if((prefilter((uint32_t)(position_index[v>>1]), (uint32_t)(position_index[rId]), + if(asm_opt.purge_level_primary <= 2 && + (prefilter((uint32_t)(position_index[v>>1]), (uint32_t)(position_index[rId]), xReads->n, yReads->n, 0, Hap_rate, seedOcc)==NON_PLOID) && (prefilter((uint32_t)(position_index[v>>1]), (uint32_t)(position_index[rId]), xReads->n, yReads->n, 1, Hap_rate, seedOcc)==NON_PLOID)) @@ -2962,64 +3032,88 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) qsort(u_can->a.a, u_can->a.n, sizeof(hap_candidates), cmp_hap_candidates); - Get_match(hap_can) = Get_total(hap_can) = 0; - memset(&hap_align, 0, sizeof(hap_overlaps)); - + memset(&hap_align, 0, sizeof(hap_overlaps)); + m = all_ovlp->x[xUid].a.n; + max_weight = 0; max_score = 0; for (k = 0; k < u_can->a.n; k++) { - is_update = 0; - if(u_can->a.a[k].weight < Get_match(hap_can)*Hap_rate) continue; - + if(u_can->a.a[k].weight < max_weight*0.33) continue; 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) + xReads, yReads, sources, reverse_sources, read_g, ruIndex, coverage_cut, Hap_rate, + (asm_opt.purge_level_primary<=2? 0:1), max_hang, min_ovlp, cov_threshold, u_buffer, + score_vc, prevIndex_vec, cov, &r_x_pos_beg, &r_x_pos_end, &r_y_pos_beg, &r_y_pos_end)!=PLOID) { continue; } - - - if(Get_match(hap_can) < Get_match(u_can->a.a[k])) + ///max_weight == 0 means the first matched chain + if(max_weight == 0 || max_score < u_can->a.a[k].score) max_score = u_can->a.a[k].score; + if(max_weight < u_can->a.a[k].weight) max_weight = u_can->a.a[k].weight; + ///if one is positive and another one is negative, it is wrong + if(!filter_secondary_chain(max_score, u_can->a.a[k].score, CHAIN_FILTER_RATE)) continue; + + hap_align.rev = Get_rev(u_can->a.a[k]); + hap_align.type = Get_type(u_can->a.a[k]); + hap_align.x_beg_id = Get_x_beg(u_can->a.a[k]); + hap_align.x_end_id = Get_x_end(u_can->a.a[k]) + 1; + hap_align.y_beg_id = Get_y_beg(u_can->a.a[k]); + hap_align.y_end_id = Get_y_end(u_can->a.a[k]) + 1; + hap_align.weight = Get_match(u_can->a.a[k]); + hap_align.score = u_can->a.a[k].score; + hap_align.x_beg_pos = r_x_pos_beg; + hap_align.x_end_pos = r_x_pos_end + 1; + if(hap_align.rev == 0) { - is_update = 1; + hap_align.y_beg_pos = r_y_pos_beg; + hap_align.y_end_pos = r_y_pos_end + 1; } - else if(Get_match(hap_can) == Get_match(u_can->a.a[k]) && - Get_total(hap_can) > Get_total(u_can->a.a[k])) + else { - is_update = 1; + hap_align.y_beg_pos = yReads->len - r_y_pos_end - 1; + hap_align.y_end_pos = yReads->len - r_y_pos_beg - 1 + 1; } - - if(is_update) + hap_align.xUid = xUid; + hap_align.yUid = yUid; + hap_align.status = SELF_EXIST; + kv_push(hap_overlaps, all_ovlp->x[hap_align.xUid].a, hap_align); + } + /** + ///chains with same xUid && yUid + for (k = m; k < all_ovlp->x[xUid].a.n; k++) + { + if(!filter_secondary_chain(max_score, + all_ovlp->x[xUid].a.a[k].score, CHAIN_FILTER_RATE)) { - hap_can = u_can->a.a[k]; - hap_align.rev = Get_rev(hap_can); - hap_align.type = Get_type(hap_can); - hap_align.x_beg_id = Get_x_beg(hap_can); - hap_align.x_end_id = Get_x_end(hap_can) + 1; - hap_align.y_beg_id = Get_y_beg(hap_can); - hap_align.y_end_id = Get_y_end(hap_can) + 1; - hap_align.weight = Get_match(hap_can); - hap_align.x_beg_pos = r_x_pos_beg; - hap_align.x_end_pos = r_x_pos_end + 1; - if(hap_align.rev == 0) + continue; + } + all_ovlp->x[xUid].a.a[m] = all_ovlp->x[xUid].a.a[k]; + m++; + } + all_ovlp->x[xUid].a.n = m; + **/ + + hap_align_x = NULL; + for (k = m; k < all_ovlp->x[xUid].a.n; k++) + { + if(all_ovlp->x[xUid].a.a[k].score != max_score) continue; + if(hap_align_x == NULL || all_ovlp->x[xUid].a.a[k].weight > hap_align_x->weight) + { + hap_align_x = &(all_ovlp->x[xUid].a.a[k]); + } + else if(all_ovlp->x[xUid].a.a[k].weight == hap_align_x->weight) + { + if((all_ovlp->x[xUid].a.a[k].x_end_pos + - all_ovlp->x[xUid].a.a[k].x_beg_pos) < + (hap_align_x->x_end_pos - hap_align_x->x_beg_pos)) { - hap_align.y_beg_pos = r_y_pos_beg; - hap_align.y_end_pos = r_y_pos_end + 1; + hap_align_x = &(all_ovlp->x[xUid].a.a[k]); } - else - { - hap_align.y_beg_pos = yReads->len - r_y_pos_end - 1; - hap_align.y_end_pos = yReads->len - r_y_pos_beg - 1 + 1; - } - hap_align.xUid = xUid; - hap_align.yUid = yUid; - hap_align.status = SELF_EXIST; } } - - if(Get_match(hap_can) == 0 || Get_total(hap_can) == 0) continue; - - kv_push(hap_overlaps, all_ovlp->x[hap_align.xUid].a, hap_align); + if(hap_align_x) + { + all_ovlp->x[xUid].a.a[m] = (*hap_align_x); + all_ovlp->x[xUid].a.n = m + 1; + } } } @@ -3137,8 +3231,9 @@ ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex) y = &(all_ovlp->x[tn].a.a[index]); if(x->rev == y->rev && types[x->type]==y->type) continue; ///if(x->weight >= y->weight) - if((calculate_bi_weight(x, ug, read_g, reverse_sources, ruIndex)) >= - (calculate_bi_weight(y, ug, read_g, reverse_sources, ruIndex))) + // if((calculate_bi_weight(x, ug, read_g, reverse_sources, ruIndex)) >= + // (calculate_bi_weight(y, ug, read_g, reverse_sources, ruIndex))) + if(x->score >= y->score) { kv_push(hap_overlaps, back_all_ovlp->x[tn].a, (*y)); set_reverse_hap_overlap(y, x, types); @@ -3160,6 +3255,7 @@ ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex) void filter_hap_overlaps_by_length(hap_overlaps_list* all_ovlp, uint32_t minLen) { + if(minLen == 0) return; hap_overlaps *x = NULL; uint32_t v, i, m, uId; @@ -3269,6 +3365,137 @@ void print_purge_gfa(ma_ug_t *ug, asg_t *purge_g) } +long long decode_score(uint32_t h_bits, uint32_t l_bits) +{ + uint64_t x; + x = h_bits; x <<= 32; x += l_bits; + long long score = ((uint64_t)((uint64_t)x<<1)>>1); + if((x>>63) == 0) score *= -1; + return score; +} + +void encode_score(long long i_s, uint32_t *h_bits, uint32_t *l_bits) +{ + uint64_t score = (i_s >= 0? (i_s) : (i_s*(-1))); + if(i_s >= 0) score += (((uint64_t)1)<<63); + (*l_bits) = (uint32_t)score; (*h_bits)= (score>>32); +} + +uint64_t asg_bub_pop1_purge_graph(asg_t *g, uint32_t v0, int max_dist, buf_t *b) +{ + uint32_t i, n_pending = 0, n_tips, tip_end; + uint64_t n_pop = 0; + ///if this node has been deleted + if (g->seq[v0>>1].del || g->seq[v0>>1].c == ALTER_LABLE) return 0; // already deleted + ///if ((uint32_t)g->idx[v0] < 2) return 0; // no bubbles + if(get_real_length(g, v0, NULL)<2) return 0; + ///S saves nodes with all incoming edges visited + b->S.n = b->T.n = b->b.n = b->e.n = 0; + ///for each node, b->a saves all related information + b->a[v0].c = b->a[v0].d = b->a[v0].m = b->a[v0].nc = b->a[v0].np = 0; + ///b->S is the nodes with all incoming edges visited + kv_push(uint32_t, b->S, v0); + n_tips = 0; + tip_end = (uint32_t)-1; + + do { + ///v is a node that all incoming edges have been visited + ///d is the distance from v0 to v + uint32_t v = kv_pop(b->S), d = b->a[v].d; + uint32_t nv = asg_arc_n(g, v); + asg_arc_t *av = asg_arc_a(g, v); + long long t_s = decode_score(b->a[v].c, b->a[v].m), c_s; + ///why we have this assert? + ///assert(nv > 0); + ///all out-edges of v + for (i = 0; i < nv; ++i) { // loop through v's neighbors + uint32_t w = av[i].v; // v->w with length l + binfo_t *t = &b->a[w]; + ///that means there is a circle, directly terminate the whole bubble poping + ///if (w == v0) goto pop_reset; + if ((w>>1) == (v0>>1)) goto pop_reset; + ///if this edge has been deleted + if (av[i].del) continue; + c_s = decode_score((uint32_t)av[i].ul, av[i].ol); + ///push the edge + ///high 32-bit of g->idx[v] is the start point of v's edges + //so here is the point of this specfic edge + kv_push(uint32_t, b->e, (g->idx[v]>>32) + i); + ///find a too far path? directly terminate the whole bubble poping + if (d + 1 > (uint32_t)max_dist) break; // too far + + ///if this node + if (t->s == 0) { // this vertex has never been visited + kv_push(uint32_t, b->b, w); // save it for revert + ///t->p is the parent node of + ///t->s = 1 means w has been visited + ///d is len(v0->v), l is len(v->w), so t->d is len(v0->w) + t->p = v, t->s = 1, t->d = d + 1; + encode_score(t_s + c_s, &(t->c), &(t->m)); + ///incoming edges of w + ///t->r = count_out(g, w^1); + t->r = get_real_length(g, w^1, NULL); + ++n_pending; + } else { // visited before + if((t_s + c_s)> decode_score(t->c, t->m)) + { + t->p = v; + encode_score(t_s + c_s, &(t->c), &(t->m)); + } + ///it is the shortest edge + if (d + 1 < t->d) t->d = d + 1; // update dist + } + ///assert(t->r > 0); + //if all incoming edges of w have visited + //push it to b->S + if (--(t->r) == 0) { + uint32_t x = get_real_length(g, w, NULL); + /****************************may have bugs for bubble********************************/ + if(x > 0) + { + kv_push(uint32_t, b->S, w); + } + else + { + ///at most one tip + if(n_tips != 0) goto pop_reset; + n_tips++; + tip_end = w; + } + /****************************may have bugs for bubble********************************/ + --n_pending; + } + } + //if found a tip + /****************************may have bugs for bubble********************************/ + if(n_tips == 1) + { + if(tip_end != (uint32_t)-1 && n_pending == 0 && b->S.n == 0) + { + kv_push(uint32_t, b->S, tip_end); + break; + } + else + { + goto pop_reset; + } + } + /****************************may have bugs for bubble********************************/ + ///if i < nv, that means (d + l > max_dist) + if (i < nv || b->S.n == 0) goto pop_reset; + } while (b->S.n > 1 || n_pending); + + asg_bub_backtrack_primary(g, v0, b); + + n_pop = 1; +pop_reset: + for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices + binfo_t *t = &b->a[b->b.a[i]]; + t->s = t->c = t->d = t->m = t->nc = t->np = 0; + } + return n_pop; +} + // pop bubbles @@ -3291,7 +3518,7 @@ int asg_pop_bubble_purge_graph(asg_t *purge_g, int max_dist) for (i = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs if (!av[i].del) ++n_arc; if (n_arc > 1) - n_pop += asg_bub_pop1_primary_trio(purge_g, NULL, v, max_dist, &b, (uint32_t)-1, DROP, 1, NULL, NULL); + n_pop += asg_bub_pop1_purge_graph(purge_g, v, purge_g->n_seq, &b); } free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); if (n_pop) asg_cleanup(purge_g); @@ -3320,9 +3547,98 @@ int min_ovlp, asg_arc_t* t) h.bl = h.el = h.ml = h.no_l_indel = 0; r = ma_hit2arc(&h, qLen, tLen, max_hang, max_hang_rate, min_ovlp, t); + if(r < 0) return r; + uint64_t score = (hap->score >= 0? (hap->score) : (hap->score*(-1))); + if(hap->score >= 0) score += (((uint64_t)1)<<63); + t->ol = (uint32_t)score; + t->ul >>= 32; t->ul <<= 32; t->ul |= (score>>32); return r; } +typedef struct { + uint64_t eid; + uint64_t score; +}e_score; + +typedef struct { + size_t n, m; + e_score* a; +}e_score_warp; + +#define e_score_key(a) ((a).score) +KRADIX_SORT_INIT(e_score, e_score, e_score_key, member_size(e_score, score)) + +int purge_g_arc_del_short_diploid_by_score(asg_t *g, float drop_ratio) +{ + e_score_warp b; + kv_init(b); + e_score *p = NULL; + + uint32_t v, n_vtx = g->n_seq * 2; + long long n_cut = 0; + + for (v = 0; v < n_vtx; ++v) + { + if(g->seq[v>>1].c == ALTER_LABLE || g->seq[v>>1].del) continue; + asg_arc_t *av = asg_arc_a(g, v); + uint32_t nv = asg_arc_n(g, v); + if (nv < 2) continue; + uint64_t i; + for (i = 0; i < nv; ++i) + { + kv_pushp(e_score, b, &p); + p->eid = av - g->arc + i; + p->score = (uint32_t)av[i].ul; + p->score <<= 32; + p->score |= av[i].ol; + } + } + + radix_sort_e_score(b.a, b.a + b.n); + + uint64_t k; + for (k = 0; k < b.n; k++) + { + asg_arc_t *a = &g->arc[b.a[k].eid]; + ///v is self id, w is the id of another end + uint32_t i, v = (a->ul)>>32; + uint32_t nv = asg_arc_n(g, v), kv; + long long ovlp_max = 0, ovlp; + asg_arc_t *av = NULL; + ///nv must be >= 2 + if (nv <= 1) continue; + av = asg_arc_a(g, v); + + ///calculate the longest edge for v and w + for (i = 0, kv = 0; i < nv; ++i) { + if (av[i].del) continue; + ovlp = decode_score((uint32_t)av[i].ul, av[i].ol); + if (kv == 0 || ovlp_max < ovlp) ovlp_max = ovlp; + ++kv; + } + + if (kv <= 1) continue; + ovlp = decode_score((uint32_t)a->ul, a->ol); + if (kv >= 2) + { + if(ovlp >= 0 && ovlp_max >= 0 && ovlp > ovlp_max * drop_ratio) continue; + } + + a->del = 1; + asg_arc_del(g, a->v^1, av->ul>>32^1, 1); + ++n_cut; + } + + kv_destroy(b); + if (n_cut) + { + asg_cleanup(g); + asg_symm(g); + } + + + return n_cut; +} void clean_purge_graph(asg_t *purge_g, int max_dist, float drop_ratio) @@ -3332,168 +3648,13 @@ void clean_purge_graph(asg_t *purge_g, int max_dist, float drop_ratio) { operation = 0; operation += asg_pop_bubble_purge_graph(purge_g, max_dist); - operation += unitig_arc_del_short_diploid_by_length(purge_g, drop_ratio); + operation += purge_g_arc_del_short_diploid_by_score(purge_g, drop_ratio); } - unitig_arc_del_short_diploid_by_length(purge_g, 1); + purge_g_arc_del_short_diploid_by_score(purge_g, 1); } -void get_node_boundary(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, -asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads, -uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, -long long yEndIndex, uint32_t dir, uint32_t rev, asg_arc_t* reture_t_f, asg_arc_t* reture_t_r) -{ - long long k, j, offset; - ma_hit_t_alloc *xR = NULL; - ma_hit_t *h = NULL; - ma_sub_t *sq = NULL, *st = NULL; - int r, index; - asg_arc_t t_f, t_r; - uint32_t rId, Hap_uId, is_Unitig, v, w, v_dir, w_dir, is_found = 0, oLen = 0; - reture_t_f->del = reture_t_r->del = 1; - if(dir == 1) - { - for (k = xEndIndex; k >= xBegIndex; k--) - { - xR = &(reverse_sources[xReads->a[k]>>33]); - is_found = 0; - for (j = 0; j < xR->length; j++) - { - h = &(xR->buffer[j]); - sq = &(coverage_cut[Get_qn(*h)]); - st = &(coverage_cut[Get_tn(*h)]); - if(st->del || read_g->seq[Get_tn(*h)].del) continue; - r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, - asm_opt.max_hang_rate, min_ovlp, &t_f); - ///if it is a contained overlap, skip - if(r < 0) continue; - - rId = t_f.v>>1; - if(read_g->seq[rId].del == 1) continue; - ///there are two cases: - ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id - ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1 - get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig); - if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; - if(Hap_uId != yUid) continue; - - v = xReads->a[k]>>32; - get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig); - if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; - if(Hap_uId != xUid) continue; - if((uint32_t)(position_index[v>>1]) != k) continue; - - w = (yReads->a[(uint32_t)(position_index[rId])])>>32; - v_dir = ((t_f.ul>>32)==v)?1:0; - w_dir = (t_f.v == w)?1:0; - - if(rev == 0 && v_dir != w_dir) continue; - if(rev == 1 && v_dir == w_dir) continue; - if(v_dir == 1) continue; - - /****************************may have bugs********************************/ - offset = (uint32_t)(position_index[rId]); - if(offset < yBegIndex || offset > yEndIndex) continue; - /****************************may have bugs********************************/ - - /************************get reverse edge*************************/ - index = get_specific_overlap(&(reverse_sources[Get_tn(*h)]), Get_tn(*h), Get_qn(*h)); - if(index == -1) continue; - h = &(reverse_sources[Get_tn(*h)].buffer[index]); - sq = &(coverage_cut[Get_qn(*h)]); - st = &(coverage_cut[Get_tn(*h)]); - if(st->del || read_g->seq[Get_tn(*h)].del) continue; - r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, - asm_opt.max_hang_rate, min_ovlp, &t_r); - if(r < 0) continue; - /************************get reverse edge*************************/ - - if(is_found == 0 || t_f.ol > oLen) - { - (*reture_t_f) = t_f; - (*reture_t_r) = t_r; - oLen = t_f.ol; - } - - is_found = 1; - } - if(is_found) return; - } - } - else - { - for (k = xBegIndex; k <= xEndIndex; k++) - { - xR = &(reverse_sources[xReads->a[k]>>33]); - is_found = 0; oLen = 0; - for (j = 0; j < xR->length; j++) - { - h = &(xR->buffer[j]); - sq = &(coverage_cut[Get_qn(*h)]); - st = &(coverage_cut[Get_tn(*h)]); - if(st->del || read_g->seq[Get_tn(*h)].del) continue; - - r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, - asm_opt.max_hang_rate, min_ovlp, &t_f); - ///if it is a contained overlap, skip - if(r < 0) continue; - - rId = t_f.v>>1; - if(read_g->seq[rId].del == 1) continue; - ///there are two cases: - ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id - ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1 - get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig); - if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; - if(Hap_uId != yUid) continue; - - v = xReads->a[k]>>32; - get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig); - if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; - if(Hap_uId != xUid) continue; - if((uint32_t)(position_index[v>>1]) != k) continue; - - w = (yReads->a[(uint32_t)(position_index[rId])])>>32; - - v_dir = ((t_f.ul>>32)==v)?1:0; - w_dir = (t_f.v == w)?1:0; - if(rev == 0 && v_dir != w_dir) continue; - if(rev == 1 && v_dir == w_dir) continue; - if(v_dir == 0) continue; - - /****************************may have bugs********************************/ - offset = (uint32_t)(position_index[rId]); - if(offset < yBegIndex || offset > yEndIndex) continue; - /****************************may have bugs********************************/ - - /************************get reverse edge*************************/ - index = get_specific_overlap(&(reverse_sources[Get_tn(*h)]), Get_tn(*h), Get_qn(*h)); - if(index == -1) continue; - h = &(reverse_sources[Get_tn(*h)].buffer[index]); - sq = &(coverage_cut[Get_qn(*h)]); - st = &(coverage_cut[Get_tn(*h)]); - if(st->del || read_g->seq[Get_tn(*h)].del) continue; - r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, - asm_opt.max_hang_rate, min_ovlp, &t_r); - if(r < 0) continue; - /************************get reverse edge*************************/ - - if(is_found == 0 || t_f.ol > oLen) - { - (*reture_t_f) = t_f; - (*reture_t_r) = t_r; - oLen = t_f.ol; - } - - is_found = 1; - } - if(is_found) return; - } - } - -} - void get_node_boundary_advance(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads, uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, @@ -3686,7 +3847,7 @@ uint32_t is_circle, uint64_t* rLen) if(k == edge->a.n) { - fprintf(stderr, "####ERROR1: i: %u, v>>1: %u, v&1: %u, w>>1: %u, w&1: %u\n", + fprintf(stderr, "####ERROR1-fill: i: %u, v>>1: %u, v&1: %u, w>>1: %u, w&1: %u\n", i, v>>1, v&1, w>>1, w&1); } } @@ -3731,7 +3892,7 @@ uint32_t is_circle, uint64_t* rLen) if(k == edge->a.n) { - fprintf(stderr, "####ERROR2: i: %u, v>>1: %u, v&1: %u, w>>1: %u, w&1: %u\n", + fprintf(stderr, "####ERROR2-fill: i: %u, v>>1: %u, v&1: %u, w>>1: %u, w&1: %u\n", i, v>>1, v&1, w>>1, w&1); } } @@ -3754,10 +3915,101 @@ uint32_t is_circle, uint64_t* rLen) } +void collect_trans_purge_cov(hap_cov_t *cov, ma_ug_t *ug, hap_overlaps* x, uint32_t is_keep_X) +{ + if(ug->u.a[x->xUid].n == 0 || ug->u.a[x->yUid].n == 0) return; + uint64_t *pri = NULL, pri_n, *aux = NULL, aux_n, i, rId, uCov = 0, uLen = 0; + + if(is_keep_X) + { + pri = ug->u.a[x->xUid].a + x->x_beg_id; + pri_n = x->x_end_id - x->x_beg_id; + + aux = ug->u.a[x->yUid].a + x->y_beg_id; + aux_n = x->y_end_id - x->y_beg_id; + } + else + { + pri = ug->u.a[x->yUid].a + x->y_beg_id; + pri_n = x->y_end_id - x->y_beg_id; + + aux = ug->u.a[x->xUid].a + x->x_beg_id; + aux_n = x->x_end_id - x->x_beg_id; + } + + + uCov = uLen = 0; + for (i = 0; i < aux_n; i++) + { + rId = aux[i]>>33; + uCov += cov->cov[rId]; + } + + for (i = 0; i < pri_n; i++) + { + rId = pri[i]>>33; + uLen += cov->read_g->seq[rId].len; + } + + uCov = (uLen == 0? 0 : uCov / uLen); + + + for (i = 0; i < pri_n; i++) + { + rId = pri[i]>>33; + cov->cov[rId] += (uCov * cov->read_g->seq[rId].len); + } +} + + +void collect_trans_purge_joint_cov(hap_cov_t *cov, ma_ug_t *ug, hap_overlaps* x) +{ + if(ug->u.a[x->xUid].n == 0 || ug->u.a[x->yUid].n == 0) return; + uint64_t *a[2], a_n[2], uCov[2], uLen[2], uDepth[2], i, rId; + + a[0] = ug->u.a[x->xUid].a + x->x_beg_id; + a_n[0] = x->x_end_id - x->x_beg_id; + uCov[0] = uLen[0] = 0; + for (i = 0; i < a_n[0]; i++) + { + rId = a[0][i]>>33; + uCov[0] += cov->cov[rId]; + uLen[0] += cov->read_g->seq[rId].len; + } + + a[1] = ug->u.a[x->yUid].a + x->y_beg_id; + a_n[1] = x->y_end_id - x->y_beg_id; + uCov[1] = uLen[1] = 0; + for (i = 0; i < a_n[1]; i++) + { + rId = a[1][i]>>33; + uCov[1] += cov->cov[rId]; + uLen[1] += cov->read_g->seq[rId].len; + } + + uDepth[0] = (uLen[0] == 0? 0 : uCov[1] / uLen[0]); + uDepth[1] = (uLen[1] == 0? 0 : uCov[0] / uLen[1]); + + for (i = 0; i < a_n[0]; i++) + { + rId = a[0][i]>>33; + cov->cov[rId] += (uDepth[0] * cov->read_g->seq[rId].len); + } + + for (i = 0; i < a_n[1]; i++) + { + rId = a[1][i]>>33; + cov->cov[rId] += (uDepth[1] * cov->read_g->seq[rId].len); + } +} + + + void purge_merge(asg_t *purge_g, ma_ug_t *ug, hap_overlaps_list* all_ovlp, buf_t* b_0, R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, uint64_t* position_index, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, -kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit) +kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit, +hap_cov_t *cov) { uint32_t i, nv, k, v, w, x_beg_index, x_end_index, y_beg_index, y_end_index, cut_beg, cut_end, begIndex, endIndex, keepUid; hap_overlaps *x = NULL/**, *y = NULL**/; @@ -3780,7 +4032,6 @@ kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edg { for (k = 0; k < xReads->n; k++) { - ///aim[query->n - j - 1] = (query->a[j])^(uint64_t)(0x100000000); kv_push(uint64_t, buffer, (xReads->a[xReads->n - k - 1])^(uint64_t)(0x100000000)); } } @@ -3819,9 +4070,6 @@ kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edg endIndex = x->x_end_id-1; if(cut_end < endIndex) endIndex = cut_end; - // get_node_boundary(ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, - // min_ovlp, xReads, yReads, v>>1, w>>1, begIndex, endIndex, x->y_beg_id, x->y_end_id-1, v&1, - // x->rev, &t_forward, &t_backward); get_node_boundary_advance(ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, min_ovlp, xReads, yReads, v>>1, w>>1, begIndex, endIndex, x->y_beg_id, x->y_end_id-1, v&1, x->rev, u_buffer, tailIndex, prevIndex, &t_forward, &t_backward); @@ -3874,10 +4122,33 @@ kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edg } purge_g->seq[w>>1].c = ALTER_LABLE; + collect_trans_purge_joint_cov(cov, ug, x); + + // if(buffer.n > 1) + // { + // for (k = 0; k < buffer.n - 1; k++) + // { + // if((buffer.a[k]>>32) == 854769 && (buffer.a[k+1]>>32) == 64486) + // { + // fprintf(stderr, "+++++++v: %u, w: %u, xReads->n: %u, yReads->n: %u\n", + // v, w, xReads->n, yReads->n); + // fprintf(stderr, "x->rev: %u, x->x_beg_id: %u, x->x_end_id: %u, x->y_beg_id: %u, x->y_end_id: %u\n", + // x->rev, x->x_beg_id, x->x_end_id, x->y_beg_id, x->y_end_id); + // fprintf(stderr, "t_forward.ul>>32: %u, t_forward.v: %u, y_beg_index: %u, y_end_index: %u\n", + // t_forward.ul>>32, t_forward.v, y_beg_index, y_end_index); + // fprintf(stderr, "type: %u, x->x_beg_pos: %u, x->x_end_pos: %u, xReads->len: %u\n", + // x->type, x->x_beg_pos, x->x_end_pos, xReads->len); + // fprintf(stderr, "x->y_beg_pos: %u, x->y_end_pos: %u, yReads->len: %u\n", + // x->y_beg_pos, x->y_end_pos, yReads->len); + // } + // } + // } } - + // fprintf(stderr, "+keepUid: %u, i: %u, b_0->b.n: %u, buffer.n: %u\n", + // keepUid, i, (uint32_t)b_0->b.n, (uint32_t)buffer.n); fill_unitig(buffer.a, buffer.n, read_g, edge, 0, &totalLen); + ///fprintf(stderr, "-keepUid: %u\n", keepUid); xReads = &(ug->u.a[keepUid]); free(xReads->a); @@ -3922,11 +4193,8 @@ kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edg if(purge_g->seq[v>1].c != ALTER_LABLE) continue; asg_seq_drop(purge_g, v>1); } - } - - 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; @@ -3953,7 +4221,6 @@ void collect_reverse_unitig_pair(hc_links* link, ma_ug_t *ug, hap_overlaps* t) push_hc_edge(&(link->a.a[pre_1]), pre_0, 1, 1, &d); } } - } @@ -3975,7 +4242,7 @@ void link_unitigs(asg_t *purge_g, ma_ug_t *ug, hap_overlaps_list* all_ovlp, R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, uint64_t* position_index, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit, -hc_links* link) +hc_links* link, hap_cov_t *cov) { uint32_t v, n_vtx = purge_g->n_seq * 2, beg, end; long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen; @@ -3998,7 +4265,7 @@ hc_links* link) if(link) collect_reverse_unitigs_purge(&b_0, link, ug, all_ovlp); purge_merge(purge_g, ug, all_ovlp, &b_0, ruIndex, reverse_sources, coverage_cut, - read_g, position_index, u_buffer, tailIndex, prevIndex,max_hang, min_ovlp, edge, visit); + read_g, position_index, u_buffer, tailIndex, prevIndex,max_hang, min_ovlp, edge, visit, cov); } free(b_0.b.a); } @@ -4219,29 +4486,143 @@ uint32_t minLen, double purge_threshold) return 0; } +int cmp_chain_score(const void * a, const void * b) +{ + if((*(hap_overlaps*)a).score < (*(hap_overlaps*)b).score) return 1; + if((*(hap_overlaps*)a).score > (*(hap_overlaps*)b).score) return -1; + + return 0; +} +long long get_ovlp_len(long long a_beg, long long a_end, long long b_beg, long long b_end) +{ + long long ovlp = (long long)(MIN(a_end, b_end)) - (long long)(MAX(a_beg, b_beg)) + 1; + return ovlp <= 0? 0 : ovlp; +} +void sort_hap_chain(hap_overlaps_list* all_ovlp) +{ + hap_overlaps *x = NULL, *p = NULL; + uint32_t v, i, k, uId; + long long ovlp, xLen, pLen; + kvec_t(hap_overlaps) pri; kv_init(pri); + kvec_t(hap_overlaps) alt; kv_init(alt); + + for (v = 0; v < all_ovlp->num; v++) + { + uId = v; + qsort(all_ovlp->x[uId].a.a, all_ovlp->x[uId].a.n, sizeof(hap_overlaps), cmp_chain_score); + pri.n = alt.n = 0; + for (i = 0; i < all_ovlp->x[uId].a.n; i++) + { + x = &(all_ovlp->x[uId].a.a[i]); + xLen = x->x_end_pos - x->x_beg_pos; + for (k = 0; k < pri.n; k++) + { + p = &(pri.a[k]); + pLen = p->x_end_pos - p->x_beg_pos; + ovlp = get_ovlp_len(x->x_beg_pos, x->x_end_pos-1, p->x_beg_pos, p->x_end_pos-1); + if(ovlp == 0) continue; + if(ovlp >= (MIN(xLen, pLen))*0.5) break; + } + + if(k < pri.n) + { + x->xUid = k; + kv_push(hap_overlaps, alt, *x); + } + else + { + kv_push(hap_overlaps, pri, *x); + } + } + + } + + kv_destroy(pri); kv_destroy(alt); +} + +void remove_contained_haplotig(hap_overlaps_list* all_ovlp, ma_ug_t *ug, asg_t* nsg, asg_t *purge_g, hc_links* link, hap_cov_t *cov) +{ + uint32_t v, i, uId, xUid; + hap_overlaps *p = NULL; + for (v = 0; v < all_ovlp->num; v++) + { + uId = v; p = NULL; + for (i = 0; i < all_ovlp->x[uId].a.n; i++) + { + if(p == NULL || p->score < all_ovlp->x[uId].a.a[i].score) + { + p = &(all_ovlp->x[uId].a.a[i]); + } + } + + for (i = 0; i < all_ovlp->x[uId].a.n; i++) + { + if(all_ovlp->x[uId].a.a[i].type == YCX) + { + if(!filter_secondary_chain(p->score, all_ovlp->x[uId].a.a[i].score, 0.95)) + { + continue; + } + + xUid = all_ovlp->x[uId].a.a[i].xUid; + + nsg->seq[xUid].c = ALTER_LABLE; + purge_g->seq[xUid].c = ALTER_LABLE; + purge_g->seq[xUid].del = 1; + + all_ovlp->x[uId].a.a[i].status = DELETE; + if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp->x[uId].a.a[i])); + collect_trans_purge_cov(cov, ug, &(all_ovlp->x[uId].a.a[i]), 0); + } + + ///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i])); + } + } + + // for (v = 0; v < all_ovlp.num; v++) + // { + // uId = v; + // for (i = 0; i < all_ovlp.x[uId].a.n; i++) + // { + // if(all_ovlp.x[uId].a.a[i].type == YCX) + // { + // nsg->seq[all_ovlp.x[uId].a.a[i].xUid].c = ALTER_LABLE; + // purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].c = ALTER_LABLE; + // purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].del = 1; + // all_ovlp.x[uId].a.a[i].status = DELETE; + // if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i])); + // collect_trans_purge_cov(cov, ug, &(all_ovlp.x[uId].a.a[i]), 0); + // } + + // if(all_ovlp.x[uId].a.a[i].type == XCY) + // { + // nsg->seq[all_ovlp.x[uId].a.a[i].yUid].c = ALTER_LABLE; + // purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].c = ALTER_LABLE; + // purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].del = 1; + // all_ovlp.x[uId].a.a[i].status = DELETE; + // if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i])); + // collect_trans_purge_cov(cov, ug, &(all_ovlp.x[uId].a.a[i]), 1); + // } + // ///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i])); + // } + // } +} + + + void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio, -uint32_t just_contain, uint32_t just_coverage, hc_links* link) +uint32_t just_contain, uint32_t just_coverage, hc_links* link, hap_cov_t *cov) { asg_t *purge_g = NULL; purge_g = asg_init(); asg_t* nsg = ug->g; uint32_t v, rId, uId, i, offset; ma_utg_t* reads = NULL; - - // kvec_t_u64_warp u_vecs; - // kv_init(u_vecs.a); - // uint8_t* visit = NULL; - // visit = (uint8_t*)malloc(sizeof(uint8_t) * nsg->n_seq); - // memset(visit, 0, nsg->n_seq); - // uint64_t* vote_counting = (uint64_t*)malloc(sizeof(uint64_t)*nsg->n_seq); - // memset(vote_counting, 0, sizeof(uint64_t)*nsg->n_seq); - // kvec_asg_arc_t_offset u_buffer; - // kv_init(u_buffer.a); - // kvec_hap_candidates u_can; - // kv_init(u_can.a); - uint64_t* position_index = (uint64_t*)malloc(sizeof(uint64_t)*read_g->n_seq); + uint64_t* position_index = NULL; + if(cov) position_index = cov->pos_idx; + else position_index = (uint64_t*)malloc(sizeof(uint64_t)*read_g->n_seq); memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq); hap_overlaps_list all_ovlp; @@ -4249,8 +4630,7 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link) hap_overlaps_list back_all_ovlp; init_hap_overlaps_list(&back_all_ovlp, nsg->n_seq); ///uint32_t junk_cov, hap_cov, dip_cov, junk_occ, repeat_occ, single_cov; - asg_arc_t t; - asg_arc_t* p = NULL; + asg_arc_t t, *p = NULL; int r; hap_alignment_struct_pip hap_buf; long long k_mer_only, coverage_only; @@ -4295,7 +4675,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); + 0.1, &all_ovlp, cov); if(hap_buf.cov_threshold < 0) { @@ -4315,56 +4695,18 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link) fprintf(stderr, "[M::%s] purge duplication coverage threshold: %lld\n", __func__, hap_buf.cov_threshold); if(just_coverage) goto end_coverage; - ///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); ///if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp); - - - // for (v = 0; v < nsg->n_seq; v++) - // { - // uId = v; - // if(nsg->seq[uId].del || nsg->seq[uId].c == ALTER_LABLE) continue; - - // hap_alignment(ug, read_g, reverse_sources, ruIndex, coverage_cut, position_index, - // vote_counting, visit, &u_vecs, &u_buffer, &u_can, uId, density, max_hang, min_ovlp, - // 0.05, &all_ovlp); - // } - - filter_hap_overlaps_by_length(&all_ovlp, purege_minLen); ///normalize_hap_overlaps(&all_ovlp, &back_all_ovlp); normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex); ///debug_hap_overlaps(&all_ovlp, &back_all_ovlp); + remove_contained_haplotig(&all_ovlp, ug, nsg, purge_g, link, cov); - for (v = 0; v < all_ovlp.num; v++) - { - uId = v; - for (i = 0; i < all_ovlp.x[uId].a.n; i++) - { - if(all_ovlp.x[uId].a.a[i].type == YCX) - { - nsg->seq[all_ovlp.x[uId].a.a[i].xUid].c = ALTER_LABLE; - purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].c = ALTER_LABLE; - purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].del = 1; - all_ovlp.x[uId].a.a[i].status = DELETE; - if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i])); - } - - if(all_ovlp.x[uId].a.a[i].type == XCY) - { - nsg->seq[all_ovlp.x[uId].a.a[i].yUid].c = ALTER_LABLE; - purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].c = ALTER_LABLE; - purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].del = 1; - all_ovlp.x[uId].a.a[i].status = DELETE; - if(link) collect_reverse_unitig_pair(link, ug, &(all_ovlp.x[uId].a.a[i])); - } - ///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i])); - } - } if(just_contain == 0) { @@ -4375,6 +4717,7 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link) for (i = 0; i < all_ovlp.x[uId].a.n; i++) { if(all_ovlp.x[uId].a.a[i].status == DELETE) continue; + ///if(all_ovlp.x[uId].a.a[i].type == ) if(purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].c == ALTER_LABLE|| purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].del|| purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].c == ALTER_LABLE|| @@ -4383,28 +4726,27 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link) continue; } + ///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i])); + r = get_hap_arch(&(all_ovlp.x[uId].a.a[i]), ug->u.a[all_ovlp.x[uId].a.a[i].xUid].len, ug->u.a[all_ovlp.x[uId].a.a[i].yUid].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); - - if (r >= 0) - { - ///push node? - p = asg_arc_pushp(purge_g); - *p = t; - } - else - { - print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i])); - fprintf(stderr, "error: uId: %u, i: %u, xUid: %u, yUid: %u\n", - uId, i, all_ovlp.x[uId].a.a[i].xUid, all_ovlp.x[uId].a.a[i].yUid); - } + + // if(all_ovlp.x[uId].a.a[i].xUid == 118 && all_ovlp.x[uId].a.a[i].yUid == 82) + // { + // fprintf(stderr, "r: %d\n", r); + // print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i])); + // } + + if(r < 0) continue; + p = asg_arc_pushp(purge_g); + *p = t; } } asg_cleanup(purge_g); asg_symm(purge_g); - + ///may need to do transitive reduction clean_purge_graph(purge_g, bubble_dist, drop_ratio); // if(debug_enable) print_purge_gfa(ug, purge_g); @@ -4412,7 +4754,7 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link) link_unitigs(purge_g, ug, &all_ovlp, ruIndex, reverse_sources, coverage_cut, read_g, position_index, &(hap_buf.buf[0].u_buffer), &(hap_buf.buf[0].u_buffer_tailIndex), &(hap_buf.buf[0].u_buffer_prevIndex), - max_hang, min_ovlp, edge, hap_buf.buf[0].visit, link); + max_hang, min_ovlp, edge, hap_buf.buf[0].visit, link, cov); } for (v = 0; v < all_ovlp.num; v++) @@ -4437,12 +4779,130 @@ uint32_t just_contain, uint32_t just_coverage, hc_links* link) destory_hap_overlaps_list(&all_ovlp); destory_hap_overlaps_list(&back_all_ovlp); asg_destroy(purge_g); - free(position_index); - // kv_destroy(u_vecs.a); - // kv_destroy(u_buffer.a); - // kv_destroy(u_can.a); - // free(vote_counting); - // free(visit); + if(cov) memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq); + else free(position_index); + destory_hap_alignment_struct_pip(&hap_buf); } +hap_cov_t* init_hap_cov_t(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* sources, R_to_U* ruIndex, +ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, int max_hang, int min_ovlp) +{ + uint32_t n_ux = ug->g->n_seq, i, k, j, v, rId, tn, is_Unitig, r_i, nv, w, C_bases; + uint8_t *set = NULL; CALLOC(set, read_g->n_seq<<1); + hap_cov_t *x = NULL; CALLOC(x, 1); + x->n = read_g->n_seq; CALLOC(x->cov, x->n); + MALLOC(x->pos_idx, x->n); memset(x->pos_idx, -1, x->n*sizeof(uint64_t)); + x->reverse_sources = reverse_sources; + x->coverage_cut = coverage_cut; + x->ruIndex = ruIndex; + x->max_hang = max_hang; + x->min_ovlp = min_ovlp; + x->read_g = read_g; + kv_init(x->u_buffer.a); + kv_init(x->tailIndex.a); + kv_init(x->prevIndex.a); + ma_utg_t* u = NULL; + asg_arc_t *av = NULL; + ma_hit_t *h = NULL; + + for (i = 0; i < n_ux; i++) + { + if(ug->g->seq[i].del) continue; + u = &(ug->u.a[i]); + for (r_i = 0; r_i < u->n; r_i++) set[u->a[r_i]>>33] = 1; + + v = i<<1; + nv = asg_arc_n(ug->g, v); + av = asg_arc_a(ug->g, v); + for (k = 0; k < nv; k++) + { + w = av[k].v; + if(av[k].del) continue; + if(ug->g->seq[w>>1].del) continue; + u = &(ug->u.a[w>>1]); + for (r_i = 0; r_i < u->n; r_i++) set[u->a[r_i]>>33] = 1; + } + + v = (i<<1)+1; + nv = asg_arc_n(ug->g, v); + av = asg_arc_a(ug->g, v); + for (k = 0; k < nv; k++) + { + w = av[k].v; + if(av[k].del) continue; + if(ug->g->seq[w>>1].del) continue; + u = &(ug->u.a[w>>1]); + for (r_i = 0; r_i < u->n; r_i++) set[u->a[r_i]>>33] = 1; + } + + + + u = &(ug->u.a[i]); + for (k = 0; k < u->n; k++) + { + C_bases = 0; + rId = u->a[k]>>33; + for (j = 0; j < (uint64_t)(sources[rId].length); j++) + { + h = &(sources[rId].buffer[j]); + tn = Get_tn((*h)); + if(read_g->seq[tn].del == 1) + { + ///get the id of read that contains it + get_R_to_U(ruIndex, tn, &tn, &is_Unitig); + if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue; + } + if(read_g->seq[tn].del == 1) continue; + if(!set[tn]) continue; + C_bases += (Get_qe((*h)) - Get_qs((*h))); + } + x->cov[rId] = MAX(C_bases, x->cov[rId]); + } + + + + u = &(ug->u.a[i]); + for (r_i = 0; r_i < u->n; r_i++) set[u->a[r_i]>>33] = 0; + + v = i<<1; + nv = asg_arc_n(ug->g, v); + av = asg_arc_a(ug->g, v); + for (k = 0; k < nv; k++) + { + w = av[k].v; + if(av[k].del) continue; + if(ug->g->seq[w>>1].del) continue; + u = &(ug->u.a[w>>1]); + for (r_i = 0; r_i < u->n; r_i++) set[u->a[r_i]>>33] = 0; + } + + v = (i<<1)+1; + nv = asg_arc_n(ug->g, v); + av = asg_arc_a(ug->g, v); + for (k = 0; k < nv; k++) + { + w = av[k].v; + if(av[k].del) continue; + if(ug->g->seq[w>>1].del) continue; + u = &(ug->u.a[w>>1]); + for (r_i = 0; r_i < u->n; r_i++) set[u->a[r_i]>>33] = 0; + } + } + + free(set); + return x; +} + +void destory_hap_cov_t(hap_cov_t **x) +{ + if(*x) + { + free((*x)->cov); + free((*x)->pos_idx); + kv_destroy((*x)->u_buffer.a); + kv_destroy((*x)->tailIndex.a); + kv_destroy((*x)->prevIndex.a); + free((*x)); + } +} \ No newline at end of file diff --git a/Purge_Dups.h b/Purge_Dups.h index 4635e22..9928624 100644 --- a/Purge_Dups.h +++ b/Purge_Dups.h @@ -11,14 +11,20 @@ #define HET_PEAK_RATE (HOM_PEAK_RATE*2) #define ALTER_COV_THRES 0.9 #define REAL_ALTER_THRES 0.1 +#define CHAIN_FILTER_RATE 0.7 void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t purege_minLen, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio, -uint32_t just_contain, uint32_t just_coverage, hc_links* link); +uint32_t just_contain, uint32_t just_coverage, hc_links* link, hap_cov_t *cov); void fill_unitig(uint64_t* buffer, uint32_t bufferLen, asg_t* read_g, kvec_asg_arc_t_warp* edge, uint32_t is_circle, uint64_t* rLen); void get_contig_length(ma_ug_t *ug, asg_t *g, uint64_t* primaryLen, uint64_t* alterLen); void enable_debug_mode(uint32_t mode); +hap_cov_t* init_hap_cov_t(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* sources, R_to_U* ruIndex, +ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, int max_hang, int min_ovlp); +void destory_hap_cov_t(hap_cov_t **x); +void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd); + #endif \ No newline at end of file diff --git a/hic.cpp b/hic.cpp index e6e0621..f35724a 100644 --- a/hic.cpp +++ b/hic.cpp @@ -38,6 +38,7 @@ typedef struct{ kvec_t(uint64_t) name_Len; kvec_t(char) r; kvec_t(uint64_t) r_Len; + uint64_t idx; } reads_t; typedef struct{ @@ -167,6 +168,18 @@ typedef struct { // global data structure for kt_pipeline() uint64_t n_thread; } pldat_t; +typedef struct { + uint64_t *a, id; + uint16_t occ1, occ2; +} pe_hit_hap; + +typedef struct { + pe_hit_hap* a; + size_t n, m; + uint64_t n_u; +} kvec_pe_hit_hap; + + typedef struct { uint64_t s, e, id; } pe_hit; @@ -190,6 +203,18 @@ KRADIX_SORT_INIT(u32, uint32_t, generic_key, 4) #define g_partition_key(x) (((x)>>1)+((x)<<63)) KRADIX_SORT_INIT(g_partition, uint64_t, g_partition_key, 8) +#define get_pe_s(x) ((x).a[0]) +#define get_pe_e(x) ((x).a[(x).occ1]) +KRADIX_SORT_INIT(pe_an1, pe_hit_hap, get_pe_s, 8) +KRADIX_SORT_INIT(pe_an2, pe_hit_hap, get_pe_e, 8) + +#define pe_occ_key_1(x) ((x).occ1) +KRADIX_SORT_INIT(pe_occ1, pe_hit_hap, pe_occ_key_1, member_size(pe_hit_hap, occ1)) +#define pe_occ_key_2(x) ((x).occ2) +KRADIX_SORT_INIT(pe_occ2, pe_hit_hap, pe_occ_key_2, member_size(pe_hit_hap, occ2)) +#define pe_occ_key_t(x) (((uint64_t)((x).occ1))+((uint64_t)((x).occ2))) +KRADIX_SORT_INIT(pe_occ_t, pe_hit_hap, pe_occ_key_t, 8) + typedef struct { // global data structure for kt_pipeline() const ha_ug_index* idx; @@ -198,7 +223,8 @@ typedef struct { // global data structure for kt_pipeline() uint64_t n_thread; uint64_t total_base; uint64_t total_pair; - kvec_pe_hit hits; + ///kvec_pe_hit hits; + kvec_pe_hit_hap hits; hc_links* link; } sldat_t; @@ -218,7 +244,8 @@ typedef struct { // data structure for each step in kt_pipeline() char **seq; ch_buf_t *buf; kvec_vote* pos_buf; - pe_hit* pos; + ///pe_hit* pos; + pe_hit_hap* pos; hc_links* link; } stepdat_t; @@ -232,6 +259,8 @@ KRADIX_SORT_INIT(hc_pos, uint64_t, hc_pos_key, 8) KRADIX_SORT_INIT(hc_s_hit_an1, s_hit, hc_s_hit_an1_key, 8) #define hc_s_hit_an2_key(a) ((uint32_t)(a).off_cnt) KRADIX_SORT_INIT(hc_s_hit_an2, s_hit, hc_s_hit_an2_key, 8) +#define hc_s_hit_off_cnt_key(a) ((a).off_cnt) +KRADIX_SORT_INIT(hc_s_hit_off_cnt, s_hit, hc_s_hit_off_cnt_key, 8) #define hc_edge_key_u(a) ((a).uID) KRADIX_SORT_INIT(hc_edge_u, hc_edge, hc_edge_key_u, 4) #define hc_edge_key_d(a) ((a).dis) @@ -870,11 +899,12 @@ uint64_t debug_hash_value(char *r, uint64_t end, uint64_t k_mer) inline uint64_t collect_votes(s_hit* a, uint64_t n) { if(n == 0) return 0; - if(n == 1) return (a[0].off_cnt>>32); + if(n == 1) return (a[0].off_cnt>>32); //seed length long long i = 0; - uint64_t cur_beg, cur_end, beg, end, ovlp = 0, tLen = 0; + uint64_t cur_beg, cur_end, beg, end, ovlp = 0, tLen = 0; cur_end = (uint32_t)a[n-1].off_cnt; cur_beg = cur_end + 1 - (a[n-1].off_cnt>>32); + for (i = n - 2; i >= 0; i--) { @@ -1164,13 +1194,198 @@ const ha_ug_index* idx, uint64_t buf_iter, uint64_t rid) /*******************************for debug************************************/ } +uint64_t get_longest_hit(char *r, uint64_t len, uint64_t k_mer, uint64_t self_p, uint64_t self_rev, kvec_vote* buf, const ha_ug_index* idx, +uint64_t *pos_list, uint64_t cnt, uint64_t* c_sfx) +{ + uint64_t max_p, map_p_occ, i, j, m, rev, ref_p, u_len, uID, k_len; + s_hit *p = NULL; + ///each k-mer at different unitigs + ///rev:uID:pos + if(c_sfx) (*c_sfx) = (uint64_t)-1; + for (j = 0; j < cnt; j++) + { + ///get + kv_pushp(s_hit, buf->a, &p); + rev = (pos_list[j]>>63) != self_rev; + ref_p = pos_list[j] & idx->pos_mode; + uID = (pos_list[j] << 1) >> (64 - idx->uID_bits); + u_len = idx->ug->u.a[uID].len; + if(rev) ref_p = u_len - 1 - (ref_p + 1 - k_mer); + p->off_cnt = self_p | ((uint64_t)k_mer << 32); ///high bits should be the legnth + + p->ref = ref_p >= self_p? (ref_p-self_p) + : (self_p-ref_p) + ((uint64_t)1 << (idx->pos_bits - 1)); + p->ref = (rev << 63)|(pos_list[j] & idx->uID_mode)|(p->ref&idx->pos_mode); + + + ///extend + k_len = check_exact_match(r, self_p + 1, len, idx->ug->u.a[uID].s, ref_p + 1, u_len, len, rev, 0); + if(c_sfx && cnt == idx->hap_cnt && k_len < (*c_sfx)) (*c_sfx) = k_len; + + p->off_cnt += ((uint64_t)k_len << 32) + k_len; + + if(self_p >= k_mer && ref_p >= k_mer) + { + k_len = check_exact_match(r, self_p - k_mer, len, idx->ug->u.a[uID].s, + ref_p - k_mer, u_len, len, rev, 1); + p->off_cnt += ((uint64_t)k_len << 32); + } + // if(cnt > 0) fprintf(stderr, "inner j: %lu, rev: %lu, uID: %lu, ref_p: %lu, self_p: %u, len: %lu\n", j, rev, uID, ref_p, (uint32_t)p->off_cnt, p->off_cnt>>32); + } + + p = buf->a.a + buf->a.n - cnt; + if(cnt > 1) radix_sort_hc_s_hit_off_cnt(p, p + cnt); + max_p = map_p_occ = 0; + for (j = 1, i = 0; j <= cnt; ++j) + { + if(j == cnt || p[j].off_cnt != p[i].off_cnt) + { + if((max_p>>32) < (p[i].off_cnt>>32)) + { + max_p = p[i].off_cnt; + map_p_occ = j - i; + } + else if(((max_p>>32) == (p[i].off_cnt>>32)) && ((j - i) > map_p_occ)) + { + max_p = p[i].off_cnt; + map_p_occ = j - i; + } + i = j;///must + } + } + + buf->a.n -= cnt; + for (j = m = 0; j < cnt; j++) + { + if(p[j].off_cnt == max_p) + { + p[m] = p[j]; + m++; + } + } + cnt = m; + buf->a.n += cnt; + + // if(cnt > 0) fprintf(stderr, "max_p_offset: %u, max_p_len: %lu, map_p_occ: %lu\n", (uint32_t)max_p, max_p>>32, map_p_occ); + return max_p; +} + +#define is_update_hit(mL, mR, cL, cR) (((mL)<(cL))||((mL)==(cL)&&(mR)<(cR))) +inline void compress_mapped_pos_advance(const ha_ug_index* idx, kvec_vote* buf, uint64_t buf_iter, uint64_t ovlp_thre) +{ + if(buf_iter >= buf->a.n) + { + buf->a.n = buf_iter; + return; + } + s_hit *p = NULL; + uint64_t rev, uID, ref_p, self_p, eLen, tLen, i, j, cnt; + uint64_t max_beg = 0, max_end = 0, max_i, max_occ, cur_beg, cur_end, ovlp; + uint64_t second_i = (uint64_t)-1, second_occ; + uint64_t max_eLen, sec_eLen; + double max_eRate, sec_eRate, eRate; + p = buf->a.a + buf_iter; + cnt = buf->a.n - buf_iter; + radix_sort_hc_s_hit_off_cnt(p, p + cnt); ///buf save all hits, here sort by offset in reads + max_eLen = 0; max_i = (uint64_t)-1; max_occ = 0; max_eRate = 0; + for (j = 1, i = 0; j <= cnt; ++j) + { + if(j == cnt || p[j].off_cnt != p[i].off_cnt) + { + ///occ = j - i; + interpret_pos((ha_ug_index*)idx, &p[i], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + eRate = (double)(eLen)/(double)(tLen); + if(is_update_hit(max_eLen, max_eRate, eLen, eRate)) + { + max_eLen = eLen; max_eRate = eRate; + max_end = self_p; max_beg = self_p + 1 - tLen; + max_i = i; max_occ = j - i; + } + // fprintf(stderr, "\n++++++[%lu, %lu] uID: %lu, ref_p: %lu, self_p: %lu, eLen: %lu, tLen: %lu, max_i: %lu\n", + // i, j, uID, ref_p, self_p, eLen, tLen, max_i); + i = j;///must + } + } + + + sec_eLen = 0; second_i = (uint64_t)-1; second_occ = 0; sec_eRate = 0; + for (j = 1, i = 0; j <= cnt; ++j) + { + if(j == cnt || p[j].off_cnt != p[i].off_cnt) + { + if(i != max_i) + { + interpret_pos((ha_ug_index*)idx, &p[i], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + eRate = (double)(eLen)/(double)(tLen); + cur_end = self_p; + cur_beg = self_p + 1 - tLen; + // fprintf(stderr, "\n----[%lu, %lu] uID: %lu, ref_p: %lu, self_p: %lu, eLen: %lu, tLen: %lu, max_i: %lu\n", + // i, j, uID, ref_p, self_p, eLen, tLen, max_i); + // fprintf(stderr, "max_beg: %lu, max_end: %lu, cur_beg: %lu, cur_end: %lu\n", + // max_beg, max_end, cur_beg, cur_end); + ///overlap with max interval + if(MAX(cur_beg, max_beg) <= MIN(cur_end, max_end)) + { + ovlp = MIN(cur_end, max_end) - MAX(cur_beg, max_beg) + 1; + if(ovlp == MIN(max_end+1-max_end, tLen)) + { + i = j;///must + continue;///fully contain + } + + if(ovlp > ((max_end+1-max_end)*0.8) && eLen > (max_eLen*0.8))///best is not unique + { + buf->a.n = buf_iter; + return; + } + if(ovlp > ((max_end+1-max_end)*0.15)) + { + i = j;///must + continue;///fully contain + } + } + + if(is_update_hit(sec_eLen, sec_eRate, eLen, eRate)) + { + sec_eLen = eLen; sec_eRate = eRate; + second_i = i; second_occ = j - i; + } + } + i = j;///must + } + } + + // fprintf(stderr, "max_i: %lu, max_occ: %lu, second_i: %lu, second_occ: %lu\n", + // max_i, max_occ, second_i, second_occ); + + if(second_i == (uint64_t)-1) + { + i = 0; + for (j = max_i; j < max_i + max_occ; j++, i++) p[i] = p[j]; + } + else ///be carful about overwritten + { + i = 0; + if(max_i <= second_i) + { + for (j = max_i; j < max_i + max_occ; j++, i++) p[i] = p[j]; + for (j = second_i; j < second_i + second_occ; j++, i++) p[i] = p[j]; + } + else + { + for (j = second_i; j < second_i + second_occ; j++, i++) p[i] = p[j]; + for (j = max_i; j < max_i + max_occ; j++, i++) p[i] = p[j]; + } + } + + buf->a.n = buf_iter + max_occ + second_occ; +} void get_alignment(char *r, uint64_t len, uint64_t k_mer, kvec_vote* buf, const ha_ug_index* idx, uint64_t buf_iter, uint64_t rid) { - uint64_t i, j, l = 0, skip, *pos_list = NULL, cnt, rev, self_p, ref_p, u_len, uID; + uint64_t i, j, k, l = 0, k_len, c_sfx, m, skip, *pos_list = NULL, cnt, rev, self_p, ref_p, uID; uint64_t x[4], mask = (1ULL<a.n = 0; for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) { int c = seq_nt4_table[(uint8_t)r[i]]; @@ -1186,144 +1401,75 @@ const ha_ug_index* idx, uint64_t buf_iter, uint64_t rid) { hash = hc_hash_long(x, &skip, k_mer); if(skip == (uint64_t)-1) continue; - /*******************************for debug************************************/ - // if(debug_hash_value(r, i, k_mer) != hash) - // { - // fprintf(stderr, "ERROR\n"); - // } - /*******************************for debug************************************/ cnt = get_hc_pt1_count((ha_ug_index*)idx, hash, &pos_list); - if(cnt > idx->hap_cnt) continue; - if(cnt != 1) continue; ///might be able to be disabled in future + if(cnt > idx->hap_cnt || cnt < 0) continue; - ///rev:uID:pos - for (j = 0; j < cnt; j++) + if(cnt > 1) { - kv_pushp(s_hit, buf->a, &p); - rev = (pos_list[j]>>63) != skip; - self_p = i; - ref_p = pos_list[j] & idx->pos_mode; - uID = (pos_list[j] << 1) >> (64 - idx->uID_bits); - u_len = idx->ug->u.a[uID].len; - if(rev) ref_p = u_len - 1 - (ref_p + 1 - k_mer); - p->off_cnt = self_p | ((uint64_t)k_mer << 32); ///high bits should be the legnth + for (j = 0; j < cnt; j++) + { + uID = (pos_list[j] << 1) >> (64 - idx->uID_bits); + for (k = j + 1; k < cnt; k++) + { + if(uID == ((pos_list[k] << 1) >> (64 - idx->uID_bits))) break; + } + if(k < cnt) break; + } - p->ref = ref_p >= self_p? (ref_p-self_p) - : (self_p-ref_p) + ((uint64_t)1 << (idx->pos_bits - 1)); - p->ref = (rev << 63)|(pos_list[j] & idx->uID_mode)|(p->ref&idx->pos_mode); - - - /*******************************for debug************************************/ - // if(check_exact_match(r, i + 1 - k_mer, len, - // idx->ug->u.a[uID].s, ref_p + 1 - k_mer, u_len, k_mer, rev, 0) != k_mer - // || - // check_exact_match(r, i, len, - // idx->ug->u.a[uID].s, ref_p, u_len, k_mer, rev, 1) != k_mer) - // { - // fprintf(stderr, "ERROR\n"); - // } - /*******************************for debug************************************/ + if(j < cnt) continue; } - - if(cnt == 1) + // if(cnt > 0) fprintf(stderr, "+i: %lu, l: %lu, cnt: %lu\n", i, l, cnt); + get_longest_hit(r, len, k_mer, i, skip, buf, idx, pos_list, cnt, &c_sfx); + // if(cnt > 0) fprintf(stderr, "c_sfx: %lu\n", c_sfx); + if(c_sfx != (uint64_t)-1) { - ///uint64_t debug_right = 0, debug_left = 0, debug_len; - - j = check_exact_match(r, self_p + 1, len, idx->ug->u.a[uID].s, ref_p + 1, u_len, len, rev, 0); - - ///debug_right = j; - ///if(j == 0) continue; - if((j + 1) >= k_mer) + k_len = c_sfx; + if((k_len + 1) >= k_mer) { l = 0, x[0] = x[1] = x[2] = x[3] = 0; - i = i + j - (k_mer - 1); + i = i + k_len - (k_mer - 1); } else { - ///l = i - (i + j - (k_mer - 1)); - l = k_mer - j -1; - } - buf->a.a[buf->a.n-1].off_cnt += ((uint64_t)j << 32) + j; - - if(self_p >= k_mer && ref_p >= k_mer) - { - j = check_exact_match(r, self_p - k_mer, len, idx->ug->u.a[uID].s, - ref_p - k_mer, u_len, len, rev, 1); - buf->a.a[buf->a.n-1].off_cnt += ((uint64_t)j << 32); - ///debug_left = j; - } - - - // debug_len = check_exact_match(r, self_p + debug_right, len, idx->ug->u.a[uID].s, - // ref_p + debug_right, u_len, len, rev, 1); - // if(debug_len!= (debug_left + debug_right + k_mer)) - // { - // fprintf(stderr, "debug_len: %lu, debug_left: %lu, debug_right: %lu\n", - // debug_len, debug_left, debug_right); - // } + ///l = i - (i + k_len - (k_mer - 1)); + l = k_mer - k_len - 1; + } } - + // if(cnt > 0) fprintf(stderr, "-i: %lu, l: %lu\n", i, l); } } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart } - ///if(buf->a.n - buf_iter <= 1) return; if(buf->a.n - buf_iter == 0) return; if(buf->a.n - buf_iter > 1) radix_sort_hc_s_hit_an1(buf->a.a + buf_iter, buf->a.a + buf->a.n); - - /*******************************for debug************************************/ - // print_pos_list(idx, buf->a.a+buf_iter, buf->a.n - buf_iter, rid, (buf_iter != 0)); - // fprintf(stderr, "len0:%lu\n", buf->a.n - buf_iter); - // for (i = buf_iter; i < buf->a.n; i++) - // { - // interpret_pos(idx, &buf->a.a[i], &rev, &uID, &ref_p, &self_p, &cnt, NULL); - // fprintf(stderr, "(%lu) rev: %lu, uID: %lu, ref_p: %lu, self_p: %lu, len: %lu\n", - // i, rev, uID, ref_p, self_p, cnt); - // } - /*******************************for debug************************************/ - - - - - - uint64_t cur_ref_p, thres = (len * HIC_R_E_RATE) + 1, m, index_beg, ovlp, maxLen = 0, max_i = (uint64_t)-1; + uint64_t cur_ref_p, thres = (len * HIC_R_E_RATE) + 1, index_beg, ovlp; i = m = buf_iter; while (i < buf->a.n) { interpret_pos(idx, &buf->a.a[i], &rev, &uID, &ref_p, &self_p, &cnt, NULL); - /*******************************for debug************************************/ - // if(check_exact_match(r, self_p, len, idx->ug->u.a[uID].s, - // ref_p, idx->ug->u.a[uID].len, cnt, rev, 1) != cnt) - // { - // fprintf(stderr, "ERROR\n"); - // } - /*******************************for debug************************************/ - // if(self_p > ref_p) - // { - // i++; - // continue; ///fix this in future - // } + ///fprintf(stderr, "after-i: %lu, uID: %lu, ref_p: %lu, self_p: %lu\n", i, uID, ref_p, self_p); cur_ref_p = buf->a.a[i].ref; index_beg = i; - while ((i < buf->a.n) && ((buf->a.a[i].ref>>idx->pos_bits) == (cur_ref_p>>idx->pos_bits)) && + ///ref>>(idx->pos_bits-1) = (rev:1):(uID:uID-bits):(ref_pos>=self_pos:1) + while ((i < buf->a.n) && + ((buf->a.a[i].ref>>(idx->pos_bits-1)) == (cur_ref_p>>(idx->pos_bits-1))) && (buf->a.a[i].ref - cur_ref_p <= thres)) { i++; } if(i - index_beg > 1) { - radix_sort_hc_s_hit_an2(buf->a.a + index_beg, buf->a.a + i); + radix_sort_hc_s_hit_an2(buf->a.a + index_beg, buf->a.a + i);//sort by self_p } ovlp = collect_votes(buf->a.a + index_beg, i - index_beg); + ///fprintf(stderr, "i-1: %lu, self_p: %u\n", i-1, (uint32_t)buf->a.a[i - 1].off_cnt); buf->a.a[m] = buf->a.a[i - 1]; buf->a.a[m].off_cnt = (buf->a.a[m].off_cnt << 32)>>32; buf->a.a[m].off_cnt += ((uint64_t)ovlp<<32); - - if(maxLen < (ovlp&((uint64_t)65535))) maxLen = (ovlp&((uint64_t)65535)), max_i = m; - + ///fprintf(stderr, "m: %lu, self_p: %u\n", m, (uint32_t)buf->a.a[m].off_cnt); m++; } buf->a.n = m; @@ -1349,9 +1495,7 @@ const ha_ug_index* idx, uint64_t buf_iter, uint64_t rid) // i, rev, uID, ref_p, self_p, eLen, tLen); // } /*******************************for debug************************************/ - - compress_mapped_pos(idx, buf, buf_iter, max_i, thres); - + compress_mapped_pos_advance(idx, buf, buf_iter, (k_mer * 0.1) > 0? (k_mer * 0.1) : 1); /*******************************for debug************************************/ // fprintf(stderr, "len2:%lu, max_i: %lu\n", buf->a.n - buf_iter, max_i); // for (i = buf_iter; i < buf->a.n; i++) @@ -1394,6 +1538,7 @@ inline int is_unreliable_hits(long long rev, long long ref_p, long long tLen, ui return 0; } + inline void set_pe_pos(ha_ug_index* idx, s_hit *l1, uint64_t occ1, s_hit *l2, uint64_t occ2, pe_hit* x, uint64_t rid, hc_links* link) { @@ -1442,31 +1587,162 @@ pe_hit* x, uint64_t rid, hc_links* link) } } +void get_5_3_list(ha_ug_index* idx, s_hit* p, uint64_t cnt, s_hit** l5, uint64_t* l5_occ, +s_hit** l3, uint64_t* l3_occ) +{ + (*l5) = (*l3) = NULL; + (*l5_occ) = (*l3_occ) = 0; + uint64_t i, j, rev, uID, ref_p, self_p, eLen, tLen, cur_beg, num; + uint64_t beg_5 = (uint64_t)-1; + for (j = 1, i = 0, num = 0; j <= cnt; ++j) + { + if(j == cnt || p[j].off_cnt != p[i].off_cnt) + { + interpret_pos((ha_ug_index*)idx, &p[i], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + ///cur_end = self_p; + cur_beg = self_p + 1 - tLen; + num++; + if(cur_beg <= beg_5) + { + (*l3_occ) = (*l5_occ); (*l3) = (*l5); + beg_5 = cur_beg; (*l5_occ) = j - i; (*l5) = p + i; + } + else + { + (*l3_occ) = j - i; (*l3) = p + i; + } + i = j;///must + } + } + + ///if(num > 2) fprintf(stderr, "ERROR: get_5_3_list\n"); +} +inline void set_pe_pos_hap(ha_ug_index* idx, s_hit *l1, uint64_t occ1, s_hit *l2, uint64_t occ2, +pe_hit_hap* x, uint64_t rid, hc_links* link) +{ + if(occ1 == 0 || occ2 == 0) return; + uint64_t rev, uID, ref_p, self_p, eLen, tLen, i, is_unreliable = 0; + s_hit *l1_5 = NULL, *l1_3 = NULL, *l2_5 = NULL, *l2_3 = NULL; + uint64_t l1_5_occ = 0, l1_3_occ = 0, l2_5_occ = 0, l2_3_occ = 0; + + /***************************for debug******************************/ + // fprintf(stderr, "\nrid: %lu, occ1: %lu, occ2: %lu\n", rid, occ1, occ2); + // for (i = 0; i < occ1; i++) + // { + // interpret_pos(idx, &l1[i], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + // fprintf(stderr, "***-1-rev: %lu, uID: %lu, ref_p: %lu, self_p: %lu\n", + // rev, uID, ref_p, self_p); + // } + + // for (i = 0; i < occ2; i++) + // { + // interpret_pos(idx, &l2[i], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + // fprintf(stderr, "***-2-rev: %lu, uID: %lu, ref_p: %lu, self_p: %lu\n", + // rev, uID, ref_p, self_p); + // } + /***************************for debug******************************/ + + + get_5_3_list(idx, l1, occ1, &l1_5, &l1_5_occ, &l1_3, &l1_3_occ); + get_5_3_list(idx, l2, occ2, &l2_5, &l2_5_occ, &l2_3, &l2_3_occ); + if(l1_5_occ == 0 || l2_5_occ == 0) return; + x->id = rid; + MALLOC(x->a, l1_5_occ + l2_5_occ); + + x->occ1 = 0; + for (i = 0; i < l1_5_occ; i++) + { + interpret_pos(idx, &l1_5[i], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + if(ref_p < self_p) continue; + ref_p -= self_p; + if(rev) ref_p = idx->ug->u.a[uID].len - 1 - ref_p; + if(link && (is_unreliable_hits(rev, ref_p, tLen, uID, link))) + { + is_unreliable = 1; + continue; + } + x->a[x->occ1++] = (rev<<63) | ((uID << (64-idx->uID_bits))>>1) | (ref_p & idx->pos_mode); + } + + x->occ2 = x->occ1; + for (i = 0; i < l2_5_occ; i++) + { + interpret_pos(idx, &l2_5[i], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + if(ref_p < self_p) continue; + ref_p -= self_p; + if(rev) ref_p = idx->ug->u.a[uID].len - 1 - ref_p; + if(link && (is_unreliable_hits(rev, ref_p, tLen, uID, link))) + { + is_unreliable = 1; + continue; + } + x->a[x->occ2++] = (rev<<63) | ((uID << (64-idx->uID_bits))>>1) | (ref_p & idx->pos_mode); + } + x->occ2 -= x->occ1; + + if(x->occ1 == 0 || x->occ2 == 0 || is_unreliable) + { + free(x->a); x->occ1 = x->occ2 = 0; x->a = NULL; x->id = (uint64_t)-1; + return; + } + + + if(x->occ1 > 1) radix_sort_hc64(x->a, x->a + x->occ1); + if(x->occ2 > 1) radix_sort_hc64(x->a + x->occ1, x->a + x->occ1 + x->occ2); + + + /***************************for debug******************************/ + // fprintf(stderr, "-------------saved: x->occ1: %u, x->occ2: %u-------------\n", x->occ1, x->occ2); + // for (i = 0; i < x->occ1; i++) + // { + // fprintf(stderr, "###-1-rev: %lu, uID: %lu, ref_p: %lu\n", + // x->a[i]>>63, (x->a[i]<<1)>>(64-idx->uID_bits), x->a[i] & idx->pos_mode); + // } + + // for (i = 0; i < x->occ2; i++) + // { + // fprintf(stderr, "###-2-rev: %lu, uID: %lu, ref_p: %lu\n", + // x->a[i+x->occ1]>>63, (x->a[i+x->occ1]<<1)>>(64-idx->uID_bits), x->a[i+x->occ1] & idx->pos_mode); + // } + // fprintf(stderr, "-------------get_pe_s-rev: %lu, uID: %lu, ref_p: %lu-------------\n", + // get_pe_s(*x)>>63, (get_pe_s(*x)<<1)>>(64-idx->uID_bits), get_pe_s(*x) & idx->pos_mode); + // fprintf(stderr, "-------------get_pe_e-rev: %lu, uID: %lu, ref_p: %lu-------------\n", + // get_pe_e(*x)>>63, (get_pe_e(*x)<<1)>>(64-idx->uID_bits), get_pe_e(*x) & idx->pos_mode); + /***************************for debug******************************/ +} + +uint64_t if_debug_read(uint64_t rid) +{ + if(rid == 177 || rid == 439 || rid == 97 || rid == 114) + { + return 1; + } + return 0; +} + static void worker_for_alignment(void *data, long i, int tid) // callback for kt_for() { stepdat_t *s = (stepdat_t*)data; - s->pos[i].id = s->pos[i].s = s->pos[i].e = (uint64_t)-1; + s->pos[i].id = (uint64_t)-1; s->pos[i].occ1 = s->pos[i].occ2 = 0; s->pos[i].a = NULL; + + /*******************************for debug************************************/ + // if(!if_debug_read(s->id+i)) return; + // fprintf(stderr, "work-rid: %lu\n", (uint64_t)(s->id+i)); + /*******************************for debug************************************/ + uint64_t len1 = s->len[i]>>32, len2 = (uint32_t)s->len[i], occ1, occ2; char *r1 = s->seq[i], *r2 = s->seq[i] + len1; - /*******************************for debug************************************/ - // if(memcmp(r1, R1.r.a + R1.r_Len.a[s->id+i], len1) != 0) - // { - // fprintf(stderr, "haha1\n"); - // } - // if(memcmp(r2, R2.r.a + R2.r_Len.a[s->id+i], len2) != 0) - // { - // fprintf(stderr, "haha2\n"); - // } - /*******************************for debug************************************/ + // fprintf(stderr, "**********R1**********\n"); s->pos_buf[tid].a.n = 0; get_alignment(r1, len1, s->idx->k, &s->pos_buf[tid], s->idx, 0, s->id+i); occ1 = s->pos_buf[tid].a.n; if(occ1 == 0) return; + // fprintf(stderr, "**********R2**********\n"); get_alignment(r2, len2, s->idx->k, &s->pos_buf[tid], s->idx, occ1, s->id+i); occ2 = s->pos_buf[tid].a.n - occ1; if(occ2 == 0) return; - set_pe_pos((ha_ug_index*)s->idx, s->pos_buf[tid].a.a, occ1, s->pos_buf[tid].a.a + occ1, occ2, &(s->pos[i]), s->id+i, s->link); + set_pe_pos_hap((ha_ug_index*)s->idx, s->pos_buf[tid].a.a, occ1, s->pos_buf[tid].a.a + occ1, occ2, &(s->pos[i]), s->id+i, s->link); /*******************************for debug************************************/ // if(memcmp(r1, R1.r.a + R1.r_Len.a[s->id+i], len1) != 0) @@ -1532,7 +1808,7 @@ static void *worker_pipeline(void *data, int step, void *in) // callback for kt_ else if (step == 1) { // step 2: alignment stepdat_t *s = (stepdat_t*)in; CALLOC(s->pos_buf, p->n_thread); - MALLOC(s->pos, s->n); + CALLOC(s->pos, s->n); int i; kt_for(p->n_thread, worker_for_alignment, s, s->n); for (i = 0; i < s->n; ++i) { @@ -1551,8 +1827,8 @@ static void *worker_pipeline(void *data, int step, void *in) // callback for kt_ stepdat_t *s = (stepdat_t*)in; int i; for (i = 0; i < s->n; ++i) { - if(s->pos[i].s == (uint64_t)-1) continue; - kv_push(pe_hit, p->hits.a, s->pos[i]); + if(s->pos[i].a == NULL) continue; + kv_push(pe_hit_hap, p->hits, s->pos[i]); } free(s->pos); free(s); @@ -1561,38 +1837,53 @@ static void *worker_pipeline(void *data, int step, void *in) // callback for kt_ } -void load_reads(reads_t* x, const char *fn) +int load_reads(reads_t* x, const enzyme *fn1, const enzyme *fn2) { kv_init(x->name); kv_init(x->name_Len); kv_init(x->r); kv_init(x->r_Len); - gzFile fp; - kseq_t *ks; + int ret; uint64_t name_tot, base_total; - + int i; name_tot = base_total = 0; - if ((fp = gzopen(fn, "r")) == 0) return; - ks = kseq_init(fp); - while (((ret = kseq_read(ks)) >= 0)) + + for (i = 0; i < fn1->n && i < fn2->n; i++) { - kv_push(uint64_t, x->name_Len, name_tot); - kv_resize(char, x->name, name_tot + ks->name.l); - memcpy(x->name.a + name_tot, ks->name.s, ks->name.l); - name_tot += ks->name.l; + gzFile fp; + if ((fp = gzopen(fn1->a[i], "r")) == 0) + { + kv_destroy(x->name); + kv_destroy(x->name_Len); + kv_destroy(x->r); + kv_destroy(x->r_Len); + return 0; + } + + kseq_t *ks; + ks = kseq_init(fp); - kv_push(uint64_t, x->r_Len, base_total); - kv_resize(char, x->r, base_total + ks->seq.l); - memcpy(x->r.a + base_total, ks->seq.s, ks->seq.l); - base_total += ks->seq.l; + while (((ret = kseq_read(ks)) >= 0)) + { + kv_push(uint64_t, x->name_Len, name_tot); + kv_resize(char, x->name, name_tot + ks->name.l); + memcpy(x->name.a + name_tot, ks->name.s, ks->name.l); + name_tot += ks->name.l; + + kv_push(uint64_t, x->r_Len, base_total); + kv_resize(char, x->r, base_total + ks->seq.l); + memcpy(x->r.a + base_total, ks->seq.s, ks->seq.l); + base_total += ks->seq.l; + } + + kseq_destroy(ks); + gzclose(fp); } - kv_push(uint64_t, x->name_Len, name_tot); kv_push(uint64_t, x->r_Len, base_total); - - kseq_destroy(ks); - gzclose(fp); + x->idx = 0; + return 1; } @@ -1625,11 +1916,11 @@ void destory_reads(reads_t* x) kv_destroy(x->r_Len); } -void print_hits(ha_ug_index* idx, kvec_pe_hit* hits, const char *fn) +void print_hits(ha_ug_index* idx, kvec_pe_hit* hits, const enzyme *fn1, const enzyme *fn2) { uint64_t k, shif = 64 - idx->uID_bits; reads_t r1; - load_reads(&r1, fn); + load_reads(&r1, fn1, fn2); char dir[2] = {'+', '-'}; for (k = 0; k < hits->a.n; ++k) { @@ -1643,33 +1934,115 @@ void print_hits(ha_ug_index* idx, kvec_pe_hit* hits, const char *fn) destory_reads(&r1); } -void dedup_hits(kvec_pe_hit* hits) +inline void swap_pe_hit_hap(pe_hit_hap* x, pe_hit_hap* y) +{ + pe_hit_hap tmp; + tmp = (*x); (*x) = (*y); (*y) = tmp; +} + +void dedup_hits(kvec_pe_hit_hap* hits, const ha_ug_index* idx) { double index_time = yak_realtime(); - uint64_t k, l, m = 0, cur; - radix_sort_pe_hit_an1(hits->a.a, hits->a.a + hits->a.n); - for (k = 1, l = 0; k <= hits->a.n; ++k) - { - if (k == hits->a.n || hits->a.a[k].s != hits->a.a[l].s) + uint64_t k, l, m = 0, cur = (uint64_t)-1; + radix_sort_pe_an1(hits->a, hits->a + hits->n); + /***************************for debug******************************/ + // for (k = 0; k < hits->n; ++k) + // { + // for (l = k + 1; l < hits->n; l++) + // { + // if(get_pe_s(hits->a[k]) == get_pe_s(hits->a[l]) && + // get_pe_e(hits->a[k]) == get_pe_e(hits->a[l])) + // { + // fprintf(stderr, "DUP: k_id=%lu, l_id=%lu\n", hits->a[k].id, hits->a[l].id); + // } + // } + // } + + /** + fprintf(stderr, "\n\n\n\n\n\n\n\n\n\n*********************dedup_hits*********************\n"); + for (k = 0; k < hits->n; ++k) + { + pe_hit_hap *x = &(hits->a[k]); + fprintf(stderr, "\nsorted-rid: %lu, occ1: %u, occ2: %u\n", x->id, x->occ1, x->occ2); + fprintf(stderr, "---get_pe_s-rev: %lu, uID: %lu, ref_p: %lu---\n", + get_pe_s(*x)>>63, (get_pe_s(*x)<<1)>>(64-idx->uID_bits), get_pe_s(*x) & idx->pos_mode); + fprintf(stderr, "---get_pe_e-rev: %lu, uID: %lu, ref_p: %lu---\n", + get_pe_e(*x)>>63, (get_pe_e(*x)<<1)>>(64-idx->uID_bits), get_pe_e(*x) & idx->pos_mode); + uint64_t i; + for (i = 0; i < x->occ1; i++) { - if (k - l > 1) radix_sort_pe_hit_an2(hits->a.a + l, hits->a.a + k); - cur = (uint64_t)-1; + fprintf(stderr, "###-1-rev: %lu, uID: %lu, ref_p: %lu\n", + x->a[i]>>63, (x->a[i]<<1)>>(64-idx->uID_bits), x->a[i] & idx->pos_mode); + } + + for (i = 0; i < x->occ2; i++) + { + fprintf(stderr, "###-2-rev: %lu, uID: %lu, ref_p: %lu\n", + x->a[i+x->occ1]>>63, (x->a[i+x->occ1]<<1)>>(64-idx->uID_bits), x->a[i+x->occ1] & idx->pos_mode); + } + } + **/ + /***************************for debug******************************/ + for (k = 1, l = 0; k <= hits->n; ++k) + { + if (k == hits->n || get_pe_s(hits->a[k]) != get_pe_s(hits->a[l])) + { + if (k - l > 1) radix_sort_pe_an2(hits->a + l, hits->a + k); + ////fprintf(stderr, "\nl: %lu, k: %lu, %s\n", l, k, k - l > 1? "Found":"NONE"); + + cur = (uint64_t)-1; while (l < k) { - if(hits->a.a[l].e != cur) + if(get_pe_e(hits->a[l]) != cur) { - cur = hits->a.a[l].e; - hits->a.a[m++] = hits->a.a[l]; + cur = get_pe_e(hits->a[l]); + if(m != l) swap_pe_hit_hap(&hits->a[m], &hits->a[l]); + m++; } l++; } l = k; } } - hits->a.n = m; - fprintf(stderr, "[M::%s::%.3f] ==> Dedup\n", __func__, yak_realtime()-index_time); + + for (k = m; k < hits->n; k++) + { + hits->a[k].id = (uint64_t)-1; + hits->a[k].occ1 = hits->a[k].occ2 = 0; + free(hits->a[k].a); hits->a[k].a = NULL; + } + + radix_sort_pe_occ_t(hits->a, hits->a + m); + for (k = 0, hits->n_u = 0; k < m; k++) + { + if(hits->a[k].occ1 == 1 && hits->a[k].occ2 == 1) hits->n_u++; + } + + + fprintf(stderr, "[M::%s::%.3f] ==> Dedup (# dup: %lu, # non-dup: %lu, # non-dup-unique: %lu)\n", + __func__, yak_realtime()-index_time, (uint64_t)(hits->n - m), m, hits->n_u); + hits->n = m; } +void int_kvec_pe_hit_hap(kvec_pe_hit_hap* x) +{ + x->m = x->n = x->n_u = 0; + x->a = NULL; +} + +void destory_kvec_pe_hit_hap(kvec_pe_hit_hap* x) +{ + uint64_t k; + for (k = 0; k < x->n; k++) + { + x->a[k].id = (uint64_t)-1; + x->a[k].occ1 = x->a[k].occ2 = 0; + free(x->a[k].a); x->a[k].a = NULL; + } + free(x->a); +} + + void sort_hits(kvec_pe_hit* hits) { double index_time = yak_realtime(); @@ -1849,7 +2222,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, hc_links* link) if(ug->g->seq[v>>1].del) continue; if(asg_arc_n(ug->g, v) < 2) continue; if((bub->index[v]&(uint32_t)3) != 0) continue; - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL)) + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL)) { //beg is v, end is b.S.a[0] //note b.b include end, does not include beg @@ -1870,7 +2243,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, hc_links* link) for (v = 0; v < n_vtx; ++v) { if((bub->index[v]&(uint32_t)3) !=2) continue; - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL)) + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL)) { //note b.b include end, does not include beg i = b.b.n + 1; @@ -1903,7 +2276,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, hc_links* link) if((bub->num.a[k]>>31) == 0) bub->s_bub++; v = (bub->num.a[k]<<1)>>1; bub->num.a[k] = bub->list.n; - if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL)) + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL)) { kv_push(uint64_t, bub->pathLen, pathLen); //beg is v, end is b.S.a[0] @@ -2026,7 +2399,7 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, hc_links* link) -void print_bubbles(ma_ug_t* ug, bubble_type* bub, kvec_pe_hit* hits, hc_links* link, ha_ug_index* idx) +void print_bubbles(ma_ug_t* ug, bubble_type* bub, kvec_pe_hit_hap* hits, hc_links* link, ha_ug_index* idx) { uint64_t tLen, t_utg, i, k; uint32_t beg, sink, n, *a; @@ -2070,10 +2443,10 @@ void print_bubbles(ma_ug_t* ug, bubble_type* bub, kvec_pe_hit* hits, hc_links* l uint64_t s_uid, e_uid, shif = 64 - idx->uID_bits; if(hits) { - for (k = 0; k < hits->a.n; ++k) + for (k = 0; k < hits->n_u; ++k) { - s_uid = ((hits->a.a[k].s<<1)>>shif); - e_uid = ((hits->a.a[k].e<<1)>>shif); + s_uid = ((get_pe_s(hits->a[k])<<1)>>shif); + e_uid = ((get_pe_e(hits->a[k])<<1)>>shif); if(bub->index[s_uid] == (uint32_t)-1 || bub->index[e_uid] == (uint32_t)-1) continue; if(IF_BUB(s_uid, *bub) && IF_BUB(e_uid, *bub)) { @@ -2803,14 +3176,14 @@ void destory_MT(MT* M) kv_destroy(M->matrix); } -void collect_hc_links(const ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, MT* M) +void collect_hc_links(const ha_ug_index* idx, kvec_pe_hit_hap* hits, hc_links* link, bubble_type* bub, MT* M) { double index_time = yak_realtime(); uint64_t k, i, shif = 64 - idx->uID_bits, beg, end, t_d; - for (k = 0; k < hits->a.n; ++k) + for (k = 0; k < hits->n_u; ++k) { - beg = ((hits->a.a[k].s<<1)>>shif); - end = ((hits->a.a[k].e<<1)>>shif); + beg = ((get_pe_s(hits->a[k])<<1)>>shif); + end = ((get_pe_e(hits->a[k])<<1)>>shif); if(beg == end) continue; if(IF_HOM(beg, *bub)) continue; @@ -3092,22 +3465,179 @@ int load_hc_links(hc_links* link, const char *fn) } -void write_hc_hits(kvec_pe_hit* hits, const char *fn) +void write_hc_hits(kvec_pe_hit_hap* hits, const char *fn) { char *buf = (char*)calloc(strlen(fn) + 25, 1); sprintf(buf, "%s.hic.lk.bin", fn); FILE* fp = fopen(buf, "w"); - fwrite(&hits->a.n, sizeof(hits->a.n), 1, fp); - fwrite(hits->a.a, sizeof(pe_hit), hits->a.n, fp); - + uint64_t k; + fwrite(&hits->n_u, sizeof(hits->n_u), 1, fp); + fwrite(&hits->n, sizeof(hits->n), 1, fp); + for (k = 0; k < hits->n; k++) + { + fwrite(&hits->a[k].id, sizeof(hits->a[k].id), 1, fp); + fwrite(&hits->a[k].occ1, sizeof(hits->a[k].occ1), 1, fp); + fwrite(&hits->a[k].occ2, sizeof(hits->a[k].occ2), 1, fp); + fwrite(hits->a[k].a, sizeof(uint64_t), hits->a[k].occ1 + hits->a[k].occ2, fp); + } + fclose(fp); free(buf); } -int load_hc_hits(kvec_pe_hit* hits, const char *fn) +void write_hc_hits_v14(kvec_pe_hit_hap* i_hits, const char *fn) +{ + char *buf = (char*)calloc(strlen(fn) + 25, 1); + sprintf(buf, "%s.v14.hic.lk.bin", fn); + FILE* fp = fopen(buf, "w"); + kvec_pe_hit hits; + kv_init(hits.a); + uint64_t i, m_u = (uint64_t)-1, m_m = (uint64_t)-1; + pe_hit* p = NULL; + for (i = 0; i < i_hits->n; i++) + { + if(i_hits->a[i].occ1 == 1 && i_hits->a[i].occ2 == 1) + { + kv_pushp(pe_hit, hits.a, &p); + p->id = i_hits->a[i].id; + p->s = i_hits->a[i].a[0]; + p->e = i_hits->a[i].a[1]; + m_u = i; + } + else + { + if(m_m == (uint64_t)-1) m_m = i; + } + } + fprintf(stderr, "m_u: %lu, m_m: %lu, n_u: %lu\n", m_u, m_m, i_hits->n_u); + + fwrite(&hits.a.n, sizeof(hits.a.n), 1, fp); + fwrite(hits.a.a, sizeof(pe_hit), hits.a.n, fp); + + kv_destroy(hits.a); + fclose(fp); + free(buf); + exit(1); +} + +#define pe_hit_hap_id_key(x) ((x).id) +KRADIX_SORT_INIT(pe_hit_hap_id, pe_hit_hap, pe_hit_hap_id_key, member_size(pe_hit_hap, id)) + +#define pe_hit_id_key(x) ((x).id) +KRADIX_SORT_INIT(pe_hit_id, pe_hit, pe_hit_id_key, member_size(pe_hit, id)) + +void debug_hc_hits_v14(kvec_pe_hit_hap* i_hits, const char *fn, const ha_ug_index* idx) { uint64_t flag = 0; + char *buf = (char*)calloc(strlen(fn) + 25, 1); + sprintf(buf, "%s.v14.hic.lk.bin", fn); + kvec_pe_hit hits; + kv_init(hits.a); + FILE* fp = NULL; + fp = fopen(buf, "r"); + + kv_init(hits.a); + flag += fread(&hits.a.n, sizeof(hits.a.n), 1, fp); + hits.a.m = hits.a.n; MALLOC(hits.a.a, hits.a.n); + flag += fread(hits.a.a, sizeof(pe_hit), hits.a.n, fp); + + radix_sort_pe_hit_id(hits.a.a, hits.a.a + hits.a.n); + radix_sort_pe_hit_hap_id(i_hits->a, i_hits->a + i_hits->n_u); + + fprintf(stderr, "i_hits->n_u: %lu, hits.a.n: %lu\n", (uint64_t)i_hits->n_u, (uint64_t)hits.a.n); + + uint64_t i, k; + uint64_t i_beg_utg, i_beg_pos, i_beg_rev; + uint64_t i_end_utg, i_end_pos, i_end_rev; + uint64_t k_beg_utg, k_beg_pos, k_beg_rev; + uint64_t k_end_utg, k_end_pos, k_end_rev; + uint64_t i_id, k_id; + uint64_t same_occ = 0, diff_occ = 0, miss_occ = 0; + for (i = 0, k = 0; i < i_hits->n_u; i++) + { + i_beg_rev = get_pe_s(i_hits->a[i])>>63; + i_beg_utg = ((get_pe_s(i_hits->a[i])<<1)>>(64 - idx->uID_bits)); + i_beg_pos = get_pe_s(i_hits->a[i]) & idx->pos_mode; + + i_end_rev = get_pe_e(i_hits->a[i])>>63; + i_end_utg = ((get_pe_e(i_hits->a[i])<<1)>>(64 - idx->uID_bits)); + i_end_pos = get_pe_e(i_hits->a[i]) & idx->pos_mode; + + i_id = i_hits->a[i].id; + for (; k < hits.a.n; k++) + { + k_beg_rev = hits.a.a[k].s>>63; + k_beg_utg = ((hits.a.a[k].s<<1)>>(64 - idx->uID_bits)); + k_beg_pos = hits.a.a[k].s & idx->pos_mode; + + k_end_rev = hits.a.a[k].e>>63; + k_end_utg = ((hits.a.a[k].e<<1)>>(64 - idx->uID_bits)); + k_end_pos = hits.a.a[k].e & idx->pos_mode; + + k_id = hits.a.a[k].id; + + if(k_id > i_id) + { + miss_occ++; + fprintf(stderr, "\n[MISS]rid=%lu\n", i_id); + fprintf(stderr, "********v0.15********\n"); + fprintf(stderr, "beg_rev: %lu, beg_utg: %lu, beg_pos: %lu\n", + i_beg_rev, i_beg_utg, i_beg_pos); + fprintf(stderr, "end_rev: %lu, end_utg: %lu, end_pos: %lu\n", + i_end_rev, i_end_utg, i_end_pos); + break; + } + + if(k_id == i_id) + { + if(get_pe_s(i_hits->a[i]) == hits.a.a[k].s && get_pe_e(i_hits->a[i]) == hits.a.a[k].e) + { + same_occ++; + // fprintf(stderr, "\n[SAME]rid=%lu\n", i_id); + // fprintf(stderr, "********v0.15********\n"); + // fprintf(stderr, "beg_rev: %lu, beg_utg: %lu, beg_pos: %lu\n", + // i_beg_rev, i_beg_utg, i_beg_pos); + // fprintf(stderr, "end_rev: %lu, end_utg: %lu, end_pos: %lu\n", + // i_end_rev, i_end_utg, i_end_pos); + // fprintf(stderr, "********v0.14********\n"); + // fprintf(stderr, "beg_rev: %lu, beg_utg: %lu, beg_pos: %lu\n", + // k_beg_rev, k_beg_utg, k_beg_pos); + // fprintf(stderr, "end_rev: %lu, end_utg: %lu, end_pos: %lu\n", + // k_end_rev, k_end_utg, k_end_pos); + } + else + { + diff_occ++; + fprintf(stderr, "\n[DIFF]rid=%lu\n", i_id); + fprintf(stderr, "********v0.15********\n"); + fprintf(stderr, "beg_rev: %lu, beg_utg: %lu, beg_pos: %lu\n", + i_beg_rev, i_beg_utg, i_beg_pos); + fprintf(stderr, "end_rev: %lu, end_utg: %lu, end_pos: %lu\n", + i_end_rev, i_end_utg, i_end_pos); + fprintf(stderr, "********v0.14********\n"); + fprintf(stderr, "beg_rev: %lu, beg_utg: %lu, beg_pos: %lu\n", + k_beg_rev, k_beg_utg, k_beg_pos); + fprintf(stderr, "end_rev: %lu, end_utg: %lu, end_pos: %lu\n", + k_end_rev, k_end_utg, k_end_pos); + } + break; + } + } + } + + fprintf(stderr, "same_occ: %lu, diff_occ: %lu, miss_occ: %lu", same_occ, diff_occ, miss_occ); + + + kv_destroy(hits.a); + fclose(fp); + free(buf); + exit(1); +} + +int load_hc_hits(kvec_pe_hit_hap* hits, const char *fn) +{ + uint64_t flag = 0, k; char *buf = (char*)calloc(strlen(fn) + 25, 1); sprintf(buf, "%s.hic.lk.bin", fn); @@ -3115,10 +3645,19 @@ int load_hc_hits(kvec_pe_hit* hits, const char *fn) fp = fopen(buf, "r"); if(!fp) return 0; - kv_init(hits->a); - flag += fread(&hits->a.n, sizeof(hits->a.n), 1, fp); - hits->a.m = hits->a.n; MALLOC(hits->a.a, hits->a.n); - flag += fread(hits->a.a, sizeof(pe_hit), hits->a.n, fp); + kv_init(*hits); + flag += fread(&hits->n_u, sizeof(hits->n_u), 1, fp); + flag += fread(&hits->n, sizeof(hits->n), 1, fp); + hits->m = hits->n; MALLOC(hits->a, hits->n); + + for (k = 0; k < hits->n; k++) + { + flag += fread(&hits->a[k].id, sizeof(hits->a[k].id), 1, fp); + flag += fread(&hits->a[k].occ1, sizeof(hits->a[k].occ1), 1, fp); + flag += fread(&hits->a[k].occ2, sizeof(hits->a[k].occ2), 1, fp); + MALLOC(hits->a[k].a, hits->a[k].occ1 + hits->a[k].occ2); + flag += fread(hits->a[k].a, sizeof(uint64_t), hits->a[k].occ1 + hits->a[k].occ2, fp); + } fclose(fp); free(buf); @@ -4414,12 +4953,12 @@ G_partition* clean_bubbles(hc_links* link, bubble_type* bub, min_cut_t* m, const -uint64_t get_hic_distance(pe_hit* hit, hc_links* link, const ha_ug_index* idx) +uint64_t get_hic_distance(pe_hit_hap* hit, hc_links* link, const ha_ug_index* idx) { uint64_t s_uid, s_dir, e_uid, e_dir, u_dis, k; long long s_pos, e_pos; - s_uid = ((hit->s<<1)>>(64 - idx->uID_bits)); s_pos = hit->s & idx->pos_mode; - e_uid = ((hit->e<<1)>>(64 - idx->uID_bits)); e_pos = hit->e & idx->pos_mode; + s_uid = ((get_pe_s(*hit)<<1)>>(64 - idx->uID_bits)); s_pos = get_pe_s(*hit) & idx->pos_mode; + e_uid = ((get_pe_e(*hit)<<1)>>(64 - idx->uID_bits)); e_pos = get_pe_e(*hit) & idx->pos_mode; if(s_uid == e_uid) return MAX(s_pos, e_pos) - MIN(s_pos, e_pos); hc_linkeage* t = &(link->a.a[s_uid]); for (k = 0; k < t->e.n; k++) @@ -4677,7 +5216,7 @@ void LeastSquare_advance(trans_idx* dis, ha_ug_index* idx, uint64_t med) } -void weight_edges(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub) +void weight_edges(ha_ug_index* idx, kvec_pe_hit_hap* hits, hc_links* link, bubble_type* bub) { uint64_t k, i, shif = 64 - idx->uID_bits, beg, end, t_d; hc_edge *e1 = NULL, *e2 = NULL; @@ -4692,16 +5231,16 @@ void weight_edges(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_ty } } - for (k = 0; k < hits->a.n; ++k) + for (k = 0; k < hits->n_u; ++k) { - beg = ((hits->a.a[k].s<<1)>>shif); - end = ((hits->a.a[k].e<<1)>>shif); + beg = ((get_pe_s(hits->a[k])<<1)>>shif); + end = ((get_pe_e(hits->a[k])<<1)>>shif); if(beg == end) continue; if(IF_HOM(beg, *bub)) continue; if(IF_HOM(end, *bub)) continue; - t_d = get_hic_distance(&(hits->a.a[k]), link, idx); + t_d = get_hic_distance(&(hits->a[k]), link, idx); if(t_d == (uint64_t)-1) continue; e1 = get_hc_edge(link, beg, end, 0); @@ -4718,7 +5257,7 @@ void weight_edges(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_ty } -void weight_edges_advance(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, trans_idx* dis) +void weight_edges_advance(ha_ug_index* idx, kvec_pe_hit_hap* hits, hc_links* link, bubble_type* bub, trans_idx* dis) { uint64_t k, i, shif = 64 - idx->uID_bits, beg, end, t_d; hc_edge *e1 = NULL, *e2 = NULL; @@ -4733,16 +5272,16 @@ void weight_edges_advance(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, b } } - for (k = 0; k < hits->a.n; ++k) + for (k = 0; k < hits->n_u; ++k) { - beg = ((hits->a.a[k].s<<1)>>shif); - end = ((hits->a.a[k].e<<1)>>shif); + beg = ((get_pe_s(hits->a[k])<<1)>>shif); + end = ((get_pe_e(hits->a[k])<<1)>>shif); if(beg == end) continue; if(IF_HOM(beg, *bub)) continue; if(IF_HOM(end, *bub)) continue; - t_d = get_hic_distance(&(hits->a.a[k]), link, idx); + t_d = get_hic_distance(&(hits->a[k]), link, idx); if(t_d == (uint64_t)-1) continue; e1 = get_hc_edge(link, beg, end, 0); @@ -7747,7 +8286,7 @@ void get_forward_distance(uint32_t src, uint32_t dest, asg_t *sg, hc_links* link -int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, MT* M, H_partition* hap, trans_idx* dis) +int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit_hap* hits, hc_links* link, bubble_type* bub, MT* M, H_partition* hap, trans_idx* dis) { kvec_t(uint64_t) buf, buf_idx; kv_init(buf); @@ -7766,16 +8305,16 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, buf.n = 0; - for (k = 0; k < hits->a.n; ++k) + for (k = 0; k < hits->n_u; ++k) { - beg = ((hits->a.a[k].s<<1)>>(64 - idx->uID_bits)); - end = ((hits->a.a[k].e<<1)>>(64 - idx->uID_bits)); + beg = ((get_pe_s(hits->a[k])<<1)>>(64 - idx->uID_bits)); + end = ((get_pe_e(hits->a[k])<<1)>>(64 - idx->uID_bits)); if(IF_HOM(beg, *bub)) continue; if(IF_HOM(end, *bub)) continue; - t_d = get_hic_distance(&(hits->a.a[k]), link, idx); + t_d = get_hic_distance(&(hits->a[k]), link, idx); if(t_d == (uint64_t)-1) continue; if(beg == end) { @@ -7809,7 +8348,6 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, if((buf.a[k]&1) == 1) f_idx = k; } buf.n = MIN(r_idx, f_idx); - for (k = 0; k < buf.n; k++) { if((buf.a[k]&1) == 1) @@ -7817,7 +8355,6 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, kv_push(uint64_t, buf_idx, buf.a[k]>>1); } } - uint64_t cutoff = buf_idx.n * 0.9, t = buf_idx.n * 0.005, pre, step; k = 0; if(cutoff >= t) k = cutoff - t; @@ -7829,14 +8366,12 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, pre = buf_idx.a[k]; i++; } - if(t_d == 0 || i == 0 || t == 0) { kv_destroy(buf); kv_destroy(buf_idx); return 0; } - step = (t_d/i)*20; if(step == 0) @@ -7845,7 +8380,6 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, kv_destroy(buf_idx); return 0; } - trans_p_t* p = NULL; dis->n = 0; uint64_t step_s = 0, step_e = step; @@ -7871,8 +8405,9 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, cnt[0] = cnt[1] = 0; } } + // fprintf(stderr, "-k: %lu, buf.n: %lu, buf.a[k]: %lu, step_s: %lu, step_e: %lu\n", + // k, (uint64_t)buf.n, (buf.a[k]>>1), step_s, step_e); } - if(cnt[0] > 0 || cnt[1] > 0) { kv_pushp(trans_p_t, *dis, &p); @@ -7882,7 +8417,6 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, p->cnt_1 = cnt[1]; } - uint64_t smooth_step = 20, k_i, cnt_0; if(dis->n > 0) med = dis->a[dis->n-1].end; for (k = 0; k+smooth_step < dis->n; k++) @@ -7898,7 +8432,6 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, break; } } - long long b_k = 0, b_i = 0, b_j, pass = 0; ///for (b_k = b_i = 0; b_k < (long long)dis->n; b_k++) @@ -7964,7 +8497,6 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, pass = 0; break; } - dis->n = b_i; if(dis->n == 0 || pass == 0) { @@ -7980,9 +8512,7 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, // dis->a[i].beg, dis->a[i].end, dis->a[i].cnt_0, dis->a[i].cnt_1, (double)(dis->a[i].cnt_1)/(double)(dis->a[i].cnt_1 + dis->a[i].cnt_0)); // } - LeastSquare_advance(dis, idx, med); - // fprintf(stderr, "idx->a: %f, idx->b: %f, idx->frac: %f, med: %lu\n", // (double)idx->a, (double)idx->b, (double)idx->frac, med); @@ -7991,7 +8521,6 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, kv_destroy(buf); kv_destroy(buf_idx); - if(idx->a < 0) idx->a = 0; if(idx->a == 0) { @@ -8002,25 +8531,24 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, idx->b = ((double)(dis->a[dis->n-1].cnt_1))/((double)(dis->a[dis->n-1].cnt_0 + dis->a[dis->n-1].cnt_1)); } - // fprintf(stderr, "idx->a: %f, idx->b: %f, idx->frac: %f, med: %lu\n", // (double)idx->a, (double)idx->b, (double)idx->frac, med); return 1; } -void init_hic_p(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, +void init_hic_p(ha_ug_index* idx, kvec_pe_hit_hap* hits, hc_links* link, bubble_type* bub, kvec_hc_edge* back_hc_edge, MT* M, H_partition* hap, uint32_t ignore_dis) { uint64_t k, i, m, uID, is_comples_weight = 0; trans_idx dis; kv_init(dis); - + if(bub->round_id > 0 && ignore_dis == 0) { is_comples_weight = get_trans_rate_function(idx, hits, link, bub, M, hap, &dis); } - + hc_edge *e = NULL; for (i = 0; i < link->a.n; i++) @@ -8082,10 +8610,8 @@ kvec_hc_edge* back_hc_edge, MT* M, H_partition* hap, uint32_t ignore_dis) } link->a.a[i].e.n = m; } - weight_edges_advance(idx, hits, link, bub, is_comples_weight == 1? &dis : NULL); - for (i = 0; i < link->a.n; i++) { for (k = 0; k < link->a.a[i].e.n; k++) @@ -11340,7 +11866,7 @@ void print_bubble_chain(bubble_type* bub) } } -void init_contig_H_partition(bubble_type* bub, ha_ug_index* idx, kvec_pe_hit* hits, H_partition* hap) +void init_contig_H_partition(bubble_type* bub, ha_ug_index* idx, kvec_pe_hit_hap* hits, H_partition* hap) { uint32_t i, k_i, k_j, uID, *a = NULL, n, *h0, h0_n, *h1, h1_n; destory_G_partition(&(hap->group_g_p)); memset(&(hap->group_g_p), 0, sizeof(G_partition)); @@ -11416,15 +11942,15 @@ void init_contig_H_partition(bubble_type* bub, ha_ug_index* idx, kvec_pe_hit* hi label_unitigs(&(hap->group_g_p), idx->ug); } -void cluster_contigs(bubble_type* bub, ha_ug_index* idx, kvec_pe_hit* hits, MT* M, H_partition* hap) +void cluster_contigs(bubble_type* bub, ha_ug_index* idx, kvec_pe_hit_hap* hits, MT* M, H_partition* hap) { uint64_t k, i, shif = 64 - idx->uID_bits, beg, end, t_d; hc_links* link = idx->link; for (i = 0; i < link->a.n; i++) link->a.a[i].e.n = 0; - for (k = 0; k < hits->a.n; ++k) + for (k = 0; k < hits->n_u; ++k) { - beg = ((hits->a.a[k].s<<1)>>shif); - end = ((hits->a.a[k].e<<1)>>shif); + beg = ((get_pe_s(hits->a[k])<<1)>>shif); + end = ((get_pe_e(hits->a[k])<<1)>>shif); if(beg == end) continue; if(IF_HOM(beg, *bub)) continue; @@ -11468,6 +11994,7 @@ void reset_H_partition(H_partition* hap, uint32_t is_init) int alignment_worker_pipeline(sldat_t* sl, const enzyme *fn1, const enzyme *fn2) { + double index_time = yak_realtime(); int i; for (i = 0; i < fn1->n && i < fn2->n; i++) { @@ -11479,21 +12006,228 @@ int alignment_worker_pipeline(sldat_t* sl, const enzyme *fn1, const enzyme *fn2) kt_pipeline(3, worker_pipeline, sl, 3); - kseq_destroy(sl->ks1); kseq_destroy(sl->ks2); gzclose(fp1); gzclose(fp2); } + fprintf(stderr, "[M::%s::%.3f] ==> Qualification\n", __func__, yak_realtime()-index_time); - - dedup_hits(&(sl->hits)); - - ///fprintf(stderr, "-sl->hits.a.n: %u\n", (uint32_t)sl->hits.a.n); - + dedup_hits(&(sl->hits), sl->idx); return 1; } +/** +typedef struct{ + FILE* fp; + kvec_t(char) buf; + kvec_t(char) name; + kvec_t_u64_warp pos; +}pe_aln_t; + +int init_pe_aln_t(pe_aln_t* x, const char* aln) +{ + memset(x, 0, sizeof(*x)); + if(!strcmp(aln,"-")) x->fp = stdin; + else if ((x->fp = fopen(aln, "r")) == 0) return 0; + kv_malloc(x->buf, 10); x->buf.n = 0; + kv_malloc(x->name, 10); x->name.n = 0; + kv_init(x->pos.a); + return 1; +} + +void destory_pe_aln_t(pe_aln_t* x) +{ + fclose(x->fp); + kv_destroy(x->buf); + kv_destroy(x->name); + kv_destroy(x->pos.a); +} + +char* get_alnLine(pe_aln_t* x) +{ + uint64_t len; + uint64_t b_size = x->buf.m; + char* b = x->buf.a; + while (fgets(b, b_size, x->fp) != NULL) + { + len = strlen(x->buf.a); + if(x->buf.a[len - 1] == '\n') + { + x->buf.a[len - 1] = '\0'; + return x->buf.a; + } + kv_resize(char, x->buf, x->buf.m<<1); + b = x->buf.a + len; b_size = x->buf.m - len; + } + return NULL; +} + +uint64_t get_read_id_by_name(char* name, uint64_t name_len, reads_t* r1) +{ + uint64_t size = r1->r_Len.n - 1, r_len; + uint64_t end_idx = r1->idx; + char* r_char = NULL; + while(1) + { + r_len = r1->name_Len.a[r1->idx + 1] - r1->name_Len.a[r1->idx]; + r_char = r1->name.a + r1->name_Len.a[r1->idx]; + if(name_len == r_len && memcmp(r_char, name, r_len) == 0) return r1->idx; + r1->idx++; + if(r1->idx >= size) r1->idx = 0; + if(r1->idx == end_idx) break; + } + return (uint64_t)-1; +} +uint64_t get_utg_id_by_name(char* u_name) +{ + uint64_t i, len = strlen(u_name), id; + char c = u_name[len - 1]; + u_name[len - 1] = '\0'; + for (i = 3; i < len; i++) + { + if(u_name[i] != '0') break; + } + id = atoi(u_name + i); + u_name[len - 1] = c; + return id; +} + +uint64_t adjust_pos(uint64_t pos, uint64_t rev, char* cigar) +{ + long long i, occ = strlen(cigar); + + if(rev == 0) + { + for (i = 0; i < occ; i++) + { + if(cigar[i] < '0' || cigar[i] > '9') + { + break; + } + } + + if(cigar[i] == 'S') + { + cigar[i] = '\0'; + pos = pos + atoll(cigar); + cigar[i] = 'S'; + } + } + else + { + if(cigar[occ-1] == 'S') + { + cigar[occ-1] = '\0'; + for (i = occ-2; i >= 0; i--) + { + if(cigar[i] < '0' || cigar[i] > '9') + { + break; + } + } + pos = pos - atoll(cigar+i+1); + cigar[occ-1] = 'S'; + } + } + + return pos; +} + +uint64_t parse_sam(char *x, char** name, uint64_t* flag, uint64_t* uid, kvec_t_u64_warp* pos) +{ + uint64_t p_pos, p_err, n_len; + p->a.n = 0; + char *t = NULL; + + t = strtok (a, "\t\0");///name + n_len = strlen(t); + kv_resize(char, x->name, n_len+1); + memcpy(x->name, t, n_len+1); + (*name) = x->name; + + (*flag) = atoll(strtok (NULL, "\t\0"));//flag + if(!((*flag)&1) || ((*flag)&4) || ((*flag)&256) || ((*flag)&2048)) return 0; + + (*uid) = get_utg_id_by_name(strtok(NULL, "\t\0")); ///utg name + + p_pos = atoll(strtok(NULL, "\t\0")) - 1;//primary pos + + strtok(NULL, "\t\0");///MAPQ + + p_pos = adjust_pos(p_pos, !!((*flag)&16), strtok(NULL, "\t\0")); ///cigar + + strtok(NULL, "\t\0"); + strtok(NULL, "\t\0"); + strtok(NULL, "\t\0"); + strtok(NULL, "\t\0"); + strtok(NULL, "\t\0"); + + p_err = atoll(strtok(NULL, "\t\0") + 5); //NM:i: + + t = strtok(NULL, "\t\0"); + while (t != NULL) + { + n_len = strlen(t); + if(n_len > 5 && t[0] == 'X' && t[1] == 'A' && t[2] == ':' && t[3] == 'Z' && t[4] == ':') + { + break; + } + t = strtok(NULL, "\t\0"); + } +} + +uint64_t get_sam(pe_aln_t* x) +{ + char *a = x->buf.a, *t = NULL; + uint64_t n_len; + while (1) + { + if(x->buf.n == 0) + { + a = get_alnLine(x); + if(a == NULL) break; + x->buf.n = 1; + } + + ///parse_sam(char *x, char** name, uint64_t* flag, uint64_t* uid, kvec_t_u64_warp* pos) + + + } + + return 0; +} + + + +int debug_hits_sam(ha_ug_index* idx, kvec_pe_hit_hap* hits, const enzyme *fn1, const enzyme *fn2, +const char* aln) +{ + uint64_t k, id, uid, shif = 64 - idx->uID_bits, b_size = 100000, flag; + reads_t r1; + load_reads(&r1, fn1, fn2); r1.idx = 0; + pe_aln_t p; + init_pe_aln_t(&p, aln); + + + + + while (get_alnLine(p) != NULL) + { + str = strtok (buffer, "\t");///name + id = get_read_id_by_name(str, strlen(str), &r1); + if(id == (uint64_t)-1) fprintf(stderr, "ERROR\n"); + flag = atoi(strtok (NULL, "\t"));//flag + if(!(flag&1) || (flag&4) || (flag&256) || (flag&2048)) continue; + uid = get_utg_id_by_name(strtok (NULL, "\t")); ///utg name + + } + + destory_reads(&r1); + destory_pe_aln_t(&p); +} +**/ + int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) { double index_time = yak_realtime(); @@ -11506,7 +12240,8 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) sl.n_thread = asm_opt.thread_num; sl.total_base = sl.total_pair = 0; idx->hap_cnt = asm_opt.hap_occ; - kv_init(sl.hits.a); + int_kvec_pe_hit_hap(&sl.hits); + if(!load_hc_hits(&sl.hits, asm_opt.output_file_name)) { @@ -11528,6 +12263,9 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) write_hc_hits(&sl.hits, asm_opt.output_file_name); } + ///debug_hc_hits_v14(&sl.hits, asm_opt.output_file_name, sl.idx); + ////dedup_hits(&(sl.hits), sl.idx); + ///write_hc_hits_v14(&sl.hits, asm_opt.output_file_name); ///fprintf(stderr, "u.n: %d, uID_bits: %lu, pos_bits: %lu, sl.hits.a.n: %u\n", (uint32_t)idx->ug->u.n, idx->uID_bits, idx->pos_bits, (uint32_t)sl.hits.a.n); H_partition hap; @@ -11550,7 +12288,6 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) init_contig_partition(&hap, idx, &bub); phasing_improvement(&hap, &(hap.g_p), idx, &bub); label_unitigs(&(hap.g_p), idx->ug); - ///print_hc_links(idx->link, 0, &hap); } @@ -11581,13 +12318,14 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) destory_contig_partition(&hap); kv_destroy(back_hc_edge.a); + destory_kvec_pe_hit_hap(&sl.hits); return 1; /*******************************for debug************************************/ // destory_reads(&R1); // destory_reads(&R2); /*******************************for debug************************************/ - print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, idx->link, idx); + print_bubbles(idx->ug, &bub, sl.hits.n?&sl.hits:NULL, idx->link, idx); collect_hc_reverse_links(idx->link, idx->ug, &bub); normalize_hc_links(idx->link); /*******************************for debug************************************/ @@ -11600,7 +12338,6 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) destory_min_cut_t(cut); free(cut); destory_G_partition(gp); free(gp); - kv_destroy(sl.hits.a); destory_bubbles(&bub); fprintf(stderr, "[M::%s::%.3f] processed %lu pairs; %lu bases\n", __func__, yak_realtime()-index_time, sl.total_pair, sl.total_base); @@ -11755,13 +12492,13 @@ void print_bench_idx(bench_idx* idx, ma_ug_t *ug) } -uint64_t get_hic_distance_bench(pe_hit* hit, hc_links* link, bench_idx* idx, ma_ug_t *ug, uint64_t* is_trans) +uint64_t get_hic_distance_bench(pe_hit_hap* hit, hc_links* link, bench_idx* idx, ma_ug_t *ug, uint64_t* is_trans) { (*is_trans) = (uint64_t)-1; uint64_t s_uid, e_uid; long long s_pos, e_pos; - s_uid = ((hit->s<<1)>>(64 - idx->uID_bits)); s_pos = hit->s & idx->pos_mode; - e_uid = ((hit->e<<1)>>(64 - idx->uID_bits)); e_pos = hit->e & idx->pos_mode; + s_uid = ((get_pe_s(*hit)<<1)>>(64 - idx->uID_bits)); s_pos = get_pe_s(*hit) & idx->pos_mode; + e_uid = ((get_pe_e(*hit)<<1)>>(64 - idx->uID_bits)); e_pos = get_pe_e(*hit) & idx->pos_mode; if(s_uid == e_uid) { (*is_trans) = 0; @@ -11829,14 +12566,14 @@ void init_bench_idx(bench_idx* idx, asg_t* read_g, ma_ug_t *ug) } } -void evaluate_bench_idx(bench_idx* idx, kvec_pe_hit* hits, ma_ug_t *ug) +void evaluate_bench_idx(bench_idx* idx, kvec_pe_hit_hap* hits, ma_ug_t *ug) { uint64_t k, distance, is_trans, trans[2]; kvec_t(uint64_t) buf; kv_init(buf); - for (k = trans[0] = trans[1] = 0; k < hits->a.n; ++k) + for (k = trans[0] = trans[1] = 0; k < hits->n_u; ++k) { - distance = get_hic_distance_bench(&(hits->a.a[k]), &(idx->link), idx, ug, &is_trans); + distance = get_hic_distance_bench(&(hits->a[k]), &(idx->link), idx, ug, &is_trans); if(is_trans != (uint64_t)-1) trans[is_trans]++; if(distance == (uint64_t)-1 || is_trans == (uint64_t)-1) continue; distance = (distance << 1) + is_trans; @@ -11904,7 +12641,7 @@ int hic_short_align_bench(const enzyme *fn1, const enzyme *fn2, const char *outp sl.n_thread = asm_opt.thread_num; sl.total_base = sl.total_pair = 0; idx->hap_cnt = asm_opt.hap_occ; - kv_init(sl.hits.a); + int_kvec_pe_hit_hap(&sl.hits); fprintf(stderr, "u.n: %d, uID_bits: %lu, pos_bits: %lu\n", (uint32_t)idx->ug->u.n, idx->uID_bits, idx->pos_bits); if(!load_hc_hits(&sl.hits, output_file_name)) @@ -11920,7 +12657,7 @@ int hic_short_align_bench(const enzyme *fn1, const enzyme *fn2, const char *outp evaluate_bench_idx(&bench, &sl.hits, idx->ug); destory_bench_idx(&bench); - kv_destroy(sl.hits.a); + destory_kvec_pe_hit_hap(&sl.hits); fprintf(stderr, "[M::%s::%.3f] processed %lu pairs; %lu bases\n", __func__, yak_realtime()-index_time, sl.total_pair, sl.total_base); return 1; } diff --git a/hifiasm.1 b/hifiasm.1 index 42bb851..245a7b1 100644 --- a/hifiasm.1 +++ b/hifiasm.1 @@ -286,7 +286,8 @@ 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. In default, [2] for non-trio assembly, [0] for trio assembly. +2 to purge all types of haplotigs, 3 to purge all types of haplotigs in most aggressive way. +In default, [2] for non-trio assembly, [0] for trio assembly. For trio assembly, only level 0 and level 1 are allowed. .TP