POA

Assembly.cpp

@@ -6,6 +6,7 @@
 #include "CommandLines.h"
 #include "kmer.h"
 #include "Hash_Table.h"
+#include "POA.h"
 
 Total_Count_Table TCB;
 Total_Pos_Table PCB;
@@ -679,7 +680,7 @@ void* Overlap_calculate_heap_merge(void* arg)
 {
 
     int thr_ID = *((int*)arg);
-
+    uint64_t POA_i;
     long long i = 0;
     int avalible_k = 0;
 
@@ -695,7 +696,9 @@ void* Overlap_calculate_heap_merge(void* arg)
 
     Candidates_list l;
     //Candidates_list debug_l;
+    Graph POA_Graph;
 
+    init_Graph(&POA_Graph);
     init_Candidates_list(&l);
     //init_Candidates_list(&debug_l);
 
@@ -711,6 +714,13 @@ void* Overlap_calculate_heap_merge(void* arg)
 
     for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num)
     {
+        /**
+        if (thr_ID == 0 && i % 1000 == 0)
+        {
+            fprintf(stderr, "i: %llu\n", i);
+        }
+        **/
+        
         
         clear_Heap(&heap);
         clear_Candidates_list(&l);
@@ -756,8 +766,30 @@ void* Overlap_calculate_heap_merge(void* arg)
 
 
 
+
+        
+
+
+
         ///reverse complement strand
         reverse_complement(g_read.seq, g_read.length);
+        /**
+        UC_Read rc_read;
+        init_UC_Read(&rc_read);
+        recover_UC_Read_RC(&rc_read, &R_INF, i);
+        uint64_t j = 0;
+        for (j = 0; j < g_read.length; j++)
+        {
+            if (g_read.seq[j] != rc_read.seq[j])
+            {
+                fprintf(stderr, "j error: %llu, i: %llu\n", j, i);
+                fprintf(stderr, "g_read.seq[j]: %c\n", g_read.seq[j]);
+                fprintf(stderr, "rc_read.seq[j]: %c\n", rc_read.seq[j]);
+            }
+            
+        }
+        destory_UC_Read(&rc_read);
+        **/
         init_HPC_seq(&HPC_read, g_read.seq, g_read.length);
         init_Hash_code(&k_code);
         avalible_k = 0;
@@ -811,6 +843,37 @@ void* Overlap_calculate_heap_merge(void* arg)
 
         ///debug_merge_result(&l, &debug_l);
 
+        clear_Graph(&POA_Graph);
+
+
+        Perform_POA(&POA_Graph, &overlap_list, &R_INF, &g_read);
+
+        ///fprintf(stderr, "i: %u\n", i);
+
+
+        /**
+        POA_i = 0;
+
+        if (overlap_list.length > 0)
+        {
+            
+        }
+        
+
+        for (POA_i = 1; POA_i < overlap_list.length; POA_i++)
+        {
+            if(overlap_list.list[POA_i].x_pos_strand == 1)
+            {
+                fprintf(stderr, "sbsbsbsbs\n");
+            }
+
+        }
+        **/
+        
+
+
+
+
     }
     
     destory_Candidates_list(&l);
@@ -819,6 +882,9 @@ void* Overlap_calculate_heap_merge(void* arg)
 
     destory_Heap(&heap);
     destory_k_mer_pos_list_alloc(&array_list);
+
+    destory_Graph(&POA_Graph);
+    destory_UC_Read(&g_read);
 }
 
 
@@ -832,6 +898,8 @@ void Overlap_calculate_multipe_thr()
 {
     double start_time = Get_T();
 
+    fprintf(stdout, "R_INF.total_reads: %llu\n", R_INF.total_reads);
+
 
     fprintf(stdout, "Begin Overlap Calculate ...... \n");
 

CommandLines.cpp

@@ -9,7 +9,8 @@ char* read_file_name = NULL;
 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 = 9;
+int k_mer_min_freq = 3;
 int k_mer_max_freq = 66;
 int load_index_from_disk = 0;
 int write_index_to_disk = 0;

Hash_Table.cpp

@@ -256,6 +256,7 @@ void append_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp
 
 
 
+
     if(tmp->x_pos_s <= tmp->y_pos_s)
     {
         tmp->y_pos_s = tmp->y_pos_s - tmp->x_pos_s;
@@ -283,18 +284,36 @@ void append_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp
     
 
 
-    
-    
-    
-    list->list[list->length].x_id = tmp->x_id;
-    list->list[list->length].x_pos_e = tmp->x_pos_e;
-    list->list[list->length].x_pos_s = tmp->x_pos_s;
-    list->list[list->length].x_pos_strand = tmp->x_pos_strand;
 
-    list->list[list->length].y_id = tmp->y_id;
-    list->list[list->length].y_pos_e = tmp->y_pos_e;
-    list->list[list->length].y_pos_s = tmp->y_pos_s;
-    list->list[list->length].y_pos_strand = tmp->y_pos_strand;
+    ///y的区间方向永远是0
+    ///x有可能是1
+    ///如果x是1，那我们把它逆过来
+    if (tmp->x_pos_strand == 1)
+    {
+       list->list[list->length].x_id = tmp->x_id;
+        list->list[list->length].x_pos_e = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_s - 1;
+        list->list[list->length].x_pos_s = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_e - 1;
+        list->list[list->length].x_pos_strand = 0;
+
+        list->list[list->length].y_id = tmp->y_id;
+        list->list[list->length].y_pos_e = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_s - 1;
+        list->list[list->length].y_pos_s = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_e - 1;
+        list->list[list->length].y_pos_strand = 1;
+    }
+    else
+    {
+        list->list[list->length].x_id = tmp->x_id;
+        list->list[list->length].x_pos_e = tmp->x_pos_e;
+        list->list[list->length].x_pos_s = tmp->x_pos_s;
+        list->list[list->length].x_pos_strand = tmp->x_pos_strand;
+
+        list->list[list->length].y_id = tmp->y_id;
+        list->list[list->length].y_pos_e = tmp->y_pos_e;
+        list->list[list->length].y_pos_s = tmp->y_pos_s;
+        list->list[list->length].y_pos_strand = tmp->y_pos_strand;
+    }
+    
+    
 
     list->list[list->length].shared_seed = tmp->shared_seed;
 
@@ -536,7 +555,13 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF)
             break;
         }
 
-        append_overlap_region_alloc(overlap_list, &tmp_region, R_INF);
+        ///自己和自己重叠的要排除
+        if (tmp_region.x_id != tmp_region.y_id)
+        {
+            append_overlap_region_alloc(overlap_list, &tmp_region, R_INF);
+        }
+        
+        
     }
 
 

Hash_Table.h

@@ -398,9 +398,7 @@ void clear_Heap(HeapSq* HBT);
 void init_overlap_region_alloc(overlap_region_alloc* list);
 void clear_overlap_region_alloc(overlap_region_alloc* list);
 void destory_overlap_region_alloc(overlap_region_alloc* list);
-void append_overlap_region_alloc(overlap_region_alloc* list,
-uint64_t x_id, uint64_t x_pos_s, uint64_t x_pos_e, uint64_t x_pos_strand,
-uint64_t y_id, uint64_t y_pos_s, uint64_t y_pos_e, uint64_t y_pos_strand);
+void append_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, All_reads* R_INF);
 void calculate_overlap_region(Candidates_list* candidates, overlap_region_alloc* overlap_list, 
 uint64_t readID, uint64_t readLength, All_reads* R_INF);
 

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
+SOURCES = main.cpp CommandLines.cpp Process_Read.cpp Assembly.cpp kmer.cpp Hash_Table.cpp POA.cpp
 OBJECTS = $(SOURCES:.c=.o)
 EXECUTABLE = ccs_assembly
 

POA.cpp

@@ -0,0 +1,464 @@
+#include "POA.h"
+#include <stdlib.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)
+{
+    if (list->list == NULL)
+    {
+        list->size = INIT_EDGE_SIZE;
+        list->length = 0;
+        list->list = (Edge*)malloc(sizeof(Edge)*list->size);
+    }
+    else
+    {
+        list->length = 0;
+    }
+    
+}
+
+void clear_Edge_alloc(Edge_alloc* list)
+{
+    list->length = 0;
+}
+
+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)
+{
+    if (list->length + 1 > list->size)
+    {
+        list->size = list->size + INCREASE_EDGE_SIZE;
+        list->list = (Edge*)realloc(list->list, sizeof(Edge)*list->size);
+    }
+
+    list->list[list->length].in_node = in_node;
+    list->list[list->length].out_node = out_node;
+    list->list[list->length].weight = weight;
+
+    list->length++;
+}
+
+
+
+
+
+
+
+
+void init_Node_alloc(Node_alloc* list)
+{
+    list->size = INIT_NODE_SIZE;
+    list->length = 0;
+    list->list = (Node*)malloc(sizeof(Node)*list->size);
+    list->sort.size = 0;
+    list->sort.list = NULL;
+    list->sort.visit = NULL;
+
+    list->sort.iterative_buffer = NULL;
+    list->sort.iterative_buffer_visit = NULL;
+
+
+    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;
+    }
+    
+}
+
+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);
+    }
+    free(list->list);
+    free(list->sort.list);
+    free(list->sort.visit);
+    free(list->sort.iterative_buffer);
+    free(list->sort.iterative_buffer_visit);
+    ///free(list->topo_order);
+}
+
+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);
+    }
+    list->length = 0;
+}
+
+
+uint64_t append_Node_alloc(Node_alloc* list, char base)
+{
+   
+    if (list->length + 1 > list->size)
+    {
+        long long i = list->size;
+
+        ///list->topo_order这里用不到，所以不用先分配空间
+        ///但是还是一起分配了吧，免得麻烦
+        list->size = list->size * 2;
+        list->list = (Node*)realloc(list->list, sizeof(Node)*list->size);
+        ///list->topo_order = (uint64_t*)realloc(list->topo_order, sizeof(uint64_t)*list->size);
+        
+        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[list->length].ID = list->length;
+    list->list[list->length].base = base;
+    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->length++;
+
+    return list->length - 1;
+}
+
+
+
+
+
+
+
+
+
+
+
+void init_Graph(Graph* g)
+{
+    init_Node_alloc(&g->g_nodes);
+    g->g_n_edges = 0;
+    g->g_n_nodes = 0;
+    g->g_next_nodeID = 0;
+    g->s_end_nodeID = 0;
+    g->s_start_nodeID = 0;
+    g->seq = NULL;
+    g->seqID = (uint64_t)-1;
+}
+
+void destory_Graph(Graph* g)
+{
+    destory_Node_alloc(&g->g_nodes);
+}
+
+void clear_Graph(Graph* g)
+{
+    clear_Node_alloc(&g->g_nodes);
+
+    g->g_n_edges = 0;
+    g->g_n_nodes = 0;
+    g->g_next_nodeID = 0;
+    g->s_end_nodeID = 0;
+    g->s_start_nodeID = 0;
+    g->seq = NULL;
+    g->seqID = (uint64_t)-1;
+}
+
+
+
+
+
+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;
+    firstID = -1;
+    lastID = -1;
+
+    if(g_read_length == 0)
+        return;
+    
+    for (i = 0; i < g_read_length; i++)
+    {
+        nodeID = add_Node_Graph(g, g_read_seq[i]);
+
+        ////fprintf(stderr, "nodeID: %llu\n", nodeID);
+        
+        if (firstID == -1)
+        {
+            firstID = nodeID;
+        }
+        if (lastID != -1)
+        {
+            add_Edge_Graph(g, lastID, nodeID, 1);
+        }
+
+        lastID = nodeID; 
+    }
+
+    *startID = firstID;
+    *endID = lastID;
+    
+}
+
+
+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;
+    
+}

POA.h

@@ -0,0 +1,104 @@
+#ifndef __POA_PARSER__
+#define __POA_PARSER__
+#include <stdint.h>
+#include "Hash_Table.h"
+#include "Process_Read.h"
+
+/**
+ 1. 单个节点信息
+    (1) ID
+    (2) base
+    (3) 入边信息
+    (4) 出边信息
+    (5) 比对到什么节点
+ 2. 各个节点信息，用数组下标组织，数组下标就是节点ID; 还要存拓扑排序后的下标和节点ID的对应关系
+ 3. 边
+    (1) 边的起始
+    (2) 边的结束节点
+    (3) 过这条边的序列的label，也就是名称
+ 4. 各个序列信息 
+    (1) 这个序列本身 
+    (2) 这个序列的name或者ID 
+    (3) 这个序列的在图中对应的起始和结束节点ID
+ 5. 有两个回溯矩阵，一个是graph的，一个是seq的
+ **/
+
+typedef struct
+{
+    uint64_t in_node;
+    uint64_t out_node;
+    uint64_t weight;
+} Edge;
+
+typedef struct
+{
+    Edge* list;
+    uint64_t size;
+    uint64_t length;
+} Edge_alloc;
+
+typedef struct
+{
+    uint64_t ID;
+    char base;
+    Edge_alloc income_edges;
+    Edge_alloc outcome_edges;
+    Edge_alloc alignedTo_Nodes;
+} 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
+{
+    Node* list;
+    topo_Sorting_buffer sort;
+    uint64_t size;
+    uint64_t length;
+} Node_alloc;
+
+typedef struct
+{
+    uint64_t g_n_nodes;
+    uint64_t g_n_edges;
+    uint64_t g_next_nodeID;
+    Node_alloc g_nodes;
+
+
+
+    char* seq;
+    uint64_t seqID;
+    uint64_t s_start_nodeID;
+    uint64_t s_end_nodeID;
+} Graph;
+
+
+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 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 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);
+
+
+#endif

Process_Read.cpp

@@ -204,6 +204,12 @@ void malloc_All_reads(All_reads* r)
 	
 }
 
+void destory_UC_Read(UC_Read* r)
+{
+	free(r->seq);
+}
+
+
 void init_UC_Read(UC_Read* r)
 {
 	r->length = 0;
@@ -219,6 +225,11 @@ void init_UC_Read(UC_Read* r)
 			bit_t_seq_table[i][1] = s_H[((i >> 4)&(uint64_t)3)];
 			bit_t_seq_table[i][2] = s_H[((i >> 2)&(uint64_t)3)];
 			bit_t_seq_table[i][3] = s_H[(i&(uint64_t)3)];
+
+			bit_t_seq_table_rc[i][0] = RC_CHAR(bit_t_seq_table[i][3]);
+			bit_t_seq_table_rc[i][1] = RC_CHAR(bit_t_seq_table[i][2]);
+			bit_t_seq_table_rc[i][2] = RC_CHAR(bit_t_seq_table[i][1]);
+			bit_t_seq_table_rc[i][3] = RC_CHAR(bit_t_seq_table[i][0]);
 		}
 		
 	}
@@ -256,6 +267,46 @@ void recover_UC_Read(UC_Read* r, All_reads* R_INF, uint64_t ID)
 		
 }
 
+void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID)
+{
+	r->length = Get_READ_LENGTH((*R_INF), ID);
+	uint8_t* src = Get_READ((*R_INF), ID);
+
+	if (r->length + 4 > r->size)
+	{
+		r->size = r->length + 4;
+		r->seq = (char*)realloc(r->seq,sizeof(char)*(r->size));
+	}
+
+	long long last_chr = r->length % 4;
+	long long i = r->length / 4 - 1 + (last_chr != 0);
+	long long index = 0;
+
+	if(last_chr!=0)
+	{
+		memcpy(r->seq + index, bit_t_seq_table_rc[src[i]] + 4 - last_chr, last_chr);
+		index = last_chr;
+		i--;
+	}
+
+	while (i >= 0)
+	{
+		memcpy(r->seq + index, bit_t_seq_table_rc[src[i]], 4);
+		i--;
+		index = index + 4;
+	}
+
+
+	if (R_INF->N_site[ID])
+	{
+		for (i = 1; i <= R_INF->N_site[ID][0]; i++)
+		{
+			r->seq[r->length - R_INF->N_site[ID][i] - 1] = 'N';
+		}
+	}
+		
+}
+
 
 
 #define COMPRESS_BASE {c = seq_nt6_table[src[i]];\

Process_Read.h

@@ -45,6 +45,7 @@ static uint8_t seq_nt6_table[256] = {
 };
 
 static char bit_t_seq_table[256][4] = {0};
+static char bit_t_seq_table_rc[256][4] = {0};
 static char s_H[4] = {'A', 'C', 'G', 'T'};
 static char rc_Table[4] = {'T', 'G', 'C', 'A'};
 
@@ -109,6 +110,8 @@ int get_reads_mul_thread(R_buffer_block* curr_sub_block);
 void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ);
 void init_UC_Read(UC_Read* r);
 void recover_UC_Read(UC_Read* r, All_reads* R_INF, uint64_t ID);
+void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID);
+void destory_UC_Read(UC_Read* r);
 void reverse_complement(char* pattern, uint64_t length);
 void write_All_reads(All_reads* r, char* read_file_name);
 int load_All_reads(All_reads* r, char* read_file_name);