diff --git a/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch b/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch index b2273db..1e53d27 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 eba1c97..75e1fe4 100644 Binary files a/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch and b/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch differ diff --git a/Assembly.cpp b/Assembly.cpp index 0779b18..1ae84bb 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -575,6 +575,7 @@ void* Overlap_calculate(void* arg) k_mer_pos_list_alloc array_list; init_k_mer_pos_list_alloc(&array_list); + for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num) ///for (i = thr_ID; i < R_INF.total_reads/50; i = i + thread_num) { @@ -674,6 +675,144 @@ void* Overlap_calculate(void* arg) } +void* Overlap_calculate_heap_merge(void* arg) +{ + + int thr_ID = *((int*)arg); + + long long i = 0; + int avalible_k = 0; + + UC_Read g_read; + init_UC_Read(&g_read); + HPC_seq HPC_read; + Hash_code k_code; + uint64_t code; + uint64_t end_pos; + k_mer_pos* list; + uint64_t list_length; + uint64_t sub_ID; + + Candidates_list l; + //Candidates_list debug_l; + + init_Candidates_list(&l); + //init_Candidates_list(&debug_l); + + k_mer_pos_list_alloc array_list; + init_k_mer_pos_list_alloc(&array_list); + + HeapSq heap; + + Init_Heap(&heap); + + for (i = thr_ID; i < R_INF.total_reads; i = i + thread_num) + { + + clear_Heap(&heap); + clear_Candidates_list(&l); + ///clear_Candidates_list(&debug_l); + + clear_k_mer_pos_list_alloc(&array_list); + recover_UC_Read(&g_read, &R_INF, i); + + + ///forward strand + init_HPC_seq(&HPC_read, g_read.seq, g_read.length); + init_Hash_code(&k_code); + avalible_k = 0; + + while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6) + { + if(code < 4) + { + k_mer_append(&k_code,code,k_mer_length); + avalible_k++; + if (avalible_k>=k_mer_length) + { + list_length = locate_Total_Pos_Table(&PCB, &k_code, &list, k_mer_length, &sub_ID); + + if (list_length != 0) + { + append_k_mer_pos_list_alloc(&array_list, list, list_length, end_pos, 0); + } + ///merge_Candidates_list(&l, list, list_length, end_pos, 0); + //merge_Candidates_list_version(&debug_l, list, list_length, end_pos, 0); + } + } + else + { + avalible_k = 0; + init_Hash_code(&k_code); + } + + ///HPC_base++; + } + + + + ///reverse complement strand + reverse_complement(g_read.seq, g_read.length); + init_HPC_seq(&HPC_read, g_read.seq, g_read.length); + init_Hash_code(&k_code); + avalible_k = 0; + + while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6) + { + if(code < 4) + { + k_mer_append(&k_code,code,k_mer_length); + avalible_k++; + if (avalible_k>=k_mer_length) + { + list_length = locate_Total_Pos_Table(&PCB, &k_code, &list, k_mer_length, &sub_ID); + if (list_length != 0) + { + append_k_mer_pos_list_alloc(&array_list, list, list_length, end_pos, 1); + } + ///merge_Candidates_list(&l, list, list_length, end_pos, 1); + //merge_Candidates_list_version(&debug_l, list, list_length, end_pos, 1); + } + } + else + { + avalible_k = 0; + init_Hash_code(&k_code); + } + + ///HPC_base++; + } + + ///merge_k_mer_pos_list_alloc(&array_list, &l); + merge_k_mer_pos_list_alloc_heap_sort(&array_list, &l, &heap); + /** + if (array_list.length < 3) + { + merge_k_mer_pos_list_alloc(&array_list, &l); + } + else + { + merge_k_mer_pos_list_alloc_heap_sort_advance(&array_list, &l, &heap); + } + **/ + + + ///merge_k_mer_pos_list_alloc_heap_sort_advance(&array_list, &l, &heap); + + + + ///debug_merge_result(&l, &debug_l); + + } + + destory_Candidates_list(&l); + //destory_Candidates_list(&debug_l); + + destory_Heap(&heap); + destory_k_mer_pos_list_alloc(&array_list); +} + + @@ -698,7 +837,8 @@ void Overlap_calculate_multipe_thr() int *arg = (int*)malloc(sizeof(*arg)); *arg = i; - pthread_create(_r_threads + i, NULL, Overlap_calculate, (void*)arg); + pthread_create(_r_threads + i, NULL, Overlap_calculate_heap_merge, (void*)arg); + //pthread_create(_r_threads + i, NULL, Overlap_calculate, (void*)arg); } diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 245a403..36dbb19 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -4,6 +4,147 @@ #include "Hash_Table.h" +void Init_Heap(HeapSq* HBT) +{ + HBT->MaxSize = 1000; + HBT->heap = (ElemType*)malloc(HBT->MaxSize*sizeof(ElemType)); + HBT->index_i = (uint64_t*)malloc(HBT->MaxSize*sizeof(uint64_t)); + HBT->len = 0; +} + +void destory_Heap(HeapSq* HBT) +{ + free(HBT->heap); + free(HBT->index_i); +} + +void clear_Heap(HeapSq* HBT) +{ + HBT->len = 0; +} + + +inline int cmp_ElemType(ElemType* x, ElemType* y) +{ + if (x->node.strand < y->node.strand) + { + return 1; + } + else if (x->node.strand > y->node.strand) + { + return 2; + } + else + { + if (x->node.readID < y->node.readID) + { + return 1; + } + else if (x->node.readID > y->node.readID) + { + return 2; + } + else + { + if (x->node.offset < y->node.offset) + { + return 1; + } + else if (x->node.offset > y->node.offset) + { + return 2; + } + else + { + if (x->node.self_offset < y->node.self_offset) + { + return 1; + } + else if (x->node.self_offset > y->node.self_offset) + { + return 2; + } + else + { + return 0; + } + + } + } + } + + + + +} + + +inline void Insert_Heap(HeapSq* HBT, ElemType* x) +{ + long long i, j; + if (HBT->len == HBT->MaxSize) //若堆满,将数组空间扩展为原来的2倍 + { + HBT->MaxSize = 2*HBT->MaxSize; + HBT->heap = (ElemType*)realloc(HBT->heap, HBT->MaxSize*sizeof(ElemType)); + HBT->index_i = (uint64_t*)realloc(HBT->index_i,HBT->MaxSize*sizeof(uint64_t)); + } + HBT->heap[HBT->len] = *x; //向堆尾添加新元素 + HBT->len++; //堆长度加1 + i = HBT->len - 1; //i指向待调整元素的位置,即其数组下标,初始指向新元素所在的堆尾位置 + while (i != 0) + { + j = (i - 1) / 2; //j指向下标为i的元素的双亲 + ///if (x >= HBT->heap[j]) //若新元素大于待调整元素的双亲,则比较调整结束,退出循环 + ///1: xheap[j])!=1) + break; + HBT->heap[i] = HBT->heap[j]; //将双亲元素下移到待调整元素的位置 + i = j; //使待调整位置变为其双亲位置,进行下一次循环 + } + HBT->heap[i] = *x;//把新元素调整到最终位置 +} + + +inline int DeleteHeap(HeapSq* HBT, ElemType* get) +{ + ElemType temp, x; + int i, j; + if (HBT->len == 0) + { + return 0; + } + temp = HBT->heap[0]; //暂存堆顶元素 + HBT->len--; + if (HBT->len == 0) //若删除操作后堆为空则返回 + { + *get = temp; + return 2; + } + + x = HBT->heap[HBT->len]; //将待调整的原堆尾元素暂存x中,以便放入最终位置 + i = 0; //用i指向待调整元素的位置,初始指向堆顶位置 + j = 2 * i + 1;//用j指向i的左孩子位置,初始指向下标为1的位置 + while (j <= HBT->len - 1)//寻找待调整元素的最终位置,每次使孩子元素上移一层,调整到孩子为空时止 + { + ///if (j < HBT->len - 1 && HBT->heap[j] > HBT->heap[j+1])//若存在右孩子且较小,使j指向右孩子 + if (j < HBT->len - 1 && cmp_ElemType(&HBT->heap[j], &HBT->heap[j + 1]) == 2) + j++; + ///if (x <= HBT->heap[j]) //若x比其较小的孩子还小,则调整结束,退出循环 + if (cmp_ElemType(&x, &HBT->heap[j])!=2) + break; + HBT->heap[i] = HBT->heap[j];//否则,将孩子元素移到双亲位置 + i = j; //将待调整位置变为其较小的孩子位置 + j = 2 * i + 1;//将j变为新的待调整位置的左孩子位置,继续下一次循环 + } + HBT->heap[i] = x; //把x放到最终位置 + + //返回原堆顶元素 + *get = temp; + return 1; +} + + + void init_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) @@ -59,6 +200,396 @@ int cmp_k_mer_pos_list(const void * a, const void * b) } +inline void append_pos_to_Candidates_list(Candidates_list* candidates, ElemType* x) +{ + candidates->list[candidates->length] = x->node; + candidates->length++; +} + + +void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, uint64_t strand) +{ + uint64_t i; + + int j = 0; + + for (i = 0; i < n_lengh; i++) + { + + for (; j < candidates->length; j++) + { + if ( + n_list[i].offset == candidates->list[j].offset + && + n_list[i].readID == candidates->list[j].readID + && + end_pos == candidates->list[j].self_offset + && + strand == candidates->list[j].strand + ) + { + break; + } + } + + if (j == candidates->length) + { + fprintf(stderr, "ERROR 4\n"); + } + } + +} + +void verify_merge_result(k_mer_pos_list_alloc* list, Candidates_list* candidates) +{ + uint64_t total_length = 0; + uint64_t i = 0; + for (i = 0; i < list->length; i++) + { + total_length = total_length + list->list[i].length; + } + + if (total_length!=candidates->length) + { + fprintf(stderr, "ERROR length & size.\n"); + } + + for (i = 1; i < candidates->length; i++) + { + if (candidates->list[i].strand < candidates->list[i-1].strand) + { + fprintf(stderr, "ERROR -1\n"); + } + else if (candidates->list[i].strand == candidates->list[i-1].strand) + { + if (candidates->list[i].readID < candidates->list[i-1].readID) + { + fprintf(stderr, "ERROR 0\n"); + } + else if (candidates->list[i].readID == candidates->list[i-1].readID) + { + if (candidates->list[i].offset < candidates->list[i-1].offset) + { + fprintf(stderr, "ERROR 1\n"); + } + else if (candidates->list[i].offset == candidates->list[i-1].offset) + { + if (candidates->list[i].self_offset < candidates->list[i-1].self_offset) + { + fprintf(stderr, "ERROR 2\n"); + } + } + } + } + } + + + uint64_t j = 0; + for (i = 0; i < list->length; i++) + { + test_single_list(candidates, list->list[i].list, list->list[i].length, list->list[i].end_pos, list->list[i].direction); + } + + + +} + + + +void merge_k_mer_pos_list_alloc_heap_sort_back(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT) +{ + clear_Heap(HBT); + + uint64_t total_length = 0; + uint64_t i; + ElemType x, y; + ///所有list的长度都不是0 + ///把各个表第一个元素加入到堆中 + for (i = 0; i < list->length; i++) + { + x.ID = i; + x.node.offset = list->list[i].list[0].offset; + x.node.readID = list->list[i].list[0].readID; + x.node.self_offset = list->list[i].end_pos; + x.node.strand = list->list[i].direction; + + Insert_Heap(HBT, &x); + + HBT->index_i[i] = 1; + + total_length = total_length + list->list[i].length; + } + + + candidates->length = 0; + if(total_length > candidates->size) + { + candidates->size = total_length; + candidates->list = (k_mer_hit*)realloc(candidates->list, sizeof(k_mer_hit)*candidates->size); + candidates->tmp = (k_mer_hit*)realloc(candidates->tmp, sizeof(k_mer_hit)*candidates->size); + } + + uint64_t ID; + while (DeleteHeap(HBT, &x)) + { + append_pos_to_Candidates_list(candidates, &x); + ///x.ID说明是从第x.ID个列表中的这个节点已经从堆里出来了 + ///HBT->index_i[x.ID]是第x.ID个列表的当前元素的下标 + i = HBT->index_i[x.ID]; + ID = x.ID; + ///fprintf(stderr, "x.ID: %llu, i: %llu, length: %llu\n", x.ID, i, list->list[x.ID].length); + + + if (i < list->list[x.ID].length) + { + y.ID = ID; + y.node.offset = list->list[x.ID].list[i].offset; + y.node.readID = list->list[x.ID].list[i].readID; + y.node.self_offset = list->list[x.ID].end_pos; + y.node.strand = list->list[x.ID].direction; + Insert_Heap(HBT, &y); + HBT->index_i[ID]++; + } + } + + ///verify_merge_result(list, candidates); + + +} + + +void merge_k_mer_pos_list_alloc_heap_sort(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT) +{ + clear_Heap(HBT); + + uint64_t total_length = 0; + uint64_t i; + ElemType x, y; + ///所有list的长度都不是0 + ///把各个表第一个元素加入到堆中 + for (i = 0; i < list->length; i++) + { + + x.ID = i; + x.node.offset = list->list[i].list[0].offset; + x.node.readID = list->list[i].list[0].readID; + x.node.self_offset = list->list[i].end_pos; + x.node.strand = list->list[i].direction; + + Insert_Heap(HBT, &x); + + HBT->index_i[i] = 1; + + + total_length = total_length + list->list[i].length; + } + + + candidates->length = 0; + if(total_length > candidates->size) + { + candidates->size = total_length; + candidates->list = (k_mer_hit*)realloc(candidates->list, sizeof(k_mer_hit)*candidates->size); + candidates->tmp = (k_mer_hit*)realloc(candidates->tmp, sizeof(k_mer_hit)*candidates->size); + } + + uint64_t ID; + int flag; + while (flag = DeleteHeap(HBT, &x)) + { + append_pos_to_Candidates_list(candidates, &x); + ///x.ID说明是从第x.ID个列表中的这个节点已经从堆里出来了 + ///HBT->index_i[x.ID]是第x.ID个列表的当前元素的下标 + i = HBT->index_i[x.ID]; + ID = x.ID; + + if (flag == 2) + { + for (; i < list->list[x.ID].length; i++) + { + y.ID = ID; + y.node.offset = list->list[x.ID].list[i].offset; + y.node.readID = list->list[x.ID].list[i].readID; + y.node.self_offset = list->list[x.ID].end_pos; + y.node.strand = list->list[x.ID].direction; + append_pos_to_Candidates_list(candidates, &y); + } + + break; + } + + if (i < list->list[x.ID].length) + { + y.ID = ID; + y.node.offset = list->list[x.ID].list[i].offset; + y.node.readID = list->list[x.ID].list[i].readID; + y.node.self_offset = list->list[x.ID].end_pos; + y.node.strand = list->list[x.ID].direction; + Insert_Heap(HBT, &y); + HBT->index_i[ID]++; + } + } + + + ///verify_merge_result(list, candidates); + + +} + + +///有bug,找时间排一下 +void merge_k_mer_pos_list_alloc_heap_sort_advance(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT) +{ + + uint64_t total_length = 0; + uint64_t i; + ElemType x, y; + + /************************************forward*************************************************/ + clear_Heap(HBT); + + ///所有list的长度都不是0 + ///把各个表第一个元素加入到堆中 + for (i = 0; i < list->length; i++) + { + if (list->list[i].direction == 0) + { + x.ID = i; + x.node.offset = list->list[i].list[0].offset; + x.node.readID = list->list[i].list[0].readID; + x.node.self_offset = list->list[i].end_pos; + x.node.strand = list->list[i].direction; + + Insert_Heap(HBT, &x); + + HBT->index_i[i] = 1; + + total_length = total_length + list->list[i].length; + } + } + + + candidates->length = 0; + if(total_length > candidates->size) + { + candidates->size = total_length; + candidates->list = (k_mer_hit*)realloc(candidates->list, sizeof(k_mer_hit)*candidates->size); + candidates->tmp = (k_mer_hit*)realloc(candidates->tmp, sizeof(k_mer_hit)*candidates->size); + } + + + uint64_t ID; + int flag; + while (flag = DeleteHeap(HBT, &x)) + { + append_pos_to_Candidates_list(candidates, &x); + ///x.ID说明是从第x.ID个列表中的这个节点已经从堆里出来了 + ///HBT->index_i[x.ID]是第x.ID个列表的当前元素的下标 + i = HBT->index_i[x.ID]; + ID = x.ID; + ///fprintf(stderr, "x.ID: %llu, i: %llu, length: %llu\n", x.ID, i, list->list[x.ID].length); + /** + if (flag == 2) + { + for (; i < list->list[x.ID].length; i++) + { + y.ID = ID; + y.node.offset = list->list[x.ID].list[i].offset; + y.node.readID = list->list[x.ID].list[i].readID; + y.node.self_offset = list->list[x.ID].end_pos; + y.node.strand = list->list[x.ID].direction; + append_pos_to_Candidates_list(candidates, &y); + } + + break; + } + **/ + + + if (i < list->list[x.ID].length) + { + y.ID = ID; + y.node.offset = list->list[x.ID].list[i].offset; + y.node.readID = list->list[x.ID].list[i].readID; + y.node.self_offset = list->list[x.ID].end_pos; + y.node.strand = list->list[x.ID].direction; + Insert_Heap(HBT, &y); + HBT->index_i[ID]++; + } + } + + + + /************************************reverse complement*************************************************/ + + clear_Heap(HBT); + ///所有list的长度都不是0 + ///把各个表第一个元素加入到堆中 + for (i = 0; i < list->length; i++) + { + if (list->list[i].direction == 1) + { + x.ID = i; + x.node.offset = list->list[i].list[0].offset; + x.node.readID = list->list[i].list[0].readID; + x.node.self_offset = list->list[i].end_pos; + x.node.strand = list->list[i].direction; + + Insert_Heap(HBT, &x); + + HBT->index_i[i] = 1; + + total_length = total_length + list->list[i].length; + } + } + + if(total_length > candidates->size) + { + candidates->size = total_length; + candidates->list = (k_mer_hit*)realloc(candidates->list, sizeof(k_mer_hit)*candidates->size); + candidates->tmp = (k_mer_hit*)realloc(candidates->tmp, sizeof(k_mer_hit)*candidates->size); + } + + while (flag = DeleteHeap(HBT, &x)) + { + append_pos_to_Candidates_list(candidates, &x); + ///x.ID说明是从第x.ID个列表中的这个节点已经从堆里出来了 + ///HBT->index_i[x.ID]是第x.ID个列表的当前元素的下标 + i = HBT->index_i[x.ID]; + ID = x.ID; + ///fprintf(stderr, "x.ID: %llu, i: %llu, length: %llu\n", x.ID, i, list->list[x.ID].length); + /** + if (flag == 2) + { + for (; i < list->list[x.ID].length; i++) + { + y.ID = ID; + y.node.offset = list->list[x.ID].list[i].offset; + y.node.readID = list->list[x.ID].list[i].readID; + y.node.self_offset = list->list[x.ID].end_pos; + y.node.strand = list->list[x.ID].direction; + append_pos_to_Candidates_list(candidates, &y); + } + + break; + } + **/ + + + if (i < list->list[x.ID].length) + { + y.ID = ID; + y.node.offset = list->list[x.ID].list[i].offset; + y.node.readID = list->list[x.ID].list[i].readID; + y.node.self_offset = list->list[x.ID].end_pos; + y.node.strand = list->list[x.ID].direction; + Insert_Heap(HBT, &y); + HBT->index_i[ID]++; + } + } + +} + void merge_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list, Candidates_list* candidates) { diff --git a/Hash_Table.h b/Hash_Table.h index 2ad3bb1..3d001f1 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -62,6 +62,24 @@ typedef struct uint8_t strand; } k_mer_hit; + +typedef struct +{ + k_mer_hit node; + uint64_t ID; +} ElemType; + + +typedef struct +{ + ElemType* heap; + uint64_t* index_i; + int len; + int MaxSize; +} HeapSq; + + + typedef struct { k_mer_hit* list; @@ -347,6 +365,12 @@ void clear_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list); void append_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length, uint64_t n_end_pos, uint8_t n_direction); void merge_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list, Candidates_list* candidates); +void merge_k_mer_pos_list_alloc_heap_sort(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT); +void merge_k_mer_pos_list_alloc_heap_sort_advance(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT); + +void Init_Heap(HeapSq* HBT); +void destory_Heap(HeapSq* HBT); +void clear_Heap(HeapSq* HBT);