diff --git a/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch b/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch index 1c7d690..f38974f 100644 Binary files a/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch and b/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch differ diff --git a/.vscode/ipch/756c62c8d5e4771d/COMMANDLINES.ipch b/.vscode/ipch/756c62c8d5e4771d/COMMANDLINES.ipch index 3caee36..25872b9 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 5e6bbb1..9114f47 100644 Binary files a/.vscode/ipch/856c69690c514632/KMER.ipch and b/.vscode/ipch/856c69690c514632/KMER.ipch differ diff --git a/.vscode/ipch/920fff44ce64d7c2/mmap_address.bin b/.vscode/ipch/920fff44ce64d7c2/mmap_address.bin new file mode 100644 index 0000000..862b842 Binary files /dev/null and b/.vscode/ipch/920fff44ce64d7c2/mmap_address.bin differ diff --git a/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch b/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch index 510b9a5..99f78ac 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 c11c872..f8248ca 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 fc39927..b61119d 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 1caed77..c4a693b 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -8,73 +8,8 @@ #include "Hash_Table.h" Total_Count_Table TCB; - - - -/********************************for debug***************************************/ -void Verify_Counting() -{ - - - kseq_t *seq = (kseq_t*)calloc(1, sizeof(kseq_t)); - - - long long read_number = 0; - - HPC_seq HPC_read; - uint64_t code; - Hash_code k_code; - int avalible_k = 0; - - fprintf(stdout, "Start Verifying ...\n"); - - while (get_read(seq)) - { - - init_HPC_seq(&HPC_read, seq->seq.s, seq->seq.l); - init_Hash_code(&k_code); - - avalible_k = 0; - - while ((code = get_HPC_code(&HPC_read)) != 6) - { - if(code < 4) - { - k_mer_append(&k_code,code,k_mer_length); - avalible_k++; - if (avalible_k>=k_mer_length) - { - if(verify_Total_Count_Table(&TCB, &k_code, k_mer_length) == -1) - { - fprintf(stderr, "ERROR when subtracting!\n"); - } - } - - } - else - { - avalible_k = 0; - init_Hash_code(&k_code); - } - - } - - - - - read_number++; - } - - fprintf(stdout, "read_number: %lld\n",read_number); - - fprintf(stdout, "Start Traversing ...\n"); - - Traverse_Total_Count_Table(&TCB); - - fprintf(stdout, "Finish Traversing!\n"); - - -} +Total_Pos_Table PCB; +All_reads R_INF; @@ -83,9 +18,6 @@ void* Perform_Counting(void* arg) { int thr_ID = *((int*)arg); - ///debug_mode(101, thr_ID, thread_num); - - int i = 0; HPC_seq HPC_read; @@ -157,9 +89,6 @@ void* Perform_Counting(void* arg) destory_R_buffer_block(&curr_sub_block); free(arg); - fprintf(stdout, "thr_ID: %d, read_number: %lld, select_k_mer_number: %lld, k_mer_number: %lld\n", - thr_ID, read_number, select_k_mer_number, k_mer_number); - } @@ -167,7 +96,94 @@ void* Perform_Counting(void* arg) +void* Build_hash_table(void* arg) +{ + int thr_ID = *((int*)arg); + int i = 0; + HPC_seq HPC_read; + + R_buffer_block curr_sub_block; + + init_R_buffer_block(&curr_sub_block); + + int file_flag = 1; + + uint64_t code; + + long long HPC_base; + + Hash_code k_code; + + int avalible_k = 0; + + char* dest; + + while (file_flag != 0) + { + + file_flag = get_reads_mul_thread(&curr_sub_block); + + + for (i = 0; i < curr_sub_block.num; i++) + { + + dest = R_INF.read+R_INF.index[curr_sub_block.read[i].ID]; + memcpy(dest, curr_sub_block.read[i].seq.s, curr_sub_block.read[i].seq.l); + + dest = R_INF.name+R_INF.name_index[curr_sub_block.read[i].ID]; + memcpy(dest, curr_sub_block.read[i].name.s, curr_sub_block.read[i].name.l); + + + + init_HPC_seq(&HPC_read, curr_sub_block.read[i].seq.s, curr_sub_block.read[i].seq.l); + init_Hash_code(&k_code); + + avalible_k = 0; + + HPC_base = 0; + + while ((code = get_HPC_code(&HPC_read)) != 6) + { + if(code < 4) + { + k_mer_append(&k_code,code,k_mer_length); + avalible_k++; + if (avalible_k>=k_mer_length) + { + + ///选取的k-mer满足两个要求 + ///1. hash(k-mer) % 101 <= 3 + ///2. occ(k-mer)要满足范围 + ///TCB表中的元素仅满足第一个要求，而PCB表中的元素满足两个要求 + ///所以如果当前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); + } + + } + else + { + avalible_k = 0; + init_Hash_code(&k_code); + } + + HPC_base++; + } + + + + } + + + } + + destory_R_buffer_block(&curr_sub_block); + free(arg); + + + +} @@ -180,20 +196,24 @@ void* Perform_Counting(void* arg) void Counting_multiple_thr() { + double start_time = Get_T(); + + fprintf(stdout, "Begin Counting ...... \n"); + + init_kseq(read_file_name); + + init_All_reads(&R_INF); + init_Total_Count_Table(k_mer_length, &TCB); - fprintf(stdout, "TCB.prefix_bits: %d\n", TCB.prefix_bits); - fprintf(stdout, "TCB.suffix_bits: %d\n", TCB.suffix_bits); - fprintf(stdout, "TCB.size: %d\n", TCB.size); - - - pthread_t inputReadsHandle; init_R_buffer(thread_num); - pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, NULL); + int *is_insert = (int*)malloc(sizeof(*is_insert)); + *is_insert = 1; + pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert); pthread_t *_r_threads; @@ -219,11 +239,297 @@ void Counting_multiple_thr() free(_r_threads); - destory_R_buffer(); - - - destory_Total_Count_Table(&TCB); + ///destory_R_buffer(); + ///destory_Total_Count_Table(&TCB); + + destory_kseq(); + + fprintf(stdout, "Finish Counting ...... \n"); + + fprintf(stdout, "%-30s%18.2f\n\n", "Counting time:", Get_T() - start_time); + + +} + + + + +void Build_hash_table_multiple_thr() +{ + double start_time = Get_T(); + + + fprintf(stdout, "Begin Building hash table ...... \n"); + + init_Total_Pos_Table(&PCB, &TCB); + + double T_start_time = Get_T(); + + Traverse_Counting_Table(&TCB, &PCB, k_mer_min_freq, k_mer_max_freq); + + fprintf(stdout, "%-30s%18.2f\n\n", "Traverse time:", Get_T() - T_start_time); + + + + init_kseq(read_file_name); + + clear_R_buffer(); + + malloc_All_reads(&R_INF); + + pthread_t inputReadsHandle; + + int *is_insert = (int*)malloc(sizeof(*is_insert)); + *is_insert = 0; + + pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert); + + + + + 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, Build_hash_table, (void*)arg); + + } + + pthread_join(inputReadsHandle, NULL); + + for (i = 0; iseq.s, seq->seq.l); + init_Hash_code(&k_code); + + avalible_k = 0; + + while ((code = get_HPC_code(&HPC_read)) != 6) + { + if(code < 4) + { + k_mer_append(&k_code,code,k_mer_length); + avalible_k++; + if (avalible_k>=k_mer_length) + { + if(verify_Total_Count_Table(&TCB, &k_code, k_mer_length) == -1) + { + fprintf(stderr, "ERROR when subtracting!\n"); + } + } + + } + else + { + avalible_k = 0; + init_Hash_code(&k_code); + } + + } + + + + + read_number++; + } + + fprintf(stdout, "read_number: %lld\n",read_number); + + fprintf(stdout, "Start Traversing ...\n"); + + Traverse_Total_Count_Table(&TCB); + + fprintf(stdout, "Finish Traversing!\n"); + + +} + + + +/********************************for debug*****************************************/ +void verify_Position_hash_table() +{ + + init_kseq(read_file_name); + + kseq_t *seq = (kseq_t*)calloc(1, sizeof(kseq_t)); + + + long long read_number = 0; + long long HPC_base; + + HPC_seq HPC_read; + uint64_t code; + Hash_code k_code; + int avalible_k = 0; + k_mer_pos* list; + uint64_t sub_ID; + + fprintf(stdout, "Start Verifying Position Table...\n"); + + while (get_read(seq)) + { + + if (seq->seq.l + != Get_READ_LENGTH(R_INF, read_number)) + { + fprintf(stderr, "seq error\n"); + } + + if(memcmp(seq->seq.s, Get_READ(R_INF, read_number), seq->seq.l)) + { + fprintf(stderr, "seq error\n"); + } + + if (seq->name.l + != Get_NAME_LENGTH(R_INF, read_number)) + { + fprintf(stderr, "name error\n"); + } + + + if(memcmp(seq->name.s, Get_NAME(R_INF, read_number), seq->name.l)) + { + fprintf(stderr, "name error\n"); + } + + + + init_HPC_seq(&HPC_read, seq->seq.s, seq->seq.l); + init_Hash_code(&k_code); + + avalible_k = 0; + HPC_base = 0; + + while ((code = get_HPC_code(&HPC_read)) != 6) + { + if(code < 4) + { + k_mer_append(&k_code,code,k_mer_length); + avalible_k++; + if (avalible_k>=k_mer_length) + { + + uint64_t count1 = get_Total_Count_Table(&TCB, &k_code, k_mer_length); + uint64_t count2 = count_Total_Pos_Table(&PCB, &k_code, k_mer_length); + if(count1>=k_mer_min_freq && count1<= k_mer_max_freq && count1!=count2) + { + fprintf(stderr, "count1: %lld\n",count1); + fprintf(stderr, "count2: %lld\n",count2); + } + + + if(locate_Total_Pos_Table(&PCB, &k_code, &list, k_mer_length, &sub_ID) == count2) + { + int j = 0; + for(j=1; j #include #include "ketopt.h" +#include + 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_max_freq = 66; + +double Get_T(void) +{ + struct timeval t; + gettimeofday(&t, NULL); + return t.tv_sec+t.tv_usec/1000000.0; +} void Print_H() { diff --git a/CommandLines.h b/CommandLines.h index f583693..2c3fe58 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -8,8 +8,12 @@ extern char* read_file_name; extern char* output_file_name; extern int thread_num; extern int k_mer_length; +extern int k_mer_min_freq; +extern int k_mer_max_freq; + int CommandLine_process (int argc, char *argv[]); +double Get_T(void); #endif \ No newline at end of file diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 6180dc9..6e7ed70 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -9,6 +9,11 @@ void init_Count_Table(Count_Table** table) *table = kh_init(COUNT64); } +void init_Pos_Table(Pos_Table** table) +{ + *table = kh_init(POS64); +} + void init_Total_Count_Table(int k, Total_Count_Table* TCB) { if(k>64) @@ -34,6 +39,7 @@ void init_Total_Count_Table(int k, Total_Count_Table* TCB) TCB->size = (1ULL<prefix_bits); TCB->sub_h = (Count_Table**)malloc(sizeof(Count_Table*)*TCB->size); TCB->sub_h_lock = (Hash_table_spin_lock*)malloc(sizeof(Hash_table_spin_lock)*TCB->size); + memset(TCB->sub_h_lock, 0, sizeof(Hash_table_spin_lock)*TCB->size); int i = 0; for (i = 0; i < TCB->size; i++) @@ -41,11 +47,35 @@ void init_Total_Count_Table(int k, Total_Count_Table* TCB) init_Count_Table(&(TCB->sub_h[i])); TCB->sub_h_lock[i].lock = 0; } - - - + TCB->non_unique_k_mer = 0; } + + +void init_Total_Pos_Table(Total_Pos_Table* TCB, Total_Count_Table* pre_TCB) +{ + + TCB->prefix_bits = pre_TCB->prefix_bits; + TCB->suffix_bits = pre_TCB->suffix_bits; + TCB->suffix_mode = pre_TCB->suffix_mode; + TCB->size = pre_TCB->size; + TCB->useful_k_mer = 0; + TCB->total_occ = 0; + TCB->k_mer_index = NULL; + TCB->sub_h_lock = (Hash_table_spin_lock*)malloc(sizeof(Hash_table_spin_lock)*TCB->size); + memset(TCB->sub_h_lock, 0, sizeof(Hash_table_spin_lock)*TCB->size); + TCB->sub_h = (Pos_Table**)malloc(sizeof(Pos_Table*)*TCB->size); + TCB->pos = NULL; + + int i = 0; + for (i = 0; i < TCB->size; i++) + { + init_Pos_Table(&(TCB->sub_h[i])); + TCB->sub_h_lock[i].lock = 0; + } +} + + void destory_Total_Count_Table(Total_Count_Table* TCB) { int i; @@ -58,6 +88,116 @@ void destory_Total_Count_Table(Total_Count_Table* TCB) } +void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k_mer_min_freq, int k_mer_max_freq) +{ + int i; + Count_Table* h; + khint_t k; + uint64_t sub_key; + uint64_t sub_ID; + PCB->useful_k_mer = 0; + PCB->total_occ = 0; + + init_Total_Pos_Table(PCB, TCB); + + khint_t t; ///这就是个迭代器 + int absent; + + for (i = 0; i < TCB->size; i++) + { + h = TCB->sub_h[i]; + for (k = kh_begin(h); k != kh_end(h); ++k) + { + if (kh_exist(h, k)) // test if a bucket contains data + { + ///只有符合频率范围要求的k-mer，才会被加入到pos table中 + if (kh_value(h, k)>=k_mer_min_freq && kh_value(h, k)<=k_mer_max_freq) + { + + sub_ID = i; + sub_key = kh_key(h, k); + + t = kh_put(POS64, PCB->sub_h[sub_ID], sub_key, &absent); + + if (absent) + { + ///kh_value(PCB->sub_h[sub_ID], t) = useful_k_mer + total_occ; + kh_value(PCB->sub_h[sub_ID], t) = PCB->useful_k_mer; + } + else ///哈希表中已有的元素 + { + ///kh_value(PCB->sub_h[sub_ID], t)++; + fprintf(stderr, "ERROR\n"); + } + + + PCB->useful_k_mer++; + PCB->total_occ = PCB->total_occ + kh_value(h, k); + } + } + } + } + + + fprintf(stdout, "useful_k_mer: %lld\n",PCB->useful_k_mer); + fprintf(stdout, "total_occ: %lld\n",PCB->total_occ); + + PCB->k_mer_index = (uint64_t*)malloc(sizeof(uint64_t)*(PCB->useful_k_mer+1)); + + PCB->k_mer_index[0] = 0; + + PCB->total_occ = 0; + PCB->useful_k_mer = 0; + + for (i = 0; i < TCB->size; i++) + { + h = TCB->sub_h[i]; + for (k = kh_begin(h); k != kh_end(h); ++k) + { + if (kh_exist(h, k)) // test if a bucket contains data + { + if (kh_value(h, k)>=k_mer_min_freq && kh_value(h, k)<=k_mer_max_freq) + { + PCB->useful_k_mer++; + PCB->total_occ = PCB->total_occ + kh_value(h, k); + PCB->k_mer_index[PCB->useful_k_mer] = PCB->total_occ; + } + } + } + } + + PCB->pos = (k_mer_pos*)malloc(sizeof(k_mer_pos)*PCB->total_occ); + memset(PCB->pos, 0, sizeof(k_mer_pos)*PCB->total_occ); + + +} + + + + +int cmp_k_mer_pos(const void * a, const void * b) +{ + if ((*(k_mer_pos*)a).readID != (*(k_mer_pos*)b).readID) + { + return (*(k_mer_pos*)a).readID > (*(k_mer_pos*)b).readID ? 1 : -1; + } + else + { + if ((*(k_mer_pos*)a).offset != (*(k_mer_pos*)b).offset) + { + return (*(k_mer_pos*)a).offset > (*(k_mer_pos*)b).offset ? 1 : -1; + } + else + { + return 0; + } + } +} + + + + + /********************************for debug***************************************/ void debug_mode(uint64_t d, uint64_t thread_ID, uint64_t thread_num) { @@ -130,4 +270,4 @@ void test_COUNT64() kh_destroy(COUNT64, h); -} +} \ No newline at end of file diff --git a/Hash_Table.h b/Hash_Table.h index 6f46d5c..3cac6a4 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -4,13 +4,16 @@ #include "kmer.h" KHASH_MAP_INIT_INT64(COUNT64, int) - typedef khash_t(COUNT64) Count_Table; +KHASH_MAP_INIT_INT64(POS64, uint64_t) +typedef khash_t(POS64) Pos_Table; + #define PREFIX_BITS 16 #define MAX_SUFFIX_BITS 64 #define MODE_VALUE 101 + typedef struct { volatile int lock; @@ -25,8 +28,30 @@ typedef struct int suffix_bits; int size; uint64_t suffix_mode; + uint64_t non_unique_k_mer; } Total_Count_Table; +typedef struct +{ + uint64_t offset; + uint64_t readID; +} k_mer_pos; + +typedef struct +{ + Pos_Table** sub_h; + Hash_table_spin_lock* sub_h_lock; + int prefix_bits; + int suffix_bits; + int size; + uint64_t suffix_mode; + k_mer_pos* pos; + uint64_t useful_k_mer; + uint64_t total_occ; + uint64_t* k_mer_index; +} Total_Pos_Table; + + /********************************for debug***************************************/ inline void print_64bit(uint64_t x) { @@ -54,7 +79,9 @@ inline uint64_t mod_d(uint64_t h_key, uint64_t low_key, uint64_t d) -inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, Total_Count_Table* TCB, Hash_code* code, int k) +///inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, Total_Count_Table* TCB, Hash_code* code, int k) +inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, uint64_t suffix_mode, int suffix_bits, +Hash_code* code, int k) { uint64_t h_key, low_key; ///k有可能是64，所以可能会有问题 @@ -73,8 +100,8 @@ inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, Total_Coun ///前一个右移不安全，因为TCB->suffix_bits有可能为64 ///后一个左移安全，因为TCB->suffix_bits不可能为0 //uint64_t sub_ID = (low_key >> TCB->suffix_bits) | (h_key << (64 - TCB->suffix_bits)); - uint64_t sub_ID = (low_key >> SAFE_SHIFT(TCB->suffix_bits)) | (h_key << (64 - TCB->suffix_bits)); - uint64_t sub_key = (low_key & TCB->suffix_mode); + uint64_t sub_ID = (low_key >> SAFE_SHIFT(suffix_bits)) | (h_key << (64 - suffix_bits)); + uint64_t sub_key = (low_key & suffix_mode); *get_sub_ID = sub_ID; *get_sub_key = sub_key; @@ -87,7 +114,7 @@ inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, Total_Coun inline int insert_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k) { uint64_t sub_ID, sub_key; - if(!get_sub_table(&sub_ID, &sub_key, TCB, code, k)) + if(!get_sub_table(&sub_ID, &sub_key, TCB->suffix_mode, TCB->suffix_bits, code, k)) { return 0; } @@ -121,7 +148,7 @@ inline int get_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k) { uint64_t sub_ID, sub_key; - if(!get_sub_table(&sub_ID, &sub_key, TCB, code, k)) + if(!get_sub_table(&sub_ID, &sub_key, TCB->suffix_mode, TCB->suffix_bits, code, k)) { return 0; } @@ -138,18 +165,192 @@ inline int get_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k) } else { - return -1; + return 0; } - return 1; } + + + +inline uint64_t get_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k, uint64_t* r_sub_ID) +{ + + uint64_t sub_ID, sub_key; + if(!get_sub_table(&sub_ID, &sub_key, PCB->suffix_mode, PCB->suffix_bits, code, k)) + { + return (uint64_t)-1; + } + + khint_t t; ///这就是个迭代器 + int absent; + + ///查询哈希表，key为k + t = kh_get(POS64, PCB->sub_h[sub_ID], sub_key); + + if (t != kh_end(PCB->sub_h[sub_ID])) + { + *r_sub_ID = sub_ID; + return kh_value(PCB->sub_h[sub_ID], t); + } + else + { + return (uint64_t)-1; + } + + +} + +inline uint64_t count_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k) +{ + uint64_t sub_ID; + uint64_t ret = get_Total_Pos_Table(PCB, code, k, &sub_ID); + if(ret != (uint64_t)-1) + { + return PCB->k_mer_index[ret + 1] - PCB->k_mer_index[ret]; + } + else + { + return 0; + } +} + + +inline uint64_t locate_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, k_mer_pos** list, int k, uint64_t* r_sub_ID) +{ + uint64_t ret = get_Total_Pos_Table(PCB, code, k, r_sub_ID); + if(ret != (uint64_t)-1) + { + *list = PCB->k_mer_index[ret] + PCB->pos; + return PCB->k_mer_index[ret + 1] - PCB->k_mer_index[ret]; + } + else + { + *list = NULL; + return 0; + } +} + +int cmp_k_mer_pos(const void * a, const void * b); + + +inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k, uint64_t readID, uint64_t pos) +{ + k_mer_pos* list; + int flag = 0; + uint64_t sub_ID; + uint64_t occ = locate_Total_Pos_Table(PCB, code, &list, k, &sub_ID); + + if (occ) + { + + while (__sync_lock_test_and_set(&PCB->sub_h_lock[sub_ID].lock, 1)) + { + while (PCB->sub_h_lock[sub_ID].lock); + } + + if (list[0].offset + 1 < occ) + { + list[0].offset++; + list[list[0].offset].readID = readID; + list[list[0].offset].offset = pos; + } + else + { + list[0].readID = readID; + list[0].offset = pos; + flag = 1; + } + + __sync_lock_release(&PCB->sub_h_lock[sub_ID].lock); + + ///当所有位置都存好后，不会再有其他线程修改该list + ///所以可以在临界区外排序 + if (flag && occ>1) + { + qsort(list, occ, sizeof(k_mer_pos), cmp_k_mer_pos); + } + + + return 1; + } + else + { + return 0; + } + + +} + + +void init_Total_Count_Table(int k, Total_Count_Table* TCB); +void init_Total_Pos_Table(Total_Pos_Table* TCB, Total_Count_Table* pre_TCB); +void destory_Total_Count_Table(Total_Count_Table* TCB); + +void init_Count_Table(Count_Table** table); +void init_Pos_Table(Count_Table** pre_table, Pos_Table** table); + +void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k_mer_min_freq, int k_mer_max_freq); + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + /********************************for debug***************************************/ inline int verify_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k) { uint64_t sub_ID, sub_key; - if(!get_sub_table(&sub_ID, &sub_key, TCB, code, k)) + if(!get_sub_table(&sub_ID, &sub_key, TCB->suffix_mode, TCB->suffix_bits, code, k)) { return 0; } @@ -178,6 +379,11 @@ inline int verify_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int } } + + + + + /********************************for debug***************************************/ inline int Traverse_Total_Count_Table(Total_Count_Table* TCB) { @@ -208,11 +414,6 @@ inline int Traverse_Total_Count_Table(Total_Count_Table* TCB) } -void init_Total_Count_Table(int k, Total_Count_Table* TCB); -void destory_Total_Count_Table(Total_Count_Table* TCB); - -void init_Count_Table(Count_Table** table); - /********************************for debug***************************************/ void test_COUNT64(); diff --git a/Process_Read.cpp b/Process_Read.cpp index 32cc014..c773616 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -9,6 +9,7 @@ gzFile fp; kseq_t *seq; R_buffer RDB; +static uint64_t total_reads; pthread_mutex_t i_readinputMutex; pthread_mutex_t i_queueMutex; @@ -20,6 +21,52 @@ pthread_cond_t i_readinputstallCond; pthread_mutex_t i_doneMutex; +void init_All_reads(All_reads* r) +{ + r->index_size = READ_INIT_NUMBER; + r->index = (uint64_t*)malloc(sizeof(uint64_t)*r->index_size); + r->index[0] = 0; + r->read = NULL; + r->total_reads_bases = 0; + + + r->name_index_size = READ_INIT_NUMBER; + r->name_index = (uint64_t*)malloc(sizeof(uint64_t)*r->name_index_size); + r->name_index[0] = 0; + r->name = NULL; + r->total_name_length = 0; + + r->total_reads = 0; + +} + + +inline void insert_read(All_reads* r, kstring_t* read, kstring_t* name) +{ + r->total_reads++; + r->total_reads_bases = r->total_reads_bases + read->l; + r->total_name_length = r->total_name_length + name->l; + + ///必须要+1 + if (r->index_size < r->total_reads + 2) + { + r->index_size = r->index_size * 2 + 2; + r->index = (uint64_t*)realloc(r->index,sizeof(uint64_t)*(r->index_size)); + + r->name_index_size = r->name_index_size * 2 + 2; + r->name_index = (uint64_t*)realloc(r->name_index,sizeof(uint64_t)*(r->name_index_size)); + } + r->index[r->total_reads] = r->index[r->total_reads-1] + read->l; + r->name_index[r->total_reads] = r->name_index[r->total_reads-1] + name->l; + +} + +void malloc_All_reads(All_reads* r) +{ + r->read = (char*)malloc(sizeof(char)*r->total_reads_bases); + r->name = (char*)malloc(sizeof(char)*r->total_name_length); +} + void init_kseq(char* file) { fp = gzopen(file, "r"); @@ -69,7 +116,11 @@ void init_R_buffer_block(R_buffer_block* curr_sub_block) curr_sub_block->num = 0; } - +void clear_R_buffer() +{ + RDB.all_read_end = 0; + RDB.num = 0; +} void init_R_buffer(int thread_num) { RDB.all_read_end = 0; @@ -111,7 +162,7 @@ void destory_R_buffer() inline void load_read_block(R_buffer_block* read_batch, int batch_read_size, - int* return_file_flag) + int* return_file_flag, int is_insert) { int inner_i = 0; int file_flag = 1; @@ -126,6 +177,14 @@ inline void load_read_block(R_buffer_block* read_batch, int batch_read_size, if (file_flag == 1) { + read_batch->read[inner_i].ID = total_reads; + total_reads++; + + if (is_insert) + { + insert_read(&R_INF, &read_batch->read[inner_i].seq, &read_batch->read[inner_i].name); + } + inner_i++; } else if (file_flag == 0) @@ -184,8 +243,13 @@ inline void pop_R_block(R_buffer_block* curr_sub_block) } -void* input_reads_muti_threads(void*) + +void* input_reads_muti_threads(void* arg) { + int is_insert = *((int*)arg); + + + total_reads = 0; int i = 0; @@ -202,7 +266,7 @@ void* input_reads_muti_threads(void*) - load_read_block(&tmp_buf, RDB.block_inner_size, &file_flag); + load_read_block(&tmp_buf, RDB.block_inner_size, &file_flag, is_insert); if (file_flag == 0) { @@ -233,6 +297,13 @@ void* input_reads_muti_threads(void*) destory_R_buffer_block(&tmp_buf); + fprintf(stdout, "total_reads: %llu\n",total_reads); + ///fprintf(stdout, "R_INF.total_reads: %llu\n",R_INF.total_reads); + ///fprintf(stdout, "R_INF.index[R_INF.total_reads]: %llu\n",R_INF.index[R_INF.total_reads]); + fprintf(stdout, "R_INF.total_reads_bases: %llu\n",R_INF.total_reads_bases); + ///fprintf(stdout, "R_INF.name_index[R_INF.total_reads]: %llu\n",R_INF.name_index[R_INF.total_reads]); + fprintf(stdout, "R_INF.total_name_length: %llu\n",R_INF.total_name_length); + } @@ -293,7 +364,7 @@ void Counting_block() load_read_block(&tmp_buf, RDB.block_inner_size, - &file_flag); + &file_flag, 0); if (file_flag == 0) diff --git a/Process_Read.h b/Process_Read.h index 4734264..a02d75b 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -7,11 +7,17 @@ #include #include "kseq.h" +#define READ_INIT_NUMBER 1000 + #define READ_BLOCK_SIZE 64 #define READ_BLOCK_NUM_PRE_THR 100 #define IS_FULL(buffer) ((buffer.num >= buffer.size)?1:0) #define IS_EMPTY(buffer) ((buffer.num == 0)?1:0) +#define Get_READ_LENGTH(R_INF, ID) (R_INF.index[ID+1] - R_INF.index[ID]) +#define Get_NAME_LENGTH(R_INF, ID) (R_INF.name_index[ID+1] - R_INF.name_index[ID]) +#define Get_READ(R_INF, ID) (R_INF.read + R_INF.index[ID]) +#define Get_NAME(R_INF, ID) (R_INF.name + R_INF.name_index[ID]) @@ -23,6 +29,24 @@ int get_read(kseq_t *s); +typedef struct +{ + char* read; + char* name; + uint64_t* index; + uint64_t index_size; + uint64_t* name_index; + uint64_t name_index_size; + uint64_t total_reads; + uint64_t total_reads_bases; + uint64_t total_name_length; + +} All_reads; + +extern All_reads R_INF; + +void malloc_All_reads(All_reads* r); + typedef struct { kseq_t* read; @@ -41,6 +65,7 @@ typedef struct } R_buffer; void init_R_buffer(int thread_num); +void init_All_reads(All_reads* r); void* input_reads_muti_threads(void*); void init_R_buffer_block(R_buffer_block* curr_sub_block); int get_reads_mul_thread(R_buffer_block* curr_sub_block); @@ -48,5 +73,7 @@ int get_reads_mul_thread(R_buffer_block* curr_sub_block); void Counting_block(); void destory_R_buffer_block(R_buffer_block* curr_sub_block); void destory_R_buffer(); +void clear_R_buffer(); + #endif diff --git a/kseq.h b/kseq.h index 2f94a64..6dea15b 100644 --- a/kseq.h +++ b/kseq.h @@ -222,6 +222,7 @@ typedef struct __kstring_t { kstring_t name, comment, seq, qual; \ int last_char; \ kstream_t *f; \ + uint64_t ID; \ } kseq_t; #define KSEQ_INIT2(SCOPE, type_t, __read) \ diff --git a/main.cpp b/main.cpp index 7307e84..971d529 100644 --- a/main.cpp +++ b/main.cpp @@ -21,14 +21,18 @@ int main(int argc, char *argv[]) if (!CommandLine_process(argc, argv)) return 1; - init_kseq(read_file_name); + ///init_kseq(read_file_name); fprintf(stdout, "k-mer length: %d\n",k_mer_length); ///Counting(); Counting_multiple_thr(); - destory_kseq(); + Build_hash_table_multiple_thr(); + + ///verify_Position_hash_table(); + + ///destory_kseq(); ///debug_Counting();