diff --git a/Assembly.cpp b/Assembly.cpp index d28e238..e75b287 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1190,7 +1190,6 @@ int ha_assemble(void) build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf, R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round, asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, !ovlp_loaded); - if (asm_opt.verbose_gfa) R_INF.paf = R_INF.reverse_paf = 0; // with --dbg-gfa, R_INF.paf and R_INF.reverse_paf are freed in load_debug_graph() destory_All_reads(&R_INF); return 0; } diff --git a/Overlaps.cpp b/Overlaps.cpp index fa3deb6..e938e40 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -22194,6 +22194,8 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex) return 0; } + R_INF.paf = (*sources); R_INF.reverse_paf = (*reverse_sources); + return 1; } @@ -22276,7 +22278,8 @@ kvec_asg_arc_t_warp* new_rtg_edges) - + + 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); @@ -22286,28 +22289,25 @@ kvec_asg_arc_t_warp* new_rtg_edges) - rescue_missing_overlaps_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang, min_ovlp, 0, 0, 1, NULL); - renew_utg(ug, read_g, new_rtg_edges); rescue_contained_reads_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang, min_ovlp, 0, 10, 0, 1, NULL, NULL); - renew_utg(ug, read_g, 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(); + 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); renew_utg(ug, read_g, new_rtg_edges); - ///set_drop_trio_flag(*ug); n_vtx = read_g->n_seq; for (v = 0; v < n_vtx; v++) diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index 9ce909d..b976330 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -5,6 +5,7 @@ #include "Purge_Dups.h" #include "Overlaps.h" #include "Correct.h" +#include "kthread.h" #define Cal_Off(OFF) ((long long)((uint32_t)((OFF)>>32)) - (long long)((uint32_t)((OFF)))) #define Get_match(x) ((x).weight) @@ -98,7 +99,7 @@ void init_hap_alignment_struct(hap_alignment_struct* x, uint32_t size) kv_init(x->u_can.a); } -void destory_hap_alignment_struct(hap_alignment_struct* x, uint32_t size) +void destory_hap_alignment_struct(hap_alignment_struct* x) { free(x->vote_counting); free(x->visit); @@ -107,6 +108,61 @@ void destory_hap_alignment_struct(hap_alignment_struct* x, uint32_t size) kv_destroy(x->u_can.a); } + +typedef struct { + hap_alignment_struct* buf; + uint32_t num_threads; + + 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; +}hap_alignment_struct_pip; + +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) +{ + uint32_t i; + x->num_threads = num_threads; + x->buf = (hap_alignment_struct*)malloc(sizeof(hap_alignment_struct)*x->num_threads); + for (i = 0; i < x->num_threads; i++) + { + init_hap_alignment_struct(&(x->buf[i]), n_seq); + } + + x->ug = ug; + x->read_g = read_g; + x->reverse_sources = reverse_sources; + x->ruIndex = ruIndex; + x->coverage_cut = coverage_cut; + x->position_index = position_index; + x->Hap_rate = Hap_rate; + x->max_hang = max_hang; + x->min_ovlp = min_ovlp; + x->chain_rate = chain_rate; + x->all_ovlp = all_ovlp; +} + + +void destory_hap_alignment_struct_pip(hap_alignment_struct_pip* x) +{ + uint32_t i; + for (i = 0; i < x->num_threads; i++) + { + destory_hap_alignment_struct(&(x->buf[i])); + } + + free(x->buf); +} + void init_hap_overlaps_list(hap_overlaps_list* x, uint32_t num) { uint32_t i = 0; @@ -118,9 +174,9 @@ void init_hap_overlaps_list(hap_overlaps_list* x, uint32_t num) } } -void enable_debug_mode() +void enable_debug_mode(uint32_t mode) { - debug_enable = 1; + debug_enable = mode; } void destory_hap_overlaps_list(hap_overlaps_list* x) @@ -411,6 +467,7 @@ ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match, { max_count = 0; min_count = Len; + break; } qn = readIDs[i]>>33; if(reverse_sources[qn].length > 0) min_count++; @@ -446,6 +503,60 @@ ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match, (*Match) = max_count; (*Total) = min_count; } +/** +void get_pair_hap_similarity_deduplicate(uint64_t* readIDs, uint32_t Len, uint32_t target_uId, +ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match, double* Total) +{ + get_pair_hap_similarity(readIDs, Len, target_uId, reverse_sources, read_g, ruIndex, Match, Total); + return; + + #define CUTOFF_THRES 100 + uint32_t i, j, qn, tn, is_Unitig, uId, min_count = 0, max_count = 0, cutoff = 0, is_found; + for (i = 0; i < Len; i++) + { + if(cutoff > CUTOFF_THRES) + { + max_count = 0; + min_count = Len; + break; + } + qn = readIDs[i]>>33; + is_found = 0; + + for (j = 0; j < reverse_sources[qn].length; j++) + { + tn = Get_tn(reverse_sources[qn].buffer[j]); + if(read_g->seq[tn].del == 1) + { + get_R_to_U(ruIndex, tn, &tn, &is_Unitig); + if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue; + } + + + get_R_to_U(ruIndex, tn, &uId, &is_Unitig); + if(uId!=(uint32_t)-1 && is_Unitig == 1 && uId == target_uId) + { + max_count++; + } + min_count++; + is_found = 1; + } + + //means there is a match + if(is_found) + { + cutoff = 0; + } + else + { + cutoff++; + } + } + + (*Match) = max_count; + (*Total) = min_count; +} +**/ inline void check_hap_match(uint32_t qn, uint32_t targetBeg, uint32_t targetEnd, uint32_t targetID, uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, uint32_t* is_found, uint32_t* is_match) @@ -477,6 +588,42 @@ uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U } +/** +inline void check_hap_match_deduplicate(uint32_t qn, uint32_t targetBeg, uint32_t targetEnd, uint32_t targetID, +uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, uint32_t* is_found, uint32_t* is_match) +{ + check_hap_match(qn, targetBeg, targetEnd, targetID, position_index, reverse_sources, read_g, + ruIndex, is_found, is_match); + return; + + + uint32_t j, tn, uId, is_Unitig, offset; + (*is_found) = (*is_match) = 0; + + for (j = 0; j < reverse_sources[qn].length; j++) + { + tn = Get_tn(reverse_sources[qn].buffer[j]); + if(read_g->seq[tn].del == 1) + { + get_R_to_U(ruIndex, tn, &tn, &is_Unitig); + if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue; + } + + + get_R_to_U(ruIndex, tn, &uId, &is_Unitig); + if(uId!=(uint32_t)-1 && is_Unitig == 1 && uId == targetID) + { + offset = (uint32_t)(position_index[tn]); + if(offset >= targetBeg && offset <= targetEnd) + { + (*is_match)++; + } + } + (*is_found)++; + } + +} +**/ void determin_hap_alignment_boundary_single_side(uint64_t* readIDs, long long queryLen, long long targetBeg, long long targetEnd, long long targetID, long long eMatch, long long eTotal, long long dir, @@ -496,6 +643,7 @@ R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_co { check_hap_match(readIDs[i]>>33, targetBeg, targetEnd, targetID, position_index, reverse_sources, read_g, ruIndex, &is_found, &is_match); + min_count += is_found; max_count += is_match; if(max_count > min_count*Hap_rate) maxId = i; @@ -506,6 +654,8 @@ R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_co { check_hap_match(readIDs[i]>>33, targetBeg, targetEnd, targetID, position_index, reverse_sources, read_g, ruIndex, &is_found, &is_match); + + ///if(is_found > 0 && is_match > 0 && is_match > is_found*Hap_rate) if(is_found == 1 && is_match == 1) { break; @@ -522,6 +672,7 @@ R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_co { check_hap_match(readIDs[i]>>33, targetBeg, targetEnd, targetID, position_index, reverse_sources, read_g, ruIndex, &is_found, &is_match); + min_count += is_found; max_count += is_match; if(max_count > min_count*Hap_rate) maxId = i; @@ -531,6 +682,8 @@ R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_co { check_hap_match(readIDs[i]>>33, targetBeg, targetEnd, targetID, position_index, reverse_sources, read_g, ruIndex, &is_found, &is_match); + + ///if(is_found > 0 && is_match > 0 && is_match > is_found*Hap_rate) if(is_found == 1 && is_match == 1) { break; @@ -765,10 +918,10 @@ uint64_t* position_index, ma_utg_t* xReads, ma_utg_t* yReads) void get_base_boundary(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 begIndex, long long endIndex, uint32_t dir, uint32_t rev, -uint32_t* x_off, uint32_t* y_off) +uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, long long yEndIndex, +uint32_t dir, uint32_t rev, uint32_t* x_off, uint32_t* y_off) { - long long k, j; + long long k, j, offset; ma_hit_t_alloc *xR = NULL; ma_hit_t *h = NULL; ma_sub_t *sq = NULL, *st = NULL; @@ -779,7 +932,7 @@ uint32_t* x_off, uint32_t* y_off) (*x_off) = (*y_off) = (uint32_t)-1; if(dir == 1) { - for (k = endIndex; k >= begIndex; k--) + for (k = xEndIndex; k >= xBegIndex; k--) { xR = &(reverse_sources[xReads->a[k]>>33]); is_found = 0; oLen = 0; @@ -817,6 +970,11 @@ uint32_t* x_off, uint32_t* y_off) 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; @@ -835,7 +993,7 @@ uint32_t* x_off, uint32_t* y_off) } else { - for (k = begIndex; k <= endIndex; k++) + for (k = xBegIndex; k <= xEndIndex; k++) { xR = &(reverse_sources[xReads->a[k]>>33]); is_found = 0; oLen = 0; @@ -872,8 +1030,11 @@ uint32_t* x_off, uint32_t* y_off) 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; @@ -905,11 +1066,12 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon uint32_t x_pos_beg, y_pos_beg, x_pos_end, y_pos_end; /*************************x***************************/ get_base_boundary(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), 0, - Get_rev(*hap_can), &x_pos_beg, &y_pos_beg); + 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), 0, Get_rev(*hap_can), &x_pos_beg, &y_pos_beg); + get_base_boundary(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), 1, - Get_rev(*hap_can), &x_pos_end, &y_pos_end); + 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), 1, Get_rev(*hap_can), &x_pos_end, &y_pos_end); /*************************x***************************/ if(x_pos_beg == (uint32_t)-1 || y_pos_beg == (uint32_t)-1 @@ -975,6 +1137,7 @@ long long* r_y_pos_beg, long long* r_y_pos_end) ///flag = classify_hap_overlap(xBasePos, xBasePos, xReads->len, yBasePos, yBasePos, yReads->len); flag = vote_overlap_type(u_buffer, hap_can, position_index, xReads, yReads); + if(flag == XCY) { get_pair_hap_similarity(yReads->a, yLen, xUid, reverse_sources, read_g, ruIndex, @@ -1007,7 +1170,6 @@ long long* r_y_pos_beg, long long* r_y_pos_end) 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; @@ -1042,6 +1204,8 @@ long long* r_y_pos_beg, long long* r_y_pos_end) hap_can->index_end = determine_hap_overlap_type(hap_can, xReads, yReads, ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, min_ovlp, xUid, yUid, 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; @@ -1053,9 +1217,11 @@ long long* r_y_pos_beg, long long* r_y_pos_end) void print_hap_paf(ma_ug_t *ug, hap_overlaps* ovlp) { - fprintf(stderr, "utg%.6d%c\t%u\t%u\t%u\t%c\tutg%.6d%c\t%u\t%u\t%u\t%u\t%u\n", - ovlp->xUid+1, "lc"[ug->u.a[ovlp->xUid].circ], ug->u.a[ovlp->xUid].len, ovlp->x_beg_pos, ovlp->x_end_pos, "+-"[ovlp->rev], - ovlp->yUid+1, "lc"[ug->u.a[ovlp->yUid].circ], ug->u.a[ovlp->yUid].len, ovlp->y_beg_pos, ovlp->y_end_pos, ovlp->type, (uint32_t)ovlp->weight); + fprintf(stderr, "utg%.6d%c\t%u(%u)\t%u(%u)\t%u(%u)\t%c\tutg%.6d%c\t%u(%u)\t%u(%u)\t%u(%u)\t%u\t%u\n", + ovlp->xUid+1, "lc"[ug->u.a[ovlp->xUid].circ], ug->u.a[ovlp->xUid].len, ug->u.a[ovlp->xUid].n, + ovlp->x_beg_pos, ovlp->x_beg_id, ovlp->x_end_pos, ovlp->x_end_id, "+-"[ovlp->rev], + ovlp->yUid+1, "lc"[ug->u.a[ovlp->yUid].circ], ug->u.a[ovlp->yUid].len, ug->u.a[ovlp->yUid].n, + ovlp->y_beg_pos, ovlp->y_beg_id, ovlp->y_end_pos, ovlp->y_end_id, ovlp->type, (uint32_t)ovlp->weight); } void hap_alignment(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources, @@ -1139,18 +1305,291 @@ hap_overlaps_list* all_ovlp) xReads = &(ug->u.a[xUid]); yReads = &(ug->u.a[yUid]); u_buffer->a.n = 0; - - //if(xUid == 404 && yUid == 307) - // if(xReads->n <= 100 && yReads->n <= 100) - // { - // debug_enable = 1; - // } - // else - // { - // debug_enable = 0; - // } - // debug_enable = 1; + 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; + + qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment); + k = 0; + u_can->a.n = 0; + while (k < u_buffer->a.n) + { + hap_can.rev = u_buffer->a.a[k].x.el; + hap_can.index_beg = k; + hap_can.index_end = k; + hap_can.weight = u_buffer->a.a[k].weight; + hap_can.x_beg_pos = hap_can.x_end_pos = (uint32_t)(u_buffer->a.a[k].Off>>32); + hap_can.y_beg_pos = hap_can.y_end_pos = (uint32_t)(u_buffer->a.a[k].Off); + cur_offset = Cal_Off(u_buffer->a.a[k].Off); + interval_len = get_hap_overlapLen(hap_can.x_beg_pos, hap_can.x_end_pos, xReads->len, + hap_can.y_beg_pos, hap_can.y_end_pos, yReads->len, NULL, NULL, NULL, NULL); + + + k++; + while (k < u_buffer->a.n) + { + new_offset = Cal_Off(u_buffer->a.a[k].Off); + if(u_buffer->a.a[k].x.el != hap_can.rev) break; + if((new_offset - cur_offset)>(interval_len*chain_rate)) break; + + + hap_can.index_end = k; + hap_can.weight += u_buffer->a.a[k].weight; + + is_update = 0; + xPos = (uint32_t)(u_buffer->a.a[k].Off>>32); + yPos = (uint32_t)(u_buffer->a.a[k].Off); + if(xPos < hap_can.x_beg_pos) + { + hap_can.x_beg_pos = xPos; + is_update = 1; + } + + if(xPos > hap_can.x_end_pos) + { + hap_can.x_end_pos = xPos; + is_update = 1; + } + + if(yPos < hap_can.y_beg_pos) + { + hap_can.y_beg_pos = yPos; + is_update = 1; + } + + if(yPos > hap_can.y_end_pos) + { + hap_can.y_end_pos = yPos; + is_update = 1; + } + + if(new_offset == cur_offset) is_update = 0; + + if(is_update) + { + interval_len = get_hap_overlapLen(hap_can.x_beg_pos, hap_can.x_end_pos, xReads->len, + hap_can.y_beg_pos, hap_can.y_end_pos, yReads->len, NULL, NULL, NULL, NULL); + } + + k++; + } + + kv_push(hap_candidates, u_can->a, hap_can); + } + + 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(u_buffer, &(u_can->a.a[k]), position_index, xUid, yUid, + xReads, yReads, reverse_sources, read_g, ruIndex, coverage_cut, Hap_rate, max_hang, + min_ovlp, &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); + } + +} + + +static 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; + 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); + + + 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, xPos, yPos, is_update; + uint64_t tmp; + long long cur_offset, new_offset, interval_len; + 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++) { @@ -1344,7 +1783,6 @@ hap_overlaps_list* all_ovlp) } - int inline get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn) { uint32_t i; @@ -1587,10 +2025,10 @@ void clean_purge_graph(asg_t *purge_g, int max_dist, float drop_ratio) void get_node_boundary(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 begIndex, long long endIndex, uint32_t dir, uint32_t rev, -asg_arc_t* reture_t_f, asg_arc_t* reture_t_r) +uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, +long long yEndIndex, uint32_t dir, uint32_t rev, asg_arc_t* reture_t_f, asg_arc_t* reture_t_r) { - long long k, j; + long long k, j, offset; ma_hit_t_alloc *xR = NULL; ma_hit_t *h = NULL; ma_sub_t *sq = NULL, *st = NULL; @@ -1600,7 +2038,7 @@ asg_arc_t* reture_t_f, asg_arc_t* reture_t_r) reture_t_f->del = reture_t_r->del = 1; if(dir == 1) { - for (k = endIndex; k >= begIndex; k--) + for (k = xEndIndex; k >= xBegIndex; k--) { xR = &(reverse_sources[xReads->a[k]>>33]); is_found = 0; @@ -1637,6 +2075,12 @@ asg_arc_t* reture_t_f, asg_arc_t* reture_t_r) if(rev == 0 && v_dir != w_dir) continue; if(rev == 1 && v_dir == w_dir) continue; if(v_dir == 1) 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; @@ -1663,7 +2107,7 @@ asg_arc_t* reture_t_f, asg_arc_t* reture_t_r) } else { - for (k = begIndex; k <= endIndex; k++) + for (k = xBegIndex; k <= xEndIndex; k++) { xR = &(reverse_sources[xReads->a[k]>>33]); is_found = 0; oLen = 0; @@ -1701,6 +2145,11 @@ asg_arc_t* reture_t_f, asg_arc_t* reture_t_r) if(rev == 0 && v_dir != w_dir) continue; if(rev == 1 && v_dir == w_dir) continue; if(v_dir == 0) 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)); @@ -1900,7 +2349,8 @@ uint64_t* position_index, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, if(cut_end < endIndex) endIndex = cut_end; get_node_boundary(ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, - min_ovlp, xReads, yReads, v>>1, w>>1, begIndex, endIndex, v&1, x->rev, &t_forward, &t_backward); + 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); if(t_forward.del || t_backward.del) break; kv_push(asg_arc_t, edge->a, t_forward); @@ -2027,41 +2477,57 @@ uint64_t* position_index, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge, } free(b_0.b.a); } + +void print_all_purge_ovlp(ma_ug_t *ug, hap_overlaps_list* all_ovlp) +{ + uint32_t v, uId, i; + for (v = 0; v < all_ovlp->num; v++) + { + uId = v; + if(uId != 96 && uId != 272) continue; + for (i = 0; i < all_ovlp->x[uId].a.n; i++) + { + print_hap_paf(ug, &(all_ovlp->x[uId].a.a[i])); + } + } + +} void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t bi_graph_Len, uint32_t long_hap_overlap, float lable_match_rate, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio, uint32_t just_contain) { - fprintf(stderr, "*****************\n"); asg_t *purge_g = NULL; purge_g = asg_init(); - kvec_t_u64_warp u_vecs; - kv_init(u_vecs.a); asg_t* nsg = ug->g; uint32_t v, rId, uId, i, k, offset; ma_utg_t* reads = NULL; - uint8_t* visit = NULL; - visit = (uint8_t*)malloc(sizeof(uint8_t) * nsg->n_seq); - ///int flag; - memset(visit, 0, nsg->n_seq); + + // kvec_t_u64_warp u_vecs; + // kv_init(u_vecs.a); + // uint8_t* visit = NULL; + // visit = (uint8_t*)malloc(sizeof(uint8_t) * nsg->n_seq); + // memset(visit, 0, nsg->n_seq); + // uint64_t* vote_counting = (uint64_t*)malloc(sizeof(uint64_t)*nsg->n_seq); + // memset(vote_counting, 0, sizeof(uint64_t)*nsg->n_seq); + // kvec_asg_arc_t_offset u_buffer; + // kv_init(u_buffer.a); + // kvec_hap_candidates u_can; + // kv_init(u_can.a); uint64_t* position_index = (uint64_t*)malloc(sizeof(uint64_t)*read_g->n_seq); memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq); - uint64_t* vote_counting = (uint64_t*)malloc(sizeof(uint64_t)*nsg->n_seq); - memset(vote_counting, 0, sizeof(uint64_t)*nsg->n_seq); - kvec_asg_arc_t_offset u_buffer; - kv_init(u_buffer.a); - kvec_hap_candidates u_can; - kv_init(u_can.a); + hap_overlaps_list all_ovlp; init_hap_overlaps_list(&all_ovlp, nsg->n_seq); hap_overlaps_list back_all_ovlp; init_hap_overlaps_list(&back_all_ovlp, nsg->n_seq); ///uint32_t junk_cov, hap_cov, dip_cov, junk_occ, repeat_occ, single_cov; - ma_sub_t* purge_cut = NULL; - purge_cut = (ma_sub_t*)malloc(sizeof(ma_sub_t)*nsg->n_seq); asg_arc_t t; asg_arc_t* p = NULL; int r; + hap_alignment_struct_pip hap_buf; + + /****************************may have bugs********************************/ for (v = 0; v < nsg->n_seq; ++v) @@ -2099,24 +2565,29 @@ uint32_t just_contain) offset += (uint32_t)reads->a[i]; } - purge_cut[uId].del = 0; - purge_cut[uId].s = 0; - purge_cut[uId].e = offset; - purge_cut[uId].c = PRIMARY_LABLE; asg_seq_set(purge_g, uId, offset, 0); purge_g->seq[uId].c = PRIMARY_LABLE; } - for (v = 0; v < nsg->n_seq; v++) - { - uId = v; - if(nsg->seq[uId].del || nsg->seq[uId].c == ALTER_LABLE) continue; + init_hap_alignment_struct_pip(&hap_buf, asm_opt.thread_num, nsg->n_seq, ug, read_g, + reverse_sources, ruIndex, coverage_cut, position_index, density, max_hang, min_ovlp, + 0.05, &all_ovlp); - hap_alignment(ug, read_g, reverse_sources, ruIndex, coverage_cut, position_index, - vote_counting, visit, &u_vecs, &u_buffer, &u_can, uId, density, max_hang, min_ovlp, - 0.05, &all_ovlp); - } + kt_for(asm_opt.thread_num, hap_alignment_worker, &hap_buf, nsg->n_seq); + + ///if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp); + + + // for (v = 0; v < nsg->n_seq; v++) + // { + // uId = v; + // if(nsg->seq[uId].del || nsg->seq[uId].c == ALTER_LABLE) continue; + + // hap_alignment(ug, read_g, reverse_sources, ruIndex, coverage_cut, position_index, + // vote_counting, visit, &u_vecs, &u_buffer, &u_can, uId, density, max_hang, min_ovlp, + // 0.05, &all_ovlp); + // } normalize_hap_overlaps(&all_ovlp, &back_all_ovlp); ///debug_hap_overlaps(&all_ovlp, &back_all_ovlp); @@ -2186,9 +2657,12 @@ uint32_t just_contain) asg_symm(purge_g); clean_purge_graph(purge_g, bubble_dist, drop_ratio); - ///if(debug_enable) print_purge_gfa(ug, purge_g); + + // if(debug_enable) print_purge_gfa(ug, purge_g); + // 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, - max_hang, min_ovlp, edge, visit); + max_hang, min_ovlp, edge, hap_buf.buf[0].visit); } for (v = 0; v < all_ovlp.num; v++) @@ -2208,16 +2682,15 @@ uint32_t just_contain) } asg_cleanup(nsg); - kv_destroy(u_vecs.a); - kv_destroy(u_buffer.a); - kv_destroy(u_can.a); destory_hap_overlaps_list(&all_ovlp); destory_hap_overlaps_list(&back_all_ovlp); asg_destroy(purge_g); - free(purge_cut); free(position_index); - free(vote_counting); - free(visit); - fprintf(stderr, "#################\n"); + // kv_destroy(u_vecs.a); + // kv_destroy(u_buffer.a); + // kv_destroy(u_can.a); + // free(vote_counting); + // free(visit); + destory_hap_alignment_struct_pip(&hap_buf); } diff --git a/Purge_Dups.h b/Purge_Dups.h index e9e3152..07fa0d8 100644 --- a/Purge_Dups.h +++ b/Purge_Dups.h @@ -14,6 +14,6 @@ uint32_t just_contain); 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); -void enable_debug_mode(); +void enable_debug_mode(uint32_t mode); #endif \ No newline at end of file