Update position hash table

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; i<thread_num; i++)
+		pthread_join(_r_threads[i], NULL);
+
+    free(_r_threads);
+
+    destory_kseq();
+
+    fprintf(stdout, "Finish Building hash table ...... \n");
+
+    fprintf(stdout, "%-30s%18.2f\n\n", "Build hash table time:", Get_T() - start_time);
+
+    ///destory_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");
+
+
+}
+
+
+
+/********************************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<count2; j++)
+                        {
+                            if(list[j].readID < list[j - 1].readID)
+                            {
+                                fprintf(stderr, "locate error\n");
+                            }
+                            else if(list[j].readID == list[j - 1].readID)
+                            {
+                                if(list[j].offset < list[j - 1].offset)
+                                {
+                                    fprintf(stderr, "locate error\n");
+                                }
+                            }
+                        }
+
+                        if(count2 != 0)
+                        {
+                            for(j=0; j<count2; j++)
+                            {   
+                                if(list[j].readID == read_number && list[j].offset == HPC_base - k_mer_length + 1)
+                                {
+                                    break;
+                                }
+                            }
+
+                            if(j == count2)
+                            {
+                                fprintf(stderr, "locate error, read_number: %llu, pos: %llu\n", 
+                                read_number, HPC_base - k_mer_length + 1);
+
+                                for(j=0; j<count2; j++)
+                                {   
+
+                                    fprintf(stderr, "j: %llu, readID: %llu, offset: %llu\n", j, list[j].readID, list[j].offset);
+                            
+                                }
+                            }
+                        }
+                        
+                        
+                    }
+                    else
+                    {
+                        fprintf(stderr, "locate error\n");
+                    }
+
+                }
+                
+            }
+            else
+            {
+                avalible_k = 0;
+                init_Hash_code(&k_code);
+            }
+
+            HPC_base++;
+            
+        }
+
+        read_number++;
+    }
+
+    destory_kseq();
+
+    fprintf(stdout, "Finish Verifying Position Table!\n");
 
 }
 
@@ -235,17 +541,6 @@ void Counting_multiple_thr()
 
 
 
-
-
-
-
-
-
-
-
-
-
-
 
 
 

Assembly.h

@@ -4,7 +4,11 @@
 
 
 void Counting_multiple_thr();
+void Build_hash_table_multiple_thr();
+
 /********************************for debug***************************************/
 void Verify_Counting();
 
+/********************************for debug*****************************************/
+void verify_Position_hash_table();
 #endif

CommandLines.cpp

@@ -2,11 +2,22 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include "ketopt.h"
+#include <sys/time.h>
+
 
 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()
 {

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

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<<TCB->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);
     
 
-}
+}

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();
 

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)

Process_Read.h

@@ -7,11 +7,17 @@
 #include <zlib.h>
 #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

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)		\

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();