diff --git a/Overlaps.cpp b/Overlaps.cpp index a8dca53..a155c7c 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -22279,7 +22279,7 @@ kvec_asg_arc_t_warp* new_rtg_edges) - enable_debug_mode(1); + ///enable_debug_mode(1); purge_dups(*ug, read_g, coverage_cut, reverse_sources, ruIndex, new_rtg_edges, 0.75, 50, 50, 0.5, max_hang, min_ovlp, bubble_dist, drop_ratio, 0); @@ -22301,7 +22301,7 @@ kvec_asg_arc_t_warp* new_rtg_edges) ///debug_purge_dup = 1; ///deduplicate_advance(*ug, read_g, coverage_cut, sources, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, 0); ///delete_useless_nodes(ug); - enable_debug_mode(0); + ///enable_debug_mode(0); purge_dups(*ug, read_g, coverage_cut, reverse_sources, ruIndex, new_rtg_edges, 0.75, 50, 50, 0.5, max_hang, min_ovlp, bubble_dist, drop_ratio, 0); delete_useless_nodes(ug); diff --git a/Overlaps.h b/Overlaps.h index 03c5d09..0bb802e 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -365,11 +365,21 @@ typedef struct { uint32_t pre_n_seq, seqID; } C_graph; +typedef struct { + kvec_t(uint8_t) a; + uint32_t i; +} kvec_t_u8_warp; + typedef struct { kvec_t(uint32_t) a; uint32_t i; } kvec_t_u32_warp; +typedef struct { + kvec_t(int32_t) a; + uint32_t i; +} kvec_t_i32_warp; + typedef struct { kvec_t(uint64_t) a; uint64_t i; diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index b976330..6588e38 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -8,6 +8,8 @@ #include "kthread.h" #define Cal_Off(OFF) ((long long)((uint32_t)((OFF)>>32)) - (long long)((uint32_t)((OFF)))) +#define Get_xOff(OFF) ((long long)((uint32_t)((OFF)>>32))) +#define Get_yOff(OFF) ((long long)((uint32_t)((OFF)))) #define Get_match(x) ((x).weight) #define Get_total(x) ((x).index_beg) #define Get_type(x) ((x).index_end) @@ -40,6 +42,7 @@ typedef struct { uint32_t index_beg; uint32_t index_end; uint8_t rev; + asg_arc_t t; }hap_candidates; typedef struct { @@ -83,31 +86,13 @@ typedef struct { uint8_t* visit; kvec_t_u64_warp u_vecs; kvec_asg_arc_t_offset u_buffer; + kvec_t_i32_warp u_buffer_tailIndex; + kvec_t_i32_warp u_buffer_prevIndex; + kvec_t_u8_warp u_buffer_flag; + kvec_t_i32_warp u_buffer_beg; kvec_hap_candidates u_can; }hap_alignment_struct; -void init_hap_alignment_struct(hap_alignment_struct* x, uint32_t size) -{ - x->vote_counting = (uint64_t*)malloc(sizeof(uint64_t)*size); - memset(x->vote_counting, 0, sizeof(uint64_t)*size); - - x->visit = (uint8_t*)malloc(sizeof(uint8_t)*size); - memset(x->visit, 0, size); - - kv_init(x->u_vecs.a); - kv_init(x->u_buffer.a); - kv_init(x->u_can.a); -} - -void destory_hap_alignment_struct(hap_alignment_struct* x) -{ - free(x->vote_counting); - free(x->visit); - kv_destroy(x->u_vecs.a); - kv_destroy(x->u_buffer.a); - kv_destroy(x->u_can.a); -} - typedef struct { hap_alignment_struct* buf; @@ -126,6 +111,37 @@ typedef struct { hap_overlaps_list* all_ovlp; }hap_alignment_struct_pip; + +void init_hap_alignment_struct(hap_alignment_struct* x, uint32_t size) +{ + x->vote_counting = (uint64_t*)malloc(sizeof(uint64_t)*size); + memset(x->vote_counting, 0, sizeof(uint64_t)*size); + + x->visit = (uint8_t*)malloc(sizeof(uint8_t)*size); + memset(x->visit, 0, size); + + kv_init(x->u_vecs.a); + kv_init(x->u_buffer.a); + kv_init(x->u_buffer_tailIndex.a); + kv_init(x->u_buffer_prevIndex.a); + kv_init(x->u_buffer_beg.a); + kv_init(x->u_buffer_flag.a); + kv_init(x->u_can.a); +} + +void destory_hap_alignment_struct(hap_alignment_struct* x) +{ + free(x->vote_counting); + free(x->visit); + kv_destroy(x->u_vecs.a); + kv_destroy(x->u_buffer.a); + kv_destroy(x->u_buffer_tailIndex.a); + kv_destroy(x->u_buffer_prevIndex.a); + kv_destroy(x->u_buffer_beg.a); + kv_destroy(x->u_buffer_flag.a); + kv_destroy(x->u_can.a); +} + void init_hap_alignment_struct_pip(hap_alignment_struct_pip* x, uint32_t num_threads, uint32_t n_seq, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, ma_sub_t *coverage_cut, uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chain_rate, hap_overlaps_list* all_ovlp) @@ -386,6 +402,21 @@ int cmp_hap_alignment(const void * a, const void * b) return 0; } + +int cmp_hap_alignment_chaining(const void * a, const void * b) +{ + if((*(asg_arc_t_offset*)a).x.el > (*(asg_arc_t_offset*)b).x.el) return 1; + if((*(asg_arc_t_offset*)a).x.el < (*(asg_arc_t_offset*)b).x.el) return -1; + + if(((*(asg_arc_t_offset*)a).Off>>32) > ((*(asg_arc_t_offset*)b).Off>>32)) return 1; + if(((*(asg_arc_t_offset*)a).Off>>32) < ((*(asg_arc_t_offset*)b).Off>>32)) return -1; + + if((uint32_t)((*(asg_arc_t_offset*)a).Off) > (uint32_t)((*(asg_arc_t_offset*)b).Off)) return 1; + if((uint32_t)((*(asg_arc_t_offset*)a).Off) < (uint32_t)((*(asg_arc_t_offset*)b).Off)) return -1; + + return 0; +} + int cmp_hap_candidates(const void * a, const void * b) { if((*(hap_candidates*)a).weight < (*(hap_candidates*)b).weight) return 1; @@ -1058,6 +1089,218 @@ uint32_t dir, uint32_t rev, uint32_t* x_off, uint32_t* y_off) } +void print_asg_arc_t_offset(asg_arc_t_offset* x, long long n, const char* info) +{ + fprintf(stderr,"\n\n(%s)n: %lld\n", info, n); + long long i, x_off, y_off; + for (i = 0; i < n; i++) + { + x_off = (long long)(x[i].Off>>32); + y_off = (long long)((uint32_t)x[i].Off); + fprintf(stderr, "i: %lld, x_off: %lld, y_off: %lld, weight: %lu, rev: %u, ol: %u\n", + i, x_off, y_off, x[i].weight, x[i].x.el, x[i].x.ol); + } +} + + +// Binary search +inline int GetCeilIndex(asg_arc_t_offset* arr, kvec_t_i32_warp* T, int l, int r, uint32_t key) +{ + while (r - l > 1) { + int m = l + (r - l) / 2; + if (Get_yOff(arr[T->a.a[m]].Off) >= key) + r = m; + else + l = m; + } + + return r; +} + + +void quick_LIS(asg_arc_t_offset* x, uint32_t n, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex) +{ + tailIndex->a.n = prevIndex->a.n = 0; + if(n == 0) return; + + + kv_resize(int32_t, tailIndex->a, n); + kv_resize(int32_t, prevIndex->a, n); + + long long len = 1, i, pos, m; + tailIndex->a.a[0] = 0; + prevIndex->a.a[0] = -1; + + ///x has already sorted by x_pos + for(i = 1; i < (long long)n; i++) + { + if(Get_yOff(x[i].Off) < Get_yOff(x[tailIndex->a.a[0]].Off)) + { + // new smallest value + tailIndex->a.a[0] = i; + } + else if(Get_yOff(x[i].Off) > Get_yOff(x[tailIndex->a.a[len - 1]].Off)) + { + // arr[i] wants to extend largest subsequence + prevIndex->a.a[i] = tailIndex->a.a[len - 1]; + tailIndex->a.a[len++] = i; + } + else + { + // arr[i] wants to be a potential condidate of + // future subsequence + // It will replace ceil value in tailIndices + pos = GetCeilIndex(x, tailIndex, -1, len - 1, Get_yOff(x[i].Off)); + + prevIndex->a.a[i] = tailIndex->a.a[pos - 1]; + tailIndex->a.a[pos] = i; + } + } + + + for (m = 0, i = tailIndex->a.a[len - 1]; m < len; i = prevIndex->a.a[i], m++) + { + tailIndex->a.a[len-m-1] = i; + } + + tailIndex->a.n = len; +} + +void get_base_boundary_advance(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, +asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads, +uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, long long yEndIndex, +uint32_t rev, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, +uint32_t* xBeg, uint32_t* xEnd, uint32_t* yBeg, uint32_t* yEnd) +{ + long long k, j, offset, m; + ma_hit_t_alloc *xR = NULL; + ma_hit_t *h = NULL; + ma_sub_t *sq = NULL, *st = NULL; + int32_t r; + asg_arc_t t; + uint32_t rId, Hap_uId, is_Unitig, v, w, v_dir, w_dir; + uint64_t tmp; + asg_arc_t_offset t_offset; + u_buffer->a.n = 0; + (*xBeg) = (*xEnd) = (*yBeg) = (*yEnd) = (uint32_t)-1; + for (k = xBegIndex; k <= xEndIndex; k++) + { + xR = &(reverse_sources[xReads->a[k]>>33]); + for (j = 0; j < xR->length; j++) + { + h = &(xR->buffer[j]); + sq = &(coverage_cut[Get_qn(*h)]); + st = &(coverage_cut[Get_tn(*h)]); + if(st->del || read_g->seq[Get_tn(*h)].del) continue; + + r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, + asm_opt.max_hang_rate, min_ovlp, &t); + ///if it is a contained overlap, skip + if(r < 0) continue; + + rId = t.v>>1; + if(read_g->seq[rId].del == 1) continue; + ///there are two cases: + ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id + ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1 + get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig); + if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; + if(Hap_uId != yUid) continue; + + v = xReads->a[k]>>32; + get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig); + if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; + if(Hap_uId != xUid) continue; + if((uint32_t)(position_index[v>>1]) != k) continue; + + w = (yReads->a[(uint32_t)(position_index[rId])])>>32; + + v_dir = ((t.ul>>32)==v)?1:0; + w_dir = (t.v == w)?1:0; + if(rev == 0 && v_dir != w_dir) continue; + if(rev == 1 && v_dir == w_dir) continue; + + /****************************may have bugs********************************/ + offset = (uint32_t)(position_index[rId]); + if(offset < yBegIndex || offset > yEndIndex) continue; + /****************************may have bugs********************************/ + + tmp = get_xy_pos(read_g, &t, v, w, xReads->len, yReads->len, position_index, &(t.el)); + if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue; + + t_offset.Off = tmp; + t_offset.x = t; + t_offset.weight = 1; + kv_push(asg_arc_t_offset, u_buffer->a, t_offset); + } + } + if(u_buffer->a.n == 0) return; + + qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment_chaining); + + ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "before"); + + for (k = 1, m = 1; k < (long long)u_buffer->a.n; k++) + { + if(u_buffer->a.a[m-1].Off == u_buffer->a.a[k].Off) + { + u_buffer->a.a[m-1].weight += u_buffer->a.a[k].weight; + if(u_buffer->a.a[k].x.ol > u_buffer->a.a[m-1].x.ol) + { + u_buffer->a.a[m-1].x = u_buffer->a.a[k].x; + } + continue; + } + u_buffer->a.a[m] = u_buffer->a.a[k]; + m++; + } + u_buffer->a.n = m; + + ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n, "after"); + + quick_LIS(u_buffer->a.a, u_buffer->a.n, tailIndex, prevIndex); + + if(tailIndex->a.n == 0) return; + + (*xBeg) = Get_xOff(u_buffer->a.a[tailIndex->a.a[0]].Off); + (*yBeg) = Get_yOff(u_buffer->a.a[tailIndex->a.a[0]].Off); + (*xEnd) = Get_xOff(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].Off); + (*yEnd) = Get_yOff(u_buffer->a.a[tailIndex->a.a[tailIndex->a.n-1]].Off); +} + + +uint32_t determine_hap_overlap_type_advance(hap_candidates* hap_can, ma_utg_t *xReads, ma_utg_t *yReads, +R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, uint64_t* position_index, +int max_hang, int min_ovlp, uint32_t xUid, uint32_t yUid, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, +kvec_t_i32_warp* prevIndex, long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end) +{ + uint32_t x_pos_beg, y_pos_beg, x_pos_end, y_pos_end; + /*************************x***************************/ + get_base_boundary_advance(ruIndex, reverse_sources, coverage_cut, read_g, position_index, + max_hang, min_ovlp, xReads, yReads, xUid, yUid, Get_x_beg(*hap_can), Get_x_end(*hap_can), + Get_y_beg(*hap_can), Get_y_end(*hap_can), Get_rev(*hap_can), u_buffer, tailIndex, prevIndex, + &x_pos_beg, &x_pos_end, &y_pos_beg, &y_pos_end); + /*************************x***************************/ + + if(x_pos_beg == (uint32_t)-1 || y_pos_beg == (uint32_t)-1 + || x_pos_end == (uint32_t)-1 || y_pos_end == (uint32_t)-1) + { + return (uint32_t)-1; + } + if(x_pos_beg > x_pos_end || y_pos_beg > y_pos_end) return (uint32_t)-1; + + /** + #define X2Y 0 + #define Y2X 1 + #define XCY 2 + #define YCX 3 + **/ + return classify_hap_overlap(x_pos_beg, x_pos_end, xReads->len, y_pos_beg, y_pos_end, yReads->len, + r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end); +} + + + uint32_t determine_hap_overlap_type(hap_candidates* hap_can, ma_utg_t *xReads, ma_utg_t *yReads, R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, uint32_t xUid, uint32_t yUid, @@ -1214,6 +1457,119 @@ long long* r_y_pos_beg, long long* r_y_pos_end) } +uint32_t calculate_pair_hap_similarity_advance(hap_candidates* hap_can, +uint64_t* position_index, uint32_t xUid, uint32_t yUid, ma_utg_t* xReads, ma_utg_t* yReads, +ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, ma_sub_t *coverage_cut, +float Hap_rate, int max_hang, int min_ovlp, kvec_asg_arc_t_offset* u_buffer, +kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, long long* r_x_pos_beg, +long long* r_x_pos_end, long long* r_y_pos_beg, long long* r_y_pos_end) +{ + uint32_t max_count = 0, min_count = 0, flag; + uint32_t xLen = xReads->n, xIndex; + uint32_t yLen = yReads->n, yIndex; + uint32_t xLeftBeg, xLeftLen, yLeftBeg, yLeftLen; + uint32_t xRightBeg, xRightLen, yRightBeg, yRightLen; + double xLeftMatch = 0, xLeftTotal = 0, yLeftMatch = 0, yLeftTotal = 0; + double xRightMatch = 0, xRightTotal = 0, yRightMatch = 0, yRightTotal = 0; + asg_arc_t* arch = NULL; + + + arch = &(hap_can->t); + xIndex = (uint32_t)(position_index[arch->ul>>33]); + yIndex = (uint32_t)(position_index[arch->v>>1]); + + if(hap_can->rev == 0) + { + xLeftBeg = 0; xLeftLen = xIndex; xRightBeg = xIndex; xRightLen = xLen - xRightBeg; + yLeftBeg = 0; yLeftLen = yIndex; yRightBeg = yIndex; yRightLen = yLen - yRightBeg; + } + else + { + xLeftBeg = 0; xLeftLen = xIndex; xRightBeg = xIndex; xRightLen = xLen - xRightBeg; + + yLeftBeg = yIndex + 1; yLeftLen = yLen - yLeftBeg; + yRightBeg = 0; yRightLen = yIndex + 1; + } + + flag = Get_type(*hap_can); + + + if(flag == XCY) + { + get_pair_hap_similarity(yReads->a, yLen, xUid, reverse_sources, read_g, ruIndex, + &yLeftMatch, &yLeftTotal); + max_count = yLeftMatch; + min_count = yLeftTotal; + } + else if(flag == YCX) + { + get_pair_hap_similarity(xReads->a, xLen, yUid, reverse_sources, read_g, ruIndex, + &xLeftMatch, &xLeftTotal); + max_count = xLeftMatch; + min_count = xLeftTotal; + } + else if(flag == X2Y) + { + get_pair_hap_similarity(yReads->a+yLeftBeg, yLeftLen, xUid, reverse_sources, read_g, ruIndex, + &yLeftMatch, &yLeftTotal); + get_pair_hap_similarity(xReads->a+xRightBeg, xRightLen, yUid, reverse_sources, read_g, ruIndex, + &xRightMatch, &xRightTotal); + max_count = yLeftMatch + xRightMatch; + min_count = yLeftTotal + xRightTotal; + } + else if(flag == Y2X) + { + get_pair_hap_similarity(xReads->a+xLeftBeg, xLeftLen, yUid, reverse_sources, read_g, ruIndex, + &xLeftMatch, &xLeftTotal); + get_pair_hap_similarity(yReads->a+yRightBeg, yRightLen, xUid, reverse_sources, read_g, ruIndex, + &yRightMatch, &yRightTotal); + max_count = xLeftMatch + yRightMatch; + min_count = xLeftTotal + yRightTotal; + } else abort(); + + hap_can->weight = hap_can->index_beg = 0; + if(min_count == 0) return NON_PLOID; + if(max_count > min_count*Hap_rate) + { + long long r_x_interval_beg, r_x_interval_end, r_y_interval_beg, r_y_interval_end; + + ///for containment, don't need to do anything + get_hap_alignment_boundary(xReads, yReads, flag, xLeftMatch, xLeftTotal, + yLeftMatch, yLeftTotal, xRightMatch, xRightTotal, yRightMatch, yRightTotal, + xLeftBeg, xLeftLen, yLeftBeg, yLeftLen, xRightBeg, xRightLen, yRightBeg, yRightLen, + xUid, yUid, Hap_rate, position_index, reverse_sources, read_g, ruIndex, hap_can->rev, + &r_x_interval_beg, &r_x_interval_end, &r_y_interval_beg, &r_y_interval_end); + + if(r_x_interval_beg < 0 || r_x_interval_end < 0 || r_y_interval_beg < 0 || r_y_interval_end < 0) + { + return NON_PLOID; + } + + get_pair_hap_similarity(xReads->a + r_x_interval_beg, r_x_interval_end + 1 - r_x_interval_beg, + yUid, reverse_sources, read_g, ruIndex, &xLeftMatch, &xLeftTotal); + if(xLeftMatch == 0 || xLeftTotal == 0) return NON_PLOID; + + hap_can->weight = xLeftMatch; + hap_can->index_beg = xLeftTotal; + hap_can->index_end = flag; + hap_can->x_beg_pos = r_x_interval_beg; + hap_can->x_end_pos = r_x_interval_end; + hap_can->y_beg_pos = r_y_interval_beg; + hap_can->y_end_pos = r_y_interval_end; + + hap_can->index_end = determine_hap_overlap_type_advance(hap_can, xReads, yReads, + ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, min_ovlp, + xUid, yUid, u_buffer, tailIndex, prevIndex, r_x_pos_beg, r_x_pos_end, r_y_pos_beg, + r_y_pos_end); + if(hap_can->index_end == XCY && yReads->len > (xReads->len*2)) return NON_PLOID; + if(hap_can->index_end == YCX && xReads->len > (yReads->len*2)) return NON_PLOID; + if(hap_can->index_end == (uint32_t)-1) return NON_PLOID; + + return PLOID; + } + return NON_PLOID; +} + void print_hap_paf(ma_ug_t *ug, hap_overlaps* ovlp) { @@ -1493,7 +1849,285 @@ hap_overlaps_list* all_ovlp) } -static void hap_alignment_worker(void *_data, long eid, int tid) + +inline long long get_max_index(asg_arc_t_offset* x, int32_t* Scores, uint8_t* Flag, long long n, +long long x_readLen, long long y_readLen) +{ + long long i = 0, max_result = -1, max_i = -1, min_xLen = x_readLen * 2 + 2, x_off, y_off, tmp_xLen; + for (i = 0; i < n; i++) + { + if(Flag[i] != 0) continue; + x_off = (long long)(x[i].Off>>32); + y_off = (long long)((uint32_t)x[i].Off); + if(Scores[i] > max_result) + { + max_result = Scores[i]; + max_i = i; + min_xLen = get_hap_overlapLen(x_off, x_off, x_readLen, y_off, y_off, y_readLen, + NULL, NULL, NULL, NULL); + } + else if(Scores[i] == max_result) + { + tmp_xLen = get_hap_overlapLen(x_off, x_off, x_readLen, y_off, y_off, y_readLen, + NULL, NULL, NULL, NULL); + + if(tmp_xLen < min_xLen) + { + max_result = Scores[i]; + max_i = i; + min_xLen = tmp_xLen; + } + } + } + + return max_i; +} + + +inline void get_chain_details(int32_t* Pres, int32_t* Results, uint8_t* Flag, long long max_i, +long long* chainLen, long long* dup) +{ + long long i = max_i; + (*chainLen) = 0; + (*dup) = 0; + + while (i >= 0) + { + if(Flag[i] == 1) (*dup)++; + Results[(*chainLen)] = i; + i = Pres[i]; + (*chainLen)++; + } +} + +inline void push_hap_can(asg_arc_t_offset* x, kvec_hap_candidates* u_can, int32_t* Results, +long long chainLen, long long x_readLen, long long y_readLen) +{ + if(chainLen <= 0) return; + hap_candidates hap_can; + hap_can.rev = x[Results[0]].x.el; + hap_can.x_beg_pos = hap_can.x_end_pos = (uint32_t)(x[Results[0]].Off>>32); + hap_can.y_beg_pos = hap_can.y_end_pos = (uint32_t)(x[Results[0]].Off); + hap_can.weight = 0; + long long i = 0; + uint64_t totalWeigth = 0; + ///fprintf(stderr, "^^^chainLen: %lld\n", chainLen); + for (i = 0; i < chainLen; i++) + { + ///fprintf(stderr, "i: %lld, Results[i]: %d\n", i, Results[i]); + hap_can.x_beg_pos = (uint32_t)(x[Results[i]].Off>>32); + hap_can.y_beg_pos = (uint32_t)(x[Results[i]].Off); + hap_can.weight += x[Results[i]].weight; + } + + for (i = 0; i < chainLen; i++) + { + totalWeigth += x[Results[i]].weight; + if(totalWeigth >= (hap_can.weight/2)) break; + } + + if(i >= chainLen) i = chainLen-1; + ///Get_total(hap_can) = Results[i]; + hap_can.t = x[Results[i]].x; + Get_type(hap_can) = classify_hap_overlap(hap_can.x_beg_pos, hap_can.x_end_pos, + x_readLen, hap_can.y_beg_pos, hap_can.y_end_pos, y_readLen, NULL, NULL, NULL, NULL); + kv_push(hap_candidates, u_can->a, hap_can); + if(hap_can.x_beg_pos > hap_can.x_end_pos || hap_can.y_beg_pos > hap_can.y_end_pos) + { + fprintf(stderr, "ERROR\n"); + } +} + + + +void print_chain_data(int32_t* Scores, int32_t* Pres, int32_t* Begs, long long n) +{ + fprintf(stderr,"*****\nn_chain: %lld\n", n); + long long i; + for (i = 0; i < n; i++) + { + fprintf(stderr, "i: %lld, Scores: %d, Pres: %d, Begs: %d\n", + i, Scores[i], Pres[i], Begs[i]); + } +} + + + +void hap_chaining(asg_arc_t_offset* x, uint32_t n, kvec_t_i32_warp* score_vc, kvec_t_i32_warp* prevIndex_vec, +kvec_t_i32_warp* begIndex_vec, kvec_t_u8_warp* flag_vec, float band_width_threshold, long long max_skip, +long long x_readLen, long long y_readLen, kvec_hap_candidates* u_can) +{ + #define DUP_OVLP_RATE 0.75 + score_vc->a.n = prevIndex_vec->a.n = begIndex_vec->a.n = flag_vec->a.n = 0; + if(n == 0) return; + kv_resize(int32_t, score_vc->a, n); + kv_resize(int32_t, prevIndex_vec->a, n); + kv_resize(int32_t, begIndex_vec->a, n); + kv_resize(uint8_t, flag_vec->a, n); + + int32_t* Scores = score_vc->a.a; + int32_t* Pres = prevIndex_vec->a.a; + int32_t* Begs = begIndex_vec->a.a; + uint8_t* Flag = flag_vec->a.a; + long long i, j, n_max_skip, x_off, y_off, max_beg, max_j = -1, max_score, score; + long long distance_x, distance_y, total_distance_x, total_distance_y, distance_gap; + float gap_rate, band_width_penalty = 1 / band_width_threshold; + long long max_result, max_i, min_xLen, tmp_xLen, chainLen = 0, dup = 0; + max_result = max_i = -1; min_xLen = x_readLen * 2 + 2; + for (i = 0; i < n; i++) + { + n_max_skip = 0; + + x_off = (long long)(x[i].Off>>32); + y_off = (long long)((uint32_t)x[i].Off); + max_j = -1; + max_score = x[i].weight; + max_beg = i; //i itself + ///may have a pre-cut condition for j + for (j = i - 1; j >= 0; --j) + { + distance_x = x_off - (long long)(x[j].Off>>32); + distance_y = y_off - (long long)((uint32_t)x[j].Off); + ///x has been sorted by x_off + if(distance_x <= 0 || distance_y <= 0) continue; + + total_distance_x = x_off - (long long)(x[Begs[j]].Off>>32); + total_distance_y = y_off - (long long)((uint32_t)x[Begs[j]].Off); + + distance_gap = total_distance_x - total_distance_y; + if(distance_gap < 0) distance_gap = -distance_gap; + if(distance_gap > band_width_threshold * total_distance_x) + { + continue; + } + + score = x[i].weight; + gap_rate = (float)((float)(distance_gap)/(float)(total_distance_x)); + score -= (long long)(score * gap_rate * band_width_penalty); + score += Scores[j]; + + ///find a new max score + if (score > max_score) { + max_score = score; + max_j = j; + max_beg = Begs[j]; + n_max_skip = 0; + } + else + { + if (++n_max_skip > max_skip) break; + } + } + + Scores[i] = max_score; + Pres[i] = max_j; + Begs[i] = max_beg; + + if(Scores[i] > max_result) + { + max_result = Scores[i]; + max_i = i; + min_xLen = get_hap_overlapLen(x_off, x_off, x_readLen, y_off, y_off, y_readLen, + NULL, NULL, NULL, NULL); + } + else if(Scores[i] == max_result) + { + tmp_xLen = get_hap_overlapLen(x_off, x_off, x_readLen, y_off, y_off, y_readLen, + NULL, NULL, NULL, NULL); + + if(tmp_xLen < min_xLen) + { + max_result = Scores[i]; + max_i = i; + min_xLen = tmp_xLen; + } + } + Flag[i] = 0; + } + + // print_asg_arc_t_offset(x, n); + // print_chain_data(Scores, Pres, Begs, n); + while (max_i != -1) + { + get_chain_details(Pres, Begs, Flag, max_i, &chainLen, &dup); + if(chainLen == 0) break; + if(dup > chainLen*DUP_OVLP_RATE) + { + for (i = 0; i < chainLen; i++) + { + if(Flag[Begs[i]] == 1) continue; + Flag[Begs[i]] = 2; + } + + } + else + { + push_hap_can(x, u_can, Begs, chainLen, x_readLen, y_readLen); + + for (i = 0; i < chainLen; i++) + { + Flag[Begs[i]] = 1; + } + } + + max_i = get_max_index(x, Scores, Flag, n, x_readLen, y_readLen); + } +} + +void get_candidate_hap_alignment(kvec_hap_candidates* u_can, kvec_asg_arc_t_offset* u_buffer, +kvec_t_i32_warp* score_vc, kvec_t_i32_warp* prevIndex_vec, kvec_t_i32_warp* begIndex_vec, +kvec_t_u8_warp* flag_vec, float band_width_threshold, long long max_skip, long long x_readLen, +long long y_readLen) +{ + u_can->a.n = 0; + if(u_buffer->a.n == 0) return; + uint32_t i = 0, anchor_i = 0, m = 1, break_point = (uint32_t)-1, is_merge; + + qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment_chaining); + + ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n); + + for (i = 1; i < u_buffer->a.n; i++) + { + is_merge = 0; + if(u_buffer->a.a[m-1].x.el == u_buffer->a.a[i].x.el) + { + if(u_buffer->a.a[m-1].Off == u_buffer->a.a[i].Off) is_merge = 1; + if(is_merge == 0 && (Get_xOff(u_buffer->a.a[m-1].Off)==Get_xOff(u_buffer->a.a[i].Off))) + { + if((Get_yOff(u_buffer->a.a[i].Off)-(Get_yOff(u_buffer->a.a[m-1].Off))) == + (i-anchor_i)) + { + is_merge = 1; + } + } + + if(is_merge) + { + u_buffer->a.a[m-1].weight += u_buffer->a.a[i].weight; + continue; + } + } + u_buffer->a.a[m] = u_buffer->a.a[i]; + anchor_i = i; + if(u_buffer->a.a[m].x.el != u_buffer->a.a[m-1].x.el) break_point = m; + m++; + } + u_buffer->a.n = m; + if(break_point > u_buffer->a.n) break_point = u_buffer->a.n; + + ///print_asg_arc_t_offset(u_buffer->a.a, u_buffer->a.n); + + hap_chaining(u_buffer->a.a, break_point, score_vc, prevIndex_vec, begIndex_vec, flag_vec, + band_width_threshold, max_skip, x_readLen, y_readLen, u_can); + + hap_chaining(u_buffer->a.a + break_point, u_buffer->a.n - break_point, score_vc, prevIndex_vec, + begIndex_vec, flag_vec, band_width_threshold, max_skip, x_readLen, y_readLen, u_can); +} + + +///static void hap_alignment_worker(void *_data, long eid, int tid) +void hap_alignment_worker(void *_data, long eid, int tid) { hap_alignment_struct_pip* hap_buf = (hap_alignment_struct_pip*)_data; ma_ug_t *ug = hap_buf->ug; @@ -1783,6 +2417,230 @@ static void hap_alignment_worker(void *_data, long eid, int tid) } + +static void hap_alignment_advance_worker(void *_data, long eid, int tid) +{ + hap_alignment_struct_pip* hap_buf = (hap_alignment_struct_pip*)_data; + ma_ug_t *ug = hap_buf->ug; + asg_t *read_g = hap_buf->read_g; + ma_hit_t_alloc* reverse_sources = hap_buf->reverse_sources; + R_to_U* ruIndex = hap_buf->ruIndex; + ma_sub_t *coverage_cut = hap_buf->coverage_cut; + uint64_t* position_index = hap_buf->position_index; + float Hap_rate = hap_buf->Hap_rate; + int max_hang = hap_buf->max_hang; + int min_ovlp = hap_buf->min_ovlp; + float chain_rate = hap_buf->chain_rate; + hap_overlaps_list* all_ovlp = hap_buf->all_ovlp; + uint32_t Input_uId = eid; + uint64_t* vote_counting = hap_buf->buf[tid].vote_counting; + uint8_t* visit = hap_buf->buf[tid].visit; + kvec_t_u64_warp* u_vecs = &(hap_buf->buf[tid].u_vecs); + kvec_asg_arc_t_offset* u_buffer = &(hap_buf->buf[tid].u_buffer); + kvec_hap_candidates* u_can = &(hap_buf->buf[tid].u_can); + kvec_t_i32_warp* score_vc = &(hap_buf->buf[tid].u_buffer_tailIndex); + kvec_t_i32_warp* prevIndex_vec = &(hap_buf->buf[tid].u_buffer_prevIndex); + kvec_t_i32_warp* begIndex_vec = &(hap_buf->buf[tid].u_buffer_beg); + kvec_t_u8_warp* flag_vec = &(hap_buf->buf[tid].u_buffer_flag); + + ma_utg_t *xReads = NULL, *yReads = NULL; + ma_hit_t_alloc *xR = NULL; + ma_hit_t *h = NULL; + ma_sub_t *sq = NULL, *st = NULL; + asg_t* nsg = ug->g; + uint32_t i, j, v, rId, k, is_Unitig, Hap_uId, xUid, yUid, seedOcc, is_update; + uint64_t tmp; + long long r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end; + int32_t r; + asg_arc_t t; + asg_arc_t_offset t_offset; + hap_candidates hap_can; + hap_overlaps hap_align; + xUid = Input_uId; + if(nsg->seq[xUid].del || nsg->seq[xUid].c == ALTER_LABLE) return; + memset(vote_counting, 0, sizeof(uint64_t)*nsg->n_seq); + memset(visit, 0, nsg->n_seq); + u_vecs->a.n = 0; + u_can->a.n = 0; + + xReads = &(ug->u.a[xUid]); + for (i = 0; i < xReads->n; i++) + { + xR = &(reverse_sources[xReads->a[i]>>33]); + + for (k = 0; k < xR->length; k++) + { + rId = Get_tn(xR->buffer[k]); + + if(read_g->seq[rId].del == 1) + { + ///get the id of read that contains it + get_R_to_U(ruIndex, rId, &rId, &is_Unitig); + if(rId == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[rId].del == 1) continue; + } + + ///there are two cases: + ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id + ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1 + get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig); + if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; + ///here rId is the id of the read coming from the different haplotype + ///Hap_cId is the id of the corresponding contig (note here is the contig, instead of untig) + if(visit[Hap_uId]!=0) continue; + visit[Hap_uId] = 1; + if(vote_counting[Hap_uId] < UINT64_MAX) vote_counting[Hap_uId]++; + } + + clean_visit_flag(visit, read_g, ruIndex, nsg->n_seq, xR); + } + + + + u_vecs->a.n = 0; + for (i = 0; i < nsg->n_seq; i++) + { + if(i == xUid) continue; + if(vote_counting[i] == 0) continue; + tmp = vote_counting[i]; tmp = tmp << 32; tmp = tmp | (uint64_t)i; + kv_push(uint64_t, u_vecs->a, tmp); + } + + if(u_vecs->a.n == 0) return; + sort_kvec_t_u64_warp(u_vecs, 1); + + + ///scan each candidate unitig + for (i = 0; i < u_vecs->a.n; i++) + { + yUid = (uint32_t)u_vecs->a.a[i]; + seedOcc = u_vecs->a.a[i]>>32; + xReads = &(ug->u.a[xUid]); + yReads = &(ug->u.a[yUid]); + u_buffer->a.n = 0; + + for (k = 0; k < xReads->n; k++) + { + xR = &(reverse_sources[xReads->a[k]>>33]); + for (j = 0; j < xR->length; j++) + { + h = &(xR->buffer[j]); + sq = &(coverage_cut[Get_qn(*h)]); + st = &(coverage_cut[Get_tn(*h)]); + if(st->del || read_g->seq[Get_tn(*h)].del) continue; + + r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, + asm_opt.max_hang_rate, min_ovlp, &t); + ///if it is a contained overlap, skip + if(r < 0) continue; + + rId = t.v>>1; + if(read_g->seq[rId].del == 1) continue; + ///there are two cases: + ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id + ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1 + get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig); + if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; + if(Hap_uId != yUid) continue; + + v = xReads->a[k]>>32; + get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig); + if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; + if(Hap_uId != xUid) continue; + if((uint32_t)(position_index[v>>1]) != k) continue; + + + if((prefilter((uint32_t)(position_index[v>>1]), (uint32_t)(position_index[rId]), + xReads->n, yReads->n, 0, Hap_rate, seedOcc)==NON_PLOID) && + (prefilter((uint32_t)(position_index[v>>1]), (uint32_t)(position_index[rId]), + xReads->n, yReads->n, 1, Hap_rate, seedOcc)==NON_PLOID)) + { + continue; + } + + t_offset.Off = get_xy_pos(read_g, &t, v, (yReads->a[(uint32_t)(position_index[rId])])>>32, + xReads->len, yReads->len, position_index, &(t.el)); + if(((t_offset.Off>>32) == (uint32_t)-1) || (((uint32_t)t_offset.Off) == (uint32_t)-1)) continue; + + t_offset.x = t; + t_offset.weight = 1; + + kv_push(asg_arc_t_offset, u_buffer->a, t_offset); + } + + deduplicate_edge(u_buffer); + } + + if(u_buffer->a.n == 0) continue; + + + get_candidate_hap_alignment(u_can, u_buffer, score_vc, prevIndex_vec, begIndex_vec, + flag_vec, chain_rate, 50, xReads->len, yReads->len); + + if(u_can->a.n == 0) continue; + + qsort(u_can->a.a, u_can->a.n, sizeof(hap_candidates), cmp_hap_candidates); + + Get_match(hap_can) = Get_total(hap_can) = 0; + memset(&hap_align, 0, sizeof(hap_overlaps)); + + for (k = 0; k < u_can->a.n; k++) + { + is_update = 0; + if(u_can->a.a[k].weight < Get_match(hap_can)*Hap_rate) continue; + + if(calculate_pair_hap_similarity_advance(&(u_can->a.a[k]), position_index, xUid, yUid, + xReads, yReads, reverse_sources, read_g, ruIndex, coverage_cut, Hap_rate, max_hang, + min_ovlp, u_buffer, score_vc, prevIndex_vec, &r_x_pos_beg, &r_x_pos_end, &r_y_pos_beg, + &r_y_pos_end)!=PLOID) + { + continue; + } + + if(Get_match(hap_can) < Get_match(u_can->a.a[k])) + { + is_update = 1; + } + else if(Get_match(hap_can) == Get_match(u_can->a.a[k]) && + Get_total(hap_can) > Get_total(u_can->a.a[k])) + { + is_update = 1; + } + + if(is_update) + { + hap_can = u_can->a.a[k]; + hap_align.rev = Get_rev(hap_can); + hap_align.type = Get_type(hap_can); + hap_align.x_beg_id = Get_x_beg(hap_can); + hap_align.x_end_id = Get_x_end(hap_can) + 1; + hap_align.y_beg_id = Get_y_beg(hap_can); + hap_align.y_end_id = Get_y_end(hap_can) + 1; + hap_align.weight = Get_match(hap_can); + hap_align.x_beg_pos = r_x_pos_beg; + hap_align.x_end_pos = r_x_pos_end + 1; + if(hap_align.rev == 0) + { + hap_align.y_beg_pos = r_y_pos_beg; + hap_align.y_end_pos = r_y_pos_end + 1; + } + else + { + hap_align.y_beg_pos = yReads->len - r_y_pos_end - 1; + hap_align.y_end_pos = yReads->len - r_y_pos_beg - 1 + 1; + } + hap_align.xUid = xUid; + hap_align.yUid = yUid; + hap_align.status = SELF_EXIST; + } + } + + if(Get_match(hap_can) == 0 || Get_total(hap_can) == 0) continue; + + kv_push(hap_overlaps, all_ovlp->x[hap_align.xUid].a, hap_align); + } + +} + int inline get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn) { uint32_t i; @@ -1861,6 +2719,62 @@ void normalize_hap_overlaps(hap_overlaps_list* all_ovlp, hap_overlaps_list* back } } +inline uint64_t calculate_bi_weight(hap_overlaps *x, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, +R_to_U* ruIndex) +{ + double xMatch, xTotal; + uint64_t weight = x->weight; + ma_utg_t *yReads = &(ug->u.a[x->yUid]); + get_pair_hap_similarity(yReads->a + x->y_beg_id, x->y_end_id - x->y_beg_id, + x->xUid, reverse_sources, read_g, ruIndex, &xMatch, &xTotal); + weight += xMatch; + return weight; +} + +void normalize_hap_overlaps_advance(hap_overlaps_list* all_ovlp, hap_overlaps_list* back_all_ovlp, +ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex) +{ + hap_overlaps *x = NULL, *y = NULL; + uint32_t v, i, uId, qn, tn; + uint32_t types[4]; + + types[X2Y] = Y2X; types[Y2X] = X2Y; types[XCY] = YCX; types[YCX] = XCY; + int index; + for (v = 0; v < all_ovlp->num; v++) + { + uId = v; + for (i = 0; i < all_ovlp->x[uId].a.n; i++) + { + qn = all_ovlp->x[uId].a.a[i].xUid; + tn = all_ovlp->x[uId].a.a[i].yUid; + x = &(all_ovlp->x[uId].a.a[i]); + index = get_specific_hap_overlap(&(all_ovlp->x[tn]), tn, qn); + if(index != -1) + { + y = &(all_ovlp->x[tn].a.a[index]); + if(x->rev == y->rev && types[x->type]==y->type) continue; + ///if(x->weight >= y->weight) + if((calculate_bi_weight(x, ug, read_g, reverse_sources, ruIndex)) >= + (calculate_bi_weight(y, ug, read_g, reverse_sources, ruIndex))) + { + kv_push(hap_overlaps, back_all_ovlp->x[tn].a, (*y)); + set_reverse_hap_overlap(y, x, types); + } + else + { + kv_push(hap_overlaps, back_all_ovlp->x[qn].a, (*x)); + set_reverse_hap_overlap(x, y, types); + } + } + else + { + kv_pushp(hap_overlaps, all_ovlp->x[tn].a, &y); + set_reverse_hap_overlap(y, x, types); + } + } + } +} + void debug_hap_overlaps(hap_overlaps_list* all_ovlp, hap_overlaps_list* back_all_ovlp) { @@ -2178,6 +3092,160 @@ long long yEndIndex, uint32_t dir, uint32_t rev, asg_arc_t* reture_t_f, asg_arc_ } +void get_node_boundary_advance(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, +asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads, +uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, +long long yEndIndex, uint32_t dir, uint32_t rev, kvec_asg_arc_t_offset* u_buffer, +kvec_t_i32_warp* tailIndex, kvec_t_i32_warp* prevIndex, asg_arc_t* reture_t_f, asg_arc_t* reture_t_r) +{ + long long k, j, offset, m; + ma_hit_t_alloc *xR = NULL; + ma_hit_t *h = NULL; + ma_sub_t *sq = NULL, *st = NULL; + int r, index; + asg_arc_t t_f, t_r; + uint32_t rId, Hap_uId, is_Unitig, v, w, v_dir, w_dir; + uint64_t tmp; + asg_arc_t_offset t_offset; + reture_t_f->del = reture_t_r->del = 1; + u_buffer->a.n = 0; + + + for (k = xBegIndex; k <= xEndIndex; k++) + { + xR = &(reverse_sources[xReads->a[k]>>33]); + for (j = 0; j < xR->length; j++) + { + h = &(xR->buffer[j]); + sq = &(coverage_cut[Get_qn(*h)]); + st = &(coverage_cut[Get_tn(*h)]); + if(st->del || read_g->seq[Get_tn(*h)].del) continue; + + r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, + asm_opt.max_hang_rate, min_ovlp, &t_f); + ///if it is a contained overlap, skip + if(r < 0) continue; + + rId = t_f.v>>1; + if(read_g->seq[rId].del == 1) continue; + ///there are two cases: + ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id + ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1 + get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig); + if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; + if(Hap_uId != yUid) continue; + + v = xReads->a[k]>>32; + get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig); + if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue; + if(Hap_uId != xUid) continue; + if((uint32_t)(position_index[v>>1]) != k) continue; + + w = (yReads->a[(uint32_t)(position_index[rId])])>>32; + + v_dir = ((t_f.ul>>32)==v)?1:0; + w_dir = (t_f.v == w)?1:0; + if(rev == 0 && v_dir != w_dir) continue; + if(rev == 1 && v_dir == w_dir) continue; + if(dir == v_dir) continue; + + /****************************may have bugs********************************/ + offset = (uint32_t)(position_index[rId]); + if(offset < yBegIndex || offset > yEndIndex) continue; + /****************************may have bugs********************************/ + + /************************get reverse edge*************************/ + index = get_specific_overlap(&(reverse_sources[Get_tn(*h)]), Get_tn(*h), Get_qn(*h)); + if(index == -1) continue; + h = &(reverse_sources[Get_tn(*h)].buffer[index]); + sq = &(coverage_cut[Get_qn(*h)]); + st = &(coverage_cut[Get_tn(*h)]); + if(st->del || read_g->seq[Get_tn(*h)].del) continue; + r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, + asm_opt.max_hang_rate, min_ovlp, &t_r); + if(r < 0) continue; + /************************get reverse edge*************************/ + + tmp = get_xy_pos(read_g, &t_f, v, w, xReads->len, yReads->len, position_index, &(t_f.el)); + if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue; + + t_offset.Off = tmp; + t_offset.x = t_f; + t_offset.weight = 1; + kv_push(asg_arc_t_offset, u_buffer->a, t_offset); + + } + } + + if(u_buffer->a.n == 0) return; + + qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment_chaining); + + for (k = 1, m = 1; k < (long long)u_buffer->a.n; k++) + { + if(u_buffer->a.a[m-1].Off == u_buffer->a.a[k].Off) + { + u_buffer->a.a[m-1].weight += u_buffer->a.a[k].weight; + if(u_buffer->a.a[k].x.ol > u_buffer->a.a[m-1].x.ol) + { + u_buffer->a.a[m-1].x = u_buffer->a.a[k].x; + } + continue; + } + u_buffer->a.a[m] = u_buffer->a.a[k]; + m++; + } + u_buffer->a.n = m; + + quick_LIS(u_buffer->a.a, u_buffer->a.n, tailIndex, prevIndex); + + if(tailIndex->a.n == 0) return; + + if(dir == 0) + { + for (k = 0; k < (long long)tailIndex->a.n; k++) + { + v = u_buffer->a.a[tailIndex->a.a[k]].x.v>>1; + w = u_buffer->a.a[tailIndex->a.a[k]].x.ul>>33; + index = get_specific_overlap(&(reverse_sources[v]), v, w); + if(index == -1) continue; + h = &(reverse_sources[v].buffer[index]); + sq = &(coverage_cut[Get_qn(*h)]); + st = &(coverage_cut[Get_tn(*h)]); + if(st->del || read_g->seq[Get_tn(*h)].del) continue; + r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, + asm_opt.max_hang_rate, min_ovlp, &t_r); + if(r < 0) continue; + + (*reture_t_f) = u_buffer->a.a[tailIndex->a.a[k]].x; + (*reture_t_r) = t_r; + return; + } + } + else + { + for (k = tailIndex->a.n-1; k >= 0; k--) + { + v = u_buffer->a.a[tailIndex->a.a[k]].x.v>>1; + w = u_buffer->a.a[tailIndex->a.a[k]].x.ul>>33; + index = get_specific_overlap(&(reverse_sources[v]), v, w); + if(index == -1) continue; + h = &(reverse_sources[v].buffer[index]); + sq = &(coverage_cut[Get_qn(*h)]); + st = &(coverage_cut[Get_tn(*h)]); + if(st->del || read_g->seq[Get_tn(*h)].del) continue; + r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, + asm_opt.max_hang_rate, min_ovlp, &t_r); + if(r < 0) continue; + + (*reture_t_f) = u_buffer->a.a[tailIndex->a.a[k]].x; + (*reture_t_r) = t_r; + return; + } + } +} + + void fill_unitig(uint64_t* buffer, uint32_t bufferLen, asg_t* read_g, kvec_asg_arc_t_warp* edge, uint32_t is_circle, uint64_t* rLen) { @@ -2286,7 +3354,8 @@ uint32_t is_circle, uint64_t* rLen) void purge_merge(asg_t *purge_g, ma_ug_t *ug, hap_overlaps_list* all_ovlp, buf_t* b_0, R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, -uint64_t* position_index, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit) +uint64_t* position_index, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, +kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit) { uint32_t i, nv, k, v, w, x_beg_index, x_end_index, y_beg_index, y_end_index, cut_beg, cut_end, begIndex, endIndex, keepUid; hap_overlaps *x = NULL/**, *y = NULL**/; @@ -2348,9 +3417,13 @@ uint64_t* position_index, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, endIndex = x->x_end_id-1; if(cut_end < endIndex) endIndex = cut_end; - get_node_boundary(ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, + // get_node_boundary(ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, + // min_ovlp, xReads, yReads, v>>1, w>>1, begIndex, endIndex, x->y_beg_id, x->y_end_id-1, v&1, + // x->rev, &t_forward, &t_backward); + get_node_boundary_advance(ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, min_ovlp, xReads, yReads, v>>1, w>>1, begIndex, endIndex, x->y_beg_id, x->y_end_id-1, v&1, - x->rev, &t_forward, &t_backward); + x->rev, u_buffer, tailIndex, prevIndex, &t_forward, &t_backward); + if(t_forward.del || t_backward.del) break; kv_push(asg_arc_t, edge->a, t_forward); @@ -2452,7 +3525,8 @@ uint64_t* position_index, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, void link_unitigs(asg_t *purge_g, ma_ug_t *ug, hap_overlaps_list* all_ovlp, R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, asg_t *read_g, -uint64_t* position_index, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit) +uint64_t* position_index, kvec_asg_arc_t_offset* u_buffer, kvec_t_i32_warp* tailIndex, +kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, uint8_t* visit) { uint32_t v, n_vtx = purge_g->n_seq * 2, beg, end; long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen; @@ -2473,7 +3547,7 @@ uint64_t* position_index, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, continue; } purge_merge(purge_g, ug, all_ovlp, &b_0, ruIndex, reverse_sources, coverage_cut, - read_g, position_index, max_hang, min_ovlp, edge, visit); + read_g, position_index, u_buffer, tailIndex, prevIndex,max_hang, min_ovlp, edge, visit); } free(b_0.b.a); } @@ -2574,7 +3648,9 @@ uint32_t just_contain) reverse_sources, ruIndex, coverage_cut, position_index, density, max_hang, min_ovlp, 0.05, &all_ovlp); - kt_for(asm_opt.thread_num, hap_alignment_worker, &hap_buf, nsg->n_seq); + ///kt_for(asm_opt.thread_num, hap_alignment_worker, &hap_buf, nsg->n_seq); + kt_for(asm_opt.thread_num, hap_alignment_advance_worker, &hap_buf, nsg->n_seq); + ///if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp); @@ -2589,7 +3665,8 @@ uint32_t just_contain) // 0.05, &all_ovlp); // } - normalize_hap_overlaps(&all_ovlp, &back_all_ovlp); + ///normalize_hap_overlaps(&all_ovlp, &back_all_ovlp); + normalize_hap_overlaps_advance(&all_ovlp, &back_all_ovlp, ug, read_g, reverse_sources, ruIndex); ///debug_hap_overlaps(&all_ovlp, &back_all_ovlp); @@ -2662,6 +3739,7 @@ uint32_t just_contain) // if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp); link_unitigs(purge_g, ug, &all_ovlp, ruIndex, reverse_sources, coverage_cut, read_g, position_index, + &(hap_buf.buf[0].u_buffer), &(hap_buf.buf[0].u_buffer_tailIndex), &(hap_buf.buf[0].u_buffer_prevIndex), max_hang, min_ovlp, edge, hap_buf.buf[0].visit); }