diff --git a/Assembly.cpp b/Assembly.cpp index 42bf649..bc719ea 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1784,8 +1784,9 @@ int ha_assemble(void) ha_extract_print_list(&R_INF, asm_opt.extract_iter, asm_opt.extract_list); exit(0); } - if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN) && !(asm_opt.flag & HA_F_VERBOSE_GFA)) ha_triobin(&asm_opt); - ///if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt), ovlp_loaded = 2; + // if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN) && !(asm_opt.flag & HA_F_VERBOSE_GFA)) ha_triobin(&asm_opt); + if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt); + // if (!(asm_opt.flag & HA_F_SKIP_TRIOBIN)) ha_triobin(&asm_opt), ovlp_loaded = 2; if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads(); if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF(); if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2; diff --git a/CommandLines.h b/CommandLines.h index b0b583d..1fd6948 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.17.2-r432" +#define HA_VERSION "0.17.2-r433" #define VERBOSE 0 diff --git a/Correct.cpp b/Correct.cpp index b930f07..b86f597 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -15508,6 +15508,24 @@ void push_sec_aln(overlap_region *z, int64_t s, int64_t e, int64_t sec_err) p->x_start = s; p->x_end = e; p->clen += sec_err; } +void push_sec_aln_robust(overlap_region *z, int64_t s, int64_t e, int64_t sec_err) +{ + window_list *p; + if(z->align_length > 0) { + p = z->w_list.a + z->w_list.n + z->align_length - 1; + if((p->x_end == s) && (p->clen == 0) && ((!!(p->clen)) == (!!sec_err))) { + p->x_end = e; p->clen += sec_err; + return; + } + } + if((z->w_list.n+z->align_length)==z->w_list.m) { + z->w_list.m = z->w_list.m? z->w_list.m<<1 : 2; + z->w_list.a = (window_list*)realloc(z->w_list.a, sizeof(window_list)*z->w_list.m); + } + p = &(z->w_list.a[z->w_list.n+z->align_length]); z->align_length++; memset(p, 0, sizeof((*p))); + p->x_start = s; p->x_end = e; p->clen += sec_err; +} + // #define id_mm ((uint64_t)0x7fffffffffffffff) #define id_set ((uint64_t)0x8000000000000000) #define id_get(a) ((uint32_t)(a)) @@ -15564,6 +15582,8 @@ uint64_t gen_region_phase(overlap_region* ol, uint64_t *id_a, uint64_t id_n, uin // fprintf(stderr, "---[M::%s::utg%.6dl] wid::%u, xoff::%u, coff::%u, err::%ld\n", __func__, // (int32_t)ol[ovlp_id(*p)].y_id+1, ovlp_cur_wid(*p), ovlp_cur_xoff(*p), ovlp_cur_coff(*p), err); + fprintf(stderr, "---[M::%s::utg%.6dl] xoff::[%lu, %lu), err::%ld\n", __func__, + (int32_t)ol[ovlp_id(*p)].y_id+1, s, e, err); assert(err >= 0); if(err < msc) { msc = err; msc_k = k; msc_n = 1; @@ -15604,7 +15624,7 @@ uint64_t gen_region_phase(overlap_region* ol, uint64_t *id_a, uint64_t id_n, uin // (int32_t)ol[oid].y_id+1, s, e); } - for (k = 0; k < mn; k++) { + for (k = 0; k < mn; k++) {///best alignment z = &(ol[id_get(buf->a[k])]); reassign_sec_err(ol, ovidx, buf, k); push_sec_aln(z, s, e, 0); @@ -15634,6 +15654,111 @@ uint64_t gen_region_phase(overlap_region* ol, uint64_t *id_a, uint64_t id_n, uin } +///[s, e) +uint64_t gen_region_phase_robust(overlap_region* ol, uint64_t *id_a, uint64_t id_n, uint64_t s, uint64_t e, uint64_t dp, ul_ov_t *c_idx, asg64_v *buf) +{ + if(!id_n) return id_n; + uint64_t k, m, mn, q[2], buf_n, rm_n, oid; int64_t err, msc, msc_k, msc_n; + overlap_region *z; ul_ov_t *p; buf->n = 0; kv_resize(uint64_t, *buf, dp); + for (k = buf_n = rm_n = 0; k < id_n; k++) { + p = &(c_idx[id_a[k]]); + q[0] = ol[ovlp_id(*p)].w_list.a[ovlp_min_wid(*p)].x_start; + q[1] = ol[ovlp_id(*p)].w_list.a[ovlp_max_wid(*p)].x_end+1; + if(q[0]<=s && q[1]>=e) { + kv_push(uint64_t, *buf, id_a[k]); + // buf[buf_n++] = id_a[k]; + } + if(q[1] < e) rm_n++; + } + buf_n = buf->n; + assert(buf_n == dp);//not right + + if(buf_n > 0) { + for (k = 0, msc = INT32_MAX, msc_k = -1, msc_n = 0; k < buf_n; k++) { + p = &(c_idx[(uint32_t)buf->a[k]]); z = &(ol[ovlp_id(*p)]); + // fprintf(stderr, "+++[M::%s::utg%.6dl] wid::%u, xoff::%u, coff::%u\n", __func__, + // (int32_t)ol[ovlp_id(*p)].y_id+1, ovlp_cur_wid(*p), ovlp_cur_xoff(*p), ovlp_cur_coff(*p)); + err = extract_sub_cigar_err(z, s, e, p); + // int64_t debug_err = extract_sub_cigar_err_debug(z, s, e); + // assert(err == debug_err); + // fprintf(stderr, "[M::%s::] err::%ld, debug_err::%ld\n", __func__, err, debug_err); + + // fprintf(stderr, "---[M::%s::utg%.6dl] wid::%u, xoff::%u, coff::%u, err::%ld\n", __func__, + // (int32_t)ol[ovlp_id(*p)].y_id+1, ovlp_cur_wid(*p), ovlp_cur_xoff(*p), ovlp_cur_coff(*p), err); + // fprintf(stderr, "---[M::%s::utg%.6dl] xoff::[%lu, %lu), err::%ld\n", __func__, + // (int32_t)ol[ovlp_id(*p)].y_id+1, s, e, err); + assert(err >= 0); + if(err < msc) { + msc = err; msc_k = k; msc_n = 1; + } else if(err == msc) { + msc_n++; + } + buf->a[k] |= (((uint64_t)err)<<32); + } + + if(msc_n == 1) { + p = &(c_idx[(uint32_t)buf->a[msc_k]]); + z = &(ol[ovlp_id(*p)]); mn = 1; + if(msc_k != 0) { + m = buf->a[msc_k]; + buf->a[msc_k] = buf->a[0]; + buf->a[0] = m; + } + } else { + for (k = mn = 0; k < buf_n && (int64_t)mn < msc_n; k++) { + p = &(c_idx[(uint32_t)buf->a[k]]); + z = &(ol[ovlp_id(*p)]); + if((buf->a[k]>>32) == (uint64_t)msc) { + if(mn != k) { + m = buf->a[k]; + buf->a[k] = buf->a[mn]; + buf->a[mn] = m; + } + mn++; + } + } + } + // fprintf(stderr, "[M::%s] buf_n::%ld, msc_n::%ld, mn::%lu\n", __func__, buf_n, msc_n, mn); + for (k = 0; k < buf_n; k++) { + oid = ovlp_id((c_idx[(uint32_t)buf->a[k]])); + buf->a[k] >>= 32; buf->a[k] <<= 32; buf->a[k] |= oid; + // ol[oid].overlapLen = k;//no need to set + // if(s == 158482) fprintf(stderr, "k->%ld::oid->%ld[M::%s::utg%.6dl] pos::[%lu, %lu)\n", k, oid, __func__, + // (int32_t)ol[oid].y_id+1, s, e); + } + + for (k = 0; k < mn; k++) {///best alignment + z = &(ol[id_get(buf->a[k])]); + // reassign_sec_err(ol, ovidx, buf, k); + // push_sec_aln(z, s, e, 0); + push_sec_aln_robust(z, s, e, 0); + } + + for (k = mn; k < buf_n; k++) { + z = &(ol[id_get(buf->a[k])]); + // push_sec_aln(z, s, e, ((buf->a[k]&id_set)?(0):(err_get(buf->a[k])-msc))); + push_sec_aln_robust(z, s, e, (err_get(buf->a[k])-msc)); + } + + // for (k = 0; k < buf_n; k++) { + // ol[id_get(buf->a[k])].overlapLen = (uint32_t)-1; + // } + } + + + if(rm_n) { + for (k = m = 0; k < id_n; k++) { + p = &(c_idx[id_a[k]]); + q[1] = ol[ovlp_id(*p)].w_list.a[ovlp_max_wid(*p)].x_end+1; + if(q[1] < e) continue; + id_a[m++] = id_a[k]; + } + id_n = m; + } + return id_n; +} + + int64_t infer_rovlp(ul_ov_t *li, ul_ov_t *lj, uc_block_t *bi, uc_block_t *bj, All_reads *ridx, ma_ug_t *ug) { int64_t in, is, ie, irev, iqs, iqe, jn, js, je, jrev, jqs, jqe, ir, jr, ts, te, max_s, min_e, s_shift, e_shift; @@ -15830,7 +15955,7 @@ void gen_gov_idx(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t // } } -void prt_overlap_region_stat(overlap_region *z) +void prt_overlap_region_stat(overlap_region *z, int64_t sid) { uint64_t k = 0, aln = 0, ualn = 0, err = 0; for (k = 0; k < z->w_list.n; k++) { @@ -15841,21 +15966,20 @@ void prt_overlap_region_stat(overlap_region *z) err += z->w_list.a[k].error; } } - fprintf(stderr, "[M::%s::utg%.6dl::%c] q::[%d, %d), t::[%d, %d), aln::%lu, ualn::%lu, err::%lu, flen::%u, blen::%u, sec_err::%u\n", __func__, - (int32_t)z->y_id+1, "+-"[z->y_pos_strand], z->x_pos_s, z->x_pos_e+1, z->y_pos_s, z->y_pos_e+1, aln, ualn, err, - z->overlapLen, z->align_length, z->non_homopolymer_errors); + fprintf(stderr, "[M::%s::utg%.6dl::%c] sid::%ld, q::[%d, %d), t::[%d, %d), aln::%lu, ualn::%lu, err::%lu\n", + __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], sid, z->x_pos_s, z->x_pos_e+1, + z->y_pos_s, z->y_pos_e+1, aln, ualn, err); } -void prt_overlap_region_phase_stat(overlap_region *z) +void prt_overlap_region_phase_stat(overlap_region *z, int64_t sid) { uint64_t k = 0; - fprintf(stderr, "[M::%s::utg%.6dl::%c] q::[%d, %d), t::[%d, %d), best::%u, sec::%u\n", __func__, - (int32_t)z->y_id+1, "+-"[z->y_pos_strand], z->x_pos_s, z->x_pos_e+1, z->y_pos_s, z->y_pos_e+1, - z->align_length, z->non_homopolymer_errors); + fprintf(stderr, "[M::%s::utg%.6dl::%c] sid::%ld, q::[%d, %d), t::[%d, %d)\n", __func__, + (int32_t)z->y_id+1, "+-"[z->y_pos_strand], sid, z->x_pos_s, z->x_pos_e+1, z->y_pos_s, z->y_pos_e+1); for (k = 0; k < z->w_list.n; k++) { - fprintf(stderr, "[k::%lu] q::[%d, %d), sec::%u\n", k, + fprintf(stderr, "[k::%lu] sid::%ld, q::[%d, %d), sec::%u\n", k, sid, z->w_list.a[k].x_start, z->w_list.a[k].x_end, z->w_list.a[k].clen); } } @@ -15868,6 +15992,7 @@ void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t kv_resize(ul_ov_t, *c_idx, ol->length); for (k = idx->n = c_idx->n = 0; k < on; k++) { z = &(ol->list[k]); zwn = z->w_list.n; + // prt_overlap_region_stat(z, ulid); // if(ulid == 35437 && z->y_id == 109111) { // fprintf(stderr, "+0+[M::%s::utg%.6dl::%c] ulid::%ld, all::%u, non-best::%ld, best::%u\n", __func__, // (int32_t)z->y_id+1, "+-"[z->y_pos_strand], ulid, z->overlapLen, zwn, z->align_length); @@ -15927,7 +16052,8 @@ void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t // } } radix_sort_bc64(idx->a, idx->a+idx->n); - gen_gov_idx(ol, uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, buf1); + //this is used with gen_region_phase, now give up + //gen_gov_idx(ol, uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, buf1); // for (m = 0; m < c_idx->n; m++) { // fprintf(stderr, "+++[M::%s::utg%.6dl] q[%d, %d), t[%d, %d), wn::%d\n", __func__, // (int32_t)ol->list[ovlp_id(c_idx->a[m])].y_id+1, @@ -15954,7 +16080,8 @@ void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t if((end > beg) && (old_dp >= 2)) { // fprintf(stderr, "\n[M::%s::] beg::%ld, end::%ld, old_dp::%ld\n", __func__, beg, end, old_dp); // kv_resize(uint64_t, *buf, ((uint32_t)old_dp)<<1); - idx->n = srt_n + gen_region_phase(ol->list, idx->a+srt_n, idx->n-srt_n, beg, end, old_dp, c_idx->a, buf, buf1); + // idx->n = srt_n + gen_region_phase(ol->list, idx->a+srt_n, idx->n-srt_n, beg, end, old_dp, c_idx->a, buf, buf1); + idx->n = srt_n + gen_region_phase_robust(ol->list, idx->a+srt_n, idx->n-srt_n, beg, end, old_dp, c_idx->a, buf); } beg = end; } @@ -15990,6 +16117,7 @@ void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t // } assert(zwn <= z->overlapLen); z->align_length = z->overlapLen - zwn; + // prt_overlap_region_phase_stat(z, ulid); // fprintf(stderr, "[M::%s::utg%.6dl::%c] all::%u, non-best::%ld, best::%u\n", __func__, // (int32_t)z->y_id+1, "+-"[z->y_pos_strand], z->overlapLen, zwn, z->align_length); // prt_overlap_region_phase_stat(&(ol->list[k])); @@ -16328,6 +16456,13 @@ void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur ol->list[i] = t; } ol->list[k].is_match = 1; ol->list[k].non_homopolymer_errors = re; + // fprintf(stderr, "+[M::%s] on::%lu\n", __func__, ol->length); + /****for debug****/ + // z = &(ol->list[k]); + // fprintf(stderr, "[M::%s::utg%.6dl::%c] sid::%ld, q::[%d, %d), t::[%d, %d), err::%u\n", + // __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], sid, z->x_pos_s, z->x_pos_e+1, z->y_pos_s, z->y_pos_e+1, + // z->non_homopolymer_errors); + /****for debug****/ k++; } } diff --git a/Overlaps.cpp b/Overlaps.cpp index d048900..c1911cb 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -17798,6 +17798,7 @@ char *f_prefix, uint8_t *kpt_buf, kvec_asg_arc_t_warp *r_edges) return NULL; } + void filter_set_kug(uint8_t* trio_flag, asg_t *rg, uint8_t *rf, kvec_asg_arc_t_warp *r_edges, float f_rate, ma_ug_t **ug) { asg_t* nsg = (*ug)->g; ma_utg_t *u = NULL; @@ -24843,14 +24844,16 @@ int load_ruIndex(R_to_U* ruIndex, char* read_file_name) (ruIndex)->index = (uint32_t*)malloc(sizeof(uint32_t)*(ruIndex)->len); f_flag += fread((ruIndex)->index, sizeof((ruIndex)->index[0]), (ruIndex)->len, fp); - R_INF.trio_flag = (uint8_t*)malloc(sizeof(uint8_t)*(ruIndex)->len); - f_flag += fread(R_INF.trio_flag, sizeof(R_INF.trio_flag[0]), (ruIndex)->len, fp); + if(!(asm_opt.ar)) { + R_INF.trio_flag = (uint8_t*)malloc(sizeof(uint8_t)*(ruIndex)->len); + f_flag += fread(R_INF.trio_flag, sizeof(R_INF.trio_flag[0]), (ruIndex)->len, fp); + } // CALLOC(ruIndex->is_het, ruIndex->len); // f_flag += fread(ruIndex->is_het, 1, ruIndex->len, fp); - f_flag += fread(&(asm_opt.hom_global_coverage_set), sizeof(asm_opt.hom_global_coverage_set), 1, fp); - f_flag += fread(&(asm_opt.hom_global_coverage), sizeof(asm_opt.hom_global_coverage), 1, fp); + // f_flag += fread(&(asm_opt.hom_global_coverage_set), sizeof(asm_opt.hom_global_coverage_set), 1, fp); + // f_flag += fread(&(asm_opt.hom_global_coverage), sizeof(asm_opt.hom_global_coverage), 1, fp); free(index_name); fflush(fp); @@ -25187,7 +25190,7 @@ uint32_t is_collect_trans) MALLOC(x->pos_idx, x->n); memset(x->pos_idx, -1, x->n*sizeof(uint64_t)); CALLOC(set, read_g->n_seq<<1); CALLOC(x->cov, x->n); - for (i = 0; i < n_ux; i++) + for (i = 0; i < n_ux; i++)//get the coverage for each read { if(ug->g->seq[i].del) continue; u = &(ug->u.a[i]); @@ -31767,10 +31770,11 @@ ma_sub_t **coverage_cut_ptr, int debug_g) debug_gfa:; gen_ug_opt_t(&uopt, sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, &b_mask_t); + set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, max_hang_length, mini_overlap_length); } if(asm_opt.ar) { ul_realignment_gfa(&uopt, sg, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, - asm_opt.max_short_tip, &b_mask_t, ha_opt_triobin(&asm_opt)); + asm_opt.max_short_tip, &b_mask_t, ha_opt_triobin(&asm_opt), o_file); } // print_debug_gfa(sg, NULL, coverage_cut, "UL.debug", sources, ruIndex, max_hang_length, mini_overlap_length, 0, 0, 0); /** diff --git a/Overlaps.h b/Overlaps.h index 3f7996d..3a05afa 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -1141,6 +1141,10 @@ int64_t count_edges_v_w(asg_t *g, uint32_t v, uint32_t w); void renew_utg(ma_ug_t **ug, asg_t* read_g, kvec_asg_arc_t_warp* edge); void merge_unitig_content(ma_utg_t* collection, ma_ug_t* ug, asg_t* read_g, kvec_asg_arc_t_warp* edge); +// void break_ug_contig(ma_ug_t **ug, asg_t *read_g, All_reads *RNF, ma_sub_t *coverage_cut, +// ma_hit_t_alloc* sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, +// int* b_low_cov, int* b_high_cov, double m_rate); + #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 8e4ecd9..6b6b4fe 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -6217,8 +6217,8 @@ void gen_integer_normalize(ul_resolve_t *uidx) void clip_integer_chimeric(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *o, ul2ul_idx_t *ul2, integer_t *buf, int64_t min_dp) { - uint64_t k, is_srt = 0, is_del = 0, rid, zs, ze; assert(o->id == qid); - ul_str_t *str = &(uidx->pstr.str.a[qid]); int64_t z, z_n; + uint64_t k, is_srt = 0, is_del = 0/**, rid, zs, ze**/; assert(o->id == qid); + ul_str_t *str = &(uidx->pstr.str.a[qid]); ///int64_t z, z_n; for (k = 0; k < o->cn; k++) { if(o->a[k].is_del) break; if(k > 0 && o->a[k].hid < o->a[k-1].hid) break; @@ -6249,8 +6249,10 @@ void clip_integer_chimeric(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *o, ul } for (k = 0, buf->u.n = 0; k < o->cn; k++) { - kv_push(uint64_t, buf->u, (o->a[k].qs<<1)); - kv_push(uint64_t, buf->u, (o->a[k].qe<<1)|1); + // kv_push(uint64_t, buf->u, (o->a[k].qs<<1)); + kv_push(uint64_t, buf->u, (o->a[k].qs_k<<1)); + // kv_push(uint64_t, buf->u, (o->a[k].qe<<1)|1); + kv_push(uint64_t, buf->u, (o->a[k].qe_k<<1)|1); } radix_sort_srt64(buf->u.a, buf->u.a + buf->u.n); @@ -6278,20 +6280,20 @@ void clip_integer_chimeric(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *o, ul start = buf->u.a[k]>>32; end = (uint32_t)buf->u.a[k]; if(start == 0) is_left = 1; else is_middle = 1; - if(end == uidx->idx->a[qid].rlen) is_right = 1; + if(end == str->cn/**uidx->idx->a[qid].rlen**/) is_right = 1; else is_middle = 1; } if(is_left == 0 && is_right == 0) is_del = 1; if(is_left && is_right && is_middle) is_del = 1; } - + /** if(!is_del) { for (k = 0; k < str->cn; k++) { rid = (((uint32_t)str->a[k])>>1); if(!IF_HOM(rid, *(uidx->bub))) break; } - if(k < str->cn) {///at least a hom node covered + if(k < str->cn) {///at least a het node covered for (k = 0, buf->u.n = 0; k < o->cn; k++) { z = o->a[k].qs_k; z_n = o->a[k].qe_k; for (; z < z_n; z++) { @@ -6338,7 +6340,7 @@ void clip_integer_chimeric(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *o, ul start = buf->u.a[k]>>32; end = (uint32_t)buf->u.a[k]; if(start == 0) is_left = 1; else is_middle = 1; - if(end == uidx->idx->a[qid].rlen) is_right = 1; + if(end == str->cn) is_right = 1; else is_middle = 1; } if(is_left == 0 && is_right == 0) is_del = 1; @@ -6346,6 +6348,7 @@ void clip_integer_chimeric(ul_resolve_t *uidx, uint32_t qid, ul2ul_item_t *o, ul } } } + **/ if(is_del) { o->is_del = 1; @@ -10862,11 +10865,22 @@ uint32_t check_hybrid_connect(usg_t *ng, uint32_t i_uid, uint32_t v, uint32_t vi return 0; } +uint32_t usg_arc_occ(usg_t *ng, uint32_t v) +{ + usg_arc_t *p; uint32_t nv, kv, i; + p = usg_arc_a(ng, v); nv = usg_arc_n(ng, v); + for (i = kv = 0; i < nv; i++) { + if(p[i].del) continue; + kv++; + } + return kv; +} + void integer_realign_g(ul_resolve_t *uidx, usg_t *ng, uinfo_srt_warp_t *seq, uint32_t seq_id, integer_t *buf) { if(seq->n < 2) return; // if(seq_id != 1137) return; - uint64_t *srt, *track, k, z, m, t, i, l, seq_n, j, sc, csc, mm_sc, mm_idx, n_v, n_u, n_v0, *p; + uint64_t *srt, *track, k, z, m, t, i, l, seq_n, j, sc, csc, mm_sc, mm_idx, n_v, n_u, n_v0, *p, lin, pn; uint32_t vi, vj; mmap_t *zm, *zt; ///ng->map: the nodes in the new graph that are mapped to the initial HiFi graph for (i = seq_n = 0; i < seq->n; ++i) seq_n += ng->mp.a[seq->a[i].v>>1].n; @@ -10966,15 +10980,23 @@ void integer_realign_g(ul_resolve_t *uidx, usg_t *ng, uinfo_srt_warp_t *seq, uin buf->u.a[n_u++] = buf->u.a[k]; } buf->u.n = n_u; buf->res_dump.n = l; - for (k = 0; k < buf->u.n; k++) { + for (k = n_u = 0; k < buf->u.n; k++) { // t = (uint32_t)-1; t <<= 32; t = seq_id; t <<= 32; t |= ((uint64_t)0x8000000000000000); - t |= (((uint32_t)buf->u.a[k])-(buf->u.a[k]>>32)); - kv_push(uint64_t, buf->res_dump, t); + t |= (((uint32_t)buf->u.a[k])-(buf->u.a[k]>>32)); + pn = buf->res_dump.n; kv_push(uint64_t, buf->res_dump, t); n_v0 = buf->u.a[k]>>32; n_v = (uint32_t)buf->u.a[k]; - for (z = n_v0; z < n_v; z++) { + for (z = n_v0, lin = 1; z < n_v; z++) { + if(lin && (z > n_v0)) { + if((usg_arc_occ(ng, ((uint32_t)track[z])^1) > 1) || + (usg_arc_occ(ng, ((uint32_t)track[z-1])) > 1)) { + lin = 0; + } + } kv_push(uint64_t, buf->res_dump, (uint32_t)track[z]); } + + if(lin == 1) buf->res_dump.n = pn; } } @@ -11321,6 +11343,10 @@ void update_usg_t_threading_0(usg_t *ng, uint64_t *a, uint64_t a_n, uint32_t *oc } occ[a[k]]--; occ[a[k]^1]--; } + // if(occ[a[a_n-1]^1] != 1) { + // fprintf(stderr, "[M::%s] utg%.6dl(%c)(occ::%u)\n", __func__, + // (int32_t)((a[a_n-1]^1)>>1)+1, "+-"[(a[a_n-1]^1)&1], occ[a[a_n-1]^1]); + // } assert(occ[a[a_n-1]^1] == 1); kv_pushp(uint64_t, *b, &p); (*p) = a[a_n-1]^1; (*p) <<= 32; (*p) |= a[a_n-1]^1; occ[a[a_n-1]^1]--; @@ -11398,6 +11424,11 @@ void update_usg_t_threading(ul_resolve_t *uidx, usg_t *ng, uint64_t *arcs, uint6 occ[integ_seq[k]]++; occ[integ_seq[k]^1]++; } occ[integ_seq[s]]++; occ[integ_seq[e]^1]++; + // if(occ[integ_seq[s]] > 1 || occ[integ_seq[e]^1] > 1) { + // fprintf(stderr, "s::utg%.6dl(%c)(occ::%u), e::utg%.6dl(%c)(occ::%u)\n", + // (int32_t)(integ_seq[s]>>1)+1, "+-"[integ_seq[s]&1], occ[integ_seq[s]], + // (int32_t)((integ_seq[e]^1)>>1)+1, "+-"[((integ_seq[e]^1)&1)], occ[integ_seq[e]^1]); + // } } for (i = 0; i < nvtx; i++) { @@ -11722,13 +11753,438 @@ void prt_intg_info(uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, const void prt_usg_t(ul_resolve_t *uidx, usg_t *ng, const char *cmd); -uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, uint32_t max_ext) + +#define occ_m(x) ((x)&((uint32_t)0x7fffffff)) +#define c_unqiue_m(v, occ) (occ_m((occ)[(v)]) == 1 && occ_m((occ)[(v)^1]) <= 1) +#define reli_pass(x) (!((x)&((uint64_t)0x8000000000000000))) + +uint64_t get_ext_tip(usg_t *g, uint64_t *seq, uint64_t s, uint64_t e, asg64_v *buf, asg64_v *set, uint8_t *f, uint64_t max_ext) +{ + uint64_t bn = buf->n, sn = set->n, v, k, nv, z, n_ext = 0, i, kv; usg_arc_t *av; + for (k = s + 1; k < e; k++) {///note: here is [s, e] + v = seq[k]; av = usg_arc_a(g, v); nv = usg_arc_n(g, v); + for (z = 0; z < nv; z++) { + if(av[z].del || f[av[z].v^1]) continue; + kv_push(uint64_t, *buf, av[z].v); + } + + v = seq[k]^1; av = usg_arc_a(g, v); nv = usg_arc_n(g, v); + for (z = 0; z < nv; z++) { + if(av[z].del || f[av[z].v^1]) continue; + kv_push(uint64_t, *buf, av[z].v); + } + } + + v = seq[s]; av = usg_arc_a(g, v); nv = usg_arc_n(g, v); + for (z = 0; z < nv; z++) { + if(av[z].del || f[av[z].v^1]) continue; + kv_push(uint64_t, *buf, av[z].v); + } + + v = seq[e]^1; av = usg_arc_a(g, v); nv = usg_arc_n(g, v); + for (z = 0; z < nv; z++) { + if(av[z].del || f[av[z].v^1]) continue; + kv_push(uint64_t, *buf, av[z].v); + } + + + while (buf->n > bn && n_ext < max_ext) { + v = buf->a[--buf->n]; if(f[v]) continue; + av = usg_arc_a(g, v^1); nv = usg_arc_n(g, v^1); + for (i = kv = 0; i < nv && kv < 1; i++) { + if (av[i].del || f[av[i].v^1]) continue; + kv++; + } + if(kv > 0) continue; + n_ext += g->a[v>>1].occ; + if(!(f[v])) { + f[v] = 1; kv_push(uint64_t, *set, v); + } + if(!(f[v^1])) { + f[v^1] = 1; kv_push(uint64_t, *set, (v^1)); + } + + av = usg_arc_a(g, v); nv = usg_arc_n(g, v); + for (i = 0; i < nv; i++) { + if (av[i].del || f[av[i].v^1]) continue; + kv_push(uint64_t, *buf, av[i].v); + } + } + buf->n = bn; + for (k = sn; k < set->n; k++) f[set->a[k]] = 0; + return n_ext; +} + +void select_unqiue_path_on_fly(usg_t *g, uint64_t s, uint64_t e, asg64_v *b64, asg64_v *ub64, +uint64_t *integer_seq, uint32_t *ng_occ, uint64_t *is_reli, uint8_t *f, uint64_t max_ext, asg64_v *res) +{ + uint64_t k, sk = (uint64_t)-1, ek = (uint64_t)-1, *p = NULL; ub64->n = 0; + assert(occ_m(ng_occ[integer_seq[s]]) == 1); + assert(occ_m(ng_occ[integer_seq[e]^1]) == 1); + + if(reli_pass(is_reli[integer_seq[s]])) { + sk = s; + } else { + for (k = s + 1; k < e; k++) { + if((reli_pass(is_reli[integer_seq[k]])) && (occ_m(ng_occ[integer_seq[k]]) == 1)) { + sk = k; + break; + } + } + } + if((sk == ((uint64_t)-1)) || (sk >= e)) return; + + for (k = sk + 1; k <= e; k++) { + if(sk != ((uint64_t)-1)) { + if(reli_pass(is_reli[(integer_seq[k]^1)])) { + if(occ_m((ng_occ[integer_seq[k]]^1)) == 1) { + ek = k; ///check and extend [sk, ek] + if((ek > sk) && (get_ext_tip(g, integer_seq, sk, ek, b64, ub64, f, max_ext) < max_ext)) { + p = ((res->n > 0)? (&(res->a[res->n-1])):NULL); + if(p && (((uint32_t)(*p)) == sk)) { + (*p) >>= 32; (*p) <<= 32; (*p) |= ek; + } else { + kv_pushp(uint64_t, *res, &p); *p = sk; (*p) <<= 32; (*p) |= ek; + } + } + sk = ek = ((uint64_t)-1); + } + } else { + sk = ek = ((uint64_t)-1); + } + } + + ///case2: sk != (uint64_t)-1 && ek == (uint64_t)-1 -> wait for an available ek + ///case3: sk == (uint64_t)-1 && ek == (uint64_t)-1 -> not available + if(k < e) { + if((!(reli_pass(is_reli[integer_seq[k]])))) { + sk = ek = (uint64_t)-1; continue; + } + if((sk == (uint64_t)-1) && (occ_m(ng_occ[integer_seq[k]]) == 1)) { + sk = k; + } + } + } +} + +void select_unqiue_path_on_fly_raw(usg_t *g, uint64_t s, uint64_t e, asg64_v *b64, asg64_v *ub64, +uint64_t *integer_seq, uint32_t *ng_occ, uint64_t *is_reli, uint8_t *f, uint64_t max_ext, asg64_v *res) +{ + uint64_t k, sk = (uint64_t)-1, ek = (uint64_t)-1, *p = NULL; ub64->n = 0; + assert(occ_m(ng_occ[integer_seq[s]]) == 1); + assert(occ_m(ng_occ[integer_seq[e]^1]) == 1); + sk = s; ek = e; + if(get_ext_tip(g, integer_seq, sk, ek, b64, ub64, f, max_ext) < max_ext) { + kv_pushp(uint64_t, *res, &p); *p = sk; (*p) <<= 32; (*p) |= ek; + // fprintf(stderr, "+ext_k::[%lu, %lu]\ts::%lu\te::%lu\n", sk, ek, s, e); + // fprintf(stderr, "occ[sk]::%u\tocc[ek]::%u\tocc[s]::%u\tocc[e]::%u\n", + // occ_m(ng_occ[integer_seq[sk]]), occ_m(ng_occ[integer_seq[ek]^1]), + // occ_m(ng_occ[integer_seq[s]]), occ_m(ng_occ[integer_seq[e]^1])); + } + return; + + + + for (k = sk + 1; k <= e; k++) { + if(sk != ((uint64_t)-1)) { + if(reli_pass(is_reli[(integer_seq[k]^1)])) { + if(occ_m((ng_occ[integer_seq[k]]^1)) == 1) { + ek = k; ///check and extend [sk, ek] + if((ek > sk) && (get_ext_tip(g, integer_seq, sk, ek, b64, ub64, f, max_ext) < max_ext)) { + p = ((res->n > 0)? (&(res->a[res->n-1])):NULL); + if(p && (((uint32_t)(*p)) == sk)) { + (*p) >>= 32; (*p) <<= 32; (*p) |= ek; + } else { + kv_pushp(uint64_t, *res, &p); *p = sk; (*p) <<= 32; (*p) |= ek; + } + } + sk = ek = ((uint64_t)-1); + } + } else { + sk = ek = ((uint64_t)-1); + } + } + + ///case2: sk != (uint64_t)-1 && ek == (uint64_t)-1 -> wait for an available ek + ///case3: sk == (uint64_t)-1 && ek == (uint64_t)-1 -> not available + if(k < e) { + if((!(reli_pass(is_reli[integer_seq[k]])))) { + sk = ek = (uint64_t)-1; continue; + } + if((sk == (uint64_t)-1) && (occ_m(ng_occ[integer_seq[k]]) == 1)) { + sk = k; + } + } + } +} + +void update_thread_path(usg_t *g, asg64_v *b64, asg64_v *ub64, uint64_t g_s, uint64_t g_e, uint64_t *integer_seq, +uint32_t *ng_occ, uint64_t *reliable_idx, uint8_t *f, uint64_t max_ext, asg64_v *res) +{ + uint64_t i, k, s, e; + for (i = g_s; i < g_e; i++) {///available intervals within the same cluster + s = b64->a[((uint32_t)b64->a[i])]>>32; e = ((uint32_t)b64->a[((uint32_t)b64->a[i])]); assert(s < e); + for (k = s + 1; k < e; k++) {///note: here is [s, e] + f[integer_seq[k]] = f[integer_seq[k]^1] = 1; + } + f[integer_seq[s]] = f[integer_seq[e]^1] = 1; + } + + for (i = g_s; i < g_e; i++) {///available intervals within the same cluster + s = b64->a[((uint32_t)b64->a[i])]>>32; e = ((uint32_t)b64->a[((uint32_t)b64->a[i])]); assert(s < e); + // fprintf(stderr, "[M::%s::] i::%lu, chain_id::%u\n", __func__, i, ((uint32_t)b64->a[i])); + // select_unqiue_path_on_fly(g, s, e, b64, ub64, integer_seq, ng_occ, reliable_idx, f, max_ext, res); + select_unqiue_path_on_fly_raw(g, s, e, b64, ub64, integer_seq, ng_occ, reliable_idx, f, max_ext, res); + } + + for (i = g_s; i < g_e; i++) {///available intervals within the same cluster + s = b64->a[((uint32_t)b64->a[i])]>>32; e = ((uint32_t)b64->a[((uint32_t)b64->a[i])]); assert(s < e); + for (k = s + 1; k < e; k++) {///note: here is [s, e] + f[integer_seq[k]] = f[integer_seq[k]^1] = 0; + } + f[integer_seq[s]] = f[integer_seq[e]^1] = 0; + } +} + + +typedef struct { + ul_resolve_t *uidx; + uint64_t *integer_seq; + uint64_t *gidx; + uint64_t *interval_idx; + uint64_t gidx_n; + uint8_t *f; + usg_t *ng; + uint64_t *int_idx; + uint64_t int_idx_n; + uint64_t *int_a; + uint64_t int_an; + uint64_t *ridx_a; + uint64_t *ridx; +} unique_bridge_check_t; + +uint64_t get_arc_support(ul_resolve_t *uidx, uint64_t v, uint64_t w) +{ + uint64_t *hid_a, hid_n; ul_str_idx_t *str_idx = &(uidx->pstr); + uint64_t z, vz, occ = 0; ul_str_t *str; int64_t s, s_n; + hid_a = str_idx->occ.a + str_idx->idx.a[v>>1]; + hid_n = str_idx->idx.a[(v>>1)+1] - str_idx->idx.a[v>>1]; + + for (z = occ = 0; z < hid_n; z++) { + str = &(str_idx->str.a[hid_a[z]>>32]); s_n = str->cn; + if(s_n < 2) continue; + vz = (uint32_t)(str->a[(uint32_t)hid_a[z]]); + assert((v>>1) == (vz>>1)); + + if(v == vz) { + for(s = ((uint32_t)hid_a[z]) + 1; s < s_n; s++) { + if(((uint32_t)(str->a[s])) == w) { + occ++; break; + } + } + } else { + for(s = ((int32_t)((uint32_t)hid_a[z]))-1; s >= 0; s--) { + if(((uint32_t)(str->a[s])) == (w^1)) { + occ++; break; + } + } + } + } + return occ; +} +#define unique_bridge_occ 2 +#define unique_bridge_rate 0.499999 + +// static void worker_unique_bridge_check(void *data, long i, int tid) // callback for kt_for() +// { +// unique_bridge_check_t *uaux = (unique_bridge_check_t *)data; +// uint64_t *integer_seq = uaux->integer_seq, s, e, v, w, nse, self_k = i, k, kv; +// uint64_t *gidx = uaux->gidx, *interval = uaux->interval_idx; +// uint8_t *f = uaux->f; + +// s = interval[((uint32_t)gidx[i])]>>32; v = integer_seq[s]; +// e = ((uint32_t)interval[((uint32_t)gidx[i])]); w = integer_seq[e]^1; +// assert(s < e); + +// nse = get_arc_support(uaux->uidx, v, w^1); +// if(nse < unique_bridge_occ) { +// f[v] = f[w] = 1; return; +// } + +// for (k = 0; k < uaux->gidx_n; k++) { +// if(k == self_k) continue; +// s = interval[((uint32_t)gidx[k])]>>32; +// e = ((uint32_t)interval[((uint32_t)gidx[k])]); +// kv = get_arc_support(uaux->uidx, v, integer_seq[s]^1); +// if((kv > 0) && (kv >= (nse*unique_bridge_rate))) { +// f[v] = f[w] = 1; return; +// } + +// kv = get_arc_support(uaux->uidx, v, integer_seq[e]); +// if((kv > 0) && (kv >= (nse*unique_bridge_rate))) { +// f[v] = f[w] = 1; return; +// } + +// kv = get_arc_support(uaux->uidx, w, integer_seq[s]^1); +// if((kv > 0) && (kv >= (nse*unique_bridge_rate))) { +// f[v] = f[w] = 1; return; +// } + +// kv = get_arc_support(uaux->uidx, w, integer_seq[e]); +// if((kv > 0) && (kv >= (nse*unique_bridge_rate))) { +// f[v] = f[w] = 1; return; +// } +// } +// } + +static void worker_unique_bridge_check_s(void *data, long i, int tid) // callback for kt_for() +{ + unique_bridge_check_t *uaux = (unique_bridge_check_t *)data; + uint64_t *integer_seq = uaux->integer_seq, s, e, v, w, nse, k; + uint64_t *gidx = uaux->gidx, *interval = uaux->interval_idx; + bubble_type *bub = uaux->uidx->bub; + + s = interval[((uint32_t)gidx[i])]>>32; + e = ((uint32_t)interval[((uint32_t)gidx[i])]); + assert(s < e); + for (k = s, v = w = (uint32_t)-1; k <= e; k++) { + if(k > s) { + nse = get_arc_support(uaux->uidx, integer_seq[k-1], integer_seq[k]); + if(nse < unique_bridge_occ) { + // fprintf(stderr, "+utg%.6dl(%c)\tutg%.6dl(%c)\tnse::%lu\n", + // ((int32_t)(integer_seq[k-1]>>1))+1, "+-"[integer_seq[k-1]&1], + // ((int32_t)(integer_seq[k]>>1))+1, "+-"[integer_seq[k]&1], nse); + gidx[i] |= ((uint64_t)0x8000000000000000); return; + } + } + if(IF_HOM((integer_seq[k]>>1), *bub)) continue; + w = integer_seq[k]; + if(v != (uint32_t)-1) { + nse = get_arc_support(uaux->uidx, v, w); + if(nse < unique_bridge_occ) { + // fprintf(stderr, "-utg%.6dl(%c)\tutg%.6dl(%c)\tnse::%lu\n", + // ((int32_t)(v>>1))+1, "+-"[v&1], + // ((int32_t)(w>>1))+1, "+-"[w&1], nse); + gidx[i] |= ((uint64_t)0x8000000000000000); return; + } + } + v = w; + } +} + +uint32_t ava_pass_unique_bridge_cov(ul_resolve_t *uidx, usg_t *g, asg64_v *b64, uint64_t g_s, uint64_t g_e, uint64_t *integer_seq, uint64_t tip_l) +{ + if((tip_l == 0) && (g_e - g_s == 1)) return 1; ///if only one path, go through in anyway + uint64_t i, is_del = 0; unique_bridge_check_t uaux; + + uaux.uidx = uidx; uaux.integer_seq = integer_seq; + uaux.gidx = b64->a + g_s; uaux.gidx_n = g_e - g_s; + uaux.interval_idx = b64->a; + kt_for(uidx->str_b.n_thread, worker_unique_bridge_check_s, (&uaux), g_e-g_s);///seq->n > 1 + + for (i = g_s; i < g_e; i++) {///available intervals within the same cluster + if(b64->a[i]&((uint64_t)0x8000000000000000)) { + b64->a[i] -= ((uint64_t)0x8000000000000000); is_del = 1; + } + } + return (!is_del); +} + + + +uint64_t old_path_ext(ul_resolve_t *uidx, usg_t *ng, asg64_v *b64, asg64_v *ub64, uint64_t a_n, uint64_t n_clus, uint64_t *i_idx, uint64_t *int_a, +uint8_t *ff, uint32_t *ng_occ, uint32_t max_ext) +{ + uint64_t k, i, z, /**s, e,**/ mm, tip_l; + /** + for (k = a_n + 1, i = a_n, mm = a_n; k <= b64->n; k++) { + if(k == b64->n || (b64->a[k]>>32) != (b64->a[i]>>32)) { + for (z = i; z < k; z++) {///all intger seqs within the same cluster + s = b64->a[((uint32_t)b64->a[z])]>>32; + e = ((uint32_t)b64->a[((uint32_t)b64->a[z])]); assert(e > s); + ///[s, e]:: available interval + if(!ava_pass_unique_bridge(i_idx, int_a, s, e)) break; + } + if(z >= k) {///all arcs in this cluster is fine -> each of arch is reliable + for (z = i; z < k; z++) b64->a[mm++] = b64->a[z]; + } + i = k; n_clus--; + } + } + assert(n_clus == 0); + b64->n = mm; + **/ + // prt_thread_info(b64.a, a_n, NULL, 0, "tt5"); + // fprintf(stderr, "**4**[M::%s::] a_n::%lu, n_clus::%lu\n", __func__, a_n, n_clus); + n_clus = 0; + for (k = a_n + 1, i = a_n, mm = a_n; k <= b64->n; k++) { + if(k == b64->n || (b64->a[k]>>32) != (b64->a[i]>>32)) { + tip_l = ava_pass_unique_bridge_tips(ng, b64, i, k, int_a, ff, max_ext); + if(tip_l < max_ext) {///no long tip + if(ava_pass_unique_bridge_cov(uidx, ng, b64, i, k, int_a, tip_l)) + { + for (z = i; z < k; z++) { + b64->a[mm++] = b64->a[z]; + } + n_clus++; + } + } + i = k; + } + } + // prt_thread_info(b64.a, a_n, NULL, 0, "tt6"); + b64->n = mm; + // fprintf(stderr, "**5**[M::%s::] a_n::%lu, n_clus::%lu, ng->n::%lu\n", __func__, a_n, n_clus, ng->n); + // prt_thread_info(b64.a, a_n, b64.a + a_n, b64.n - a_n, "thred"); + if(n_clus > 0) { + update_usg_t_threading(uidx, ng, b64->a, b64->a + a_n, b64->n - a_n, int_a, ng_occ, ub64); + } + // fprintf(stderr, "**6**[M::%s::] a_n::%lu, n_clus::%lu, ng->n::%lu\n", __func__, a_n, n_clus, ng->n); + return n_clus; +} + + + +void new_path_ext(ul_resolve_t *uidx, usg_t *ng, asg64_v *b64, asg64_v *ub64, uint64_t a_n, uint64_t n_clus, uint64_t *i_idx, uint64_t *int_a, +uint8_t *ff, uint32_t *ng_occ, uint32_t max_ext) +{ + asg64_v res; kv_init(res); uint64_t k, i, s, e, nvtx = ng->n<<1; + for (k = a_n + 1, i = a_n, res.n = 0; k <= b64->n; k++) { + if(k == b64->n || (b64->a[k]>>32) != (b64->a[i]>>32)) { + update_thread_path(ng, b64, ub64, i, k, int_a, ng_occ, i_idx, ff, max_ext, &res); + i = k; n_clus--; + } + } + assert(n_clus == 0); + + + if(res.n > 0) { + memset(ng_occ, 0, sizeof((*ng_occ))*nvtx); + for (i = 0; i < res.n; i++) { + s = res.a[i]>>32; e = ((uint32_t)res.a[i]); assert(s < e); + for (k = s + 1; k < e; k++) {///note: here is [s, e] + ng_occ[int_a[k]]++; ng_occ[int_a[k]^1]++; + } + ng_occ[int_a[s]]++; ng_occ[int_a[e]^1]++; + } + for (i = 0; i < nvtx; i++) { + if(!ng_occ[i]) ng_occ[i] = (uint32_t)-1; + } + + for (i = 0, ub64->n = 0; i < res.n; i++) { + s = res.a[i]>>32; e = ((uint32_t)res.a[i]); + update_usg_t_threading_0(ng, int_a + s, e + 1 - s, ng_occ, ub64); + } + } + kv_destroy(res); +} + +uint64_t gen_unique_g_adv_old(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, uint32_t max_ext) { - #define occ_m(x) ((x)&((uint32_t)0x7fffffff)) - #define c_unqiue_m(v, occ) (occ_m((occ)[(v)]) == 1 && occ_m((occ)[(v)^1]) <= 1) uint64_t pi, ei, *pz, bn, a_n, ua_n, *r_a, r_n; uint8_t *ff; CALLOC(ff, ng->n<<1); - uint32_t i, k, n_vtx = ng->n<<1, s, e, z, zs, ze, mm, v, b64_n; - uint32_t *ng_occ; CALLOC(ng_occ, n_vtx); asg64_v b64, ub64; kv_init(b64); kv_init(ub64); + uint32_t i, k, n_vtx = ng->n<<1, s, e, zs, ze, v, b64_n; + uint32_t *ng_occ; CALLOC(ng_occ, n_vtx); + asg64_v b64, ub64; kv_init(b64); kv_init(ub64); for (k = 0; k < int_idx_n; k++) {///scan all integer contigs r_a = int_a + (int_idx[k]>>32); r_n = (uint32_t)int_idx[k]; @@ -11841,6 +12297,12 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint b64.n = a_n; assert(ua_n == 0); for (i = 0; i < a_n; i++) { ///available intervals s = b64.a[i]>>32; e = (uint32_t)b64.a[i]; assert(e > s);///[s, e] -> coordinates within int_a[] + if(i == 3201 || i == 3202) { + fprintf(stderr, "[M::%s::] s::%u, e::%u, int_a[s]::%lu(occ::%u), int_a[e]^1::%lu(occ::%u)\n", + __func__, s, e, int_a[s], occ_m(ng_occ[int_a[s]]), + int_a[e]^1, occ_m(ng_occ[int_a[e]^1])); + } + assert(!(b64.a[i]&((uint64_t)0x8000000000000000))); for (k = s + 1; k < e; k++) {///note: here is [s, e]; s && e are unique, but [s+1, e-1] are not unique kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[k]]++; @@ -11871,49 +12333,155 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint ///b64.a[a_n, b64.n]:: (raw unitig/non-unqiue node id)|(resolvable path id) // n_clus = usg_unique_arcs_cluster(&b64, a_n, i_idx, int_a);///this function might be wrong n_clus = usg_unique_arcs_cluster_adv(&b64, a_n, i_idx, int_a, &ub64); - - // fprintf(stderr, "**3**[M::%s::] a_n::%lu, n_clus::%lu\n", __func__, a_n, n_clus); + fprintf(stderr, "**3**[M::%s::] a_n::%lu, n_clus::%lu\n", __func__, a_n, n_clus); // prt_thread_info(b64.a, a_n, NULL, 0, "tt4"); assert(b64.n == (a_n<<1)); - for (k = a_n + 1, i = a_n, mm = a_n; k <= b64.n; k++) { - if(k == b64.n || (b64.a[k]>>32) != (b64.a[i]>>32)) { - for (z = i; z < k; z++) {///all intger seqs within the same cluster - s = b64.a[((uint32_t)b64.a[z])]>>32; e = ((uint32_t)b64.a[((uint32_t)b64.a[z])]); assert(e > s); - ///[s, e]:: available interval - if(!ava_pass_unique_bridge(i_idx, int_a, s, e)) break; - } - if(z >= k) {///all arcs in this cluster is fine -> each of arch is reliable - for (z = i; z < k; z++) b64.a[mm++] = b64.a[z]; - } - i = k; n_clus--; - } - } - assert(n_clus == 0); - // prt_thread_info(b64.a, a_n, NULL, 0, "tt5"); - // fprintf(stderr, "**4**[M::%s::] a_n::%lu, n_clus::%lu\n", __func__, a_n, n_clus); - b64.n = mm; n_clus = 0; - for (k = a_n + 1, i = a_n, mm = a_n; k <= b64.n; k++) { - if(k == b64.n || (b64.a[k]>>32) != (b64.a[i]>>32)) { - if(ava_pass_unique_bridge_tips(ng, &b64, i, k, int_a, ff, max_ext) < max_ext) {///no long tip - for (z = i; z < k; z++) b64.a[mm++] = b64.a[z]; - n_clus++; - } - i = k; - } - } - // prt_thread_info(b64.a, a_n, NULL, 0, "tt6"); - b64.n = mm; - // fprintf(stderr, "**5**[M::%s::] a_n::%lu, n_clus::%lu, ng->n::%lu\n", __func__, a_n, n_clus, ng->n); - // prt_thread_info(b64.a, a_n, b64.a + a_n, b64.n - a_n, "thred"); - if(n_clus > 0) { - update_usg_t_threading(uidx, ng, b64.a, b64.a + a_n, b64.n - a_n, int_a, ng_occ, &ub64); - } - // fprintf(stderr, "**6**[M::%s::] a_n::%lu, n_clus::%lu, ng->n::%lu\n", __func__, a_n, n_clus, ng->n); + old_path_ext(uidx, ng, &b64, &ub64, a_n, n_clus, i_idx, int_a, ff, ng_occ, max_ext); + // new_path_ext(uidx, ng, &b64, &ub64, a_n, n_clus, i_idx, int_a, ff, ng_occ, max_ext); - free(ng_occ); free(i_idx); free(ff); kv_destroy(b64); kv_destroy(ub64); + free(ng_occ); free(i_idx); free(ff); kv_destroy(b64); kv_destroy(ub64); return b64.n - a_n; } +uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, uint32_t max_ext) +{ + uint64_t pi, ei, *pz, bn, a_n, ua_n, *r_a, r_n; uint8_t *ff; CALLOC(ff, ng->n<<1); + uint32_t i, k, n_vtx = ng->n<<1, s, e, zs, ze, v, b64_n; + uint32_t *ng_occ; CALLOC(ng_occ, n_vtx); + asg64_v b64, ub64; kv_init(b64); kv_init(ub64); + uint64_t *i_idx, n_clus; CALLOC(i_idx, ng->n<<1); + + for (k = 0; k < int_idx_n; k++) {///scan all integer contigs + r_a = int_a + (int_idx[k]>>32); r_n = (uint32_t)int_idx[k]; + assert(r_n >= 2); + for (i = 1; i + 1 < r_n; i++) { + ng_occ[r_a[i]]++; ng_occ[r_a[i]^1]++; + } + ng_occ[r_a[0]]++; ng_occ[r_a[r_n-1]^1]++; + } + + ma_ug_t *un_g = ma_ug_hybrid_gen(ng); ma_utg_t *u; + // print_debug_gfa(uidx->sg, un_g, uidx->uopt->coverage_cut, "iig0", uidx->uopt->sources, + // uidx->uopt->ruIndex, uidx->uopt->max_hang, uidx->uopt->min_ovlp, 0, 0, 0); + int32_t ui, un; + for (k = 0; k < un_g->u.n; k++) {///all unitigs of raw utg + u = &(un_g->u.a[k]); + zs = ze = (uint32_t)-1; un = u->n; + for (ui = 0; ui < un; ui++) { + if(occ_m(ng_occ[(u->a[ui]>>32)^1]) == 1) { + zs = ui; break; + } + } + + for (ui = ((int32_t)un)-1; ui >= 0; ui--) { + if(occ_m(ng_occ[u->a[ui]>>32]) == 1) { + ze = ui; break; + } + } + if(zs != (uint32_t)-1 && ze != (uint32_t)-1 && zs > ze) continue; + if(zs != (uint32_t)-1) ng_occ[(u->a[zs]>>32)^1] |= ((uint32_t)0x80000000); + if(ze != (uint32_t)-1) ng_occ[(u->a[ze]>>32)] |= ((uint32_t)0x80000000); + } + ma_ug_destroy(un_g); + + // fprintf(stderr, ">>>>>>[M::%s::] int_idx[0]::%lu, int_idx[1]::%lu\n", __func__, int_idx[0], int_idx[1]); + // prt_intg_info(int_idx, int_idx_n, int_a, "intg"); + for (i = b64.n = ua_n = a_n = 0; i < int_idx_n; i++) {///scan all integer contigs + s = int_idx[i]>>32; e = s + ((uint32_t)int_idx[i]); + assert(e > s + 1);//the length is at least 2 + for (k = s, bn = b64.n, pi = (uint64_t)-1; k < e; k++) { + v = int_a[k]; + if((!(ng_occ[v]&((uint32_t)0x80000000)))&&(!(ng_occ[v^1]&((uint32_t)0x80000000)))) { + continue;///it must be a unique node + } + if((pi != (uint64_t)-1) && (ng_occ[v^1]&((uint32_t)0x80000000))) { + // fprintf(stderr, "+++[M::%s::i->%u::s->%u::e->%u] pi::%lu, k::%u\n", __func__, i, s, e, pi, k); + pz = (b64.n>0)? &(b64.a[b64.n-1]):(NULL); + if((!pz) || (((*pz)>>32) != pi)) {//keep the shortest pi<->k + kv_pushp(uint64_t, b64, &pz); + *pz = pi; (*pz) <<= 32; (*pz) |= k; + a_n++; + } + } + if(ng_occ[v]&((uint32_t)0x80000000)) pi = k;//keep the shortest pi<->k + } + + b64_n = b64.n; + for (k = bn, pi = s; k < b64_n; k++) { + ei = b64.a[k]>>32; + if(pi < ei) { + kv_pushp(uint64_t, b64, &pz); ua_n++; + (*pz) = pi; (*pz) <<= 32; (*pz) |= ei; (*pz) |= ((uint64_t)0x8000000000000000); + } + pi = (uint32_t)b64.a[k]; + } + + ei = e - 1; + if(pi < ei) { + kv_pushp(uint64_t, b64, &pz); ua_n++; + (*pz) = pi; (*pz) <<= 32; (*pz) |= ei; (*pz) |= ((uint64_t)0x8000000000000000); + } + } + assert(a_n + ua_n == b64.n); + // prt_thread_info(b64.a, a_n+ua_n, NULL, 0, "tt_minus"); + // fprintf(stderr, "**0**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); + radix_sort_srt64(b64.a, b64.a + b64.n);///keeps the coordinates within int_a[] + + for (i = a_n; i < b64.n; i++) {///unavailable intervals; mask all unavailable nodes + assert(b64.a[i]&((uint64_t)0x8000000000000000)); ua_n--; + s = (b64.a[i]<<1)>>33; e = (uint32_t)b64.a[i]; assert(e > s);//[s, e] + for (k = s + 1; k < e; k++) {///note: here is [s, e] + i_idx[int_a[k]] |= ((uint64_t)0x8000000000000000); + i_idx[int_a[k]^1] |= ((uint64_t)0x8000000000000000); + } + i_idx[int_a[s]] |= ((uint64_t)0x8000000000000000); + i_idx[int_a[e]^1] |= ((uint64_t)0x8000000000000000); + } + // prt_thread_info(b64.a, a_n, NULL, 0, "tt1"); + ///unavailable intervals are useless + b64.n = a_n; assert(ua_n == 0); + for (i = 0; i < a_n; i++) { ///available intervals + s = b64.a[i]>>32; e = (uint32_t)b64.a[i]; assert(e > s);///[s, e] -> coordinates within int_a[] + assert(!(b64.a[i]&((uint64_t)0x8000000000000000))); + + for (k = s + 1; k < e; k++) {///note: here is [s, e]; s && e are unique, but [s+1, e-1] are not unique + kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[k]]++; + (*pz) = int_a[k]; (*pz) <<= 32; (*pz) |= i;///(raw unitig/non-unqiue node id)|(integer contig id) + + kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[k]^1]++; + (*pz) = int_a[k]^1; (*pz) <<= 32; (*pz) |= i; + } + kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[s]]++; + (*pz) = int_a[s]; (*pz) <<= 32; (*pz) |= i; + + kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[e]^1]++; + (*pz) = int_a[e]^1; (*pz) <<= 32; (*pz) |= i; + } + // prt_thread_info(b64.a, a_n, NULL, 0, "tt2"); + // fprintf(stderr, "**1**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); + ///index + radix_sort_srt64(b64.a + a_n, b64.a + b64.n);///(raw unitig node id)|(integer contig id) + for (k = a_n + 1, i = a_n; k <= b64.n; k++) { + if(k == b64.n || (b64.a[k]>>32) != (b64.a[i]>>32)) { + i_idx[b64.a[i]>>32] |= (((uint64_t)i)<<32)|((uint64_t)k);///b64.a[i]>>32 appear once (unique ends)/multipe times + i = k; + } + } + // prt_thread_info(b64.a, a_n, NULL, 0, "tt3"); + // fprintf(stderr, "**2**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); + ///b64.a[0, a_n]:: all resolvable paths with unique beg && end + ///b64.a[a_n, b64.n]:: (raw unitig/non-unqiue node id)|(resolvable path id) + // n_clus = usg_unique_arcs_cluster(&b64, a_n, i_idx, int_a);///this function might be wrong + n_clus = usg_unique_arcs_cluster_adv(&b64, a_n, i_idx, int_a, &ub64); + // fprintf(stderr, "**3**[M::%s::] a_n::%lu, n_clus::%lu\n", __func__, a_n, n_clus); + // prt_thread_info(b64.a, a_n, NULL, 0, "tt4"); + assert(b64.n == (a_n<<1)); + old_path_ext(uidx, ng, &b64, &ub64, a_n, n_clus, i_idx, int_a, ff, ng_occ, max_ext); + // new_path_ext(uidx, ng, &b64, &ub64, a_n, n_clus, i_idx, int_a, ff, ng_occ, max_ext);//wrong + + free(ng_occ); free(i_idx); free(ff); kv_destroy(b64); kv_destroy(ub64); + return b64.n - a_n; +} void u2g_hybrid_detan(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, asg64_v *in, asg64_v *ib) { @@ -12017,7 +12585,697 @@ void u2g_hybrid_detan(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, asg64_v * if(!in) free(tx.a); if(!ib) free(tb.a); } +void u2g_hybrid_aln(ul_resolve_t *uidx, usg_t *ng, asg64_v *ob, asg64_v *ub) +{ + uint64_t k, m, i, x, *tmp, sn; + ob->n = ub->n = 0; + for (k = 0; k < uidx->str_b.n_thread; k++) { + uidx->str_b.buf[k].res_dump.n = uidx->str_b.buf[k].u.n = uidx->str_b.buf[k].o.n = 0; + } + kt_for(uidx->str_b.n_thread, worker_integer_realign_g, uidx, uidx->uovl.i_ug->u.n); + + for (k = ob->n = ub->n = m = 0; k < uidx->str_b.n_thread; k++) { + for (i = 0; i < uidx->str_b.buf[k].res_dump.n; i++) { + x = uidx->str_b.buf[k].res_dump.a[i]; + kv_push(uint64_t, *ob, x);//aln details + if(x&((uint64_t)0x8000000000000000)) { + x -= ((uint64_t)0x8000000000000000); x >>= 32; x <<= 32;//seq_id + x |= ob->n;//offset + kv_push(uint64_t, *ub, x);///idx:: seq_id|offset_in_ob + } else { + m++; + } + } + } + + radix_sort_srt64(ub->a, ub->a + ub->n); + kv_resize(uint64_t, *ub, ub->n+m); tmp = ub->a + ub->n; m = 0; + for (k = 0; k < ub->n; k++) { + sn = ((uint32_t)(ob->a[((uint32_t)ub->a[k])-1])); + memcpy(tmp + m, ob->a + ((uint32_t)ub->a[k]), sn*sizeof((*tmp))); + ub->a[k] = m; ub->a[k] <<= 32; ub->a[k] |= sn;//offset_in_ob|occ + m += sn; + } + assert(m <= ob->n); + memcpy(ob->a, tmp, m*sizeof((*tmp))); ob->n = m; +} + +uint32_t ug_ext(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, uint32_t max_ext, +uint8_t *ff, uint32_t *ng_occ, uint64_t *i_idx, asg64_v *b64, asg64_v *ub64) +{ + uint32_t n_vtx = ng->n<<1, k, i; uint64_t pi, ei, *pz, v, b64_n; + uint64_t *r_a, r_n, zs, ze, s, e, a_n, ua_n, bn, n_clus; + memset(ff, 0, sizeof((*ff))*n_vtx); + memset(ng_occ, 0, sizeof((*ng_occ))*n_vtx); + memset(i_idx, 0, sizeof((*i_idx))*n_vtx); + b64->n = ub64->n = 0; + + for (k = 0; k < int_idx_n; k++) {///scan all integer contigs + r_a = int_a + (int_idx[k]>>32); r_n = (uint32_t)int_idx[k]; + assert(r_n >= 2); + for (i = 1; i + 1 < r_n; i++) { + ng_occ[r_a[i]]++; ng_occ[r_a[i]^1]++; + } + ng_occ[r_a[0]]++; ng_occ[r_a[r_n-1]^1]++; + } + + ma_ug_t *un_g = ma_ug_hybrid_gen(ng); + int32_t ui, un; ma_utg_t *u = NULL; + for (k = 0; k < un_g->u.n; k++) {///all unitigs of raw utg + u = &(un_g->u.a[k]); + zs = ze = (uint32_t)-1; un = u->n; + for (ui = 0; ui < un; ui++) { + if(occ_m(ng_occ[(u->a[ui]>>32)^1]) == 1) { + zs = ui; break; + } + } + + for (ui = ((int32_t)un)-1; ui >= 0; ui--) { + if(occ_m(ng_occ[u->a[ui]>>32]) == 1) { + ze = ui; break; + } + } + if(zs != (uint32_t)-1 && ze != (uint32_t)-1 && zs > ze) continue; + if(zs != (uint32_t)-1) ng_occ[(u->a[zs]>>32)^1] |= ((uint32_t)0x80000000); + if(ze != (uint32_t)-1) ng_occ[(u->a[ze]>>32)] |= ((uint32_t)0x80000000); + } + ma_ug_destroy(un_g); + + // fprintf(stderr, ">>>>>>[M::%s::] int_idx[0]::%lu, int_idx[1]::%lu\n", __func__, int_idx[0], int_idx[1]); + // prt_intg_info(int_idx, int_idx_n, int_a, "intg"); + for (i = b64->n = ua_n = a_n = 0; i < int_idx_n; i++) {///scan all integer contigs + s = int_idx[i]>>32; e = s + ((uint32_t)int_idx[i]); + assert(e > s + 1);//the length is at least 2 + for (k = s, bn = b64->n, pi = (uint64_t)-1; k < e; k++) { + v = int_a[k]; + if((!(ng_occ[v]&((uint32_t)0x80000000)))&&(!(ng_occ[v^1]&((uint32_t)0x80000000)))) { + continue;///it must be a unique node + } + if((pi != (uint64_t)-1) && (ng_occ[v^1]&((uint32_t)0x80000000))) { + pz = (b64->n>0)? &(b64->a[b64->n-1]):(NULL); + if((!pz) || (((*pz)>>32) != pi)) {//keep the shortest pi<->k + kv_pushp(uint64_t, *b64, &pz); + *pz = pi; (*pz) <<= 32; (*pz) |= k; + a_n++; + } + } + if(ng_occ[v]&((uint32_t)0x80000000)) pi = k;//keep the shortest pi<->k + } + + b64_n = b64->n; + for (k = bn, pi = s; k < b64_n; k++) { + ei = b64->a[k]>>32; + if(pi < ei) { + kv_pushp(uint64_t, *b64, &pz); ua_n++; + (*pz) = pi; (*pz) <<= 32; (*pz) |= ei; (*pz) |= ((uint64_t)0x8000000000000000); + } + pi = (uint32_t)b64->a[k]; + } + + ei = e - 1; + if(pi < ei) { + kv_pushp(uint64_t, *b64, &pz); ua_n++; + (*pz) = pi; (*pz) <<= 32; (*pz) |= ei; (*pz) |= ((uint64_t)0x8000000000000000); + } + } + assert(a_n + ua_n == b64->n); + // prt_thread_info(b64.a, a_n+ua_n, NULL, 0, "tt_minus"); + // fprintf(stderr, "**0**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); + radix_sort_srt64(b64->a, b64->a + b64->n);///keeps the coordinates within int_a[] + + + for (i = a_n; i < b64->n; i++) {///unavailable intervals; mask all unavailable nodes + assert(b64->a[i]&((uint64_t)0x8000000000000000)); ua_n--; + s = (b64->a[i]<<1)>>33; e = (uint32_t)b64->a[i]; assert(e > s);//[s, e] + for (k = s + 1; k < e; k++) {///note: here is [s, e] + i_idx[int_a[k]] |= ((uint64_t)0x8000000000000000); + i_idx[int_a[k]^1] |= ((uint64_t)0x8000000000000000); + } + i_idx[int_a[s]] |= ((uint64_t)0x8000000000000000); + i_idx[int_a[e]^1] |= ((uint64_t)0x8000000000000000); + } + // prt_thread_info(b64.a, a_n, NULL, 0, "tt1"); + ///unavailable intervals are useless + b64->n = a_n; assert(ua_n == 0); + for (i = 0; i < a_n; i++) { ///available intervals + s = b64->a[i]>>32; e = (uint32_t)b64->a[i]; assert(e > s);///[s, e] -> coordinates within int_a[] + assert(!(b64->a[i]&((uint64_t)0x8000000000000000))); + + for (k = s + 1; k < e; k++) {///note: here is [s, e]; s && e are unique, but [s+1, e-1] are not unique + kv_pushp(uint64_t, *b64, &pz); //i_idx[int_a[k]]++; + (*pz) = int_a[k]; (*pz) <<= 32; (*pz) |= i;///(raw unitig/non-unqiue node id)|(integer contig id) + + kv_pushp(uint64_t, *b64, &pz); //i_idx[int_a[k]^1]++; + (*pz) = int_a[k]^1; (*pz) <<= 32; (*pz) |= i; + } + kv_pushp(uint64_t, *b64, &pz); //i_idx[int_a[s]]++; + (*pz) = int_a[s]; (*pz) <<= 32; (*pz) |= i; + + kv_pushp(uint64_t, *b64, &pz); //i_idx[int_a[e]^1]++; + (*pz) = int_a[e]^1; (*pz) <<= 32; (*pz) |= i; + } + // prt_thread_info(b64.a, a_n, NULL, 0, "tt2"); + // fprintf(stderr, "**1**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); + ///index + radix_sort_srt64(b64->a + a_n, b64->a + b64->n);///(raw unitig node id)|(integer contig id) + for (k = a_n + 1, i = a_n; k <= b64->n; k++) { + if(k == b64->n || (b64->a[k]>>32) != (b64->a[i]>>32)) { + i_idx[b64->a[i]>>32] |= (((uint64_t)i)<<32)|((uint64_t)k);///b64.a[i]>>32 appear once (unique ends)/multipe times + i = k; + } + } + // prt_thread_info(b64.a, a_n, NULL, 0, "tt3"); + // fprintf(stderr, "**2**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); + ///b64.a[0, a_n]:: all resolvable paths with unique beg && end + ///b64.a[a_n, b64.n]:: (raw unitig/non-unqiue node id)|(resolvable path id) + // n_clus = usg_unique_arcs_cluster(&b64, a_n, i_idx, int_a);///this function might be wrong + n_clus = usg_unique_arcs_cluster_adv(b64, a_n, i_idx, int_a, ub64); + // fprintf(stderr, "**3**[M::%s::] a_n::%lu, n_clus::%lu\n", __func__, a_n, n_clus); + // prt_thread_info(b64.a, a_n, NULL, 0, "tt4"); + assert(b64->n == (a_n<<1)); + old_path_ext(uidx, ng, b64, ub64, a_n, n_clus, i_idx, int_a, ff, ng_occ, max_ext); + // new_path_ext(uidx, ng, &b64, &ub64, a_n, n_clus, i_idx, int_a, ff, ng_occ, max_ext);//wrong + return b64->n - a_n; +} + void merge_hybrid_utg_content(ma_utg_t* cc, ma_ug_t* raw, asg_t* rg, usg_t *ng, kvec_asg_arc_t_warp* edge); +void renew_usg_t_bub(ul_resolve_t *uidx, usg_t *ng, uint32_t *id_map, uint8_t *ff, uint32_t rocc_cut) +{ + ma_ug_t *ug = ma_ug_hybrid_gen(ng); + uint32_t i, k, v, p[2]; ma_utg_t *u; p[0] = 1; p[1] = 2; + memset(id_map, -1, sizeof((*id_map))*ng->n); + memset(ff, 0, sizeof((*ff))*(ng->n<<1)); + kvec_asg_arc_t_warp e; kv_init(e.a); e.i = 0; + // for (i = dn = 0; i < ng->n; i++) { + // if(ng->a[i].del) continue; + // dn++; + // } + // fprintf(stderr, "ng->n::%u, ug->u.n::%u, dn::%u\n", (uint32_t)ng->n, (uint32_t)ug->u.n, dn); + // dn = 0 + for (i = 0; i < ug->u.n; i++) { + ug->g->seq[i].c = PRIMARY_LABLE; + u = &(ug->u.a[i]); + if(u->m == 0) continue; + for (k = 0; k < u->n; k++) { + v = u->a[k]>>32; + id_map[v>>1] = i<<2; + if(k == 0) id_map[v>>1] |= p[(v^1)&1]; + if(k + 1 == u->n) id_map[v>>1] |= p[(v)&1]; + // dn++; + } + merge_hybrid_utg_content(u, uidx->l1_ug, uidx->sg, ng, &e); + ug->g->seq[i].len = u->len; + } + kv_destroy(e.a); + + destory_bubbles(uidx->bub); free(uidx->bub); CALLOC(uidx->bub, 1); + // identify_bubbles(ug, uidx->bub, uidx->r_het, NULL); + identify_bubbles_recal(uidx->sg, ug, uidx->bub, uidx->r_het, uidx->uopt->sources, uidx->uopt->ruIndex, NULL); + // fprintf(stderr, "0[M::%s::] f[51]::%u\n", __func__, ff[51]); + for (i = 0; i < ng->n; i++) { + if(id_map[i] != (uint32_t)-1) { + k = id_map[i]>>2; + // fprintf(stderr, "k::%u, ng->n::%u, ug->u.n::%u, dn::%u\n", k, (uint32_t)ng->n, (uint32_t)ug->u.n, dn); + // if(i == 7075 || i == 28174 || i == 77111 || i == 3826 || i == 12150 || i == 58312 || i == 59190 || i == 72134) { + // fprintf(stderr, "[M::%s::k->%u] utg%.6dl(%c), ug->u.a[k].n::%u, rocc_cut::%u, is_hom::%u\n", __func__, k, + // i+1, "+-"[0], (uint32_t)ug->u.a[k].n, rocc_cut, IF_HOM(k, (*(uidx->bub)))); + // } + if((ug->u.a[k].n >= rocc_cut) && (!(IF_HOM(k, (*(uidx->bub)))))) { + if(id_map[i]&p[0]) { + ff[i<<1] = 2; + // fprintf(stderr, "[M::%s::] utg%.6dl(%c), ug->u.a[k].n::%u, rocc_cut::%u\n", __func__, + // i+1, "+-"[0], (uint32_t)ug->u.a[k].n, rocc_cut); + } + if(id_map[i]&p[1]) { + ff[(i<<1)+1] = 2; + // fprintf(stderr, "[M::%s::] utg%.6dl(%c), ug->u.a[k].n::%u, rocc_cut::%u\n", __func__, + // i+1, "+-"[1], (uint32_t)ug->u.a[k].n, rocc_cut); + } + } + id_map[i] = uidx->bub->index[k]; + } + } + // fprintf(stderr, "1[M::%s::] f[51]::%u\n", __func__, ff[51]); + free(uidx->bub->index); MALLOC(uidx->bub->index, ng->n); + memcpy(uidx->bub->index, id_map, sizeof((*id_map))*ng->n); + ma_ug_destroy(ug); +} + +void gen_hybrid_aln_idx(usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, uint64_t int_an, asg64_v *b64, uint64_t *ridx) +{ + uint64_t k, i, l, m, v, s, e, *a, a_n; + kv_resize(uint64_t, *b64, int_an); b64->n = int_an; + memset(ridx, 0, sizeof((*ridx))*ng->n); + for (k = 0; k < int_idx_n; k++) {///scan all integer contigs + s = int_idx[k]>>32; e = s + ((uint32_t)int_idx[k]); assert(e > s + 1);//the length is at least 2 + for (i = s; i < e; i++) ridx[int_a[i]>>1]++; + } + + for (k = l = 0; k < ng->n; k++) { + m = ridx[k]; + ridx[k] = l; ridx[k] <<= 32; ridx[k] |= m; + l += m; + } + + for (k = 0; k < ng->n; k++) { + a = b64->a + (ridx[k]>>32); a_n = (uint32_t)ridx[k]; + if(a_n) a[a_n-1] = 0; + } + + for (k = 0; k < int_idx_n; k++) { + s = int_idx[k]>>32; e = s + ((uint32_t)int_idx[k]); assert(e > s + 1);//the length is at least 2 + for (i = s; i < e; i++) { + v = int_a[i]>>1; + a = b64->a + (ridx[v]>>32); + a_n = (uint32_t)ridx[v]; + if(a_n) { + if(a[a_n-1] == a_n-1) { + a[a_n-1] = (k<<32)|i; + } else { + a[a[a_n-1]++] = (k<<32)|i; + } + } + } + } +} + +void gen_sub_integer_path(uint64_t *int_a, int64_t s, int64_t e, int64_t it, uint64_t v, asg64_v *res) +{ + int64_t k; res->n = 0; assert((int_a[it]>>1) == (v>>1)); + if(int_a[it] == v) {///forward + kv_resize(uint64_t, *res, (uint64_t)(e - it)); + for (k = it; k < e; k++) res->a[res->n++] = int_a[k]; + } else {//reverse + kv_resize(uint64_t, *res, (uint64_t)(it - s)); + for (k = it; k >= s; k--) res->a[res->n++] = int_a[k]^1; + } +} + +void prt_sub_integer_path(asg64_v *res, uint64_t it, uint64_t w0, uint64_t w1, uint64_t z_n, uint64_t z) +{ + uint64_t k; + fprintf(stderr, "(it::%lu) res->n::%u, w0::%lu, w1::%lu, z_n::%lu, z::%lu\n", it, (uint32_t)res->n, + w0, w1, z_n, z); + for (k = 0; k < res->n; k++) { + fprintf(stderr, "(it::%lu) utg%.6dl,", it, (int32_t)(res->a[k]>>1)+1); + } + fprintf(stderr, "\n"); +} + +uint64_t is_best_path(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t *int_a, +uint64_t s, uint64_t e, uint64_t it, uint64_t v, uint64_t *ridx_a, uint64_t *ridx, +asg64_v *b0, asg64_v *b1, double cutoff) +{ + uint64_t *arc_a, arc_n, k, is, ie, z, zn, w0, w1, min_w0, min_w1, alt_n = 0, is_contain = 1; + gen_sub_integer_path(int_a, s, e, it, v, b0); + + arc_a = ridx_a + (ridx[v>>1]>>32); + arc_n = (uint32_t)ridx[v>>1]; + + + // uint64_t is_debug = 0; + // if(((v>>1) == 10600)/** && (b0->n > 1 && (b0->a[b0->n-1]>>1) == 34944)**/) { + // is_debug = 1; + // fprintf(stderr, "[M::%s::] utg%.6dl(%c)\n", __func__, + // (int32_t)(v>>1)+1, "+-"[v&1]); + // prt_sub_integer_path(b0, it, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1); + // } + + + for (k = 0, min_w0 = min_w1 = (uint64_t)-1, is_contain = 1; k < arc_n; k++) { + if(it == ((uint32_t)arc_a[k])) continue; + is = int_idx[arc_a[k]>>32]>>32; alt_n++; + ie = is + ((uint32_t)(int_idx[arc_a[k]>>32])); + gen_sub_integer_path(int_a, is, ie, ((uint32_t)arc_a[k]), v, b1); + zn = MIN(b0->n, b1->n); + for (z = 0; z < zn && b0->a[z] == b1->a[z]; z++); ///z: first raw unitig that is different between two paths + assert(z > 0); w0 = w1 = (uint64_t)-1; + if(z < b0->n) get_integer_seq_ovlps(uidx, b0->a, b0->n, z - 1, 0, NULL, &w0); + if(z < b1->n) get_integer_seq_ovlps(uidx, b1->a, b1->n, z - 1, 0, NULL, &w1); + // if(is_debug) { + // prt_sub_integer_path(b0, it, w0, w1, zn, z); + // prt_sub_integer_path(b1, it, w0, w1, zn, z); + // } + if(w0 == (uint64_t)-1) w0 = 0; + if(w1 == (uint64_t)-1) w1 = 0; + if(b0->n == zn) return 0;///b0 is shorter + if(b1->n == zn) continue;///b1 is shorter + if((min_w0 == (uint64_t)-1) || (z == zn) || (min_w0 > w0) || (min_w0 == w0 && min_w1 < w1)) { + min_w0 = w0; min_w1 = w1; + } + is_contain = 0; + } + // fprintf(stderr, "[M::%s::] utg%.6dl(%c)->utg%.6dl(%c)\n", __func__, + // (int32_t)(v>>1)+1, "+-"[v&1], (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1]); + if(alt_n == 0) return 1; + if(is_contain) return 1; + + if(min_w0 == (uint64_t)-1) min_w0 = 0; + if(min_w1 == (uint64_t)-1) min_w1 = 0; + if((min_w0 > min_w1) && (min_w1 <= (min_w0*cutoff))) return 1; + return 0; +} + + +void get_best_path(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t *int_a, +uint8_t *f, uint64_t s, uint64_t e, asg64_v *b0, asg64_v *b1, uint64_t *ridx_a, +uint64_t *ridx, asg64_v *res) +{ + uint64_t pi, v, k, res_n = res->n, *pz; + b0->n = b1->n = 0; + for (k = s, pi = (uint64_t)-1; k < e; k++) { + v = int_a[k]; + if((!f[v])&&(!f[v^1])) continue; + if((pi != (uint64_t)-1) && (f[v^1])) { + // fprintf(stderr, "+[M::%s::] utg%.6dl(%c), f[v^1]::%u ,v^1::%lu\n", + // __func__, (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1], f[v^1], v^1); + if(is_best_path(uidx, ng, int_idx, int_a, s, e, k, v^1, ridx_a, ridx, b0, b1, 0.51)) { + // fprintf(stderr, "[M::%s::] utg%.6dl(%c)->utg%.6dl(%c)\n", __func__, + // (int32_t)(int_a[pi]>>1)+1, "+-"[int_a[pi]&1], (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1]); + pz = (res->n > res_n)? &(res->a[res->n-1]):(NULL); + if((!pz) || (((*pz)>>32) != pi)) {//keep the shortest pi<->k + kv_pushp(uint64_t, *res, &pz); + *pz = pi; (*pz) <<= 32; (*pz) |= k; + } + } else { + pi = (uint64_t)-1; + } + } + if(f[v]) { + // fprintf(stderr, "-[M::%s::] utg%.6dl(%c), f[v]::%u, v::%lu\n", + // __func__, (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1], f[v], v); + if(is_best_path(uidx, ng, int_idx, int_a, s, e, k, v, ridx_a, ridx, b0, b1, 0.51)) { + pi = k;//keep the shortest pi<->k + } else { + pi = (uint64_t)-1; + } + } + } +} + +static void worker_ul_aln_path(void *data, long i, int tid) // callback for kt_for() +{ + unique_bridge_check_t *u_aux = (unique_bridge_check_t*)data; + ul_resolve_t *uidx = u_aux->uidx; + integer_t *buf = &(uidx->str_b.buf[tid]); + uint64_t s, e; + // uint64_t *x = &(uidx->uovl.iug_tra->a[i]); + // asg_arc_t *ve = &(uidx->uovl.i_ug->g->arc[*x]); + + asg64_v b_v, b_r, res; + b_v.a = buf->u.a; b_v.n = buf->u.n; b_v.m = buf->u.m; + b_r.a = buf->o.a; b_r.n = buf->o.n; b_r.m = buf->o.m; + res.a = buf->res_dump.a; res.n = buf->res_dump.n; res.m = buf->res_dump.m; + + b_v.n = b_r.n = 0; + s = u_aux->int_idx[i]>>32; ///the i-th integer contig/path + e = s + ((uint32_t)(u_aux->int_idx[i])); + assert(e > s + 1);//the length is at least 2 + + get_best_path(uidx, u_aux->ng, u_aux->int_idx, u_aux->int_a, u_aux->f, s, e, &b_v, &b_r, u_aux->ridx_a, u_aux->ridx, &res); + // get_ul_arc_supports(uidx, ve, &b_v, &b_r, 1, &w_v, &w_r); + + buf->u.a = b_v.a; buf->u.n = b_v.n; buf->u.m = b_v.m; + buf->o.a = b_r.a; buf->o.n = b_r.n; buf->o.m = b_r.m; + buf->res_dump.a = res.a; buf->res_dump.n = res.n; buf->res_dump.m = res.m; +} + +uint32_t ug_ext_0(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, uint32_t max_ext, +uint8_t *ff, uint32_t *ng_occ, uint64_t *i_idx, asg64_v *b64, asg64_v *ub64, uint64_t a_n) +{ + uint32_t n_vtx = ng->n<<1, k, i; uint64_t *pz; + uint64_t s, e, n_clus; + memset(ff, 0, sizeof((*ff))*n_vtx); + // memset(ng_occ, 0, sizeof((*ng_occ))*n_vtx); + memset(i_idx, 0, sizeof((*i_idx))*n_vtx); + ub64->n = 0; + + for (i = 0; i < a_n; i++) { ///available intervals + s = b64->a[i]>>32; e = (uint32_t)b64->a[i]; assert(e > s);///[s, e] -> coordinates within int_a[] + assert(!(b64->a[i]&((uint64_t)0x8000000000000000))); + // fprintf(stderr, "\n[M::%s::] occ::%lu\n", __func__, e - s); + // for (k = s; k <= e; k++) { + // // fprintf(stderr, "utg%.6dl(%c),", (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1]); + // fprintf(stderr, "utg%.6dl,", (int32_t)(int_a[k]>>1)+1); + // } + // fprintf(stderr, "\n"); + + for (k = s + 1; k < e; k++) {///note: here is [s, e]; s && e are unique, but [s+1, e-1] are not unique + kv_pushp(uint64_t, *b64, &pz); //i_idx[int_a[k]]++; + (*pz) = int_a[k]; (*pz) <<= 32; (*pz) |= i;///(raw unitig/non-unqiue node id)|(integer contig id) + + kv_pushp(uint64_t, *b64, &pz); //i_idx[int_a[k]^1]++; + (*pz) = int_a[k]^1; (*pz) <<= 32; (*pz) |= i; + } + kv_pushp(uint64_t, *b64, &pz); //i_idx[int_a[s]]++; + (*pz) = int_a[s]; (*pz) <<= 32; (*pz) |= i; + + kv_pushp(uint64_t, *b64, &pz); //i_idx[int_a[e]^1]++; + (*pz) = int_a[e]^1; (*pz) <<= 32; (*pz) |= i; + } + // prt_thread_info(b64.a, a_n, NULL, 0, "tt2"); + // fprintf(stderr, "**1**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); + ///index + radix_sort_srt64(b64->a + a_n, b64->a + b64->n);///(raw unitig node id)|(integer contig id) + for (k = a_n + 1, i = a_n; k <= b64->n; k++) { + if(k == b64->n || (b64->a[k]>>32) != (b64->a[i]>>32)) { + i_idx[b64->a[i]>>32] |= (((uint64_t)i)<<32)|((uint64_t)k);///b64.a[i]>>32 appear once (unique ends)/multipe times + i = k; + } + } + + n_clus = usg_unique_arcs_cluster_adv(b64, a_n, i_idx, int_a, ub64); + // fprintf(stderr, "**3**[M::%s::] a_n::%lu, n_clus::%lu\n", __func__, a_n, n_clus); + // prt_thread_info(b64.a, a_n, NULL, 0, "tt4"); + assert(b64->n == (a_n<<1)); + + // new_path_ext(uidx, ng, &b64, &ub64, a_n, n_clus, i_idx, int_a, ff, ng_occ, max_ext);//wrong + return old_path_ext(uidx, ng, b64, ub64, a_n, n_clus, i_idx, int_a, ff, ng_occ, max_ext); +} + +void debug_prt_renew_aln(usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, uint64_t int_an, asg64_v *b64, uint64_t *ridx) +{ + uint64_t k, s, e, i; + for (i = 0; i < int_idx_n; i++) {///scan all integer contigs + s = int_idx[i]>>32; e = s + ((uint32_t)int_idx[i]); assert(e > s + 1);//the length is at least 2 + fprintf(stderr, "\n[M::%s::] occ::%lu\n", __func__, e - s); + for (k = s; k <= e; k++) { + // fprintf(stderr, "utg%.6dl(%c),", (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1]); + fprintf(stderr, "utg%.6dl,", (int32_t)(int_a[k]>>1)+1); + } + fprintf(stderr, "\n"); + } +} + +uint32_t ug_ext_free(asg64_v *ob, asg64_v *ub, ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, +uint8_t **ff, uint32_t **ng_occ, uint64_t **i_idx, asg64_v *b64, asg64_v *ub64, uint32_t rocc_cut) +{ + fprintf(stderr, "[M::%s::] rocc_cut::%u\n", __func__, rocc_cut); + uint32_t k, n_vtx = ng->n<<1, a_n; unique_bridge_check_t u_aux; + REALLOC((*ff), n_vtx); REALLOC((*ng_occ), n_vtx); REALLOC((*i_idx), n_vtx); + renew_usg_t_bub(uidx, ng, *ng_occ, *ff, rocc_cut); + u2g_hybrid_aln(uidx, ng, ob, ub); + gen_hybrid_aln_idx(ng, ub->a, ub->n, ob->a, ob->n, b64, *i_idx); + // if(rocc_cut == 10) + // { + // fprintf(stderr, "\n[M::%s::] rocc_cut::%u\n", __func__, rocc_cut); + // debug_prt_renew_aln(ng, ub->a, ub->n, ob->a, ob->n, b64, *i_idx); + // } + + u_aux.f = *ff; u_aux.ng = ng; u_aux.int_idx = ub->a; u_aux.int_idx_n = ub->n; u_aux.uidx = uidx; + u_aux.int_a = ob->a; u_aux.int_an = ob->n; u_aux.ridx_a = b64->a; u_aux.ridx = *i_idx; + for (k = 0; k < uidx->str_b.n_thread; k++) { + uidx->str_b.buf[k].res_dump.n = uidx->str_b.buf[k].u.n = uidx->str_b.buf[k].o.n = 0; + } + kt_for(uidx->str_b.n_thread, worker_ul_aln_path, &u_aux, u_aux.int_idx_n); + for (k = b64->n = a_n = 0; k < uidx->str_b.n_thread; k++) { + a_n += uidx->str_b.buf[k].res_dump.n; + kv_resize(uint64_t, *b64, a_n); + memcpy(b64->a+b64->n, uidx->str_b.buf[k].res_dump.a, uidx->str_b.buf[k].res_dump.n*(sizeof(*(b64->a)))); + b64->n = a_n; + } + radix_sort_srt64(b64->a, b64->a + b64->n);///keeps the coordinates within int_a[] + // fprintf(stderr, "[M::%s::] a_n::%u, int_idx_n::%u, int_n::%u\n", + // __func__, a_n, (uint32_t)ub->n, (uint32_t)ob->n); + a_n = ug_ext_0(uidx, ng, ub->a, ub->n, ob->a, max_ext, *ff, *ng_occ, *i_idx, b64, ub64, a_n); + // if(a_n) usg_cleanup(ng); + return a_n; +} + +uint32_t ug_ext_strict(asg64_v *ob, asg64_v *ub, ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, +uint8_t **ff, uint32_t **ng_occ, uint64_t **i_idx, asg64_v *b64, asg64_v *ub64) +{ + fprintf(stderr, "[M::%s::]\n", __func__); + uint32_t k, i, n_vtx = ng->n<<1, a_n; + REALLOC((*ff), n_vtx); REALLOC((*ng_occ), n_vtx); REALLOC((*i_idx), n_vtx); + u2g_hybrid_aln(uidx, ng, ob, ub); b64->n = ub64->n = 0; + uint64_t *int_idx = ub->a, int_idx_n = ub->n, *int_a = ob->a; + uint64_t *r_a, r_n, zs, ze, s, e, pi, v, *pz; + + memset((*ng_occ), 0, sizeof((*(*ng_occ)))*n_vtx); + for (k = 0; k < int_idx_n; k++) {///scan all integer contigs + r_a = int_a + (int_idx[k]>>32); r_n = (uint32_t)int_idx[k]; + assert(r_n >= 2); + for (i = 1; i + 1 < r_n; i++) { + (*ng_occ)[r_a[i]]++; (*ng_occ)[r_a[i]^1]++; + } + (*ng_occ)[r_a[0]]++; (*ng_occ)[r_a[r_n-1]^1]++; + } + + ma_ug_t *un_g = ma_ug_hybrid_gen(ng); + int32_t ui, un; ma_utg_t *u = NULL; + for (k = 0; k < un_g->u.n; k++) {///all unitigs of raw utg + u = &(un_g->u.a[k]); + zs = ze = (uint32_t)-1; un = u->n; + for (ui = 0; ui < un; ui++) { + if(occ_m((*ng_occ)[(u->a[ui]>>32)^1]) == 1) { + zs = ui; break; + } + } + + for (ui = ((int32_t)un)-1; ui >= 0; ui--) { + if(occ_m((*ng_occ)[u->a[ui]>>32]) == 1) { + ze = ui; break; + } + } + if(zs != (uint32_t)-1 && ze != (uint32_t)-1 && zs > ze) continue; + if(zs != (uint32_t)-1) (*ng_occ)[(u->a[zs]>>32)^1] |= ((uint32_t)0x80000000); + if(ze != (uint32_t)-1) (*ng_occ)[(u->a[ze]>>32)] |= ((uint32_t)0x80000000); + } + ma_ug_destroy(un_g); + + // fprintf(stderr, ">>>>>>[M::%s::] int_idx[0]::%lu, int_idx[1]::%lu\n", __func__, int_idx[0], int_idx[1]); + // prt_intg_info(int_idx, int_idx_n, int_a, "intg"); + for (i = b64->n = a_n = 0; i < int_idx_n; i++) {///scan all integer contigs + s = int_idx[i]>>32; e = s + ((uint32_t)int_idx[i]); + assert(e > s + 1);//the length is at least 2 + for (k = s, pi = (uint64_t)-1; k < e; k++) { + v = int_a[k]; + if((!((*ng_occ)[v]&((uint32_t)0x80000000)))&&(!((*ng_occ)[v^1]&((uint32_t)0x80000000)))) { + continue;///it must be a unique node + } + if((pi != (uint64_t)-1) && ((*ng_occ)[v^1]&((uint32_t)0x80000000))) { + pz = (b64->n>0)? &(b64->a[b64->n-1]):(NULL); + if((!pz) || (((*pz)>>32) != pi)) {//keep the shortest pi<->k + kv_pushp(uint64_t, *b64, &pz); + *pz = pi; (*pz) <<= 32; (*pz) |= k; + a_n++; + } + } + if((*ng_occ)[v]&((uint32_t)0x80000000)) pi = k;//keep the shortest pi<->k + } + } + // prt_thread_info(b64.a, a_n+ua_n, NULL, 0, "tt_minus"); + // fprintf(stderr, "**0**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); + radix_sort_srt64(b64->a, b64->a + b64->n);///keeps the coordinates within int_a[] + + a_n = ug_ext_0(uidx, ng, ub->a, ub->n, ob->a, max_ext, *ff, *ng_occ, *i_idx, b64, ub64, a_n); + // if(a_n) usg_cleanup(ng); + return a_n; +} + +void u2g_hybrid_detan_iter(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, uint32_t clean_round, asg64_v *in, asg64_v *ib) +{ + uint32_t n_vtx = ng->n<<1, ncut = 0, k; + uint8_t *ff; CALLOC(ff, n_vtx); + uint32_t *ng_occ; CALLOC(ng_occ, n_vtx); + uint64_t *i_idx; CALLOC(i_idx, n_vtx); + asg64_v b64, ub64; kv_init(b64); kv_init(ub64); + asg64_v tx = {0,0,0}, tb = {0,0,0}, *ob = NULL, *ub = NULL; + ob = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); ob->n = ub->n = 0; + + for (k = 0; k < clean_round; k++) { + ncut += ug_ext_strict(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64); + // prt_usg_t(uidx, ng, "ng0"); + ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 48); + ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 16); + } + // // ug_ext_strict(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64); + // ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 50); + // // prt_usg_t(uidx, ng, "ng0"); + // ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 10); + // ncut += ug_ext_strict(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64); + if(ncut) { + usg_cleanup(ng); usg_arc_cut_tips(ng, max_ext, 1, ub); + } + + // u2g_hybrid_aln(uidx, ng, ob, ub); + // if(ug_ext(uidx, ng, ub->a, ub->n, ob->a, max_ext, ff, ng_occ, i_idx, &b64, &ub64)) { + // // debug_sysm_usg_t(ng, __func__); + // usg_cleanup(ng); + // // debug_sysm_usg_t(ng, __func__); + // } + + // prt_usg_t(uidx, ng, "ng0"); + // if(ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 50)) { + // usg_cleanup(ng); + // } + // prt_usg_t(uidx, ng, "ng1"); + + if(!in) free(tx.a); if(!ib) free(tb.a); + free(ng_occ); free(i_idx); free(ff); kv_destroy(b64); kv_destroy(ub64); +} + + +/** + +void u2g_hybrid_detan(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, asg64_v *in, asg64_v *ib) +{ + uint64_t k, i, x, m, *tmp, sn; asg64_v tx = {0,0,0}, tb = {0,0,0}, *ob = NULL, *ub = NULL; + ob = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); ob->n = ub->n = 0; + + for (k = 0; k < uidx->str_b.n_thread; k++) { + uidx->str_b.buf[k].res_dump.n = uidx->str_b.buf[k].u.n = uidx->str_b.buf[k].o.n = 0; + } + kt_for(uidx->str_b.n_thread, worker_integer_realign_g, uidx, uidx->uovl.i_ug->u.n); + + for (k = ob->n = ub->n = m = 0; k < uidx->str_b.n_thread; k++) { + for (i = 0; i < uidx->str_b.buf[k].res_dump.n; i++) { + x = uidx->str_b.buf[k].res_dump.a[i]; + kv_push(uint64_t, *ob, x);//aln details + if(x&((uint64_t)0x8000000000000000)) { + x -= ((uint64_t)0x8000000000000000); x >>= 32; x <<= 32;//seq_id + x |= ob->n;//offset + kv_push(uint64_t, *ub, x);///idx:: seq_id|offset_in_ob + } else { + m++; + } + } + } + + radix_sort_srt64(ub->a, ub->a + ub->n); + kv_resize(uint64_t, *ub, ub->n+m); tmp = ub->a + ub->n; m = 0; + for (k = 0; k < ub->n; k++) { + sn = ((uint32_t)(ob->a[((uint32_t)ub->a[k])-1])); + memcpy(tmp + m, ob->a + ((uint32_t)ub->a[k]), sn*sizeof((*tmp))); + ub->a[k] = m; ub->a[k] <<= 32; ub->a[k] |= sn;//offset_in_ob|occ + m += sn; + } + assert(m <= ob->n); + memcpy(ob->a, tmp, m*sizeof((*tmp))); ob->n = m; + + debug_sysm_usg_t(ng, __func__); + // prt_usg_t(uidx, ng, "ng4"); + if(gen_unique_g_adv(uidx, ng, ub->a, ub->n, ob->a, max_ext)) { + // usg_arc_t *z = get_usg_arc(ng, 2, 576), *q = get_usg_arc(ng, 577, 3); + // fprintf(stderr, "xxxx0xxx[M::%s::] p->del::%u, q->del::%u\n", + // __func__, z?z->del:1, q?q->del:1); + ///debug + debug_sysm_usg_t(ng, __func__); + // z = get_usg_arc(ng, 2, 576); q = get_usg_arc(ng, 577, 3); + // fprintf(stderr, "xxxx1xxx[M::%s::] p->del::%u, q->del::%u\n", + // __func__, z?z->del:1, q?q->del:1); + // fprintf(stderr, "+[M::%s::] ng->n::%u\n", __func__, (uint32_t)ng->n); + usg_cleanup(ng); + ///debug + debug_sysm_usg_t(ng, __func__); + // fprintf(stderr, "-[M::%s::] ng->n::%u\n", __func__, (uint32_t)ng->n); + } + // prt_usg_t(uidx, ng, "ng_dbg"); + if(!in) free(tx.a); if(!ib) free(tb.a); +} +**/ + ma_ug_t *gen_debug_hybrid_ug(ul_resolve_t *uidx, usg_t *ng) { ma_ug_t *ug = NULL; uint32_t k, nv, z; usg_arc_t *av; asg_arc_t *p; @@ -12084,8 +13342,8 @@ void u2g_hybrid_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, usg_t *ng, asg64_v * for (ss = 1; ss <= 1/**6**/; ss++) { mm_tip = ulopt->max_tip_hifi*ss; - fprintf(stderr, "\n[M::%s::] ss::%ld, mm_tip::%ld, ulopt->clean_round::%ld\n", - __func__, ss, mm_tip, ulopt->clean_round); + // fprintf(stderr, "\n[M::%s::] ss::%ld, mm_tip::%ld, ulopt->clean_round::%ld\n", + // __func__, ss, mm_tip, ulopt->clean_round); for (i = 0, drop = ulopt->min_ovlp_drop_ratio; i < ulopt->clean_round; i++, drop += step) { if(drop > ulopt->max_ovlp_drop_ratio) drop = ulopt->max_ovlp_drop_ratio; // fprintf(stderr, "-0-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop); @@ -12135,11 +13393,12 @@ void u2g_hybrid_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, usg_t *ng, asg64_v * debug_sysm_usg_t(ng, __func__); /******for debug******/ - prt_usg_t(uidx, ng, "ng_dbg"); + // prt_usg_t(uidx, ng, "ng_dbg"); /******for debug******/ // u2g_hybrid_extend(ng, NULL, b, ub); - u2g_hybrid_detan(uidx, ng, mm_tip, b, ub); + // u2g_hybrid_detan(uidx, ng, mm_tip, b, ub); + u2g_hybrid_detan_iter(uidx, ng, mm_tip, ulopt->clean_round, b, ub); } void merge_hybrid_utg_content(ma_utg_t* cc, ma_ug_t* raw, asg_t* rg, usg_t *ng, kvec_asg_arc_t_warp* edge) @@ -12879,12 +14138,15 @@ ul2ul_idx_t *gen_ul2ul(ul_resolve_t *uidx, ug_opt_t *uopt, ulg_opt_t *ulopt, uin kt_for(uidx->str_b.n_thread, worker_integer_postprecess, uidx, uls->n); gen_integer_normalize(uidx); + chimeric_integer_deal(uidx); append_utg_es(uidx); + + kt_for(uidx->str_b.n_thread, worker_integert_debug_sym, uidx, z->tot);///all ul + ug print_integert_ovlp_stat(z); - remove_integert_containment(uidx, keep_raw_utg); + kt_for(uidx->str_b.n_thread, worker_integert_debug_sym, uidx, z->tot);///all ul + ug print_integert_ovlp_stat(z); print_uls_seq(uidx, asm_opt.output_file_name); @@ -12954,8 +14216,63 @@ double max_ovlp_drop_ratio, double hom_check_drop_rate, int64_t max_tip, int64_t z->is_trio = is_trio; } +ma_ug_t* output_trio_unitig_graph_ul(ug_opt_t *uopt, ul_resolve_t *uidx, char* ou, uint8_t flag) +{ + char* gfa_name; MALLOC(gfa_name, strlen(ou)+100); + sprintf(gfa_name, "%s.%s.p_ctg.gfa", ou, (flag==FATHER?"hap1":"hap2")); + FILE* output_file = fopen(gfa_name, "w"); + + ma_ug_t *ug = copy_untig_graph(uidx->uovl.hybrid_ug); + kvec_asg_arc_t_warp ne; kv_init(ne.a); + + adjust_utg_by_trio(&ug, uidx->sg, flag, TRIO_THRES, uopt->sources, uopt->reverse_sources, + uopt->coverage_cut, uopt->tipsLen, uopt->tip_drop_ratio, uopt->stops_threshold, uopt->ruIndex, + uopt->chimeric_rate, uopt->drop_ratio, uopt->max_hang, uopt->min_ovlp, &ne, uopt->b_mask_t); + + // if(asm_opt.b_low_cov > 0) { + // break_ug_contig(&ug, uidx->sg, &R_INF, uopt->coverage_cut, uopt->sources, uopt->ruIndex, &ne, + // uopt->max_hang, uopt->min_ovlp, &asm_opt.b_low_cov, NULL, asm_opt.m_rate); + // } + + // if(asm_opt.b_high_cov > 0) + // { + // break_ug_contig(&ug, uidx->sg, &R_INF, uopt->coverage_cut, uopt->sources, uopt->ruIndex, &ne, + // uopt->max_hang, uopt->min_ovlp, NULL, &asm_opt.b_high_cov, asm_opt.m_rate); + // } + + fprintf(stderr, "Writing %s to disk... \n", gfa_name); + ma_ug_seq(ug, uidx->sg, uopt->coverage_cut, uopt->sources, &ne, uopt->max_hang, uopt->min_ovlp, 0, 1); + + ma_ug_print(ug, uidx->sg, uopt->coverage_cut, uopt->sources, uopt->ruIndex, (flag==FATHER?"h1tg":"h2tg"), output_file); + fclose(output_file); + + sprintf(gfa_name, "%s.%s.p_ctg.noseq.gfa", ou, (flag==FATHER?"hap1":"hap2")); + output_file = fopen(gfa_name, "w"); + ma_ug_print_simple(ug, uidx->sg, uopt->coverage_cut, uopt->sources, uopt->ruIndex, (flag==FATHER?"h1tg":"h2tg"), output_file); + fclose(output_file); + // if(asm_opt.bed_inconsist_rate != 0) + // { + // sprintf(gfa_name, "%s.%s.p_ctg.lowQ.bed", output_file_name, f_prefix?f_prefix:(flag==FATHER?"hap1":"hap2")); + // output_file = fopen(gfa_name, "w"); + // ma_ug_print_bed(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, + // max_hang, min_ovlp, asm_opt.bed_inconsist_rate, (flag==FATHER?"h1tg":"h2tg"), output_file, NULL); + // fclose(output_file); + // } + + free(gfa_name); + ma_ug_destroy(ug); + kv_destroy(ne.a); + return NULL; +} + +void gen_ul_trio_graph(ug_opt_t *uopt, ul_resolve_t *uidx, char *o_file) +{ + output_trio_unitig_graph_ul(uopt, uidx, o_file, FATHER); + output_trio_unitig_graph_ul(uopt, uidx, o_file, MOTHER); +} + void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio, -double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio) +double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file) { uint64_t i; uint8_t *r_het = NULL; bubble_type *bub = NULL; ulg_opt_t uu; for (i = 0; i < sg->n_seq; ++i) { @@ -12976,7 +14293,7 @@ double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_ print_raw_uls_seq(uidx, asm_opt.output_file_name); ul_re_correct(uidx, 3); init_ulg_opt_t(&uu, uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.55, max_tip, max_tip<<1, b_mask_t, is_trio); - print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug0", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1); + // print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug0", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1); /**ul2ul_idx_t *u2o = **/gen_ul2ul(uidx, uopt, &uu, 0); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-3"); @@ -12985,7 +14302,12 @@ double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_ // resolve_dip_bub_chains(uidx); // free(r_het); destory_bubbles(bub); free(bub); + print_debug_gfa(sg, uidx->uovl.hybrid_ug, uopt->coverage_cut, "hybrid_ug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1); print_debug_gfa(sg, uidx->uovl.hybrid_ug, uopt->coverage_cut, "hybrid_ug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0); print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 0); + + if(is_trio) gen_ul_trio_graph(uopt, uidx, o_file); + + exit(0); } \ No newline at end of file diff --git a/gfa_ut.h b/gfa_ut.h index 7304c1a..761cc71 100644 --- a/gfa_ut.h +++ b/gfa_ut.h @@ -24,7 +24,7 @@ void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float o uint32_t asg_cut_large_indel(asg_t *g, asg64_v *in, int32_t max_ext, float ou_rat, uint32_t is_ou); uint32_t asg_cut_semi_circ(asg_t *g, uint32_t lim_len, uint32_t is_clean); void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg, int64_t clean_round, double min_ovlp_drop_ratio, -double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio); +double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_trio, char *o_file); void recover_contain_g(asg_t *g, ma_hit_t_alloc *src, R_to_U* ruIndex, int64_t max_hang, int64_t min_ovlp, int64_t ul_occ); void normalize_gou(asg_t *g); void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd); diff --git a/hic.cpp b/hic.cpp index cc06ae0..a1d51ea 100644 --- a/hic.cpp +++ b/hic.cpp @@ -2519,6 +2519,268 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_t // fprintf(stderr, "-bub->index[18759]: %u, bub->num.n: %u\n", (uint32_t)bub->index[18759], bub->num.n); } +uint32_t get_unitig_het_fly(ma_ug_t* ug, uint32_t uid, asg_t* sg, int64_t het_cov_thres, +ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag, uint32_t m_het_occ, uint32_t m_het_label, +uint32_t p_het_label, uint32_t n_het_label) +{ + ma_utg_t *u = &(ug->u.a[uid]); + uint32_t k, i, j, rId, nv, tn, is_Unitig; + asg_arc_t *av = NULL; ma_hit_t *h; + int64_t R_bases = 0, C_bases = 0, cov; + + ///set + u = &(ug->u.a[uid]); + for (k = 0; k < u->n; k++) { + rId = u->a[k]>>33; + r_flag[rId] = 1; + } + for (i = 0; i < 2; i++) { + nv = asg_arc_n(ug->g, (uid<<1)+i); + av = asg_arc_a(ug->g, (uid<<1)+i); + for (j = 0; j < nv; j++) { + u = &(ug->u.a[av[j].v>>1]); + for (k = 0; k < u->n; k++) { + rId = u->a[k]>>33; + r_flag[rId] = 2; + } + } + } + + u = &(ug->u.a[uid]); + for (k = 0; k < u->n; k++) { + if(u->a[k] == (uint64_t)-1) continue; + rId = u->a[k]>>33; + R_bases += sg->seq[rId].len; + for (j = 0; j < (uint64_t)(sources[rId].length); j++) { + h = &(sources[rId].buffer[j]); + if(h->el != 1) continue; + tn = Get_tn((*h)); + if(sg->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 || sg->seq[tn].del == 1) continue; + } + if(sg->seq[tn].del == 1) continue; + if(r_flag[tn] == 0) continue; + if(r_flag[tn] == 1) { + C_bases += (Get_qe((*h)) - Get_qs((*h))); + } + if(r_flag[tn] == 2) { + C_bases += ((Get_qe((*h)) - Get_qs((*h)))/2); + } + } + } + + ///reset + u = &(ug->u.a[uid]); + for (k = 0; k < u->n; k++) { + rId = u->a[k]>>33; + r_flag[rId] = 0; + } + for (i = 0; i < 2; i++) { + nv = asg_arc_n(ug->g, (uid<<1)+i); + av = asg_arc_a(ug->g, (uid<<1)+i); + for (j = 0; j < nv; j++) { + u = &(ug->u.a[av[j].v>>1]); + for (k = 0; k < u->n; k++) { + rId = u->a[k]>>33; + r_flag[rId] = 0; + } + } + } + + u = &(ug->u.a[uid]); cov = 0; + if(R_bases > 0) cov = C_bases/R_bases; + // if(uid == 9253) { + // fprintf(stderr, "[M::%s::uid->%u] u->n::%u, C_bases::%ld, R_bases::%ld, het_cov_thres::%ld, m_het_occ::%u\n", + // __func__, uid, (uint32_t)u->n, C_bases, R_bases, het_cov_thres, m_het_occ); + // } + if((cov <= (het_cov_thres*1.333333)) && (u->n >= m_het_occ)) return m_het_label; ///must het + if((cov >= (het_cov_thres*1.6))) return n_het_label; ///hom + return p_het_label; ///potential het +} + +void identify_bubbles_recal(asg_t* sg, ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, ma_hit_t_alloc* sources, R_to_U* ruIndex, +kv_u_trans_t *ref) +{ + asg_cleanup(ug->g); + if (!ug->g->is_symm) asg_symm(ug->g); + uint32_t v, n_vtx = ug->g->n_seq * 2, i, k, mode = (((uint32_t)-1)<<2); + uint32_t beg, sink, n, *a, n_occ; + uint64_t pathLen; + bub->ug = ug; + bub->b_bub = bub->b_end_bub = bub->tangle_bub = bub->cross_bub = bub->mess_bub = 0; + if(bub->round_id == 0) + { + buf_t b; memset(&b, 0, sizeof(buf_t)); b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); + uint64_t tLen = get_bub_pop_max_dist_advance(ug->g, &b); + kv_init(bub->list); kv_init(bub->num); kv_init(bub->pathLen); + kv_init(bub->b_s_idx); kv_malloc(bub->b_s_idx, ug->g->n_seq); + bub->b_ug = NULL; kv_init(bub->chain_weight); + bub->b_s_idx.n = ug->g->n_seq; + memset(bub->b_s_idx.a, -1, bub->b_s_idx.n * sizeof(uint64_t)); + CALLOC(bub->index, n_vtx); + for (i = 0; i < ug->g->n_seq; i++) ug->g->seq[i].c = 0; + 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, NULL, 0, 0, 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; + } + } + + kvec_t_u32_warp stack, result; + kv_init(stack.a); kv_init(result.a); + for (v = 0; v < n_vtx; ++v) + { + if((bub->index[v]&(uint32_t)3) !=2) continue; + if(asg_bub_pop1_primary_trio(ug->g, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL, 0, 0, NULL)) + { + //note b.b include end, does not include beg + i = b.b.n + 1; + if(b.b.n == 2 || b.b.n == 3 || b.b.n == 5) + { + for (i = 0; i < b.b.n; i++) + { + if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue; + dfs_bubble(ug->g, &stack, &result, b.b.a[i]>>1, v>>1, b.S.a[0]>>1); + if((result.a.n + 3) != b.b.n && (result.a.n + 2) != b.b.n) break; + } + } + + if(i == b.b.n) + { + kv_push(uint32_t, bub->num, v); + } + else + { + kv_push(uint32_t, bub->num, v + (1<<31)); + } + } + } + kv_destroy(stack.a); kv_destroy(result.a); + radix_sort_u32(bub->num.a, bub->num.a + bub->num.n); + bub->s_bub = 0; + for (k = 0; k < bub->num.n; k++) + { + 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, ug, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0, &pathLen, NULL, NULL, 0, 0, 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->s_bub = bub->num.n - 1; + + + uint64_t dip_thre_max; + memset(r_het_flag, 0, sizeof((*r_het_flag))*sg->n_seq); + if(asm_opt.hom_global_coverage_set) { + dip_thre_max = asm_opt.hom_global_coverage; + } else { + dip_thre_max = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE); + } + // dip_thre_max = (double)(dip_thre_max) - (((double)(dip_thre_max)*0.5)/asm_opt.polyploidy); + dip_thre_max = (double)(dip_thre_max) - ((double)(dip_thre_max)/asm_opt.polyploidy); + for (i = 0; i < ug->g->n_seq; i++) { + bub->index[i] = get_unitig_het_fly(ug, i, sg, dip_thre_max, sources, ruIndex, + r_het_flag, 20, M_het(*bub), P_het(*bub), (uint32_t)-1); + } + + for (i = 0; i < bub->f_bub; i++) + { + get_bubbles(bub, i, &beg, &sink, &a, &n, &pathLen); + for (v = n_occ = 0; v < n; v++) + { + bub->index[(a[v]>>1)] = i; + n_occ += ug->u.a[a[v]>>1].n; + } + + // if((pathLen*2) >= ug->g->seq[beg>>1].len && (pathLen*2) >= ug->g->seq[sink>>1].len) + // { + // bub->index[(beg>>1)] = (uint32_t)-1; + // bub->index[(sink>>1)] = (uint32_t)-1; + // } + + if(n_occ > 3) + { + if(bub->index[(beg>>1)] != M_het(*bub)) bub->index[(beg>>1)] = (uint32_t)-1; + if(bub->index[(sink>>1)] != M_het(*bub)) bub->index[(sink>>1)] = (uint32_t)-1; + } + + + v = beg>>1; + if(bub->b_s_idx.a[v] == (uint64_t)-1) + { + bub->b_s_idx.a[v] <<= 32; + bub->b_s_idx.a[v] |= i; + } + else if((bub->b_s_idx.a[v] & 0xffffffff00000000) == 0xffffffff00000000) + { + bub->b_s_idx.a[v] <<= 32; + bub->b_s_idx.a[v] |= i; + } + + + v = sink>>1; + if(bub->b_s_idx.a[v] == (uint64_t)-1) + { + bub->b_s_idx.a[v] <<= 32; + bub->b_s_idx.a[v] |= i; + } + else if((bub->b_s_idx.a[v] & 0xffffffff00000000) == 0xffffffff00000000) + { + bub->b_s_idx.a[v] <<= 32; + bub->b_s_idx.a[v] |= i; + } + } + for (i = 0; i < ug->g->n_seq; i++) + { + if(bub->index[i] == M_het(*bub)) bub->index[i] = P_het(*bub); + } + } + else + { + bub->num.n = bub->f_bub + 1; + bub->pathLen.n = bub->f_bub; + bub->list.n = bub->num.a[bub->num.n-1]; + update_bub_b_s_idx(bub); + bub->check_het = 0; + asg_destroy(bub->b_g); bub->b_g = NULL; + ma_ug_destroy(bub->b_ug); bub->b_ug = NULL; + kv_destroy(bub->chain_weight); kv_init(bub->chain_weight); + } + bub->b_g = NULL; + bub->b_ug = NULL; + build_bub_graph(ug, bub); + // fprintf(stderr, "-bub->index[18759]: %u, bub->num.n: %u\n", (uint32_t)bub->index[18759], bub->num.n); +} + void print_bubbles(ma_ug_t* ug, bubble_type* bub, kvec_pe_hit* hits, hc_links* link, ha_ug_index* idx) { uint64_t tLen, t_utg, i, k; diff --git a/hic.h b/hic.h index c0f8c82..0f2ffda 100644 --- a/hic.h +++ b/hic.h @@ -90,6 +90,8 @@ int load_hc_links(hc_links* link, const char *fn); void write_hc_links(hc_links* link, const char *fn); void destory_bubbles(bubble_type* bub); void identify_bubbles(ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, kv_u_trans_t *ref); +void identify_bubbles_recal(asg_t* sg, ma_ug_t* ug, bubble_type* bub, uint8_t *r_het_flag, ma_hit_t_alloc* sources, R_to_U* ruIndex, +kv_u_trans_t *ref); void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub); uint32_t connect_bub_occ(bubble_type* bub, uint32_t root_id, uint32_t check_het); void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, uint32_t check_het); diff --git a/inter.cpp b/inter.cpp index 746c52e..baa7638 100644 --- a/inter.cpp +++ b/inter.cpp @@ -4839,8 +4839,8 @@ void update_ul_vec_t_ug(const ul_idx_t *uref, ul_vec_t *rch, vec_mg_lchain_t *uc z->qs = a[m].qs; z->qe = a[m].qe; z->te = a[m].re; z->ts = a[m].rs; z->pidx = k + 1 + m; z->pdis = z->aidx = (uint32_t)-1; - // fprintf(stderr, "[M::%s::k->%ld] m->%ld, utg%.6dl(%c)\n", - // __func__, k, m, (int32_t)z->hid+1, "+-"[z->rev]); + // fprintf(stderr, "[M::%s::k->%ld] m->%ld, utg%.6dl(%c), q::[%u, %u), t::[%u, %u)\n", + // __func__, k, m, (int32_t)z->hid+1, "+-"[z->rev], z->qs, z->qe, z->ts, z->te); } } @@ -5058,7 +5058,7 @@ int64_t debug_i, int64_t tid, void *km) } // #define aln_sc(a, w) (((int64_t)((a).sec))-((int64_t)(((a).qe-(a).qs-(a).sec)*(w)))) -#define aln_sc(a, w) (((int64_t)((a).align_length))-((int64_t)(((a).overlapLen-(a).align_length)*(w)))-(((int64_t)((a).non_homopolymer_errors))*UG_TRANS_ERR_W)) +#define aln_sc(a, trans_w, err_w) (((int64_t)((a).align_length))-((int64_t)(((a).overlapLen-(a).align_length)*(trans_w)))-(((int64_t)((a).non_homopolymer_errors))*(err_w))) void gen_gl_aln(overlap_region_alloc* olist, const ul_idx_t *uref, kv_ul_ov_t *res) { @@ -5090,7 +5090,7 @@ void gen_gg_aln(overlap_region_alloc* olist, const ul_idx_t *uref, int64_t trans for (k = 0; k < olist->length; k++) { p = &(res->a[res->n++]); memset(p, 0, sizeof((*p))); p->v = ((olist->list[k].y_id<<1)|(olist->list[k].y_pos_strand)); p->off = k; - p->score = aln_sc((olist->list[k]), (trans_sc)); + p->score = aln_sc((olist->list[k]), (trans_sc), UG_TRANS_ERR_W); p->qs = olist->list[k].x_pos_s; p->qe = olist->list[k].x_pos_e+1; if((p->v&1)) { p->rs = uref->ug->u.a[p->v>>1].len - (olist->list[k].y_pos_e+1); @@ -5177,16 +5177,69 @@ int64_t get_overlap_region_sub_err(overlap_region *o, rtrace_iter *it, int64_t q return it->werr; } -int64_t cal_gl_chain_lin_sc(ul_ov_t *li, ul_ov_t *lj, rtrace_iter *tc, overlap_region *ol, All_reads *ridx, ma_ug_t *ug, -const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, uint64_t mode, int64_t trans_sc, int64_t sec_sec) + +int64_t get_ecov_el(const ul_idx_t *uref, const ug_opt_t *uopt, uint32_t v, uint32_t w, int64_t bw, double diff_ec_ul, int64_t dq, uint64_t mode, int64_t *el) { + int64_t dt = -1, dif, mm; (*el) = 0; + uint32_t nv, i; asg_arc_t *av = NULL; ///ma_hit_t *x = NULL; + if(!mode) { + const asg_t *g = uref?uref->ug->g:NULL; + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + for (i = 0; i < nv; i++) { + if(av[i].del || av[i].v != w) continue; + dt = av[i].ol; (*el) = av[i].el; + // if(v==1772 && w==1769) fprintf(stderr, "+++v:%u, w:%u, ou:%u\n", v, w, av[i].ou); + // if((v>>1) == 3012 && (w>>1) == 3011) fprintf(stderr, "******************\n"); + break; + } + }else { + ma_hit_t_alloc* src = uopt->sources; + int64_t min_ovlp = uopt->min_ovlp; + int64_t max_hang = uopt->max_hang; + uint64_t z, qn, tn, x = v>>1; int32_t r = 1; asg_arc_t e; + for (z = 0; z < src[x].length; z++) { + qn = Get_qn(src[x].buffer[z]); tn = Get_tn(src[x].buffer[z]); + if(tn != (w>>1)) continue; + r = ma_hit2arc(&(src[x].buffer[z]), Get_READ_LENGTH(R_INF, qn), Get_READ_LENGTH(R_INF, tn), max_hang, asm_opt.max_hang_rate, min_ovlp, &e); + if(r < 0) continue; + if((e.ul>>32) != v || e.v != w) continue; + dt = e.ol; (*el) = src[x].buffer[z].el; + break; + } + } + if(dt < 0) return 0; + dif = (dq>dt? dq-dt:dt-dq); + mm = MAX(dq, dt); mm *= diff_ec_ul; if(mm < bw) mm = bw; + // if((v>>1) == 1163 && (w>>1) == 1168) fprintf(stderr, ">>>>>>dis_q:%ld, dis_t:%ld, dif:%ld, mm:%ld\n", dis_q, dis_t, dif, mm); + if(dif <= mm) return 1; + return 0; +} + +//ai > aj +int64_t cal_gl_chain_lin_sc(ul_ov_t *a, int32_t ai, int32_t aj, rtrace_iter *tc, overlap_region *ol, All_reads *ridx, ma_ug_t *ug, +const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, uint64_t mode, int64_t trans_sc, int64_t sec_sec, +int32_t *f) +{ + ul_ov_t *li = &(a[ai]), *lj = &(a[aj]); ///li is the suffix of lj if(lj->qs >= li->qs) return INT32_MIN; uint32_t li_v = (li->tn<<1)|li->rev, lj_v = (lj->tn<<1)|lj->rev; - int64_t qo = infer_rovlp(li, lj, NULL, NULL, ridx, ug), trans_l = 0, sec_err = 0, sc; ///overlap length in query (UL read) - if(/**li_v != lj_v &&**/ get_ecov_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, mode, NULL)) { + int64_t qo = infer_rovlp(li, lj, NULL, NULL, ridx, ug), trans_l = 0, sec_err = 0, sc, sc0, el; ///overlap length in query (UL read) + ///li_v == lj_v is possiable + if(get_ecov_el(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, mode, &el)) { trans_l = get_overlap_region_sub_err(&(ol[li->qn]), tc, lj->qe, &sec_err); - + sc = f[aj] + ((int64_t)(li->qe - lj->qe)) - (trans_l*trans_sc) - (sec_err*sec_sec); + // fprintf(stderr, "+[M::utg%.6dl] utg%.6dl, liq::[%u, %u), ljq::[%u, %u), trans_l::%ld, sec_err::%ld, sc::%ld, f[aj]::%d\n", + // (int32_t)li->tn+1, (int32_t)lj->tn+1, li->qs, li->qe, lj->qs, lj->qe, trans_l, sec_err, sc, f[aj]); + if((el == 0) && ((trans_l > 0) || (sec_err > 0)) && (lj->qe > li->qs)) { + rtrace_iter tr; tr.k = INT32_MAX; sc0 = sc; + sc = f[aj] + aln_sc(ol[(*li).qn], trans_sc, sec_sec); + trans_l = get_overlap_region_sub_err(&(ol[lj->qn]), &tr, li->qs, &sec_err); + sc -= (((int64_t)(lj->qe - li->qs)) - (trans_l*trans_sc) - (sec_err*sec_sec)); + if(sc < sc0) sc = sc0; + // fprintf(stderr, "-[M::utg%.6dl] utg%.6dl, liq::[%u, %u), ljq::[%u, %u), trans_l::%ld, sec_err::%ld, sc::%ld, f[aj]::%d, aln_sc::%ld\n", + // (int32_t)li->tn+1, (int32_t)lj->tn+1, li->qs, li->qe, lj->qs, lj->qe, trans_l, sec_err, sc, f[aj], aln_sc(ol[(*li).qn], trans_sc, sec_sec)); + } // int64_t trans_l_debug, sec_err_debug; // trans_l_debug = get_overlap_region_sub_err_debug(&(ol[li->qn]), lj->qe, &sec_err_debug); // if(!(trans_l_debug == trans_l && sec_err_debug == sec_err)) { @@ -5194,8 +5247,6 @@ const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, uint6 // __func__, lj->qe, trans_l, trans_l_debug, sec_err, sec_err_debug); // } // assert(trans_l_debug == trans_l && sec_err_debug == sec_err); - - sc = (li->qe - lj->qe) - (trans_l*trans_sc) - (sec_err*sec_sec); // if(li->tn == 308 || li->tn == 311 || lj->tn == 305 || lj->tn == 304) { // fprintf(stderr, "[M::%s::utg%.6dl] utg%.6dl, liq::[%u, %u), ljq::[%u, %u), trans_l::%ld, sec_err::%ld, sc::%ld\n", __func__, // (int32_t)li->tn+1, (int32_t)lj->tn+1, li->qs, li->qe, lj->qs, lj->qe, trans_l, sec_err, sc); @@ -5235,16 +5286,19 @@ uint64_t mode, All_reads *ridx, ma_ug_t *ug, int64_t need_srt) x += li->qs + mm_ovlp; if (x > qlen+1) x = qlen+1; x = find_ul_ov_max(i, res->a, x+G_CHAIN_INDEL); - csc = aln_sc(ol[(*li).qn], trans_sc); + csc = aln_sc(ol[(*li).qn], trans_sc, UG_TRANS_ERR_W); + // fprintf(stderr, "[M::%s::utg%.6dl] i::%ld, csc::%ld, q::[%u, %u), aln::%u, ol::%u, sec_e::%u\n", + // __func__, (int32_t)li->tn+1, i, csc, li->qs, li->qe, + // (ol[(*li).qn]).align_length, (ol[(*li).qn]).overlapLen, + // (ol[(*li).qn]).non_homopolymer_errors); mm_sc = csc; mm_idx = -1; n_skip = 0; end_j = -1; tc.k = INT32_MAX; if ((x-st) > max_iter) st = x-max_iter; for (j = x; j >= st; --j) { // collect potential destination vertices lj = &(res->a[j]); if(lj->qe+G_CHAIN_INDEL <= li->qs) break;//even this pair has a overlap, its length will be very small; just ignore - sc = cal_gl_chain_lin_sc(li, lj, &tc, ol, ridx, ug, uref, uopt, bw, diff_ec_ul, mode, trans_sc, UG_TRANS_ERR_W); + sc = cal_gl_chain_lin_sc(res->a, i, j, &tc, ol, ridx, ug, uref, uopt, bw, diff_ec_ul, mode, trans_sc, UG_TRANS_ERR_W, f); if(sc == INT32_MIN) continue; - sc += f[j]; if(sc > mm_sc) { mm_sc = sc, mm_idx = j; if (n_skip > 0) --n_skip; @@ -5268,9 +5322,8 @@ uint64_t mode, All_reads *ridx, ma_ug_t *ug, int64_t need_srt) if (max_ii >= 0 && max_ii < end_j) {///just have a try with a[i]<->a[max_ii] lj = &(res->a[max_ii]); if(lj->qe+G_CHAIN_INDEL > li->qs && lj->qs < li->qs) { - sc = cal_gl_chain_lin_sc(li, lj, &tc, ol, ridx, ug, uref, uopt, bw, diff_ec_ul, mode, trans_sc, UG_TRANS_ERR_W); + sc = cal_gl_chain_lin_sc(res->a, i, max_ii, &tc, ol, ridx, ug, uref, uopt, bw, diff_ec_ul, mode, trans_sc, UG_TRANS_ERR_W, f); if(sc != INT32_MIN) { - sc += f[max_ii]; if(sc > mm_sc) { mm_sc = sc; mm_idx = max_ii; } @@ -5747,11 +5800,12 @@ st_mt_t *bf, Chain_Data* dp, int64_t max_skip, int64_t max_iter, int64_t max_dis -inline int32_t cal_gchain_sc_adv(overlap_region *ol, const mg_path_dst_t *dj, const mg_lchain_t *li, const mg_lchain_t *lc, int64_t *f, int64_t b_w, float diff_thre, float chn_pen_gap, -rtrace_iter *tc, int64_t trans_sc, int64_t sec_sec) +inline int32_t cal_gchain_sc_adv(const ma_ug_t *ug, const ul_idx_t *uref, const ug_opt_t *uopt, +overlap_region *ol, const mg_path_dst_t *dj, const mg_lchain_t *li, ul_ov_t *ui, mg_lchain_t *lc, int64_t *f, +int64_t b_w, float diff_thre, float chn_pen_gap, rtrace_iter *tc, int64_t trans_sc, int64_t sec_sec) { // const mg_lchain_t *lj; - int32_t gap, sc; + int32_t gap; float lin_pen, log_pen; if (dj->n_path == 0) return INT32_MIN; gap = dj->dist - dj->target_dist; @@ -5760,8 +5814,23 @@ rtrace_iter *tc, int64_t trans_sc, int64_t sec_sec) if ((gap > ((dj->target_dist)*diff_thre)) && (gap > b_w)) return INT32_MIN; // if (lj->qe <= li->qs) sc = li->score; // else sc = (int32_t)((double)(li->qe - lj->qe) / (li->qe - li->qs) * li->score + .499); // dealing with overlap on query - int64_t trans_l = 0, sec_err = 0; - trans_l = get_overlap_region_sub_err(&(ol[li->off]), tc, lc[dj->meta].qe, &sec_err); + int64_t trans_l = 0, sec_err = 0, qo, el = 0, sc, sc0; + mg_lchain_t *lj = &(lc[dj->meta]); ul_ov_t uj; + + trans_l = get_overlap_region_sub_err(&(ol[li->off]), tc, lj->qe, &sec_err); + sc = (li->qe - lj->qe) - (trans_l*trans_sc) - (sec_err*sec_sec); sc += f[dj->meta]; + if(((trans_l > 0) || (sec_err > 0)) && (lj->qe > li->qs)) { + set_ul_ov_t_by_mg_lchain_t(&uj, lj); + qo = infer_rovlp(ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug); + if((!get_ecov_el(uref, uopt, li->v^1, lj->v^1, b_w, N_GCHAIN_RATE, qo, 0, &el)) || (el == 0)) { + rtrace_iter tr; tr.k = INT32_MAX; sc0 = sc; + sc = f[dj->meta] + li->score; + trans_l = get_overlap_region_sub_err(&(ol[lj->off]), &tr, li->qs, &sec_err); + sc -= (((int64_t)(lj->qe - li->qs)) - (trans_l*trans_sc) - (sec_err*sec_sec)); + if(sc < sc0) sc = sc0; + } + } + // int64_t trans_l_debug, sec_err_debug; // trans_l_debug = get_overlap_region_sub_err_debug(&(ol[li->off]), lc[dj->meta].qe, &sec_err_debug); @@ -5771,14 +5840,13 @@ rtrace_iter *tc, int64_t trans_sc, int64_t sec_sec) // } // assert(trans_l_debug == trans_l && sec_err_debug == sec_err); - sc = (li->qe - lc[dj->meta].qe) - (trans_l*trans_sc) - (sec_err*sec_sec); + // sc = li->score; //sc += dj->mlen; // TODO: is this line the right thing to do? // if (dj->is_0) sc += ref_bonus; lin_pen = chn_pen_gap * (float)gap; log_pen = gap >= 2? mg_log2(gap) : 0.0f; sc -= (int32_t)(lin_pen + log_pen); - sc += f[dj->meta]; return sc; } @@ -5900,7 +5968,7 @@ int64_t need_srt) for (j = 0; j < dst_n; ++j) { dj = &dst->a[j]; if (dj->n_path == 0) continue; // unreachable - sc = cal_gchain_sc_adv(ol, dj, li, lc->a, f, bw, diff_thre, W_CHN_PEN_GAP, &tc, trans_sc, sec_sec); + sc = cal_gchain_sc_adv(ug, uref, uopt, ol, dj, li, &ui, lc->a, f, bw, diff_thre, W_CHN_PEN_GAP, &tc, trans_sc, sec_sec); if (sc == INT32_MIN) continue; // out of band // if (sc < 0) continue;// negative score @@ -8811,11 +8879,19 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call // void *km = s->buf?(s->buf[tid]?s->buf[tid]->km:NULL):NULL; // if(s->id+i!=3046/** && s->id+i!=3111**/) return; // if((s->id+i!=871) && (s->id+i!=963) && (s->id+i!=980)) return; - // if(s->id+i!=963) return; + // if(s->id+i!=944) return; // if(s->id+i != 35437) return; + // if((s->id+i != 7086) && (s->id+i != 51705) && (s->id+i != 266022) && (s->id+i != 353608) + // && (s->id+i != 399416) && (s->id+i != 403014) && (s->id+i != 420915) && (s->id+i != 603855) + // && (s->id+i != 680134) && (s->id+i != 766261) && (s->id+i != 794527)) { + // return; + // } // fprintf(stderr, "\n[M::%s] rid::%ld, len::%lu, name::%.*s\n", __func__, s->id+i, s->len[i], // (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a); + // fprintf(stderr, ">%.*s\n%.*s\n", (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a, + // (int32_t)s->len[i], s->seq[i]); + // if (memcmp(UL_INF.nid.a[s->id+i].a, "d0aab024-b3a7-40fb-83cc-22c3d6d951f8", UL_INF.nid.a[s->id+i].n-1)) return; // fprintf(stderr, "[M::%s::] ==> len: %lu\n", __func__, s->len[i]); // ha_get_ul_candidates_interface(b->abl, i, s->seq[i], s->len[i], s->opt->w, s->opt->k, s->uu, &b->olist, &b->olist_hp, &b->clist, s->opt->bw_thres, @@ -14690,7 +14766,7 @@ ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_c ///for debug interval if(!load_all_ul_t(&UL_INF, gfa_name, &R_INF, ug)) { gen_UL_reovlps(&sl, ug, sg, gfa_name, cutoff); - exit(1); + // exit(1); write_all_ul_t(&UL_INF, gfa_name, ug); } else if(double_check_cache){ if(drenew_UL_reovlps(&sl, ug, sg, gfa_name, cutoff)) { diff --git a/inter.h b/inter.h index 1684581..c51f8ba 100644 --- a/inter.h +++ b/inter.h @@ -17,8 +17,10 @@ #define UG_SKIP_N 100 #define UG_ITER_N 5000 #define UG_DIS_N 50000 +// #define UG_TRANS_W 2 #define UG_TRANS_W 2 -#define UG_TRANS_ERR_W 512 +// #define UG_TRANS_ERR_W 512 +#define UG_TRANS_ERR_W 64 #define G_CHAIN_TRANS_RATE 0.25 #define G_CHAIN_TRANS_WEIGHT -1 #define G_CHAIN_INDEL 128