diff --git a/.vscode/ipch/2549b0c830b4882d/POA.ipch b/.vscode/ipch/2549b0c830b4882d/POA.ipch new file mode 100644 index 0000000..a84cd9b Binary files /dev/null and b/.vscode/ipch/2549b0c830b4882d/POA.ipch differ diff --git a/.vscode/ipch/2549b0c830b4882d/mmap_address.bin b/.vscode/ipch/2549b0c830b4882d/mmap_address.bin new file mode 100644 index 0000000..862b842 Binary files /dev/null and b/.vscode/ipch/2549b0c830b4882d/mmap_address.bin differ diff --git a/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch b/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch index 11f3e3e..a33202d 100644 Binary files a/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch and b/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch differ diff --git a/.vscode/ipch/856c69690c514632/KMER.ipch b/.vscode/ipch/856c69690c514632/KMER.ipch index 0f60803..23c0c45 100644 Binary files a/.vscode/ipch/856c69690c514632/KMER.ipch and b/.vscode/ipch/856c69690c514632/KMER.ipch differ diff --git a/.vscode/ipch/a64815cd9ef3d0a4/mmap_address.bin b/.vscode/ipch/a64815cd9ef3d0a4/mmap_address.bin new file mode 100644 index 0000000..862b842 Binary files /dev/null and b/.vscode/ipch/a64815cd9ef3d0a4/mmap_address.bin differ diff --git a/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch b/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch index c01351c..b11e08a 100644 Binary files a/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch and b/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch differ diff --git a/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch b/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch index 21270d7..fcf8a67 100644 Binary files a/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch and b/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch differ diff --git a/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch b/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch index 8a465a3..bb86d9c 100644 Binary files a/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch and b/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch differ diff --git a/Assembly.cpp b/Assembly.cpp index 2e607ea..128ac7b 100644 --- a/Assembly.cpp +++ b/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"); diff --git a/CommandLines.cpp b/CommandLines.cpp index dc70ddc..0b2e7a2 100644 --- a/CommandLines.cpp +++ b/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; diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 3161939..9ddae36 100644 --- a/Hash_Table.cpp +++ b/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); + } + + } diff --git a/Hash_Table.h b/Hash_Table.h index 0ac5782..bb0c68e 100644 --- a/Hash_Table.h +++ b/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); diff --git a/Makefile b/Makefile index 383a714..fc0fe01 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ CC=g++ CFLAGS = -w -c -msse4.2 -mpopcnt -fomit-frame-pointer -Winline -O3 -lz LDFLAGS = -lm -lz -lpthread -O3 -mpopcnt -msse4.2 -lz -w -SOURCES = main.cpp CommandLines.cpp Process_Read.cpp Assembly.cpp kmer.cpp Hash_Table.cpp +SOURCES = main.cpp CommandLines.cpp Process_Read.cpp Assembly.cpp kmer.cpp Hash_Table.cpp POA.cpp OBJECTS = $(SOURCES:.c=.o) EXECUTABLE = ccs_assembly diff --git a/POA.cpp b/POA.cpp new file mode 100644 index 0000000..a65f8bf --- /dev/null +++ b/POA.cpp @@ -0,0 +1,464 @@ +#include "POA.h" +#include +#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; + +} diff --git a/POA.h b/POA.h new file mode 100644 index 0000000..9571cf3 --- /dev/null +++ b/POA.h @@ -0,0 +1,104 @@ +#ifndef __POA_PARSER__ +#define __POA_PARSER__ +#include +#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 \ No newline at end of file diff --git a/Process_Read.cpp b/Process_Read.cpp index eafa1fe..d3d4e03 100644 --- a/Process_Read.cpp +++ b/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]];\ diff --git a/Process_Read.h b/Process_Read.h index da75bf0..9623601 100644 --- a/Process_Read.h +++ b/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);