final DP-based assembly

POA.cpp

@@ -15,11 +15,13 @@ void init_Edge_alloc(Edge_alloc* list)
     {
         list->size = INIT_EDGE_SIZE;
         list->length = 0;
+        list->delete_length = 0;
         list->list = (Edge*)malloc(sizeof(Edge)*list->size);
     }
     else
     {
         list->length = 0;
+        list->delete_length = 0;
     }
     
 }
@@ -27,6 +29,7 @@ void init_Edge_alloc(Edge_alloc* list)
 void clear_Edge_alloc(Edge_alloc* list)
 {
     list->length = 0;
+    list->delete_length = 0;
 }
 
 void destory_Edge_alloc(Edge_alloc* list)
@@ -47,17 +50,138 @@ void append_Edge_alloc(Edge_alloc* list,  uint64_t in_node, uint64_t out_node, u
     list->list[list->length].weight = weight;
     list->list[list->length].length = length;
     list->list[list->length].num_insertions = 0;
+    list->list[list->length].self_edge_ID = list->length;
 
     list->length++;
 }
 
 
+int add_and_check_bi_direction_edge(Graph* graph, Node* in_node, Node* out_node, uint64_t weight, uint64_t flag)
+{
+    Edge* e_forward;
+    Edge* e_backward;
+
+    //if there are no edge from in_node to out_node
+    if(!get_bi_Edge(graph, in_node, out_node, &e_forward, &e_backward))
+    {
+        append_Edge_alloc(&(Output_Edges((*in_node))), (*in_node).ID, (*out_node).ID, weight, flag);
+        append_Edge_alloc(&(Input_Edges((*out_node))), (*in_node).ID, (*out_node).ID, weight, flag);
+
+        Output_Edges((*in_node)).list[Output_Edges((*in_node)).length - 1].reverse_edge_ID 
+        = Input_Edges((*out_node)).length - 1;
+
+        Input_Edges((*out_node)).list[Input_Edges((*out_node)).length - 1].reverse_edge_ID
+        = Output_Edges((*in_node)).length - 1;
+
+        return 1;
+    }
+    else//if there is an edge from in_node to out_node, do nothing
+    {
+        return 0;
+    }
+}
+
+
+void add_bi_direction_edge(Graph* graph, Node* in_node, Node* out_node, uint64_t weight, uint64_t flag)
+{
+
+    append_Edge_alloc(&(Output_Edges((*in_node))), (*in_node).ID, (*out_node).ID, weight, flag);
+    append_Edge_alloc(&(Input_Edges((*out_node))), (*in_node).ID, (*out_node).ID, weight, flag);
+
+    Output_Edges((*in_node)).list[Output_Edges((*in_node)).length - 1].reverse_edge_ID 
+    = Input_Edges((*out_node)).length - 1;
+
+    Input_Edges((*out_node)).list[Input_Edges((*out_node)).length - 1].reverse_edge_ID
+    = Output_Edges((*in_node)).length - 1;
+}
+
+
+
+int remove_and_check_bi_direction_edge_from_nodes(Graph* graph, Node* in_node, Node* out_node)
+{
+    Edge* e_forward;
+    Edge* e_backward;
+
+    //if there are no edge from in_node to out_node
+    //1. remove these two edges
+    //2. increase the edge_list.delete_length in both in_node and out_node
+    if(get_bi_Edge(graph, in_node, out_node, &e_forward, &e_backward))
+    {
+        e_forward->in_node = (uint64_t)-1;
+        e_forward->out_node = (uint64_t)-1;
+        e_forward->weight = (uint64_t)-1;
+        e_forward->length = (uint64_t)-1;
+        e_forward->num_insertions = (uint64_t)-1;
+        e_forward->self_edge_ID = (uint64_t)-1;
+        e_forward->reverse_edge_ID = (uint64_t)-1;
+
+
+        e_backward->in_node = (uint64_t)-1;
+        e_backward->out_node = (uint64_t)-1;
+        e_backward->weight = (uint64_t)-1;
+        e_backward->length = (uint64_t)-1;
+        e_backward->num_insertions = (uint64_t)-1;
+        e_backward->self_edge_ID = (uint64_t)-1;
+        e_backward->reverse_edge_ID = (uint64_t)-1;
+
+        Output_Edges(*in_node).delete_length++;
+        Input_Edges((*out_node)).delete_length++;
+
+        return 1;
+    }
+    else//if there is an edge from in_node to out_node, do nothing
+    {
+        return 0;
+    }
+}
+
+
+
+int remove_and_check_bi_direction_edge_from_edge(Graph* graph, Edge* e)
+{
+    Edge* e_forward;
+    Edge* e_backward;
+
+    if(If_Edge_Exist(*e))
+    {
+        get_bi_direction_edges(graph, e, &e_forward, &e_backward);
+        Output_Edges(G_Node(*graph, e_forward->in_node)).delete_length++;
+        Input_Edges(G_Node(*graph, e_forward->out_node)).delete_length++;
+
+        e_forward->in_node = (uint64_t)-1;
+        e_forward->out_node = (uint64_t)-1;
+        e_forward->weight = (uint64_t)-1;
+        e_forward->length = (uint64_t)-1;
+        e_forward->num_insertions = (uint64_t)-1;
+        e_forward->self_edge_ID = (uint64_t)-1;
+        e_forward->reverse_edge_ID = (uint64_t)-1;
+
+
+        e_backward->in_node = (uint64_t)-1;
+        e_backward->out_node = (uint64_t)-1;
+        e_backward->weight = (uint64_t)-1;
+        e_backward->length = (uint64_t)-1;
+        e_backward->num_insertions = (uint64_t)-1;
+        e_backward->self_edge_ID = (uint64_t)-1;
+        e_backward->reverse_edge_ID = (uint64_t)-1;
+
+        return 1;
+    }
+    else
+    {
+        return 0;
+    }
+    
+    
+}
+
 
 
 void init_Node_alloc(Node_alloc* list)
 {
     list->size = INIT_NODE_SIZE;
     list->length = 0;
+    list->delete_length = 0;
     list->list = (Node*)malloc(sizeof(Node)*list->size);
     list->sort.size = 0;
     list->sort.list = NULL;
@@ -105,6 +229,7 @@ void clear_Node_alloc(Node_alloc* list)
     }
 
     list->length = 0;
+    list->delete_length = 0;
 }
 
 
@@ -162,11 +287,14 @@ void init_Graph(Graph* g)
     g->s_start_nodeID = 0;
     g->seq = NULL;
     g->seqID = (uint64_t)-1;
+
+    init_Queue(&(g->node_q));
 }
 
 void destory_Graph(Graph* g)
 {
     destory_Node_alloc(&g->g_nodes);
+    destory_Queue(&(g->node_q));
 }
 
 void clear_Graph(Graph* g)
@@ -180,6 +308,8 @@ void clear_Graph(Graph* g)
     g->s_start_nodeID = 0;
     g->seq = NULL;
     g->seqID = (uint64_t)-1;
+
+    clear_Queue(&(g->node_q));
 }
 
 
@@ -237,6 +367,282 @@ void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length,
 
 
 
+inline void add_insertionEdge_weight_print(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);
+
+    /******************************for homopolymer*************************/
+    long long i = 0;
+    char hom;
+    if (insert_length > 0)
+    {
+        hom = insert[0];
+    }
+
+    for (i = 0; i < insert_length; i++)
+    {
+        if(insert[i] != hom)
+        {
+            break;
+        }
+    }
+
+    fprintf(stderr, "###insert_length: %d\n", insert_length);
+
+    ///if it is a homopolymer
+    if(i == insert_length)
+    {
+        ///single base
+        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);
+        }
+
+        ///multiple bases
+        for (i = 1; i < insert_length; i++)
+        {
+            edgeID = get_insertion_Edges(g, edge, i + 1, insert);
+            if (edgeID != -1)
+            {
+                ///这条路均只有一个出度
+                edge->list[edgeID].weight++;
+            }
+            else
+            {
+                create_insertion_Edges(g, alignNodeID, i + 1, insert);
+            }
+        }
+
+        return;
+    }
+    /******************************for homopolymer*************************/
+
+    fprintf(stderr, "###not homopolymer: %d\n", insert_length);
+
+    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]);
+        fprintf(stderr, "edgeID[0]: %d, length: %d\n", edgeID, edge->list[edgeID].length);
+        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);
+        }
+
+        fprintf(stderr, "edge->length: %d\n", edge->length);
+        
+        /*******************第0个字符********************* */
+
+        /*******************第1个字符********************* */
+        if (insert[1] != insert[0])
+        {
+            edgeID = getEdge(g, edge, 1, insert[1]);
+            fprintf(stderr, "edgeID[1]: %d, length: %d\n", edgeID, edge->list[edgeID].length);
+            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);
+            }
+
+            fprintf(stderr, "edge->length: %d\n", edge->length);
+        }
+        /*******************第1个字符********************* */
+
+        /**********************两个字符******************* */
+        
+        edgeID = get_insertion_Edges(g, edge, 2, insert);
+        fprintf(stderr, "edgeID[len2]: %d, length: %d\n", edgeID, edge->list[edgeID].length);
+        if (edgeID != -1)
+        {
+            ///这条路均只有一个出度
+            edge->list[edgeID].weight++;
+        }
+        else
+        {
+            create_insertion_Edges(g, alignNodeID, insert_length, insert);
+        }
+
+        fprintf(stderr, "edge->length: %d\n", edge->length);
+
+        for (i = 0; i < edge->length; i++)
+        {
+            fprintf(stderr, "edgeID[%d].length: %d\n", i, edge->list[i].length);
+        }
+        
+        
+        /**********************两个字符******************* */
+    }
+    else if (insert_length > 2)
+    {
+        ////fprintf(stderr, "too long insertion\n");
+        /*************************大于2个字符************************** */
+
+        edgeID = get_insertion_Edges(g, edge, insert_length, insert);
+        if (edgeID != -1)
+        {
+            ///这条路均只有一个出度
+            edge->list[edgeID].weight++;
+        }
+        else
+        {
+            create_insertion_Edges(g, alignNodeID, insert_length, insert);
+        }
+    }
+    
+    
+    
+}
+
+
+
+void addmatchedSeqToGraph_print(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;
+    int last_operation = -1;
+
+    fprintf(stderr,"*******\n");
+    ///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 (cigar_i < cigar->length)
+    {
+        operation = cigar->C_C[cigar_i];
+        operationLen = cigar->C_L[cigar_i];
+
+        // fprintf(stderr, "operation: %d, operationLen: %d\n", 
+        // operation, operationLen);
+
+        ///这种情况代表匹配和mismatch
+        if (operation == 0 || operation == 1)
+        {
+            
+            for (i = 0; i < operationLen; i++)
+            {
+                //backbone->g_nodes.list[currentNodeID].weight++;
+                ///前面是插入，后面有可能是误配，也有可能是匹配
+                add_mismatchEdge_weight(backbone, currentNodeID, y_string[y_i], last_operation);    
+                x_i++;
+                y_i++;
+                currentNodeID++;
+            }
+        }///insertion
+        else if (operation == 2)
+        {
+            ///cigar的起始和结尾不可能是2，所以这里-1没问题
+            ///if (operationLen <= CORRECT_INDEL_LENGTH)
+            {
+                add_insertionEdge_weight_print(backbone, currentNodeID, y_string + y_i, operationLen);
+                backbone->g_nodes.list[currentNodeID].num_insertions++;
+            }
+
+            ///fprintf(stderr, "y_string: %.*s\n",  operationLen, y_string+y_i);
+            y_i += operationLen;
+        }
+        else if (operation == 3)
+        {
+            ///3是y缺字符（x多字符），也就是backbone多字符
+            ///这个相当于在backbone对应字符处变成了‘——’
+            ///因此可以用mismatch类似的方法处理
+            ///if (operationLen <= CORRECT_INDEL_LENGTH)
+            {
+                ///add_deletion_to_backbone(backbone, &currentNodeID, operationLen);
+                ///在编辑距离中，前面是个插入，后面是个删除，这种情况是不存在的
+                ///为了保险要不还给他加上吧
+                ///先不加
+                add_deletionEdge_weight(backbone, currentNodeID, operationLen);
+            }
+
+            
+            currentNodeID += operationLen;
+            x_i += operationLen;
+        }
+
+        last_operation = operation;
+        
+        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);
+    }
+    **/
+    
+}
 
 
 
@@ -311,28 +717,6 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri
         
         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);
-    }
-    **/
-    
 }