diff --git a/Assembly.cpp b/Assembly.cpp index a87345c..75d706f 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -299,30 +299,6 @@ void push_overlaps(ma_hit_t_alloc* paf, overlap_region_alloc* overlap_list, int } } -int if_exact_match(char* x, long long xLen, char* y, long long yLen, long long xBeg, long long xEnd, long long yBeg, long long yEnd) -{ - long long overlapLen = xEnd - xBeg + 1; - - if(yEnd - yBeg + 1 == overlapLen) - { - long long i; - - for (i = 0; i < overlapLen; i++) - { - if(x[xBeg + i] != y[yBeg + i]) - { - break; - } - } - - if(i == overlapLen) - { - return 1; - } - } - - return 0; -} long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_list, overlap_region_alloc* overlap_list, int flag) { diff --git a/Correct.h b/Correct.h index f51a5df..452cd29 100644 --- a/Correct.h +++ b/Correct.h @@ -1175,6 +1175,20 @@ long long get_affine_gap_score(overlap_region* ovc, UC_Read* g_read, UC_Read* ov uint8_t* y_num, uint64_t EstimateXOlen, uint64_t EstimateYOlen); int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u8_warp* k_flag, float hp_rate, int rlen, FILE* fp); +inline int if_exact_match(char* x, long long xLen, char* y, long long yLen, long long xBeg, long long xEnd, long long yBeg, long long yEnd) +{ + long long overlapLen = xEnd - xBeg + 1; + + if(yEnd - yBeg + 1 == overlapLen) + { + if(memcmp(x+xBeg, y+yBeg, overlapLen)==0) return 1; + } + + return 0; +} + + + #define FORWARD_KSW 0 #define BACKWARD_KSW 1 #define MATCH_SCORE_KSW 2 diff --git a/Overlaps.cpp b/Overlaps.cpp index 46f1fa7..7cbf948 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -9045,12 +9045,172 @@ ma_sub_t *coverage_cut, int max_hang, int min_ovlp) return 1; } +void reduce_ma_utg_t(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) +{ + asg_arc_t t_f; + asg_arc_t* av = NULL; + uint32_t m = 0, k, i, nv; + 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, v, 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; + } +} + uint32_t polish_unitig(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 < 3) return 0; - uint32_t i, k, v, pre, afte, nv, exactLen, inexactLen, tmp_exactLen, tmp_inexactLen, m = 0, skip = 0; + uint32_t i, k, v, pre, afte, nv, exactLen, inexactLen, tmp_exactLen, tmp_inexactLen, skip = 0; asg_arc_t* av = NULL; asg_arc_t *pE = NULL, *aE = NULL; asg_arc_t t_f, t_b; @@ -9064,14 +9224,6 @@ 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++) @@ -9124,15 +9276,6 @@ 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; if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, pre, @@ -9147,31 +9290,12 @@ 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; } @@ -9183,450 +9307,179 @@ 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; } - - /****************************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; - } + reduce_ma_utg_t(collection, read_g, sources, coverage_cut, edge, max_hang, min_ovlp); return 1; } +int get_consensus_rate(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, +int* r_match, int* r_total) +{ + (*r_match) = (*r_total) = 0; + if(cur_i < 1) return -1; + asg_arc_t *t = NULL, i_t; + uint32_t v, w, k, v_beg, v_end, w_beg, w_end; + int j; + if(collection->a[cur_i] == (uint64_t)-1 || collection->a[next_i] == (uint64_t)-1) return 0; + + v = (uint64_t)(collection->a[cur_i])>>32; + w = (uint64_t)(collection->a[next_i])>>32; + + if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, v, w, &i_t) == 0) + { + 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) + { + t = &(edge->a.a[k]); + break; + } + } + } + else + { + t = &i_t; + } + + if(t == NULL) return -1; + v_beg = 0; v_end = asg_arc_len(*t) - 1; r_read->length = 0; + ///cur_i must >= 1 + for (j = cur_i-1; j >= 0; j--) + { + if(collection->a[j] == (uint64_t)-1) continue; + + w = (uint64_t)(collection->a[j])>>32; + t = NULL; + + if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, w, v, &i_t) == 0) + { + for (k = 0; k < edge->a.n; k++) + { + if(edge->a.a[k].del) continue; + if((edge->a.a[k].ul>>32) == w && edge->a.a[k].v == v) + { + t = &(edge->a.a[k]); + break; + } + } + } + else + { + t = &i_t; + } + + if(t == NULL) break; + + w_beg = asg_arc_len(*t); + w_end = MIN(((int)(w_beg + v_end)), ((int)(read_g->seq[w>>1].len-1))); + // if(t->el == 1) + // { + // if(r_read->length == 0) recover_UC_sub_Read(r_read, v_beg, v_end - v_beg + 1, v&1, RNF, v>>1); + // ///for debug: check if v[v_beg, v_end] == w[w_beg, w_end] + // recover_UC_sub_Read(q_read, w_beg, w_end - w_beg +1, w&1, RNF, w>>1); + // ///q_read->length<=r_read->length + // if(if_exact_match(r_read->seq, r_read->length, q_read->seq, q_read->length, + // 0, q_read->length-1, 0, q_read->length-1) == 0) + // { + // fprintf(stderr, "ERROR: cur_i: %u, j: %d, r_len: %lld, q_len: %lld\n", cur_i, j, r_read->length, q_read->length); + // } + // } + + //filter overlaps which cannot cover the whole [v_beg, v_end] + if(w_end - w_beg != v_end - v_beg) break; + + (*r_total)++; + + if(t->el == 1) + { + (*r_match)++; + } + else + { + if(r_read->length == 0) recover_UC_sub_Read(r_read, v_beg, v_end - v_beg + 1, v&1, RNF, v>>1); + recover_UC_sub_Read(q_read, w_beg, w_end - w_beg +1, w&1, RNF, w>>1); + ///q_read->length<=r_read->length + if(if_exact_match(r_read->seq, r_read->length, q_read->seq, q_read->length, + 0, q_read->length-1, 0, q_read->length-1) == 1) + { + (*r_match)++; + } + } + } -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) + return 1; +} + +uint32_t polish_unitig_advance(ma_utg_t* collection, 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) { 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(collection->n < 3) return 0; + uint32_t i, skip = 0; + int match_v, total_v, max_i, match_max, k; + double match_rate, match_rate_max; - if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, pre, - afte, &t_f) == 0) + ///we should be able to handle i = collection->n-1 + for (i = 1; i < collection->n-1; i++) + { + ///in practice, collection->a[i] and collection->a[i+1] must be available + ///collection->a[index] might be unavailable only if index < i + if(get_consensus_rate(collection, i, i+1, read_g, RNF, sources, coverage_cut, + 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 collection[i], so it is right + if(match_v >= total_v * 0.5 && total_v > 0 && match_v > 0) continue; - if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, afte^1, - pre^1, &t_b) == 0) + max_i = i; match_max = match_v; match_rate_max = match_rate; + for (k = i - 1; k >= 0; k--) { - continue; - } - - - 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); - 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; - } - - - 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) + if(collection->a[k] == (uint64_t)-1) continue; + ///collection->a[k] might be unavailable, while collection->a[i+1] must be available + ///return -1 means there is no overlap from k to i+1 + if(get_consensus_rate(collection, k, i+1, read_g, RNF, sources, coverage_cut, + edge, r_read, q_read, max_hang, min_ovlp, &match_v, &total_v) < 0) { - l = asg_arc_len(av[k]); break; } - } + ///no read support k to i+1 + if(total_v == 0) break; - if(k == nv) - { - for (k = 0; k < edge->a.n; k++) + 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)) { - 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); + max_i = k; match_max = match_v; match_rate_max = match_rate; } } - 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) + + ///set [max_i+1, i] to be unavailable + for (k = max_i+1; k <= (int)i; k++) { - 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; + if(collection->a[k] == (uint64_t)-1) continue; + collection->a[k] = (uint64_t)-1; + skip++; } } - collection->len = totalLen; - if(!collection->circ) - { - collection->start = collection->a[0]>>32; - collection->end = (collection->a[collection->n-1]>>32)^1; - } + if(skip == 0) return 1; + + reduce_ma_utg_t(collection, read_g, sources, coverage_cut, edge, max_hang, min_ovlp); return 1; } -// generate unitig sequences - // generate unitig sequences @@ -9653,6 +9506,8 @@ ma_hit_t_alloc* sources, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp) polish_unitig(u, read_g, sources, coverage_cut, edge, max_hang, min_ovlp); + polish_unitig_advance(u, read_g, RNF, sources, coverage_cut, edge, &g_read, &tmp, max_hang, min_ovlp); + uint32_t l = 0; u->s = (char*)calloc(1, u->len + 1); memset(u->s, 'N', u->len); @@ -9734,8 +9589,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; - C_bases += (Get_qe((*h)) - Get_qs((*h))); if(r_flag[tn] != 1) continue; + C_bases += (Get_qe((*h)) - Get_qs((*h))); } } @@ -23196,13 +23051,7 @@ 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"); @@ -27354,9 +27203,8 @@ long long bubble_dist, int read_graph, int write) &R_INF, output_file_name); } - debug_info_of_specfic_read("m64062_190803_042216/177341795/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)) { diff --git a/Process_Read.cpp b/Process_Read.cpp index cb23a98..16855a7 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -338,93 +338,6 @@ void init_UC_Read(UC_Read* r) } } -void recover_UC_Read_sub_region_begin_end(char* r, long long start_pos, long long length, uint8_t strand, - All_reads* R_INF, long long ID, int extra_begin, int extra_end) -{ - long long readLen = Get_READ_LENGTH((*R_INF), ID); - uint8_t* src = Get_READ((*R_INF), ID); - - long long i; - long long copyLen; - long long end_pos = start_pos + length - 1; - - if (strand == 0) - { - i = start_pos; - copyLen = 0; - - long long initLen = start_pos % 4; - - if (initLen != 0) - { - memcpy(r, bit_t_seq_table[src[i>>2]] + initLen, 4 - initLen); - copyLen = copyLen + 4 - initLen; - i = i + copyLen; - } - while (copyLen < length) - { - memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4); - copyLen = copyLen + 4; - i = i + 4; - } - - if (R_INF->N_site[ID]) - { - for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) - { - if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)R_INF->N_site[ID][i] <= end_pos) - { - r[R_INF->N_site[ID][i] - start_pos] = 'N'; - } - else if((long long)R_INF->N_site[ID][i] > end_pos) - { - break; - } - } - } - } - else - { - start_pos = readLen - start_pos - 1; - end_pos = readLen - end_pos - 1; - - ///start_pos > end_pos - i = start_pos; - copyLen = 0; - long long initLen = (start_pos + 1) % 4; - - if (initLen != 0) - { - memcpy(r, bit_t_seq_table_rc[src[i>>2]] + 4 - initLen, initLen); - copyLen = copyLen + initLen; - i = i - initLen; - } - - while (copyLen < length) - { - memcpy(r+copyLen, bit_t_seq_table_rc[src[i>>2]], 4); - copyLen = copyLen + 4; - i = i - 4; - } - - if (R_INF->N_site[ID]) - { - long long offset = readLen - start_pos - 1; - - for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) - { - if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos) - { - r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N'; - } - else if((long long)R_INF->N_site[ID][i] > start_pos) - { - break; - } - } - } - } -} void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID) { @@ -515,6 +428,104 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, } +void recover_UC_sub_Read(UC_Read* i_r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID) +{ + i_r->length = length;i_r->RID = ID; + if (i_r->length + 8 > i_r->size) + { + i_r->size = i_r->length + 4; + i_r->seq = (char*)realloc(i_r->seq,sizeof(char)*(i_r->size)); + } + char* r = i_r->seq; + long long readLen = Get_READ_LENGTH((*R_INF), ID); + uint8_t* src = Get_READ((*R_INF), ID); + + long long i; + long long copyLen; + long long end_pos = start_pos + length - 1; + + if (strand == 0) + { + i = start_pos; + copyLen = 0; + + long long initLen = start_pos % 4; + + if (initLen != 0) + { + memcpy(r, bit_t_seq_table[src[i>>2]] + initLen, 4 - initLen); + copyLen = copyLen + 4 - initLen; + i = i + copyLen; + } + + while (copyLen < length) + { + memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4); + copyLen = copyLen + 4; + i = i + 4; + } + + if (R_INF->N_site[ID]) + { + for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) + { + if ((long long)R_INF->N_site[ID][i] >= start_pos && (long long)R_INF->N_site[ID][i] <= end_pos) + { + r[R_INF->N_site[ID][i] - start_pos] = 'N'; + } + else if((long long)R_INF->N_site[ID][i] > end_pos) + { + break; + } + } + } + } + else + { + start_pos = readLen - start_pos - 1; + end_pos = readLen - end_pos - 1; + + ///start_pos > end_pos + i = start_pos; + copyLen = 0; + long long initLen = (start_pos + 1) % 4; + + if (initLen != 0) + { + memcpy(r, bit_t_seq_table_rc[src[i>>2]] + 4 - initLen, initLen); + copyLen = copyLen + initLen; + i = i - initLen; + } + + while (copyLen < length) + { + memcpy(r+copyLen, bit_t_seq_table_rc[src[i>>2]], 4); + copyLen = copyLen + 4; + i = i - 4; + } + + if (R_INF->N_site[ID]) + { + long long offset = readLen - start_pos - 1; + + for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) + { + if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos) + { + r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N'; + } + else if((long long)R_INF->N_site[ID][i] > start_pos) + { + break; + } + } + } + } +} + + + + void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID) { r->length = Get_READ_LENGTH((*R_INF), ID); diff --git a/Process_Read.h b/Process_Read.h index 01cd42c..4b27d81 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -169,4 +169,5 @@ void destory_All_reads(All_reads* r); int destory_read_bin(All_reads* r); void init_Debug_reads(Debug_reads* x, const char* file); void destory_Debug_reads(Debug_reads* x); +void recover_UC_sub_Read(UC_Read* i_r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID); #endif