From 90195293f175b92db9c927da2385adb4430756aa Mon Sep 17 00:00:00 2001 From: Haoyu Cheng Date: Thu, 12 Dec 2019 20:12:52 -0500 Subject: [PATCH] for backup --- Correct.cpp | 295 +++++++++++++++++++++++++++++++++++++++++++++------- Correct.h | 4 +- 2 files changed, 259 insertions(+), 40 deletions(-) diff --git a/Correct.cpp b/Correct.cpp index 9e04f41..bd2d824 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -7115,7 +7115,6 @@ Cigar_record* current_cigar, long long uncorrected_window_start, Round2_alignmen backbone_length = corrected_window_end - corrected_window_start + 1; addUnmatchedSeqToGraph(g, backbone, backbone_length, &startNodeID, &endNodeID); - ///fprintf(stderr, "startNodeID: %d, endNodeID: %d\n", startNodeID, endNodeID); long long correct_x_pos_s; long long matched_coverage = 0; @@ -7125,7 +7124,7 @@ Cigar_record* current_cigar, long long uncorrected_window_start, Round2_alignmen correct_x_pos_s = (overlap_list->list[overlapID].x_pos_s / WINDOW) * WINDOW; windowID = (uncorrected_window_start - correct_x_pos_s) / WINDOW; - ///如果这个window不匹配，跳过 + ///skip if window is unmatched if (overlap_list->list[overlapID].w_list[windowID].y_end == -1) { continue; @@ -7139,7 +7138,7 @@ Cigar_record* current_cigar, long long uncorrected_window_start, Round2_alignmen /** * There are total 3 cases: - * 1. this window of x is overlapped totally by y + * 1. this window of x is overlapped fully by y * x: ------|------|--------- * y: ------|------|--------- * in this case, x_start == uncorrected_window_start, x_length == WINDOW @@ -7626,24 +7625,26 @@ CIGAR* cigar, haplotype_evdience_alloc* hap) int operation; int operationLen; int i; + ///mismatches based on the offset of x long long inner_offset = x_total_start - window_offset; ///note that node 0 is the start node ///0 is match, 1 is mismatch, 2 is up, 3 is left - ///2是x缺字符（y多字符），而3是y缺字符（x多字符） + ///2 represents thre are more bases at y + ///3 represents thre are more bases at x while (cigar_i < cigar->length) { operation = cigar->C_C[cigar_i]; operationLen = cigar->C_L[cigar_i]; - ///这种情况代表匹配和mismatch + ///match if (operation == 0) { x_i += operationLen; y_i += operationLen; } - else if(operation == 1) + else if(operation == 1) ///mismatch { for (i = 0; i < operationLen; i++) { @@ -7667,8 +7668,6 @@ CIGAR* cigar, haplotype_evdience_alloc* hap) cigar_i++; } - - } @@ -7693,17 +7692,19 @@ CIGAR* cigar, haplotype_evdience_alloc* hap, int snp_threshold) ///note that node 0 is the start node ///0 is match, 1 is mismatch, 2 is up, 3 is left - ///2是x缺字符（y多字符），而3是y缺字符（x多字符） + ///2 represents thre are more bases at y + ///3 represents thre are more bases at x while (cigar_i < cigar->length) { operation = cigar->C_C[cigar_i]; operationLen = cigar->C_L[cigar_i]; - ///这种情况代表匹配和mismatch + ///matches if (operation == 0) { for (i = 0; i < operationLen; i++) { + ///should be at least 2 mismatches if(hap->flag[inner_offset] > snp_threshold) { ev.misBase = y_string[y_i]; @@ -7747,7 +7748,7 @@ CIGAR* cigar, haplotype_evdience_alloc* hap, int snp_threshold) } else if (operation == 3) { - + /****************************may have bugs********************************/ for (i = 0; i < operationLen; i++) { if(hap->flag[inner_offset] > snp_threshold) @@ -7763,6 +7764,7 @@ CIGAR* cigar, haplotype_evdience_alloc* hap, int snp_threshold) inner_offset++; x_i++; } + /****************************may have bugs********************************/ } cigar_i++; @@ -7779,10 +7781,10 @@ CIGAR* cigar, haplotype_evdience_alloc* hap, int snp_threshold) - void cluster(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) { + ///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; long long x_length; @@ -7796,23 +7798,13 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, haplot long long correct_x_pos_s; long long inner_window_offset; int snp_threshold; - /** - if(overlap_list->mapped_overlaps > Coverage_Threshold(coverage)) - { - snp_threshold = 0; - } - else**/ - { - snp_threshold = 1; - } - + snp_threshold = 1; - - ///与当前window重叠的所有overlap + ///all overlaps related to the current window [window_start, window_end] ///first mark all snp pos for (i = 0; i < dumy->length; i++) { - ///这个是那个overlap的ID，而不是overlap里对应窗口的ID + ///overlap id, instead of the window id or the y id overlapID = dumy->overlapID[i]; ///overlap_list->list[overlapID].x_pos_s is the begining of the whole overlap @@ -7820,13 +7812,14 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, haplot ///window_start is the begining of this window in the whole x_read windowID = (window_start - correct_x_pos_s) / WINDOW; - ///如果这个window不匹配，跳过 + ///skip if this window is not matched if (overlap_list->list[overlapID].w_list[windowID].y_end == -1) { continue; } - ///both x_start and y_start are the begining of the whole x_read and y_read + ///both x_start and y_start are the offsets of the whole x_read and y_read + ///instead of the offsets of window 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; @@ -7836,33 +7829,35 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, haplot - overlap_list->list[overlapID].w_list[windowID].y_start + 1; - markSNP(window_start, x_start, x_length, y_start, y_length, &(overlap_list->list[overlapID].w_list[windowID].cigar), - hap); + markSNP(window_start, x_start, x_length, y_start, y_length, + &(overlap_list->list[overlapID].w_list[windowID].cigar), hap); } /****************************may have bugs********************************/ - long long last_snp = useful_length - 1; + long long last_snp = -1; long long first_snp = -1; for (i = 0; i < useful_length; i++) { if(hap->flag[i] != 0) { last_snp = i; - if(first_snp != -1) + if(first_snp == -1) { first_snp = i; } } - + ///for a real snp, the coverage should be at least 2 if(hap->flag[i] > snp_threshold) { hap->snp++; } } - if(first_snp == -1) + ///if there are any >0 elements, both first_snp and last_snp should be != -1 + if(first_snp == -1 || last_snp == -1) { first_snp = 0; + last_snp = -1; } /****************************may have bugs********************************/ @@ -7871,7 +7866,7 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, haplot ///add the information related to snp to haplotype_evdience_alloc for (i = 0; i < dumy->length; i++) { - ///这个是那个overlap的ID，而不是overlap里对应窗口的ID + ///overlap ID, instead of the window ID overlapID = dumy->overlapID[i]; ///overlap_list->list[overlapID].x_pos_s is the begining of the whole overlap @@ -7879,13 +7874,14 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, haplot ///window_start is the begining of this window in the whole x_read windowID = (window_start - correct_x_pos_s) / WINDOW; - ///如果这个window不匹配，跳过 + ///skip if this window is not matched if (overlap_list->list[overlapID].w_list[windowID].y_end == -1) { continue; } - ///both x_start and y_start are the begining of the whole x_read and y_read + ///both x_start and y_start are the offsets of the whole x_read and y_read + ///instead of the offsets of window 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; @@ -7907,7 +7903,133 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, haplot hap, snp_threshold); } - RsetInitHaplotypeEvdienceFlag(hap + first_snp, last_snp + 1); + RsetInitHaplotypeEvdienceFlag(hap, first_snp, last_snp + 1 - first_snp); +} + + +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) +{ + ///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; + long long x_length; + char* x_string; + char* y_string; + long long i; + long long y_start, y_length; + long long overlapID, windowID; + long long startNodeID, endNodeID, currentNodeID; + + long long correct_x_pos_s; + long long inner_window_offset; + int snp_threshold; + snp_threshold = 1; + + ///all overlaps related to the current window [window_start, window_end] + ///first mark all snp pos + for (i = 0; i < dumy->length; i++) + { + ///overlap id, instead of the window id or the y id + overlapID = dumy->overlapID[i]; + + ///overlap_list->list[overlapID].x_pos_s is the begining of the whole overlap + correct_x_pos_s = (overlap_list->list[overlapID].x_pos_s / WINDOW) * WINDOW; + ///window_start is the begining of this window in the whole x_read + windowID = (window_start - correct_x_pos_s) / WINDOW; + + ///skip if this window is not matched + if (overlap_list->list[overlapID].w_list[windowID].y_end == -1) + { + continue; + } + + ///both x_start and y_start are the offsets of the whole x_read and y_read + ///instead of the offsets of window + 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; + + 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; + + + markSNP(window_start, x_start, x_length, y_start, y_length, + &(overlap_list->list[overlapID].w_list[windowID].cigar), hap); + } + + + /****************************may have bugs********************************/ + long long last_snp = -1; + long long first_snp = -1; + for (i = 0; i < useful_length; i++) + { + if(hap->flag[i] != 0) + { + last_snp = i; + if(first_snp == -1) + { + first_snp = i; + } + } + ///for a real snp, the coverage should be at least 2 + if(hap->flag[i] > snp_threshold) + { + hap->snp++; + } + } + ///if there are any >0 elements, both first_snp and last_snp should be != -1 + if(first_snp == -1 || last_snp == -1) + { + first_snp = 0; + last_snp = -1; + } + /****************************may have bugs********************************/ + + + + ///add the information related to snp to haplotype_evdience_alloc + for (i = 0; i < dumy->length; i++) + { + ///overlap ID, instead of the window ID + overlapID = dumy->overlapID[i]; + + ///overlap_list->list[overlapID].x_pos_s is the begining of the whole overlap + correct_x_pos_s = (overlap_list->list[overlapID].x_pos_s / WINDOW) * WINDOW; + ///window_start is the begining of this window in the whole x_read + windowID = (window_start - correct_x_pos_s) / WINDOW; + + ///skip if this window is not matched + if (overlap_list->list[overlapID].w_list[windowID].y_end == -1) + { + continue; + } + + ///both x_start and y_start are the offsets of the whole x_read and y_read + ///instead of the offsets of window + 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; + + 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; + + + 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; + + + 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); + } + + RsetInitHaplotypeEvdienceFlag(hap, first_snp, last_snp + 1 - first_snp); } int cmp_haplotype_evdience(const void * a, const void * b) @@ -11245,7 +11367,103 @@ 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) +{ + ResizeInitHaplotypeEvdience(hap); + long long i, j, overlap_length; + long long window_start, window_end; + + long long num_availiable_win = 0; + + Window_Pool w_inf; + init_Window_Pool(&w_inf, g_read->length, WINDOW, TAIL_LENGTH); + + int flag = 0; + while(get_Window(&w_inf, &window_start, &window_end) && flag != -2) + { + dumy->length = 0; + dumy->lengthNT = 0; + ///flag返回的是重叠数量 + ///dumy->length返回的是有效完全重叠的数量 + ///dumy->lengthNT返回的是有效不完全重叠的数量 + ///return overlaps that is overlaped with [window_start, window_end] + flag = get_available_interval(window_start, window_end, overlap_list, dumy); + switch (flag) + { + case 1: ///找到匹配 + break; + case 0: ///没找到匹配 + break; + case -2: ///下一个window也不会存在匹配, 直接跳出 + break; + } + + ///这个是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); + } + + + + ///very time-consuming + qsort(hap->list, hap->length, sizeof(haplotype_evdience), cmp_haplotype_evdience); + + + + + ///debug_hap_information(overlap_list, R_INF, g_read, hap, dumy); + + SetSnpMatrix(hap, hap->snp, overlap_list->length); + + + uint64_t pre_site = (uint64_t)-1; + uint64_t num_of_snps = 0; + long long pre_i = -1; + long long sub_length; + haplotype_evdience* sub_list; + + ///long long debug_total_length = 0; + + ////split reads + for (i = 0; i < hap->length; i++) + { + if(pre_site != hap->list[i].site) + { + if(i != 0) + { + sub_list = hap->list + pre_i; + sub_length = i - pre_i; + ///debug_total_length = debug_total_length + sub_length; + split_sub_list(hap, sub_list, sub_length, hap->snp, overlap_list, R_INF, g_read); + } + num_of_snps++; + pre_site = hap->list[i].site; + pre_i = i; + } + } + + if(pre_i != -1) + { + sub_list = hap->list + pre_i; + sub_length = i - pre_i; + ///debug_total_length = debug_total_length + sub_length; + split_sub_list(hap, sub_list, sub_length, hap->snp, overlap_list, R_INF, g_read); + } + + ///debug_snp_matrix(hap); + generate_haplotypes_DP(hap, overlap_list, R_INF, g_read->length, force_repeat); + ///generate_haplotypes_naive(hap, overlap_list, R_INF, g_read->length, force_repeat); + + lable_large_indels(overlap_list, R_INF, g_read->length, dumy); + + + ///debug_snp_matrix(hap); +} @@ -11545,7 +11763,8 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, - partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat); + ////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); // print_overlap("m64016_190918_162737/174131552/ccs", diff --git a/Correct.h b/Correct.h index 3ef973a..13b0630 100644 --- a/Correct.h +++ b/Correct.h @@ -1152,10 +1152,10 @@ inline void ResizeInitHaplotypeEvdience(haplotype_evdience_alloc* h) /****************************may have bugs********************************/ } -inline void RsetInitHaplotypeEvdienceFlag(haplotype_evdience_alloc* h, long long useful_length) +inline void RsetInitHaplotypeEvdienceFlag(haplotype_evdience_alloc* h, long long beg, long long useful_length) { /****************************may have bugs********************************/ - memset(h->flag, 0, useful_length * sizeof(uint8_t)); + memset(h->flag + beg, 0, useful_length * sizeof(uint8_t)); /****************************may have bugs********************************/ }