From 1ead2d3081ff472c136f5898b670673cf5811fdb Mon Sep 17 00:00:00 2001 From: Haoyu Cheng Date: Mon, 23 Dec 2019 09:37:26 -0500 Subject: [PATCH] M50 improved --- Assembly.cpp | 130 +++++++ Correct.cpp | 817 +++++++++++++++++++++++++++++++--------- Correct.h | 4 + Hash_Table.h | 2 +- Overlaps.cpp | 944 ++++++++++++++++++++++++++++++++++++++++++++++- Overlaps.h | 7 + Process_Read.cpp | 21 ++ Process_Read.h | 2 +- 8 files changed, 1751 insertions(+), 176 deletions(-) diff --git a/Assembly.cpp b/Assembly.cpp index e6c56b3..2872496 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -3321,6 +3321,132 @@ void Output_PAF() } +int check_cluster(uint64_t* list, long long listLen, ma_hit_t_alloc* paf, float threshold) +{ + long long i, k; + uint32_t qn, tn; + long long T_edges, A_edges; + T_edges = A_edges = 0; + for (i = 0; i < listLen; i++) + { + qn = (uint32_t)list[i]; + for (k = i + 1; k < listLen; k++) + { + tn = (uint32_t)list[k]; + if(get_specific_overlap(&(paf[qn]), qn, tn) != -1) + { + A_edges++; + } + + if(get_specific_overlap(&(paf[tn]), tn, qn) != -1) + { + A_edges++; + } + + T_edges = T_edges + 2; + } + + } + + if(A_edges >= (T_edges*threshold)) + { + return 1; + } + else + { + return 0; + } +} + + + + +void rescue_edges(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, +long long readNum, long long rescue_threshold, float cluster_threshold) +{ + double startTime = Get_T(); + long long i, j, revises = 0; + uint32_t qn, tn; + + kvec_t(uint64_t) edge_vector; + kv_init(edge_vector); + kvec_t(uint64_t) edge_vector_index; + kv_init(edge_vector_index); + uint64_t flag; + int index; + + for (i = 0; i < readNum; i++) + { + edge_vector.n = 0; + edge_vector_index.n = 0; + for (j = 0; j < paf[i].length; j++) + { + qn = Get_qn(paf[i].buffer[j]); + tn = Get_tn(paf[i].buffer[j]); + index = get_specific_overlap(&(rev_paf[tn]), tn, qn); + if(index != -1) + { + flag = tn; + flag = flag << 32; + flag = flag | (uint64_t)(index); + kv_push(uint64_t, edge_vector, flag); + kv_push(uint64_t, edge_vector_index, j); + } + } + + ///the read itself has these overlaps, but all related reads do not have + ///we need to remove all overlaps from paf[i], and then add all overlaps to rev_paf[i] + if(edge_vector.n >= rescue_threshold && + check_cluster(edge_vector.a, edge_vector.n, paf, cluster_threshold) == 1) + { + fprintf(stderr,"\nremove following %d edges...\n", edge_vector.n); + print_revise_edges(&(paf[i]), edge_vector_index.a, edge_vector_index.n); + + add_overlaps(&(paf[i]), &(rev_paf[i]), edge_vector_index.a, edge_vector_index.n); + remove_overlaps(&(paf[i]), edge_vector_index.a, edge_vector_index.n); + revises = revises + edge_vector.n; + } + + edge_vector.n = 0; + edge_vector_index.n = 0; + for (j = 0; j < rev_paf[i].length; j++) + { + qn = Get_qn(rev_paf[i].buffer[j]); + tn = Get_tn(rev_paf[i].buffer[j]); + index = get_specific_overlap(&(paf[tn]), tn, qn); + if(index != -1) + { + flag = tn; + flag = flag << 32; + flag = flag | (uint64_t)(index); + kv_push(uint64_t, edge_vector, flag); + kv_push(uint64_t, edge_vector_index, j); + } + } + + ///the read itself do not have these overlaps, but all related reads have + ///we need to remove all overlaps from rev_paf[i], and then add all overlaps to paf[i] + if(edge_vector.n >= rescue_threshold && + check_cluster(edge_vector.a, edge_vector.n, paf, cluster_threshold) == 1) + { + fprintf(stderr,"\nadd following %d edges...\n", edge_vector.n); + print_revise_edges(&(rev_paf[i]), edge_vector_index.a, edge_vector_index.n); + + + remove_overlaps(&(rev_paf[i]), edge_vector_index.a, edge_vector_index.n); + add_overlaps_from_different_sources(paf, &(paf[i]), edge_vector.a, edge_vector.n); + revises = revises + edge_vector.n; + } + } + + + kv_destroy(edge_vector); + kv_destroy(edge_vector_index); + + fprintf(stderr, "[M::%s] took %0.2fs, revise edges #: %lld\n\n", __func__, Get_T()-startTime, revises); +} + + void generate_overlaps(int last_round) { @@ -3350,6 +3476,8 @@ void generate_overlaps(int last_round) pthread_join(_r_threads[i], NULL); free(_r_threads); + ///rescue_edges(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads, 4, 0.985); + fprintf(stdout, "Final overlaps have been calculated.\n"); fprintf(stdout, "%-30s%18.2f\n\n", "Final overlaps calculation time:", Get_T() - start_time); @@ -3369,6 +3497,8 @@ void generate_overlaps(int last_round) + + void Correct_Reads(int last_round) { diff --git a/Correct.cpp b/Correct.cpp index d4bb5bc..148b73b 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -601,6 +601,65 @@ int* r_extra_begin, int* r_extra_end, long long* r_y_start, long long* r_y_lengt } +int verify_single_window(long long x_start, long long x_end, +long long overlap_x_s, long long overlap_y_s, int x_id, +int y_id, int x_strand, char* x_buffer, char* y_buffer, +All_reads* R_INF) +{ + char* x_string = NULL; + char* y_string = NULL; + int extra_begin, extra_end, x_len, threshold; + long long y_start; + long long Window_Len, o_len; + int end_site; + unsigned int error; + + + + x_len = x_end - x_start + 1; + threshold = x_len * THRESHOLD_RATE; + /****************************may have bugs********************************/ + threshold = Adjust_Threshold(threshold, x_len); + /****************************may have bugs********************************/ + + y_start = (x_start - overlap_x_s) + overlap_y_s; + + Window_Len = x_len + (threshold << 1); + + if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF, + &extra_begin, &extra_end, &y_start, &o_len)) + { + return 0; + } + + ///use unusual direction here + /** + fill_subregion(y_buffer, y_start, o_len, y_strand, R_INF, y_id, extra_begin, extra_end); + ///x is always forward strand + recover_UC_Read_sub_region(x_buffer, x_start, x_len, 0, R_INF, x_id); + **/ + ///use unusual direction here, here y is always forward strand + fill_subregion(y_buffer, y_start, o_len, 0, R_INF, y_id, extra_begin, extra_end); + recover_UC_Read_sub_region(x_buffer, x_start, x_len, x_strand, R_INF, x_id); + + + x_string = x_buffer; + y_string = y_buffer; + + end_site = Reserve_Banded_BPM(y_string, Window_Len, x_string, x_len, threshold, &error); + + // if(y_id == 439960 && (x_id == 439950 || x_id == 5598234)) + // { + // fprintf(stderr, "y_id: %d, x_id: %d, error: %llu", y_id, x_id, error); + // } + + if (error!=(unsigned int)-1) + { + return 1; + } + + return 0; +} void verify_window(long long window_start, long long window_end, overlap_region_alloc* overlap_list,Correct_dumy* dumy, All_reads* R_INF, char* r_string) @@ -4227,7 +4286,8 @@ UC_Read* g_read) int scanLen = 10; long long boundaryLen = 200; long long single_sideLen = boundaryLen/2; - long long useless_side = single_sideLen/2; + long long force_useless_side = single_sideLen/2; + long long L_useless_side, R_useless_side; int alpha = 1; long long y_start, x_start; long long y_end, x_end; @@ -4242,22 +4302,32 @@ UC_Read* g_read) unsigned int error; int real_y_start; sub_list->boundary_cigars.length = sub_list->w_list_length - 1; + ///the (i)-th boundary between the (i)-th window and the (i+1)-th window + ///that means it includes (the tail of (i)-th window) and (the header of (i+1)-th window) + ///note the (i)-th boundary is calculated at the (i)-th window for (i = 0; i + 1 < sub_list->w_list_length; i++) { - ///if two windows are not aligned - if(sub_list->w_list[i].y_end == -1 && sub_list->w_list[i+1].y_end == -1) + ///if both of the two windows are not aligned + ///it is not necessary to calculate the boundary + ///if(sub_list->w_list[i].y_end == -1 && sub_list->w_list[i+1].y_end == -1) + if(sub_list->w_list[i].y_end == -1 || sub_list->w_list[i+1].y_end == -1) { sub_list->boundary_cigars.buffer[i].error = -1; sub_list->boundary_cigars.buffer[i].y_end = -1; continue; } + ///note if w_list[i+1].y_start or sub_list->w_list[i].y_end is -1 + ///y_distance might have some problems at the last of this function + ///we need to deal with it carefully y_distance = sub_list->w_list[i+1].y_start - sub_list->w_list[i].y_end - 1; ///if two windows are aligned if(sub_list->w_list[i].y_end != -1 && sub_list->w_list[i+1].y_end != -1 && y_distance == 0) { + ///scan backward scan_cigar(&(sub_list->w_list[i].cigar), &b_err, scanLen, 1); + ///scan forward scan_cigar(&(sub_list->w_list[i+1].cigar), &f_err, scanLen, 0); if(b_err == 0 && f_err == 0) { @@ -4273,12 +4343,12 @@ UC_Read* g_read) { y_start = sub_list->w_list[i].y_end; x_start = sub_list->w_list[i].x_end; - } + }///if the (i)-th window is not matched, have a look at the (i+1)-th window else if(sub_list->w_list[i+1].y_end != -1) { y_start = sub_list->w_list[i+1].y_start; x_start = sub_list->w_list[i+1].x_start; - } + }///if both of these two windows are not matched, directly skip else { sub_list->boundary_cigars.buffer[i].error = -1; @@ -4286,30 +4356,34 @@ UC_Read* g_read) continue; } - + ///it seems we don't need to record x_start and y_start sub_list->boundary_cigars.buffer[i].extra_begin = x_start; sub_list->boundary_cigars.buffer[i].extra_end = y_start; ///leftLen and rightLen are used for x ///x should be at [sub_list->w_list[i].x_start, sub_list->w_list[i+1].x_end] ///y shouldn't have limitation + ///note that the x_start and x_end should not be -1 in any case + ///up to now, x_start and y_start are not -1 + ///leftLen does not include x_start itself, rightLen does + ///gnerally speaking, rightLen should be always larger than leftLen leftLen = MIN(MIN((x_start - sub_list->w_list[i].x_start), y_start), single_sideLen); rightLen = MIN(MIN((sub_list->w_list[i+1].x_end + 1 - x_start), y_readLen - y_start), single_sideLen); - + ///xLen should be the sum length of two windows xLen = leftLen + rightLen; x_start = x_start - leftLen; x_end = x_start + xLen - 1; y_start = y_start - leftLen; - ///leftLen, rightLen - if(leftLen <= useless_side || rightLen <= useless_side) - { - sub_list->boundary_cigars.buffer[i].error = -1; - sub_list->boundary_cigars.buffer[i].y_end = -1; - continue; - } + ///if we don't have enough leftLen and rightLen + // if(leftLen <= useless_side || rightLen <= useless_side) + // { + // sub_list->boundary_cigars.buffer[i].error = -1; + // sub_list->boundary_cigars.buffer[i].y_end = -1; + // continue; + // } threshold = xLen * THRESHOLD_RATE; @@ -4341,6 +4415,7 @@ UC_Read* g_read) end_site = Reserve_Banded_BPM_PATH(y_string, yLen, x_string, xLen, threshold, &error, &real_y_start, &(dumy->path_length), dumy->matrix_bit, dumy->path, -1, -1); + ///means this window is matched if (error!=(unsigned int)-1) { sub_list->boundary_cigars.buffer[i].x_start = x_start; @@ -4348,23 +4423,51 @@ UC_Read* g_read) generate_cigar(dumy->path, dumy->path_length, &(sub_list->boundary_cigars.buffer[i]), &real_y_start, &end_site, &error, x_string, xLen, y_string); - + ///should not adjust cigar here, adjust cigar may cause problem + ///that is not what we want ///y_distance can be less than 0, or larger than 0 - if(y_distance<0) y_distance = y_distance * (-1); + ///please if one of the two windows is not matched, + ///y_distance may have potential problems + if(y_distance < 0) y_distance = y_distance * (-1); ///leftLen, rightLen - if(leftLen <= useless_side || rightLen <= useless_side) + // if(leftLen <= useless_side || rightLen <= useless_side) + // { + // sub_list->boundary_cigars.buffer[i].error = -1; + // sub_list->boundary_cigars.buffer[i].y_end = -1; + // continue; + // } + L_useless_side = R_useless_side = force_useless_side; + + ///first window + if((i == 0) && (x_start == sub_list->w_list[0].x_start)) + { + L_useless_side = 0; + } + ///last window + if((i == sub_list->w_list_length - 2) && + (x_end == sub_list->w_list[sub_list->w_list_length - 1].x_end)) + { + R_useless_side = 0; + } + + if(leftLen <= L_useless_side || rightLen <= R_useless_side) { sub_list->boundary_cigars.buffer[i].error = -1; sub_list->boundary_cigars.buffer[i].y_end = -1; continue; } - + + + ///up to now, if we require (i)-th window and (i+1)-th window are matched + ///boundary_cigars.buffer[i].cigar, w_list[i].cigar and w_list[i+1].cigar are avaiable + ///get the error excluding the first and the last useless_side bases scan_cigar_interval(&(sub_list->boundary_cigars.buffer[i].cigar), &m_error, - useless_side, xLen-useless_side-1); - scan_cigar(&(sub_list->w_list[i].cigar), &b_err, leftLen-useless_side, 1); - scan_cigar(&(sub_list->w_list[i+1].cigar), &f_err, rightLen-useless_side, 0); + L_useless_side, xLen-R_useless_side-1); + scan_cigar(&(sub_list->w_list[i].cigar), &b_err, leftLen-L_useless_side, 1); + scan_cigar(&(sub_list->w_list[i+1].cigar), &f_err, rightLen-R_useless_side, 0); + if(f_err + b_err + y_distance + alpha < m_error) { sub_list->boundary_cigars.buffer[i].error = -1; @@ -4378,9 +4481,9 @@ UC_Read* g_read) sub_list->boundary_cigars.buffer[i].x_start = x_start; sub_list->boundary_cigars.buffer[i].x_end = x_end; - sub_list->boundary_cigars.buffer[i].error_threshold = useless_side; - - + ///sub_list->boundary_cigars.buffer[i].error_threshold = useless_side; + sub_list->boundary_cigars.buffer[i].extra_begin = L_useless_side; + sub_list->boundary_cigars.buffer[i].extra_end = R_useless_side; } else { @@ -4393,6 +4496,104 @@ UC_Read* g_read) } + +void debug_window_cigar(overlap_region_alloc* overlap_list, UC_Read* g_read, Correct_dumy* dumy, +All_reads* R_INF, int test_window, int test_boundary) +{ + long long i, j, y_id, y_strand; + char* x_string; + char* y_string; + long long x_start; + long long x_end; + long long x_len; + long long y_start; + long long y_end; + long long y_len; + + for (j = 0; j < overlap_list->length; j++) + { + y_id = overlap_list->list[j].y_id; + y_strand = overlap_list->list[j].y_pos_strand; + if(overlap_list->list[j].is_match == 1) + { + + if(test_window == 1) + { + for (i = 0; i < overlap_list->list[j].w_list_length; i++) + { + if(overlap_list->list[j].w_list[i].y_end != -1) + { + ///there is no problem for x + x_start = overlap_list->list[j].w_list[i].x_start; + x_end = overlap_list->list[j].w_list[i].x_end; + x_len = x_end - x_start + 1; + + x_string = g_read->seq + x_start; + + y_start = overlap_list->list[j].w_list[i].y_start; + y_end = overlap_list->list[j].w_list[i].y_end; + y_len = y_end - y_start + 1; + + recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, R_INF, y_id); + y_string = dumy->overlap_region; + + + if(verify_cigar(x_string, x_len, y_string, y_len, &overlap_list->list[j].w_list[i].cigar, + overlap_list->list[j].w_list[i].error)) + { + fprintf(stderr, "j: %d, i: %d, y_id: %d, y_start: %d, y_end: %d\n", j, i, y_id, y_start, y_end); + } + } + } + } + + + + if(test_boundary == 1) + { + for (i = 0; i < overlap_list->list[j].boundary_cigars.length; i++) + { + if(overlap_list->list[j].boundary_cigars.buffer[i].y_end != -1) + { + x_start = overlap_list->list[j].boundary_cigars.buffer[i].x_start; + x_end = overlap_list->list[j].boundary_cigars.buffer[i].x_end; + x_len = x_end - x_start + 1; + x_string = g_read->seq + x_start; + + y_start = overlap_list->list[j].boundary_cigars.buffer[i].y_start; + y_end = overlap_list->list[j].boundary_cigars.buffer[i].y_end; + y_len = y_end - y_start + 1; + + recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, + R_INF, y_id); + y_string = dumy->overlap_region; + + + if(verify_cigar(x_string, x_len, y_string, y_len, + &overlap_list->list[j].boundary_cigars.buffer[i].cigar, + overlap_list->list[j].boundary_cigars.buffer[i].error)) + { + fprintf(stderr, "j: %d, i: %d, y_id: %d, y_start: %d, y_end: %d\n", j, i, y_id, y_start, y_end); + } + } + } + } + + + if(test_window == 1 && test_boundary == 1) + { + if(overlap_list->list[j].w_list_length != + overlap_list->list[j].boundary_cigars.length + 1) + { + fprintf(stderr, "w_list_length: %d, boundary_cigars.length: %d\n", + overlap_list->list[j].w_list_length, overlap_list->list[j].boundary_cigars.length); + } + } + + } + } +} + inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read) { @@ -4821,9 +5022,10 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea ///error_rate = trim_error_rate(overlap_list, j); error_rate = non_trim_error_rate(overlap_list, j, R_INF, dumy, g_read); - + ///if(error_rate <= 0.015) - if(error_rate <= 0.03) + ///if(error_rate <= 0.03) + if(error_rate <= FINAL_OVERLAP_ERROR_RATE) { overlap_list->mapped_overlaps_length += overlap_length; overlap_list->list[j].is_match = 1; @@ -4836,74 +5038,7 @@ inline void recalcate_window_advance(overlap_region_alloc* overlap_list, All_rea } } - - /** - for (j = 0; j < overlap_list->length; j++) - { - y_id = overlap_list->list[j].y_id; - y_strand = overlap_list->list[j].y_pos_strand; - if(overlap_list->list[j].is_match == 1) - { - for (i = 0; i < overlap_list->list[j].w_list_length; i++) - { - if(overlap_list->list[j].w_list[i].y_end != -1) - { - ///there is no problem for x - x_start = overlap_list->list[j].w_list[i].x_start; - x_end = overlap_list->list[j].w_list[i].x_end; - x_len = x_end - x_start + 1; - - x_string = g_read->seq + x_start; - - y_start = overlap_list->list[j].w_list[i].y_start; - y_end = overlap_list->list[j].w_list[i].y_end; - y_len = y_end - y_start + 1; - - recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, R_INF, y_id); - y_string = dumy->overlap_region; - - - if(verify_cigar(x_string, x_len, y_string, y_len, &overlap_list->list[j].w_list[i].cigar, - overlap_list->list[j].w_list[i].error)) - { - fprintf(stderr, "j: %d, i: %d, y_id: %d, y_start: %d, y_end: %d\n", j, i, y_id, y_start, y_end); - } - } - - } - - - - for (i = 0; i < overlap_list->list[j].boundary_cigars.length; i++) - { - if(overlap_list->list[j].boundary_cigars.buffer[i].y_end != -1) - { - x_start = overlap_list->list[j].boundary_cigars.buffer[i].x_start; - x_end = overlap_list->list[j].boundary_cigars.buffer[i].x_end; - x_len = x_end - x_start + 1; - x_string = g_read->seq + x_start; - - y_start = overlap_list->list[j].boundary_cigars.buffer[i].y_start; - y_end = overlap_list->list[j].boundary_cigars.buffer[i].y_end; - y_len = y_end - y_start + 1; - - recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_len, y_strand, - R_INF, y_id); - y_string = dumy->overlap_region; - - - if(verify_cigar(x_string, x_len, y_string, y_len, - &overlap_list->list[j].boundary_cigars.buffer[i].cigar, - overlap_list->list[j].boundary_cigars.buffer[i].error)) - { - fprintf(stderr, "j: %d, i: %d, y_id: %d, y_start: %d, y_end: %d\n", j, i, y_id, y_start, y_end); - } - } - } - } - } - **/ - + ///debug_window_cigar(overlap_list, g_read, dumy, R_INF, 1, 1); } @@ -8696,11 +8831,12 @@ CIGAR* cigar, haplotype_evdience_alloc* hap) } -///[xBeg, xEnd] +///mark SNPs at [xBeg, xEnd], note we need to deal with flag_offset carefully void markSNP_detail(CIGAR* cigar_record, uint8_t* flag, long long xBeg, long long xEnd, long long flag_offset) { if(xBeg > xEnd) return; + int operation; int operationLen; long long x_i, y_i, cigar_i, i; @@ -8715,7 +8851,6 @@ long long xBeg, long long xEnd, long long flag_offset) break; } - ///match if (operation == 0) { @@ -8726,7 +8861,7 @@ long long xBeg, long long xEnd, long long flag_offset) { for (i = 0; i < operationLen; i++) { - ///flag_offset + /// note we need to deal with flag_offset carefully if(flag[x_i - flag_offset] < 127 && x_i >= xBeg && x_i <= xEnd) { flag[x_i - flag_offset]++; @@ -8735,7 +8870,7 @@ long long xBeg, long long xEnd, long long flag_offset) x_i++; y_i++; } - }///insertion + }///insertion, that means y has more bases than x else if (operation == 2) { y_i += operationLen; @@ -8749,38 +8884,54 @@ long long xBeg, long long xEnd, long long flag_offset) } } + +///window_offset is still the x-based offset +///x_total_start and y_total_start are global positions, instead of local positions void markSNP_advance( long long window_offset, long long x_total_start, long long x_length, long long y_total_start, long long y_length, -window_list* current_cigar, window_list* beg_cigar, window_list* end_cigar, haplotype_evdience_alloc* hap) +window_list* current_cigar, window_list* beg_cigar, window_list* end_cigar, +haplotype_evdience_alloc* hap) { long long x_total_end = x_total_start + x_length - 1; ///mismatches based on the offset of x long long inner_offset = x_total_start - window_offset; CIGAR* cigar_record; - long long useless_side, xleftLen, xrightLen, x_interval_beg, x_interval_end; + ///long long useless_side, xleftLen, xrightLen, x_interval_beg, x_interval_end; + long long L_useless_side, R_useless_side, xleftLen, xrightLen, x_interval_beg, x_interval_end; long long current_cigar_beg, current_cigar_end; ///for current_cigar, [current_cigar_beg, current_cigar_end] current_cigar_beg = 0; current_cigar_end = x_length - 1; - - if(beg_cigar!=NULL && beg_cigar->y_end!=-1) + ///if the beg_cigar is available + if(beg_cigar != NULL && beg_cigar->y_end!=-1) { - useless_side = beg_cigar->error_threshold; + ///useless_side = beg_cigar->error_threshold; + L_useless_side = beg_cigar->extra_begin; + R_useless_side = beg_cigar->extra_end; cigar_record = &(beg_cigar->cigar); - + ///again, xleftLen does not include x_total_start itself, but includes beg_cigar->x_start xleftLen = x_total_start - beg_cigar->x_start; + ///xrightLen includes both x_total_start and beg_cigar->x_end xrightLen = beg_cigar->x_end - x_total_start + 1; ///actually xleftLen could be no larger than useless_side - ///but such window has already been filtered out in calculate_boundary_cigars - if(xleftLen > useless_side && xrightLen > useless_side) + ///but such window has already been filtered out at calculate_boundary_cigars + ///if(xleftLen > useless_side && xrightLen > useless_side) + if(xleftLen > L_useless_side && xrightLen > R_useless_side) { + ///[x_interval_beg, x_interval_end] are the offsets to beg_cigar->x_start + ///they are local postions, instead of global positions + x_interval_beg = xleftLen; - x_interval_end = x_interval_beg + (xrightLen - useless_side) - 1; - ///this is the offset of the current cigar - current_cigar_beg = xrightLen - useless_side; + //x_interval_end = x_interval_beg + (xrightLen - useless_side) - 1; + x_interval_end = x_interval_beg + (xrightLen - R_useless_side) - 1; + ///current_cigar_beg is the offset of the current cigar + ///that is the beg of current_cigar_beg + ///current_cigar_beg = xrightLen - useless_side; + current_cigar_beg = xrightLen - R_useless_side; + markSNP_detail(cigar_record, hap->flag + inner_offset, x_interval_beg, x_interval_end, x_interval_beg); } @@ -8788,31 +8939,53 @@ window_list* current_cigar, window_list* beg_cigar, window_list* end_cigar, hapl if(end_cigar!=NULL && end_cigar->y_end!=-1) { - useless_side = end_cigar->error_threshold; + ///useless_side = end_cigar->error_threshold; + L_useless_side = end_cigar->extra_begin; + R_useless_side = end_cigar->extra_end; cigar_record = &(end_cigar->cigar); + ///again, xleftLen does not include x_total_end, but includes end_cigar->x_start + ///it seems to be not what we want xleftLen = x_total_end - end_cigar->x_start; + ///xrightLen includes both x_total_end and end_cigar->x_end + ///it is also not what we want xrightLen = end_cigar->x_end - x_total_end + 1; + ///we hope that x_total_end should be included in xleftLen, instead of xrightLen + ///that means xleftLen should + 1, while xrightLen should -1 + ///but it is fine here ///actually xrightLen could be no larger than useless_side ///but such window has already been filtered out in calculate_boundary_cigars - if(xleftLen > useless_side && xrightLen > useless_side) + ///if(xleftLen > useless_side && xrightLen > useless_side) + if(xleftLen > L_useless_side && xrightLen > R_useless_side) { + ///[x_interval_beg, x_interval_end] are the offsets to beg_cigar->x_start + ///they are local postions, instead of global positions x_interval_end = xleftLen; - x_interval_beg = x_interval_end - (xleftLen + 1 - useless_side) + 1; - current_cigar_end = (x_length - 1) - (xleftLen + 1 - useless_side); + ///the real left part length is (xleftLen + 1) + ///so the useful left part length is ((xleftLen + 1) - useless_side) + ///x_interval_beg = x_interval_end - ((xleftLen + 1) - useless_side) + 1; + x_interval_beg = x_interval_end - ((xleftLen + 1) - L_useless_side) + 1; + ///current_cigar_end = (x_length - 1) - ((xleftLen + 1) - useless_side); + current_cigar_end = (x_length - 1) - ((xleftLen + 1) - L_useless_side); + markSNP_detail(cigar_record, hap->flag + end_cigar->x_start - window_offset, x_interval_beg, x_interval_end, 0); } } - - - - // markSNP_detail(&(current_cigar->cigar), hap->flag + inner_offset, 0, x_length/2); - // markSNP_detail(&(current_cigar->cigar), hap->flag + inner_offset, x_length/2+1, x_length - 1); - ///markSNP_detail(&(current_cigar->cigar), hap->flag + inner_offset, 0, x_length - 1, 0); markSNP_detail(&(current_cigar->cigar), hap->flag + inner_offset, current_cigar_beg, current_cigar_end, 0); + + // end_mark: + // if(x_length >= 4) + // { + // markSNP_detail(&(current_cigar->cigar), hap->flag + inner_offset, 0, x_length/2, 0); + // markSNP_detail(&(current_cigar->cigar), hap->flag + inner_offset, x_length/2+1, x_length - 1, 0); + // } + // else + // { + // markSNP_detail(&(current_cigar->cigar), hap->flag + inner_offset, 0, x_length - 1, 0); + // } } @@ -8922,6 +9095,239 @@ CIGAR* cigar, haplotype_evdience_alloc* hap, int snp_threshold) +///mark SNPs at [xBeg, xEnd], note we need to deal with flag_offset carefully +void addSNPtohaplotype_details(CIGAR* cigar_record, uint8_t* flag, +char* x_string, char* y_string, long long x_total_start, long long y_total_start, +long long xBeg, long long xEnd, int overlapID, long long flag_offset, +haplotype_evdience_alloc* hap, long long snp_threshold) +{ + if(xBeg > xEnd) return; + int operation; + int operationLen; + long long x_i, y_i, cigar_i, i; + i = cigar_i = x_i = y_i = 0; + + haplotype_evdience ev; + + ///note that node 0 is the start node + ///0 is match, 1 is mismatch, 2 is up, 3 is left + ///2 represents thre are more bases at y + ///3 represents thre are more bases at x + while (cigar_i < cigar_record->length) + { + operation = cigar_record->C_C[cigar_i]; + operationLen = cigar_record->C_L[cigar_i]; + if(x_i > xEnd) + { + break; + } + + ///matches + if (operation == 0) + { + for (i = 0; i < operationLen; i++) + { + ///should be at least 2 mismatches + /// note we need to deal with flag_offset carefully + if(flag[x_i - flag_offset] > snp_threshold && x_i >= xBeg && x_i <= xEnd) + { + ev.misBase = y_string[y_i]; + ev.overlapID = overlapID; + ev.site = x_total_start + x_i; + ev.overlapSite = y_total_start + y_i; + ev.type = 0; + addHaplotypeEvdience(hap, &ev); + } + ///inner_offset++; + x_i++; + y_i++; + } + } + else if(operation == 1) + { + for (i = 0; i < operationLen; i++) + { + + /// should be at least 2 mismatches + /// note we need to deal with flag_offset carefully + if(flag[x_i - flag_offset] > snp_threshold && x_i >= xBeg && x_i <= xEnd) + { + ev.misBase = y_string[y_i]; + ev.overlapID = overlapID; + ev.site = x_total_start + x_i; + ev.overlapSite = y_total_start + y_i; + ev.type = 1; + addHaplotypeEvdience(hap, &ev); + } + + ///inner_offset++; + x_i++; + y_i++; + } + }///insertion, 2 represents thre are more bases at y + else if (operation == 2) + { + y_i += operationLen; + }///3 represents thre are more bases at x + else if (operation == 3) + { + /****************************may have bugs********************************/ + for (i = 0; i < operationLen; i++) + { + ///if(hap->flag[inner_offset] > snp_threshold) + /// should be at least 2 mismatches + /// note we need to deal with flag_offset carefully + if(flag[x_i - flag_offset] > snp_threshold && x_i >= xBeg && x_i <= xEnd) + { + ev.misBase = 'N'; + ev.overlapID = overlapID; + ev.site = x_total_start + x_i; + ev.overlapSite = y_total_start + y_i; + ev.type = 2; + addHaplotypeEvdience(hap, &ev); + } + + ///inner_offset++; + x_i++; + } + /****************************may have bugs********************************/ + } + cigar_i++; + } +} + + +void addSNPtohaplotype_advance( +long long window_offset, int overlapID, +long long x_total_start, long long x_length, +long long y_total_start, long long y_length, +window_list* current_cigar, window_list* beg_cigar, window_list* end_cigar, +haplotype_evdience_alloc* hap, int snp_threshold, char* x_T_string, char* y_T_string) +{ + long long x_total_end = x_total_start + x_length - 1; + long long inner_offset = x_total_start - window_offset; + CIGAR* cigar_record; + ///long long useless_side, xleftLen, xrightLen, x_interval_beg, x_interval_end; + long long L_useless_side, R_useless_side, xleftLen, xrightLen, x_interval_beg, x_interval_end; + long long current_cigar_beg, current_cigar_end; + + ///for current_cigar, [current_cigar_beg, current_cigar_end] + current_cigar_beg = 0; + current_cigar_end = x_length - 1; + ///if the beg_cigar is available + if(beg_cigar != NULL && beg_cigar->y_end!=-1) + { + ///useless_side = beg_cigar->error_threshold; + L_useless_side = beg_cigar->extra_begin; + R_useless_side = beg_cigar->extra_end; + cigar_record = &(beg_cigar->cigar); + ///again, xleftLen does not include x_total_start itself, but includes beg_cigar->x_start + xleftLen = x_total_start - beg_cigar->x_start; + ///xrightLen includes both x_total_start and beg_cigar->x_end + xrightLen = beg_cigar->x_end - x_total_start + 1; + ///actually xleftLen could be no larger than useless_side + ///but such window has already been filtered out at calculate_boundary_cigars + ///if(xleftLen > useless_side && xrightLen > useless_side) + if(xleftLen > L_useless_side && xrightLen > R_useless_side) + { + ///[x_interval_beg, x_interval_end] are the offsets to beg_cigar->x_start + ///they are local postions, instead of global positions + + x_interval_beg = xleftLen; + ///x_interval_end = x_interval_beg + (xrightLen - useless_side) - 1; + x_interval_end = x_interval_beg + (xrightLen - R_useless_side) - 1; + ///current_cigar_beg is the offset of the current cigar + ///that is the beg of current_cigar_beg + ///current_cigar_beg = xrightLen - useless_side; + current_cigar_beg = xrightLen - R_useless_side; + + // markSNP_detail(cigar_record, hap->flag + inner_offset, x_interval_beg, + // x_interval_end, x_interval_beg); + addSNPtohaplotype_details(cigar_record, hap->flag + inner_offset, + x_T_string + beg_cigar->x_start, y_T_string + beg_cigar->y_start, + beg_cigar->x_start, beg_cigar->y_start, x_interval_beg, x_interval_end, + overlapID, x_interval_beg, hap, snp_threshold); + } + } + + + if(end_cigar!=NULL && end_cigar->y_end!=-1) + { + ///useless_side = end_cigar->error_threshold; + L_useless_side = end_cigar->extra_begin; + R_useless_side = end_cigar->extra_end; + cigar_record = &(end_cigar->cigar); + ///again, xleftLen does not include x_total_end, but includes end_cigar->x_start + ///it seems to be not what we want + xleftLen = x_total_end - end_cigar->x_start; + ///xrightLen includes both x_total_end and end_cigar->x_end + ///it is also not what we want + xrightLen = end_cigar->x_end - x_total_end + 1; + ///we hope that x_total_end should be included in xleftLen, instead of xrightLen + ///that means xleftLen should + 1, while xrightLen should -1 + ///but it is fine here + + ///actually xrightLen could be no larger than useless_side + ///but such window has already been filtered out in calculate_boundary_cigars + ///if(xleftLen > useless_side && xrightLen > useless_side) + if(xleftLen > L_useless_side && xrightLen > R_useless_side) + { + ///[x_interval_beg, x_interval_end] are the offsets to beg_cigar->x_start + ///they are local postions, instead of global positions + x_interval_end = xleftLen; + ///the real left part length is (xleftLen + 1) + ///so the useful left part length is ((xleftLen + 1) - useless_side) + ///x_interval_beg = x_interval_end - ((xleftLen + 1) - useless_side) + 1; + x_interval_beg = x_interval_end - ((xleftLen + 1) - L_useless_side) + 1; + ///current_cigar_end = (x_length - 1) - ((xleftLen + 1) - useless_side); + current_cigar_end = (x_length - 1) - ((xleftLen + 1) - L_useless_side); + + // markSNP_detail(cigar_record, hap->flag + end_cigar->x_start - window_offset, + // x_interval_beg, x_interval_end, 0); + addSNPtohaplotype_details(cigar_record, hap->flag + end_cigar->x_start - window_offset, + x_T_string + end_cigar->x_start, y_T_string + end_cigar->y_start, + end_cigar->x_start, end_cigar->y_start, x_interval_beg, x_interval_end, + overlapID, 0, hap, snp_threshold); + } + } + + // markSNP_detail(&(current_cigar->cigar), hap->flag + inner_offset, current_cigar_beg, + // current_cigar_end, 0); + addSNPtohaplotype_details(&(current_cigar->cigar), hap->flag + inner_offset, + x_T_string + current_cigar->x_start, y_T_string + current_cigar->y_start, + current_cigar->x_start, current_cigar->y_start, current_cigar_beg, + current_cigar_end, overlapID, 0, hap, snp_threshold); + + /** + end_add: + if(x_length >= 4) + { + addSNPtohaplotype_details(&(current_cigar->cigar), hap->flag + inner_offset, + x_T_string + current_cigar->x_start, y_T_string + current_cigar->y_start, + current_cigar->x_start, current_cigar->y_start, 0, + x_length/2, overlapID, 0, hap, snp_threshold); + + addSNPtohaplotype_details(&(current_cigar->cigar), hap->flag + inner_offset, + x_T_string + current_cigar->x_start, y_T_string + current_cigar->y_start, + current_cigar->x_start, current_cigar->y_start, x_length/2 + 1, + x_length - 1, overlapID, 0, hap, snp_threshold); + } + else + { + addSNPtohaplotype_details(&(current_cigar->cigar), hap->flag + inner_offset, + x_T_string + current_cigar->x_start, y_T_string + current_cigar->y_start, + current_cigar->x_start, current_cigar->y_start, 0, + x_length - 1, overlapID, 0, hap, snp_threshold); + } + **/ +} + + + + + + + @@ -9063,17 +9469,17 @@ window_list** end_cigar) void cluster_advance(char* r_string, long long window_start, long long window_end, overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, -haplotype_evdience_alloc* hap) +haplotype_evdience_alloc* hap, UC_Read* overlap_read) { window_list* beg_cigar; window_list* end_cigar; ///window_start, window_end, and useful_length correspond to x, instead of y long long useful_length = window_end - window_start + 1; - long long x_start, x_end, x_length; + long long x_start, x_length; char* x_string; char* y_string; long long i; - long long y_start, y_end, y_length; + long long y_start, y_length; long long overlapID, windowID; long long startNodeID, endNodeID, currentNodeID; @@ -9082,7 +9488,6 @@ haplotype_evdience_alloc* hap) int snp_threshold; snp_threshold = 1; - ///resize_window_list_alloc(boundary_cigars, dumy->length * 2); ///all overlaps related to the current window [window_start, window_end] ///first mark all snp pos @@ -9107,26 +9512,22 @@ haplotype_evdience_alloc* hap) x_start = overlap_list->list[overlapID].w_list[windowID].x_start; x_length = overlap_list->list[overlapID].w_list[windowID].x_end - overlap_list->list[overlapID].w_list[windowID].x_start + 1; - x_end = overlap_list->list[overlapID].w_list[windowID].x_end; - y_start = overlap_list->list[overlapID].w_list[windowID].y_start; y_length = overlap_list->list[overlapID].w_list[windowID].y_end - overlap_list->list[overlapID].w_list[windowID].y_start + 1; - y_end = overlap_list->list[overlapID].w_list[windowID].y_end; + beg_cigar = end_cigar = NULL; + + if(windowID >= 1) + { + beg_cigar = &(overlap_list->list[overlapID].boundary_cigars.buffer[windowID-1]); + } - beg_cigar = end_cigar = NULL; - - if(windowID >= 1) - { - beg_cigar = &(overlap_list->list[overlapID].boundary_cigars.buffer[windowID-1]); - } - - if(windowID < overlap_list->list[overlapID].w_list_length - 1) - { - end_cigar = &(overlap_list->list[overlapID].boundary_cigars.buffer[windowID]); - } + if(windowID < overlap_list->list[overlapID].w_list_length - 1) + { + end_cigar = &(overlap_list->list[overlapID].boundary_cigars.buffer[windowID]); + } markSNP_advance(window_start, x_start, x_length, y_start, y_length, @@ -9191,11 +9592,17 @@ haplotype_evdience_alloc* hap) - overlap_list->list[overlapID].w_list[windowID].y_start + 1; - recover_UC_Read_sub_region(dumy->overlap_region, y_start, y_length, overlap_list->list[overlapID].y_pos_strand, - R_INF, overlap_list->list[overlapID].y_id); - - x_string = r_string + x_start; - y_string = dumy->overlap_region; + if(overlap_list->list[overlapID].y_pos_strand == 0) + { + recover_UC_Read(overlap_read, R_INF, overlap_list->list[overlapID].y_id); + } + else + { + recover_UC_Read_RC(overlap_read, R_INF, overlap_list->list[overlapID].y_id); + } + + x_string = r_string; + y_string = overlap_read->seq; beg_cigar = end_cigar = NULL; @@ -9209,9 +9616,9 @@ haplotype_evdience_alloc* hap) } - addSNPtohaplotype(window_start, overlapID, x_string, x_start, x_length, - y_string, y_start, y_length, &(overlap_list->list[overlapID].w_list[windowID].cigar), - hap, snp_threshold); + addSNPtohaplotype_advance(window_start, overlapID, x_start, x_length, + y_start, y_length, &(overlap_list->list[overlapID].w_list[windowID]), + beg_cigar, end_cigar, hap, snp_threshold, x_string, y_string); } RsetInitHaplotypeEvdienceFlag(hap, first_snp, last_snp + 1 - first_snp); @@ -11706,7 +12113,7 @@ int force_repeat) return 0; } - // debug_print_snp_stat("m64016_190918_162737/174131552/ccs", hap, overlap_list, R_INF); + ///debug_print_snp_stat("m64013_190324_024932/23660629/ccs", hap, overlap_list, R_INF); ///if hap->available_snp == 1, the following codes would have bugs @@ -11929,7 +12336,7 @@ int force_repeat) } } - // debug_print_snp_stat("m64016_190918_162737/174131552/ccs", hap, overlap_list, R_INF); + ///debug_print_snp_stat("m64013_190324_024932/23660629/ccs", hap, overlap_list, R_INF); //if(hap->dp.max_snp_num > 0) @@ -12553,8 +12960,8 @@ void partition_overlaps(overlap_region_alloc* overlap_list, All_reads* R_INF, void partition_overlaps_advance(overlap_region_alloc* overlap_list, All_reads* R_INF, - UC_Read* g_read, Correct_dumy* dumy, haplotype_evdience_alloc* hap, - int force_repeat) + UC_Read* g_read, UC_Read* overlap_read, Correct_dumy* dumy, + haplotype_evdience_alloc* hap, int force_repeat) { ResizeInitHaplotypeEvdience(hap); @@ -12587,10 +12994,11 @@ void partition_overlaps_advance(overlap_region_alloc* overlap_list, All_reads* R } ///这个是available overlap里所有window的数量... - ///num_availiable_win = num_availiable_win + dumy->length + dumy->lengthNT; num_availiable_win = num_availiable_win + dumy->length; - - cluster_advance(g_read->seq, window_start, window_end, overlap_list, dumy, R_INF, hap); + + ///need to deal with + cluster_advance(g_read->seq, window_start, window_end, overlap_list, + dumy, R_INF, hap, overlap_read); } @@ -12617,6 +13025,13 @@ void partition_overlaps_advance(overlap_region_alloc* overlap_list, All_reads* R ////split reads for (i = 0; i < hap->length; i++) { + // if(overlap_list->list[0].x_id == 109837) + // { + // fprintf(stderr, "hap->list[%d].site: %d, type: %d, misBase: %c, overlapID: %d, %.*s\n", + // i, hap->list[i].site, hap->list[i].type, hap->list[i].misBase, hap->list[i].overlapID, + // Get_NAME_LENGTH((*R_INF), overlap_list->list[hap->list[i].overlapID].y_id), + // Get_NAME((*R_INF), overlap_list->list[hap->list[i].overlapID].y_id)); + // } if(pre_site != hap->list[i].site) { if(i != 0) @@ -12748,7 +13163,7 @@ void correct_overlap_back(overlap_region_alloc* overlap_list, All_reads* R_INF, void print_overlap(char* name, long long readID, -overlap_region_alloc* overlap_list, All_reads* R_INF, int output_reads) +overlap_region_alloc* overlap_list, All_reads* R_INF, int output_reads, int output_cigar) { if(memcmp(name, Get_NAME((*R_INF), readID), Get_NAME_LENGTH((*R_INF),readID)) == 0) @@ -12804,8 +13219,8 @@ overlap_region_alloc* overlap_list, All_reads* R_INF, int output_reads) if(overlap_list->list[i].is_match == 4) { fprintf(stderr, "%.*s\n", Get_NAME_LENGTH((*R_INF),overlap_list->list[i].y_id), - Get_NAME((*R_INF),overlap_list->list[i].y_id)); - fprintf(stderr, "alignLen: %d, x_s: %d, x_e: %d, y_s: %d, y_e: %d, y_dir: %d, strong: %d\n", + Get_NAME((*R_INF),overlap_list->list[i].y_id)); + fprintf(stderr, "alignLen: %d, x_s: %d, x_e: %d, y_s: %d, y_e: %d, y_dir: %d, strong: %d\n", overlap_list->list[i].align_length, overlap_list->list[i].x_pos_s, overlap_list->list[i].x_pos_e, @@ -12841,9 +13256,11 @@ overlap_region_alloc* overlap_list, All_reads* R_INF, int output_reads) fprintf(stderr, "Has already output all related reads\n\n"); destory_UC_Read(&g_read); + } - + if(output_cigar) + { for (i = 0; i < overlap_list->length; i++) { @@ -12859,7 +13276,7 @@ overlap_region_alloc* overlap_list, All_reads* R_INF, int output_reads) overlap_list->list[i].strong, overlap_list->list[i].is_match); - + fprintf(stderr, "window cigar: \n"); for (j = 0; j < overlap_list->list[i].w_list_length; j++) { fprintf(stderr, "************************\ncigar_j: %d, x_s: %d, x_e: %d, y_s: %d, y_end: %d\n", @@ -12882,20 +13299,50 @@ overlap_region_alloc* overlap_list, All_reads* R_INF, int output_reads) operationLen = cigar->C_L[cigar_i]; fprintf(stderr, "oper: %d, Len: %d\n", operation, operationLen); } - - } - } - - + + + + fprintf(stderr, "boundary cigar: \n"); + for (j = 0; j < overlap_list->list[i].boundary_cigars.length; j++) + { + fprintf(stderr, "###################\ncigar_j: %d, x_s: %d, x_e: %d, y_s: %d, y_end: %d\n", + j, overlap_list->list[i].boundary_cigars.buffer[j].x_start, + overlap_list->list[i].boundary_cigars.buffer[j].x_end, + overlap_list->list[i].boundary_cigars.buffer[j].y_start, + overlap_list->list[i].boundary_cigars.buffer[j].y_end); + + + if(overlap_list->list[i].boundary_cigars.buffer[j].y_end == -1) + { + fprintf(stderr, "not match\n"); + } + else + { + int cigar_i, operation, operationLen; + CIGAR* cigar = &(overlap_list->list[i].boundary_cigars.buffer[j].cigar); + fprintf(stderr, "length: %d\n", cigar->length); + for (cigar_i = 0; cigar_i < cigar->length; cigar_i++) + { + operation = cigar->C_C[cigar_i]; + operationLen = cigar->C_L[cigar_i]; + fprintf(stderr, "oper: %d, Len: %d\n", operation, operationLen); + } + } + } } } + + + fprintf(stderr, "End\n\n\n\n\n"); + } } + void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, Graph* g, Graph* DAGCon, long long* matched_overlap_0, long long* matched_overlap_1, @@ -12944,16 +13391,50 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, verify_window(window_start, window_end, overlap_list, dumy, R_INF, g_read->seq); } + /** recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read); - ///recalcate_window_advance(overlap_list, R_INF, g_read, dumy, overlap_read); - - partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat); - ///partition_overlaps_advance(overlap_list, R_INF, g_read, dumy, hap, force_repeat); + **/ + + recalcate_window_advance(overlap_list, R_INF, g_read, dumy, overlap_read); + ///partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat); + partition_overlaps_advance(overlap_list, R_INF, g_read, overlap_read, dumy, hap, force_repeat); + + + + + + + // print_overlap("m64013_190410_223304/61081096/ccs", + // overlap_list->list[0].x_id, overlap_list, R_INF, 0, 1); + + // print_overlap("m64013_190412_043951/83756835/ccs", + // overlap_list->list[0].x_id, overlap_list, R_INF, 0, 1); + + // print_overlap("m64013_190322_203854/152832751/ccs", + // overlap_list->list[0].x_id, overlap_list, R_INF, 0, 1); + + // print_overlap("m64011_190329_072846/175047040/ccs", + // overlap_list->list[0].x_id, overlap_list, R_INF, 0, 1); + + // print_overlap("m64013_190322_203854/13174300/ccs", + // overlap_list->list[0].x_id, overlap_list, R_INF, 0, 1); + + // print_overlap("m64013_190322_203854/177145456/ccs", + // overlap_list->list[0].x_id, overlap_list, R_INF, 0, 1); + + // print_overlap("m64013_190322_203854/120717695/ccs", + // overlap_list->list[0].x_id, overlap_list, R_INF, 0, 1); + + // print_overlap("m64011_190329_072846/76548728/ccs", + // overlap_list->list[0].x_id, overlap_list, R_INF, 0, 1); + + // print_overlap("m64013_190324_024932/23660629/ccs", + // overlap_list->list[0].x_id, overlap_list, R_INF, 0, 1); - // print_overlap("m64016_190918_162737/174131552/ccs", - // overlap_list->list[0].x_id, overlap_list, R_INF, 1); + + if(is_consensus) { generate_consensus(overlap_list, R_INF, g_read, dumy, g, DAGCon, current_cigar, second_round); diff --git a/Correct.h b/Correct.h index 95b068c..8aff96e 100644 --- a/Correct.h +++ b/Correct.h @@ -1447,6 +1447,10 @@ void append_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list, k_mer_pos uint64_t n_end_pos, uint8_t n_direction, UC_Read* g_read, All_reads* R_INF, Correct_dumy* dumy); /**********************for prefilter************************ */ +int verify_single_window(long long x_start, long long x_end, +long long overlap_x_s, long long overlap_y_s, int x_id, +int y_id, int y_strand, char* x_buffer, char* y_buffer, +All_reads* R_INF); void init_Cigar_record_alloc(Cigar_record_alloc* x); void resize_Cigar_record_alloc(Cigar_record_alloc* x, long long new_size); diff --git a/Hash_Table.h b/Hash_Table.h index 175084b..dff6500 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -28,7 +28,7 @@ typedef khash_t(POS64) Pos_Table; #define OVERLAP_THRESHOLD_FILTER 0.9 #define WINDOW_MAX_SIZE WINDOW + TAIL_LENGTH + 3 #define THRESHOLD_MAX_SIZE 31 - +#define FINAL_OVERLAP_ERROR_RATE 0.03 #define GROUP_SIZE 4 ///最长是10M10D10M10D10M这种 diff --git a/Overlaps.cpp b/Overlaps.cpp index 7458830..035d376 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -6,6 +6,8 @@ #include "ksort.h" #include "Process_Read.h" #include "CommandLines.h" +#include "Hash_Table.h" +#include "Correct.h" KDQ_INIT(uint64_t) @@ -55,6 +57,109 @@ void asg_arc_sort(asg_t *g) radix_sort_asg(g->arc, g->arc + g->n_arc); } + +void add_overlaps(ma_hit_t_alloc* source_paf, ma_hit_t_alloc* dest_paf, uint64_t* source_index, long long listLen) +{ + long long i; + ma_hit_t ele; + ma_hit_t* tmp; + for (i = 0; i < listLen; i++) + { + tmp = &(source_paf->buffer[(uint32_t)(source_index[i])]); + /** + ele.rev = tmp->rev; + ele.qns = Get_tn((*tmp)); + ele.qns = ele.qns << 32; + ele.qns = ele.qns | (uint64_t)(Get_ts((*tmp))); + ele.qe = Get_te((*tmp)); + + ele.tn = Get_qn((*tmp)); + ele.ts = Get_qs((*tmp)); + ele.te = Get_qe((*tmp)); + + ele.bl = R_INF.read_length[ele.tn]; + ele.ml = tmp->ml; + ele.el = tmp->el; + ele.no_l_indel = tmp->no_l_indel; + + add_ma_hit_t_alloc(dest_paf, &ele); + **/ + add_ma_hit_t_alloc(dest_paf, tmp); + } +} + + +void remove_overlaps(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen) +{ + long long i, m; + ma_hit_t ele; + ma_hit_t* tmp; + for (i = 0; i < listLen; i++) + { + source_paf->buffer[(uint32_t)(source_index[i])].qns = (uint64_t)(-1); + } + + m = 0; + for (i = 0; i < source_paf->length; i++) + { + if(source_paf->buffer[i].qns != (uint64_t)(-1)) + { + source_paf->buffer[m] = source_paf->buffer[i]; + m++; + } + } + source_paf->length = m; +} + + +void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf, +uint64_t* source_index, long long listLen) +{ + long long i; + ma_hit_t ele; + ma_hit_t* tmp; + uint32_t source_n, source_i; + for (i = 0; i < listLen; i++) + { + source_n = source_index[i] >> 32; + source_i = (uint32_t)(source_index[i]); + tmp = &(source_paf_list[source_n].buffer[source_i]); + + ele.rev = tmp->rev; + ele.qns = Get_tn((*tmp)); + ele.qns = ele.qns << 32; + ele.qns = ele.qns | (uint64_t)(Get_ts((*tmp))); + ele.qe = Get_te((*tmp)); + + ele.tn = Get_qn((*tmp)); + ele.ts = Get_qs((*tmp)); + ele.te = Get_qe((*tmp)); + + ele.bl = R_INF.read_length[ele.tn]; + ele.ml = tmp->ml; + ele.el = tmp->el; + ele.no_l_indel = tmp->no_l_indel; + + add_ma_hit_t_alloc(dest_paf, &ele); + } +} + +void print_revise_edges(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen) +{ + long long i, m; + ma_hit_t ele; + ma_hit_t* tmp; + for (i = 0; i < listLen; i++) + { + tmp = &(source_paf->buffer[(uint32_t)(source_index[i])]); + + fprintf(stderr, "%.*s(%d) ---(+)--> %.*s(%d), Len: %d\n", + Get_NAME_LENGTH(R_INF, Get_qn((*tmp))), Get_NAME(R_INF, Get_qn((*tmp))), Get_qn((*tmp)), + Get_NAME_LENGTH(R_INF, Get_tn((*tmp))), Get_NAME(R_INF, Get_tn((*tmp))), Get_tn((*tmp)), + Get_qe((*tmp)) - Get_qs((*tmp))); + } +} + void ma_ug_destroy(ma_ug_t *ug) { uint32_t i; @@ -347,6 +452,7 @@ inline void set_reverse_overlap(ma_hit_t* dest, ma_hit_t* source) /****************************may have bugs********************************/ + /** if(dest->ml == 0 || source->ml == 0) { dest->ml = source->ml = 0; @@ -365,7 +471,9 @@ inline void set_reverse_overlap(ma_hit_t* dest, ma_hit_t* source) { dest->no_l_indel = source->no_l_indel = 1; } - + **/ + dest->ml = source->ml; + dest->no_l_indel = source->no_l_indel; /****************************may have bugs********************************/ dest->bl = Get_qe(*dest) - Get_qs(*dest); } @@ -970,6 +1078,16 @@ long long n_read, uint64_t* readLen, ma_sub_t* coverage_cut) } else { + // if(coverage_cut[i].c == 2) + // { + // fprintf(stderr, "\ni: %d, readLen: %d\n", i, readLen[i]); + // for (j = 0; j < sources[i].length; j++) + // { + // fprintf(stderr, "j: %d, qs: %d, qe: %d, ts: %d, te: %d\n", + // j, Get_qs(sources[i].buffer[j]), Get_qe(sources[i].buffer[j]), + // Get_ts(sources[i].buffer[j]), Get_te(sources[i].buffer[j])); + // } + // } coverage_cut[i].c = 0; } } @@ -981,6 +1099,487 @@ long long n_read, uint64_t* readLen, ma_sub_t* coverage_cut) +int boundary_verify_back(uint32_t interval_s, uint32_t interval_e, ma_hit_t* map, +char* x_buffer, char* y_buffer, All_reads* R_INF) +{ + uint32_t xs, xe, ys, ye, dir, x_id, y_id, xLen; + dir = (*map).rev; + xs = Get_qs((*map)); + xe = Get_qe((*map)) - 1; + x_id = Get_qn((*map)); + y_id = Get_tn((*map)); + long long yLen = Get_READ_LENGTH((*R_INF), y_id); + + if(dir == 1) + { + ys = yLen - (Get_te((*map)) - 1) - 1; + ye = yLen - Get_ts((*map)) - 1; + + } + else + { + ys = Get_ts((*map)); + ye = Get_te((*map)) - 1; + } + + xLen = interval_e - interval_s; + + if(xLen <= WINDOW) + { + return verify_single_window(interval_s, interval_e-1, xs, ys, x_id, y_id, + dir, x_buffer, y_buffer, R_INF); + } + else + { + if(verify_single_window(interval_s, interval_s + WINDOW - 1, + xs, ys, x_id, y_id, dir, x_buffer, y_buffer, R_INF) == 0) + { + return 0; + } + + if(verify_single_window(interval_e - WINDOW, interval_e-1, + xs, ys, x_id, y_id, dir, x_buffer, y_buffer, R_INF) == 0) + { + return 0; + } + + return 1; + } +} + +int boundary_verify(uint32_t x_interval_s, uint32_t x_interval_e, ma_hit_t* map, +char* x_buffer, char* y_buffer, All_reads* R_INF) +{ + uint32_t xs, xe, ys, ye, dir, x_id, y_id, x_interval_Len, y_interval_Len, y_interval_s, y_interval_e; + dir = (*map).rev; + xs = Get_qs((*map)); + xe = Get_qe((*map)) - 1; + x_id = Get_qn((*map)); + y_id = Get_tn((*map)); + long long yLen = Get_READ_LENGTH((*R_INF), y_id); + + if(dir == 1) + { + ys = yLen - (Get_te((*map)) - 1) - 1; + ye = yLen - Get_ts((*map)) - 1; + + } + else + { + ys = Get_ts((*map)); + ye = Get_te((*map)) - 1; + } + ///[x_interval_s, x_interval_e) + x_interval_Len = x_interval_e - x_interval_s; + + ///[y_interval_s, y_interval_e] + y_interval_s = (x_interval_s - xs) + ys; + if(y_interval_s >= yLen) + { + return 0; + } + y_interval_e = y_interval_s + x_interval_Len - 1; + if(y_interval_e >= yLen) + { + y_interval_e = yLen - 1; + } + + if(y_interval_e < y_interval_s) + { + return 0; + } + + y_interval_Len = y_interval_e - y_interval_s + 1; + + + if(y_interval_Len <= WINDOW) + { + return verify_single_window(y_interval_s, y_interval_e, ys, xs, y_id, x_id, + dir, y_buffer, x_buffer, R_INF); + } + else + { + + if(verify_single_window(y_interval_s, y_interval_s + WINDOW - 1, ys, xs, y_id, x_id, + dir, y_buffer, x_buffer, R_INF) == 0) + { + return 0; + } + + if(verify_single_window(y_interval_e - WINDOW + 1, y_interval_e, ys, xs, y_id, x_id, + dir, y_buffer, x_buffer, R_INF) == 0) + { + return 0; + } + + return 1; + } +} + +void collect_sides(ma_hit_t_alloc* paf, uint64_t rLen, ma_sub_t* max_left, ma_sub_t* max_right) +{ + long long j; + uint32_t qs, qe; + for (j = 0; j < paf->length; j++) + { + qs = Get_qs(paf->buffer[j]); + qe = Get_qe(paf->buffer[j]); + + + ///overlaps from left side + if(qs == 0) + { + if(qs < max_left->s) max_left->s = qs; + if(qe > max_left->e) max_left->e = qe; + } + + ///overlaps from right side + if(qe == rLen) + { + if(qs < max_right->s) max_right->s = qs; + if(qe > max_right->e) max_right->e = qe; + } + + ///note: if (qs == 0 && qe == rLen) + ///this overlap would be added to both b_left and b_right + ///that is what we want + } +} + +void collect_contain(ma_hit_t_alloc* paf, uint64_t rLen, ma_sub_t* max_left, ma_sub_t* max_right) +{ + long long j; + uint32_t qs, qe; + for (j = 0; j < paf->length; j++) + { + qs = Get_qs(paf->buffer[j]); + qe = Get_qe(paf->buffer[j]); + ///check contained overlaps + if(qs != 0 && qe != rLen) + { + ///[qs, qe), [max_left.s, max_left.e) + if(qs < max_left->e && qe > max_left->e) + { + if(qe > max_left->e) max_left->e = qe; + } + + ///[qs, qe), [max_right.s, max_right.e) + if(qs < max_right->s && qe > max_right->s) + { + if(qs < max_right->s) max_right->s = qs; + } + } + } +} + + + +int intersection_check(ma_hit_t_alloc* paf, uint64_t rLen, uint32_t interval_s, uint32_t interval_e) +{ + long long j, cov = 0; + uint32_t qs, qe; + + for (j = 0; j < paf->length; j++) + { + qs = Get_qs(paf->buffer[j]); + qe = Get_qe(paf->buffer[j]); + ///[interval_s, interval_e) must be at least contained at one of the [qs, qe) + if(qs<=interval_s && qe>=interval_e) + { + cov++; + } + } + + return cov; +} + + +int intersection_check_by_base(ma_hit_t_alloc* paf, uint64_t rLen, uint32_t interval_s, uint32_t interval_e, +char* bq, char* bt) +{ + long long j; + uint32_t qs, qe; + + for (j = 0; j < paf->length; j++) + { + qs = Get_qs(paf->buffer[j]); + qe = Get_qe(paf->buffer[j]); + ///[interval_s, interval_e) must be at least contained at one of the [qs, qe) + if(qs<=interval_s && qe>=interval_e) + { + if(boundary_verify(interval_s, interval_e, &(paf->buffer[j]), bq, bt, &R_INF) == 0) + { + return 1; + } + } + } + + return 0; +} + + + +void detect_chimeric_reads(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, +long long n_read, uint64_t* readLen, ma_sub_t* coverage_cut, float shift_rate) +{ + double startTime = Get_T(); + init_aux_table(); + long long i, j, rLen, /**cov,**/ n_simple_remove = 0, n_complex_remove = 0, n_complex_remove_real = 0; + uint32_t qs, qe; + uint32_t interval_s, interval_e; + ma_sub_t max_left, max_right; + kvec_t(char) b_q = {0,0,0}; + kvec_t(char) b_t = {0,0,0}; + for (i = 0; i < n_read; ++i) + { + coverage_cut[i].c = 0; + rLen = readLen[i]; + + + max_left.s = max_right.s = rLen; + max_left.e = max_right.e = 0; + + + collect_sides(&(paf[i]), rLen, &max_left, &max_right); + collect_sides(&(rev_paf[i]), rLen, &max_left, &max_right); + ///that means this read is an end node + if(max_left.s == rLen || max_right.s == rLen) + { + continue; + } + + + collect_contain(&(paf[i]), rLen, &max_left, &max_right); + collect_contain(&(rev_paf[i]), rLen, &max_left, &max_right); + ////shift_rate should be (FINAL_OVERLAP_ERROR_RATE*2) + ///this read is a normal read + if(max_left.e > max_right.s && + (max_left.e - max_right.s >= rLen * shift_rate)) + { + ///coverage_cut[i].c = 0; + continue; + } + + ///simple chimeric reads + if(max_left.e <= max_right.s) + { + ///coverage_cut[i].c = 2; + coverage_cut[i].c = 1; + coverage_cut[i].del = 1; + paf[i].length = 0; + n_simple_remove++; + continue; + } + + ///now max_left.e > max_right.s && max_left.e - max_right.s is small enough + //[interval_s, interval_e) + interval_s = max_right.s; + interval_e = max_left.e; + + /** + cov = 0; + cov += intersection_check(&(paf[i]), rLen, interval_s, interval_e); + cov += intersection_check(&(rev_paf[i]), rLen, interval_s, interval_e); + if(interval_e - interval_s < WINDOW && cov <= 2) + { + coverage_cut[i].c = 1; + coverage_cut[i].del = 1; + paf[i].length = 0; + n_complex_remove_real++; + } + else**/ + { + kv_resize(char, b_q, WINDOW*4+20); + kv_resize(char, b_t, WINDOW*4+20); + if(intersection_check_by_base(&(paf[i]), rLen, interval_s, interval_e, b_q.a, b_t.a) + || + intersection_check_by_base(&(rev_paf[i]), rLen, interval_s, interval_e, b_q.a, b_t.a)) + { + coverage_cut[i].c = 1; + coverage_cut[i].del = 1; + paf[i].length = 0; + n_complex_remove_real++; + } + } + + n_complex_remove++; + } + + free(b_q.a); + free(b_t.a); + + fprintf(stderr, "[M::%s] takes %0.2f s, n_simple_remove: %d, n_complex_remove: %d/%d\n\n", + __func__, Get_T()-startTime, n_simple_remove, n_complex_remove_real, n_complex_remove); +} + + +void detect_chimeric_reads_back(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, +long long n_read, uint64_t* readLen, ma_sub_t* coverage_cut, float shift_rate) +{ + double startTime = Get_T(); + long long i, j, cov, rLen, n_simple_remove = 0, n_complex_remove = 0, n_complex_remove_real = 0; + uint32_t qs, qe; + uint32_t interval_s, interval_e; + ma_sub_t max_left, max_right; + kvec_t(char) b_q = {0,0,0}; + kvec_t(char) b_t = {0,0,0}; + for (i = 0; i < n_read; ++i) + { + coverage_cut[i].c = 0; + rLen = readLen[i]; + + + max_left.s = max_right.s = rLen; + max_left.e = max_right.e = 0; + + for (j = 0; j < paf[i].length; j++) + { + qs = Get_qs(paf[i].buffer[j]); + qe = Get_qe(paf[i].buffer[j]); + + + ///overlaps from left side + if(qs == 0) + { + if(qs < max_left.s) max_left.s = qs; + if(qe > max_left.e) max_left.e = qe; + } + + ///overlaps from right side + if(qe == rLen) + { + if(qs < max_right.s) max_right.s = qs; + if(qe > max_right.e) max_right.e = qe; + } + + ///note: if (qs == 0 && qe == rLen) + ///this overlap would be added to both b_left and b_right + ///that is what we want + } + + + ///that means this read is an end node + ///if(max_left.e == 0 || max_right.e == 0) + if(max_left.s == rLen || max_right.s == rLen) + { + continue; + } + + + for (j = 0; j < paf[i].length; j++) + { + qs = Get_qs(paf[i].buffer[j]); + qe = Get_qe(paf[i].buffer[j]); + ///check contained overlaps + if(qs != 0 && qe != rLen) + { + ///[qs, qe), [max_left.s, max_left.e) + if(qs < max_left.e && qe > max_left.e) + { + if(qe > max_left.e) max_left.e = qe; + } + + ///[qs, qe), [max_right.s, max_right.e) + if(qs < max_right.s && qe > max_right.s) + { + if(qs < max_right.s) max_right.s = qs; + } + } + } + + ////shift_rate should be (FINAL_OVERLAP_ERROR_RATE*2) + ///this read is a normal read + if(max_left.e > max_right.s && + (max_left.e - max_right.s >= rLen * shift_rate)) + { + ///coverage_cut[i].c = 0; + continue; + } + + ///simple chimeric reads + if(max_left.e <= max_right.s) + { + ///coverage_cut[i].c = 2; + coverage_cut[i].c = 1; + coverage_cut[i].del = 1; + paf[i].length = 0; + n_simple_remove++; + continue; + } + + ///now max_left.e > max_right.s && max_left.e - max_right.s is small enough + //[interval_s, interval_e) + interval_s = max_right.s; + interval_e = max_left.e; + + kv_resize(char, b_q, rLen+10); + kv_resize(char, b_t, (rLen*2+10)); + + cov = 0; + for (j = 0; j < paf[i].length; j++) + { + qs = Get_qs(paf[i].buffer[j]); + qe = Get_qe(paf[i].buffer[j]); + ///[interval_s, interval_e) must be at least contained at one of the [qs, qe) + if(qs<=interval_s && qe>=interval_e) + { + // if(i == 439960) + // { + // fprintf(stderr, "?????? interval_s: %d, interval_e: %d, qn: %d, tn: %d, j: %d, qs: %d, qe: %d, ts: %d, te: %d, dir: %d\n", + // interval_s, interval_e, + // Get_qn(paf[i].buffer[j]), Get_tn(paf[i].buffer[j]), + // j, Get_qs(paf[i].buffer[j]), Get_qe(paf[i].buffer[j]), + // Get_ts(paf[i].buffer[j]), Get_te(paf[i].buffer[j]), + // paf[i].buffer[j].rev); + + // fprintf(stderr, "%.*s\n", Get_NAME_LENGTH(R_INF, Get_tn(paf[i].buffer[j])), + // Get_NAME(R_INF, Get_tn(paf[i].buffer[j]))); + // } + ///num++; + cov++; + if(boundary_verify(interval_s, interval_e, &(paf[i].buffer[j]), + b_q.a, b_t.a, &R_INF) == 0) + { + coverage_cut[i].c = 1; + coverage_cut[i].del = 1; + paf[i].length = 0; + n_complex_remove_real++; + break; + } + } + } + + if(coverage_cut[i].c == 0 && interval_e - interval_s < WINDOW && cov <= 2) + { + coverage_cut[i].c = 1; + coverage_cut[i].del = 1; + paf[i].length = 0; + n_complex_remove_real++; + } + + + + + + ///if all of them are wrong + // if(j == paf[i].length) + // { + // coverage_cut[i].c = 1; + // coverage_cut[i].del = 1; + // paf[i].length = 0; + // n_complex_remove_real++; + // } + + + n_complex_remove++; + } + + free(b_q.a); + free(b_t.a); + + fprintf(stderr, "[M::%s] takes %0.2f s, n_simple_remove: %d, n_complex_remove: %d/%d\n\n", + __func__, Get_T()-startTime, n_simple_remove, n_complex_remove_real, n_complex_remove); +} @@ -8199,8 +8798,8 @@ ma_hit_t_alloc* reverse_sources, int id, char* command) { fprintf(stderr, "\n\n\nafter %s\n", command); - fprintf(stderr, "****************ma_hit_t ref_read: %.*s****************\n", - Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + fprintf(stderr, "****************ma_hit_t (%d)ref_read: %.*s****************\n", + i, Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); fprintf(stderr, "sources Len: %d, is_fully_corrected: %d\n", @@ -9558,6 +10157,326 @@ int output_tips(asg_t *g, const All_reads *RNF) } } +void collect_abnormal_edges(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long readNum) +{ + double startTime = Get_T(); + long long T_edges, T_Single_Dir_Edges_0, T_Single_Dir_Edges_1, T_Conflict_Equal_Edges, T_Conflict_Strong_Edges; + T_edges = T_Single_Dir_Edges_0 = T_Single_Dir_Edges_1 = T_Conflict_Equal_Edges = T_Conflict_Strong_Edges = 0; + long long T_Single_Dir_Edges_1_1000 = 0; + long long related_reads = 0; + long long related_overlaps = 0; + long long i, j; + uint32_t qn, tn; + int is_equal_f, is_strong_f, is_exist_f; + int is_equal_b, is_strong_b, is_exist_b; + + kvec_t(uint64_t) edge_vector; + kv_init(edge_vector); + + for (i = 0; i < readNum; i++) + { + for (j = 0; j < paf[i].length; j++) + { + qn = Get_qn(paf[i].buffer[j]); + tn = Get_tn(paf[i].buffer[j]); + T_edges++; + + is_equal_f = paf[i].buffer[j].el; + is_strong_f = paf[i].buffer[j].ml; + is_exist_f = 1; + + is_exist_b = get_specific_overlap(&(paf[tn]), tn, qn); + if(is_exist_b == -1) + { + is_exist_b = get_specific_overlap(&(rev_paf[tn]), tn, qn); + if(is_exist_b != -1) + { + T_Single_Dir_Edges_0++; + kv_push(uint64_t, edge_vector, qn); + kv_push(uint64_t, edge_vector, tn); + ///related_overlaps += paf[qn].length + rev_paf[qn].length + paf[tn].length + rev_paf[tn].length; + // fprintf(stderr, "%.*s(%d) ---(+)--> %.*s(%d), Len: %d\n", + // Get_NAME_LENGTH(R_INF, qn), Get_NAME(R_INF, qn), qn, + // Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn), tn, + // Get_qe(paf[i].buffer[j]) - Get_qs(paf[i].buffer[j])); + + // fprintf(stderr, "%.*s(%d) ---(-)--> %.*s(%d), Len: %d\n\n", + // Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn), tn, + // Get_NAME_LENGTH(R_INF, qn), Get_NAME(R_INF, qn), qn, + // Get_qe(rev_paf[tn].buffer[is_exist_b]) - Get_qs(rev_paf[tn].buffer[is_exist_b])); + } + else + { + T_Single_Dir_Edges_1++; + if(Get_qe(paf[i].buffer[j]) - Get_qs(paf[i].buffer[j]) >= 1000) + { + T_Single_Dir_Edges_1_1000++; + // fprintf(stderr, "%.*s(%d) ---(%d)--> %.*s(%d), Len: %d\n\n", + // Get_NAME_LENGTH(R_INF, qn), Get_NAME(R_INF, qn), qn, + // is_strong_f, + // Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn), tn, + // Get_qe(paf[i].buffer[j]) - Get_qs(paf[i].buffer[j])); + } + } + + related_reads = related_reads + 2; + } + else + { + is_equal_b = paf[tn].buffer[is_exist_b].el; + is_strong_b = paf[tn].buffer[is_exist_b].ml; + + if(is_equal_f != is_equal_b) + { + T_Conflict_Equal_Edges++; + + // fprintf(stderr, "%.*s(%d) ---(%d)--> %.*s(%d), Len: %d\n", + // Get_NAME_LENGTH(R_INF, qn), Get_NAME(R_INF, qn), qn, + // is_equal_f, + // Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn), tn, + // Get_qe(paf[i].buffer[j]) - Get_qs(paf[i].buffer[j])); + + // fprintf(stderr, "%.*s(%d) ---(%d)--> %.*s(%d), Len: %d\n\n", + // Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn), tn, + // is_equal_b, + // Get_NAME_LENGTH(R_INF, qn), Get_NAME(R_INF, qn), qn, + // Get_qe(paf[tn].buffer[is_exist_b]) - Get_qs(paf[tn].buffer[is_exist_b])); + } + + if(is_strong_f != is_strong_b) + { + T_Conflict_Strong_Edges++; + // fprintf(stderr, "%.*s(%d) ---(%d)--> %.*s(%d), Len: %d\n", + // Get_NAME_LENGTH(R_INF, qn), Get_NAME(R_INF, qn), qn, + // is_strong_f, + // Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn), tn, + // Get_qe(paf[i].buffer[j]) - Get_qs(paf[i].buffer[j])); + + // fprintf(stderr, "%.*s(%d) ---(%d)--> %.*s(%d), Len: %d\n\n", + // Get_NAME_LENGTH(R_INF, tn), Get_NAME(R_INF, tn), tn, + // is_strong_b, + // Get_NAME_LENGTH(R_INF, qn), Get_NAME(R_INF, qn), qn, + // Get_qe(paf[tn].buffer[is_exist_b]) - Get_qs(paf[tn].buffer[is_exist_b])); + } + + if(is_equal_f != is_equal_b || is_strong_f != is_strong_b) + { + related_reads++; + } + } + } + } + + radix_sort_arch64(edge_vector.a, edge_vector.a + edge_vector.n); + uint64_t pre = (uint64_t)-1; + long long mn = 0; + for (i = 0; i < edge_vector.n; i++) + { + if(pre != edge_vector.a[i]) + { + mn++; + pre = edge_vector.a[i]; + related_overlaps += paf[pre].length + rev_paf[pre].length; + } + + if(i>0 && edge_vector.a[i] < edge_vector.a[i-1]) fprintf(stderr, "hehe\n"); + } + + + + fprintf(stdout, "****************statistic for abnormal overlaps****************\n"); + fprintf(stdout, "overlaps #: %lld\n", T_edges); + fprintf(stdout, "one direction overlaps (different phasing)#: %lld\n", T_Single_Dir_Edges_0); + fprintf(stdout, "one direction overlaps (missing)#: %lld\n", T_Single_Dir_Edges_1); + fprintf(stdout, "one direction overlaps (missing) >= 1000#: %lld\n", T_Single_Dir_Edges_1_1000); + fprintf(stdout, "conflict strong/weak overlaps #: %lld\n", T_Conflict_Strong_Edges); + fprintf(stdout, "conflict exact/inexact overlaps #: %lld\n", T_Conflict_Equal_Edges); + fprintf(stdout, "related_reads #: %lld/%lld\n", related_reads, mn); + fprintf(stdout, "related_overlaps #: %lld\n", related_overlaps); + + fprintf(stdout, "****************statistic for abnormal overlaps****************\n"); + + fprintf(stderr, "[M::%s] took %0.2fs\n\n", __func__, Get_T()-startTime); + + kv_destroy(edge_vector); +} + +///if we don't have this function, we just simply remove all one-direction edges +///by utilizing this function, some one-direction edges can be recovered as two-direction edges +void try_rescue_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long readNum, +long long rescue_threshold) +{ + double startTime = Get_T(); + long long i, j, revises = 0; + uint32_t qn, tn, qs, qe; + + kvec_t(uint64_t) edge_vector; + kv_init(edge_vector); + kvec_t(uint64_t) edge_vector_index; + kv_init(edge_vector_index); + kvec_t(uint32_t) b; + kv_init(b); + uint64_t flag; + int index; + + for (i = 0; i < readNum; i++) + { + edge_vector.n = 0; + edge_vector_index.n = 0; + for (j = 0; j < rev_paf[i].length; j++) + { + qn = Get_qn(rev_paf[i].buffer[j]); + tn = Get_tn(rev_paf[i].buffer[j]); + index = get_specific_overlap(&(paf[tn]), tn, qn); + if(index != -1) + { + flag = tn; + flag = flag << 32; + flag = flag | (uint64_t)(index); + kv_push(uint64_t, edge_vector, flag); + kv_push(uint64_t, edge_vector_index, j); + } + } + + ///based on qn, all edges at edge_vector/edge_vector_index come from different haplotype + ///but at another direction, all these edges come from the same haplotype + //here we want to recover these edges + if(edge_vector_index.n >= rescue_threshold) + { + kv_resize(uint32_t, b, edge_vector_index.n); + b.n = 0; + for (j = 0; j < edge_vector_index.n; j++) + { + qs = Get_qs(rev_paf[i].buffer[edge_vector_index.a[j]]); + qe = Get_qe(rev_paf[i].buffer[edge_vector_index.a[j]]); + kv_push(uint32_t, b, qs<<1); + kv_push(uint32_t, b, qe<<1|1); + } + + ks_introsort_uint32_t(b.n, b.a); + int dp, start, max_dp = 0; + ma_sub_t max_interval; + for (j = 0, dp = 0; j < b.n; ++j) + { + int old_dp = dp; + ///if a[j] is qe + if (b.a[j]&1) + { + --dp; + } + else + { + ++dp; + } + + /** + there are two cases: + 1. old_dp = dp + 1 (b.a[j] is qe); 2. old_dp = dp - 1 (b.a[j] is qs); + **/ + ///if (old_dp < min_dp && dp >= min_dp) ///old_dp < dp, b.a[j] is qs + if(old_dp < dp) ///b.a[j] is qs + { + ///case 2, a[j] is qs + //here should use dp >= max_dp, instead of dp > max_dp + if(dp >= max_dp) + { + start = b.a[j]>>1; + max_dp = dp; + } + } + ///else if (old_dp >= min_dp && dp < min_dp) ///old_dp > min_dp, b.a[j] is qe + else if (old_dp > dp) ///old_dp > min_dp, b.a[j] is qe + { + if(old_dp == max_dp) + { + max_interval.s = start; + max_interval.e = b.a[j]>>1; + } + } + // else + // { + // fprintf(stderr, "error\n"); + // } + + } + + if(max_dp>= rescue_threshold) + { + long long m = 0; + for (j = 0; j < edge_vector_index.n; j++) + { + qs = Get_qs(rev_paf[i].buffer[edge_vector_index.a[j]]); + qe = Get_qe(rev_paf[i].buffer[edge_vector_index.a[j]]); + if(qs <= max_interval.s && qe >= max_interval.e) + { + edge_vector_index.a[m] = edge_vector_index.a[j]; + edge_vector.a[m] = edge_vector.a[j]; + m++; + } + } + edge_vector_index.n = m; + edge_vector.n = m; + + // if(max_dp != edge_vector_index.n) + // { + // fprintf(stderr, "error\n"); + // } + + ///the read itself do not have these overlaps, but all related reads have + ///we need to remove all overlaps from rev_paf[i], and then add all overlaps to paf[i] + // fprintf(stderr,"\nadd following %d edges...\n", edge_vector.n); + // print_revise_edges(&(rev_paf[i]), edge_vector_index.a, edge_vector_index.n); + remove_overlaps(&(rev_paf[i]), edge_vector_index.a, edge_vector_index.n); + add_overlaps_from_different_sources(paf, &(paf[i]), edge_vector.a, edge_vector.n); + revises = revises + edge_vector.n; + + ///for debug + /** + for (j = 0; j < edge_vector.n; j++) + { + int index_qn, index_tn; + qn = i; + tn = (uint32_t)(edge_vector.a[j] >> 32); + + index_qn = get_specific_overlap(&(paf[qn]), qn, tn); + if(index_qn == -1) fprintf(stderr, "error1\n"); + + index_tn = get_specific_overlap(&(paf[tn]), tn, qn); + if(index_tn == -1) fprintf(stderr, "error2\n"); + + if(paf[qn].buffer[index_qn].el != paf[tn].buffer[index_tn].el || + paf[qn].buffer[index_qn].ml != paf[tn].buffer[index_tn].ml || + paf[qn].buffer[index_qn].no_l_indel != paf[tn].buffer[index_tn].no_l_indel || + paf[qn].buffer[index_qn].rev != paf[tn].buffer[index_tn].rev || + Get_qs(paf[qn].buffer[index_qn]) != Get_ts(paf[tn].buffer[index_tn]) || + Get_qe(paf[qn].buffer[index_qn]) != Get_te(paf[tn].buffer[index_tn]) || + Get_qn(paf[qn].buffer[index_qn]) != Get_tn(paf[tn].buffer[index_tn]) || + Get_tn(paf[qn].buffer[index_qn]) != Get_qn(paf[tn].buffer[index_tn]) + ) + { + fprintf(stderr, "error full\n"); + } + + index_qn = get_specific_overlap(&(rev_paf[qn]), qn, tn); + if(index_qn != -1) fprintf(stderr, "error3\n"); + + + } + **/ + } + } + } + + + kv_destroy(edge_vector); + kv_destroy(edge_vector_index); + kv_destroy(b); + + fprintf(stderr, "[M::%s] took %0.2fs, rescue edges #: %lld\n\n", __func__, Get_T()-startTime, revises); +} + + + void build_string_graph_without_clean( int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long n_read, uint64_t* readLen, @@ -9572,11 +10491,19 @@ char* output_file_name, long long bubble_dist, int read_graph, int write) &R_INF, output_file_name); } + + ///collect_abnormal_edges(sources, reverse_sources, n_read); - // debug_info_of_specfic_read("m64016_190918_162737/179635219/ccs", + try_rescue_overlaps(sources, reverse_sources, n_read, 4); + + ///collect_abnormal_edges(sources, reverse_sources, n_read); + + + + // debug_info_of_specfic_read("m64016_190918_162737/76482648/ccs", // sources, reverse_sources, -1, "init"); - // debug_info_of_specfic_read("m64016_190918_162737/130811282/ccs", + // debug_info_of_specfic_read("m64016_190918_162737/133368200/ccs", // sources, reverse_sources, -1, "init"); @@ -9617,8 +10544,13 @@ char* output_file_name, long long bubble_dist, int read_graph, int write) ma_hit_sub(min_dp, sources, n_read, readLen, mini_overlap_length, &coverage_cut); + + + detect_chimeric_reads(sources, reverse_sources, n_read, readLen, coverage_cut, + FINAL_OVERLAP_ERROR_RATE*2); + ////这个会断开 - ma_hit_chimeric(1, sources, reverse_sources, n_read, readLen, coverage_cut); + ///ma_hit_chimeric(1, sources, reverse_sources, n_read, readLen, coverage_cut); ma_hit_cut(min_dp, sources, n_read, readLen, mini_overlap_length, &coverage_cut); ///it seems we do not need ma_hit_flt diff --git a/Overlaps.h b/Overlaps.h index 9538ee8..c28625a 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -330,4 +330,11 @@ int write); void debug_info_of_specfic_read(char* name, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int id, char* command); +void collect_abnormal_edges(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long readNum); + +void add_overlaps(ma_hit_t_alloc* source_paf, ma_hit_t_alloc* dest_paf, uint64_t* source_index, long long listLen); +void remove_overlaps(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen); +void add_overlaps_from_different_sources(ma_hit_t_alloc* source_paf_list, ma_hit_t_alloc* dest_paf, +uint64_t* source_index, long long listLen); +void print_revise_edges(ma_hit_t_alloc* source_paf, uint64_t* source_index, long long listLen); #endif \ No newline at end of file diff --git a/Process_Read.cpp b/Process_Read.cpp index 32101b6..710752c 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -329,6 +329,27 @@ void destory_UC_Read(UC_Read* r) free(r->seq); } +void init_aux_table() +{ + if (bit_t_seq_table[0][0] == 0) + { + uint64_t i = 0; + + for (i = 0; i < 256; i++) + { + bit_t_seq_table[i][0] = s_H[((i >> 6)&(uint64_t)3)]; + bit_t_seq_table[i][1] = s_H[((i >> 4)&(uint64_t)3)]; + bit_t_seq_table[i][2] = s_H[((i >> 2)&(uint64_t)3)]; + bit_t_seq_table[i][3] = s_H[(i&(uint64_t)3)]; + + bit_t_seq_table_rc[i][0] = RC_CHAR(bit_t_seq_table[i][3]); + bit_t_seq_table_rc[i][1] = RC_CHAR(bit_t_seq_table[i][2]); + bit_t_seq_table_rc[i][2] = RC_CHAR(bit_t_seq_table[i][1]); + bit_t_seq_table_rc[i][3] = RC_CHAR(bit_t_seq_table[i][0]); + } + + } +} void init_UC_Read(UC_Read* r) { diff --git a/Process_Read.h b/Process_Read.h index fdc6b3a..a9426ec 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -56,7 +56,7 @@ static char rc_Table[5] = {'T', 'G', 'C', 'A', 'N'}; #define RC_CHAR(x) rc_Table[seq_nt6_table[(uint8_t)x]] - +void init_aux_table(); void init_kseq(char* file); void destory_kseq(); int get_read(kseq_t *s, int adapterLen);