diff --git a/.vscode/settings.json b/.vscode/settings.json index a09f8ba..edb48a6 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -6,6 +6,7 @@ "functional": "cpp", "deque": "cpp", "vector": "cpp", - "initializer_list": "cpp" + "initializer_list": "cpp", + "*.tcc": "cpp" } } \ No newline at end of file diff --git a/Assembly.cpp b/Assembly.cpp index e201ac6..0a2f486 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -783,6 +783,8 @@ void* Overlap_calculate_heap_merge(void* arg) for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num) { + + clear_Heap(&heap); clear_Candidates_list(&l); diff --git a/CommandLines.cpp b/CommandLines.cpp index 0b2e7a2..a90da2d 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -10,7 +10,8 @@ char* output_file_name = NULL; int thread_num = 1; int k_mer_length = 40; //int k_mer_min_freq = 9; -int k_mer_min_freq = 3; +//int k_mer_min_freq = 3; +int k_mer_min_freq = 2; int k_mer_max_freq = 66; int load_index_from_disk = 0; int write_index_to_disk = 0; diff --git a/Correct.cpp b/Correct.cpp index 55ef66e..daeb78b 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -1593,6 +1593,22 @@ inline void recalcate_window(overlap_region_alloc* overlap_list, All_reads* R_IN } + +inline void add_base_to_correct_read_directly(Correct_dumy* dumy, char base) +{ + + 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'; + +} + inline void add_base_to_correct_read(Correct_dumy* dumy, char base, int is_error) { ///deletion就不要管 @@ -1632,6 +1648,223 @@ inline void add_segment_to_correct_read(Correct_dumy* dumy, char* segment, long dumy->corrected_read[dumy->corrected_read_length] = '\0'; } + + + + +///返回下一个backbone节点上的ID +long long inline add_path_to_correct_read(Graph* backbone, Correct_dumy* dumy, long long currentNodeID, +long long type, long long edgeID) +{ + //long long i; + long long nodeID; + + ///包括匹配和误配两种情况 + if (type == MISMATCH) + { + ///这是match的情况 + if(backbone->g_nodes.list[currentNodeID].mismatch_edges.list[edgeID].length == 0) + { + nodeID = backbone->g_nodes.list[currentNodeID].mismatch_edges.list[edgeID].out_node; + add_base_to_correct_read_directly(dumy, backbone->g_nodes.list[nodeID].base); + ///match所以dumy->corrected_base不要+1 + + /***********需要注释掉********* */ + if (nodeID != currentNodeID + 1) + { + fprintf(stderr, "error match\n"); + } + /***********需要注释掉********* */ + + return nodeID; + } + else ///这是mismatch的情况 + { + nodeID = backbone->g_nodes.list[currentNodeID].mismatch_edges.list[edgeID].out_node; + add_base_to_correct_read_directly(dumy, backbone->g_nodes.list[nodeID].base); + dumy->corrected_base++; + + ///这种中间节点只有一个元素，所以直接list[0] + nodeID = backbone->g_nodes.list[nodeID].mismatch_edges.list[0].out_node; + /***********需要注释掉********* */ + if (nodeID != currentNodeID + 1) + { + fprintf(stderr, "error mismatch\n"); + } + /***********需要注释掉********* */ + + return nodeID; + } + } + else if (type == DELETION) + { + nodeID = backbone->g_nodes.list[currentNodeID].deletion_edges.list[edgeID].out_node; + + /***********需要注释掉********* */ + if (!(nodeID >= backbone->s_start_nodeID && nodeID <= backbone->s_end_nodeID)) + { + fprintf(stderr, "error deletion 1\n"); + } + if (nodeID <= currentNodeID) + { + fprintf(stderr, "error deletion 2\n"); + } + + /***********需要注释掉********* */ + dumy->corrected_base += nodeID - currentNodeID; + return nodeID; + } + else if (type == INSERTION) + { + ///这个一定要变成0 + backbone->g_nodes.list[currentNodeID].num_insertions = 0; + + nodeID = backbone->g_nodes.list[currentNodeID].insertion_edges.list[edgeID].out_node; + long long step = backbone->g_nodes.list[currentNodeID].insertion_edges.list[edgeID].length; + long long i; + for (i = 0; i < step; i++) + { + add_base_to_correct_read_directly(dumy, backbone->g_nodes.list[nodeID].base); + ///只有一条边 + nodeID = backbone->g_nodes.list[nodeID].insertion_edges.list[0].out_node; + } + dumy->corrected_base += step; + + /***********需要注释掉********* */ + if (nodeID != currentNodeID) + { + fprintf(stderr, "error insertion\n"); + } + /***********需要注释掉********* */ + return nodeID; + } + else + { + fprintf(stderr, "error type\n"); + } + +} + +void get_seq_from_Graph(Graph* backbone, 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_edge; + long long total_count; + long long nodeID; + char current_base; + long long current_weight; + + currentNodeID = backbone->s_start_nodeID; + + while (currentNodeID != backbone->s_end_nodeID) + { + total_count = 0; + max_count = -1; + max_type = -1; + max_edge = -1; + + + ///假如这是个backbone节点 + ///有三类出边 + ///1. mismatch_edges 2. insertion_edges 3. deletion_edges + if (currentNodeID >= backbone->s_start_nodeID && currentNodeID <= backbone->s_end_nodeID) + { + + ///mismatch_edges + for (i = 0; i < backbone->g_nodes.list[currentNodeID].mismatch_edges.length; i++) + { + if (backbone->g_nodes.list[currentNodeID].num_insertions != 0) + { + current_weight = backbone->g_nodes.list[currentNodeID].mismatch_edges.list[i].weight - + backbone->g_nodes.list[currentNodeID].mismatch_edges.list[i].num_insertions; + } + else + { + current_weight = backbone->g_nodes.list[currentNodeID].mismatch_edges.list[i].weight; + } + + + total_count = total_count + current_weight; + + ///match + ///match要处理插入的情况 + ///如果这里有insertion, 这个节点会过两遍 + ///第一遍num_insertions > 0, 第二遍num_insertions=0 + if (current_weight > max_count) + { + max_count = current_weight; + max_edge = i; + max_type = MISMATCH; + } + } + + ///insertion_edges + if (backbone->g_nodes.list[currentNodeID].num_insertions != 0) + { + for (i = 0; i < backbone->g_nodes.list[currentNodeID].insertion_edges.length; i++) + { + total_count = total_count + backbone->g_nodes.list[currentNodeID].insertion_edges.list[i].weight; + + if (backbone->g_nodes.list[currentNodeID].insertion_edges.list[i].weight > max_count) + { + max_count = backbone->g_nodes.list[currentNodeID].insertion_edges.list[i].weight; + max_edge = i; + max_type = INSERTION; + } + } + } + + + + ///deletion_edges + for (i = 0; i < backbone->g_nodes.list[currentNodeID].deletion_edges.length; i++) + { + total_count = total_count + backbone->g_nodes.list[currentNodeID].deletion_edges.list[i].weight; + + if (backbone->g_nodes.list[currentNodeID].deletion_edges.list[i].weight > max_count) + { + max_count = backbone->g_nodes.list[currentNodeID].deletion_edges.list[i].weight; + max_edge = i; + max_type = DELETION; + } + } + + ///这种情况下矫正 + if(max_count >= total_count*CORRECT_THRESHOLD) + { + currentNodeID = add_path_to_correct_read(backbone, dumy, currentNodeID, max_type, max_edge); + } + else ///不矫正, 直接取下一个backbone节点 + { + currentNodeID++; + add_base_to_correct_read_directly(dumy, backbone->g_nodes.list[currentNodeID].base); + } + + + ///fprintf(stderr, "currentNodeID: %d, max_type: %d\n", currentNodeID, max_type); + + + } + else ///非backbone节点就会出错了 + { + fprintf(stderr, "error\n"); + } + + } + +} + + + + +/** ///从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) { @@ -1836,9 +2069,232 @@ void get_seq_from_Graph(Graph* backbone, long long backbone_start, long long bac add_base_to_correct_read(dumy, current_base, max_type); } +**/ + + + + + + + +/** +///从backbone_start遍历到backbone_end节点，生成出来的seq要接着放到dumy->corrected_read中 +void get_seq_from_Graph_Len2(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; + char buffer[2]; + + 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) + { + ///fprintf(stderr, "error\n"); + 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节点，就无脑看出边 + { + if (backbone->g_nodes.list[currentNodeID].outcome_edges.length!=1) + { + fprintf(stderr, "error000\n"); + } + + ///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节点上不达标很麻烦...,只能选最大的了 + { + ///因为现在每个insert节点只有一个出边，且这个出边到backbone + fprintf(stderr, "error111\n"); + } + } + + + if (max_count <= 0) + { + fprintf(stderr, "error\n"); + } + + if (current_base < 'A') + { + buffer[0] = s_H[(current_base >> 2) & ((uint8_t)3)]; + buffer[1] = s_H[current_base & ((uint8_t)3)]; + add_base_to_correct_read(dumy, buffer[0], max_type); + add_base_to_correct_read(dumy, buffer[1], max_type); + } + else + { + 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) { @@ -1902,31 +2358,19 @@ overlap_region_alloc* overlap_list, Correct_dumy* dumy, All_reads* R_INF, Graph* y_string, y_length, &(overlap_list->list[overlapID].w_list[windowID].cigar), startNodeID, endNodeID); } + get_seq_from_Graph(g, dumy); - get_seq_from_Graph(g, startNodeID, endNodeID, dumy); + ///get_seq_from_Graph(g, startNodeID, endNodeID, dumy); + ///get_seq_from_Graph_Len2(g, startNodeID, endNodeID, dumy); + + + + ///debug_graph(g, backbone_length); + + /** - 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 @@ -2160,11 +2604,14 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, generate_consensus(overlap_list, R_INF, g_read, dumy, g); - ///fprintf(stderr, "length: %lld, corrected_base: %lld\n", g_read->length, dumy->corrected_base); /** + ///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) { @@ -2193,6 +2640,15 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, } edlibFreeAlignResult(result); + + for (i = 0; i < dumy->corrected_read_length; i++) + { + if (dumy->corrected_read[i] != 'A' && dumy->corrected_read[i] != 'C' && dumy->corrected_read[i] != 'G' && dumy->corrected_read[i] != 'T') + { + fprintf(stderr, "error\n"); + } + + } **/ /************需要注释掉********* */ diff --git a/Correct.h b/Correct.h index a83ccf6..f727669 100644 --- a/Correct.h +++ b/Correct.h @@ -5,10 +5,12 @@ #include "Levenshtein_distance.h" #include "POA.h" -#define CORRECT_THRESHOLD 0.6 +#define CORRECT_THRESHOLD 0.65 #define MIN_COVERAGE_THRESHOLD 4 - - +#define CORRECT_INDEL_LENGTH 2 +#define MISMATCH 1 +#define INSERTION 2 +#define DELETION 3 #define MAX(x, y) ((x >= y)?x:y) #define MIN(x, y) ((x <= y)?x:y) @@ -47,6 +49,7 @@ void pre_filter_by_nearby(k_mer_pos* new_n_list, k_mer_pos* old_n_list, uint64_t All_reads* R_INF, Correct_dumy* dumy, uint64_t* new_n_length); void pre_filter_by_nearby_single(k_mer_pos* new_n_list, k_mer_pos* old_n_list, uint64_t n_length, uint64_t n_end_pos, UC_Read* g_read, All_reads* R_INF, Correct_dumy* dumy, uint64_t* new_n_length); +void get_seq_from_Graph(Graph* backbone, Correct_dumy* dumy); /**********************for prefilter************************ */ void destory_k_mer_pos_list_alloc_prefilter(k_mer_pos_list_alloc* list); diff --git a/Hash_Table.h b/Hash_Table.h index 78219e1..52b2fdf 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -15,14 +15,22 @@ typedef khash_t(POS64) Pos_Table; ///#define WINDOW 350 ///#define THRESHOLD 14 + #define WINDOW 375 #define THRESHOLD 15 -///#define THRESHOLD_RATE 0.04 #define THRESHOLD_RATE 0.04 +#define OVERLAP_THRESHOLD 0.9 + +/** +#define WINDOW 500 +#define THRESHOLD 15 +#define THRESHOLD_RATE 0.03 +#define OVERLAP_THRESHOLD 0.95 +**/ + #define GROUP_SIZE 4 ///最长是10M10D10M10D10M这种 #define CIGAR_MAX_LENGTH THRESHOLD*2+2 -#define OVERLAP_THRESHOLD 0.9 typedef struct { diff --git a/POA.cpp b/POA.cpp index c649292..3f1819f 100644 --- a/POA.cpp +++ b/POA.cpp @@ -1,79 +1,12 @@ #include "POA.h" #include #include "Correct.h" +#include "Process_Read.h" #define INIT_EDGE_SIZE 50 #define INCREASE_EDGE_SIZE 5 #define INIT_NODE_SIZE 16000 -void check_addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length, long long startID, long long endID) -{ - long long reverse_startID, reverse_endID; - reverse_startID = startID; - reverse_endID = endID; - long long i = 0; - if (endID - startID + 1 != g_read_length) - { - fprintf(stderr, "ERROR Length...\n"); - fprintf(stderr, "startID: %lld, endID: %lld\n", startID, endID); - fprintf(stderr, "g_read_length: %lld\n", g_read_length); - } - - if (g_read_length == 0) - { - return; - } - - - while (1) - { - if(g->g_nodes.list[startID].base != g_read_seq[i]) - { - fprintf(stderr, "i: %llu, ERROR Node Base...\n", i); - - } - - if(g->g_nodes.list[startID].outcome_edges.length == 0) - { - break; - } - - startID = g->g_nodes.list[startID].outcome_edges.list[0].out_node; - - i++; - } - - if (startID != endID) - { - fprintf(stderr, "ERROR End Node Base\n"); - } - - - i = g_read_length - 1; - while (1) - { - if(g->g_nodes.list[reverse_endID].base != g_read_seq[i]) - { - fprintf(stderr, "i: %llu, g_read_length: %llu, ERROR Node Base...\n", i, g_read_length); - - } - - if(g->g_nodes.list[reverse_endID].income_edges.length == 0) - { - break; - } - - reverse_endID = g->g_nodes.list[reverse_endID].income_edges.list[0].in_node; - - i--; - } - - if (reverse_startID != reverse_endID) - { - fprintf(stderr, "ERROR Start Node Base\n"); - } - -} void init_Edge_alloc(Edge_alloc* list) @@ -101,7 +34,7 @@ void destory_Edge_alloc(Edge_alloc* list) free(list->list); } -void append_Edge_alloc(Edge_alloc* list, uint64_t in_node, uint64_t out_node, uint64_t weight) +void append_Edge_alloc(Edge_alloc* list, uint64_t in_node, uint64_t out_node, uint64_t weight, uint64_t length) { if (list->length + 1 > list->size) { @@ -112,6 +45,8 @@ void append_Edge_alloc(Edge_alloc* list, uint64_t in_node, uint64_t out_node, u list->list[list->length].in_node = in_node; list->list[list->length].out_node = out_node; list->list[list->length].weight = weight; + list->list[list->length].length = length; + list->list[list->length].num_insertions = 0; list->length++; } @@ -119,10 +54,6 @@ void append_Edge_alloc(Edge_alloc* list, uint64_t in_node, uint64_t out_node, u - - - - void init_Node_alloc(Node_alloc* list) { list->size = INIT_NODE_SIZE; @@ -139,15 +70,10 @@ void init_Node_alloc(Node_alloc* list) long long i; for (i = 0; i < list->size; i++) { - list->list[i].income_edges.list=NULL; - 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; - + list->list[i].insertion_edges.list=NULL; + list->list[i].mismatch_edges.list=NULL; + list->list[i].deletion_edges.list=NULL; + } } void destory_Node_alloc(Node_alloc* list) @@ -155,13 +81,10 @@ void destory_Node_alloc(Node_alloc* list) uint64_t i =0; for (i = 0; i < list->length; i++) { - destory_Edge_alloc(&list->list[i].income_edges); - destory_Edge_alloc(&list->list[i].outcome_edges); - destory_Edge_alloc(&list->list[i].alignedTo_Nodes); + destory_Edge_alloc(&list->list[i].deletion_edges); + destory_Edge_alloc(&list->list[i].insertion_edges); + destory_Edge_alloc(&list->list[i].mismatch_edges); } - 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); @@ -176,13 +99,11 @@ void clear_Node_alloc(Node_alloc* list) uint64_t i =0; for (i = 0; i < list->length; i++) { - clear_Edge_alloc(&list->list[i].income_edges); - clear_Edge_alloc(&list->list[i].outcome_edges); - clear_Edge_alloc(&list->list[i].alignedTo_Nodes); + clear_Edge_alloc(&list->list[i].insertion_edges); + clear_Edge_alloc(&list->list[i].mismatch_edges); + clear_Edge_alloc(&list->list[i].deletion_edges); } - 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; } @@ -202,18 +123,19 @@ uint64_t append_Node_alloc(Node_alloc* list, char base) for (; i < list->size; i++) { - list->list[i].income_edges.list=NULL; - list->list[i].outcome_edges.list=NULL; - list->list[i].alignedTo_Nodes.list=NULL; + list->list[i].deletion_edges.list=NULL; + list->list[i].insertion_edges.list=NULL; + list->list[i].mismatch_edges.list=NULL; } } 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); + list->list[list->length].num_insertions = 0; + init_Edge_alloc(&list->list[list->length].deletion_edges); + init_Edge_alloc(&list->list[list->length].insertion_edges); + init_Edge_alloc(&list->list[list->length].mismatch_edges); list->length++; @@ -264,26 +186,8 @@ void clear_Graph(Graph* g) -uint64_t inline add_Node_Graph(Graph* g, char base) -{ - return append_Node_alloc(&g->g_nodes, base); -} -void inline add_Edge_Graph(Graph* g, uint64_t start, uint64_t end, uint64_t weight) -{ - if (start >= g->g_nodes.length || end >= g->g_nodes.length) - { - fprintf(stderr, "Not existing nodes ..."); - exit(0); - } - - ///对起始节点加出边 - append_Edge_alloc(&g->g_nodes.list[start].outcome_edges, start, end, weight); - ///对结束节点加入边 - append_Edge_alloc(&g->g_nodes.list[end].income_edges, start, end, weight); -} - void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length, long long* startID, long long* endID) { long long firstID, lastID, nodeID, i; @@ -292,7 +196,13 @@ void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length, if(g_read_length == 0) return; - + + ///start node + nodeID = add_Node_Graph(g, 'S'); + firstID = nodeID; + lastID = nodeID; + + for (i = 0; i < g_read_length; i++) { nodeID = add_Node_Graph(g, g_read_seq[i]); @@ -305,8 +215,13 @@ void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length, } if (lastID != -1) { + /** ///0是match边 add_Edge_Graph(g, lastID, nodeID, 0); + **/ + ///只有match边长度是0 + ///mismatch边长度都是1 + append_Edge_alloc(&(g->g_nodes.list[lastID].mismatch_edges), lastID, nodeID, 1, 0); } lastID = nodeID; @@ -314,212 +229,23 @@ void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length, *startID = firstID; *endID = lastID; - -} -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; + g->s_start_nodeID = firstID; + g->s_end_nodeID = lastID; } -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; @@ -527,42 +253,41 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri int operation; int operationLen; int i; + int last_operation = -1; - + + ///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多字符） - ///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) + ///这种情况代表匹配和mismatch + if (operation == 0 || operation == 1) { for (i = 0; i < operationLen; i++) { - backbone->g_nodes.list[currentNodeID].weight++; - + //backbone->g_nodes.list[currentNodeID].weight++; + ///前面是插入，后面有可能是误配，也有可能是匹配 + add_mismatchEdge_weight(backbone, currentNodeID, y_string[y_i], last_operation); 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; - } + }///insertion 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); + if (operationLen <= CORRECT_INDEL_LENGTH) + { + add_insertionEdge_weight(backbone, currentNodeID, y_string + y_i, operationLen); + backbone->g_nodes.list[currentNodeID].num_insertions++; + } y_i += operationLen; } else if (operation == 3) @@ -570,12 +295,25 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri ///3是y缺字符（x多字符），也就是backbone多字符 ///这个相当于在backbone对应字符处变成了‘——’ ///因此可以用mismatch类似的方法处理 - add_deletion_to_backbone(backbone, ¤tNodeID, operationLen); + if (operationLen <= CORRECT_INDEL_LENGTH) + { + ///add_deletion_to_backbone(backbone, ¤tNodeID, operationLen); + ///在编辑距离中，前面是个插入，后面是个删除，这种情况是不存在的 + ///为了保险要不还给他加上吧 + ///先不加 + add_deletionEdge_weight(backbone, currentNodeID, operationLen); + } + + + currentNodeID += operationLen; x_i += operationLen; } + + last_operation = operation; cigar_i++; } + /** @@ -597,11 +335,140 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri } **/ +} + + +void debug_graph(Graph* g, long long backbone_length) +{ + long long i = 0; + + if (g->s_start_nodeID != 0 || g->s_end_nodeID != backbone_length) + { + fprintf(stderr, "error\n"); + } + + + for (i = g->s_start_nodeID; i <= g->s_end_nodeID; i++) + { + if(g->g_nodes.list[i].weight != 1) + { + fprintf(stderr, "error node weight\n"); + } + + if(g->g_nodes.list[i].mismatch_edges.length > 4) + { + fprintf(stderr, "error mismatch_edges\n"); + } + + if(g->g_nodes.list[i].mismatch_edges.length < 1 && i != g->s_end_nodeID) + { + fprintf(stderr, "i: %d, error mismatch_edges: %d\n", i, g->g_nodes.list[i].mismatch_edges.length); + } + } + + + for (i = 0; i < g->g_nodes.length; i++) + { + if(g->g_nodes.list[i].ID < g->s_start_nodeID || g->g_nodes.list[i].ID > g->s_end_nodeID) + { + if (g->g_nodes.list[i].deletion_edges.length + + g->g_nodes.list[i].insertion_edges.length + + g->g_nodes.list[i].mismatch_edges.length + != 1) + { + fprintf(stderr, "g->s_start_nodeID: %lld\n", + g->s_start_nodeID); + fprintf(stderr, "g->s_end_nodeID: %lld\n", + g->s_end_nodeID); + fprintf(stderr, "deletion_edges_length: %lld\n", + g->g_nodes.list[i].deletion_edges.length); + fprintf(stderr, "insertion_edges_length: %lld, \n", + g->g_nodes.list[i].insertion_edges.length); + fprintf(stderr, "g->g_nodes.list[i].insertion_edges.list[0].length: %lld, \n", + g->g_nodes.list[i].insertion_edges.list[0].length); + fprintf(stderr, "g->g_nodes.list[i].insertion_edges.list[1].length: %lld, \n", + g->g_nodes.list[i].insertion_edges.list[1].length); + + fprintf(stderr, "mismatch_edges_length: %lld\n", + g->g_nodes.list[i].mismatch_edges.length); + } + else + { + ///不是0肯定是1 + if (g->g_nodes.list[i].deletion_edges.length != 0) + { + long long step = g->g_nodes.list[i].deletion_edges.list[0].length; + long long nodeID = i; + for (int j = 0; j < step; j++) + { + nodeID = g->g_nodes.list[nodeID].deletion_edges.list[0].out_node; + } + + nodeID = g->g_nodes.list[nodeID].deletion_edges.list[0].out_node; + + if (nodeID < g->s_start_nodeID || nodeID > g->s_end_nodeID) + { + fprintf(stderr, "error\n"); + } + } + + if (g->g_nodes.list[i].insertion_edges.length != 0) + { + + long long step = g->g_nodes.list[i].insertion_edges.list[0].length; + long long nodeID = i; + + + for (int j = 0; j < step; j++) + { + nodeID = g->g_nodes.list[nodeID].insertion_edges.list[0].out_node; + } + + nodeID = g->g_nodes.list[nodeID].insertion_edges.list[0].out_node; + + if ((nodeID < g->s_start_nodeID || nodeID > g->s_end_nodeID)) + { + fprintf(stderr, "error: step: %d\n", step); + } + + } + + + if (g->g_nodes.list[i].mismatch_edges.length != 0) + { + + long long step = g->g_nodes.list[i].mismatch_edges.list[0].length; + long long nodeID = i; + + for (int j = 0; j < step; j++) + { + nodeID = g->g_nodes.list[nodeID].mismatch_edges.list[0].out_node; + } + + nodeID = g->g_nodes.list[nodeID].mismatch_edges.list[0].out_node; + + if (nodeID < g->s_start_nodeID || nodeID > g->s_end_nodeID) + { + fprintf(stderr, "error\n"); + } + + } + + } + + + } + } + + + + } 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; @@ -746,168 +613,10 @@ void Graph_debug(Graph* backbone, long long currentNodeID, char* x_string, long { 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) -{ - long long startNodeID, endNodeID; - ///第一条序列是read本身，g_read里存的是反向互补，所以首先要恢复回正向 - ///这一步可以优化掉 - reverse_complement(g_read->seq, g_read->length); - addUnmatchedSeqToGraph(g, g_read->seq, g_read->length, &startNodeID, &endNodeID); - - /** - if (startNodeID!=0||endNodeID!=g_read->length-1) - { - fprintf(stderr, "Error startNodeID or endNodeID ...\n"); - } - check_addUnmatchedSeqToGraph(g, g_read->seq, g_read->length, startNodeID, endNodeID); **/ - - get_Topo_Sort_Order(&g->g_nodes, 0); - - uint64_t i = 0; - for ( i = 0; i < overlap_list->length; i++) - { - /** - if(overlap_list->list[i].x_id == overlap_list->list[i].y_id) - { - fprintf(stderr, "Error x_id or y_id ...\n"); - } - **/ - - if (overlap_list->list[i].y_pos_strand) - { - recover_UC_Read_RC(g_read, R_INF, overlap_list->list[i].y_id); - } - else - { - recover_UC_Read(g_read, R_INF, overlap_list->list[i].y_id); - } - - } - -} - -void topologicalSortDFS(Node_alloc* list, uint64_t nodeID) -{ - list->sort.visit[nodeID] = 1; - - long long i; - uint64_t out_nodeID; - - for (i = 0; i < list->list[nodeID].outcome_edges.length; i++) - { - out_nodeID = list->list[nodeID].outcome_edges.list[i].out_node; - if (list->sort.visit[out_nodeID] == 0) - { - topologicalSortDFS(list, out_nodeID); - } - } - - list->sort.length--; - list->sort.list[list->sort.length] = nodeID; -} - -#define INIT_STACK(stack) stack.iterative_i = 0; -#define PUSH(stack, nodeID, time) stack.iterative_buffer[stack.iterative_i]=nodeID;\ -stack.iterative_buffer_visit[stack.iterative_i++]=time; -#define IF_EMPTY(stack) (stack.iterative_i == 0) -#define POP(stack, nodeID, time) --stack.iterative_i;nodeID = stack.iterative_buffer[stack.iterative_i];\ -time = stack.iterative_buffer_visit[stack.iterative_i]; - - -void topologicalSortDFS_Iterative(Node_alloc* list, uint64_t nodeID) -{ - long long i; - uint64_t out_nodeID; - int flag; - - INIT_STACK(list->sort); - PUSH(list->sort, nodeID, 0); - - while (!IF_EMPTY(list->sort)) - { - POP(list->sort, nodeID, flag); - ///flag == 1说明是第二次访问; flag == 0说明是第一次访问 - if (flag) - { - list->sort.length--; - list->sort.list[list->sort.length] = nodeID; - continue; - } - - list->sort.visit[nodeID] = 1; - PUSH(list->sort, nodeID, 1); - - for (i = 0; i < list->list[nodeID].outcome_edges.length; i++) - { - out_nodeID = list->list[nodeID].outcome_edges.list[i].out_node; - if (list->sort.visit[out_nodeID] == 0) - { - PUSH(list->sort, out_nodeID, 0); - } - } - - } -} - -uint64_t* get_Topo_Sort_Order(Node_alloc* list, int need_sort) -{ - long long i = 0; - list->sort.length = list->length; - if (list->length > list->sort.size) - { - list->sort.size = list->length; - list->sort.list = (uint64_t*)realloc(list->sort.list, sizeof(uint64_t)*list->sort.size); - list->sort.visit = (uint8_t*)realloc(list->sort.visit, sizeof(uint8_t)*list->sort.size); - - list->sort.iterative_buffer - = (uint64_t*)realloc(list->sort.iterative_buffer, sizeof(uint64_t)*list->sort.size); - list->sort.iterative_buffer_visit - = (uint8_t*)realloc(list->sort.iterative_buffer_visit, sizeof(uint8_t)*list->sort.size); - } - - if (!need_sort) - { - ///可以循环展开, 作用微乎其微 - for (i = 0; i < list->length; i++) - { - list->sort.list[i] = i; - } - - } - else - { - memset(list->sort.visit, 0 , list->length); - for (i = 0; i < list->length; i++) - { - if(list->sort.visit[i] == 0) - { - ///topologicalSortDFS(list, i); - topologicalSortDFS_Iterative(list, i); - } - } - - - /** - for (i = 0; i < list->length; i++) - { - if (list->sort.list[i] != i) - { - fprintf(stderr, "ERROR Sort ....\n"); - } - } - **/ - - - } - - return list->sort.list; } + + + + diff --git a/POA.h b/POA.h index a392341..0e2e39e 100644 --- a/POA.h +++ b/POA.h @@ -29,6 +29,9 @@ typedef struct uint64_t out_node; ///0是match，1是mismatch，2是x缺字符（y多字符），而3是y缺字符（x多字符） uint64_t weight; + uint64_t num_insertions; + ///这条路径上到backbone节点之前总共有多少节点 + uint64_t length; } Edge; typedef struct @@ -42,10 +45,13 @@ typedef struct { uint64_t ID; uint64_t weight; + ///记录的是以当前节点为尾的deletion个数 + uint64_t num_insertions; char base; - Edge_alloc income_edges; - Edge_alloc outcome_edges; - Edge_alloc alignedTo_Nodes; + Edge_alloc mismatch_edges; + Edge_alloc deletion_edges; + Edge_alloc insertion_edges; + } Node; typedef struct @@ -62,7 +68,7 @@ typedef struct typedef struct { - Node total_start; + ///has a indivial start node 0 Node* list; topo_Sorting_buffer sort; uint64_t size; @@ -88,7 +94,7 @@ typedef struct void init_Edge_alloc(Edge_alloc* list); void clear_Edge_alloc(Edge_alloc* list); void destory_Edge_alloc(Edge_alloc* list); -void append_Edge_alloc(Edge_alloc* list, uint64_t in_node, uint64_t out_node, uint64_t weight); +void append_Edge_alloc(Edge_alloc* list, uint64_t in_node, uint64_t out_node, uint64_t weight, uint64_t length); void init_Node_alloc(Node_alloc* list); void destory_Node_alloc(Node_alloc* list); @@ -108,5 +114,285 @@ void Perform_POA(Graph* g, overlap_region_alloc* overlap_list, All_reads* R_INF, 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); +void debug_graph(Graph* g, long long backbone_length); + +uint64_t inline add_Node_Graph(Graph* g, char base) +{ + return append_Node_alloc(&g->g_nodes, base); +} + + +///仅仅用于误配边 +inline void add_mismatchEdge_weight(Graph* g, uint64_t in_node, char base, int last_operation) +{ + long long i = 0; + long long nodeID; + Edge_alloc* edge = &(g->g_nodes.list[in_node].mismatch_edges); + + for (i = 0; i < edge->length; i++) + { + nodeID = edge->list[i].out_node; + if(g->g_nodes.list[nodeID].base == base) + { + edge->list[i].weight++; + ///如果上一个操作是insertion + if (last_operation == 2) + { + edge->list[i].num_insertions++; + } + + break; + } + } + + ///说明不存在这么一条边 + if (i == edge->length) + { + nodeID = add_Node_Graph(g, base); + + ///只有match边长度是0 + ///mismatch边长度都是1 + append_Edge_alloc(edge, in_node, nodeID, 1, 1); + ///如果上一个操作是insertion + if (last_operation == 2) + { + edge->list[edge->length - 1].num_insertions++; + } + + ///将新节点的mismatch_edges连到backbone上 + append_Edge_alloc(&(g->g_nodes.list[nodeID].mismatch_edges), nodeID, in_node + 1, 1, 0); + } + ///获得节点的mismatch_edges长度为1，其他均为0 + +} + + + +inline void add_single_deletionEdge_weight(Graph* g, long long alignNodeID, long long nextNodeID, uint64_t edge_length) +{ + long long i = 0; + long long nodeID; + Edge_alloc* edge = &(g->g_nodes.list[alignNodeID].deletion_edges); + + for (i = 0; i < edge->length; i++) + { + nodeID = edge->list[i].out_node; + if(nodeID == nextNodeID) + { + edge->list[i].weight++; + break; + } + } + + ///说明不存在这么一条边 + if (i == edge->length) + { + append_Edge_alloc(edge, alignNodeID, nextNodeID, 1, edge_length); + } +} + +inline void add_deletionEdge_weight(Graph* g, long long alignNodeID, long long deletion_length) +{ + if (deletion_length == 1) + { + add_single_deletionEdge_weight(g, alignNodeID, alignNodeID + 1, 0); + } + else if (deletion_length == 2) + { + add_single_deletionEdge_weight(g, alignNodeID, alignNodeID + 1, 0); + add_single_deletionEdge_weight(g, alignNodeID, alignNodeID + 2, 0); + add_single_deletionEdge_weight(g, alignNodeID + 1, alignNodeID + 2, 0); + } + else + { + fprintf(stderr, "too long deletion!\n"); + } +} + + +inline int getEdge(Graph* g, Edge_alloc* edge, uint64_t edge_length, char base) +{ + long long i = 0; + long long nodeID; + + for (i = 0; i < edge->length; i++) + { + if (edge->list[i].length == edge_length) + { + nodeID = edge->list[i].out_node; + if(g->g_nodes.list[nodeID].base == base) + { + return i; + } + } + } + + return -1; +} + + +inline int get_insertion_Edges(Graph* g, Edge_alloc* edge, uint64_t edge_length, char* bases) +{ + long long i = 0; + long long nodeID; + long long edgeID; + + if (edge_length < 1) + { + return -1; + } + + + edgeID = getEdge(g, edge, edge_length, bases[0]); + + long long return_edgeID = edgeID; + + if(edgeID == -1) + { + return -1; + } + + + Edge_alloc* new_edge = edge; + + for (i = 1; i < edge_length; i++) + { + nodeID = new_edge->list[edgeID].out_node; + new_edge = &(g->g_nodes.list[nodeID].insertion_edges); + edgeID = getEdge(g, new_edge, edge_length - i, bases[i]); + if(edgeID == -1) + { + return -1; + } + } + + return edgeID; +} + + + +inline int create_insertion_Edges(Graph* g, long long alignNodeID, uint64_t edge_length, char* bases) +{ + long long i = 0; + long long nodeID; + ///最后应该连回原节点 + ///long long backboneID = alignNodeID + 1; + long long backboneID = alignNodeID; + + + if (edge_length < 1) + { + return -1; + } + + + nodeID = add_Node_Graph(g, bases[0]); + ///将新加入的节点通过insertion_edges接到alignNodeID上 + append_Edge_alloc(&(g->g_nodes.list[alignNodeID].insertion_edges), alignNodeID, nodeID, 1, edge_length); + + alignNodeID = nodeID; + + for (i = 1; i < edge_length; i++) + { + nodeID = add_Node_Graph(g, bases[i]); + ///将新加入的节点通过insertion_edges接到alignNodeID上 + append_Edge_alloc(&(g->g_nodes.list[alignNodeID].insertion_edges), alignNodeID, nodeID, 1, edge_length - i); + alignNodeID = nodeID; + } + + append_Edge_alloc(&(g->g_nodes.list[alignNodeID].insertion_edges), alignNodeID, backboneID, 1, 0); + +} + +inline void add_insertionEdge_weight(Graph* g, long long alignNodeID, char* insert, long long insert_length) +{ + long long i = 0; + long long nodeID; + long long edgeID; + Edge_alloc* edge = &(g->g_nodes.list[alignNodeID].insertion_edges); + + if (insert_length == 1) + { + edgeID = getEdge(g, edge, 1, insert[0]); + if (edgeID != -1) + { + ///这条路均只有一个出度 + edge->list[edgeID].weight++; + } + else ///不存在这么一条边 + { + nodeID = add_Node_Graph(g, insert[0]); + append_Edge_alloc(edge, alignNodeID, nodeID, 1, 1); + ///将新加入的节点通过insertion_edges接回backbone上 + ///应该连回到原节点，而不是原节点的下一个节点 + ///append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID + 1, 1, 0); + append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID, 1, 0); + } + } + else if (insert_length == 2) + { + /*******************第0个字符********************* */ + edgeID = getEdge(g, edge, 1, insert[0]); + if (edgeID != -1) + { + ///这条路均只有一个出度 + edge->list[edgeID].weight++; + } + else ///不存在这么一条边 + { + nodeID = add_Node_Graph(g, insert[0]); + append_Edge_alloc(edge, alignNodeID, nodeID, 1, 1); + ///将新加入的节点通过insertion_edges接回backbone上 + ///应该连回到原节点，而不是原节点的下一个节点 + ///append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID + 1, 1, 0); + append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID, 1, 0); + } + /*******************第0个字符********************* */ + + /*******************第1个字符********************* */ + if (insert[1] != insert[0]) + { + edgeID = getEdge(g, edge, 1, insert[1]); + if (edgeID != -1) + { + ///这条路均只有一个出度 + edge->list[edgeID].weight++; + } + else ///不存在这么一条边 + { + nodeID = add_Node_Graph(g, insert[1]); + append_Edge_alloc(edge, alignNodeID, nodeID, 1, 1); + ///将新加入的节点通过insertion_edges接回backbone上 + ///应该连回到原节点，而不是原节点的下一个节点 + ///append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID + 1, 1, 0); + append_Edge_alloc(&(g->g_nodes.list[nodeID].insertion_edges), nodeID, alignNodeID, 1, 0); + } + } + /*******************第1个字符********************* */ + + /**********************两个字符******************* */ + + edgeID = get_insertion_Edges(g, edge, 2, insert); + if (edgeID != -1) + { + ///这条路均只有一个出度 + edge->list[edgeID].weight++; + } + else + { + create_insertion_Edges(g, alignNodeID, insert_length, insert); + } + + /**********************两个字符******************* */ + } + else + { + fprintf(stderr, "too long insertion\n"); + } + + + +} + #endif \ No newline at end of file