diff --git a/.vscode/settings.json b/.vscode/settings.json index 18a7993..a09f8ba 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -3,6 +3,9 @@ "files.associations": { "limits": "cpp", "random": "cpp", - "functional": "cpp" + "functional": "cpp", + "deque": "cpp", + "vector": "cpp", + "initializer_list": "cpp" } } \ No newline at end of file diff --git a/Assembly.cpp b/Assembly.cpp index e871bcb..e201ac6 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -8,6 +8,7 @@ #include "Hash_Table.h" #include "POA.h" #include "Correct.h" +#include "Output.h" Total_Count_Table TCB; Total_Pos_Table PCB; @@ -18,6 +19,9 @@ long long total_matched_overlap_0 = 0; long long total_matched_overlap_1 = 0; long long total_potiental_matched_overlap_0 = 0; long long total_potiental_matched_overlap_1 = 0; +long long total_num_read_base = 0; +long long total_num_correct_base = 0; + long long complete_threads = 0; @@ -682,6 +686,37 @@ void* Overlap_calculate(void* arg) destory_k_mer_pos_list_alloc(&array_list); } +inline void output_read_to_buffer(long long readID, All_reads* R_INF, char* corrected_read, long long correct_read_length, +Output_buffer_sub_block* current_sub_buffer) +{ + /** + if (seq->name.l + != Get_NAME_LENGTH(R_INF, read_number)) + { + fprintf(stderr, "name error\n"); + } + + + if(memcmp(seq->name.s, Get_NAME(R_INF, read_number), seq->name.l)) + { + fprintf(stderr, "name error\n"); + } + **/ + + + ///先清空 + current_sub_buffer->length = 0; + + ///一个是>一个是\n + add_base_to_sub_buffer(current_sub_buffer, '>'); + add_segment_to_sub_buffer(current_sub_buffer, Get_NAME((*R_INF), readID), Get_NAME_LENGTH((*R_INF), readID)); + add_base_to_sub_buffer(current_sub_buffer, '\n'); + add_segment_to_sub_buffer(current_sub_buffer, corrected_read, correct_read_length); + add_base_to_sub_buffer(current_sub_buffer, '\n'); + push_results_to_buffer(current_sub_buffer); + +} + void* Overlap_calculate_heap_merge(void* arg) { @@ -694,6 +729,8 @@ void* Overlap_calculate_heap_merge(void* arg) long long matched_overlap_1 = 0; long long potiental_matched_overlap_0 = 0; long long potiental_matched_overlap_1 = 0; + long long num_read_base = 0; + long long num_correct_base = 0; long long j; int thr_ID = *((int*)arg); @@ -739,10 +776,14 @@ void* Overlap_calculate_heap_merge(void* arg) Correct_dumy correct; init_Correct_dumy(&correct); + + Output_buffer_sub_block current_sub_buffer; + + init_buffer_sub_block(¤t_sub_buffer); + for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num) { - clear_Heap(&heap); clear_Candidates_list(&l); ///clear_Candidates_list(&debug_l); @@ -875,10 +916,17 @@ void* Overlap_calculate_heap_merge(void* arg) ///以x_pos_e，即结束位置为主元排序 calculate_overlap_region(&l, &overlap_list, i, g_read.length, &R_INF); + ///clear_Graph(&POA_Graph); - correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, + correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1); + num_read_base = num_read_base + g_read.length; + num_correct_base = num_correct_base + correct.corrected_base; + + output_read_to_buffer(i, &R_INF, correct.corrected_read, correct.corrected_read_length, ¤t_sub_buffer); + ///output_read_to_buffer(i, &R_INF, g_read.seq, g_read.length, ¤t_sub_buffer); + /** POA_i = 0; @@ -963,7 +1011,9 @@ void* Overlap_calculate_heap_merge(void* arg) // fprintf(stderr, "filtered_debug_overlap: %llu\n", filtered_debug_overlap); /************需要注释掉**********/ - + finish_output_buffer(); + + destory_buffer_sub_block(¤t_sub_buffer); @@ -992,6 +1042,9 @@ void* Overlap_calculate_heap_merge(void* arg) total_matched_overlap_1 += matched_overlap_1; total_potiental_matched_overlap_0 += potiental_matched_overlap_0; total_potiental_matched_overlap_1 += potiental_matched_overlap_1; + total_num_read_base += num_read_base; + total_num_correct_base += num_correct_base; + complete_threads++; if(complete_threads == thread_num) { @@ -999,6 +1052,8 @@ void* Overlap_calculate_heap_merge(void* arg) fprintf(stderr, "total_matched_overlap_1: %llu\n", total_matched_overlap_1); fprintf(stderr, "total_potiental_matched_overlap_0: %llu\n", total_potiental_matched_overlap_0); fprintf(stderr, "total_potiental_matched_overlap_1: %llu\n", total_potiental_matched_overlap_1); + fprintf(stderr, "total_num_read_base: %llu\n", total_num_read_base); + fprintf(stderr, "total_num_correct_base: %llu\n", total_num_correct_base); } pthread_mutex_unlock(&statistics); } @@ -1014,6 +1069,14 @@ void Overlap_calculate_multipe_thr() { double start_time = Get_T(); + init_output_buffer(thread_num); + + pthread_t outputResultSinkHandle; + + pthread_create(&outputResultSinkHandle, NULL, pop_buffer, NULL); + + + fprintf(stdout, "R_INF.total_reads: %llu\n", R_INF.total_reads); @@ -1039,8 +1102,12 @@ void Overlap_calculate_multipe_thr() for (i = 0; istart_i; i < overlap_list->length; i++) + { + ///只会发生在这个interval比list里所有元素都小的情况 + ///这种情况下一个interval需要从0开始 + if (window_end < overlap_list->list[i].x_pos_s) + { + dumy->start_i = 0; + dumy->length = 0; + dumy->lengthNT = 0; + return 0; + } + else ///只要window_end >= overlap_list->list[i].x_pos_s，就有可能重叠 + { + dumy->start_i = i; + break; + } + } + + ///只会发生在这个window比list里所有元素都大的情况 + ///这种情况下一个window也无需遍历了 + if (i >= overlap_list->length) + { + dumy->start_i = overlap_list->length; + dumy->length = 0; + dumy->lengthNT = 0; + return -2; + } + + dumy->length = 0; + dumy->lengthNT = 0; + + + + long long fake_length = 0; + + for (; i < overlap_list->length; i++) + { + ///是否重叠 + if((Len = OVERLAP(window_start, window_end, overlap_list->list[i].x_pos_s, overlap_list->list[i].x_pos_e)) > 0) + { + ///重叠数量 + fake_length++; + + ///重叠是否有效 + overlap_length = overlap_list->list[i].x_pos_e - overlap_list->list[i].x_pos_s + 1; + if (overlap_length * OVERLAP_THRESHOLD <= overlap_list->list[i].align_length) + { + dumy->overlapID[dumy->length] = i; + dumy->length++; + } + } + + if(overlap_list->list[i].x_pos_s > window_end) + { + break; + } + } + + ///fake_length是重叠的数量，而不是有效重叠的数量 + if (fake_length == 0) + { + return 0; + } + else + { + return 1; + } +} + ///Len = OVERLAP(window_start, window_end, overlap_list->list[i].x_pos_s, overlap_list->list[i].x_pos_e)) void print_string(char* s, int l) @@ -648,6 +725,7 @@ char* r_string) currentID = dumy->overlapID[reverse_i--]; x_start = MAX(window_start, overlap_list->list[currentID].x_pos_s); x_end = MIN(window_end, overlap_list->list[currentID].x_pos_e); + ///这个是和当前窗口重叠的长度 x_len = x_end - x_start + 1; threshold = x_len * THRESHOLD_RATE; @@ -766,7 +844,7 @@ void debug_stats(overlap_region_alloc* overlap_list, All_reads* R_INF, { Len_x = overlap_list->list[j].x_pos_e - overlap_list->list[j].x_pos_s + 1; - if (Len_x * 0.90 <= overlap_list->list[j].align_length) + if (Len_x * OVERLAP_THRESHOLD <= overlap_list->list[j].align_length) { if (overlap_list->list[j].y_pos_strand == 0) { @@ -998,6 +1076,8 @@ int verify_cigar(char* x, int x_len, char* y, int y_len, CIGAR* cigar, int error } + + /** ///记住等于0也要重算，虽然意义不那么大就是了 #define NEED_START_POS(x) (x.error<=0) @@ -1255,7 +1335,7 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN - /** + ///j负责遍历整个overlap list for (j = 0; j < overlap_list->length; j++) { @@ -1264,7 +1344,7 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN 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; - ///if (overlap_length * 0.90 <= overlap_list->list[j].align_length) + if (overlap_length * OVERLAP_THRESHOLD <= overlap_list->list[j].align_length) { for (i = 0; i < overlap_list->list[j].w_list_length; i++) { @@ -1306,7 +1386,7 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN } } } - **/ + @@ -1513,9 +1593,505 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN } +inline void add_base_to_correct_read(Correct_dumy* dumy, char base, int is_error) +{ + ///deletion就不要管 + if (base != 'D') + { + + if (dumy->corrected_read_length + 2 > dumy->corrected_read_size) + { + dumy->corrected_read_size = dumy->corrected_read_size * 2; + dumy->corrected_read = (char*)realloc(dumy->corrected_read, dumy->corrected_read_size); + } + + dumy->corrected_read[dumy->corrected_read_length] = base; + dumy->corrected_read_length++; + dumy->corrected_read[dumy->corrected_read_length] = '\0'; + } + + if (is_error) + { + dumy->corrected_base++; + } + +} + + +inline void add_segment_to_correct_read(Correct_dumy* dumy, char* segment, long long segment_length) +{ + + if (dumy->corrected_read_length + segment_length + 2 > dumy->corrected_read_size) + { + dumy->corrected_read_size = dumy->corrected_read_length + segment_length + 2; + dumy->corrected_read = (char*)realloc(dumy->corrected_read, dumy->corrected_read_size); + } + + memcpy(dumy->corrected_read + dumy->corrected_read_length, segment, segment_length); + dumy->corrected_read_length += segment_length; + dumy->corrected_read[dumy->corrected_read_length] = '\0'; +} + +///从backbone_start遍历到backbone_end节点，生成出来的seq要接着放到dumy->corrected_read中 +void get_seq_from_Graph(Graph* backbone, long long backbone_start, long long backbone_end, Correct_dumy* dumy) +{ + long long new_seq_length = 0; + long long currentNodeID; + long long i; + // 总共有以下几种情况: + // 1. match 2. mismatch (A, C, G, T, N) 3. deletion 4. insertion (A, C, G, T) + // 其实就是 1. 自己本身的weight 2. alignToNode的weight 3. insertion节点的weight + long long max_count; + int max_type; + long long max_node; + long long total_count; + long long nodeID; + char current_base; + + currentNodeID = backbone_start; + while (currentNodeID != backbone_end) + { + total_count = 0; + max_count = -1; + ///图上能够被遍历到的有两种节点 + ///1. backbone节点 2. insertion节点 + ///backbone节点才有match/mismatch/deletion + ///insertion这些都没有，就是无脑看出边 + + ///假如这是个backbone节点 + if (currentNodeID >= backbone_start && currentNodeID <= backbone_end) + { + ///match + total_count += backbone->g_nodes.list[currentNodeID].weight; + max_count = backbone->g_nodes.list[currentNodeID].weight; + max_type = 0; + max_node = currentNodeID; + ///mismatch和deletion (A, C, G, T, N, D, 除了自己的那个字符) + for (i = 0; i < backbone->g_nodes.list[currentNodeID].alignedTo_Nodes.length; i++) + { + nodeID = backbone->g_nodes.list[currentNodeID].alignedTo_Nodes.list[i].out_node; + total_count += backbone->g_nodes.list[nodeID].weight; + + if (backbone->g_nodes.list[nodeID].weight > max_count) + { + max_count = backbone->g_nodes.list[nodeID].weight; + max_type = 1; + max_node = nodeID; + } + } + ///insertion (A, C, G, T, N) + ///注意这个出边还得避开下一个backbone节点 + for (i = 0; i < backbone->g_nodes.list[currentNodeID].outcome_edges.length; i++) + { + if (backbone->g_nodes.list[currentNodeID].outcome_edges.list[i].weight == 2) + { + + nodeID = backbone->g_nodes.list[currentNodeID].outcome_edges.list[i].out_node; + total_count += backbone->g_nodes.list[nodeID].weight; + + if (backbone->g_nodes.list[nodeID].weight > max_count) + { + max_count = backbone->g_nodes.list[nodeID].weight; + max_type = 2; + max_node = nodeID; + } + + ///拿到的nodeID应该一定不是backbone上的，如果是就错了 + if (nodeID >= backbone_start && nodeID <= backbone_end) + { + fprintf(stderr, "error\n"); + } + + } + } + } + else ///如果是insertion节点，就无脑看出边 + { + ///insertion (A, C, G, T, N) + ///注意这个出边不用避开下一个backbone节点 + for (i = 0; i < backbone->g_nodes.list[currentNodeID].outcome_edges.length; i++) + { + if (backbone->g_nodes.list[currentNodeID].outcome_edges.list[i].weight == 2) + { + + nodeID = backbone->g_nodes.list[currentNodeID].outcome_edges.list[i].out_node; + total_count += backbone->g_nodes.list[nodeID].weight; + + if (backbone->g_nodes.list[nodeID].weight > max_count) + { + max_count = backbone->g_nodes.list[nodeID].weight; + max_type = 2; + max_node = nodeID; + } + + } + else ///如果边不是2就不对了 + { + fprintf(stderr, "error\n"); + } + + } + } + + + + + + + + if(max_count >= total_count*CORRECT_THRESHOLD) + { + current_base = backbone->g_nodes.list[max_node].base; + ///说明是insertion + if (max_type == 2) + { + currentNodeID = max_node; + } + else ///其他情况依然沿着backbone向前 + { + currentNodeID++; + } + } + else + { + ///假如这是个backbone节点, 不矫正 + if (currentNodeID >= backbone_start && currentNodeID <= backbone_end) + { + current_base = backbone->g_nodes.list[currentNodeID].base; + currentNodeID++; + } + else///如果在insertion节点上不达标很麻烦...,只能选最大的了 + { + current_base = backbone->g_nodes.list[max_node].base; + ///说明是insertion + if (max_type == 2) + { + currentNodeID = max_node; + } + else ///insertion节点不可能出现这种情况 + { + fprintf(stderr, "error\n"); + } + } + } + + + if (max_count <= 0) + { + fprintf(stderr, "error\n"); + } + + + add_base_to_correct_read(dumy, current_base, max_type); + + } + + + + + ///最后还要处理backbone_end这个节点 + total_count = 0; + max_count = -1; + ///这个节点肯定是backbone上的节点啊 + ///match + total_count += backbone->g_nodes.list[currentNodeID].weight; + max_count = backbone->g_nodes.list[currentNodeID].weight; + max_type = 0; + max_node = currentNodeID; + ///mismatch和deletion (A, C, G, T, N, D, 除了自己的那个字符) + for (i = 0; i < backbone->g_nodes.list[currentNodeID].alignedTo_Nodes.length; i++) + { + nodeID = backbone->g_nodes.list[currentNodeID].alignedTo_Nodes.list[i].out_node; + total_count += backbone->g_nodes.list[nodeID].weight; + + if (backbone->g_nodes.list[nodeID].weight > max_count) + { + max_count = backbone->g_nodes.list[nodeID].weight; + max_type = 1; + max_node = nodeID; + } + } + ///这个节点不应该有任何出边了 + if(backbone->g_nodes.list[currentNodeID].outcome_edges.length) + { + fprintf(stderr, "haha\n"); + } + + + if(max_count >= total_count*CORRECT_THRESHOLD) + { + current_base = backbone->g_nodes.list[max_node].base; + } + else + { + current_base = backbone->g_nodes.list[currentNodeID].base; + } + + + if (max_count <= 0) + { + fprintf(stderr, "error\n"); + } + + add_base_to_correct_read(dumy, current_base, max_type); + +} + + + +void window_consensus(char* r_string, long long window_start, long long window_end, +overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, Graph* g) +{ + clear_Graph(g); + + long long x_start; + long long x_length; + char* x_string; + char* y_string; + char* backbone; + long long backbone_length; + long long i; + long long y_start, y_length; + long long overlapID, windowID; + long long startNodeID, endNodeID, currentNodeID; + + ///这个和前面算alignment还不一样 + ///那个时候x_start和x_end是当前窗口内的overlap的起始和结束位置 + ///这个window就是要做consensus啊，所以起始和结束就是window本身，固定的 + backbone = r_string + window_start; + backbone_length = window_end - window_start + 1; + + addUnmatchedSeqToGraph(g, backbone, backbone_length, &startNodeID, &endNodeID); + + + long long correct_x_pos_s; + ///与当前window重叠的所有overlap + for (i = 0; i < dumy->length; i++) + { + ///这个是那个overlap的ID，而不是overlap里对应窗口的ID + overlapID = dumy->overlapID[i]; + + correct_x_pos_s = (overlap_list->list[overlapID].x_pos_s / WINDOW) * WINDOW; + windowID = (window_start - correct_x_pos_s) / WINDOW; + + ///如果这个window不匹配，跳过 + if (overlap_list->list[overlapID].w_list[windowID].y_end == -1) + { + continue; + } + + + 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; + ///这个是比对上的起始base在backbone上对应的位置，也就是节点ID + currentNodeID = x_start - window_start; + ///这个是要用的cigar: overlap_list->list[overlapID].w_list[windowID].cigar; + + addmatchedSeqToGraph(g, currentNodeID, x_string, x_length, + y_string, y_length, &(overlap_list->list[overlapID].w_list[windowID].cigar), startNodeID, endNodeID); + } + + + get_seq_from_Graph(g, startNodeID, endNodeID, dumy); + + + /** + if (startNodeID != 0 || endNodeID != backbone_length - 1) + { + fprintf(stderr, "error\n"); + } + + + for (i = startNodeID; i <= backbone_length; i++) + { + if(g->g_nodes.list[i].alignedTo_Nodes.length > 5) + { + fprintf(stderr, "error 1\n"); + } + + if(g->g_nodes.list[i].outcome_edges.length > 4) + { + fprintf(stderr, "error 2\n"); + } + } + + + for (i = 0; i < dumy->length; i++) + { + ///这个是那个overlap的ID，而不是overlap里对应窗口的ID + overlapID = dumy->overlapID[i]; + + correct_x_pos_s = (overlap_list->list[overlapID].x_pos_s / WINDOW) * WINDOW; + windowID = (window_start - correct_x_pos_s) / WINDOW; + + ///如果这个window不匹配，跳过 + if (overlap_list->list[overlapID].w_list[windowID].y_end == -1) + { + continue; + } + + + 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; + ///这个是比对上的起始base在backbone上对应的位置，也就是节点ID + currentNodeID = x_start - window_start; + ///这个是要用的cigar: overlap_list->list[overlapID].w_list[windowID].cigar; + + Graph_debug(g, currentNodeID, x_string, x_length, + y_string, y_length, &(overlap_list->list[overlapID].w_list[windowID].cigar), startNodeID, endNodeID); + } + + for (i = 0; i < g->g_nodes.length; i++) + { + if (i >= startNodeID && i <= endNodeID) + { + if (g->g_nodes.list[i].weight != 1) + { + fprintf(stderr, "error 1\n"); + } + } + else + { + if (g->g_nodes.list[i].weight != 0) + { + fprintf(stderr, "error 2\n"); + } + + if (g->g_nodes.list[i].alignedTo_Nodes.length != 0) + { + fprintf(stderr, "error 3\n"); + } + + ///节点入边不为0，说明这个不是alignTO节点，而是insert节点 + if (g->g_nodes.list[i].income_edges.length != 0 && g->g_nodes.list[i].outcome_edges.length != 0) + { + fprintf(stderr, "error 4\n"); + } + } + } + **/ + + + + +} + + +void generate_consensus(overlap_region_alloc* overlap_list, All_reads* R_INF, + UC_Read* g_read, Correct_dumy* dumy, Graph* g) +{ + + long long window_num = (g_read->length + WINDOW - 1) / WINDOW; + long long i, j, overlap_length; + long long window_start, window_end; + + long long num_availiable_win = 0; + + + window_start = 0; + window_end = WINDOW - 1; + if (window_end >= g_read->length) + { + window_end = g_read->length - 1; + } + int flag; + for (i = 0; i < window_num; i++) + { + dumy->length = 0; + dumy->lengthNT = 0; + ///flag返回的是重叠数量 + ///dumy->length返回的是有效完全重叠的数量 + ///dumy->lengthNT返回的是有效不完全重叠的数量 + flag = get_available_interval(window_start, window_end, overlap_list, dumy); + switch (flag) + { + case 1: ///找到匹配 + break; + case 0: ///没找到匹配 + break; + case -2: ///下一个window也不会存在匹配, 直接跳出 + i = window_num; + break; + } + + + ///这个是available overlap里所有window的数量... + ///num_availiable_win = num_availiable_win + dumy->length + dumy->lengthNT; + num_availiable_win = num_availiable_win + dumy->length; + + ///重叠窗口数，也就是coverage大小 + if(dumy->length >= MIN_COVERAGE_THRESHOLD) + { + window_consensus(g_read->seq, window_start, window_end, overlap_list, dumy, R_INF, g); + } + else + { + add_segment_to_correct_read(dumy, g_read->seq + window_start, window_end - window_start + 1); + } + + window_start = window_start + WINDOW; + window_end = window_end + WINDOW; + if (window_end >= g_read->length) + { + window_end = g_read->length - 1; + } + + } + + if (window_start < g_read->length) + { + add_segment_to_correct_read(dumy, g_read->seq + window_start, g_read->length - window_start); + } + + + + + + + /***********************要注释掉*************************/ + // long long debug_num_availiable_win = 0; + // for (j = 0; j < overlap_list->length; j++) + // { + // overlap_length = overlap_list->list[j].x_pos_e - overlap_list->list[j].x_pos_s + 1; + + // if (overlap_length * OVERLAP_THRESHOLD <= overlap_list->list[j].align_length) + // { + // debug_num_availiable_win = debug_num_availiable_win + overlap_list->list[j].w_list_length; + // } + // } + + // if (debug_num_availiable_win != num_availiable_win) + // { + // fprintf(stderr, "error, debug_num_availiable_win: %d, num_availiable_win: %d\n", + // debug_num_availiable_win, num_availiable_win); + // } + /***********************要注释掉*************************/ +} void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, - UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, + UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, Graph* g, long long* matched_overlap_0, long long* matched_overlap_1, long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1) { @@ -1531,6 +2107,11 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, window_start = 0; window_end = WINDOW - 1; + if (window_end >= g_read->length) + { + window_end = g_read->length - 1; + } + int flag; for (i = 0; i < window_num; i++) @@ -1574,7 +2155,45 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read); debug_stats(overlap_list, R_INF, g_read, dumy, overlap_read, matched_overlap_0, matched_overlap_1); + + + generate_consensus(overlap_list, R_INF, g_read, dumy, g); + + + ///fprintf(stderr, "length: %lld, corrected_base: %lld\n", g_read->length, dumy->corrected_base); + /** + EdlibAlignResult result = edlibAlign(g_read->seq, g_read->length, dumy->corrected_read, dumy->corrected_read_length, + edlibNewAlignConfig(-1, EDLIB_MODE_NW, EDLIB_TASK_PATH, NULL, 0)); + + if (result.status == EDLIB_STATUS_OK) { + if (dumy->corrected_base == 0 && result.editDistance!= 0) + { + fprintf(stderr, "error 0\n"); + } + + + if (dumy->corrected_base < result.editDistance) + { + fprintf(stderr, "error 1\n"); + } + + // fprintf(stderr, "****\n distance: %d, alignmentLength: %d, startLocations: %d, endLocations: %d, corrected_base: %d\n", + // result.editDistance, result.alignmentLength, result.startLocations[0], result.endLocations[0], dumy->corrected_base); + + + // char* cigar = edlibAlignmentToCigar(result.alignment, result.alignmentLength, EDLIB_CIGAR_STANDARD); + // fprintf(stderr,"%s\n", cigar); + // free(cigar); + + } + else + { + fprintf(stderr, "error\n"); + } + + edlibFreeAlignResult(result); + **/ /************需要注释掉********* */ // pthread_mutex_lock(&debug_statistics); @@ -1599,11 +2218,19 @@ void init_Correct_dumy(Correct_dumy* list) { list->Peq_SSE[i] = _mm_setzero_si128(); } + + + list->corrected_read_size = 1000; + list->corrected_read_length = 0; + list->corrected_read = (char*)malloc(sizeof(char)*list->corrected_read_size); + list->corrected_base = 0; + } void destory_Correct_dumy(Correct_dumy* list) { free(list->overlapID); + free(list->corrected_read); } void clear_Correct_dumy(Correct_dumy* list, overlap_region_alloc* overlap_list) @@ -1617,6 +2244,9 @@ void clear_Correct_dumy(Correct_dumy* list, overlap_region_alloc* overlap_list) list->size = overlap_list->length; list->overlapID = (uint64_t*)realloc(list->overlapID, list->size*sizeof(uint64_t)); } + + list->corrected_read_length = 0; + list->corrected_base = 0; } diff --git a/Correct.h b/Correct.h index 2ae3438..a83ccf6 100644 --- a/Correct.h +++ b/Correct.h @@ -3,7 +3,10 @@ #include #include "Hash_Table.h" #include "Levenshtein_distance.h" +#include "POA.h" +#define CORRECT_THRESHOLD 0.6 +#define MIN_COVERAGE_THRESHOLD 4 @@ -14,6 +17,11 @@ typedef struct { + char* corrected_read; + long long corrected_read_length; + long long corrected_read_size; + long long corrected_base; + uint64_t* overlapID; uint64_t length; uint64_t lengthNT; @@ -29,7 +37,7 @@ typedef struct void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, - UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, + UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, Graph* g, long long* matched_overlap_0, long long* matched_overlap_1, long long* potiental_matched_overlap_0, long long* potiental_matched_overlap_1); void init_Correct_dumy(Correct_dumy* list); diff --git a/Hash_Table.h b/Hash_Table.h index f6a4355..78219e1 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -17,10 +17,12 @@ typedef khash_t(POS64) Pos_Table; ///#define THRESHOLD 14 #define WINDOW 375 #define THRESHOLD 15 +///#define THRESHOLD_RATE 0.04 #define THRESHOLD_RATE 0.04 #define GROUP_SIZE 4 ///最长是10M10D10M10D10M这种 #define CIGAR_MAX_LENGTH THRESHOLD*2+2 +#define OVERLAP_THRESHOLD 0.9 typedef struct { diff --git a/Makefile b/Makefile index c4fa7f3..d2e34d8 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ CC=g++ CFLAGS = -w -c -msse4.2 -mpopcnt -fomit-frame-pointer -Winline -O3 -lz LDFLAGS = -lm -lz -lpthread -O3 -mpopcnt -msse4.2 -lz -w -SOURCES = main.cpp CommandLines.cpp Process_Read.cpp Assembly.cpp kmer.cpp Hash_Table.cpp POA.cpp Correct.cpp Levenshtein_distance.cpp edlib.cpp +SOURCES = main.cpp Output.cpp CommandLines.cpp Process_Read.cpp Assembly.cpp kmer.cpp Hash_Table.cpp POA.cpp Correct.cpp Levenshtein_distance.cpp edlib.cpp OBJECTS = $(SOURCES:.c=.o) EXECUTABLE = ccs_assembly diff --git a/Output.cpp b/Output.cpp new file mode 100644 index 0000000..b294583 --- /dev/null +++ b/Output.cpp @@ -0,0 +1,234 @@ +#include "Output.h" +#include "CommandLines.h" +#include + +pthread_mutex_t o_queueMutex; +pthread_cond_t o_flushCond; +pthread_cond_t o_stallCond; +pthread_mutex_t o_doneMutex; + + +Output_buffer buffer_out; +Output_buffer_sub_block tmp_buffer_sub_block; + +void init_buffer_sub_block(Output_buffer_sub_block* sub_block) +{ + sub_block->length = 0; + sub_block->size = SUB_BLOCK_INIT_SIZE; + sub_block->buffer = (char*)malloc(sub_block->size); +} + +void destory_buffer_sub_block(Output_buffer_sub_block* sub_block) +{ + free(sub_block->buffer); +} + +void destory_output_buffer() +{ + for (int i = 0; i < buffer_out.sub_block_size; i++) + { + destory_buffer_sub_block(&(buffer_out.sub_buffer[i])); + } + + free(buffer_out.sub_buffer); +} + +void init_output_buffer(int thread_number) +{ + buffer_out.sub_block_size = OUTPUT_BUFFER_SIZE * thread_number; + buffer_out.sub_block_number = 0; + + buffer_out.sub_buffer = (Output_buffer_sub_block*)malloc(sizeof(Output_buffer_sub_block)*buffer_out.sub_block_size); + + + + for (int i = 0; i < buffer_out.sub_block_size; i++) + { + init_buffer_sub_block(&(buffer_out.sub_buffer[i])); + } + + buffer_out.all_buffer_end = 0; +} + +inline int if_empty_buffer() +{ + + if (buffer_out.sub_block_number == 0) + { + return 1; + } + else + { + return 0; + } +} + + +inline int if_full_buffer() +{ + + if (buffer_out.sub_block_number >= buffer_out.sub_block_size) + { + return 1; + } + else + { + return 0; + } +} + +inline void pop_single_buffer(Output_buffer_sub_block* curr_sub_block) +{ + buffer_out.sub_block_number--; + char *k; + k = buffer_out.sub_buffer[buffer_out.sub_block_number].buffer; + buffer_out.sub_buffer[buffer_out.sub_block_number].buffer = curr_sub_block->buffer; + curr_sub_block->buffer = k; + + + long long tmp_size; + tmp_size = curr_sub_block->size; + curr_sub_block->size = buffer_out.sub_buffer[buffer_out.sub_block_number].size; + buffer_out.sub_buffer[buffer_out.sub_block_number].size = tmp_size; + + + curr_sub_block->length = buffer_out.sub_buffer[buffer_out.sub_block_number].length; + buffer_out.sub_buffer[buffer_out.sub_block_number].length = 0; +} + + + +void add_segment_to_sub_buffer(Output_buffer_sub_block* current_sub_buffer, char* seg, long long segLen) +{ + if(current_sub_buffer->length + segLen + 2 > current_sub_buffer->size) + { + current_sub_buffer->size = current_sub_buffer->length + segLen + 2; + current_sub_buffer->buffer = (char*)realloc(current_sub_buffer->buffer, current_sub_buffer->size); + } + + memcpy(current_sub_buffer->buffer + current_sub_buffer->length, seg, segLen); + current_sub_buffer->length += segLen; + current_sub_buffer->buffer[current_sub_buffer->length] = '\0'; +} + + +void add_base_to_sub_buffer(Output_buffer_sub_block* current_sub_buffer, char base) +{ + if(current_sub_buffer->length + 2 > current_sub_buffer->size) + { + current_sub_buffer->size = current_sub_buffer->length + 2; + current_sub_buffer->buffer = (char*)realloc(current_sub_buffer->buffer, current_sub_buffer->size); + } + + current_sub_buffer->buffer[current_sub_buffer->length] = base; + current_sub_buffer->length++; + current_sub_buffer->buffer[current_sub_buffer->length] = '\0'; + +} + + +inline void push_single_buffer(Output_buffer_sub_block* curr_sub_block) +{ + + char *k; + k = buffer_out.sub_buffer[buffer_out.sub_block_number].buffer; + buffer_out.sub_buffer[buffer_out.sub_block_number].buffer = curr_sub_block->buffer; + curr_sub_block->buffer = k; + + long long tmp_size; + tmp_size = curr_sub_block->size; + curr_sub_block->size = buffer_out.sub_buffer[buffer_out.sub_block_number].size; + buffer_out.sub_buffer[buffer_out.sub_block_number].size = tmp_size; + + buffer_out.sub_buffer[buffer_out.sub_block_number].length = curr_sub_block->length; + curr_sub_block->length = 0; + + buffer_out.sub_block_number++; +} + + +void* pop_buffer(void*) +{ + + FILE* output_file = fopen(output_file_name, "w"); + + init_buffer_sub_block(&tmp_buffer_sub_block); + + + + while (buffer_out.all_buffer_end < thread_num) + { + + pthread_mutex_lock(&o_queueMutex); + + while (if_empty_buffer() && (buffer_out.all_buffer_end < thread_num)) + { + pthread_cond_signal(&o_stallCond); + pthread_cond_wait(&o_flushCond, &o_queueMutex); + } + + if (!if_empty_buffer()) + { + pop_single_buffer(&tmp_buffer_sub_block); + } + + pthread_cond_signal(&o_stallCond); + pthread_mutex_unlock(&o_queueMutex); + + + + if (tmp_buffer_sub_block.length != 0) + { + fprintf(output_file, "%s", tmp_buffer_sub_block.buffer); + } + + } + + + while (buffer_out.sub_block_number>0) + { + buffer_out.sub_block_number--; + fprintf(output_file, "%s", buffer_out.sub_buffer[buffer_out.sub_block_number].buffer); + } + + + destory_buffer_sub_block(&tmp_buffer_sub_block); + + fclose(output_file); + +} + + +void push_results_to_buffer(Output_buffer_sub_block* sub_block) +{ + + pthread_mutex_lock(&o_queueMutex); + + while (if_full_buffer()) + { + pthread_cond_signal(&o_flushCond); + pthread_cond_wait(&o_stallCond, &o_queueMutex); + } + + push_single_buffer(sub_block); + pthread_cond_signal(&o_flushCond); + pthread_mutex_unlock(&o_queueMutex); + +} + + + +void finish_output_buffer() +{ + + pthread_mutex_lock(&o_doneMutex); + buffer_out.all_buffer_end++; + + + if (buffer_out.all_buffer_end == thread_num) + { + pthread_cond_signal(&o_flushCond); + } + + pthread_mutex_unlock(&o_doneMutex); +} \ No newline at end of file diff --git a/Output.h b/Output.h new file mode 100644 index 0000000..9f95698 --- /dev/null +++ b/Output.h @@ -0,0 +1,37 @@ +#ifndef __OUTPUT__ +#define __OUTPUT__ + +#include +#include +#include +#include + +typedef struct +{ + char* buffer; + long long size; + long long length; +} Output_buffer_sub_block; + +typedef struct Output_buffer +{ + Output_buffer_sub_block* sub_buffer; + long long sub_block_size; + long long sub_block_number; + int all_buffer_end; +} Output_buffer; + +#define OUTPUT_BUFFER_SIZE 100 +#define SUB_BLOCK_INIT_SIZE 10000 + +void init_buffer_sub_block(Output_buffer_sub_block* sub_block); +void* pop_buffer(void*); +void add_segment_to_sub_buffer(Output_buffer_sub_block* current_sub_buffer, char* seg, long long segLen); +void add_base_to_sub_buffer(Output_buffer_sub_block* current_sub_buffer, char base); +void push_results_to_buffer(Output_buffer_sub_block* sub_block); +void finish_output_buffer(); +void destory_buffer_sub_block(Output_buffer_sub_block* sub_block); +void destory_output_buffer(); +void init_output_buffer(int thread_number); + +#endif diff --git a/POA.cpp b/POA.cpp index a65f8bf..c649292 100644 --- a/POA.cpp +++ b/POA.cpp @@ -1,5 +1,6 @@ #include "POA.h" #include +#include "Correct.h" #define INIT_EDGE_SIZE 50 #define INCREASE_EDGE_SIZE 5 #define INIT_NODE_SIZE 16000 @@ -142,6 +143,10 @@ void init_Node_alloc(Node_alloc* list) list->list[i].outcome_edges.list=NULL; list->list[i].alignedTo_Nodes.list=NULL; } + + list->total_start.income_edges.list = NULL; + list->total_start.outcome_edges.list = NULL; + list->total_start.alignedTo_Nodes.list = NULL; } @@ -154,6 +159,10 @@ void destory_Node_alloc(Node_alloc* list) destory_Edge_alloc(&list->list[i].outcome_edges); destory_Edge_alloc(&list->list[i].alignedTo_Nodes); } + destory_Edge_alloc(&list->total_start.income_edges); + destory_Edge_alloc(&list->total_start.outcome_edges); + destory_Edge_alloc(&list->total_start.alignedTo_Nodes); + free(list->list); free(list->sort.list); free(list->sort.visit); @@ -171,6 +180,9 @@ void clear_Node_alloc(Node_alloc* list) clear_Edge_alloc(&list->list[i].outcome_edges); clear_Edge_alloc(&list->list[i].alignedTo_Nodes); } + clear_Edge_alloc(&list->total_start.income_edges); + clear_Edge_alloc(&list->total_start.outcome_edges); + clear_Edge_alloc(&list->total_start.alignedTo_Nodes); list->length = 0; } @@ -198,6 +210,7 @@ uint64_t append_Node_alloc(Node_alloc* list, char base) list->list[list->length].ID = list->length; list->list[list->length].base = base; + list->list[list->length].weight = 1; init_Edge_alloc(&list->list[list->length].income_edges); init_Edge_alloc(&list->list[list->length].outcome_edges); init_Edge_alloc(&list->list[list->length].alignedTo_Nodes); @@ -292,7 +305,8 @@ void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length, } if (lastID != -1) { - add_Edge_Graph(g, lastID, nodeID, 1); + ///0是match边 + add_Edge_Graph(g, lastID, nodeID, 0); } lastID = nodeID; @@ -303,6 +317,441 @@ void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length, } +inline long long get_alignToNode(Graph* backbone, long long currentNodeID, char base) +{ + if(backbone->g_nodes.list[currentNodeID].alignedTo_Nodes.length == 0) + { + return -1; + } + + long long i = 0; + long long nodeID; + for (i = 0; i < backbone->g_nodes.list[currentNodeID].alignedTo_Nodes.length; i++) + { + nodeID = backbone->g_nodes.list[currentNodeID].alignedTo_Nodes.list[i].out_node; + if(backbone->g_nodes.list[nodeID].base == base) + { + return nodeID; + } + } + + return -1; + +} + + +inline void add_mismatch_to_backbone(Graph* backbone, long long* alignNodeID, char* mis_base, long long mis_base_length) +{ + long long i; + long long mismatch_nodeID; + char base; + + for (i = 0; i < mis_base_length; i++, (*alignNodeID)++) + { + base = mis_base[i]; + mismatch_nodeID = get_alignToNode(backbone, *alignNodeID, base); + + ///如果已经存在一个误配节点，那么给误配节点的权重+1 + if (mismatch_nodeID != -1) + { + backbone->g_nodes.list[mismatch_nodeID].weight++; + } + else + { + mismatch_nodeID = add_Node_Graph(backbone, base); + ///1代表是mismatch边 + append_Edge_alloc(&backbone->g_nodes.list[*alignNodeID].alignedTo_Nodes, + *alignNodeID, mismatch_nodeID, 1); + } + } + +} + +inline void add_deletion_to_backbone(Graph* backbone, long long* alignNodeID, long long deletion_length) +{ + long long i; + long long mismatch_nodeID; + char base; + + for (i = 0; i < deletion_length; i++, (*alignNodeID)++) + { + base = 'D'; + mismatch_nodeID = get_alignToNode(backbone, *alignNodeID, base); + + ///如果已经存在一个误配节点，那么给误配节点的权重+1 + if (mismatch_nodeID != -1) + { + backbone->g_nodes.list[mismatch_nodeID].weight++; + } + else + { + mismatch_nodeID = add_Node_Graph(backbone, base); + ///3代表是deletion边 + append_Edge_alloc(&backbone->g_nodes.list[*alignNodeID].alignedTo_Nodes, + *alignNodeID, mismatch_nodeID, 3); + } + } + +} + + + +///注意这里返回的有可能是新加的节点，也有可能返回的是backbone上的节点 +inline long long get_insertion_Node(Graph* backbone, long long currentNodeID, char base) +{ + ///看出边数量 + if(backbone->g_nodes.list[currentNodeID].outcome_edges.length == 0) + { + return -1; + } + + long long i = 0; + long long nodeID; + int type; + ///遍历所有出边 + for (i = 0; i < backbone->g_nodes.list[currentNodeID].outcome_edges.length; i++) + { + ///出边类型 + type = backbone->g_nodes.list[currentNodeID].outcome_edges.list[i].weight; + + ///为2的时候才是deletion边 + ///这个边有可能是match边，也就是type = 0 + ///这个似乎不需要...，加上反而坏事 + ///也不一定 + if (type == 2) + { + nodeID = backbone->g_nodes.list[currentNodeID].outcome_edges.list[i].out_node; + if(backbone->g_nodes.list[nodeID].base == base) + { + return nodeID; + } + } + } + + return -1; +} + + + +///注意这里返回的有可能是新加的节点，也有可能返回的是backbone上的节点 +inline void link_insertion_Node(Graph* backbone, long long currentNodeID, long long backboneNodeID) +{ + ///如果出边数量为0，那这就是个新节点 + if(backbone->g_nodes.list[currentNodeID].outcome_edges.length == 0) + { + + ///2代表是insertion边 + add_Edge_Graph(backbone, currentNodeID, backboneNodeID, 2); + } + else ////如果不为0，backboneNodeID应该一定在出边中 + { + long long i = 0; + long long nodeID; + int type; + ///遍历所有出边 + for (i = 0; i < backbone->g_nodes.list[currentNodeID].outcome_edges.length; i++) + { + ///出边类型 + type = backbone->g_nodes.list[currentNodeID].outcome_edges.list[i].weight; + + nodeID = backbone->g_nodes.list[currentNodeID].outcome_edges.list[i].out_node; + if(backboneNodeID == nodeID) + { + break; + } + + } + + ///如果这个节点没有被连到backboneNodeID上，就要处理 + if (i >= backbone->g_nodes.list[currentNodeID].outcome_edges.length) + { + ///2代表是insertion边 + add_Edge_Graph(backbone, currentNodeID, backboneNodeID, 2); + } + + } + +} + + +inline void add_insertion_to_backbone(Graph* backbone, long long alignNodeID, char* insertion_base, long long insertion_length, +long long backbone_start, long long backbone_end) +{ + long long i; + long long insertion_nodeID; + long long backboneID = alignNodeID + 1; + char base; + + for (i = 0; i < insertion_length; i++) + { + base = insertion_base[i]; + insertion_nodeID = get_insertion_Node(backbone, alignNodeID, base); + + ///如果已经存在一个insertion节点，那么给insertion节点的权重+1 + ///注意这里返回的有可能是新加的节点，也有可能返回的是backbone上的节点 + ///不可能，这里要是返回了backbone上的节点就错了，最后要验证下 + if (insertion_nodeID != -1) + { + /** + if (insertion_nodeID >= backbone_start && insertion_nodeID <= backbone_end) + { + fprintf(stderr, "error\n"); + } + **/ + + backbone->g_nodes.list[insertion_nodeID].weight++; + } + else + { + insertion_nodeID = add_Node_Graph(backbone, base); + ///2代表是insertion边 + add_Edge_Graph(backbone, alignNodeID, insertion_nodeID, 2); + } + + alignNodeID = insertion_nodeID; + } + + ///最后要把节点接回到backbone上去 + link_insertion_Node(backbone, alignNodeID, backboneID); +} + + + +void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_string, long long x_length, + char* y_string, long long y_length, CIGAR* cigar, long long backbone_start, long long backbone_end) +{ + int x_i, y_i, cigar_i; + x_i = 0; + y_i = 0; + cigar_i = 0; + int operation; + int operationLen; + int i; + + + + + ///0 is match, 1 is mismatch, 2 is up, 3 is left + ///2是x缺字符（y多字符），而3是y缺字符（x多字符） + ///while (x_i < x_len && y_i < y_len && cigar_i < cigar->length) + while (cigar_i < cigar->length) + { + operation = cigar->C_C[cigar_i]; + operationLen = cigar->C_L[cigar_i]; + + ///这种情况代表匹配 + if (operation == 0) + { + + for (i = 0; i < operationLen; i++) + { + backbone->g_nodes.list[currentNodeID].weight++; + + x_i++; + y_i++; + currentNodeID++; + } + } + else if (operation == 1) + { + add_mismatch_to_backbone(backbone, ¤tNodeID, y_string + y_i, operationLen); + x_i = x_i + operationLen; + y_i = y_i + operationLen; + } + else if (operation == 2) + { + ///记住要传currentNodeID - 1而不是currentNodeID + ///cigar的起始和结尾不可能是2，所以这里-1没问题 + add_insertion_to_backbone(backbone, currentNodeID - 1, y_string + y_i, operationLen, backbone_start, backbone_end); + y_i += operationLen; + } + else if (operation == 3) + { + ///3是y缺字符（x多字符），也就是backbone多字符 + ///这个相当于在backbone对应字符处变成了‘——’ + ///因此可以用mismatch类似的方法处理 + add_deletion_to_backbone(backbone, ¤tNodeID, operationLen); + x_i += operationLen; + } + + cigar_i++; + } + + + /** + ///cigar的起始和结尾不可能是2 + if (cigar->C_C[0] == 2 || cigar->C_C[cigar->length - 1] == 2) + { + fprintf(stderr, "error\n"); + } + + + if (x_i != x_length) + { + fprintf(stderr, "x_i: %d, x_length: %d\n", x_i, x_length); + } + + if (y_i != y_length) + { + fprintf(stderr, "y_i: %d, y_length: %d\n", y_i, y_length); + } + **/ + +} + +void Graph_debug(Graph* backbone, long long currentNodeID, char* x_string, long long x_length, + char* y_string, long long y_length, CIGAR* cigar, long long backbone_start, long long backbone_end) +{ + int x_i, y_i, cigar_i; + x_i = 0; + y_i = 0; + cigar_i = 0; + int operation; + int operationLen; + int i; + + + + + ///0 is match, 1 is mismatch, 2 is up, 3 is left + ///2是x缺字符（y多字符），而3是y缺字符（x多字符） + ///while (x_i < x_len && y_i < y_len && cigar_i < cigar->length) + while (cigar_i < cigar->length) + { + operation = cigar->C_C[cigar_i]; + operationLen = cigar->C_L[cigar_i]; + + ///这种情况代表匹配 + if (operation == 0) + { + + for (i = 0; i < operationLen; i++) + { + if (backbone->g_nodes.list[currentNodeID].base != y_string[y_i]) + { + fprintf(stderr, "error match\n"); + } + + backbone->g_nodes.list[currentNodeID].weight--; + + x_i++; + y_i++; + currentNodeID++; + } + } + else if (operation == 1) + { + for (i = 0; i < operationLen; i++) + { + if (backbone->g_nodes.list[currentNodeID].base == y_string[y_i]) + { + fprintf(stderr, "error mismatch 1\n"); + } + + long long mismatchID = get_alignToNode(backbone, currentNodeID, y_string[y_i]); + + + + if(mismatchID == -1) + { + fprintf(stderr, "error mismatch 2\n"); + } + else + { + backbone->g_nodes.list[mismatchID].weight--; + } + + + x_i++; + y_i++; + currentNodeID++; + } + } + else if (operation == 2) + { + long long nodeID = currentNodeID - 1; + long long mismatchID; + + for (i = 0; i < operationLen; i++) + { + mismatchID = get_insertion_Node(backbone, nodeID, y_string[y_i]); + + if (mismatchID == -1) + { + fprintf(stderr, "error insertion 1, i: %d\n", i); + } + else + { + backbone->g_nodes.list[mismatchID].weight--; + } + + nodeID = mismatchID; + + y_i++; + } + ///注意这里是x_string[x_i]而不是x_string[currentNodeID] + mismatchID = get_insertion_Node(backbone, nodeID, x_string[x_i]); + if (mismatchID == -1) + { + fprintf(stderr, "error insertion 2, i: %d, x_i: %d\n", i, x_i); + } + + + if (mismatchID != currentNodeID) + { + fprintf(stderr, "error insertion 3, i: mismatchID: %d, currentNodeID: %d\n", mismatchID, currentNodeID); + } + + + + + } + else if (operation == 3) + { + for (i = 0; i < operationLen; i++) + { + + long long mismatchID = get_alignToNode(backbone, currentNodeID, 'D'); + + if(mismatchID == -1) + { + fprintf(stderr, "error deletion 2\n"); + } + else + { + backbone->g_nodes.list[mismatchID].weight--; + } + + x_i++; + currentNodeID++; + } + } + + cigar_i++; + } + + + if (cigar->C_C[0] == 2 || cigar->C_C[cigar->length - 1] == 2) + { + fprintf(stderr, "error\n"); + } + + + if (x_i != x_length) + { + fprintf(stderr, "x_i: %d, x_length: %d\n", x_i, x_length); + } + + if (y_i != y_length) + { + fprintf(stderr, "y_i: %d, y_length: %d\n", y_i, y_length); + } + +} + + + + void Perform_POA(Graph* g, overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read) { diff --git a/POA.h b/POA.h index 9571cf3..a392341 100644 --- a/POA.h +++ b/POA.h @@ -27,6 +27,7 @@ typedef struct { uint64_t in_node; uint64_t out_node; + ///0是match，1是mismatch，2是x缺字符（y多字符），而3是y缺字符（x多字符） uint64_t weight; } Edge; @@ -40,6 +41,7 @@ typedef struct typedef struct { uint64_t ID; + uint64_t weight; char base; Edge_alloc income_edges; Edge_alloc outcome_edges; @@ -60,6 +62,7 @@ typedef struct typedef struct { + Node total_start; Node* list; topo_Sorting_buffer sort; uint64_t size; @@ -96,9 +99,14 @@ uint64_t* get_Topo_Sort_Order(Node_alloc* list, int need_sort); void init_Graph(Graph* g); void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length, long long* startID, long long* endID); +void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_string, long long x_length, + char* y_string, long long y_length, CIGAR* cigar, long long backbone_start, long long backbone_end); void destory_Graph(Graph* g); void clear_Graph(Graph* g); void Perform_POA(Graph* g, overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read); +void Graph_debug(Graph* backbone, long long currentNodeID, char* x_string, long long x_length, + char* y_string, long long y_length, CIGAR* cigar, long long backbone_start, long long backbone_end); + #endif \ No newline at end of file