diff --git a/Overlaps.cpp b/Overlaps.cpp index ca1c72d..46f1fa7 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -8759,69 +8759,6 @@ static char comp_tab[] = { // complement base 'p', 'q', 'y', 's', 'a', 'a', 'b', 'w', 'x', 'r', 'z', 123, 124, 125, 126, 127 }; -// generate unitig sequences -int ma_ug_seq_back(ma_ug_t *g, All_reads *RNF, const ma_sub_t *coverage_cut, -const long long n_read) -{ - UC_Read g_read; - init_UC_Read(&g_read); - utg_intv_t *tmp; - uint32_t i, j; - - - ///why we need n_read here? it is just beacuse one read can only be included in one untig - ///but it is not true - tmp = (utg_intv_t*)calloc(n_read, sizeof(utg_intv_t)); - ///number of unitigs - for (i = 0; i < g->u.n; ++i) { - ma_utg_t *u = &g->u.a[i]; - uint32_t l = 0; - u->s = (char*)calloc(1, u->len + 1); - memset(u->s, 'N', u->len); - for (j = 0; j < u->n; ++j) { - ///u->a[j]>>33 is the readID - utg_intv_t *t = &tmp[u->a[j]>>33]; - ///assert(t->len == 0); - t->utg = i, t->ori = u->a[j]>>32&1; - ///l is the start pos of this read at its corresponding untig - t->start = l, t->len = (uint32_t)u->a[j]; - l += t->len; - } - } - - - - int32_t id; - for (id = 0; id < n_read; id++) - { - utg_intv_t *t; - ma_utg_t *u; - if (id < 0 || tmp[id].len == 0) continue; - - t = &tmp[id]; - u = &g->u.a[t->utg]; - recover_UC_Read(&g_read, RNF, id); - - memmove(g_read.seq, g_read.seq + coverage_cut[id].s, coverage_cut[id].e - coverage_cut[id].s); - g_read.length = coverage_cut[id].e - coverage_cut[id].s; - - if (!t->ori) { // forward strand - for (i = 0; i < t->len; ++i) - u->s[t->start + i] = g_read.seq[i]; - } else { - for (i = 0; i < t->len; ++i) { - int c = (uint8_t)g_read.seq[g_read.length - 1 - i]; - u->s[t->start + i] = c >= 128? 'N' : comp_tab[c]; - } - } - } - - - free(tmp); - - destory_UC_Read(&g_read); - return 0; -} void recover_fake_read(UC_Read* result, UC_Read* tmp, ma_utg_t *u, All_reads *RNF, const ma_sub_t *coverage_cut) @@ -9127,6 +9064,14 @@ ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp) if(pre == (uint32_t)-1) continue; if(afte == (uint32_t)-1) continue; + /****************************may have bugs********************************/ + ///if(debug_purge_dup == 1 && ((v>>1)==1239235 || (v>>1)==4576917 || (v>>1)==4479766)) + if(debug_purge_dup == 1 && (v>>1)==4576917) + { + fprintf(stderr, "v>>1: %u, pre>>1: %u, afte>>1: %u, len: %u\n", v>>1, pre>>1, afte>>1, collection->len); + } + /****************************may have bugs********************************/ + av = asg_arc_a(read_g, v^1); nv = asg_arc_n(read_g, v^1); for (k = 0; k < nv; k++) @@ -9177,8 +9122,15 @@ ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp) } } if(k == edge->a.n) fprintf(stderr, "ERROR\n"); - } + } + /****************************may have bugs********************************/ + ///if(debug_purge_dup == 1 && ((v>>1)==1239235 || (v>>1)==4576917 || (v>>1)==4479766)) + if(debug_purge_dup == 1 && (v>>1)==4576917) + { + fprintf(stderr, "v>>1: %u, pE->el: %u, aE->el: %u\n", v>>1, pE->el, aE->el); + } + /****************************may have bugs********************************/ if(pE->el == 1 && aE->el == 1) continue; @@ -9195,6 +9147,304 @@ ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp) continue; } + + /****************************may have bugs********************************/ + ///if(debug_purge_dup == 1 && ((v>>1)==1239235 || (v>>1)==4576917 || (v>>1)==4479766)) + if(debug_purge_dup == 1 && (v>>1)==4576917) + { + fprintf(stderr, "v>>1: %u, t_f.el: %u, t_b.el: %u\n", v>>1, t_f.el, t_b.el); + } + /****************************may have bugs********************************/ + + if(t_f.el == 0 || t_b.el == 0) continue; + + get_overlapLen(v>>1, sources, &exactLen, &inexactLen); + if(pE->el == 0) + { + get_overlapLen(pre>>1, sources, &tmp_exactLen, &tmp_inexactLen); + + + /****************************may have bugs********************************/ + ///if(debug_purge_dup == 1 && ((v>>1)==1239235 || (v>>1)==4576917 || (v>>1)==4479766)) + if(debug_purge_dup == 1 && (v>>1)==4576917) + { + fprintf(stderr, "v>>1: %u, exactLen: %u, inexactLen: %u, tmp_exactLen: %u, tmp_inexactLen: %u\n", + v>>1, exactLen, inexactLen, tmp_exactLen, tmp_inexactLen); + } + /****************************may have bugs********************************/ + if(inexactLen < tmp_inexactLen) continue; + if(inexactLen == tmp_inexactLen && exactLen > tmp_exactLen) continue; + } + + if(aE->el == 0) + { + get_overlapLen(afte>>1, sources, &tmp_exactLen, &tmp_inexactLen); + if(inexactLen < tmp_inexactLen) continue; + if(inexactLen == tmp_inexactLen && exactLen > tmp_exactLen) continue; + } + + + /****************************may have bugs********************************/ + ///if(debug_purge_dup == 1 && ((v>>1)==1239235 || (v>>1)==4576917 || (v>>1)==4479766)) + if(debug_purge_dup == 1 && (v>>1)==4576917) + { + fprintf(stderr, "v>>1: %u, i: %u\n", v>>1, i); + } + /****************************may have bugs********************************/ + + collection->a[i] = (uint64_t)-1; + skip++; + } + + if(skip == 0) return 0; + + m = 0; + for (i = 0; i < collection->n; i++) + { + if(collection->a[i] == (uint64_t)-1) continue; + collection->a[m] = collection->a[i]; + m++; + } + collection->n = m; + + + uint32_t totalLen = 0, w, l; + for (i = 0; i < collection->n - 1; i++) + { + v = (uint64_t)(collection->a[i])>>32; + w = (uint64_t)(collection->a[i + 1])>>32; + + + + + /*******************************for debug************************************/ + l = (uint32_t)-1; + av = asg_arc_a(read_g, v); + nv = asg_arc_n(read_g, 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(sources, coverage_cut, 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, read_g->r_seq); + } + l = asg_arc_len(t_f); + } + } + if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n"); + /*******************************for debug************************************/ + + + + + + + + + + collection->a[i] = v; collection->a[i] = collection->a[i]<<32; + collection->a[i] = collection->a[i] | (uint64_t)(l); + totalLen += l; + } + + if(i < collection->n) + { + if(collection->circ) + { + v = (uint64_t)(collection->a[i])>>32; + w = (uint64_t)(collection->a[0])>>32; + + + + + /*******************************for debug************************************/ + l = (uint32_t)-1; + av = asg_arc_a(read_g, v); + nv = asg_arc_n(read_g, 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(sources, coverage_cut, 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, read_g->r_seq); + } + l = asg_arc_len(t_f); + } + } + if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n"); + /*******************************for debug************************************/ + + + + + + + + + + + + + collection->a[i] = v; collection->a[i] = collection->a[i]<<32; + collection->a[i] = collection->a[i] | (uint64_t)(l); + + totalLen += l; + } + else + { + v = (uint64_t)(collection->a[i])>>32; + l = read_g->seq[v>>1].len; + collection->a[i] = v; + collection->a[i] = collection->a[i]<<32; + collection->a[i] = collection->a[i] | (uint64_t)(l); + totalLen += l; + } + } + + collection->len = totalLen; + if(!collection->circ) + { + collection->start = collection->a[0]>>32; + collection->end = (collection->a[collection->n-1]>>32)^1; + } + + return 1; +} + + + + +uint32_t polish_unitig_advance(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) +{ + if(collection->m == 0) return 0; + if(collection->n < 2) return 0; + uint32_t i, k, v, pre, afte, nv, exactLen, inexactLen, tmp_exactLen, tmp_inexactLen, m = 0, skip = 0; + asg_arc_t* av = NULL; + asg_arc_t *pE = NULL, *aE = NULL; + asg_arc_t t_f, t_b; + + for (i = 1; i < collection->n; i++) + { + v = (uint64_t)(collection->a[i])>>32; + pre = (uint64_t)(collection->a[i-1])>>32; + afte = (uint64_t)(collection->a[i+1])>>32; + if(v == (uint32_t)-1) continue; + if(pre == (uint32_t)-1) continue; + if(afte == (uint32_t)-1) continue; + + av = asg_arc_a(read_g, v^1); + nv = asg_arc_n(read_g, v^1); + for (k = 0; k < nv; k++) + { + if(av[k].del) continue; + if(av[k].v == (pre^1)) + { + pE = &(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^1) && edge->a.a[k].v == (pre^1)) + { + pE = &(edge->a.a[k]); + break; + } + } + + if(k == edge->a.n) fprintf(stderr, "ERROR\n"); + } + + av = asg_arc_a(read_g, v); + nv = asg_arc_n(read_g, v); + for (k = 0; k < nv; k++) + { + if(av[k].del) continue; + if(av[k].v == afte) + { + aE = &(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 == afte) + { + aE = &(edge->a.a[k]); + break; + } + } + if(k == edge->a.n) fprintf(stderr, "ERROR\n"); + } + + + if(pE->el == 1 && aE->el == 1) continue; + + if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, pre, + afte, &t_f) == 0) + { + continue; + } + + if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, afte^1, + pre^1, &t_b) == 0) + { + continue; + } + + if(t_f.el == 0 || t_b.el == 0) continue; get_overlapLen(v>>1, sources, &exactLen, &inexactLen); @@ -9212,6 +9462,7 @@ ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp) if(inexactLen == tmp_inexactLen && exactLen > tmp_exactLen) continue; } + collection->a[i] = (uint64_t)-1; skip++; } @@ -9375,76 +9626,6 @@ ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp) return 1; } // generate unitig sequences -int ma_ug_seq_back(ma_ug_t *g, asg_t *read_g, All_reads *RNF, const ma_sub_t *coverage_cut) -{ - UC_Read g_read; - init_UC_Read(&g_read); - UC_Read tmp; - init_UC_Read(&tmp); - ///utg_intv_t *tmp; - uint32_t i, j, k; - uint32_t rId, /**uId,**/ori, start, eLen, readLen; - char* readS = NULL; - - - ///why we need n_read here? it is just beacuse one read can only be included in one untig - ///but it is not true - ///tmp = (utg_intv_t*)calloc(n_read, sizeof(utg_intv_t)); - ///number of unitigs - for (i = 0; i < g->u.n; ++i) { - ma_utg_t *u = &g->u.a[i]; - if(u->m == 0) continue; - - uint32_t l = 0; - u->s = (char*)calloc(1, u->len + 1); - memset(u->s, 'N', u->len); - for (j = 0; j < u->n; ++j) { - rId = u->a[j]>>33; - ///uId = i; - ori = u->a[j]>>32&1; - start = l; - eLen = (uint32_t)u->a[j]; - l += eLen; - - if(eLen == 0) continue; - if(rId < read_g->r_seq) - { - recover_UC_Read(&g_read, RNF, rId); - } - else - { - recover_fake_read(&g_read, &tmp, &(read_g->F_seq[rId-read_g->r_seq]), - RNF, coverage_cut); - } - - 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++) - { - u->s[start + k] = readS[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]; - } - } - - - } - } - - destory_UC_Read(&g_read); - destory_UC_Read(&tmp); - return 0; -} - @@ -9553,8 +9734,8 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag) if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue; } if(read_g->seq[tn].del == 1) continue; - if(r_flag[tn] != 1) continue; C_bases += (Get_qe((*h)) - Get_qs((*h))); + if(r_flag[tn] != 1) continue; } } @@ -23015,8 +23196,13 @@ R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ov adjust_utg_by_primary(&ug, sg, TRIO_THRES, sources, reverse_sources, coverage_cut, bubble_dist, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, max_hang, min_ovlp, &new_rtg_edges); - + /****************************may have bugs********************************/ + debug_purge_dup = 1; + /****************************may have bugs********************************/ ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp); + /****************************may have bugs********************************/ + debug_purge_dup = 0; + /****************************may have bugs********************************/ fprintf(stderr, "Writing primary contig GFA to disk... \n"); @@ -27168,10 +27354,9 @@ long long bubble_dist, int read_graph, int write) &R_INF, output_file_name); } - debug_info_of_specfic_read("m64062_190803_042216/120916077/ccs", sources, reverse_sources, -1, "beg"); - debug_info_of_specfic_read("m64062_190806_063919/71436446/ccs", sources, reverse_sources, -1, "beg"); - debug_info_of_specfic_read("m64062_190803_042216/82117654/ccs", sources, reverse_sources, -1, "beg"); - + debug_info_of_specfic_read("m64062_190803_042216/177341795/ccs", sources, reverse_sources, -1, "beg"); + // debug_info_of_specfic_read("m64062_190807_194840/126682874/ccs", sources, reverse_sources, -1, "beg"); + exit(1); if (!(asm_opt.flag & HA_F_BAN_ASSEMBLY)) {