diff --git a/Assembly.cpp b/Assembly.cpp index f707edd..addab11 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1919,8 +1919,9 @@ void* Overlap_calculate_heap_merge(void* arg) for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num) { - get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); - + ///get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); + get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02); + clear_Cigar_record(¤t_cigar); clear_Round2_alignment(&second_round); @@ -2119,8 +2120,8 @@ void* Output_related_reads(void* arg) && memcmp(required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { - get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); - + ////get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); + get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02); fprintf(stderr, ">%.*s\n", Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i)); @@ -2998,7 +2999,8 @@ void* Final_overlap_calculate_heap_merge(void* arg) - get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); + ////get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); + get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.001); /** correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon, &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1, @@ -3081,22 +3083,6 @@ void* Final_overlap_calculate_heap_merge(void* arg) - // if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i), - // Get_NAME_LENGTH((R_INF), i)) == 0) - // { - // fprintf(stderr, "\n2\n"); - // debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "second"); - // } - - // if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i), - // Get_NAME_LENGTH((R_INF), i)) == 0) - // { - // fprintf(stderr, "\n2\n"); - // debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "second"); - // } - - - inner_j = 0; j = 0; @@ -3125,22 +3111,6 @@ void* Final_overlap_calculate_heap_merge(void* arg) } - - // if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i), - // Get_NAME_LENGTH((R_INF), i)) == 0) - // { - // fprintf(stderr, "\n3\n"); - // debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "third"); - // } - - // if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i), - // Get_NAME_LENGTH((R_INF), i)) == 0) - // { - // fprintf(stderr, "\n3\n"); - // debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "third"); - // } - - ///recover missing exact overlaps reverse_complement(g_read.seq, g_read.length); for (j = 0; j < overlap_list.length; j++) @@ -3171,23 +3141,6 @@ void* Final_overlap_calculate_heap_merge(void* arg) } - // if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i), - // Get_NAME_LENGTH((R_INF), i)) == 0) - // { - // fprintf(stderr, "\n4\n"); - // debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "fourth"); - // } - - // if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i), - // Get_NAME_LENGTH((R_INF), i)) == 0) - // { - // fprintf(stderr, "\n4\n"); - // debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "fourth"); - // } - - - - if(R_INF.paf[i].is_fully_corrected) { potiental_matched_overlap_1++; @@ -3201,25 +3154,6 @@ void* Final_overlap_calculate_heap_merge(void* arg) push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &overlap_list, &g_read, &overlap_read, 2, 0); - - // if(memcmp("m64016_190918_162737/130811282/ccs", Get_NAME((R_INF), i), - // Get_NAME_LENGTH((R_INF), i)) == 0) - // { - // fprintf(stderr, "\n5\n"); - // debug_print_overlap("m64016_190918_162737/179635219/ccs", &overlap_list, &R_INF, "last"); - - // debug_print_overlap(NULL, &overlap_list, &R_INF, "SET"); - // } - - // if(memcmp("m64016_190918_162737/179635219/ccs", Get_NAME((R_INF), i), - // Get_NAME_LENGTH((R_INF), i)) == 0) - // { - // fprintf(stderr, "\n5\n"); - // debug_print_overlap("m64016_190918_162737/130811282/ccs", &overlap_list, &R_INF, "last"); - - // debug_print_overlap(NULL, &overlap_list, &R_INF, "SET"); - // } - } diff --git a/Correct.cpp b/Correct.cpp index 746ad3d..e3a7b36 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -537,19 +537,12 @@ int* r_extra_begin, int* r_extra_end, long long* r_y_start, long long* r_y_lengt currentIDLen = Get_READ_LENGTH((*R_INF), y_ID); ///since Window_Len == x_len + (threshold << 1) - if(currentIDLen <= y_start || + if(y_start < 0 || currentIDLen <= y_start || currentIDLen - y_start + 2 * threshold + THRESHOLD_MAX_SIZE < Window_Len) { return 0; } - /** - if(currentIDLen <= y_start) - { - return 0; - } - **/ - extra_begin = extra_end = 0; ///y maybe less than 0 y_start = y_start - threshold; @@ -1054,6 +1047,71 @@ UC_Read* g_read) return T_error; } +inline int double_error_threshold(int pre_threshold, int x_len) +{ + + pre_threshold = Adjust_Threshold(pre_threshold, x_len); + int threshold = pre_threshold * 2; + ///may have some bugs + if(x_len >= 300 && threshold < THRESHOLD_MAX_SIZE) + { + threshold = THRESHOLD_MAX_SIZE; + } + + if(threshold > THRESHOLD_MAX_SIZE) + { + threshold = THRESHOLD_MAX_SIZE; + } + + return threshold; +} + + +inline int verify_sub_window(All_reads* R_INF, Correct_dumy* dumy, UC_Read* g_read, +long long x_beg, long long xLen, long long y_beg, long long yLen, uint64_t y_id, +uint64_t y_pos_strand, int threshold, int alignment_strand, +unsigned int* get_error, int* get_y_end, int* get_x_end, int* get_aligned_xLen) +{ + (*get_aligned_xLen) = 0; + (*get_y_end) = -1; + (*get_x_end) = -1; + (*get_error) = (unsigned int)-1; + + int extra_begin, extra_end, r_x_end, r_y_end, aligned_xLen; + long long o_len; + unsigned int r_error; + if(!determine_overlap_region(threshold, y_beg, y_id, yLen, + R_INF, &extra_begin, &extra_end, &y_beg, &o_len)) + { + return 0; + } + + fill_subregion(dumy->overlap_region, y_beg, o_len, y_pos_strand, R_INF, + y_id, extra_begin, extra_end); + + char* x_string = g_read->seq + x_beg; + char* y_string = dumy->overlap_region; + + aligned_xLen = 0; + + alignment_extension(y_string, yLen, x_string, xLen, threshold, + alignment_strand, &r_error, &r_y_end, &r_x_end, &aligned_xLen); + + (*get_error) = r_error; + (*get_y_end) = r_y_end; + (*get_x_end) = r_x_end; + (*get_aligned_xLen) = aligned_xLen; + + if(aligned_xLen == 0) + { + return 0; + } + else + { + return 1; + } +} + inline double non_trim_error_rate(overlap_region_alloc* overlap_list, long long ID, All_reads* R_INF, Correct_dumy* dumy, UC_Read* g_read) { @@ -1076,13 +1134,96 @@ All_reads* R_INF, Correct_dumy* dumy, UC_Read* g_read) } else { - tError += - get_high_error(overlap_list->list[ID].w_list[i].x_start, - overlap_list->list[ID].w_list[i].x_end, overlap_list->list[ID].w_list[i].y_start, - overlap_list->list[ID].w_list[i].y_end, overlap_list->list[ID].y_id, - overlap_list->list[ID].y_pos_strand, - overlap_list->list[ID].w_list[i].error_threshold - overlap_list->list[ID].w_list[i].extra_begin, - 3, 0.3, R_INF, dumy, g_read); + int x_len = subWinLen; + int threshold = double_error_threshold(overlap_list->list[ID].w_list[i].error_threshold, x_len); + int Window_Len = x_len + (threshold << 1); + unsigned int r_error_left; + int r_x_end_left, r_y_end_left, aligned_xLen_left; + unsigned int r_error_right; + int r_x_end_right, r_y_end_right, aligned_xLen_right; + long long y_beg_left, y_beg_right; + + aligned_xLen_left = aligned_xLen_right = 0; + y_beg_left = y_beg_right = -1; + + if(overlap_list->list[ID].w_list[i].y_start == -1) + { + tError += x_len; + continue; + } + + ///if the previous window is mapped + if(i > 0 && overlap_list->list[ID].w_list[i - 1].y_end != -1) + { + y_beg_left = overlap_list->list[ID].w_list[i - 1].y_end + 1; + } + + ///if the next window is mapped + if(i < overlap_list->list[ID].w_list_length - 1 && overlap_list->list[ID].w_list[i + 1].y_end != -1) + { + y_beg_right = 1 + overlap_list->list[ID].w_list[i + 1].y_start - 1 - x_len; + } + + ///cannot apply this strategy to the leftmost and the rightmost window + if(y_beg_left == -1 && y_beg_right == -1 && i > 0 && i < subWinNum - 1) + { + y_beg_left = overlap_list->list[ID].w_list[i].y_start; + if(overlap_list->list[ID].w_list[i].extra_begin >= 0) + { + y_beg_left = y_beg_left + overlap_list->list[ID].w_list[i].error_threshold - + overlap_list->list[ID].w_list[i].extra_begin; + } + y_beg_right = y_beg_left; + } + + + if(y_beg_left != -1) + { + verify_sub_window(R_INF, dumy, g_read, overlap_list->list[ID].w_list[i].x_start, + x_len, y_beg_left, Window_Len, overlap_list->list[ID].y_id, overlap_list->list[ID].y_pos_strand, + threshold, 0, &r_error_left, &r_y_end_left, &r_x_end_left, &aligned_xLen_left); + } + + if(y_beg_right != -1) + { + verify_sub_window(R_INF, dumy, g_read, overlap_list->list[ID].w_list[i].x_start, + x_len, y_beg_right, Window_Len, overlap_list->list[ID].y_id, overlap_list->list[ID].y_pos_strand, + threshold, 1, &r_error_right, &r_y_end_right, &r_x_end_right, &aligned_xLen_right); + } + + ///here one of aligned_xLen_left and aligned_xLen_right must be 0 + if((i == 0 || i == subWinNum - 1) && + (x_len - (aligned_xLen_left + aligned_xLen_right) > 10)) + { + return 1.0; + } + + ///aligned in both direction + if(aligned_xLen_left != 0 && aligned_xLen_right != 0) + { + if(aligned_xLen_left + aligned_xLen_right <= x_len) + { + tError = tError + r_error_left + r_error_right + + (x_len - aligned_xLen_left - aligned_xLen_right); + } + else + { + float E_rate = (float)(x_len)/(float)(aligned_xLen_left + aligned_xLen_right); + tError = tError + (r_error_left + r_error_right)*E_rate; + } + }///not aligned in both direction + else if(aligned_xLen_left == 0 && aligned_xLen_right == 0) + { + tError += x_len; + }///only aligned in left + else if(aligned_xLen_left != 0) + { + tError = tError + r_error_left + (x_len - aligned_xLen_left); + }///only aligned in right + else if(aligned_xLen_right != 0) + { + tError = tError + r_error_right + (x_len - aligned_xLen_right); + } } } @@ -1912,6 +2053,7 @@ int* r_extra_begin, int* r_extra_end, unsigned int* r_error) int path_length; + ///if the start pos at the left boundary if(local_y_start == 0) { total_y_start = total_y_start + local_y_start; @@ -2012,24 +2154,7 @@ int* r_extra_begin, int* r_extra_end, unsigned int* r_error) } -inline int double_error_threshold(int pre_threshold, int x_len) -{ - pre_threshold = Adjust_Threshold(pre_threshold, x_len); - int threshold = pre_threshold * 2; - ///may have some bugs - if(x_len >= 300 && threshold < THRESHOLD_MAX_SIZE) - { - threshold = THRESHOLD_MAX_SIZE; - } - - if(threshold > THRESHOLD_MAX_SIZE) - { - threshold = THRESHOLD_MAX_SIZE; - } - - return threshold; -} inline void recalcate_window_back(overlap_region_alloc* overlap_list, All_reads* R_INF, @@ -2716,7 +2841,7 @@ inline void recalcate_window_back(overlap_region_alloc* overlap_list, All_reads* } -inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_INF, +inline void recalcate_window_simple(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read) { long long j, k, i; @@ -3275,6 +3400,497 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN + +inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_INF, + UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read) +{ + long long j, k, i; + long long Len_x; + int threshold; + long long y_len; + long long currentIDLen; + long long matches; + long long y_id; + int y_strand; + long long y_readLen; + long long x_start; + long long x_end; + long long x_len; + long long total_y_start; + long long total_y_end; + long long y_start; + long long y_end; + long long Window_Len; + char* x_string; + char* y_string; + int end_site; + unsigned int error; + int real_y_start; + long long overlap_length; + int extra_begin, extra_end; + long long o_len; + + + ///j负责遍历整个overlap list + 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; + y_readLen = Get_READ_LENGTH((*R_INF), y_id); + + //i corresponding to each window of a overlap + //utilize the the end pos of pre-window in backwards + for (i = overlap_list->list[j].w_list_length - 1; i >= 0; i--) + { + ///the first matched window + if(overlap_list->list[j].w_list[i].y_end != -1) + { + ///note!!! need notification + ///this is the actual end postion in ystring + total_y_start = overlap_list->list[j].w_list[i].y_end + - overlap_list->list[j].w_list[i].extra_begin + 1; + + + + ///k corresponding to all unmatched windows at the right side of overlap_list->list[j].w_list[i] + ///so k starts from i + 1, and end to the first matched window + for (k = i + 1; k < overlap_list->list[j].w_list_length && overlap_list->list[j].w_list[k].y_end == -1; k++) + { + extra_begin = extra_end = 0; + + ///if y_start > y_readLen, direct terminate + if (total_y_start >= y_readLen) + { + break; + } + + ///there is no problem for x + x_start = overlap_list->list[j].w_list[k].x_start; + x_end = overlap_list->list[j].w_list[k].x_end; + x_len = x_end - x_start + 1; + ///there are two potiential reasons for unmatched window: + ///1. this window has a large number of differences + ///2. DP does not start from the right offset + threshold = double_error_threshold(overlap_list->list[j].w_list[k].error_threshold, x_len); + + y_start = total_y_start; + 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)) + { + break; + } + + if(o_len + threshold < x_len) + { + break; + } + + fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, + R_INF, y_id, extra_begin, extra_end); + + x_string = g_read->seq + x_start; + y_string = dumy->overlap_region; + + ///note!!! need notification + end_site = Reserve_Banded_BPM(y_string, Window_Len, x_string, x_len, threshold, &error); + + + + // unsigned int debug_error; + // int return_y_end, return_x_end; + // Reserve_Banded_BPM_Extension(y_string, Window_Len, x_string, x_len, threshold, + // &debug_error, &return_y_end, &return_x_end); + // if(return_x_end != x_len - 1 && end_site != -1) + // { + // fprintf(stderr, "error\n"); + // } + + // if(return_x_end == x_len - 1 && end_site == -1) + // { + // fprintf(stderr, "error\n"); + // } + + // if(return_x_end == x_len - 1 && end_site != -1) + // { + // if(return_y_end != end_site || debug_error != error) + // { + // fprintf(stderr, "error\n"); + // } + // } + + + + ///error等于-1说明没匹配 + if (error!=(unsigned int)-1) + { + overlap_list->list[j].w_list[k].cigar.length = -1; + overlap_list->list[j].w_list[k].y_start = y_start; + overlap_list->list[j].w_list[k].y_end = y_start + end_site; + overlap_list->list[j].w_list[k].error = (int)error; + ///note!!! need notification + overlap_list->list[j].w_list[k].extra_begin = extra_begin; + overlap_list->list[j].w_list[k].extra_end = extra_end; + overlap_list->list[j].w_list[k].error_threshold = threshold; + + overlap_list->list[j].align_length += x_len; + } + else + { + break; + } + + ///note!!! need notification + total_y_start = y_start + end_site - extra_begin + 1; + } + + } + + } + + + + //i corresponding to each window of a overlap + //utilize the the start pos of next window in forward + for (i = 0; i < overlap_list->list[j].w_list_length; i++) + { + ///find the first matched window, which should not be the first window + ///the pre-window of this matched window must be unmatched + if(overlap_list->list[j].w_list[i].y_end != -1 && i != 0 && overlap_list->list[j].w_list[i - 1].y_end == -1) + { + ///check if the start pos of this matched window has been calculated + if(overlap_list->list[j].w_list[i].cigar.length == -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; + /****************************may have bugs********************************/ + threshold = overlap_list->list[j].w_list[i].error_threshold; + /****************************may have bugs********************************/ + /****************************may have bugs********************************/ + ///should not adjust threshold, since this window can be matched by the old threshold + ///threshold = Adjust_Threshold(threshold, x_len); + /****************************may have bugs********************************/ + Window_Len = x_len + (threshold << 1); + + + ///y_start is the real y_start + y_start = overlap_list->list[j].w_list[i].y_start; + extra_begin = overlap_list->list[j].w_list[i].extra_begin; + extra_end = overlap_list->list[j].w_list[i].extra_end; + o_len = Window_Len - extra_end - extra_begin; + fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, + R_INF, y_id, extra_begin, extra_end); + x_string = g_read->seq + x_start; + y_string = dumy->overlap_region; + + ///note!!! need notification + end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start, + &(dumy->path_length), dumy->matrix_bit, dumy->path, + overlap_list->list[j].w_list[i].error, overlap_list->list[j].w_list[i].y_end - y_start); + + + ///y_start has already been calculated + if (error != (unsigned int)-1) + { + ///this condition is always wrong + ///in best case, real_y_start = threshold, end_site = Window_Len - threshold - 1 + if (end_site == Window_Len - 1 || real_y_start == 0) + { + if(fix_boundary(x_string, x_len, threshold, y_start, real_y_start, + end_site, extra_begin, extra_end, y_id, Window_Len, R_INF, dumy, + y_strand, error, &y_start, &real_y_start, &end_site, &extra_begin, + &extra_end, &error)) + { + overlap_list->list[j].w_list[i].error = error; + overlap_list->list[j].w_list[i].extra_begin = extra_begin; + overlap_list->list[j].w_list[i].extra_end = extra_end; + } + } + + + + generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[i]), + &real_y_start, &end_site, &error, x_string, x_len, y_string); + + // if(real_y_start < extra_begin || end_site >= Window_Len - extra_end) + // { + // fprintf(stderr, "\nreal_y_start: %d, extra_begin: %d\n", + // real_y_start, extra_begin); + + // fprintf(stderr, "end_site: %d, Window_Len: %d, extra_end: %d\n", + // end_site, Window_Len, extra_end); + // } + + + ///note!!! need notification + real_y_start = y_start + real_y_start - extra_begin; + overlap_list->list[j].w_list[i].y_start = real_y_start; + ///I forget why don't reduce the extra_begin for y_end + ///it seems extra_begin will be reduced at the end of this function + overlap_list->list[j].w_list[i].y_end = y_start + end_site; + overlap_list->list[j].w_list[i].error = error; + } + else + { + fprintf(stderr, "error\n"); + } + } + else + { + real_y_start = overlap_list->list[j].w_list[i].y_start; + } + + + + ///the end pos for pre window is real_y_start - 1 + total_y_end = real_y_start - 1; + ///k遍历匹配window左侧所有不匹配的window + ///如果i匹配，则k从i-1开始 + ///直到第一个匹配的window结束 + ///因为i!=0，所以k的大小不用担心 + for (k = i - 1; k >= 0 && overlap_list->list[j].w_list[k].y_end == -1; k--) + { + ///there is no problem in x + x_start = overlap_list->list[j].w_list[k].x_start; + x_end = overlap_list->list[j].w_list[k].x_end; + x_len = x_end - x_start + 1; + ///there are two potiential reasons for unmatched window: + ///1. this window has a large number of differences + ///2. DP does not start from the right offset + threshold = double_error_threshold(overlap_list->list[j].w_list[k].error_threshold, x_len); + + Window_Len = x_len + (threshold << 1); + + if(total_y_end <= 0) + { + break; + } + + ///y_start might be less than 0 + y_start = total_y_end - x_len + 1; + if(!determine_overlap_region(threshold, y_start, y_id, Window_Len, R_INF, + &extra_begin, &extra_end, &y_start, &o_len)) + { + break; + } + + if(o_len + threshold < x_len) + { + break; + } + + fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, + R_INF, y_id, extra_begin, extra_end); + x_string = g_read->seq + x_start; + y_string = dumy->overlap_region; + + ///note!!! need notification + end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start, + &(dumy->path_length), dumy->matrix_bit, dumy->path, -1, -1); + + if (error!=(unsigned int)-1) + { + ///this condition is always wrong + ///in best case, real_y_start = threshold, end_site = Window_Len - threshold - 1 + if (end_site == Window_Len - 1 || real_y_start == 0) + { + fix_boundary(x_string, x_len, threshold, y_start, real_y_start, end_site, + extra_begin, extra_end, y_id, Window_Len, R_INF, dumy, y_strand, error, + &y_start, &real_y_start, &end_site, + &extra_begin, &extra_end, &error); + } + + generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[k]), + &real_y_start, &end_site, &error, x_string, x_len, y_string); + + // if(real_y_start < extra_begin || end_site >= Window_Len - extra_end) + // { + // fprintf(stderr, "\nreal_y_start: %d, extra_begin: %d\n", + // real_y_start, extra_begin); + + // fprintf(stderr, "end_site: %d, Window_Len: %d, extra_end: %d\n", + // end_site, Window_Len, extra_end); + // } + + ///y_start has no shift, but y_end has shift + overlap_list->list[j].w_list[k].y_start = y_start + real_y_start - extra_begin; + overlap_list->list[j].w_list[k].y_end = y_start + end_site; + overlap_list->list[j].w_list[k].error = error; + overlap_list->list[j].align_length += x_len; + overlap_list->list[j].w_list[k].extra_begin = extra_begin; + overlap_list->list[j].w_list[k].extra_end = extra_end; + overlap_list->list[j].w_list[k].error_threshold = threshold; + } + else + { + break; + } + + total_y_end = y_start + real_y_start - 1 - extra_begin; + } + } + } + } + + + + overlap_list->mapped_overlaps_length = 0; + + int pre_threshold; + long long tLen, tError; + double error_rate; + 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; + y_readLen = Get_READ_LENGTH((*R_INF), y_id); + overlap_length = overlap_list->list[j].x_pos_e - overlap_list->list[j].x_pos_s + 1; + overlap_list->list[j].is_match = 0; + + + ///only calculate cigar for high quality overlaps + if (overlap_length * OVERLAP_THRESHOLD_FILTER <= overlap_list->list[j].align_length) + { + + for (i = 0; i < overlap_list->list[j].w_list_length; i++) + { + ///first we need to check if this window is matched + if(overlap_list->list[j].w_list[i].y_end != -1) + { + ///second check if the cigar of this window has been got + if(overlap_list->list[j].w_list[i].cigar.length == -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; + /****************************may have bugs********************************/ + ///threshold = x_len * THRESHOLD_RATE; + threshold = overlap_list->list[j].w_list[i].error_threshold; + /****************************may have bugs********************************/ + /****************************may have bugs********************************/ + ///should not adjust threshold, since this window can be matched by the old threshold + ///threshold = Adjust_Threshold(threshold, x_len); + /****************************may have bugs********************************/ + Window_Len = x_len + (threshold << 1); + + + ///y_start is the real y_start + ///for the window with cigar, y_start has already reduced extra_begin + y_start = overlap_list->list[j].w_list[i].y_start; + extra_begin = overlap_list->list[j].w_list[i].extra_begin; + extra_end = overlap_list->list[j].w_list[i].extra_end; + o_len = Window_Len - extra_end - extra_begin; + fill_subregion(dumy->overlap_region, y_start, o_len, y_strand, + R_INF, y_id, extra_begin, extra_end); + x_string = g_read->seq + x_start; + y_string = dumy->overlap_region; + + + ///note!!! need notification + end_site = Reserve_Banded_BPM_PATH(y_string, Window_Len, x_string, x_len, threshold, &error, &real_y_start, + &(dumy->path_length), dumy->matrix_bit, dumy->path, + overlap_list->list[j].w_list[i].error, overlap_list->list[j].w_list[i].y_end - y_start); + + // if(error != overlap_list->list[j].w_list[i].error) + // { + // fprintf(stderr, "error\n"); + // } + + if (error != (unsigned int)-1) + { + if (end_site == Window_Len - 1 || real_y_start == 0) + { + + if(fix_boundary(x_string, x_len, threshold, y_start, real_y_start, end_site, + extra_begin, extra_end, y_id, Window_Len, R_INF, dumy, y_strand, error, + &y_start, &real_y_start, &end_site, + &extra_begin, &extra_end, &error)) + { + overlap_list->list[j].w_list[i].error = error; + overlap_list->list[j].w_list[i].extra_begin = extra_begin; + overlap_list->list[j].w_list[i].extra_end = extra_end; + } + + } + + generate_cigar(dumy->path, dumy->path_length, &(overlap_list->list[j].w_list[i]), + &real_y_start, &end_site, &error, x_string, x_len, y_string); + + + // if(real_y_start < extra_begin || end_site >= Window_Len - extra_end) + // { + // fprintf(stderr, "\nreal_y_start: %d, extra_begin: %d\n", + // real_y_start, extra_begin); + + // fprintf(stderr, "end_site: %d, Window_Len: %d, extra_end: %d\n", + // end_site, Window_Len, extra_end); + // } + + ///note!!! need notification + real_y_start = y_start + real_y_start - extra_begin; + overlap_list->list[j].w_list[i].y_start = real_y_start; + overlap_list->list[j].w_list[i].y_end = y_start + end_site - extra_begin; + overlap_list->list[j].w_list[i].error = error; + } + else + { + fprintf(stderr, "error\n"); + } + } + else + { + overlap_list->list[j].w_list[i].y_end -= overlap_list->list[j].w_list[i].extra_begin; + } + + + } + } + + ///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) + { + overlap_list->mapped_overlaps_length += overlap_length; + overlap_list->list[j].is_match = 1; + } + else if(error_rate <= 0.045) + { + overlap_list->list[j].is_match = 3; + } + } + } +} + +/** +inline void adjust_alignment_windows(overlap_region* overlap, All_reads* R_INF) +{ + long long y_id = overlap->y_id; + long long y_strand = overlap->y_pos_strand; + long long y_readLen = Get_READ_LENGTH((*R_INF), y_id); + long long overlap_length = overlap->x_pos_e - overlap->x_pos_s + 1; + long long i; + ///note here we start from i = 1, instead of i = 0 + for (i = 1; i < overlap->w_list_length; i++) + { + ///that means we have both start pos and end pos + if(overlap->w_list[i].y_end != -1 && overlap->w_list[i].cigar.length != -1) + { + ; + } + } +} +**/ + + + inline void add_base_to_correct_read_directly(Correct_dumy* dumy, char base) { diff --git a/Hash_Table.cpp b/Hash_Table.cpp index d1d6036..1b40777 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -632,6 +632,225 @@ void append_overlap_region_alloc_from_existing(overlap_region_alloc* list, overl } +int append_inexact_overlap_region_alloc_back(overlap_region_alloc* list, overlap_region* tmp, All_reads* R_INF) +{ + + if (list->length + 1 > list->size) + { + list->size = list->size * 2; + list->list = (overlap_region*)realloc(list->list, sizeof(overlap_region)*list->size); + ///新分配空间要初始化 + memset(list->list + (list->size/2), 0, sizeof(overlap_region)*(list->size/2)); + } + + if (list->length!=0 && + list->list[list->length - 1].y_id==tmp->y_id + ) + { + if(list->list[list->length - 1].shared_seed >= tmp->shared_seed) + { + return 0; + } + else + { + list->length--; + } + } + + if(tmp->x_pos_s <= tmp->y_pos_s) + { + tmp->y_pos_s = tmp->y_pos_s - tmp->x_pos_s; + tmp->x_pos_s = 0; + } + else + { + tmp->x_pos_s = tmp->x_pos_s - tmp->y_pos_s; + tmp->y_pos_s = 0; + } + + + long long x_right_length = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_e - 1; + long long y_right_length = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_e - 1; + + + + if(x_right_length <= y_right_length) + { + tmp->x_pos_e = Get_READ_LENGTH((*R_INF), tmp->x_id) - 1; + tmp->y_pos_e = tmp->y_pos_e + x_right_length; + } + else + { + tmp->x_pos_e = tmp->x_pos_e + y_right_length; + tmp->y_pos_e = Get_READ_LENGTH((*R_INF), tmp->y_id) - 1; + } + + + if (tmp->x_pos_strand == 1) + { + list->list[list->length].x_id = tmp->x_id; + list->list[list->length].x_pos_e = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_s - 1; + list->list[list->length].x_pos_s = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_e - 1; + list->list[list->length].x_pos_strand = 0; + + list->list[list->length].y_id = tmp->y_id; + list->list[list->length].y_pos_e = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_s - 1; + list->list[list->length].y_pos_s = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_e - 1; + list->list[list->length].y_pos_strand = 1; + + + + + resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2)); + add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0); + long long distance_gap; + long long pre_distance_gap = 0; + long long i = 0; + for (i = 0; i < tmp->f_cigar.length; i++) + { + distance_gap = get_fake_gap_shift(&(tmp->f_cigar), i); + if(distance_gap != pre_distance_gap) + { + pre_distance_gap = distance_gap; + add_fake_cigar(&(list->list[list->length].f_cigar), + Get_READ_LENGTH((*R_INF), tmp->x_id) - get_fake_gap_pos(&(tmp->f_cigar), i) - 1, + pre_distance_gap); + } + } + + if(get_fake_gap_pos(&(list->list[list->length].f_cigar), + list->list[list->length].f_cigar.length - 1) != list->list[list->length].x_pos_e) + { + add_fake_cigar(&(list->list[list->length].f_cigar), + list->list[list->length].x_pos_e, + get_fake_gap_shift(&(list->list[list->length].f_cigar), + list->list[list->length].f_cigar.length - 1)); + } + + + /******************************for debug********************************/ + // long long distance_self_pos = tmp->x_pos_e - tmp->x_pos_s; + // long long distance_pos = tmp->y_pos_e - tmp->y_pos_s; + // distance_gap = distance_pos - distance_self_pos; + + // if(distance_gap != + // get_fake_gap_shift(&(list->list[list->length].f_cigar), + // list->list[list->length].f_cigar.length - 1)) + // { + // fprintf(stderr, "error\n"); + // } + // long long distance_self_pos = tmp->x_pos_e - tmp->x_pos_s; + // long long distance_pos = tmp->y_pos_e - tmp->y_pos_s; + // distance_gap = distance_pos - distance_self_pos; + // fprintf(stderr, "\nx_s: %d, x_e: %d, y_s: %d, y_e: %d, distance_gap: %d, xLen: %d\n", + // list->list[list->length].x_pos_s, list->list[list->length].x_pos_e, + // list->list[list->length].y_pos_s, list->list[list->length].y_pos_e, distance_gap, + // Get_READ_LENGTH((*R_INF), tmp->x_id)); + // for (i = 0; i < list->list[list->length].f_cigar.length; i++) + // { + // fprintf(stderr, "##i: %d, gap_pos_in_x: %d, gap_shift: %d\n", + // i, get_fake_gap_pos(&(list->list[list->length].f_cigar), i), + // get_fake_gap_shift(&(list->list[list->length].f_cigar), i)); + // } + // for (i = 0; i < tmp->f_cigar.length; i++) + // { + // fprintf(stderr, "**i: %d, gap_pos_in_x: %d, gap_shift: %d\n", + // i, get_fake_gap_pos(&(tmp->f_cigar), i), + // get_fake_gap_shift(&(tmp->f_cigar), i)); + // } + /******************************for debug********************************/ + } + else + { + list->list[list->length].x_id = tmp->x_id; + list->list[list->length].x_pos_e = tmp->x_pos_e; + list->list[list->length].x_pos_s = tmp->x_pos_s; + list->list[list->length].x_pos_strand = tmp->x_pos_strand; + + list->list[list->length].y_id = tmp->y_id; + list->list[list->length].y_pos_e = tmp->y_pos_e; + list->list[list->length].y_pos_s = tmp->y_pos_s; + list->list[list->length].y_pos_strand = tmp->y_pos_strand; + + + + resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2)); + add_fake_cigar(&(list->list[list->length].f_cigar), list->list[list->length].x_pos_s, 0); + long long distance_self_pos = tmp->x_pos_e - tmp->x_pos_s; + long long distance_pos = tmp->y_pos_e - tmp->y_pos_s; + long long init_distance_gap = distance_pos - distance_self_pos; + long long pre_distance_gap = init_distance_gap; + long long distance_gap; + long long i = 0; + for (i = tmp->f_cigar.length - 1; i >= 0; i--) + { + distance_gap = get_fake_gap_shift(&(tmp->f_cigar), i); + if(distance_gap != pre_distance_gap) + { + pre_distance_gap = distance_gap; + + add_fake_cigar(&(list->list[list->length].f_cigar), + get_fake_gap_pos(&(tmp->f_cigar), i), init_distance_gap - pre_distance_gap); + } + } + + if(get_fake_gap_pos(&(list->list[list->length].f_cigar), + list->list[list->length].f_cigar.length - 1) != list->list[list->length].x_pos_e) + { + add_fake_cigar(&(list->list[list->length].f_cigar), + list->list[list->length].x_pos_e, + get_fake_gap_shift(&(list->list[list->length].f_cigar), + list->list[list->length].f_cigar.length - 1)); + } + + + + + /******************************for debug********************************/ + // distance_self_pos = tmp->x_pos_e - tmp->x_pos_s; + // distance_pos = tmp->y_pos_e - tmp->y_pos_s; + // distance_gap = distance_pos - distance_self_pos; + + // if(distance_gap != + // get_fake_gap_shift(&(list->list[list->length].f_cigar), + // list->list[list->length].f_cigar.length - 1)) + // { + // fprintf(stderr, "error\n"); + // } + // distance_self_pos = tmp->x_pos_e - tmp->x_pos_s; + // distance_pos = tmp->y_pos_e - tmp->y_pos_s; + // distance_gap = distance_pos - distance_self_pos; + // fprintf(stderr, "\nx_s: %d, x_e: %d, y_s: %d, y_e: %d, distance_gap: %d\n", + // list->list[list->length].x_pos_s, list->list[list->length].x_pos_e, + // list->list[list->length].y_pos_s, list->list[list->length].y_pos_e, distance_gap); + // for (i = 0; i < list->list[list->length].f_cigar.length; i++) + // { + // fprintf(stderr, "##i: %d, gap_pos_in_x: %d, gap_shift: %d\n", + // i, get_fake_gap_pos(&(list->list[list->length].f_cigar), i), + // get_fake_gap_shift(&(list->list[list->length].f_cigar), i)); + // } + // for (i = 0; i < tmp->f_cigar.length; i++) + // { + // fprintf(stderr, "**i: %d, gap_pos_in_x: %d, gap_shift: %d\n", + // i, get_fake_gap_pos(&(tmp->f_cigar), i), + // get_fake_gap_shift(&(tmp->f_cigar), i)); + // } + /******************************for debug********************************/ + } + + + list->list[list->length].shared_seed = tmp->shared_seed; + list->list[list->length].align_length = 0; + list->list[list->length].is_match = 0; + list->list[list->length].non_homopolymer_errors = 0; + list->list[list->length].strong = 0; + + list->length++; + + return 1; +} + + int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, All_reads* R_INF) { @@ -1278,6 +1497,164 @@ long long y_beg, long long y_end, long long yLen) return x_end - x_beg + 1; } +void chain_DP_back(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, +double band_width_threshold) +{ + long long i, j; + long long self_pos, pos, max_j, max_i, max_score, score, n_skip; + long long distance_pos, distance_self_pos, distance_gap, log_distance_gap, distance_min; + ///double band_width_threshold = 0.05; + double band_width_penalty = 1 / band_width_threshold; + long long min_score = k_mer_length; + long long max_indels, max_self_length; + double gap_rate; + long long total_indels, total_self_length; + + resize_Chain_Data(dp, a_n); + // fill the score and backtrack arrays + for (i = 0; i < a_n; ++i) + { + pos = a[i].offset; + self_pos = a[i].self_offset; + max_j = -1; + max_score = min_score; + n_skip = 0; + max_indels = 0; + max_self_length = 0; + + + ///may have a pre-cut condition for j + for (j = i - 1; j >= 0; --j) + { + distance_pos = pos - a[j].offset; + distance_self_pos = self_pos - a[j].self_offset; + ///a has been sorted by a[].offset + ///note for a, we do not have any two elements that have both equal offsets and self_offsets + ///but there maybe two elements that have equal offsets or equal self_offsets + if(distance_pos == 0 || distance_self_pos <= 0) + { + continue; + } + + distance_gap = distance_pos > distance_self_pos? distance_pos - distance_self_pos : distance_self_pos - distance_pos; + + total_indels = dp->indels[j] + distance_gap; + total_self_length = dp->self_length[j] + distance_self_pos; + if(total_indels > band_width_threshold * total_self_length) + { + continue; + } + + distance_min = distance_pos < distance_self_pos? distance_pos:distance_self_pos; + score = distance_min < min_score? distance_min : min_score; + + /** + log_distance_gap = distance_gap? ilog2_32(distance_gap) : 0; + score -= (long long)(distance_gap * 0.01 * min_score) + (log_distance_gap/2); + **/ + gap_rate = (double)((double)(total_indels)/(double)(total_self_length)); + ///if the gap rate > 0.05, score will be negative + score -= (long long)(gap_rate * score * band_width_penalty); + + score += dp->score[j]; + + if(score > max_score) + { + max_score = score; + max_j = j; + max_indels = total_indels; + max_self_length = total_self_length; + if (n_skip > 0) + { + n_skip--; + } + } + } + + dp->score[i] = max_score; + dp->pre[i] = max_j; + dp->indels[i] = max_indels; + dp->self_length[i] = max_self_length; + } + + + ///debug_chain(a, a_n, dp); + + + + max_score = -1; + max_i = -1; + for (i = 0; i < a_n; ++i) + { + if(dp->score[i] > max_score) + { + max_score = dp->score[i]; + max_i = i; + } + } + + + clear_fake_cigar(&(result->f_cigar)); + + i = max_i; + result->x_pos_e = a[i].self_offset; + result->y_pos_e = a[i].offset; + result->shared_seed = max_score; + + distance_self_pos = result->x_pos_e - a[i].self_offset; + distance_pos = result->y_pos_e - a[i].offset; + long long pre_distance_gap = distance_pos - distance_self_pos; + ///record first site + ///the length of f_cigar should be at least 1 + add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap); + long long chainLen = 0; + if(result->x_pos_strand == 1) + { + while (i >= 0) + { + distance_self_pos = result->x_pos_e - a[i].self_offset; + distance_pos = result->y_pos_e - a[i].offset; + distance_gap = distance_pos - distance_self_pos; + if(distance_gap != pre_distance_gap) + { + pre_distance_gap = distance_gap; + ///record this site + add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap); + } + + chainLen++; + result->x_pos_s = a[i].self_offset; + result->y_pos_s = a[i].offset; + i = dp->pre[i]; + } + } + else + { + + while (i >= 0) + { + distance_self_pos = result->x_pos_e - a[i].self_offset; + distance_pos = result->y_pos_e - a[i].offset; + distance_gap = distance_pos - distance_self_pos; + if(distance_gap == pre_distance_gap) + { + result->f_cigar.length--; + add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap); + } + else + { + pre_distance_gap = distance_gap; + add_fake_cigar(&(result->f_cigar), a[i].self_offset, pre_distance_gap); + } + + chainLen++; + result->x_pos_s = a[i].self_offset; + result->y_pos_s = a[i].offset; + i = dp->pre[i]; + } + } +} + ///double band_width_threshold = 0.05; void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, double band_width_threshold, int max_skip, int x_readLen, int y_readLen) @@ -1618,95 +1995,20 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh continue; } - - // if( - // ((memcmp("m64016_190918_162737/130811282/ccs", - // Get_NAME((*R_INF), tmp_region.x_id), Get_NAME_LENGTH((*R_INF), tmp_region.x_id)) == 0) - // && - // (memcmp("m64016_190918_162737/179635219/ccs", - // Get_NAME((*R_INF), tmp_region.y_id), Get_NAME_LENGTH((*R_INF), tmp_region.y_id)) == 0)) - // || - // ((memcmp("m64016_190918_162737/179635219/ccs", - // Get_NAME((*R_INF), tmp_region.x_id), Get_NAME_LENGTH((*R_INF), tmp_region.x_id)) == 0) - // && - // (memcmp("m64016_190918_162737/130811282/ccs", - // Get_NAME((*R_INF), tmp_region.y_id), Get_NAME_LENGTH((*R_INF), tmp_region.y_id)) == 0))) - // { - // fprintf(stderr, "****************x_name: %.*s****************\n", - // Get_NAME_LENGTH((*R_INF), tmp_region.x_id), Get_NAME((*R_INF), tmp_region.x_id)); - // fprintf(stderr, "****************y_name: %.*s****************\n", - // Get_NAME_LENGTH((*R_INF), tmp_region.y_id), Get_NAME((*R_INF), tmp_region.y_id)); - - - // k_mer_hit* k_list = candidates->list + sub_region_beg; - // long long k_listLen = sub_region_end - sub_region_beg + 1; - // long long k = 0; - - // fprintf(stderr, "k_listLen: %d\n", k_listLen); - - // for (k = 0; k < k_listLen; k++) - // { - // fprintf(stderr, "k: %d, readID: %d, strand: %d, offset: %d, self_offset: %d\n", - // k, k_list[k].readID, k_list[k].strand, k_list[k].offset, k_list[k].self_offset); - // } - // } - - + chain_DP(candidates->list + sub_region_beg, sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold, 50, Get_READ_LENGTH((*R_INF), tmp_region.x_id), Get_READ_LENGTH((*R_INF), tmp_region.y_id)); + // chain_DP_back(candidates->list + sub_region_beg, + // sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold); - // if( - // ((memcmp("m64016_190918_162737/130811282/ccs", - // Get_NAME((*R_INF), tmp_region.x_id), Get_NAME_LENGTH((*R_INF), tmp_region.x_id)) == 0) - // && - // (memcmp("m64016_190918_162737/179635219/ccs", - // Get_NAME((*R_INF), tmp_region.y_id), Get_NAME_LENGTH((*R_INF), tmp_region.y_id)) == 0)) - // || - // ((memcmp("m64016_190918_162737/179635219/ccs", - // Get_NAME((*R_INF), tmp_region.x_id), Get_NAME_LENGTH((*R_INF), tmp_region.x_id)) == 0) - // && - // (memcmp("m64016_190918_162737/130811282/ccs", - // Get_NAME((*R_INF), tmp_region.y_id), Get_NAME_LENGTH((*R_INF), tmp_region.y_id)) == 0))) - // { - // fprintf(stderr, "****************x_name: %.*s****************\n", - // Get_NAME_LENGTH((*R_INF), tmp_region.x_id), Get_NAME((*R_INF), tmp_region.x_id)); - // fprintf(stderr, "****************y_name: %.*s****************\n", - // Get_NAME_LENGTH((*R_INF), tmp_region.y_id), Get_NAME((*R_INF), tmp_region.y_id)); - - - // k_mer_hit* k_list = candidates->list + sub_region_beg; - // long long k_listLen = sub_region_end - sub_region_beg + 1; - // long long k = 0; - - // fprintf(stderr, "k_listLen: %d\n", k_listLen); - - // for (k = 0; k < k_listLen; k++) - // { - // fprintf(stderr, "k: %d, score: %d, pre: %d, indels: %d, self_length: %d\n", - // k, candidates->chainDP.score[k], candidates->chainDP.pre[k], - // candidates->chainDP.indels[k], candidates->chainDP.self_length[k]); - // } - // } - - - - ///自己和自己重叠的要排除 ///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1) if (tmp_region.x_id != tmp_region.y_id) { - ///append_overlap_region_alloc(overlap_list, &tmp_region, R_INF); - /** - if(append_inexact_overlap_region_alloc(overlap_list, &tmp_region, R_INF)) - { - debug_seed_offset(candidates->list + sub_region_beg, - sub_region_end - sub_region_beg + 1, - &(candidates->chainDP), &(overlap_list->list[overlap_list->length - 1]), R_INF); - } - **/ append_inexact_overlap_region_alloc(overlap_list, &tmp_region, R_INF); + ///append_inexact_overlap_region_alloc_back(overlap_list, &tmp_region, R_INF); } } diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h index ea0bd0e..b783b53 100644 --- a/Levenshtein_distance.h +++ b/Levenshtein_distance.h @@ -13,6 +13,394 @@ typedef uint64_t Word; typedef uint32_t Word_32; +typedef struct +{ + /** + * vec[num_words-1], vec[num_words-2], ..., vec[1], vec[0] + * num_bits-1, num_bits-2, num_bits-3, ..., 8, 7, 6, 5, 4, 3, 2, 1, 0 + **/ + Word* vec; + Word num_words; + Word num_bits; + ///for vec[num_words-1] + Word last_bit_mode; + Word last_bit_shift; +} bit_vectors; + +inline void init_bit_vector(bit_vectors* b, Word num_bits) +{ + b->num_bits = num_bits; + b->num_words = ((b->num_bits)>>6); + b->last_bit_mode = (Word)-1; + b->last_bit_shift = b->num_bits&63; + if(b->last_bit_shift != 0) + { + b->num_words++; + b->last_bit_mode = b->last_bit_mode >> (64 - b->last_bit_shift); + } + b->vec = (Word*)calloc(b->num_words, sizeof(Word)); +} + +inline void destory_bit_vector(bit_vectors* b, Word num_bits) +{ + free(b->vec); +} + +///& +inline int AND_bit_vector(bit_vectors* x, bit_vectors* y, bit_vectors* dest) +{ + if(x->num_bits != y->num_bits || x->num_bits != dest->num_bits) return 0; + + Word i; + for (i = 0; i < x->num_words; i++) + { + dest->vec[i] = x->vec[i] & y->vec[i]; + } + + if(dest->num_words > 0) + { + dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode; + } + + return 1; +} + +///| +inline int OR_bit_vector(bit_vectors* x, bit_vectors* y, bit_vectors* dest) +{ + if(x->num_bits != y->num_bits || x->num_bits != dest->num_bits) return 0; + + Word i; + for (i = 0; i < x->num_words; i++) + { + dest->vec[i] = x->vec[i] | y->vec[i]; + } + + if(dest->num_words > 0) + { + dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode; + } + + return 1; +} + + +///^ +inline int XOR_bit_vector(bit_vectors* x, bit_vectors* y, bit_vectors* dest) +{ + if(x->num_bits != y->num_bits || x->num_bits != dest->num_bits) return 0; + + Word i; + for (i = 0; i < x->num_words; i++) + { + dest->vec[i] = x->vec[i] ^ y->vec[i]; + } + + if(dest->num_words > 0) + { + dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode; + } + + return 1; +} + +///~ +inline int NOT_bit_vector(bit_vectors* source, bit_vectors* dest) +{ + if(source->num_bits != dest->num_bits) return 0; + + Word i; + for (i = 0; i < source->num_words; i++) + { + dest->vec[i] = ~(source->vec[i]); + } + + if(dest->num_words > 0) + { + dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode; + } + + return 1; +} + +///<< 1 +inline int L_shift_1_bit_vector(bit_vectors* source, bit_vectors* dest) +{ + /** + * vec[num_words-1], vec[num_words-2], ..., vec[1], vec[0] + * num_bits-1, num_bits-2, num_bits-3, ..., 8, 7, 6, 5, 4, 3, 2, 1, 0 + **/ + if(source->num_bits != dest->num_bits || source->num_words < 1) return 0; + Word i; + for (i = source->num_words - 1; i >= 1; i--) + { + dest->vec[i] = (source->vec[i])<<1; + dest->vec[i] = dest->vec[i] | ((source->vec[i-1])>>63); + } + + dest->vec[0] = (source->vec[0])<<1; + dest->vec[dest->num_words - 1] = dest->vec[dest->num_words - 1] & dest->last_bit_mode; + return 1; +} + + +inline void get_error(int t_length, int errthold, int init_err, Word VP, Word VN, +unsigned int* return_err, int* back_site) +{ + (*return_err) = (unsigned int)-1; + int site = t_length - 1; + int return_site = -1; + ///p_length大部分情况下应该是t_length + 2 * errthold，这是i要小于last_high = 2 * errthold + ///也就是p_length - t_length + ///那么当p_length < t_length + 2 * errthold, available_i也应该是这个值 + ///int available_i = p_length - t_length; + int available_i = 2 * errthold; + + + if ((init_err <= errthold) && (init_err <= (*return_err))) + { + (*return_err) = init_err; + return_site = site; + } + + + int i = 0; + unsigned int ungap_error = (unsigned int)-1; + + while (i < available_i) + { + init_err = init_err + ((VP >> i)&(Word)1); + init_err = init_err - ((VN >> i)&(Word)1); + ++i; + + if ((init_err <= errthold) && (init_err <= *return_err)) + { + *return_err = init_err; + return_site = site + i; + } + + /****************************may have bugs********************************/ + if(i == errthold) + { + ungap_error = init_err; + } + /****************************may have bugs********************************/ + } + + /****************************may have bugs********************************/ + if((ungap_error<=errthold) && (ungap_error == (*return_err))) + { + return_site = site + errthold; + } + /****************************may have bugs********************************/ + + (*back_site) = return_site; +} + +inline int Reserve_Banded_BPM_Extension +(char *pattern, int p_length, char *text, int t_length, unsigned short errthold, +unsigned int* return_err, int* return_p_end, int* return_t_end) +{ + (*return_err) = (unsigned int)-1; + (*return_p_end) = -1; + (*return_t_end) = -1; + + Word Peq[256]; + + unsigned int line_error = (unsigned int)-1; + int return_site; + int band_length = (errthold << 1) + 1; + int i = 0; + Word tmp_Peq_1 = (Word)1; + + Peq['A'] = (Word)0; + Peq['T'] = (Word)0; + Peq['G'] = (Word)0; + Peq['C'] = (Word)0; + + + Word Peq_A; + Word Peq_T; + Word Peq_C; + Word Peq_G; + + ///band_length = 2k + 1 + for (i = 0; i> 1; + VN = X&HP; + VP = HN | ~(X | HP); + + if (!(D0&err_mask)) + { + ++err; + if ((err - last_high)>errthold) + { + return (*return_t_end); + } + } + get_error(i + 1, errthold, err, VP, VN, &line_error, &return_site); + if(line_error != (unsigned int)-1) + { + (*return_t_end) = i; + (*return_p_end) = return_site; + (*return_err) = line_error; + } + + Peq['A'] = Peq['A'] >> 1; + Peq['C'] = Peq['C'] >> 1; + Peq['G'] = Peq['G'] >> 1; + Peq['T'] = Peq['T'] >> 1; + + + ++i; + ++i_bd; + Peq[pattern[i_bd]] = Peq[pattern[i_bd]] | Mask; + } + + + + + + X = Peq[text[i]] | VN; + D0 = ((VP + (X&VP)) ^ VP) | X; + HN = VP&D0; + HP = VN | ~(VP | D0); + X = D0 >> 1; + VN = X&HP; + VP = HN | ~(X | HP); + if (!(D0&err_mask)) + { + ++err; + if ((err - last_high)>errthold) + { + return (*return_t_end); + } + } + ///i = t_length - 1 + get_error(i + 1, errthold, err, VP, VN, &line_error, &return_site); + if(line_error != (unsigned int)-1) + { + (*return_t_end) = i; + (*return_p_end) = return_site; + (*return_err) = line_error; + } + + return (*return_t_end); +} + +inline void reverse_string(char* str, int strLen) +{ + int i, Len; + char k; + Len = strLen / 2; + for (i = 0; i < Len; i++) + { + k = str[i]; + str[i] = str[strLen - i - 1]; + str[strLen - i - 1] = k; + } +} + +inline int alignment_extension(char *pattern, int p_length, char *text, int t_length, +unsigned short errthold, int direction, unsigned int* return_err, int* return_p_end, +int* return_t_end, int* return_aligned_t_len) +{ + (*return_aligned_t_len) = 0; + + if(direction == 0) + { + Reserve_Banded_BPM_Extension(pattern, p_length, text, t_length, errthold, return_err, + return_p_end, return_t_end); + if((*return_p_end) != -1 && (*return_t_end) != -1) + { + (*return_aligned_t_len) = (*return_t_end) + 1; + return 1; + } + else + { + return -1; + } + + } + else + { + reverse_string(pattern, p_length); + reverse_string(text, t_length); + + Reserve_Banded_BPM_Extension(pattern, p_length, text, t_length, errthold, return_err, + return_p_end, return_t_end); + + reverse_string(pattern, p_length); + reverse_string(text, t_length); + + if((*return_p_end) != -1 && (*return_t_end) != -1) + { + (*return_aligned_t_len) = (*return_t_end) + 1; + (*return_p_end) = p_length - (*return_p_end); + (*return_t_end) = t_length - (*return_t_end); + return 1; + } + else + { + return -1; + } + } +} + diff --git a/Overlaps.cpp b/Overlaps.cpp index 2fd624e..5862bb1 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -9890,7 +9890,7 @@ char* output_file_name, long long bubble_dist, int read_graph, int write) // debug_info_of_specfic_node("m64016_190918_162737/141297762/ccs", sg); out: - ///output_tips(sg, &R_INF); + output_tips(sg, &R_INF); output_unitig_graph(sg, coverage_cut, output_file_name, n_read);