diff --git a/.vscode/ipch/756c62c8d5e4771d/COMMANDLINES.ipch b/.vscode/ipch/756c62c8d5e4771d/COMMANDLINES.ipch index 943f0e4..895d32c 100644 Binary files a/.vscode/ipch/756c62c8d5e4771d/COMMANDLINES.ipch and b/.vscode/ipch/756c62c8d5e4771d/COMMANDLINES.ipch differ diff --git a/.vscode/ipch/856c69690c514632/KMER.ipch b/.vscode/ipch/856c69690c514632/KMER.ipch index f25e0b1..0f60803 100644 Binary files a/.vscode/ipch/856c69690c514632/KMER.ipch and b/.vscode/ipch/856c69690c514632/KMER.ipch differ diff --git a/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch b/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch index 799bcc0..20f8933 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 08d61ef..6ef201d 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 14984c6..8a465a3 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 f9164e7..9c8b6a7 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -33,6 +33,8 @@ void* Perform_Counting(void* arg) uint64_t code; + uint64_t end_pos; + Hash_code k_code; int avalible_k = 0; @@ -54,7 +56,7 @@ void* Perform_Counting(void* arg) avalible_k = 0; - while ((code = get_HPC_code(&HPC_read)) != 6) + while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6) { if(code < 4) { @@ -90,7 +92,7 @@ void* Perform_Counting(void* arg) avalible_k = 0; - while ((code = get_HPC_code(&HPC_read)) != 6) + while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6) { if(code < 4) { @@ -148,8 +150,9 @@ void* Build_hash_table(void* arg) int file_flag = 1; uint64_t code; + uint64_t end_pos; - long long HPC_base; + ///long long HPC_base; Hash_code k_code; @@ -171,9 +174,9 @@ void* Build_hash_table(void* arg) avalible_k = 0; - HPC_base = 0; + ///HPC_base = 0; - while ((code = get_HPC_code(&HPC_read)) != 6) + while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6) { if(code < 4) { @@ -189,7 +192,7 @@ void* Build_hash_table(void* arg) ///所以如果当前k-mer在PCB表中存在,则他的位置一定要加入到候选位置中去 ///insert_Total_Pos_Table(&PCB, &k_code, k_mer_length, curr_sub_block.read[i].ID, HPC_base - k_mer_length + 1, FORWARD); insert_Total_Pos_Table(&PCB, &k_code, k_mer_length, - curr_sub_block.read[i].ID, HPC_base - k_mer_length + 1); + curr_sub_block.read[i].ID, end_pos); } } @@ -199,7 +202,7 @@ void* Build_hash_table(void* arg) init_Hash_code(&k_code); } - HPC_base++; + ///HPC_base++; } @@ -222,7 +225,7 @@ void* Build_hash_table(void* arg) HPC_base = 0; - while ((code = get_HPC_code(&HPC_read)) != 6) + while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6) { if(code < 4) { @@ -411,6 +414,145 @@ void Build_hash_table_multiple_thr() } + +void* Overlap_calculate(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; + + k_mer_pos_list* merge_list; + + init_Candidates_list(&l); + + 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) + { + /** + if (i % 1000 == 0) + { + fprintf(stderr, "i: %llu\n", i); + } + **/ + + + + clear_Candidates_list(&l); + 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); + merge_Candidates_list(&l, list, list_length, end_pos, 0); + //merge_Candidates_list_version(&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); + merge_Candidates_list(&l, list, list_length, end_pos, 1); + //merge_Candidates_list_version(&l, list, list_length, end_pos, 1); + } + } + else + { + avalible_k = 0; + init_Hash_code(&k_code); + } + + ///HPC_base++; + } + + } + + destory_Candidates_list(&l); +} + + +void Overlap_calculate_multipe_thr() +{ + double start_time = Get_T(); + + + fprintf(stdout, "Begin Overlap Calculate ...... \n"); + + pthread_t *_r_threads; + + _r_threads = (pthread_t *)malloc(sizeof(pthread_t)*thread_num); + + int i = 0; + + for (i = 0; i < thread_num; i++) + { + int *arg = (int*)malloc(sizeof(*arg)); + *arg = i; + + pthread_create(_r_threads + i, NULL, Overlap_calculate, (void*)arg); + + } + + + for (i = 0; isize = 1000; + list->length = 0; + list->list = (k_mer_pos_list*)malloc(sizeof(k_mer_pos_list)*list->size); +} + +void destory_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) +{ + free(list->list); +} + +void clear_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) +{ + list->length = 0; +} + + +void append_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length) +{ + list->length++; + if (list->length > list->size) + { + list->size = list->size * 2; + list->list = (k_mer_pos_list*)realloc(list->list, sizeof(k_mer_pos_list)*list->size); + } + + list->list[list->length - 1].list = n_list; + list->list[list->length - 1].length = n_length; +} + void init_Count_Table(Count_Table** table) { *table = kh_init(COUNT64); @@ -300,6 +333,386 @@ int cmp_k_mer_pos(const void * a, const void * b) } +void init_Candidates_list(Candidates_list* l) +{ + l->length = 0; + l->size = 0; + l->list = NULL; + l->tmp = NULL; + l->foward_pos = 0; + l->rc_pos = 0; +} + +void clear_Candidates_list(Candidates_list* l) +{ + l->length = 0; + l->foward_pos = 0; + l->rc_pos = 0; +} + +void destory_Candidates_list(Candidates_list* l) +{ + free(l->list); + free(l->tmp); +} + +///1是x小,2是y小,0是相等 +inline int cmp_pos(k_mer_hit* x, k_mer_pos* y, uint64_t y_pos, uint8_t strand) +{ + if (x->strand < strand) + { + return 1; + } + else if (x->strand > strand) + { + return 2; + } + else + { + if (x->readID < y->readID) + { + return 1; + } + else if (x->readID > y->readID) + { + return 2; + } + else + { + if (x->offset < y->offset) + { + return 1; + } + else if (x->offset > y->offset) + { + return 2; + } + else + { + if (x->self_offset < y_pos) + { + return 1; + } + else if (x->self_offset > y_pos) + { + return 2; + } + else + { + return 0; + } + + } + } + } + +} + +void debug_merge_Candidates_list(Candidates_list* l, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, uint64_t i, +uint64_t i_1, uint64_t i_2, uint64_t len1, uint64_t len2, uint64_t strand) +{ + if (i != l->length) + { + fprintf(stderr, "ERROR\n"); + fprintf(stderr, "i: %llu, l->length: %llu\n", i , l->length); + fprintf(stderr, "i_1: %llu, len1: %llu\n", i_1 , len1); + fprintf(stderr, "i_2: %llu, len2: %llu\n", i_2 , len2); + fprintf(stderr, "strand: %llu\n", strand); + + } + + for (i = 1; i < l->length; i++) + { + if (l->list[i].strand < l->list[i-1].strand) + { + fprintf(stderr, "ERROR -1\n"); + } + else if (l->list[i].strand == l->list[i-1].strand) + { + if (l->list[i].readID < l->list[i-1].readID) + { + fprintf(stderr, "ERROR 0\n"); + } + else if (l->list[i].readID == l->list[i-1].readID) + { + if (l->list[i].offset < l->list[i-1].offset) + { + fprintf(stderr, "ERROR 1\n"); + } + else if (l->list[i].offset == l->list[i-1].offset) + { + if (l->list[i].self_offset < l->list[i-1].self_offset) + { + fprintf(stderr, "ERROR 2\n"); + } + } + } + } + } + + for (i = 0; i < n_lengh; i++) + { + int j = 0; + for (j = 0; j < l->length; j++) + { + if ( + n_list[i].offset == l->list[j].offset + && + n_list[i].readID == l->list[j].readID + && + end_pos == l->list[j].self_offset + && + strand == l->list[j].strand + ) + { + break; + } + } + + if (j == l->length) + { + fprintf(stderr, "ERROR 4\n"); + } + } + + + for (i = 0; i < l->length - n_lengh; i++) + { + int j = 0; + for (j = 0; j < l->length; j++) + { + if ( + l->tmp[i].offset == l->list[j].offset + && + l->tmp[i].readID == l->list[j].readID + && + l->tmp[i].self_offset == l->list[j].self_offset + && + l->tmp[i].strand == l->list[j].strand + ) + { + break; + } + } + + if (j == l->length) + { + fprintf(stderr, "ERROR 5\n"); + } + } +} + +void merge_Candidates_list(Candidates_list* l, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, int strand) +{ + if (n_lengh == 0) ///不加这个就会巨慢 + { + return; + } + + + if(l->length + n_lengh > l->size) + { + l->size = l->length + n_lengh; + l->list = (k_mer_hit*)realloc(l->list, sizeof(k_mer_hit)*l->size); + l->tmp = (k_mer_hit*)realloc(l->tmp, sizeof(k_mer_hit)*l->size); + } + + int flag; + uint64_t len1, len2, i_1, i_2, i; + k_mer_hit *l_1; + k_mer_pos *l_2; + + l_1 = l->list; + if (strand == 0) + { + ///由于标号全是0,则标号为1的位置不受影响(数组后半部分),可以直接复制过去 + ///需要从l_1 + l->rc_pos开始,复制l->length-l->rc_pos个元素,到l->tmp + l->rc_pos + n_lengh处 + memcpy(l->tmp + l->rc_pos + n_lengh, l_1 + l->rc_pos, sizeof(k_mer_hit)*(l->length-l->rc_pos)); + ///l_1只需要从0扫描到l->rc_pos就好了 + len1 = l->rc_pos; + i_1 = 0; + ///此时目标数组没数据, + i = 0; + ///更新l->rc_pos + l->rc_pos = l->rc_pos + n_lengh; + } + else + { + ///如果strand==1,则标号为0的位置不受影响(数组前半部分),可以直接复制过去 + ///从l_1开始,复制l->rc_pos个元素,到l->tmp + memcpy(l->tmp, l_1, sizeof(k_mer_hit)*l->rc_pos); + ///l_1只需要从l->rc_pos扫描到末尾就好了 + i_1 = l->rc_pos; + len1 = l->length; + ///此时目标数组已经有了l->rc_pos + i = i_1; + ///l->rc_pos不用更新 + } + + + + + l_2 = n_list; + len2 = n_lengh; + i_2 = 0; + + + while (i_1tmp[i].readID = l_1[i_1].readID; + l->tmp[i].offset = l_1[i_1].offset; + l->tmp[i].self_offset = l_1[i_1].self_offset; + l->tmp[i].strand = l_1[i_1].strand; + i_1++; + i++; + } + else ///l_2小或者l_1 == l_2 + { + l->tmp[i].readID = l_2[i_2].readID; + l->tmp[i].offset = l_2[i_2].offset; + l->tmp[i].self_offset = end_pos; + l->tmp[i].strand = strand; + i_2++; + i++; + } + } + + + while (i_1tmp + i, l_1 + i_1, sizeof(k_mer_hit)*(len1-i_1)); + i = i + len1-i_1; + i_1 = len1; + } + + while (i_2tmp[i].readID = l_2[i_2].readID; + l->tmp[i].offset = l_2[i_2].offset; + l->tmp[i].self_offset = end_pos; + l->tmp[i].strand = strand; + i_2++; + i++; + } + + + + k_mer_hit* k; + k = l->list; + l->list = l->tmp; + l->tmp = k; + l->length = l->length + n_lengh; + + + /** + if (strand == 0) + { + i = l->length; + len1 = l->length - n_lengh; + } + debug_merge_Candidates_list(l, n_list, n_lengh, end_pos, i, i_1, i_2, len1, len2, strand); + **/ + +} + + + + + + + +/********************************for debug***************************************/ +void merge_Candidates_list_version(Candidates_list* l, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, int strand) +{ + if (n_lengh == 0) ///不加这个就会巨慢 + { + return; + } + + + if(l->length + n_lengh > l->size) + { + l->size = l->length + n_lengh; + l->list = (k_mer_hit*)realloc(l->list, sizeof(k_mer_hit)*l->size); + l->tmp = (k_mer_hit*)realloc(l->tmp, sizeof(k_mer_hit)*l->size); + } + + int flag; + uint64_t len1, len2, i_1, i_2, i; + k_mer_hit *l_1; + k_mer_pos *l_2; + + l_1 = l->list; + len1 = l->length; + i_1 = 0; + + l_2 = n_list; + len2 = n_lengh; + i_2 = 0; + + i = 0; + + while (i_1tmp[i].readID = l_1[i_1].readID; + l->tmp[i].offset = l_1[i_1].offset; + l->tmp[i].self_offset = l_1[i_1].self_offset; + l->tmp[i].strand = l_1[i_1].strand; + i_1++; + i++; + } + else ///l_2小或者l_1 == l_2 + { + l->tmp[i].readID = l_2[i_2].readID; + l->tmp[i].offset = l_2[i_2].offset; + l->tmp[i].self_offset = end_pos; + l->tmp[i].strand = strand; + i_2++; + i++; + } + } + + + while (i_1tmp[i].readID = l_1[i_1].readID; + l->tmp[i].offset = l_1[i_1].offset; + l->tmp[i].self_offset = l_1[i_1].self_offset; + l->tmp[i].strand = l_1[i_1].strand; + i_1++; + i++; + } + + while (i_2tmp[i].readID = l_2[i_2].readID; + l->tmp[i].offset = l_2[i_2].offset; + l->tmp[i].self_offset = end_pos; + l->tmp[i].strand = strand; + i_2++; + i++; + } + + k_mer_hit* k; + k = l->list; + l->list = l->tmp; + l->tmp = k; + l->length = l->length + n_lengh; + + ///debug_merge_Candidates_list(l, n_list, n_lengh, end_pos, i, i_1, i_2, len1, len2, strand); +} + + + + diff --git a/Hash_Table.h b/Hash_Table.h index 929b866..cac686c 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -37,6 +37,39 @@ typedef struct uint64_t readID; } k_mer_pos; +typedef struct +{ + k_mer_pos* list; + uint64_t length; +} k_mer_pos_list; + +typedef struct +{ + k_mer_pos_list* list; + uint64_t size; + uint64_t length; +} k_mer_pos_list_alloc; + + + +typedef struct +{ + uint64_t offset; + uint64_t readID; + uint64_t self_offset; + uint8_t strand; +} k_mer_hit; + +typedef struct +{ + k_mer_hit* list; + k_mer_hit* tmp; + long long length; + long long size; + uint64_t foward_pos; + uint64_t rc_pos; +} Candidates_list; + typedef struct { Pos_Table** sub_h; @@ -273,8 +306,8 @@ inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, in { qsort(list, occ, sizeof(k_mer_pos), cmp_k_mer_pos); } - + return 1; } else @@ -300,6 +333,10 @@ int load_Total_Pos_Table(Total_Pos_Table* TCB, char* read_file_name); void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k_mer_min_freq, int k_mer_max_freq); +void init_Candidates_list(Candidates_list* l); +void clear_Candidates_list(Candidates_list* l); +void destory_Candidates_list(Candidates_list* l); +void merge_Candidates_list(Candidates_list* l, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, int strand); @@ -429,4 +466,8 @@ void test_COUNT64(); void debug_mode(uint64_t d, uint64_t thread_ID, uint64_t thread_num); + +/********************************for debug***************************************/ +void merge_Candidates_list_version(Candidates_list* l, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, int strand); + #endif \ No newline at end of file diff --git a/kmer.h b/kmer.h index 9b1e5c7..c3cdab7 100644 --- a/kmer.h +++ b/kmer.h @@ -33,13 +33,15 @@ typedef struct } HPC_seq; -inline uint64_t get_HPC_code(HPC_seq* seq) +inline uint64_t get_HPC_code(HPC_seq* seq, uint64_t* end_pos) { if(seq->i < seq ->l) { uint8_t code = seq_nt6_table[(uint8_t)seq->str[seq->i]]; + (*end_pos) = seq->i; + for (; seq->i < seq->l; seq->i++) { ///统计N的个数 diff --git a/main.cpp b/main.cpp index ce53b95..1b55d49 100644 --- a/main.cpp +++ b/main.cpp @@ -35,11 +35,9 @@ int main(int argc, char *argv[]) Build_hash_table_multiple_thr(); } - - - ///verify_Position_hash_table(); + Overlap_calculate_multipe_thr();