removed hifiasm hash table index

Assembly.cpp

@@ -11,426 +11,13 @@
 #include "Output.h"
 #include "htab.h"
 
-Total_Count_Table TCB;
+void *ha_flt_tab;
-Total_Pos_Table PCB;
+void *ha_idx;
+
 All_reads R_INF;
 pthread_mutex_t statistics;
 
-void* Perform_Counting(void* arg)
+void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length, char* new_read, int* new_length)
-{
-    int i = 0;
-    HPC_seq HPC_read;
-
-    R_buffer_block curr_sub_block;
-
-    init_R_buffer_block(&curr_sub_block);
-
-    long long read_number = 0;
-    long long select_k_mer_number = 0 ;
-    long long k_mer_number = 0 ;
-
-	int file_flag = 1;
-
-    uint64_t code;
-
-    uint64_t end_pos;
-
-    Hash_code k_code;
-
-    int avalible_k = 0;
-
-
-	while (file_flag != 0)
-	{
-
-        file_flag = get_reads_mul_thread(&curr_sub_block);
-
-        read_number = read_number + curr_sub_block.num; 
-
-        for (i = 0; i < curr_sub_block.num; i++)
-        {
-            ///forward strand
-            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;
-
-            while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
-            {
-                if(code < 4)
-                {
-                    k_mer_append(&k_code,code, asm_opt.k_mer_length);
-                    avalible_k++;
-                    if (avalible_k >= asm_opt.k_mer_length)
-                    {
-                        if(insert_Total_Count_Table(&TCB, &k_code, asm_opt.k_mer_length))
-                        {
-                            select_k_mer_number++;
-                        }
-
-                        k_mer_number++;
-                    }
-                }
-                else
-                {
-                    avalible_k = 0;
-                    init_Hash_code(&k_code);
-                }
-            }
-        }
-    }
-
-    destory_R_buffer_block(&curr_sub_block);
-    free(arg);
-
-    return NULL;
-}
-
-void* Perform_Counting_non_first(void* arg)
-{
-    int thr_ID = *((int*)arg);
-
-    uint64_t i = 0;
-    HPC_seq HPC_read;
-
-    long long select_k_mer_number = 0 ;
-    long long k_mer_number = 0 ;
-
-    uint64_t code;
-
-    uint64_t end_pos;
-
-    Hash_code k_code;
-
-    int avalible_k = 0;
-
-    UC_Read g_read;
-    init_UC_Read(&g_read);
-
-    for (i = thr_ID; i < R_INF.total_reads; i = i + asm_opt.thread_num)
-    {
-        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, asm_opt.k_mer_length);
-                avalible_k++;
-                if (avalible_k >= asm_opt.k_mer_length)
-                {
-                    if(insert_Total_Count_Table(&TCB, &k_code, asm_opt.k_mer_length))
-                    {
-                        select_k_mer_number++;
-                    }
-                    k_mer_number++;
-                }
-            }
-            else
-            {
-                avalible_k = 0;
-                init_Hash_code(&k_code);
-            }
-        }
-    }
-    
-    destory_UC_Read(&g_read);
-    free(arg);
-
-    return NULL;
-}
-
-void* Build_hash_table_non_first(void* arg)
-{
-    int thr_ID = *((int*)arg);
-
-    uint64_t i = 0;
-    HPC_seq HPC_read;
-    uint64_t code;
-    uint64_t end_pos;
-
-    ///long long HPC_base;
-
-    Hash_code k_code;
-
-    int avalible_k = 0;
-
-    UC_Read g_read;
-    init_UC_Read(&g_read);
-
-    for (i = thr_ID; i < R_INF.total_reads; i = i + asm_opt.thread_num)
-    {
-        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;
-
-        ///HPC_base = 0;
-
-        while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
-        {
-            if(code < 4)
-            {
-                k_mer_append(&k_code,code, asm_opt.k_mer_length);
-                avalible_k++;
-                if (avalible_k >= asm_opt.k_mer_length)
-                {
-                    ///there are two requirements
-                    ///1. hash(k-mer) 
-                    ///2. occ(k-mer) 
-                    ///TCB just meet the first requirement，while PCB needs to meet both of them
-                    insert_Total_Pos_Table(&PCB, &k_code, asm_opt.k_mer_length, i, end_pos);
-                }
-            }
-            else
-            {
-                avalible_k = 0;
-                init_Hash_code(&k_code);
-            }
-        }
-    }
-
-    destory_UC_Read(&g_read);
-    free(arg);
-
-    return NULL;
-}
-
-
-void* Build_hash_table(void* 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;
-    uint64_t end_pos;
-
-    ///long long HPC_base;
-
-    Hash_code k_code;
-
-    int avalible_k = 0;
-
-	while (file_flag != 0)
-	{
-
-        file_flag = get_reads_mul_thread(&curr_sub_block);
-
-
-        for (i = 0; i < curr_sub_block.num; i++)
-        {
- 
-            ///forward strand
-            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, &end_pos)) != 6)
-            {
-                if(code < 4)
-                {
-                    k_mer_append(&k_code,code, asm_opt.k_mer_length);
-                    avalible_k++;
-                    if (avalible_k >= asm_opt.k_mer_length)
-                    {
-                       ///there are two requirements
-                       ///1. hash(k-mer) 
-                       ///2. occ(k-mer) 
-                       ///TCB just meet the first requirement，while PCB needs to meet both of them
-                        insert_Total_Pos_Table(&PCB, &k_code, asm_opt.k_mer_length, curr_sub_block.read[i].ID, end_pos);
-                    }
-                }
-                else
-                {
-                    avalible_k = 0;
-                    init_Hash_code(&k_code);
-                }
-            }
-
-            ///load read
-            ha_compress_base(Get_READ(R_INF, curr_sub_block.read[i].ID),
-            curr_sub_block.read[i].seq.s, curr_sub_block.read[i].seq.l, 
-            &R_INF.N_site[curr_sub_block.read[i].ID], HPC_read.N_occ);
-            
-            memcpy(R_INF.name+R_INF.name_index[curr_sub_block.read[i].ID], 
-            curr_sub_block.read[i].name.s, curr_sub_block.read[i].name.l);
-        }
-    }
-
-    destory_R_buffer_block(&curr_sub_block);
-    free(arg);
-
-    return NULL;
-}
-
-
-void Counting_multiple_thr()
-{
-
-    double start_time = Get_T();
-
-    fprintf(stderr, "Begin Counting... \n");
-
-    init_Total_Count_Table(asm_opt.k_mer_length, &TCB);
-
-    pthread_t inputReadsHandle;
-
-    int *is_insert = (int*)malloc(sizeof(*is_insert));
-	*is_insert = 1;
-
-    if (asm_opt.roundID == 0)
-    {
-        init_gz_files(&asm_opt);
-        init_All_reads(&R_INF);
-        init_R_buffer(asm_opt.thread_num);
-        pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert);
-    }
-    
-
-    pthread_t *_r_threads;
-
-	_r_threads = (pthread_t *)malloc(sizeof(pthread_t)*asm_opt.thread_num);
-
-    int i = 0;
-
-	for (i = 0; i < asm_opt.thread_num; i++)
-	{
-		int *arg = (int*)malloc(sizeof(*arg));
-		*arg = i;
-        if (asm_opt.roundID == 0)
-        {
-            pthread_create(_r_threads + i, NULL, Perform_Counting, (void*)arg);
-        }
-        else
-        {
-            pthread_create(_r_threads + i, NULL, Perform_Counting_non_first, (void*)arg);
-        }
-	}
-
-
-    for (i = 0; i < asm_opt.thread_num; i++)
-		pthread_join(_r_threads[i], NULL);
-
-    free(_r_threads);
-
-    ///destory_R_buffer();
-    
-    ///destory_Total_Count_Table(&TCB);
-    if (asm_opt.roundID == 0)
-    {
-        pthread_join(inputReadsHandle, NULL);
-        ///destory_kseq();
-        destory_gz_files();
-    }
-
-    fprintf(stderr, "Counting has been completed.\n");
-
-    fprintf(stderr, "%-30s%18.2f\n\n", "Counting time:", Get_T() - start_time);
-   
-    free(is_insert);
-
-}
-
-
-void Build_hash_table_multiple_thr()
-{
-    double start_time = Get_T();
-
-    fprintf(stderr, "Begin building hash table... \n");
-
-    init_Total_Pos_Table(&PCB, &TCB);
-
-    double T_start_time = Get_T();
-
-    Traverse_Counting_Table(&TCB, &PCB, asm_opt.k_mer_min_freq, asm_opt.k_mer_max_freq);
-
-    fprintf(stderr, "%-30s%18.2f\n\n", "Traverse time:", Get_T() - T_start_time);
-   
-    ///at this moment, TCB can be free
-    destory_Total_Count_Table(&TCB);
-
-    pthread_t inputReadsHandle;
-
-    int *is_insert = (int*)malloc(sizeof(*is_insert));
-	*is_insert = 0;
-
-    if (asm_opt.roundID == 0)
-    {
-        init_gz_files(&asm_opt);
-        clear_R_buffer();
-        malloc_All_reads(&R_INF);
-        pthread_create(&inputReadsHandle, NULL, input_reads_muti_threads, (void*)is_insert);
-    }
-
-
-
-    pthread_t *_r_threads;
-
-	_r_threads = (pthread_t *)malloc(sizeof(pthread_t) * asm_opt.thread_num);
-
-    int i = 0;
-
-	for (i = 0; i < asm_opt.thread_num; i++)
-	{
-		int *arg = (int*)malloc(sizeof(*arg));
-		*arg = i;
-
-        if (asm_opt.roundID == 0)
-        {
-            pthread_create(_r_threads + i, NULL, Build_hash_table, (void*)arg);
-        }
-        else
-        {
-            pthread_create(_r_threads + i, NULL, Build_hash_table_non_first, (void*)arg);
-        }
-	}
-
-    if (asm_opt.roundID == 0)
-    {
-        pthread_join(inputReadsHandle, NULL);
-    }
-
-    for (i = 0; i < asm_opt.thread_num; i++)
-		pthread_join(_r_threads[i], NULL);
-
-    free(_r_threads);
-
-    
-
-    fprintf(stderr, "Hash table has been built.\n");
-
-    fprintf(stderr, "%-30s%18.2f\n\n", "Build hash table time:", Get_T() - start_time);
-
-    if (asm_opt.roundID == 0)
-    {
-        ///destory_kseq();
-        destory_gz_files();
-        destory_R_buffer();
-    }
-    
-    free(is_insert);
-}
-
-
-void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length,
-char* new_read, int* new_length)
 {
     int i, j;
     int pre_i, new_i;
@@ -538,8 +125,7 @@ inline int get_cigar_errors(Cigar_record* cigar)
     return total_errors;
 }
 
-int debug_cigar(Cigar_record* cigar, char* pre_read, int pre_length,
+int debug_cigar(Cigar_record* cigar, char* pre_read, int pre_length, char* new_read, int new_length, int correct_base)
-char* new_read, int new_length, int correct_base)
 {
     int i;
     int total_errors = 0;
@@ -556,13 +142,10 @@ char* new_read, int new_length, int correct_base)
         fprintf(stderr, "total_errors: %d, correct_base: %d\n", total_errors, correct_base);
     }
 
-    
-
     int pre_i, new_i;
     int operation, operation_length;
     pre_i = new_i = 0;
 
-
     for (i = 0; i < (long long)cigar->length; i++)
     {
         operation = Get_Cigar_Type(cigar->record[i]);
@@ -592,8 +175,6 @@ char* new_read, int new_length, int correct_base)
 
     }
 
-
-
     if (pre_i != pre_length)
     {
         fprintf(stderr, "pre_i: %d, pre_length: %d\n", pre_i, pre_length);
@@ -607,7 +188,6 @@ char* new_read, int new_length, int correct_base)
 
     return 1;
 
-
     char* tmp_seq = (char*)malloc(new_length + pre_length);
     int tmp_length;
 
@@ -624,9 +204,6 @@ char* new_read, int new_length, int correct_base)
     }
 
 
-
-
-    
     get_uncorrected_read_from_cigar(cigar, new_read, new_length, tmp_seq, &tmp_length);
 
 
@@ -652,7 +229,6 @@ char* new_read, int new_length, int correct_base)
 
 inline void push_cigar(Compressed_Cigar_record* records, long long ID, Cigar_record* input)
 {
-
     if (input->length > records[ID].size)
     {
         records[ID].size = input->length;
@@ -673,8 +249,7 @@ inline void push_cigar(Compressed_Cigar_record* records, long long ID, Cigar_rec
 }
 
 
-void push_overlaps(ma_hit_t_alloc* paf, overlap_region_alloc* overlap_list, int flag, 
+void push_overlaps(ma_hit_t_alloc* paf, overlap_region_alloc* overlap_list, int flag, All_reads* R_INF, int if_reverse)
-All_reads* R_INF, int if_reverse)
 {
     long long i = 0, xLen, yLen;
     ma_hit_t tmp;
@@ -717,12 +292,9 @@ All_reads* R_INF, int if_reverse)
             add_ma_hit_t_alloc(paf, &tmp);
         }
     }
-    
 }
 
-
+int if_exact_match(char* x, long long xLen, char* y, long long yLen, long long xBeg, long long xEnd, long long yBeg, long long yEnd)
-int if_exact_match(char* x, long long xLen, char* y, long long yLen, 
-long long xBeg, long long xEnd, long long yBeg, long long yEnd)
 {
     long long overlapLen = xEnd - xBeg + 1;
 
@@ -747,8 +319,7 @@ long long xBeg, long long xEnd, long long yBeg, long long yEnd)
     return 0;
 }
 
-long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_list, 
+long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_list, overlap_region_alloc* overlap_list, int flag)
-overlap_region_alloc* overlap_list, int flag)
 {
     long long i = 0;
     long long available_overlaps = 0;
@@ -835,7 +406,7 @@ void get_new_candidates(long long readID, UC_Read* g_read, overlap_region_alloc*
             avalible_k++;
             if (avalible_k >= asm_opt.k_mer_length)
             {
-                list_length = locate_Total_Pos_Table(&PCB, &k_code, &list, asm_opt.k_mer_length, &sub_ID);
+                //list_length = locate_Total_Pos_Table(&PCB, &k_code, &list, asm_opt.k_mer_length, &sub_ID);
 
                 if (list_length != 0)
                 {
@@ -1246,7 +817,8 @@ void Overlap_calculate_multipe_thr()
     }
 
 
-    destory_Total_Pos_Table(&PCB);
+	ha_idx_destroy(ha_idx);
+	ha_idx = 0;
 
     fprintf(stderr, "All overlaps have been calculated.\n");
 
@@ -1813,15 +1385,13 @@ long long readNum, long long rescue_threshold, float cluster_threshold)
 }
 
 
-
 void generate_overlaps(int last_round)
 {
     double start_time = Get_T();
     asm_opt.roundID = asm_opt.number_of_round - last_round;
     fprintf(stderr, "Begin calculting final overlaps ...\n");
 
-    Counting_multiple_thr();
+	ha_idx = ha_gen_mzidx(&asm_opt, ha_flt_tab, asm_opt.roundID == 0? 0 : 1, &R_INF);
-    Build_hash_table_multiple_thr();
 
     pthread_t *_r_threads;
 
@@ -1842,7 +1412,8 @@ void generate_overlaps(int last_round)
     free(_r_threads);
 
     ///rescue_edges(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads, 4, 0.985);
-    destory_Total_Pos_Table(&PCB);
+	ha_idx_destroy(ha_idx);
+	ha_idx = 0;
 
     fprintf(stderr, "Final overlaps have been calculated.\n");
     fprintf(stderr, "%-30s%18.2f\n\n", "Final overlaps calculation time:", Get_T() - start_time); 
@@ -1860,14 +1431,12 @@ void generate_overlaps(int last_round)
 
 void Correct_Reads(int last_round)
 {
-    
+    if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name))
-    if(asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, 
-    asm_opt.output_file_name))
     {
         build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf, 
         R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round, 
         asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, 0);
-        exit(1);
+        exit(0);
     }
     else
     {
@@ -1884,9 +1453,9 @@ void Correct_Reads(int last_round)
 
     fprintf(stderr, "Error correction: Start the %d-th round ...\n", asm_opt.roundID);
 
-    Counting_multiple_thr();
+	ha_idx = ha_gen_mzidx(&asm_opt, ha_flt_tab, asm_opt.roundID == 0? 0 : 1, &R_INF);
-    Build_hash_table_multiple_thr();
     Overlap_calculate_multipe_thr();
+	ha_idx_destroy(ha_idx);
 
     fprintf(stderr, "Error correction: The %d-th round has been completed.\n", asm_opt.roundID);
 

CommandLines.h

@@ -53,7 +53,7 @@ extern hifiasm_opt_t asm_opt;
 void init_opt(hifiasm_opt_t* asm_opt);
 void destory_opt(hifiasm_opt_t* asm_opt);
 void clear_opt(hifiasm_opt_t* asm_opt, int last_round);
-int CommandLine_process (int argc, char *argv[], hifiasm_opt_t* asm_opt);
+int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt);
 double Get_T(void);
 
 #endif

Hash_Table.cpp

@@ -779,195 +779,6 @@ void init_Pos_Table(Pos_Table** table)
     *table = ha_pt_init();
 }
 
-void init_Total_Count_Table(int k, Total_Count_Table* TCB)
-{
-    if(k>64)
-    {
-        fprintf(stderr, "k-mer is too long. The length of k-mer must <= 64.");
-        fflush(stderr);
-        exit(0);
-    }
-
-    int total_bits = k * 2; 
-    TCB->prefix_bits = PREFIX_BITS;
-    TCB->suffix_bits = total_bits - TCB->prefix_bits;
-    if (TCB->suffix_bits > MAX_SUFFIX_BITS)
-    {
-        TCB->suffix_bits = MAX_SUFFIX_BITS;
-        TCB->prefix_bits = total_bits - TCB->suffix_bits;
-    }
-    ///TCB->suffix_mode = (1ULL<<TCB->suffix_bits) - 1;
-    ///right shift is safe, since TCB->suffix_bits cannot be 0
-    TCB->suffix_mode = ALL >> (64 - TCB->suffix_bits);
-    
-    ///number of small hash table
-    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++)
-    {
-        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;
-    for (i = 0; i < TCB->size; i++)
-    {
-        ha_ct_destroy(TCB->sub_h[i]);
-    }
-    free(TCB->sub_h);
-    free(TCB->sub_h_lock);
-}
-
-
-void destory_Total_Pos_Table(Total_Pos_Table* TCB)
-{
-    free(TCB->k_mer_index);
-    free(TCB->sub_h_lock);
-    free(TCB->pos);
-
-    int i;
-    for (i = 0; i < TCB->size; i++)
-    {
-        ha_pt_destroy(TCB->sub_h[i]);
-    }
-    free(TCB->sub_h);
-}
-
-/*
-void write_Total_Pos_Table(Total_Pos_Table* TCB, char* read_file_name)
-{
-    fprintf(stderr, "Writing index to disk... \n");
-    char* index_name = (char*)malloc(strlen(read_file_name)+5);
-    sprintf(index_name, "%s.idx", read_file_name);
-    FILE* fp = fopen(index_name, "w");
-    fwrite(&asm_opt.adapterLen, sizeof(asm_opt.adapterLen), 1, fp);
-    fwrite(&asm_opt.k_mer_min_freq, sizeof(asm_opt.k_mer_min_freq), 1, fp);
-    fwrite(&asm_opt.k_mer_max_freq, sizeof(asm_opt.k_mer_max_freq), 1, fp);
-    fwrite(&TCB->prefix_bits, sizeof(TCB->prefix_bits), 1, fp);
-    fwrite(&TCB->suffix_bits, sizeof(TCB->suffix_bits), 1, fp);
-    fwrite(&TCB->suffix_mode, sizeof(TCB->suffix_mode), 1, fp);
-    fwrite(&TCB->size, sizeof(TCB->size), 1, fp);
-    fwrite(&TCB->useful_k_mer, sizeof(TCB->useful_k_mer), 1, fp);
-    fwrite(&TCB->total_occ, sizeof(TCB->total_occ), 1, fp);
-    fwrite(TCB->k_mer_index, sizeof(uint64_t), TCB->useful_k_mer+1, fp);
-    fwrite(TCB->pos, sizeof(k_mer_pos), TCB->total_occ, fp);
-
-
-    int i;
-    for (i = 0; i < TCB->size; i++)
-    {
-        kh_write(POS64, TCB->sub_h[i], fp);
-    }
-
-    free(index_name);    
-    fclose(fp);
-    fprintf(stderr, "Index has been written.\n");
-}
-
-
-int load_Total_Pos_Table(Total_Pos_Table* TCB, char* read_file_name)
-{
-    fprintf(stderr, "Loading index from disk... \n");
-    char* index_name = (char*)malloc(strlen(read_file_name)+5);
-    sprintf(index_name, "%s.idx", read_file_name);
-    FILE* fp = fopen(index_name, "r");
-    if (!fp)
-    {
-        return 0;
-    }
-    int f_flag;
-    int local_adapterLen;
-    f_flag = fread(&local_adapterLen, sizeof(local_adapterLen), 1, fp);
-    if(local_adapterLen != asm_opt.adapterLen)
-    {
-        fprintf(stderr, "the adapterLen of index is: %d, but the adapterLen set by user is: %d\n", 
-        local_adapterLen, asm_opt.adapterLen);
-        exit(1);
-    }
-    f_flag += fread(&asm_opt.k_mer_min_freq, sizeof(asm_opt.k_mer_min_freq), 1, fp);
-    f_flag += fread(&asm_opt.k_mer_max_freq, sizeof(asm_opt.k_mer_max_freq), 1, fp);
-    f_flag += fread(&TCB->prefix_bits, sizeof(TCB->prefix_bits), 1, fp);
-    f_flag += fread(&TCB->suffix_bits, sizeof(TCB->suffix_bits), 1, fp);
-    f_flag += fread(&TCB->suffix_mode, sizeof(TCB->suffix_mode), 1, fp);
-    f_flag += fread(&TCB->size, sizeof(TCB->size), 1, fp);
-    f_flag += fread(&TCB->useful_k_mer, sizeof(TCB->useful_k_mer), 1, fp);
-    f_flag += fread(&TCB->total_occ, sizeof(TCB->total_occ), 1, fp);
-
-    if (TCB->useful_k_mer+1)
-    {
-        TCB->k_mer_index = (uint64_t*)malloc(sizeof(uint64_t)*(TCB->useful_k_mer+1));
-        f_flag += fread(TCB->k_mer_index, sizeof(uint64_t), TCB->useful_k_mer+1, fp);
-    }
-    else
-    {
-        TCB->k_mer_index = NULL;
-    }
-    
-    
-    if (TCB->total_occ)
-    {
-        TCB->pos = (k_mer_pos*)malloc(sizeof(k_mer_pos)*TCB->total_occ);
-        f_flag += fread(TCB->pos, sizeof(k_mer_pos), TCB->total_occ, fp);
-    }
-    else
-    {
-        TCB->pos = 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);
-
-    int i;
-    for (i = 0; i < TCB->size; i++)
-    {
-        init_Pos_Table(&(TCB->sub_h[i]));
-        TCB->sub_h_lock[i].lock = 0;
-        kh_load(POS64, TCB->sub_h[i], fp);
-    }
-
-    free(index_name);    
-    fclose(fp);
-    fprintf(stderr, "Index has been loaded.\n");
-
-    return 1;
-}
-*/
 
 typedef struct
 {
@@ -988,242 +799,6 @@ void insert_H_peaks(H_peaks* h, long long index, long long value)
     h->list[index] += value;
 }
 
-///1: a > b; -1: a < b; 0: a=b
-int cmp_Hash_code(Hash_code* a, Hash_code* b)
-{
-    if(a->x[1] > b->x[1])
-    {
-        return 1;
-    }
-    if(a->x[1] < b->x[1])
-    {
-        return -1;
-    }
-    ///a->x[1] == b->x[1]
-    if(a->x[0] > b->x[0])
-    {
-        return 1;
-    }
-    if(a->x[0] < b->x[0])
-    {
-        return -1;
-    }
-
-    return 0;
-}
-
-void get_total_freq(Total_Count_Table* TCB, uint64_t sub_ID, uint64_t sub_key, long long* T_count)
-{
-	khint_t t = ha_ct_get(TCB->sub_h[sub_ID], sub_key);
-	*T_count = kh_val(TCB->sub_h[sub_ID], t);
-}
-
-void get_peak(Total_Count_Table* TCB, long long* min, long long* max, long long* up_boundary)
-{
-    int i;
-    Count_Table* h;
-    khint_t k;
-    long long count;
-    H_peaks LH;
-    LH.list = NULL;
-    LH.length = 0;
-    uint64_t sub_ID;
-    uint64_t sub_key;
-    
-    for (i = 0; i < TCB->size; i++)
-    {
-        h = TCB->sub_h[i];
-        for (k = 0; k != kh_end(h); ++k)
-        {
-            if (kh_exist(h, k))            // test if a bucket contains data
-            {
-                sub_ID = i;
-                sub_key = kh_key(h, k);
-
-                get_total_freq(TCB, sub_ID, sub_key, &count);
-                insert_H_peaks(&LH, count, count);
-            }
-        }
-    }
-
-    (*max) = -1;
-    (*min) = -1;
-    long long max_value = -1;
-    //// seed with freq 1 is useless
-    for (i = 2; i < (long long)LH.length; i++)
-    {
-        if(LH.list[i] >= max_value)
-        {
-            max_value = LH.list[i];
-            (*max) = i;
-        }
-    }
-
-    long long opt = 4;
-    (*up_boundary) = -1;
-    for (i = (*max) + opt; i < (long long)LH.length; i++)
-    {
-        if(LH.list[i] > LH.list[i-opt])
-        {
-            long long j = i-opt;
-            for (; j < i; j++)
-            {
-                if(LH.list[j] < LH.list[j+1])
-                {
-                    (*up_boundary) = j;
-                    goto end_opt;
-                }
-            }
-            
-            (*up_boundary) = i;
-            goto end_opt;
-        }
-    }
-
-    end_opt:
-    if((*up_boundary) == -1 || (*up_boundary) > (*max) * 10)
-    {
-        (*up_boundary) = (*max) * 10;
-    }
-
-    long long min_value = max_value;
-    //// seed with freq 1 is useless
-    for (i = 2; i < (long long)LH.length && i < (*max); i++)
-    {
-        if(LH.list[i] < min_value && LH.list[i] != 0)
-        {
-            min_value = LH.list[i];
-            (*min) = i;
-        }
-    }
-
-    free(LH.list);
-}
-
-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;
-
-    long long freq_min, max, freq_up;
-    ///get_peak_debug(TCB, &freq_min, &freq_max);
-    get_peak(TCB, &freq_min, &max, &freq_up);
-    // fprintf(stdout, "freq_min: %d, freq_max: %d, freq_up:%d\n", 
-    // freq_min, max, freq_up);
-    if(freq_min < k_mer_min_freq)
-    {
-        k_mer_min_freq = freq_min;
-    }
-    if(freq_up > k_mer_max_freq)
-    {
-        k_mer_max_freq = freq_up;
-    }
-
-    // fprintf(stdout, "k_mer_min_freq: %d, k_mer_max_freq: %d\n", 
-    // k_mer_min_freq, k_mer_max_freq);
-
-
-    khint_t t;  
-    int absent;
-    long long count;
-
-    /********************************************
-     hash_table(key) ----> PCB->k_mer_index ------> PCB->pos
-     ********************************************/
-    for (i = 0; i < TCB->size; i++)
-    {
-        h = TCB->sub_h[i];
-        for (k = 0; k != kh_end(h); ++k)
-        {
-            if (kh_exist(h, k))            // test if a bucket contains data
-            {
-                sub_ID = i;
-                sub_key = kh_key(h, k);
-                get_total_freq(TCB, sub_ID, sub_key, &count);
-
-                if (count>=k_mer_min_freq && count<=k_mer_max_freq)
-                {
-                    t = ha_pt_put(PCB->sub_h[sub_ID], sub_key, &absent);
-    
-                    if (absent)
-                    {
-                        ///kh_val(PCB->sub_h[sub_ID], t) = useful_k_mer + total_occ;
-                        kh_val(PCB->sub_h[sub_ID], t) = PCB->useful_k_mer;
-                    }
-                    else   
-                    {
-                        ///kh_val(PCB->sub_h[sub_ID], t)++;
-                        fprintf(stderr, "ERROR\n");
-                    }
-                   
-                    PCB->useful_k_mer++;
-                    PCB->total_occ = PCB->total_occ + kh_val(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 = 0; k != kh_end(h); ++k)
-        {
-            if (kh_exist(h, k))            // test if a bucket contains data
-            {
-                sub_ID = i;
-                sub_key = kh_key(h, k);
-                get_total_freq(TCB, sub_ID, sub_key, &count);
-
-                ///if (kh_val(h, k)>=k_mer_min_freq && kh_val(h, k)<=k_mer_max_freq)
-                if (count>=k_mer_min_freq && count<=k_mer_max_freq)
-                {
-                    PCB->useful_k_mer++;
-                    PCB->total_occ = PCB->total_occ + kh_val(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;
-        }
-    }
-}
-
-
 void init_Chain_Data(Chain_Data* x)
 {
     x->length = 0;

Hash_Table.h

@@ -33,23 +33,6 @@ KHASHL_MAP_INIT(static inline, Pos_Table, ha_pt, uint64_t, uint64_t, kh_hash_dum
 ///#define CIGAR_MAX_LENGTH THRESHOLD*2+2
 #define CIGAR_MAX_LENGTH 31*2+4
 
-typedef struct
-{
-    volatile int lock;
-} Hash_table_spin_lock;
-
-typedef struct
-{
-	Count_Table** sub_h;
-    Hash_table_spin_lock* sub_h_lock;
-    int prefix_bits;
-    int suffix_bits;
-    ///number of subtable
-	int size;
-    uint64_t suffix_mode;
-    uint64_t non_unique_k_mer;
-} Total_Count_Table;
-
 typedef struct
 {
     uint32_t offset;
@@ -187,21 +170,6 @@ typedef struct
     Chain_Data chainDP;
 } Candidates_list;
 
-typedef struct
-{
-	Pos_Table** sub_h;
-    Hash_table_spin_lock* sub_h_lock;
-    int prefix_bits;
-    int suffix_bits;
-    ///number of subtable
-	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;
-
 ////suffix_bits = 64 in default
 inline int recover_hash_code(uint64_t sub_ID, uint64_t sub_key, Hash_code* code,
                              uint64_t suffix_mode, int suffix_bits, int k) // FIXME: not working right now
@@ -258,204 +226,29 @@ inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, uint64_t s
 	return ha_get_sub_table_long(get_sub_ID, get_sub_key, suffix_mode, suffix_bits, code, k);
 }
 
-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->suffix_mode, TCB->suffix_bits, code, k))
-    {
-        return 0;
-    }
-
-    khint_t t;  
-    int absent;
-
-    while (__sync_lock_test_and_set(&TCB->sub_h_lock[sub_ID].lock, 1))
-    {
-        while (TCB->sub_h_lock[sub_ID].lock);
-    }
-
-    t = ha_ct_put(TCB->sub_h[sub_ID], sub_key, &absent);
-    if (absent)
-    {
-        kh_val(TCB->sub_h[sub_ID], t) = 1;
-    }
-    else   
-    {
-        //kh_value(TCB->sub_h[sub_ID], t) = kh_value(TCB->sub_h[sub_ID], t) + 1;
-        kh_val(TCB->sub_h[sub_ID], t)++;
-    }
-
-    __sync_lock_release(&TCB->sub_h_lock[sub_ID].lock);
-
-    return 1;
-}
-
-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->suffix_mode, TCB->suffix_bits, code, k))
-    {
-        return 0;
-    }
-
-    khint_t t;  
-
-    t = ha_ct_get(TCB->sub_h[sub_ID], sub_key);
-
-    if (t != kh_end(TCB->sub_h[sub_ID]))
-    {
-        return kh_val(TCB->sub_h[sub_ID], t);
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-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;  
-
-    t = ha_pt_get(PCB->sub_h[sub_ID], sub_key);
-
-    if (t != kh_end(PCB->sub_h[sub_ID]))
-    {
-        *r_sub_ID = sub_ID;
-        return kh_val(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++; // if not the last k-mer, this field is reused to keep the number of inserted positions
-            list[list[0].offset].readID = readID;
-            list[list[0].offset].offset = pos;
-            list[list[0].offset].rev = ha_code2rev(code);
-        }
-        else // now comes to the last k-mer position; then save it to list[0]
-        {
-            list[0].readID = readID;
-            list[0].offset = pos;
-			list[0].rev = ha_code2rev(code);
-            flag = 1;
-        }
-
-        __sync_lock_release(&PCB->sub_h_lock[sub_ID].lock);
-
-        //if all pos has been saved, it is safe to sort
-        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 destory_Total_Pos_Table(Total_Pos_Table* TCB);
-void write_Total_Pos_Table(Total_Pos_Table* TCB, char* read_file_name);
-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 init_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list);
 void destory_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list);
 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 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_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error,
 int extra_begin, int extra_end, int error_threshold);
 
-
-
 void overlap_region_sort_y_id(overlap_region *a, long long n);
 
-
 void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, 
 uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end);
 
-
-
 void init_fake_cigar(Fake_Cigar* x);
 void destory_fake_cigar(Fake_Cigar* x);
 void clear_fake_cigar(Fake_Cigar* x);
@@ -463,6 +256,7 @@ void add_fake_cigar(Fake_Cigar* x, uint32_t gap_site, int32_t gap_shift);
 void resize_fake_cigar(Fake_Cigar* x, uint64_t size);
 int get_fake_gap_pos(Fake_Cigar* x, int index);
 int get_fake_gap_shift(Fake_Cigar* x, int index);
+
 inline long long y_start_offset(long long x_start, Fake_Cigar* o)
 {
     if(x_start == get_fake_gap_pos(o, o->length - 1))
@@ -470,7 +264,6 @@ inline long long y_start_offset(long long x_start, Fake_Cigar* o)
         return get_fake_gap_shift(o, o->length - 1);
     }
     
-    
     long long i;
     for (i = 0; i < (long long)o->length; i++)
     {
@@ -508,5 +301,4 @@ void clear_window_list_alloc(window_list_alloc* x);
 void destory_window_list_alloc(window_list_alloc* x);
 void resize_window_list_alloc(window_list_alloc* x, long long size);
 
-
 #endif

htab.h

@@ -10,10 +10,6 @@ typedef struct {
 	uint64_t rid:28, pos:27, rev:1, span:8;
 } ha_mz1_t;
 
-typedef struct {
-	uint64_t rid:28, pos:27, rev:1, span:8;
-} ha_seed_t;
-
 typedef struct { uint32_t n, m; ha_mz1_t *a; } ha_mz1_v;
 
 extern const unsigned char seq_nt4_table[256];