diff --git a/CommandLines.cpp b/CommandLines.cpp index 1fdeebb..186b952 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -294,6 +294,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->ul_pst_join = 1; asm_opt->trio_cov_het_ovlp = -1; asm_opt->ul_min_base = 0; + asm_opt->self_scaf = 0; } void destory_enzyme(enzyme* f) diff --git a/CommandLines.h b/CommandLines.h index b037651..1b4eb33 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.19.6-r595" +#define HA_VERSION "0.19.6-r597" #define VERBOSE 0 @@ -147,6 +147,7 @@ typedef struct { uint8_t dbg_ovec_cal; uint8_t hifi_pst_join, ul_pst_join; uint32_t ul_min_base; + uint8_t self_scaf; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Overlaps.cpp b/Overlaps.cpp index a43877c..9a455e7 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -75,7 +75,7 @@ void print_vw_edge(asg_t *sg, uint32_t vid, uint32_t wid, const char *cmd); void output_trio_graph_joint(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, -int min_ovlp, long long gap_fuzz, bub_label_t* b_mask_t, ma_ug_t **rhu0, ma_ug_t **rhu1); +int min_ovlp, long long gap_fuzz, bub_label_t* b_mask_t, ma_ug_t **rhu0, ma_ug_t **rhu1, ug_opt_t *opt); typedef struct { uint32_t d, tot, ma, p; @@ -143,6 +143,16 @@ typedef struct { uint64_t ridx_n, ra_n; } dedup_idx_t; +typedef struct { + uint32_t id0, id1; + uint64_t len; +} NN_t; + +typedef struct { + NN_t *a; + size_t n, m; +} kvect_N_t; + ///this value has been updated at the first line of build_string_graph_without_clean long long min_thres; @@ -9353,6 +9363,171 @@ ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, k } +void reduce_ma_utg_t_scaf(ma_utg_t *in, asg_t *rg, ma_hit_t_alloc* src, ma_sub_t *cov, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* newE) +{ + asg_arc_t t_f; + asg_arc_t* av = NULL; + uint32_t m = 0, k, i, nv; + for (i = 0; i < in->n; i++) { + if(in->a[i] == (uint64_t)-1) { + // if(newE) asg_seq_del(read_g, collection->a[i]>>33); + continue; + } + + in->a[m] = in->a[i]; + m++; + } + in->n = m; + + + uint32_t totalLen = 0, v, w, l; + for (i = 0; i + 1 < in->n; i++) { + v = (uint64_t)(in->a[i])>>32; + w = (uint64_t)(in->a[i + 1])>>32; + l = Get_READ_LENGTH(R_INF, (v>>1)); ///for Ns + + if((!IS_SCAF_READ(R_INF, v>>1)) && (!IS_SCAF_READ(R_INF, w>>1))) { + /*******************************for debug************************************/ + l = (uint32_t)-1; + av = asg_arc_a(rg, v); nv = asg_arc_n(rg, v); + for (k = 0; k < nv; k++) { + if(av[k].del) continue; + if(av[k].v == w) { + l = asg_arc_len(av[k]); + break; + } + } + + if(k == nv) { + for (k = 0; k < edge->a.n; k++) { + if(edge->a.a[k].del) continue; + if((edge->a.a[k].ul>>32) == v && edge->a.a[k].v == w) { + l = asg_arc_len(edge->a.a[k]); + break; + } + } + + if(k == edge->a.n) { + if(get_edge_from_source(src, cov, NULL, max_hang, min_ovlp, v, w, &t_f)==0) { + fprintf(stderr, "####ERROR1: v>>1: %u, v&1: %u, w>>1: %u, w&1: %u, r_seq: %u\n", + v>>1, v&1, w>>1, w&1, rg->r_seq); + } + l = asg_arc_len(t_f); + + + if(newE) { + kv_push(asg_arc_t, newE->a, t_f); + if(get_edge_from_source(src, cov, NULL, max_hang, min_ovlp, w^1, v^1, &t_f)==0) { + fprintf(stderr, "####ERROR2: v>>1: %u, v&1: %u, w>>1: %u, w&1: %u, r_seq: %u\n", + v>>1, v&1, w>>1, w&1, rg->r_seq); + } + kv_push(asg_arc_t, newE->a, t_f); + } + + } + else if(newE) { + kv_push(asg_arc_t, newE->a, edge->a.a[k]); + for (k = 0; k < edge->a.n; k++) { + if(edge->a.a[k].del) continue; + if((edge->a.a[k].ul>>32) == (w^1) && edge->a.a[k].v == (v^1)) { + l = asg_arc_len(edge->a.a[k]); + break; + } + } + kv_push(asg_arc_t, newE->a, edge->a.a[k]); + } + } + if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n"); + /*******************************for debug************************************/ + } + + + in->a[i] = v; in->a[i] = in->a[i]<<32; + in->a[i] = in->a[i] | (uint64_t)(l); + totalLen += l; + } + + if(i < in->n) { + if(in->circ) { + v = (uint64_t)(in->a[i])>>32; + w = (uint64_t)(in->a[0])>>32; + l = Get_READ_LENGTH(R_INF, (v>>1)); ///for Ns + + if((!IS_SCAF_READ(R_INF, v>>1)) && (!IS_SCAF_READ(R_INF, w>>1))) { + /*******************************for debug************************************/ + l = (uint32_t)-1; + av = asg_arc_a(rg, v); + nv = asg_arc_n(rg, v); + for (k = 0; k < nv; k++) { + if(av[k].del) continue; + if(av[k].v == w) { + l = asg_arc_len(av[k]); + break; + } + } + + if(k == nv) { + for (k = 0; k < edge->a.n; k++) { + if(edge->a.a[k].del) continue; + if((edge->a.a[k].ul>>32) == v && edge->a.a[k].v == w) { + l = asg_arc_len(edge->a.a[k]); + break; + } + } + + if(k == edge->a.n) { + if(get_edge_from_source(src, cov, NULL, max_hang, min_ovlp, v, w, &t_f)==0) { + fprintf(stderr, "####ERROR1: v>>1: %u, v&1: %u, w>>1: %u, w&1: %u, r_seq: %u\n", + v>>1, v&1, w>>1, w&1, rg->r_seq); + } + l = asg_arc_len(t_f); + + if(newE) { + kv_push(asg_arc_t, newE->a, t_f); + if(get_edge_from_source(src, cov, NULL, max_hang, min_ovlp, w^1, v^1, &t_f)==0) { + fprintf(stderr, "####ERROR2: v>>1: %u, v&1: %u, w>>1: %u, w&1: %u, r_seq: %u\n", + v>>1, v&1, w>>1, w&1, rg->r_seq); + } + kv_push(asg_arc_t, newE->a, t_f); + } + } else if(newE) { + kv_push(asg_arc_t, newE->a, edge->a.a[k]); + for (k = 0; k < edge->a.n; k++) { + if(edge->a.a[k].del) continue; + if((edge->a.a[k].ul>>32) == (w^1) && edge->a.a[k].v == (v^1)) { + l = asg_arc_len(edge->a.a[k]); + break; + } + } + kv_push(asg_arc_t, newE->a, edge->a.a[k]); + } + } + if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n"); + /*******************************for debug************************************/ + } + + in->a[i] = v; in->a[i] = in->a[i]<<32; + in->a[i] = in->a[i] | (uint64_t)(l); + + totalLen += l; + } else { + v = (uint64_t)(in->a[i])>>32; + l = rg->seq[v>>1].len; + in->a[i] = v; + in->a[i] = in->a[i]<<32; + in->a[i] = in->a[i] | (uint64_t)(l); + totalLen += l; + } + } + + in->len = totalLen; + if(!in->circ) { + in->start = in->a[0]>>32; + in->end = (in->a[in->n-1]>>32)^1; + } +} + + uint32_t detect_exact_ovec(ma_utg_t* collection, asg_t* read_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, uint32_t src, uint32_t dest_idx) @@ -9442,11 +9617,9 @@ kvec_asg_arc_t_warp* newE) uint32_t i, k, v, pre, pre_i, afte, afte_i, exactLen, inexactLen, skip = 0; uint32_t min_inexactLen, max_exactLen; asg_arc_t pE, aE; + pre = (uint64_t)(collection->a[0])>>32; pre_i = 0; afte_i = (uint32_t)-1; - - - for (i = 1; i < collection->n - 1; i++) - { + for (i = 1; i < collection->n - 1; i++) { if(collection->a[i] == (uint64_t)-1) continue; ///v and after must be available v = (uint64_t)(collection->a[i])>>32; @@ -9514,6 +9687,80 @@ kvec_asg_arc_t_warp* newE) } +uint32_t polish_unitig_scaf(ma_utg_t* in, asg_t* rg, ma_hit_t_alloc* src, ma_sub_t *cov, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* newE) +{ + if(in->m == 0) return 0; + if(in->n < 3) return 0; + uint32_t i, k, v, pre, pre_i, afte, afte_i, exactLen, inexactLen, skip = 0, z, l; + uint32_t min_inexactLen, max_exactLen; + asg_arc_t pE, aE; + + for (z = 1, l = 0; z <= in->n; z++) { + if(z == in->n || IS_SCAF_READ(R_INF, (in->a[z])>>33) ) { + if(z - l >= 3) { + pre = (uint64_t)(in->a[l])>>32; pre_i = l; afte_i = (uint32_t)-1; + for (i = l + 1; i < z - 1; i++) { + if(in->a[i] == (uint64_t)-1) continue; + ///v and after must be available + v = (uint64_t)(in->a[i])>>32; + afte = (uint64_t)(in->a[i+1])>>32; + + get_specific_edge(src, cov, NULL, edge, pre_i == i-1? rg:NULL, max_hang, min_ovlp, v^1, pre^1, &pE); + get_specific_edge(src, cov, NULL, edge, rg, max_hang, min_ovlp, v, afte, &aE); + + if(pE.el == 1 && aE.el == 1) { + pre = (uint64_t)(in->a[i])>>32; pre_i = i; + continue; + } + + ///pre must be a good read, we need to find a good after + ///update a new afte + afte_i = detect_exact_ovec(in, rg, src, cov, edge, max_hang, min_ovlp, pre, i+1); + if(afte_i == (uint32_t)-1) { + pre = (uint64_t)(in->a[i])>>32; pre_i = i; + continue; + } + afte = (uint64_t)(in->a[afte_i])>>32; + + + min_inexactLen = (uint32_t)-1;max_exactLen = 0; + for (k = i; k < afte_i; k++) { + get_overlapLen((uint64_t)(in->a[k])>>33, src, &exactLen, &inexactLen); + if(inexactLen < min_inexactLen) { + min_inexactLen = inexactLen; + max_exactLen = exactLen; + } + } + + get_overlapLen(pre>>1, src, &exactLen, &inexactLen); + if((inexactLen > min_inexactLen) || (inexactLen == min_inexactLen && exactLen <= max_exactLen)) { + pre = (uint64_t)(in->a[i])>>32; pre_i = i; + continue; + } + + get_overlapLen(afte>>1, src, &exactLen, &inexactLen); + if((inexactLen > min_inexactLen) || (inexactLen == min_inexactLen && exactLen <= max_exactLen)) { + pre = (uint64_t)(in->a[i])>>32; pre_i = i; + continue; + } + + for (k = i; k < afte_i; k++) { + in->a[k] = (uint64_t)-1; + skip++; + } + } + } + l = z; + } + } + + if(skip == 0) return 0; + reduce_ma_utg_t_scaf(in, rg, src, cov, edge, max_hang, min_ovlp, newE); + + return 1; +} + + void print_rough_inconsistent_sites(ma_utg_t* collection, uint32_t cur_i, uint32_t next_i, asg_t* read_g, All_reads *RNF, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, UC_Read* r_read, UC_Read* q_read, int max_hang, int min_ovlp, @@ -9918,6 +10165,69 @@ UC_Read* r_read, UC_Read* q_read, int max_hang, int min_ovlp, kvec_asg_arc_t_war } +uint32_t polish_unitig_advance_scaf(ma_utg_t* in, asg_t* rg, All_reads *RNF, ma_hit_t_alloc* src, ma_sub_t *cov, kvec_asg_arc_t_warp* edge, +UC_Read* r_read, UC_Read* q_read, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* newE) +{ + if(in->m == 0) return 0; + if(in->n < 3) return 0; + // uint32_t i, skip = 0; + int match_v, total_v, max_i, match_max, k, z, l, in_n = in->n, i, skip = 0; + double match_rate, match_rate_max; + + + for (z = 1, l = 0; z <= in_n; z++) { + if(z == in_n || IS_SCAF_READ(R_INF, (in->a[z])>>33) ) { + if(z - l >= 3) { + ///we should be able to handle i = z-1 + for (i = l + 1; i < z - 1; i++) { + ///in practice, in->a[i] and in->a[i+1] must be available + ///in->a[index] might be unavailable only if index < i + if(get_consensus_rate(in, i, i+1, rg, RNF, src, cov, edge, r_read, q_read, max_hang, min_ovlp, &match_v, &total_v) != 1) { + continue; + } + + match_rate = (total_v == 0)? 0:((double)(match_v)/(double)(total_v)); + ///most reads support in[i], so it is right + if(match_v >= total_v * 0.5 && total_v > 0 && match_v > 0) continue; + + max_i = i; match_max = match_v; match_rate_max = match_rate; + for (k = i - 1; k >= 0; k--) { + if(in->a[k] == (uint64_t)-1) continue; + ///in->a[k] might be unavailable, while in->a[i+1] must be available + ///return -1 means there is no overlap from k to i+1 + if(get_consensus_rate(in, k, i+1, rg, RNF, src, cov, edge, r_read, q_read, max_hang, min_ovlp, &match_v, &total_v) < 0) { + break; + } + + ///no read support k to i+1 + if(total_v == 0) break; + + match_rate = (total_v == 0)? 0:((double)(match_v)/(double)(total_v)); + if(match_rate > match_rate_max || (match_rate == match_rate_max && match_v > match_max)) { + max_i = k; match_max = match_v; match_rate_max = match_rate; + } + } + + ///set [max_i+1, i] to be unavailable + for (k = max_i+1; k <= (int)i; k++) { + if(in->a[k] == (uint64_t)-1) continue; + in->a[k] = (uint64_t)-1; + skip++; + } + } + } + l = z; + } + } + + if(skip == 0) return 1; + + reduce_ma_utg_t_scaf(in, rg, src, cov, edge, max_hang, min_ovlp, newE); + + return 1; +} + + ma_ug_t *gen_polished_ug(const ug_opt_t *uopt, asg_t *sg) { kvec_asg_arc_t_warp e, d; @@ -9984,10 +10294,9 @@ kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, kvec_asg_arc_t_warp *E, u for (i = 0; i < g->u.n; ++i) { ma_utg_t *u = &g->u.a[i]; if(u->m == 0) continue; - if(is_polish) - { - polish_unitig(u, read_g, sources, coverage_cut, edge, max_hang, min_ovlp, E); - polish_unitig_advance(u, read_g, &R_INF, sources, coverage_cut, edge, &g_read, &tmp, max_hang, min_ovlp, E); + if(is_polish) { + polish_unitig_scaf(u, read_g, sources, coverage_cut, edge, max_hang, min_ovlp, E); + polish_unitig_advance_scaf(u, read_g, &R_INF, sources, coverage_cut, edge, &g_read, &tmp, max_hang, min_ovlp, E); } g->g->seq[i].len = u->len; @@ -10004,12 +10313,9 @@ kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, kvec_asg_arc_t_warp *E, u l += eLen; if(eLen == 0) continue; - if(rId < read_g->r_seq) - { + if(rId < read_g->r_seq) { recover_UC_Read(&g_read, &R_INF, rId); - } - else - { + } else { recover_fake_read(&g_read, &tmp, &(read_g->F_seq[rId-read_g->r_seq]), &R_INF, coverage_cut); } @@ -10017,17 +10323,12 @@ kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp, kvec_asg_arc_t_warp *E, u readS = g_read.seq + coverage_cut[rId].s; readLen = coverage_cut[rId].e - coverage_cut[rId].s; - if (!ori) // forward strand - { - for (k = 0; k < eLen; k++) - { + if (!ori) {// forward strand + for (k = 0; k < eLen; k++) { u->s[start + k] = readS[k]; } - } - else - { - for (k = 0; k < eLen; k++) - { + } else { + for (k = 0; k < eLen; k++) { uint8_t c = (uint8_t)readS[readLen - 1 - k]; u->s[start + k] = c >= 128? 'N' : comp_tab[c]; } @@ -15758,7 +16059,7 @@ long long gap_fuzz, bub_label_t* b_mask_t, ug_opt_t *opt) output_trio_graph_joint(sg, coverage_cut, output_file_name, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, - 0.05, 0.9, max_hang, min_ovlp, gap_fuzz, b_mask_t, rhits?(&ug_fa):NULL, rhits?(&ug_mo):NULL); + 0.05, 0.9, max_hang, min_ovlp, gap_fuzz, b_mask_t, rhits?(&ug_fa):NULL, rhits?(&ug_mo):NULL, opt); if(rhits) { ha_aware_order(rhits, sg, ug_fa, ug_mo, cov?&(cov->t_ch->k_trans):&(t_ch->k_trans), opt, 3); @@ -17031,7 +17332,7 @@ long long gap_fuzz, ug_opt_t *opt) // 0.05, 0.9, max_hang, min_ovlp, gap_fuzz, 0, b_mask_t, NULL, NULL, NULL); output_trio_graph_joint(sg, coverage_cut, output_file_name, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, - 0.05, 0.9, max_hang, min_ovlp, gap_fuzz, b_mask_t, NULL, NULL); + 0.05, 0.9, max_hang, min_ovlp, gap_fuzz, b_mask_t, NULL, NULL, opt); } void output_bp_trio_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, @@ -17087,7 +17388,7 @@ long long gap_fuzz, ug_opt_t *opt) output_trio_graph_joint(sg, coverage_cut, output_file_name, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, - 0.05, 0.9, max_hang, min_ovlp, gap_fuzz, b_mask_t, NULL, NULL); + 0.05, 0.9, max_hang, min_ovlp, gap_fuzz, b_mask_t, NULL, NULL, opt); } ma_ug_t* merge_utg(ma_ug_t **dest, ma_ug_t **src) @@ -20835,10 +21136,183 @@ uint64_t append_miss_nid(asg_t *sg, ma_ug_t *hap0, ma_ug_t *hap1, uint8_t *ff, u return n_base; } +void prt_scaf_res_t(scaf_res_t *pa, ma_ug_t *ref, ma_ug_t *ctg) +{ + uint32_t k, i, z, a_n; ma_utg_t *rch; ul_vec_t *idx; uc_block_t *a; + for(i = 0; i < ctg->u.n; i++) { + rch = &(ctg->u.a[i]); idx = &(pa->a[i]); + fprintf(stderr, "[M::%s] rch->len::%u, rch->n::%u, idx->n::%u\n", __func__, (uint32_t)rch->len, (uint32_t)rch->n, (uint32_t)idx->bb.n); + // for (k = 0; k < idx->bb.n; k++) { + // fprintf(stderr, "[k->%u::utg%.6u%c(len->%u::n->%u)]\tq::[%u, %u)\t%c\tt::[%u, %u)\n", + // k, (idx->bb.a[k].hid) + 1, "lc"[ref->u.a[(idx->bb.a[k].hid)].circ], ref->u.a[(idx->bb.a[k].hid)].len, ref->u.a[(idx->bb.a[k].hid)].n, + // idx->bb.a[k].qs, idx->bb.a[k].qe, "+-"[idx->bb.a[k].rev], idx->bb.a[k].ts, idx->bb.a[k].te); + // } + + for (k = 0; k < idx->bb.n; k++) { + a = idx->bb.a + idx->bb.a[k].ts; a_n = idx->bb.a[k].te - idx->bb.a[k].ts; + fprintf(stderr, "q::[%u, %u)\n", idx->bb.a[k].qs, idx->bb.a[k].qe); + for (z = 0; z < a_n; z++) fprintf(stderr, "utg%.6u%c,", (a[z].hid)+1, "lc"[ref->u.a[a[z].hid].circ]); + fprintf(stderr, "\n"); + } + } +} + + +typedef struct { + uint32_t len[2], num[2], h; +} ug_res_t; + +typedef struct { + uint32_t *idx; + ug_res_t *map; + uint8_t *f; + ma_ug_t *ref; + kvec_t_u32_warp st, res; +} tangle_res_t; + + +ma_ug_t *gen_clean_ug(ma_ug_t *ref, scaf_res_t *cp0, scaf_res_t *cp1) +{ + ma_ug_t *ug = copy_untig_graph(ref); + uint32_t i, k, a_n, z, v, w; scaf_res_t *ctg = NULL; ul_vec_t *idx; uc_block_t *a; + if(cp0 || cp1) { + for (i = 0; i < ug->g->n_seq; i++) ug->g->seq[i].del = 1; + for (i = 0; i < ug->g->n_arc; i++) ug->g->arc[i].del = 1; + + ctg = cp0; + if (ctg) { + for(i = 0; i < ctg->n; i++) { + idx = &(ctg->a[i]); + for (k = 0; k < idx->bb.n; k++) { + a = idx->bb.a + idx->bb.a[k].ts; a_n = idx->bb.a[k].te - idx->bb.a[k].ts; + for (z = 0; z < a_n; z++) ug->g->seq[a[z].hid].del = 0; + for (z = 1; z < a_n; z++) { + v = a[z-1].hid<<1; v |= (uint32_t)a[z-1].rev; + w = a[z].hid<<1; w |= (uint32_t)a[z].rev; + + asg_arc_del(ug->g, v, w, 0); + asg_arc_del(ug->g, w^1, v^1, 0); + } + } + } + } + + ctg = cp1; + if (ctg) { + for(i = 0; i < ctg->n; i++) { + idx = &(ctg->a[i]); + for (k = 0; k < idx->bb.n; k++) { + a = idx->bb.a + idx->bb.a[k].ts; a_n = idx->bb.a[k].te - idx->bb.a[k].ts; + for (z = 0; z < a_n; z++) ug->g->seq[a[z].hid].del = 0; + for (z = 1; z < a_n; z++) { + v = a[z-1].hid<<1; v |= (uint32_t)a[z-1].rev; + w = a[z].hid<<1; w |= (uint32_t)a[z].rev; + + asg_arc_del(ug->g, v, w, 0); + asg_arc_del(ug->g, w^1, v^1, 0); + } + } + } + } + } + + + for (i = 0; i < ug->g->n_seq; i++) { + if(ug->g->seq[i].del == 1) { + asg_seq_del(ug->g, i); + if(ug->u.a[i].m!=0) { + ug->u.a[i].m = ug->u.a[i].n = 0; + free(ug->u.a[i].a); + ug->u.a[i].a = NULL; + } + } + } + + asg_cleanup(ug->g); + + return ug; +} + +void double_scaffold(ma_ug_t *ref, ma_ug_t *hu0, ma_ug_t *hu1, scaf_res_t *cp0, scaf_res_t *cp1, asg_t *sg, ma_sub_t* cover, ma_hit_t_alloc* src, ma_hit_t_alloc* rev, +long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, bub_label_t* b_mask_t, ug_opt_t *opt) +{ + kvec_asg_arc_t_warp new_rtg_edges; kv_init(new_rtg_edges.a); + ma_ug_t *cref = gen_clean_ug(ref, cp0, cp1); + asg_t *csg = copy_read_graph(sg); + hap_cov_t *cov = NULL; + + print_debug_gfa(sg, cref, cover, "cl.sb.utg", src, ruIndex, opt->max_hang, opt->min_ovlp, 0, 0, 0); + + new_rtg_edges.a.n = 0; + ///asm_opt.purge_overlap_len = asm_opt.purge_overlap_len_hic; + ///asm_opt.purge_simi_thres = asm_opt.purge_simi_rate_hic; + adjust_utg_by_primary(&cref, csg, TRIO_THRES, src, rev, cover, + tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, + max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1, 0); + + + ma_ug_destroy(cref); + asg_destroy(csg); + // clean_u_trans_t_idx(&(cov->t_ch->k_trans), ug, sg); + + filter_u_trans(&(cov->t_ch->k_trans), asm_opt.is_bub_trans, asm_opt.is_topo_trans, asm_opt.is_read_trans, asm_opt.is_base_trans); + clean_u_trans_t_idx_filter_adv(&(cov->t_ch->k_trans), ref, sg); + // dbg_prt_utg_trans(&(cov->t_ch->k_trans), ref, "after"); + + /**kv_u_trans_t *os =**/ gen_contig_trans(opt, sg, hu1, cp1, hu0, cp0, ref, &(cov->t_ch->k_trans)); + + + +} + + +void gen_self_scaf(ug_opt_t *opt, ma_ug_t *hu0, ma_ug_t *hu1, asg_t *sg, ma_sub_t *cov, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U *ri, +long long tipsLen, float tip_drop_ratio, long long stops_threshold, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, bub_label_t* b_mask_t) +{ + /** + dedup_idx_t *hidx0 = NULL, *hidx1 = NULL, *uidx = NULL; uint8_t *ff = NULL; if(hu0 || hu1) CALLOC(ff, sg->n_seq); + ma_ug_t *ug = NULL; uint64_t pscut = 0; kvect_N_t Ns; kv_init(Ns); + // ug = ma_ug_gen(sg); + pscut = (asm_opt.hom_global_coverage_set?(asm_opt.hom_global_coverage):(((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE))); + pscut *= PHASE_SEF; if(pscut < PHASE_SEP) pscut = PHASE_SEP; + ug = ma_ug_gen_phase(sg, pscut, PHASE_SEP_RATE); + + + uidx = gen_dedup_idx_t(ug, sg); if(hu0) hidx0 = gen_dedup_idx_t(hu0, sg); if(hu1) hidx1 = gen_dedup_idx_t(hu1, sg); + update_recover_atg_cov(); + + if(hidx0 && haploid_self_scaf(hidx0, uidx)); + if(hidx1); + if(hidx0 && hidx1); + + + if(hidx0) destroy_dedup_idx_t(hidx0); if(hidx1) destroy_dedup_idx_t(hidx1); if(uidx) destroy_dedup_idx_t(uidx); + free(ff); ma_ug_destroy(ug); kv_destroy(Ns); + **/ + ma_ug_t *ug = NULL; + //uint64_t pscut = 0; + //pscut = (asm_opt.hom_global_coverage_set?(asm_opt.hom_global_coverage):(((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE))); + //pscut *= PHASE_SEF; if(pscut < PHASE_SEP) pscut = PHASE_SEP; + //ug = ma_ug_gen_phase(sg, pscut, PHASE_SEP_RATE); + ug = ma_ug_gen(sg); + + print_debug_gfa(sg, ug, cov, "sb.utg", src, ri, opt->max_hang, opt->min_ovlp, 0, 0, 0); + + fprintf(stderr, "[M::%s] hu0\n", __func__); + scaf_res_t *cp0 = gen_contig_path(opt, sg, hu0, ug); prt_scaf_res_t(cp0, ug, hu0); + fprintf(stderr, "[M::%s] hu1\n", __func__); + scaf_res_t *cp1 = gen_contig_path(opt, sg, hu1, ug); prt_scaf_res_t(cp1, ug, hu1); + + double_scaffold(ug, hu0, hu1, cp0, cp1, sg, cov, src, rev, tipsLen, tip_drop_ratio, stops_threshold, ri, chimeric_rate, drop_ratio, max_hang, min_ovlp, b_mask_t, opt); + + ma_ug_destroy(ug); destroy_scaf_res_t(cp0); destroy_scaf_res_t(cp1); + +} + void output_trio_graph_joint(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, -int min_ovlp, long long gap_fuzz, bub_label_t* b_mask_t, ma_ug_t **rhu0, ma_ug_t **rhu1) +int min_ovlp, long long gap_fuzz, bub_label_t* b_mask_t, ma_ug_t **rhu0, ma_ug_t **rhu1, ug_opt_t *opt) { ma_ug_t *hu0 = NULL, *hu1 = NULL; kvec_asg_arc_t_warp arcs0, arcs1; memset(&arcs0, 0, sizeof(arcs0)); memset(&arcs1, 0, sizeof(arcs1)); @@ -20869,6 +21343,10 @@ int min_ovlp, long long gap_fuzz, bub_label_t* b_mask_t, ma_ug_t **rhu0, ma_ug_t renew_utg((&hu0), sg, &arcs0); renew_utg((&hu1), sg, &arcs1); fprintf(stderr, "[M::%s] dedup_base::%lu, miss_base::%lu\n", __func__, dedup_base, miss_base); + if(asm_opt.self_scaf) { + gen_self_scaf(opt, hu0, hu1, sg, coverage_cut, sources, reverse_sources, ruIndex, tipsLen, tip_drop_ratio, stops_threshold, chimeric_rate, drop_ratio, max_hang, min_ovlp, b_mask_t); + } + if(!rhu0) { output_hap_graph(hu0, sg, &arcs0, coverage_cut, output_file_name, FATHER, sources, ruIndex, max_hang, min_ovlp, NULL); ma_ug_destroy(hu0); @@ -36454,7 +36932,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, &b_mask_t, gap_fuzz, &uopt); } else { output_trio_graph_joint(sg, coverage_cut, o_file, sources, reverse_sources, (asm_opt.max_short_tip*2), - 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, gap_fuzz, &b_mask_t, NULL, NULL); + 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, gap_fuzz, &b_mask_t, NULL, NULL, &uopt); } } else if(ha_opt_hic(&asm_opt)) diff --git a/Process_Read.cpp b/Process_Read.cpp index ef99c2e..64f8933 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -541,31 +541,27 @@ void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID) r->length = Get_READ_LENGTH((*R_INF), ID); uint8_t* src = Get_READ((*R_INF), ID); - if (r->length + 4 > r->size) - { + if (r->length + 4 > r->size) { r->size = r->length + 4; r->seq = (char*)realloc(r->seq,sizeof(char)*(r->size)); } uint64_t i = 0; - while ((long long)i < r->length) - { - memcpy(r->seq+i, bit_t_seq_table[src[i>>2]], 4); - i = i + 4; - } - - - if (R_INF->N_site[ID]) - { - for (i = 1; i <= R_INF->N_site[ID][0]; i++) - { - r->seq[R_INF->N_site[ID][i]] = 'N'; + if(src) { + while ((long long)i < r->length) { + memcpy(r->seq+i, bit_t_seq_table[src[i>>2]], 4); + i = i + 4; } + + if (R_INF->N_site[ID]) { + for (i = 1; i <= R_INF->N_site[ID][0]; i++) r->seq[R_INF->N_site[ID][i]] = 'N'; + } + } else {///N + memset(r->seq, 'N', r->length); } - r->RID = ID; - + r->RID = ID; } void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID) @@ -573,8 +569,7 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID) r->length = Get_READ_LENGTH((*R_INF), ID); uint8_t* src = Get_READ((*R_INF), ID); - if (r->length + 4 > r->size) - { + if (r->length + 4 > r->size) { r->size = r->length + 4; r->seq = (char*)realloc(r->seq,sizeof(char)*(r->size)); } @@ -583,27 +578,29 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID) long long i = r->length / 4 - 1 + (last_chr != 0); long long index = 0; - if(last_chr!=0) - { - memcpy(r->seq + index, bit_t_seq_table_rc[src[i]] + 4 - last_chr, last_chr); - index = last_chr; - i--; - } - - while (i >= 0) - { - memcpy(r->seq + index, bit_t_seq_table_rc[src[i]], 4); - i--; - index = index + 4; - } - - if (R_INF->N_site[ID]) - { - for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) - { - r->seq[r->length - R_INF->N_site[ID][i] - 1] = 'N'; + if(src) { + if(last_chr!=0) { + memcpy(r->seq + index, bit_t_seq_table_rc[src[i]] + 4 - last_chr, last_chr); + index = last_chr; + i--; } + + while (i >= 0) { + memcpy(r->seq + index, bit_t_seq_table_rc[src[i]], 4); + i--; + index = index + 4; + } + + if (R_INF->N_site[ID]) { + for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) { + r->seq[r->length - R_INF->N_site[ID][i] - 1] = 'N'; + } + } + } else {///N + memset(r->seq, 'N', r->length); } + + } #define COMPRESS_BASE {c = seq_nt6_table[(uint8_t)src[i]];\ @@ -1893,4 +1890,24 @@ int64_t load_compress_base_disk(FILE *fp, uint64_t *ul_rid, char *dest, uint32_t des_i += tailLen; src_i += tailLen; } return 1; +} + + +scaf_res_t *init_scaf_res_t(uint32_t n) +{ + scaf_res_t *p = NULL; CALLOC(p, 1); + p->n = p->m = n; CALLOC(p->a, n); + return p; +} + +void destroy_scaf_res_t(scaf_res_t *p) +{ + if(p) { + uint32_t k; + for (k = 0; k < p->m; k++) { + free(p->a[k].N_site.a); free(p->a[k].r_base.a); free(p->a[k].bb.a); + } + free(p->a); + free(p); + } } \ No newline at end of file diff --git a/Process_Read.h b/Process_Read.h index 71c88a1..2e1acbd 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -25,6 +25,7 @@ #define Get_NAME(R_INF, ID) ((R_INF).name + (R_INF).name_index[(ID)]) #define CHECK_BY_NAME(R_INF, NAME, ID) (Get_NAME_LENGTH((R_INF),(ID))==strlen((NAME)) && \ memcmp((NAME), Get_NAME((R_INF), (ID)), Get_NAME_LENGTH((R_INF),(ID))) == 0) +#define IS_SCAF_READ(R_INF, ID) ((R_INF).read_sperate[(ID)] == NULL) extern uint8_t seq_nt6_table[256]; extern char bit_t_seq_table[256][4]; @@ -204,6 +205,13 @@ typedef struct // uint32_t mm; } all_ul_t; + +typedef struct { + ul_vec_t *a; + size_t n, m; + uint8_t dd; +} scaf_res_t; + extern all_ul_t UL_INF; extern all_ul_t ULG_INF; // extern uint32_t *het_cnt; @@ -239,5 +247,7 @@ uint64_t retrieve_r_cov_region(const ul_idx_t *ul, uint64_t id, uint8_t strand, void append_ul_t_back(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, int64_t str_l, ul_ov_t *o, int64_t on, float p_chain_rate); void write_compress_base_disk(FILE *fp, uint64_t ul_rid, char *str, uint32_t len, ul_vec_t *buf); int64_t load_compress_base_disk(FILE *fp, uint64_t *ul_rid, char *dest, uint32_t *len, ul_vec_t *buf); +scaf_res_t *init_scaf_res_t(uint32_t n); +void destroy_scaf_res_t(scaf_res_t *p); #endif diff --git a/hic.cpp b/hic.cpp index 28d95db..7d46e88 100644 --- a/hic.cpp +++ b/hic.cpp @@ -6934,7 +6934,7 @@ void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, u if(check_het) { get_bubbles(bub, b_id0, &beg, &sink, NULL, NULL, NULL); - if(IF_HET(beg>>1, *bub) && IF_HET(sink>>1, *bub)) b_id0 = (uint64_t)-1; + if(((beg == (uint32_t)-1) || IF_HET(beg>>1, *bub)) && ((sink == (uint32_t)-1) || IF_HET(sink>>1, *bub))) b_id0 = (uint64_t)-1; } } @@ -6944,7 +6944,7 @@ void get_bub_id(bubble_type* bub, uint32_t root, uint64_t* id0, uint64_t* id1, u if(check_het) { get_bubbles(bub, b_id1, &beg, &sink, NULL, NULL, NULL); - if(IF_HET(beg>>1, *bub) && IF_HET(sink>>1, *bub)) b_id1 = (uint64_t)-1; + if(((beg == (uint32_t)-1) || IF_HET(beg>>1, *bub)) && ((sink == (uint32_t)-1) || IF_HET(sink>>1, *bub))) b_id1 = (uint64_t)-1; } } diff --git a/inter.cpp b/inter.cpp index 41885c0..1b49004 100644 --- a/inter.cpp +++ b/inter.cpp @@ -395,8 +395,20 @@ typedef struct { // data structure for each step in kt_pipeline() // glchain_t *sec_ll; uint64_t num_bases, num_corrected_bases, num_recorrected_bases; int64_t n_thread; + scaf_res_t *rsc; } utepdat_t; +typedef struct { // global data structure for kt_pipeline() + utepdat_t *s; + ma_ug_t *qry; + scaf_res_t *qry_sc; + ma_ug_t *ref; + scaf_res_t *ref_sc; + ma_ug_t *gfa; + bubble_type *bub; + kv_u_trans_t *ta; +} ctdat_t; + #define ha_mzl_t_key(p) ((p).x) KRADIX_SORT_INIT(ha_mzl_t_srt, ha_mzl_t, ha_mzl_t_key, member_size(ha_mzl_t, x)) @@ -5958,6 +5970,53 @@ st_mt_t *bf, Chain_Data* dp, int64_t max_skip, int64_t max_iter, int64_t max_dis +int64_t ctg_chain_graph(void *km, const ul_idx_t *uref, const ma_ug_t *ug, vec_mg_lchain_t *lc, +vec_mg_lchain_t *sw, vec_mg_path_dst_t *dst, vec_sp_node_t *out, vec_mg_pathv_t *path, +int64_t qlen, const ug_opt_t *uopt, int64_t bw, double ng_diff_thre, uint64_t *srt, +st_mt_t *bf, Chain_Data* dp, int64_t max_skip, int64_t max_iter, int64_t max_dis) +{ + bf->n = 0; + if(lc->n == 0) return 0; + int64_t i, j, lc_n = lc->n, x, m_idx; mg_lchain_t *r; + + + for (i = 0; i < lc_n; i++) { + r = &lc->a[i]; r->dist_pre = -1; + srt[i] = r->qe; srt[i] <<= 32; srt[i] |= (uint64_t)i; + } + radix_sort_gfa64(srt, srt+lc_n); + for (i = 1, j = 0; i <= lc_n; i++) { + if (i == lc_n || (srt[i]>>32) != (srt[j]>>32)) { + if(i - j > 1) { + for (x = j; x < i; x++) { + srt[x] <<= 32; srt[x] >>= 32; srt[x] |= ((uint64_t)lc->a[(uint32_t)srt[x]].qs)<<32; + } + radix_sort_gfa64(srt+j, srt+i); + } + j = i; + } + } + + + kv_resize(mg_lchain_t, *sw, (uint64_t)lc_n); sw->n = lc_n; + for (i = 0; i < lc_n; i++) sw->a[i] = lc->a[(uint32_t)srt[i]]; + memcpy(lc->a, sw->a, lc_n *sizeof((*(lc->a)))); + + + resize_Chain_Data(dp, lc_n, NULL); + int32_t *p; int64_t *f, *t; + t = dp->tmp; p = dp->score; f = dp->pre; + + m_idx = gl_chain_linear(uref, ug, lc, sw, qlen, uopt, bw, ng_diff_thre, bf, + max_skip, max_iter, max_dis, p, f, t, lc->n); + + // if(m_idx >= 0 && primary_chain_check(bf->a, bf->n, lc->a)) return m_idx; + // else m_idx = -1; + // bf->n = 0; + return m_idx; +} + + 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) @@ -6356,6 +6415,149 @@ const asg_t *g, st_mt_t *dst_done, vec_sp_node_t *out, vec_mg_pathv_t *res, vec_ return n_mchain; } +inline uint32_t is_compat_contig_chain(const asg_t *g, ul_ov_t *p, ul_ov_t *m, mg_lchain_t *a, double overhead) +{ + uint64_t ovlp; + ovlp = ((MIN(m->qe, p->qe) > MAX(m->qs, p->qs))? (MIN(m->qe, p->qe) - MAX(m->qs, p->qs)):0); + if(ovlp == 0) return 1; + if(ovlp == (m->qe - m->qs)) return 0;////contain + if(ovlp == (p->qe - p->qs)) return 0;////contain + ul_ov_t *i0, *i1; + if(p->qs <= m->qs) { + i0 = p; i1 = m; + } else { + i0 = m; i1 = p; + } + + uint64_t o0 = max_ovlp(g,a[i0->te-1].v); + uint64_t o1 = max_ovlp(g,a[i1->ts].v^1); + uint64_t cut = MIN(o0, o1); cut *= (1 + overhead); + if(ovlp > cut) return 0; + return 1; +} + + +uint32_t select_max_ctg_chain(void *km, const ul_idx_t *uref, int64_t ulid, st_mt_t *idx, vec_mg_lchain_t *e, kv_ul_ov_t *raw_idx, +const asg_t *g, st_mt_t *dst_done, vec_sp_node_t *out, vec_mg_pathv_t *res, vec_mg_lchain_t *gchains) +{ + gchains->n = 0; + if(idx->n <= 0) return 0; + int64_t a_n, idx_n = idx->n, i, k, n_mchain = 0, min_sc, max_sc; uint64_t om, ok, ovlp; + ul_ov_t *m = NULL, *p = NULL, kp; mg_lchain_t *a = e->a, *g_item; int64_t raw_idx_n = raw_idx->n; + ul_ov_t *gb = NULL; int64_t gb_n = 0; + for (i = a_n = 0; i < idx_n; ++i) { + kv_pushp(ul_ov_t, *raw_idx, &p); + p->qn = ((int64_t)(idx->a[i]>>32));//score + p->ts = a_n; p->te = a_n + ((uint32_t)idx->a[i]); + p->qs = a[p->ts].qs; p->qe = a[p->te-1].qe; p->tn = 1;//(tn = 1) -> normal; (t = 0) -> duplicated chain + a_n += ((uint32_t)idx->a[i]); + } + gb = raw_idx->a + raw_idx_n; gb_n = raw_idx->n - raw_idx_n; + radix_sort_ul_ov_srt_qn(gb, gb + gb_n);//sort by scores + for (k = 0, n_mchain = gb_n>>1; k < n_mchain; k++) { + kp = gb[k]; gb[k] = gb[gb_n-k-1]; gb[gb_n-k-1] = kp; + } + + for (k = 0; k < gb_n; k++) {//filter too close chains + m = &(gb[k]); om = m->qe - m->qs; ///current chain + // fprintf(stderr, "k::%ld[M::%s::sc->%u] q::[%u, %u), set::%u\n", k, __func__, m->qn, m->qs, m->qe, m->tn); + if(m->tn == 0) continue; + for (i = k-1; i >= 0; i--) { + p = &(gb[i]); + ovlp = ((MIN(m->qe, p->qe) > MAX(m->qs, p->qs))? (MIN(m->qe, p->qe) - MAX(m->qs, p->qs)):0); + if(ovlp == 0) continue; + min_sc = MIN(p->qn, m->qn); max_sc = MAX(p->qn, m->qn); + ok = p->qe - p->qs; ok = MAX(ok, om); + if(min_sc < (max_sc*0.98)) break; + if((ovlp > GC_OFFSET_POS) && (min_sc > (max_sc*0.98)) && (ovlp > (ok*0.8))) { + // fprintf(stderr, "k::%ld[M::%s::i->%ld] min_sc::%ld, max_sc::%ld\n", + // k, __func__, i, min_sc, max_sc); + m->tn = p->tn = 0; + } + } + + for (i = k+1; i < gb_n; i++) { + p = &(gb[i]); + ovlp = ((MIN(m->qe, p->qe) > MAX(m->qs, p->qs))? (MIN(m->qe, p->qe) - MAX(m->qs, p->qs)):0); + if(ovlp == 0) continue; + min_sc = MIN(p->qn, m->qn); max_sc = MAX(p->qn, m->qn); + ok = p->qe - p->qs; ok = MAX(ok, om); + if(min_sc < (max_sc*0.98)) break; + if((ovlp > GC_OFFSET_POS) && (min_sc > (max_sc*0.98)) && (ovlp > (ok*0.8))) { + // fprintf(stderr, "k::%ld[M::%s::i->%ld] min_sc::%ld, max_sc::%ld\n", + // k, __func__, i, min_sc, max_sc); + m->tn = p->tn = 0; + } + } + } + + for (k = i = 0; k < gb_n; k++) { + m = &(gb[k]); if(m->tn == 0) continue; + gb[i++] = gb[k]; + } + // fprintf(stderr, "[M::%s::] gb_n0::%ld, gb_n::%ld\n", __func__, gb_n, i); + gb_n = i; + + for (k = n_mchain = 0; k < gb_n; k++) { + m = &(gb[k]); + for (i = 0; i < n_mchain; i++) { + p = &(gb[i]); + if(!is_compat_contig_chain(g, p, m, a, 0.333333)) break; + } + if(i < n_mchain) continue; + gb[n_mchain++] = gb[k]; + } + gb_n = n_mchain; + + if(gb_n) { + gchains->n = 0; + for (k = 0; k < gb_n; k++) { + kv_pushp(mg_lchain_t, *gchains, &g_item); + g_item->v = (uint32_t)-1; + g_item->qs = gb[k].qs; g_item->qe = gb[k].qe; + g_item->rs = gb[k].ts; g_item->re = gb[k].te; + g_item->cnt = g_item->off = 0; + gen_gchain_track(km, a + g_item->rs, g_item->re - g_item->rs, g, dst_done, out, res); + kv_resize(mg_lchain_t, *gchains, gchains->n + res->n); ///a = gchains->a + gchains->n; + for (i = 0, g_item = &(gchains->a[gchains->n-1]); i < ((int64_t)res->n); i++) { + if(res->a[i].v == (uint32_t)-1) { + gchains->a[i+gchains->n] = a[res->a[i].pre + g_item->rs]; + gchains->a[i+gchains->n].dist_pre = res->a[i].d; + // if(ulid == 14714) { + // fprintf(stderr, "+[M::%s::]\tutg%.6dl(%c)\tq::[%d, %d)\tt::[%d, %d)\n", __func__, + // (int32_t)(gchains->a[i+gchains->n].v>>1)+1, "+-"[gchains->a[i+gchains->n].v&1], + // gchains->a[i+gchains->n].qs, gchains->a[i+gchains->n].qe, + // gchains->a[i+gchains->n].rs, gchains->a[i+gchains->n].re); + // } + } else { + gchains->a[i+gchains->n].v = res->a[i].v; + gchains->a[i+gchains->n].off = -1; + gchains->a[i+gchains->n].dist_pre = res->a[i].d; + ///the nodes detected by the graph chaining should be fully covered + gchains->a[i+gchains->n].rs = 0; + gchains->a[i+gchains->n].re = uref->ug->g->seq[res->a[i].v>>1].len; + // fprintf(stderr, "aaaaaaa, ulid->%ld\n", ulid); + // fprintf(stderr, "-[M::%s::]\tutg%.6dl(%c)\n", __func__, + // (int32_t)(gchains->a[i+gchains->n].v>>1)+1, "+-"[gchains->a[i+gchains->n].v&1]); + // if(ulid == 14714) { + // fprintf(stderr, "-[M::%s::]\tutg%.6dl(%c)\tq::[%d, %d)\tt::[%d, %d)\n", __func__, + // (int32_t)(gchains->a[i+gchains->n].v>>1)+1, "+-"[gchains->a[i+gchains->n].v&1], + // gchains->a[i+gchains->n].qs, gchains->a[i+gchains->n].qe, + // gchains->a[i+gchains->n].rs, gchains->a[i+gchains->n].re); + // } + } + } + g_item->cnt = res->n; + gchains->n += res->n; + // fprintf(stderr, "sbsbsbsb, ulid->%ld\n", ulid); + // debug_gchain(km, g, gchains->a + gchains->n - res->n, res->n, dst_done, out); + } + } + + raw_idx->n = raw_idx_n; + return n_mchain; +} + // void prt_chains_vlog(ul_ov_t *l_idx, int64_t l_idx_n, ul_ov_t *l_a, uint64_t *g_idx, int64_t g_idx_n, vec_mg_lchain_t *g_a, int64_t ql, int64_t ulid) // { // int64_t k, i, s, e; char *as = NULL; @@ -15833,6 +16035,56 @@ int64_t dp_max_skip, int64_t dp_max_iter, int64_t dp_max_dis) dd_ul_vec_t(uref, swap->a, swap->n, rch, ulid); } + +void push_ctg_res(ul_vec_t *rch, vec_mg_lchain_t *uc) +{ + // fprintf(stderr, "\n++[M::%s::%.*s(id:%ld), len:%u]\n", __func__, + // UL_INF.nid.a[ulid].n, UL_INF.nid.a[ulid].a, ulid, rch->rlen); + + uint64_t k, ucn = uc->n; uint64_t tt, s, e, z, rrn; mg_lchain_t *ix; uc_block_t *p; uc_block_t *rr; mg_lchain_t *src; + for (k = 0, rch->bb.n = 0, tt = 0; k < ucn; k += ix->cnt + 1) { + ix = &(uc->a[k]); assert(ix->v == (uint32_t)-1); + kv_pushp(uc_block_t, rch->bb, &p); memset(p, 0, sizeof((*p))); + p->hid = (uint32_t)-1; p->qs = ix->qs; p->qe = ix->qe; p->ts = ix->cnt; tt += ix->cnt; + + // if(ulid == 14714) { + // fprintf(stderr, "\n[M::%s::ucn->%ld, k->%ld, kcnt->%d]\n", __func__, ucn, k, ix->cnt); + // print_debug_gchain(uref, uc->a + k + 1, ix->cnt, rch); + // } + // gen_rovlp_chain_by_ul(rg, rch, uref, raw_idx, raw_chn, uc->a + k + 1, ix->cnt, swap, dp, + // dp_max_skip, dp_max_iter, dp_max_dis, ulid); + } + radix_sort_uc_block_t_qe_srt(rch->bb.a, rch->bb.a + rch->bb.n); + + kv_resize(uc_block_t, rch->bb, rch->bb.n + tt); + memset(rch->bb.a + rch->bb.n, 0, tt * sizeof(*(rch->bb.a))); + for (k = 0, tt = rch->bb.n; k < rch->bb.n; k++) { + s = tt; e = tt + rch->bb.a[k].ts; + rch->bb.a[k].ts = s; rch->bb.a[k].te = e; + tt = e; + } + assert(tt <= rch->bb.m); + + for (k = 0, tt = 0; k < ucn; k += ix->cnt + 1, tt++) { + ix = &(uc->a[k]); assert(ix->v == (uint32_t)-1); + rr = rch->bb.a + rch->bb.a[tt].ts; rrn = rch->bb.a[tt].te - rch->bb.a[tt].ts; + src = uc->a + k + 1; assert((uint32_t)ix->cnt == rrn); + for(z = 0; z < rrn; z++) { + rr[z].hid = src[z].v>>1; + rr[z].rev = src[z].v&1; + rr[z].qs = src[z].qs; rr[z].qe = src[z].qe; + rr[z].ts = src[z].rs; rr[z].te = src[z].re; + rr[z].aidx = k; + } + } + ///up to now, given a in swap + ///x->ts and x->te are the coordinates in unitig (x->score>>1), instead of HiFi read (x->v>>1) + ///x->qs and x->qe are the coordinates in UL + ///x->dist_pre is the idx of this chain at rch + // debug_intermediate_chain(uref->ug, swap->a, swap->n); + // dd_ul_vec_t(uref, swap->a, swap->n, rch, ulid); +} + void print_ru_raw_chains(kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn, vec_mg_lchain_t *gch, ul_vec_t *rch, ma_ug_t *ug) { int64_t k, i; uint64_t ts, te, qs, qe; @@ -15952,6 +16204,713 @@ int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid, Chain_Data* dp, c } + +uint64_t rov2uov_ctg(uint64_t rid, const ul_idx_t *uref, utg_rid_dt *ru_map, uc_block_t *rovlp, ul_ov_t *res, uint32_t qidx, uint32_t adjust_rev) +{ + uint64_t ori = ru_map->u&1, ts, te; + if(!ori) { + ts = rovlp->ts; te = rovlp->te; + } else { + ts = Get_READ_LENGTH(R_INF, rid) - rovlp->te; + te = Get_READ_LENGTH(R_INF, rid) - rovlp->ts; + } + + ts += ru_map->pos; te += ru_map->pos; + if(ts >= 0 && te <= uref->ug->g->seq[ru_map->u>>1].len) { + memset(res, 0, sizeof(*res)); + res->qn = qidx; res->qs = rovlp->qs; res->qe = rovlp->qe; + res->tn = ru_map->u>>1; res->ts = ts; res->te = te; + res->el = rovlp->el; res->rev = (rovlp->rev == ori?0:1); res->sec = ru_map->off; + if(adjust_rev && res->rev) {///for linear chaining + res->ts = uref->ug->g->seq[ru_map->u>>1].len - te; + res->te = uref->ug->g->seq[ru_map->u>>1].len - ts; + res->sec = uref->ug->u.a[ru_map->u>>1].n - ru_map->off; + } + return 1; + } + return 0; +} + + +uint64_t get_unique_rctg_aln(const ul_idx_t *uref, uint64_t v, uint64_t vi, uint64_t vs, uint64_t ulen, ul_ov_t *ures, kv_ul_ov_t *rr, uint32_t sec) +{ + utg_rid_dt *a = NULL; uint64_t a_n, k, rrn = 0; uc_block_t z; ul_ov_t p; + p.tn = p.qn = ((uint32_t)-1); if(ures) ures->tn = ures->qn = ((uint32_t)-1); + if(IS_SCAF_READ(R_INF, (v>>1))) return 0;///actually no scaf nodes within the assembly graph + a = get_r_ug_region(uref->r_ug, &a_n, v>>1); + if(!a) return 0; + + + z.ts = 0; z.te = Get_READ_LENGTH(R_INF, (v>>1)); + z.qs = vs; z.qe = vs + Get_READ_LENGTH(R_INF, (v>>1)); + if(z.qs > ulen) z.qs = ulen; if(z.qe > ulen) z.qe = ulen; + z.hid = v>>1; z.rev = v&1; z.el = 1; + + + for (k = rrn = 0; k < a_n; k++) { + if(!rov2uov_ctg(v>>1, uref, &(a[k]), &z, &p, vi, 1)) continue; + p.el = 1; p.tn <<= 1; p.tn |= p.rev; ///p.qn: offset within q; p.sec: offset within t; p.tn: tid + if(sec != ((uint32_t)-1)) p.sec = sec; + if(rr) kv_push(ul_ov_t, *rr, p); + rrn++; + } + + if((ures) && (rrn == 1)) *ures = p; + return rrn; +} + +uint64_t cmp_shrink_ul_ov_t(ul_ov_t *a0, ul_ov_t *a1, uint64_t off) +{ + if((a0->tn == ((uint32_t)-1)) || (a1->tn == ((uint32_t)-1))) return 1; + if((a0->tn == a1->tn) && (a0->qn + off == a1->qn) && (a0->sec + off == a1->sec)) return 0; + return 1; +} + +void ctg_rg2ug_gen(ma_utg_t *r_cl, kv_ul_ov_t *u_cl, const ul_idx_t *uref) +{ + if(r_cl->n <= 0) return ; + uint64_t k, i, m, l, is_push, vs, ls; ul_ov_t kp, lp, kp0; + u_cl->n = 0; kp.tn = kp.qn = kp0.tn = kp0.qn = ((uint32_t)-1); + + k = 1; l = 0; ls = 0; vs = 0; + get_unique_rctg_aln(uref, r_cl->a[l]>>32, l, ls, r_cl->len, &lp, NULL, ((uint32_t)-1)); ///get lp for l = 0 + vs += (uint32_t)(r_cl->a[l]); + + for (; k <= r_cl->n; k++) { + is_push = 0; + if(k < r_cl->n) { + get_unique_rctg_aln(uref, r_cl->a[k]>>32, k, vs, r_cl->len, &kp, NULL, ((uint32_t)-1)); + if(cmp_shrink_ul_ov_t(&lp, &kp, k - l)) is_push = 1; + } else { + is_push = 1; + } + + if(is_push) { + kp0 = kp; + if(k - l > 1) { + for (i = l; i < k; i++) {///could be merged + m = get_unique_rctg_aln(uref, r_cl->a[i]>>32, i, ls, r_cl->len, &kp, NULL, ((uint32_t)-1)); + assert(m == 1); assert(kp.tn == lp.tn); assert(kp.qn == lp.qn + i - l); assert(kp.sec == lp.sec + i - l); + if(kp.qs < lp.qs) lp.qs = kp.qs; if(kp.qe > lp.qe) lp.qe = kp.qe; + if(kp.ts < lp.ts) lp.ts = kp.ts; if(kp.te > lp.te) lp.te = kp.te; + ls += (uint32_t)(r_cl->a[i]); + } + lp.sec = k - l; + kv_push(ul_ov_t, *u_cl, lp); + } else if(k - l == 1) { + m = get_unique_rctg_aln(uref, r_cl->a[l]>>32, l, ls, r_cl->len, NULL, u_cl, 1); + } + + l = k; lp = kp0; ls = vs; + } + + if(k < r_cl->n) vs += (uint32_t)(r_cl->a[k]); + } +} + +inline int64_t comput_linear_ctg_sc(ul_ov_t *li, ul_ov_t *lj, double diff_ec_ul, int64_t bw) +{ ///li is the suffix of lj + int64_t dq, dt, dd, mm; + if(lj->te > li->te) return INT32_MIN; + dq = li->qe - lj->qs; dt = li->te - lj->ts; + dd = (dq>dt? dq-dt:dt-dq); + mm = MAX(dq, dt); mm *= diff_ec_ul; if(mm < bw) mm = bw; + if(dd > mm) return INT32_MIN; + return li->sec; +} + +inline uint64_t is_available_ctg_aln(const ul_idx_t *uref, ul_ov_t *z) +{ + if((z->te > z->ts) && (uref->ug->u.a[z->tn].len >= (z->te-z->ts)) && ((uref->ug->u.a[z->tn].len) - (z->te-z->ts) <= 16) && (z->sec >= (uref->ug->u.a[z->tn].n*0.333333))) { + return 1; + } + return 0; +} + +uint64_t linear_ctg_chain_dp_adv(ul_ov_t *ch, int64_t ch_n, ul_ov_t *sv, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, +double diff_ec_ul, int64_t qlen, int64_t max_skip, int64_t max_iter, int64_t max_dis, Chain_Data* dp, ma_ug_t *ug, +int64_t chain_offset) +{ ///all in[].el must be 1 + if(ch_n == 0) return 0; + int64_t i, j, k, sc, csc, mm_sc, mm_idx, its, ite, max; + ul_ov_t *li = NULL, *lj = NULL; int64_t *p, *t, st, plus, max_ii, n_skip, end_j; int32_t *f; + resize_Chain_Data(dp, ch_n, NULL); + t = dp->tmp; f = dp->score; p = dp->pre; + + radix_sort_ul_ov_srt_qe(ch, ch + ch_n); + for (i = 1, j = 0; i <= ch_n; i++) { + if (i == ch_n || ch[i].qe != ch[j].qe) { + if(i - j > 1) radix_sort_ul_ov_srt_qs(ch+j, ch+i); + j = i; + } + } + + // fprintf(stderr, "[M::%s::] ch_n:%ld\n", __func__, ch_n); + memset(t, 0, (ch_n*sizeof((*t)))); + for (i = st = plus = 0, max_ii = -1; i < ch_n; ++i) { + li = &(ch[i]); csc = li->sec; + mm_sc = csc; mm_idx = -1; n_skip = 0; end_j = -1; + st = (i= st; --j) { + lj = &(ch[j]); + sc = comput_linear_ctg_sc(li, lj, diff_ec_ul, bw); ///should allow contain + if(sc == INT32_MIN) continue; + sc += f[j]; + if(sc > mm_sc) { + mm_sc = sc, mm_idx = j; + if (n_skip > 0) --n_skip; + } else if (t[j] == i) { + if (++n_skip > max_skip) + break; + } + if (p[j] >= 0) t[p[j]] = i; + } + end_j = j; + if (max_ii < 0 || (ch[i].qe>(ch[max_ii].qe+max_dis))) {//too long + max = INT32_MIN; max_ii = -1; + for (j = i - 1; (j >= st) && (ch[i].qe<=(max_dis+ch[j].qe)); --j) { + if (max < f[j]) { + max = f[j], max_ii = j; + } + } + } + + if (max_ii >= 0 && max_ii < end_j) {///just have a try with a[i]<->a[max_ii] + lj = &(ch[max_ii]); + sc = comput_linear_ctg_sc(li, lj, diff_ec_ul, bw); ///should allow contain + if(sc != INT32_MIN) { + sc += f[max_ii]; + if(sc > mm_sc) { + mm_sc = sc; mm_idx = max_ii; + } + } + } + + f[i] = mm_sc; p[i] = mm_idx; + if ((max_ii < 0) || ((ch[i].qe<=max_dis+ch[max_ii].qe) && (f[max_ii]sec = (mm_idx<0?0x3FFFFFFF:i-mm_idx); + sv[i] = *li; + // fprintf(stderr, "##(%ld) %u\t%u\t%c\tutg%.6d%c(%u)\t%u\t%u\tmm_idx:%ld\tmm_sc:%ld\n", i, li->qs, li->qe, "+-"[li->rev], + // (int32_t)(li->tn)+1, "lc"[uref->ug->u.a[li->tn].circ], uref->ug->u.a[li->tn].len, li->ts, li->te, mm_idx, mm_sc); + // track[i] = push_sc_pre(mm_sc, mm_idx); + // li->sec = (mm_idx<0?0x3FFFFFFF:i-mm_idx); sv[i] = *li; + } + + for (i = 0; i < ch_n; ++i) t[i] = 0; + int64_t n_u; + for (k = ch_n-1, n_u = 0; k >= 0; --k) { + if(t[k]) continue; + i = k; ch[n_u]=sv[i]; sc = f[i]; + for (;i>=0;) { + if(sv[i].qs < ch[n_u].qs) ch[n_u].qs = sv[i].qs; + if(sv[i].ts < ch[n_u].ts) ch[n_u].ts = sv[i].ts; + if(sv[i].qe > ch[n_u].qe) ch[n_u].qe = sv[i].qe; + if(sv[i].te > ch[n_u].te) ch[n_u].te = sv[i].te; + // ch[n_u].qn = i;//start idx of read alignment in chain + t[i] = 1; i = p[i]; + } + adjust_rev_tse(&(ch[n_u]), ug->g->seq[ch[n_u].tn].len, &its, &ite); + ch[n_u].ts = its; ch[n_u].te = ite; ch[n_u].sec = (sc>0x3FFFFFFF?0x3FFFFFFF:sc); + // ch[n_u].qn += chain_offset; //start idx of read alignment in chain + // ch[n_u].tn = k + chain_offset; //end idx of read alignment in chain + ch[n_u].qn = k + chain_offset; //end idx of read alignment in chain + + + n_u++; + } + for (i = 0; i < ch_n; ++i) { + adjust_rev_tse(&(sv[i]), ug->g->seq[sv[i].tn].len, &its, &ite); + sv[i].ts = its; sv[i].te = ite; + k = p[i]; sv[i].tn = k>=0?k+chain_offset:(uint32_t)-1; + } + return n_u; +} + +void gen_linear_chains_ctg(kv_ul_ov_t *res, kv_ul_ov_t *buf, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, +double diff_ec_ul, int64_t qlen, Chain_Data* dp) +{ + uint64_t k, l, z, an, m; + radix_sort_ul_ov_srt_tn(res->a, res->a + res->n); + ///after this function, res keeps unitig alignment, while buf keeps read alignments + kv_resize(ul_ov_t, *buf, res->n); buf->n = res->n; + for (k = 1, l = m = 0; k <= res->n; k++) { + if(k == res->n || res->a[k].tn != res->a[l].tn) {///qn <- (tn|rev) + for (z = l; z < k; z++) res->a[z].tn>>=1; + an = k; + if(k - l > 1) {///if there is only one alignment for a node, then no need for DP + an = l + linear_ctg_chain_dp_adv(res->a+l, k-l, buf->a+l, uref, uopt, bw, diff_ec_ul, qlen, UG_SKIP_N, UG_ITER_N, UG_DIS_N, dp, uref->ug, l); + } + + for (z = l; z < an; z++) { + if(is_available_ctg_aln(uref, &(res->a[z]))) { + res->a[z].ts = 0; res->a[z].te = uref->ug->u.a[res->a[z].tn].len; + res->a[m++] = res->a[z]; + } + } + l = k; + } + } + res->n = m; +} + + + + +///sps and hap are just vector for uint64_t; used for buffer +uint32_t direct_gchain_scaf(mg_tbuf_t *b, ma_utg_t *rch, glchain_t *ll, gdpchain_t *gdp, st_mt_t *sps, haplotype_evdience_alloc *hap, const ul_idx_t *uref, const ug_opt_t *uopt, +int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid, Chain_Data* dp, const asg_t *rg, ul_vec_t *res) +{ + res->bb.n = 0; + // if(ulid != 86660) return 0; + kv_ul_ov_t *idx = &(ll->lo), *init = &(ll->tk); int64_t max_idx; + idx->n = init->n = 0; + // gl_rg2ug_gen(rch, idx, uref, 1, 2, ulid); + ctg_rg2ug_gen(rch, idx, uref); + + // uint32_t k; + // fprintf(stderr, "[M::%s] rch->len::%u, rch->n::%u, idx->n::%u\n", __func__, (uint32_t)rch->len, (uint32_t)rch->n, (uint32_t)idx->n); + // for (k = 0; k < idx->n; k++) { + // fprintf(stderr, "[k->%u::utg%.6u%c(len->%u::n->%u)]\tq::[%u, %u)\t%c\tt::[%u, %u)\n", + // k, (idx->a[k].tn>>1) + 1, "lc"[uref->ug->u.a[(idx->a[k].tn>>1)].circ], uref->ug->u.a[(idx->a[k].tn>>1)].len, idx->a[k].sec, + // idx->a[k].qs, idx->a[k].qe, "+-"[idx->a[k].rev], idx->a[k].ts, idx->a[k].te); + // } + + // for (k = 0; k < idx->n; k++) { + // fprintf(stderr, "utg%.6u%c,", (idx->a[k].tn>>1)+1, "lc"[uref->ug->u.a[(idx->a[k].tn>>1)].circ]); + // } + // fprintf(stderr, "\n"); + + if(idx->n == 0) return 0; + ///generate linear chains + gen_linear_chains_ctg(idx, init, uref, uopt, bw, diff_ec_ul, rch->len, dp); + + if(idx->n == 0) return 0; + + dump_linear_chain(uref->ug, idx, &(gdp->l), rch->len, ulid); + if(gdp->l.n == 0) return 0; + + kv_resize(uint64_t, ll->srt.a, gdp->l.n); + max_idx = ctg_chain_graph(b->km, uref, uref->ug, &(gdp->l), &(gdp->swap), &(gdp->dst), &(gdp->out), + &(gdp->path), rch->len, uopt, G_CHAIN_BW, N_GCHAIN_RATE, ll->srt.a.a, sps, dp, UG_SKIP_GRAPH_N, UG_ITER_N, /**UG_DIS_N**/UG_DIS_N*100); + + //sps -> idx; (gdp->l) -> alignments + if(max_idx >= 0) { + max_idx = select_max_ctg_chain(b->km, uref, ulid, sps, &(gdp->l), &(ll->tk), uref->ug->g, &(gdp->dst_done), &(gdp->out), &(gdp->path), &(gdp->swap)); + if(max_idx) { + push_ctg_res(res, &(gdp->swap)); + return 1; + } + } + + return 0; +} + +uint64_t get_seq_bub_id(uint64_t v, bubble_type *bub, uint64_t *bid, uint64_t *nxt) +{ + uint64_t id0, id1, id; uint32_t beg, sink; + + (*bid) = (*nxt) = (uint64_t)-1; + get_bub_id(bub, v>>1, &id0, &id1, 1); + + id = id0; + if((id != (uint64_t)-1) && (id < bub->f_bub)) {///no need broken bubble + get_bubbles(bub, id, &beg, &sink, NULL, NULL, NULL); + if(v == beg) { + (*bid) = id; (*nxt) = sink^1; return 1; + } + + if(v == sink) { + (*bid) = id; (*nxt) = beg^1; return 1; + } + } + + id = id1; + if((id != (uint64_t)-1) && (id < bub->f_bub)) {///no need broken bubble + get_bubbles(bub, id, &beg, &sink, NULL, NULL, NULL); + if(v == beg) { + (*bid) = id; (*nxt) = sink^1; return 1; + } + + if(v == sink) { + (*bid) = id; (*nxt) = beg^1; return 1; + } + } + + return 0; +} + +uint64_t pick_bubble(uc_block_t *a, uint64_t a_n, bubble_type *bub, uint64_t *r_bid) +{ + uint64_t k = 0, v, bid = (uint64_t)-1, w = (uint64_t)-1; (*r_bid) = (uint64_t)-1; + if(a_n <= 1) return 0;///no bubble + v = (((uint32_t)a[k].hid)<<1)|((uint32_t)a[k].rev); + if(!get_seq_bub_id(v, bub, &bid, &w)) return 0; + (*r_bid) = bid; + for (k = 1; k < a_n; k++) { + v = (((uint32_t)a[k].hid)<<1)|((uint32_t)a[k].rev); + if(v == w) return k; + if(bub->index[v>>1] != bid) return 0; + } + + return 0; +} + +uint64_t qry_bub_sc(const ul_idx_t *uref, ma_ug_t *gfa, scaf_res_t *ref_sc, bubble_type *bub, uint64_t bid, uint64_t qsidx, uint64_t qeidx, ul_vec_t *qstr, kv_ul_ov_t *res) +{ + // fprintf(stderr, "qsidx::%lu(utg%.6u%c), qeidx::%lu(utg%.6u%c)\n", qsidx, (qstr->bb.a[qsidx].hid)+1, "lc"[gfa->u.a[qstr->bb.a[qsidx].hid].circ], qeidx, (qstr->bb.a[qeidx].hid)+1, "lc"[gfa->u.a[qstr->bb.a[qeidx].hid].circ]); + uint32_t beg, sink, v, w, z, si, ei, qrev = (uint32_t)-1, qv, qw, trev, nn = 0, tn; + get_bubbles(bub, bid, &beg, &sink, NULL, NULL, NULL); + if((beg == ((uint32_t)-1)) || (sink == ((uint32_t)-1))) return 0; + utg_rid_dt *a; uint64_t a_n, i, k; ul_vec_t *q; ul_ov_t p; + qv = qstr->bb.a[qsidx].hid; qv <<= 1; qv |= ((uint32_t)qstr->bb.a[qsidx].rev); + qw = qstr->bb.a[qeidx].hid; qw <<= 1; qw |= ((uint32_t)qstr->bb.a[qeidx].rev); + if((qv == beg) && (qw == (sink^1))) qrev = 0; + if((qv == sink) && (qw == (beg^1))) qrev = 1; + assert(qrev != (uint32_t)-1); + memset(&p, 0, sizeof(p)); + + v = beg; w = sink^1; trev = 0; + a = get_r_ug_region(uref->r_ug, &a_n, v>>1); + // fprintf(stderr, "+a_n::%lu\n", a_n); + for (i = 0; i < a_n; i++) { + // fprintf(stderr, "+(0)i::%lu\n", i); + if((a[i].u&1) != (v&1)) continue; + q = &(ref_sc->a[a[i].u>>1]); si = ei = a[i].off; + assert((q->bb.a[si].hid == (v>>1)) && (q->bb.a[si].rev == (v&1))); + // fprintf(stderr, "+(1)i::%lu\n", i); + for (k = tn = 0; k < q->bb.n; k++) tn += q->bb.a[k].te - q->bb.a[k].ts; + for (k = si + 1; (k < tn) && (q->bb.a[k].aidx == q->bb.a[si].aidx); k++) { + z = (((uint32_t)q->bb.a[k].hid)<<1)|((uint32_t)q->bb.a[k].rev); + if(z == w) {ei = k; break;} + if(bub->index[z>>1] != bid) break; + } + + if(ei > si) {////found a bubble + // fprintf(stderr, "+(2)i::%lu, si::%u, ei::%u, ref_id::%u\nn", i, si, ei, a[i].u>>1); + p.qn = q->bb.a[si].aidx; p.qs = qsidx; p.qe = qeidx; + p.tn = a[i].u>>1; p.ts = si; p.te = ei; + p.rev = (qrev == trev?0:1); + kv_push(ul_ov_t, *res, p); + nn++; + } + } + + + v = sink; w = beg^1; trev = 1; + a = get_r_ug_region(uref->r_ug, &a_n, v>>1); + // fprintf(stderr, "-a_n::%lu\n", a_n); + for (i = 0; i < a_n; i++) { + if((a[i].u&1) != (v&1)) continue; + q = &(ref_sc->a[a[i].u>>1]); si = ei = a[i].off; + assert((q->bb.a[si].hid == (v>>1)) && (q->bb.a[si].rev == (v&1))); + for (k = tn = 0; k < q->bb.n; k++) tn += q->bb.a[k].te - q->bb.a[k].ts; + for (k = si + 1; (k < tn) && (q->bb.a[k].aidx == q->bb.a[si].aidx); k++) { + z = (((uint32_t)q->bb.a[k].hid)<<1)|((uint32_t)q->bb.a[k].rev); + if(z == w) {ei = k; break;} + if(bub->index[z>>1] != bid) break; + } + + if(ei > si) {////found a bubble + // fprintf(stderr, "-(2)i::%lu, si::%u, ei::%u, ref_id::%u\nn", i, si, ei, a[i].u>>1); + p.qn = q->bb.a[si].aidx; p.qs = qsidx; p.qe = qeidx; + p.tn = a[i].u>>1; p.ts = si; p.te = ei; + p.rev = (qrev == trev?0:1); + kv_push(ul_ov_t, *res, p); + nn++; + } + } + + return nn; +} + +/** +uint64_t bubble_check_push(uc_block_t *a, uint64_t a_n, bubble_type *bub) +{ + uint64_t k = 0, v, bid, w; + if(a_n <= 1) return 0;///no bubble + while(k < a_n) { + v = (((uint32_t)a[k].hid)<<1)|((uint32_t)a[k].rev); + if(!get_seq_bub_id(v, bub, &bid, &w)) break; + } + + + for (k = 1, l = 0; k <= a_n; k++) { + if(k ) { + l = k; + } + } + + get_bub_id(bub, root_id, &id0, &id1, 1);; +} +**/ + + +uint64_t cc_bub_match(ul_ov_t *ch, int64_t ch_n) { + if(ch_n == 0) return 0; + int64_t i, j, k; + + radix_sort_ul_ov_srt_qs(ch, ch + ch_n); + for (i = 1, j = 0; i <= ch_n; i++) { + if (i == ch_n || ch[i].qs != ch[j].qs) { + if(i - j > 1) radix_sort_ul_ov_srt_qe(ch+j, ch+i); + j = i; + } + } + + for (i = k = 0; i < ch_n; ++i) { + if(ch[i].tn == (uint32_t)-1) continue; + for (j = i + 1; j < ch_n && ch[i].qe >= ch[j].qs; j++) { + if(ch[j].tn == (uint32_t)-1) continue; + if((ch[i].qe == ch[j].qs) && (ch[i].rev == ch[j].rev) && (ch[i].qn == ch[j].qn)) { + if(!ch[i].rev) { + if(ch[i].te == ch[j].ts) { + ch[i].qe = ch[j].qe; + ch[i].te = ch[j].te; + ch[j].tn = (uint32_t)-1; + } + } else { + if(ch[j].te == ch[i].ts) { + ch[i].qe = ch[j].qe; + ch[i].ts = ch[j].ts; + ch[j].tn = (uint32_t)-1; + } + } + } + } + + ch[k++] = ch[i]; + } + + return k; +} + +uint64_t extract_scaf_res_t(uc_block_t *in, const ul_idx_t *uref, scaf_res_t *ref_sc, kv_ul_ov_t *res) +{ + utg_rid_dt *a; uint64_t i, a_n, nn, nw; uc_block_t *ou; ul_ov_t p; memset(&p, 0, sizeof(p)); + a = get_r_ug_region(uref->r_ug, &a_n, in->hid); + for (i = nn = 0; i < a_n; i++) { + ou = &(ref_sc->a[a[i].u>>1].bb.a[a[i].off]); + p.qn = (uint32_t)-1; p.qs = in->qs; p.qe = in->qe; + p.tn = a[i].u>>1; p.ts = ou->qs; p.te = ou->qe; + p.el = 1; p.rev = ((in->rev == ou->rev)?0:1); + nw = MAX(p.te-p.ts, p.qe-p.qs); nw *= CHAIN_MATCH; + p.sec = ((nw>=0x3fffffff)?(0x3fffffff):(nw)); + kv_push(ul_ov_t, *res, p); + + nn++; + } + return nn; +} + +uint64_t gen_trans_ovlp_scaf(ma_ug_t *gfa, u_trans_t *map, uc_block_t *qin, uc_block_t *rin, ul_ov_t *res) +{ + uint64_t s, e, s0, e0; uc_block_t *z; + + z = qin; s0 = map->qs; e0 = map->qe; + if(s0 > gfa->u.a[z->hid].len || e0 > gfa->u.a[z->hid].len) return 0; + if(!(z->rev)) { + s = s0; e = e0; + } else { + s = gfa->u.a[z->hid].len - e0; + e = gfa->u.a[z->hid].len - s0; + } + s += z->qs; e += z->qs; + if(s >= e) return 0; + res->qs = s; res->qe = e; + + + z = rin; s0 = map->ts; e0 = map->te; + if(s0 > gfa->u.a[z->hid].len || e0 > gfa->u.a[z->hid].len) return 0; + if(!(z->rev)) { + s = s0; e = e0; + } else { + s = gfa->u.a[z->hid].len - e0; + e = gfa->u.a[z->hid].len - s0; + } + s += z->qs; e += z->qs; + if(s >= e) return 0; + res->ts = s; res->te = e; + + res->rev = map->rev^(qin->rev^rin->rev); + + return 1; +} + +void extract_trans_scaf_res_t(uc_block_t *qin, const ul_idx_t *uref, scaf_res_t *ref_sc, ma_ug_t *gfa, kv_u_trans_t *ta, kv_ul_ov_t *res) +{ + utg_rid_dt *a; uint64_t i, a_n; uc_block_t *rin; + u_trans_t *u; uint64_t k, u_n; + ul_ov_t p; memset(&p, 0, sizeof(p)); + + // fprintf(stderr, "\n[M::%s] utg%.6u%c, [%u, %u)\n", __func__, qin->hid+1, "lc"[gfa->u.a[qin->hid].circ], qin->qs, qin->qe); + u = u_trans_a(*ta, qin->hid); + u_n = u_trans_n(*ta, qin->hid); + for (k = 0; k < u_n; k++) { + // fprintf(stderr, "+[M::%s] utg%.6u%c->utg%.6u%c, q::[%u, %u), %c, t::[%u, %u)\n", __func__, qin->hid+1, "lc"[gfa->u.a[qin->hid].circ], u[k].tn+1, "lc"[gfa->u.a[u[k].tn].circ], u[k].qs, u[k].qe, "+-"[u[k].rev], u[k].ts, u[k].te); + a = get_r_ug_region(uref->r_ug, &a_n, u[k].tn); + for (i = 0; i < a_n; i++) { + rin = &(ref_sc->a[a[i].u>>1].bb.a[a[i].off]); + if(gen_trans_ovlp_scaf(gfa, &(u[k]), qin, rin, &p)) { + p.tn = a[i].u>>1; p.qn = u[k].f; p.el = 0; p.sec = u[k].nw; + kv_push(ul_ov_t, *res, p); + + // fprintf(stderr, "-[M::%s] q::[%u, %u), %c, t::[%u, %u), qin->qs::%u, rin->qs::%u\n", __func__, p.qs, p.qe, "+-"[p.rev], p.ts, p.te, qin->qs, rin->qs); + } + } + } +} + +void cl_trans_gen(uint32_t qid, uint32_t qlen, ul_vec_t *qstr, kv_ul_ov_t *res, const ul_idx_t *uref, ma_ug_t *ref, scaf_res_t *ref_sc, ma_ug_t *gfa, bubble_type *bub, kv_u_trans_t *ta) +{ + uint64_t k, l, m = 0, mi, z, zn, a_n, bid, nw; uc_block_t *a; + res->n = 0; + ///identify bubble chain first + for (k = 0, a = qstr->bb.a; k < qstr->bb.n; k++) { + z = qstr->bb.a[k].ts; a_n = qstr->bb.a[k].te; + while (z < a_n) { + zn = 1; + zn = pick_bubble(a + z, a_n - z, bub, &bid); + if((!zn) || (!qry_bub_sc(uref, gfa, ref_sc, bub, bid, z, z + zn, qstr, res))) {///found a bubble; a[z] -> beg; a[z + zn] -> sink + zn = 1; + } + z += zn; + } + } + + if(res->n) {///matched bubbles + for (k = 0; k < res->n; k++) res->a[k].tn = (((uint64_t)res->a[k].tn)<<1)|((uint64_t)res->a[k].rev); + radix_sort_ul_ov_srt_tn(res->a, res->a + res->n); + + for (k = 1, l = m = 0; k <= res->n; k++) { + if(k == res->n || res->a[k].tn != res->a[l].tn) {///qn <- (tn|rev) + for (z = l; z < k; z++) res->a[z].tn>>=1; + a_n = k; + if(k - l > 1) {///if there is only one alignment for a node, then no need for DP + a_n = l + cc_bub_match(res->a+l, k-l); + } + + for (z = l; z < a_n; z++) res->a[m++] = res->a[z]; + // fprintf(stderr, "l::%lu, k::%lu, m::%lu\n", l, k, m); + l = k; + } + } + res->n = m; + radix_sort_ul_ov_srt_qe(res->a, res->a + res->n); + // for (k = 0; k < res->n; k++) { + // fprintf(stderr, "***qsidx::%u(utg%.6u%c), qeidx::%u(utg%.6u%c)\n", res->a[k].qs, (qstr->bb.a[res->a[k].qs].hid)+1, "lc"[gfa->u.a[qstr->bb.a[res->a[k].qs].hid].circ], res->a[k].qe, (qstr->bb.a[res->a[k].qe].hid)+1, "lc"[gfa->u.a[qstr->bb.a[res->a[k].qe].hid].circ]); + // } + + } + + for (k = 0; k < qstr->bb.n; k++) { + z = qstr->bb.a[k].ts; a_n = qstr->bb.a[k].te; + for (; z < a_n; z++) { + ///filter by bubble + for (mi = 0; mi < m && z <= res->a[mi].qe; mi++) { + if((z >= res->a[mi].qs) && (z <= res->a[mi].qe)) break; + } + if((mi < m) && (z >= res->a[mi].qs) && (z <= res->a[mi].qe)) continue; + + ///exact match + if(extract_scaf_res_t(&(qstr->bb.a[z]), uref, ref_sc, res)) continue; + + ///trans match + extract_trans_scaf_res_t(&(qstr->bb.a[z]), uref, ref_sc, gfa, ta, res); + } + } + + for(k = 0; k < m; k++) {///reset bubble + res->a[k].el = 1; res->a[k].qn = res->a[k].qs; + + res->a[k].qs = qstr->bb.a[res->a[k].qs].qs; + res->a[k].qe = qstr->bb.a[res->a[k].qe].qe; + + res->a[k].ts = ref_sc->a[res->a[k].tn].bb.a[res->a[k].ts].qs; + res->a[k].te = ref_sc->a[res->a[k].tn].bb.a[res->a[k].te].qe; + + nw = MAX(res->a[k].te-res->a[k].ts, res->a[k].qe-res->a[k].qs); nw *= CHAIN_MATCH; + res->a[k].sec = ((nw>=0x3fffffff)?(0x3fffffff):(nw)); + } + + // for (k = 0; k < res->n; k++) { + // fprintf(stderr, "[M::%s]\tutg%.6ul(len::%u)\tq::[%u, %u)\t%c\tutg%.6ul(len::%u)\tt::[%u, %u)\tel::%u\ttp::%u\n", __func__, qid + 1, qlen, res->a[k].qs, res->a[k].qe, "+-"[res->a[k].rev], res->a[k].tn+1, ref->u.a[res->a[k].tn].len, res->a[k].ts, res->a[k].te, res->a[k].el, res->a[k].qn); + // } +} + +void ctg_trans_gp_chain(kv_ul_ov_t *res, kv_ul_ov_t *buf, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, +double diff_ec_ul, int64_t qlen, Chain_Data* dp) +{ + uint64_t k/**, l, z, an, m**/; + for (k = 0; k < res->n; k++) { + res->a[k].tn <<= 1; res->a[k].tn |= res->a[k].rev; + } + radix_sort_ul_ov_srt_tn(res->a, res->a + res->n); + + +} + + +uint32_t direct_ctg_trans_chain(mg_tbuf_t *b, uint32_t id, glchain_t *ll, gdpchain_t *gdp, st_mt_t *sps, haplotype_evdience_alloc *hap, const ul_idx_t *uref, const ug_opt_t *uopt, +int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid, Chain_Data* dp, ma_ug_t *qry, scaf_res_t *qry_sc, ma_ug_t *ref, scaf_res_t *ref_sc, ma_ug_t *gfa, bubble_type *bub, kv_u_trans_t *ta, const asg_t *rg) +{ + // res->bb.n = 0; + // if(ulid != 86660) return 0; + kv_ul_ov_t *idx = &(ll->lo), *init = &(ll->tk); ///int64_t max_idx; + idx->n = init->n = 0; + cl_trans_gen(id, qry->u.a[id].len, &(qry_sc->a[id]), idx, uref, ref, ref_sc, gfa, bub, ta); + if(idx->n == 0) return 0; + + ctg_trans_gp_chain(idx, init, uref, uopt, bw, diff_ec_ul, qry->u.a[id].len, dp); + + if(idx->n == 0) return 0; + + + /** + // gl_rg2ug_gen(rch, idx, uref, 1, 2, ulid); + ctg_rg2ug_gen(rch, idx, uref); + + // uint32_t k; + // fprintf(stderr, "[M::%s] rch->len::%u, rch->n::%u, idx->n::%u\n", __func__, (uint32_t)rch->len, (uint32_t)rch->n, (uint32_t)idx->n); + // for (k = 0; k < idx->n; k++) { + // fprintf(stderr, "[k->%u::utg%.6u%c(len->%u::n->%u)]\tq::[%u, %u)\t%c\tt::[%u, %u)\n", + // k, (idx->a[k].tn>>1) + 1, "lc"[uref->ug->u.a[(idx->a[k].tn>>1)].circ], uref->ug->u.a[(idx->a[k].tn>>1)].len, idx->a[k].sec, + // idx->a[k].qs, idx->a[k].qe, "+-"[idx->a[k].rev], idx->a[k].ts, idx->a[k].te); + // } + + // for (k = 0; k < idx->n; k++) { + // fprintf(stderr, "utg%.6u%c,", (idx->a[k].tn>>1)+1, "lc"[uref->ug->u.a[(idx->a[k].tn>>1)].circ]); + // } + // fprintf(stderr, "\n"); + + if(idx->n == 0) return 0; + ///generate linear chains + gen_linear_chains_ctg(idx, init, uref, uopt, bw, diff_ec_ul, rch->len, dp); + + if(idx->n == 0) return 0; + + dump_linear_chain(uref->ug, idx, &(gdp->l), rch->len, ulid); + if(gdp->l.n == 0) return 0; + + kv_resize(uint64_t, ll->srt.a, gdp->l.n); + max_idx = ctg_chain_graph(b->km, uref, uref->ug, &(gdp->l), &(gdp->swap), &(gdp->dst), &(gdp->out), + &(gdp->path), rch->len, uopt, G_CHAIN_BW, N_GCHAIN_RATE, ll->srt.a.a, sps, dp, UG_SKIP_GRAPH_N, UG_ITER_N, UG_DIS_N*100); + + //sps -> idx; (gdp->l) -> alignments + if(max_idx >= 0) { + max_idx = select_max_ctg_chain(b->km, uref, ulid, sps, &(gdp->l), &(ll->tk), uref->ug->g, &(gdp->dst_done), &(gdp->out), &(gdp->path), &(gdp->swap)); + if(max_idx) { + push_ctg_res(res, &(gdp->swap)); + return 1; + } + } + **/ + + return 0; +} + + uint32_t refine_rid_chain(const asg_t *rg, mg_tbuf_t *b, ul_vec_t *rch, uint64_t ulid) { if(rch->bb.n == 1 && rch->bb.a[0].base) return 1;///no alignment @@ -15996,6 +16955,25 @@ static void worker_for_ul_gchains_alignment(void *data, long i, int tid) // gl_chain_refine_advance(&b->olist, &b->correct, &b->hap, bl, s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, s->id+i, km); } +static void worker_for_ctg_gchains_alignment(void *data, long i, int tid) +{ + utepdat_t *s = (utepdat_t*)data; ma_utg_t *q = &(s->ug->u.a[i]); s->rsc->a[i].bb.n = 0; + + s->hab[tid]->num_read_base++; + s->hab[tid]->num_correct_base += direct_gchain_scaf(s->buf[tid], q, &(s->ll[tid]), &(s->gdp[tid]), &(s->sps[tid]), &(s->hab[tid]->hap), s->uu, s->uopt, G_CHAIN_BW, s->opt->diff_ec_ul, UG_SKIP, i, &(s->hab[tid]->clist.chainDP), s->rg, &(s->rsc->a[i])); +} + +static void worker_for_ctg_trans_alignment(void *data, long i, int tid) +{ + ctdat_t *c = (ctdat_t*)data; + utepdat_t *s = c->s; + // ul_vec_t *q = &(c->qry_sc->a[i]); + + s->hab[tid]->num_read_base++; + s->hab[tid]->num_correct_base += direct_ctg_trans_chain(s->buf[tid], i, &(s->ll[tid]), &(s->gdp[tid]), &(s->sps[tid]), &(s->hab[tid]->hap), s->uu, s->uopt, G_CHAIN_BW, s->opt->diff_ec_ul, UG_SKIP, i, &(s->hab[tid]->clist.chainDP), + c->qry, c->qry_sc, c->ref, c->ref_sc, c->gfa, c->bub, c->ta, s->rg); +} + uint32_t inline is_rg_connect(const asg_t *rg, uint32_t v, uint32_t w) { uint32_t nv = asg_arc_n(rg, v), i; @@ -17244,6 +18222,67 @@ uint64_t work_ul_gchains(uldat_t *sl) return s.n; } +scaf_res_t *work_ctg_path_gchains(uldat_t *sl) +{ + utepdat_t s; uint64_t i; memset(&s, 0, sizeof(s)); scaf_res_t *res = init_scaf_res_t(sl->ug->u.n); + s.id = 0; s.opt = sl->opt; s.ug = sl->ug; s.uopt = sl->uopt; s.rg = sl->rg; s.uu = sl->uu; s.rsc = res; + CALLOC(s.hab, sl->n_thread); CALLOC(s.buf, sl->n_thread); CALLOC(s.ll, sl->n_thread); + CALLOC(s.gdp, sl->n_thread); CALLOC(s.mzs, sl->n_thread); CALLOC(s.sps, sl->n_thread); + + for (i = 0; i < sl->n_thread; ++i) { + s.hab[i] = ha_ovec_init(0, 0, 1); s.buf[i] = mg_tbuf_init(); + } + + // detect_outlier_len("+++work_ul_gchains"); + + kt_for(sl->n_thread, worker_for_ctg_gchains_alignment, &s, s.ug->u.n); + + // detect_outlier_len("---work_ul_gchains"); + + for (i = 0; i < sl->n_thread; ++i) { + s.sum_len += s.hab[i]->num_read_base; s.n += s.hab[i]->num_correct_base; + ha_ovec_destroy(s.hab[i]); mg_tbuf_destroy(s.buf[i]); hc_glchain_destroy(&(s.ll[i])); + hc_gdpchain_destroy(&(s.gdp[i])); kv_destroy(s.mzs[i]); kv_destroy(s.sps[i]); + } + + free(s.hab); free(s.buf); free(s.ll); free(s.gdp); free(s.mzs); free(s.sps); + fprintf(stderr, "[M::%s::] # try:%d, # done:%d\n", __func__, s.sum_len, s.n); + return res; +} + + +kv_u_trans_t *work_ctg_path_trans(uldat_t *sl, asg_t *sg, ma_ug_t *qry, scaf_res_t *qry_sc, ma_ug_t *ref, scaf_res_t *ref_sc, ma_ug_t *gfa, kv_u_trans_t *ta) +{ + kv_u_trans_t *res = NULL; CALLOC(res, 1); + utepdat_t s; uint64_t i; memset(&s, 0, sizeof(s)); + s.id = 0; s.opt = sl->opt; s.ug = sl->ug; s.uopt = sl->uopt; s.rg = sl->rg; s.uu = sl->uu; + CALLOC(s.hab, sl->n_thread); CALLOC(s.buf, sl->n_thread); CALLOC(s.ll, sl->n_thread); + CALLOC(s.gdp, sl->n_thread); CALLOC(s.mzs, sl->n_thread); CALLOC(s.sps, sl->n_thread); + + for (i = 0; i < sl->n_thread; ++i) { + s.hab[i] = ha_ovec_init(0, 0, 1); s.buf[i] = mg_tbuf_init(); + } + + ctdat_t c; memset(&c, 0, sizeof(c)); uint8_t *bf = NULL; c.bub = gen_bubble_chain(sg, gfa, (ug_opt_t *)(sl->uopt), &bf, 0); free(bf); + c.s = &(s); c.qry = qry; c.qry_sc = qry_sc; c.ref = ref; c.ref_sc = ref_sc; c.gfa = gfa; c.ta = ta; + + // detect_outlier_len("+++work_ul_gchains"); + + kt_for(sl->n_thread, worker_for_ctg_trans_alignment, &c, c.qry_sc->n); + + // detect_outlier_len("---work_ul_gchains"); + destory_bubbles(c.bub); free(c.bub); + for (i = 0; i < sl->n_thread; ++i) { + s.sum_len += s.hab[i]->num_read_base; s.n += s.hab[i]->num_correct_base; + ha_ovec_destroy(s.hab[i]); mg_tbuf_destroy(s.buf[i]); hc_glchain_destroy(&(s.ll[i])); + hc_gdpchain_destroy(&(s.gdp[i])); kv_destroy(s.mzs[i]); kv_destroy(s.sps[i]); + } + + free(s.hab); free(s.buf); free(s.ll); free(s.gdp); free(s.mzs); free(s.sps); + fprintf(stderr, "[M::%s::] # try:%d, # done:%d\n", __func__, s.sum_len, s.n); + return res; +} + uint64_t work_ul_gchains_consensus(uldat_t *sl) { @@ -19510,20 +20549,24 @@ ul_idx_t *gen_ul_idx(const ug_opt_t *uopt, ma_ug_t *ug, asg_t *sg) -utg_rid_t *gen_r_ug_idx(ma_ug_t *ug, ul_contain *ct, asg_t *rg) +utg_rid_t *gen_r_ug_idx(ma_ug_t *ug, ul_contain *ct, asg_t *rg, uint32_t keep_offset) { - uint64_t i, k, l, m, rid, a_n, cn; utg_rid_dt *a; ma_utg_t *u = NULL; utg_ct_t *ca; - utg_rid_t *cc = NULL; CALLOC(cc, 1); CALLOC(cc->idx, rg->n_seq+1); kv_init(cc->p); cc->rg = rg; + uint64_t i, k, l, m, rid, a_n, cn, rn = R_INF.total_reads; ///R_INF.total_reads might be larger than rg->n_seq due scaffold reads + utg_rid_dt *a; ma_utg_t *u = NULL; utg_ct_t *ca; + + utg_rid_t *cc = NULL; CALLOC(cc, 1); CALLOC(cc->idx, rn+1); kv_init(cc->p); cc->rg = rg; for (i = 0; i < ug->u.n; i++) { u = &(ug->u.a[i]); for (k = 0; k < u->n; k++) cc->idx[u->a[k]>>33]++; - cn = ((uint32_t)(ct->idx.a[i])); - ca = ct->rids.a + ((ct->idx.a[i])>>32); - for (k = 0; k < cn; k++) cc->idx[ca[k].x>>1]++; + if(ct) { + cn = ((uint32_t)(ct->idx.a[i])); + ca = ct->rids.a + ((ct->idx.a[i])>>32); + for (k = 0; k < cn; k++) cc->idx[ca[k].x>>1]++; + } } - for (k = l = 0; k <= rg->n_seq; k++) { + for (k = l = 0; k <= rn; k++) { m = cc->idx[k]; cc->idx[k] = l; l += m; @@ -19538,30 +20581,32 @@ utg_rid_t *gen_r_ug_idx(ma_ug_t *ug, ul_contain *ct, asg_t *rg) if(a_n) { if(a[a_n-1].off == a_n-1) { a[a_n-1].u = (i<<1)|((u->a[k]>>32)&1); - a[a_n-1].pos = l; a[a_n-1].off = l; + a[a_n-1].pos = l; a[a_n-1].off = (keep_offset?k:l); } else { a[a[a_n-1].off].u = (i<<1)|((u->a[k]>>32)&1); - a[a[a_n-1].off].pos = l; a[a[a_n-1].off].off = l; + a[a[a_n-1].off].pos = l; a[a[a_n-1].off].off = (keep_offset?k:l); a[a_n-1].off++; } } l += (uint32_t)u->a[k]; } - cn = ((uint32_t)(ct->idx.a[i])); - ca = ct->rids.a + ((ct->idx.a[i])>>32); - for (k = 0; k < cn; k++) { - rid = ca[k].x>>1; - a = cc->p.a + cc->idx[rid]; - a_n = cc->idx[rid+1] - cc->idx[rid]; - if(a_n) { - if(a[a_n-1].off == a_n-1) { - a[a_n-1].u = (i<<1)|(ca[k].x&1); - a[a_n-1].pos = ca[k].s; a[a_n-1].off = ca[k].s; - } else { - a[a[a_n-1].off].u = (i<<1)|(ca[k].x&1); - a[a[a_n-1].off].pos = ca[k].s; a[a[a_n-1].off].off = ca[k].s; - a[a_n-1].off++; + if(ct) { + cn = ((uint32_t)(ct->idx.a[i])); + ca = ct->rids.a + ((ct->idx.a[i])>>32); + for (k = 0; k < cn; k++) { + rid = ca[k].x>>1; + a = cc->p.a + cc->idx[rid]; + a_n = cc->idx[rid+1] - cc->idx[rid]; + if(a_n) { + if(a[a_n-1].off == a_n-1) { + a[a_n-1].u = (i<<1)|(ca[k].x&1); + a[a_n-1].pos = ca[k].s; a[a_n-1].off = (keep_offset?k:ca[k].s); + } else { + a[a[a_n-1].off].u = (i<<1)|(ca[k].x&1); + a[a[a_n-1].off].pos = ca[k].s; a[a[a_n-1].off].off = (keep_offset?k:ca[k].s); + a[a_n-1].off++; + } } } } @@ -19569,18 +20614,72 @@ utg_rid_t *gen_r_ug_idx(ma_ug_t *ug, ul_contain *ct, asg_t *rg) return cc; } -ul_idx_t *gen_ul_idx_t(const ug_opt_t *uopt, asg_t *sg, uint64_t is_el, uint64_t is_del) +utg_rid_t *gen_r_ug_sc_idx(scaf_res_t *ug_sc, uint64_t ug_sc_num) +{ + uint64_t i, z, k, l, m, rid, a_n, s, e, rn = ug_sc_num; ///ug_sc_num is the number of nodes within the inital graph + utg_rid_dt *a; ul_vec_t *idx; + + utg_rid_t *cc = NULL; CALLOC(cc, 1); CALLOC(cc->idx, rn+1); kv_init(cc->p); ///cc->rg = rg; + for (i = 0; i < ug_sc->n; i++) { + idx = &(ug_sc->a[i]); + for (z = 0; z < idx->bb.n; z++) { + s = idx->bb.a[z].ts; e = idx->bb.a[z].te; + for (k = s; k < e; k++) cc->idx[idx->bb.a[k].hid]++; + } + } + + for (k = l = 0; k <= rn; k++) { + m = cc->idx[k]; + cc->idx[k] = l; + l += m; + } + cc->p.n = cc->p.m = l; CALLOC(cc->p.a, cc->p.n); + for (i = 0; i < ug_sc->n; i++) { + idx = &(ug_sc->a[i]); + for (z = 0; z < idx->bb.n; z++) { + s = idx->bb.a[z].ts; e = idx->bb.a[z].te; + for (k = s; k < e; k++) { + rid = idx->bb.a[k].hid; + a = cc->p.a + cc->idx[rid]; + a_n = cc->idx[rid+1] - cc->idx[rid]; + if(a_n) { + if(a[a_n-1].off == a_n-1) { + a[a_n-1].u = (i<<1)|(idx->bb.a[k].rev); + a[a_n-1].pos = idx->bb.a[k].qs; a[a_n-1].off = k; + } else { + a[a[a_n-1].off].u = (i<<1)|(idx->bb.a[k].rev); + a[a[a_n-1].off].pos = idx->bb.a[k].qs; a[a[a_n-1].off].off = k; + a[a_n-1].off++; + } + } + } + } + } + return cc; +} + +ul_idx_t *gen_ul_idx_t(const ug_opt_t *uopt, ma_ug_t *ug, asg_t *sg, uint64_t is_el, uint64_t is_del, uint64_t is_scaf) { uint64_t n_read = R_INF.total_reads; ma_hit_t_alloc* src = uopt->sources; int64_t min_ovlp = uopt->min_ovlp; int64_t max_hang = uopt->max_hang; // int64_t gap_fuzz = uopt->gap_fuzz; - ul_idx_t *uu = NULL; CALLOC(uu, 1); uu->ug = ma_ug_gen(sg); - uu->ct = ul_contain_gen(uu->ug, sg, src, min_ovlp, max_hang, is_el, is_del); - uu->cc = gen_cov_track(uu->ug, sg, uu->ct, src, min_ovlp, max_hang, is_el, is_del); - uu->cr = gen_r_contain(uu->ug, sg, src, n_read, min_ovlp, max_hang, asm_opt.thread_num, is_el, is_del); - uu->r_ug = gen_r_ug_idx(uu->ug, uu->ct, sg); + ul_idx_t *uu = NULL; CALLOC(uu, 1); uu->ug = (ug?(ug):(ma_ug_gen(sg))); + if(!is_scaf) { + uu->ct = ul_contain_gen(uu->ug, sg, src, min_ovlp, max_hang, is_el, is_del); + uu->cc = gen_cov_track(uu->ug, sg, uu->ct, src, min_ovlp, max_hang, is_el, is_del); + uu->cr = gen_r_contain(uu->ug, sg, src, n_read, min_ovlp, max_hang, asm_opt.thread_num, is_el, is_del); + } + uu->r_ug = gen_r_ug_idx(uu->ug, (!is_scaf)?(uu->ct):(NULL), sg, is_scaf); + return uu; +} + +ul_idx_t *gen_ul_idx_t_sc(ma_ug_t *ug, asg_t *sg, scaf_res_t *ug_sc, uint64_t ug_sc_num) +{ + // int64_t gap_fuzz = uopt->gap_fuzz; + ul_idx_t *uu = NULL; CALLOC(uu, 1); uu->ug = (ug?(ug):(ma_ug_gen(sg))); + uu->r_ug = gen_r_ug_sc_idx(ug_sc, ug_sc_num); return uu; } @@ -20071,7 +21170,7 @@ uint64_t ul_refine_alignment(const ug_opt_t *uopt, asg_t *sg) init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov); cutoff = asm_opt.max_n_chain; init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, asm_opt.ul_error_rate, asm_opt.ul_error_rate_low, asm_opt.ul_error_rate_hpc, asm_opt.ul_ec_round); - ul_idx_t *uu = gen_ul_idx_t(uopt, sg, 0, 0);///record contained reads; is_el = is_del = 0 + ul_idx_t *uu = gen_ul_idx_t(uopt, NULL, sg, 0, 0, 0);///record contained reads; is_el = is_del = 0 init_uldat_t(&sl, NULL, NULL, &opt, CHUNK_SIZE, asm_opt.thread_num, uopt, uu); sl.rg = sg; if(work_ul_gchains(&sl)) { free(UL_INF.ridx.idx.a); free(UL_INF.ridx.occ.a); memset(&(UL_INF.ridx), 0, sizeof(UL_INF.ridx)); @@ -20086,6 +21185,35 @@ uint64_t ul_refine_alignment(const ug_opt_t *uopt, asg_t *sg) } } +scaf_res_t *gen_contig_path(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *ctg, ma_ug_t *ref) +{ + mg_idxopt_t opt; uldat_t sl; int32_t cutoff; + init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov); + cutoff = asm_opt.max_n_chain; + init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, asm_opt.ul_error_rate, asm_opt.ul_error_rate_low, asm_opt.ul_error_rate_hpc, asm_opt.ul_ec_round); + ul_idx_t *uu = gen_ul_idx_t(uopt, ref, sg, 0, 0, 1); + init_uldat_t(&sl, NULL, NULL, &opt, CHUNK_SIZE, asm_opt.thread_num, uopt, uu); sl.rg = sg; sl.ug = ctg; + + scaf_res_t *res = work_ctg_path_gchains(&sl); uu->ug = NULL; + destroy_ul_idx_t(uu); + return res; +} + +kv_u_trans_t *gen_contig_trans(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *qry, scaf_res_t *qry_sc, ma_ug_t *ref, scaf_res_t *ref_sc, ma_ug_t *gfa, kv_u_trans_t *ta) +{ + mg_idxopt_t opt; uldat_t sl; int32_t cutoff; + init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov); + cutoff = asm_opt.max_n_chain; + init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, asm_opt.ul_error_rate, asm_opt.ul_error_rate, asm_opt.ul_error_rate_low, asm_opt.ul_error_rate_hpc, asm_opt.ul_ec_round); + ul_idx_t *uu = gen_ul_idx_t_sc(ref, sg, ref_sc, gfa->u.n); + init_uldat_t(&sl, NULL, NULL, &opt, CHUNK_SIZE, asm_opt.thread_num, uopt, uu); sl.rg = sg; sl.ug = qry; + + kv_u_trans_t *res = work_ctg_path_trans(&sl, sg, qry, qry_sc, ref, ref_sc, gfa, ta); uu->ug = NULL; + destroy_ul_idx_t(uu); + return res; +} + + uint32_t dd_ug(asg_t *sg, ma_ug_t *ug, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, const char* output_file_name) { fprintf(stderr, "Writing raw unitig GFA to disk... \n"); diff --git a/inter.h b/inter.h index 648157a..e7e180e 100644 --- a/inter.h +++ b/inter.h @@ -127,5 +127,7 @@ uint32_t rqs, uint32_t rqe, uint32_t *rts, uint32_t *rte); uint32_t clean_contain_g(const ug_opt_t *uopt, asg_t *sg, uint32_t push_trans); void dedup_contain_g(const ug_opt_t *uopt, asg_t *sg); void trans_base_mmhap_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res); +scaf_res_t *gen_contig_path(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *ctg, ma_ug_t *ref); +kv_u_trans_t *gen_contig_trans(const ug_opt_t *uopt, asg_t *sg, ma_ug_t *qry, scaf_res_t *qry_sc, ma_ug_t *ref, scaf_res_t *ref_sc, ma_ug_t *gfa, kv_u_trans_t *ta); #endif