#ifndef __POA_PARSER__ #define __POA_PARSER__ #define __STDC_LIMIT_MACROS #include #include "Hash_Table.h" #include "Process_Read.h" typedef struct { long long beg; ///end is the index of next input data, instead of the index of last data long long end; long long length; long long size; long long* buffer; } Queue; inline void init_Queue(Queue* q) { q->beg = 0; q->end = 0; q->length = 0; q->size = 20; q->buffer = (long long*)malloc(sizeof(long long)*q->size); } inline void clear_Queue(Queue* q) { q->beg = 0; q->end = 0; q->length = 0; } inline void destory_Queue(Queue* q) { free(q->buffer); } inline int is_empty_Queue(Queue* q) { ///end is the index of next input data, instead of the index of last data if(q->beg == q->end) { return 1; } else { return 0; } } inline int is_full_Queue(Queue* q) { ///end is the index of next input data, instead of the index of last data if(q->end < q->size) { return 0; } else { return 1; } } inline void push_to_Queue(Queue* q, long long nodeID) { if(is_full_Queue(q)) { long long move_length = q->beg; ///end is the index of next input data, instead of the index of last data long long current_length = q->end - q->beg; ///recalloc directly if(move_length == 0) { q->size = q->size * 2; q->buffer = (long long*)realloc(q->buffer, q->size*sizeof(long long)); } else { ///won't overlap if(current_length <= move_length) { memcpy(q->buffer, q->buffer+q->beg, sizeof(long long)*current_length); } else///may overlap { memmove(q->buffer, q->buffer+q->beg, sizeof(long long)*current_length); } q->beg = 0; q->end = current_length; } } q->buffer[q->end] = nodeID; q->end++; } inline int pop_from_Queue(Queue* q, long long* nodeID) { if(is_empty_Queue(q)) { (*nodeID) = -1; return 0; } else { (*nodeID) = q->buffer[q->beg]; q->beg++; return 1; } } typedef struct { uint64_t in_node; uint64_t out_node; ///0 is match,1 is mismatch,2 means y has more bases, 3 means x has more bases uint64_t weight; uint64_t num_insertions; uint64_t length; uint64_t self_edge_ID; uint64_t reverse_edge_ID; } Edge; typedef struct { Edge* list; uint64_t size; uint64_t length; uint64_t delete_length; } Edge_alloc; #define Real_Length(X) ((X).length - (X).delete_length) #define Input_Edges(Node) ((Node).insertion_edges) #define Output_Edges(Node) ((Node).deletion_edges) #define G_Node(G, Node) ((G).g_nodes.list[(Node)]) #define If_Node_Exist(Node) ((Node).base != 'D') #define If_Edge_Exist(E) ((E).out_node != (uint64_t)-1) #define Visit(E) (E).length typedef struct { long long index; } RSet; inline void clear_RSet(RSet* set) { set->index = 0; } typedef struct { uint64_t ID; uint64_t weight; ///number of deletion end with current node uint64_t num_insertions; char base; Edge_alloc mismatch_edges; Edge_alloc deletion_edges; Edge_alloc insertion_edges; } Node; typedef struct { uint64_t* list; uint8_t* visit; uint64_t size; uint64_t length; uint64_t* iterative_buffer; uint8_t* iterative_buffer_visit; uint64_t iterative_i; } topo_Sorting_buffer; typedef struct { ///has a indivial start node 0 Node* list; topo_Sorting_buffer sort; uint64_t size; uint64_t length; uint64_t delete_length; } Node_alloc; typedef struct { uint64_t g_n_nodes; uint64_t g_n_edges; uint64_t g_next_nodeID; Node_alloc g_nodes; Queue node_q; char* seq; uint64_t seqID; uint64_t s_start_nodeID; uint64_t s_end_nodeID; } Graph; int add_and_check_bi_direction_edge(Graph* graph, Node* in_node, Node* out_node, uint64_t weight, uint64_t flag); void add_bi_direction_edge(Graph* graph, Node* in_node, Node* out_node, uint64_t weight, uint64_t flag); int remove_and_check_bi_direction_edge_from_nodes(Graph* graph, Node* in_node, Node* out_node); int remove_and_check_bi_direction_edge_from_edge(Graph* graph, Edge* e); inline int Pop_Node(Graph* DAGCon, Node** node) { long long nodeID = 0; int return_flag = pop_from_Queue(&(DAGCon->node_q), &nodeID); (*node) = &(G_Node(*DAGCon, nodeID)); return return_flag; } inline int Push_Node(Graph* DAGCon, Node** node) { push_to_Queue(&(DAGCon->node_q), (**node).ID); return 1; } inline int getInputNodes(RSet* set, Graph* graph, Node* node, Node** get_Node) { if(set->index >= (long long)Input_Edges(*node).length) { return 0; } ///skip all deleted edges while ( set->index < (long long)Input_Edges(*node).length && !(If_Edge_Exist(Input_Edges(*node).list[set->index])) ) { set->index++; } if( set->index < (long long)Input_Edges(*node).length && If_Edge_Exist(Input_Edges(*node).list[set->index]) ) { (*get_Node) = &(G_Node((*graph), Input_Edges(*node).list[set->index].in_node)); set->index++; return 1; } else { return 0; } } inline int getInputEdges(RSet* set, Graph* graph, Node* node, Edge** get_Edge) { if(set->index >= (long long)Input_Edges(*node).length) { return 0; } ///skip all deleted edges while ( set->index < (long long)Input_Edges(*node).length && !(If_Edge_Exist(Input_Edges(*node).list[set->index])) ) { set->index++; } if( set->index < (long long)Input_Edges(*node).length && If_Edge_Exist(Input_Edges(*node).list[set->index]) ) { (*get_Edge) = &(Input_Edges(*node).list[set->index]); set->index++; return 1; } else { return 0; } } inline int getOutputNodes(RSet* set, Graph* graph, Node* node, Node** get_Node) { if(set->index >= (long long)Output_Edges(*node).length) { return 0; } ///skip all deleted edges while ( set->index < (long long)Output_Edges(*node).length && !(If_Edge_Exist(Output_Edges(*node).list[set->index])) ) { set->index++; } if(set->index < (long long)Output_Edges(*node).length && If_Edge_Exist(Output_Edges(*node).list[set->index])) { (*get_Node) = &(G_Node((*graph), Output_Edges(*node).list[set->index].out_node)); set->index++; return 1; } else { return 0; } } inline int getOutputEdges(RSet* set, Graph* graph, Node* node, Edge** get_Edge) { if(set->index >= (long long)Output_Edges(*node).length) { return 0; } ///skip all deleted edges while ( set->index < (long long)Output_Edges(*node).length && !(If_Edge_Exist(Output_Edges(*node).list[set->index])) ) { set->index++; } if(set->index < (long long)Output_Edges(*node).length && If_Edge_Exist(Output_Edges(*node).list[set->index])) { (*get_Edge) = &(Output_Edges(*node).list[set->index]); set->index++; return 1; } else { return 0; } } inline void get_bi_direction_edges(Graph* DAGCon, Edge* edge, Edge** e_forward, Edge** e_backward) { long long in_node = edge->in_node; long long out_node = edge->out_node; if( edge->self_edge_ID < Output_Edges(G_Node(*DAGCon, in_node)).length && (long long)Output_Edges(G_Node(*DAGCon, in_node)).list[edge->self_edge_ID].in_node == in_node && (long long)Output_Edges(G_Node(*DAGCon, in_node)).list[edge->self_edge_ID].out_node == out_node ) { (*e_forward) = &(Output_Edges(G_Node(*DAGCon, in_node)).list[edge->self_edge_ID]); (*e_backward) = &(Input_Edges(G_Node(*DAGCon, out_node)).list[edge->reverse_edge_ID]); } else { (*e_forward) = &(Output_Edges(G_Node(*DAGCon, in_node)).list[edge->reverse_edge_ID]); (*e_backward) = &(Input_Edges(G_Node(*DAGCon, out_node)).list[edge->self_edge_ID]); } } inline long long get_bi_Edge(Graph* DAGCon, Node* inNode, Node* outNode, Edge** e_forward, Edge** e_backward) { Edge* e; RSet iter; clear_RSet(&iter); if(If_Node_Exist(*inNode) && If_Node_Exist(*outNode)) { //find in-edge of outNode while(getInputEdges(&iter, DAGCon, outNode, &e)) { if(e->in_node == inNode->ID) { get_bi_direction_edges(DAGCon, e, e_forward, e_backward); return 1; } } } return 0; } inline long long get_Edge_Weight(Graph* DAGCon, Node* inNode, Node* outNode) { Edge* e_forward = NULL; Edge* e_backward = NULL; get_bi_Edge(DAGCon, inNode, outNode, &e_forward, &e_backward); return e_forward->weight; } 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, uint64_t length); void init_Node_alloc(Node_alloc* list); void destory_Node_alloc(Node_alloc* list); void clear_Node_alloc(Node_alloc* list); uint64_t append_Node_alloc(Node_alloc* list, char base); 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, window_list *cigar_idx, window_list_alloc *cigar_s, 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); uint64_t inline add_Node_Graph(Graph* g, char base) { return append_Node_alloc(&g->g_nodes, base); } inline Node* add_Node_DAGCon(Graph* g, char base) { return &(G_Node(*g, append_Node_alloc(&g->g_nodes, base))); } ///to delete a node ///1. set the corresponding base to be 'D' ///2. remove all related edges ///2. clear all related edges ///3. g_nodes.delete_length++, please do not substract g_nodes.length uint64_t inline delete_Node_DAGCon(Graph* g, Node* node) { g->g_nodes.delete_length++; g->g_nodes.list[(*node).ID].base = 'D'; g->g_nodes.list[(*node).ID].num_insertions = (uint64_t)-1; g->g_nodes.list[(*node).ID].weight = (uint64_t)-1; RSet iter; Edge* e; clear_RSet(&iter); while (getOutputEdges(&iter, g, node, &e)) { remove_and_check_bi_direction_edge_from_edge(g, e); } clear_RSet(&iter); while (getInputEdges(&iter, g, node, &e)) { remove_and_check_bi_direction_edge_from_edge(g, e); } clear_Edge_alloc(&(g->g_nodes.list[(*node).ID].insertion_edges)); clear_Edge_alloc(&(g->g_nodes.list[(*node).ID].mismatch_edges)); clear_Edge_alloc(&(g->g_nodes.list[(*node).ID].deletion_edges)); return 1; } ///just for mimatch edges 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 < (long long)edge->length; i++) { nodeID = edge->list[i].out_node; if(g->g_nodes.list[nodeID].base == base) { edge->list[i].weight++; ///if last operation is insertion if (last_operation == 2) { edge->list[i].num_insertions++; } break; } } ///there are no such edge if (i == (long long)edge->length) { nodeID = add_Node_Graph(g, base); ///the length of match edge is 0, while the length of mismatch edge is 1 append_Edge_alloc(edge, in_node, nodeID, 1, 1); ///if last operation is insertion if (last_operation == 2) { edge->list[edge->length - 1].num_insertions++; } ///add the mismatch_edges of new node to the backbone append_Edge_alloc(&(g->g_nodes.list[nodeID].mismatch_edges), nodeID, in_node + 1, 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 < (long long)edge->length; i++) { nodeID = edge->list[i].out_node; if(nodeID == nextNodeID) { edge->list[i].weight++; break; } } ///there are no such edge if (i == (long long)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) { long long i; for (i = 0; i < deletion_length; i++) { add_single_deletionEdge_weight(g, alignNodeID + i, alignNodeID + i + 1, 0); } } 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 < (long long)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 < (long long)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; } } /****************************may have bugs********************************/ return return_edgeID; /****************************may have bugs********************************/ } inline int create_insertion_Edges(Graph* g, long long alignNodeID, uint64_t edge_length, char* bases) { long long i = 0; long long nodeID; ///should link back to the intial node ///long long backboneID = alignNodeID + 1; long long backboneID = alignNodeID; if (edge_length < 1) { return -1; } nodeID = add_Node_Graph(g, bases[0]); ///add the new node to alignNodeID by insertion_edges append_Edge_alloc(&(g->g_nodes.list[alignNodeID].insertion_edges), alignNodeID, nodeID, 1, edge_length); alignNodeID = nodeID; for (i = 1; i < (long long)edge_length; i++) { nodeID = add_Node_Graph(g, bases[i]); ///add the new node to alignNodeID by insertion_edges 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); return 1; } inline void extract_path(Graph* backbone, int debug_node_in_backbone, int path_i, char* pre) { int step = G_Node(*backbone, debug_node_in_backbone).insertion_edges.list[path_i].length; int string_i = 0, preNode = 0, j = 0; if(step != 0) { string_i = 0; preNode = G_Node(*backbone, debug_node_in_backbone).insertion_edges.list[path_i].out_node; for (j = 0; j < step; j++) { pre[string_i++] = G_Node(*backbone, preNode).base; preNode = G_Node(*backbone, preNode).insertion_edges.list[0].out_node; } } pre[string_i] = '\0'; } inline int get_insertion_Edges_new(Graph* backbone, int debug_node_in_backbone, uint64_t edge_length, char* bases) { int path_i, j, step, preNode; for (path_i = 0; path_i < (long long)G_Node(*backbone, debug_node_in_backbone).insertion_edges.length; path_i++) { step = G_Node(*backbone, debug_node_in_backbone).insertion_edges.list[path_i].length; if(step != (long long)edge_length) { continue; } if(step != 0) { preNode = G_Node(*backbone, debug_node_in_backbone).insertion_edges.list[path_i].out_node; for (j = 0; j < step; j++) { if(G_Node(*backbone, preNode).base != bases[j]) { break; } preNode = G_Node(*backbone, preNode).insertion_edges.list[0].out_node; } if(j == step) { return path_i; } } } return -1; } inline void add_insertionEdge_weight(Graph* g, long long alignNodeID, char* insert, long long insert_length) { 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 ///there is no such edge { nodeID = add_Node_Graph(g, insert[0]); append_Edge_alloc(edge, alignNodeID, nodeID, 1, 1); ///add the new node to alignNodeID by insertion_edges //should link to the initial node, instead of the next node of the initial node ///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 { ///edgeID = get_insertion_Edges(g, edge, insert_length, insert); edgeID = get_insertion_Edges_new(g, alignNodeID, insert_length, insert); if (edgeID != -1) { ///just one outdegree edge->list[edgeID].weight++; } else { create_insertion_Edges(g, alignNodeID, insert_length, insert); } } } #endif