diff --git a/Assembly.cpp b/Assembly.cpp index 4afb606..d4f0921 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -11,426 +11,13 @@ #include "Output.h" #include "htab.h" -Total_Count_Table TCB; -Total_Pos_Table PCB; +void *ha_flt_tab; +void *ha_idx; + All_reads R_INF; pthread_mutex_t statistics; -void* Perform_Counting(void* arg) -{ - 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) +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; @@ -452,7 +39,7 @@ char* new_read, int* new_length) } else if (operation == 1) { - + for (j = 0; j < operation_length; j++) { new_read[new_i] = Get_MisMatch_Base(cigar->lost_base[diff_char_i]); @@ -499,7 +86,7 @@ void get_uncorrected_read_from_cigar(Cigar_record* cigar, char* new_read, int ne } else if (operation == 1) { - + for (j = 0; j < operation_length; j++) { pre_read[pre_i] = Get_Match_Base(cigar->lost_base[diff_char_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, -char* new_read, int new_length, int correct_base) +int debug_cigar(Cigar_record* cigar, char* pre_read, int pre_length, 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]); @@ -589,24 +172,21 @@ char* new_read, int new_length, int correct_base) { new_i = new_i + operation_length; } - + } - - if (pre_i != pre_length) { fprintf(stderr, "pre_i: %d, pre_length: %d\n", pre_i, pre_length); } - - + + if(new_i != new_length) { fprintf(stderr, "new_i: %d, new_length: %d\n", new_i, new_length); } return 1; - char* tmp_seq = (char*)malloc(new_length + pre_length); int tmp_length; @@ -624,12 +204,9 @@ char* new_read, int new_length, int correct_base) } - - - get_uncorrected_read_from_cigar(cigar, new_read, new_length, tmp_seq, &tmp_length); - + if(tmp_length != pre_length) { fprintf(stderr, "tmp_length: %d, pre_length: %d\n", tmp_length, pre_length); @@ -639,7 +216,7 @@ char* new_read, int new_length, int correct_base) { fprintf(stderr, "error pre string\n"); } - + free(tmp_seq); @@ -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, -All_reads* R_INF, int if_reverse) +void push_overlaps(ma_hit_t_alloc* paf, overlap_region_alloc* overlap_list, int flag, All_reads* R_INF, int if_reverse) { long long i = 0, xLen, yLen; ma_hit_t tmp; @@ -689,7 +264,7 @@ All_reads* R_INF, int if_reverse) tmp.qns = overlap_list->list[i].x_id; tmp.qns = tmp.qns << 32; tmp.tn = overlap_list->list[i].y_id; - + if(if_reverse != 0) { tmp.qns = tmp.qns | (uint64_t)(xLen - overlap_list->list[i].x_pos_s - 1); @@ -704,8 +279,8 @@ All_reads* R_INF, int if_reverse) tmp.ts = overlap_list->list[i].y_pos_s; tmp.te = overlap_list->list[i].y_pos_e; } - - + + ///for overlap_list, the x_strand of all overlaps are 0, so the tmp.rev is the same as the y_strand tmp.rev = overlap_list->list[i].y_pos_strand; @@ -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; @@ -743,12 +315,11 @@ long long xBeg, long long xEnd, long long yBeg, long long yEnd) return 1; } } - + return 0; } -long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_list, -overlap_region_alloc* overlap_list, int flag) +long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_list, overlap_region_alloc* overlap_list, int flag) { long long i = 0; long long available_overlaps = 0; @@ -790,7 +361,7 @@ overlap_region_alloc* overlap_list, int flag) ///the end pos is open tmp.te++; /**********************target***************************/ - + tmp.bl = R_INF.read_length[overlap_list->list[i].y_id]; tmp.ml = overlap_list->list[i].strong; tmp.no_l_indel = overlap_list->list[i].without_large_indel; @@ -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) { @@ -950,7 +521,7 @@ void* Overlap_calculate_heap_merge(void* arg) destory_Correct_dumy(&correct); destoryHaplotypeEvdience(&hap); destory_Round2_alignment(&second_round); - + pthread_mutex_lock(&statistics); @@ -988,7 +559,7 @@ void* Output_related_reads(void* arg) init_Graph(&DAGCon); init_Graph(&POA_Graph); - + init_Candidates_list(&l); //init_Candidates_list(&debug_l); @@ -1021,7 +592,7 @@ void* Output_related_reads(void* arg) long long required_read_name_length = strlen(asm_opt.required_read_name); for (i = thr_ID; i < (long long)R_INF.total_reads; i = i + asm_opt.thread_num) { - + if(required_read_name_length == (long long)Get_NAME_LENGTH((R_INF),i) && memcmp(asm_opt.required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) @@ -1043,7 +614,7 @@ void* Output_related_reads(void* arg) recover_UC_Read(&g_read, &R_INF, overlap_list.list[k].y_id); fprintf(stderr, "%.*s\n", (int)g_read.length, g_read.seq); } - + } } @@ -1061,7 +632,7 @@ void* Output_related_reads(void* arg) destory_Correct_dumy(&correct); destoryHaplotypeEvdience(&hap); destory_Round2_alignment(&second_round); - + free(arg); @@ -1103,7 +674,7 @@ void* Save_corrected_reads(void* arg) char* new_read; int new_read_length; - + for (i = thr_ID; i < (long long)R_INF.total_reads; i = i + asm_opt.thread_num) { recover_UC_Read(&g_read, &R_INF, i); @@ -1121,7 +692,7 @@ void* Save_corrected_reads(void* arg) cigar.lost_base = R_INF.cigars[i].lost_base; get_corrected_read_from_cigar(&cigar, g_read.seq, g_read.length, first_round_read, &first_round_read_length); - + /********************************1 round******************************/ /********************************2 round******************************/ @@ -1136,16 +707,16 @@ void* Save_corrected_reads(void* arg) cigar.lost_base = R_INF.second_round_cigar[i].lost_base; get_corrected_read_from_cigar(&cigar, first_round_read, first_round_read_length, second_round_read, &second_round_read_length); - + /********************************2 round******************************/ - + new_read = second_round_read; new_read_length = second_round_read_length; - + if (asm_opt.roundID != asm_opt.number_of_round - 1) { ///need modification @@ -1157,10 +728,10 @@ void* Save_corrected_reads(void* arg) reverse_complement(new_read, new_read_length); } - + N_occ = get_N_occ(new_read, new_read_length); - + if((long long)R_INF.read_size[i] < new_read_length) { R_INF.read_size[i] = new_read_length; @@ -1168,7 +739,7 @@ void* Save_corrected_reads(void* arg) } R_INF.read_length[i] = new_read_length; - + ha_compress_base(Get_READ(R_INF, i), new_read, new_read_length, @@ -1210,9 +781,9 @@ void Output_corrected_reads() void Overlap_calculate_multipe_thr() { double start_time = Get_T(); - + fprintf(stderr, "Begin calculating overlaps... \n"); - + pthread_t *_r_threads; _r_threads = (pthread_t *)malloc(sizeof(pthread_t)*asm_opt.thread_num); @@ -1232,7 +803,7 @@ void Overlap_calculate_multipe_thr() pthread_create(_r_threads + i, NULL, Output_related_reads, (void*)arg); } } - + for (i = 0; i < asm_opt.thread_num; i++) pthread_join(_r_threads[i], NULL); @@ -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"); @@ -1269,14 +841,14 @@ void Overlap_calculate_multipe_thr() fprintf(stderr, "%-30s%18.2f\n\n", "Corrected read saving time:", Get_T() - start_time); - + ///only the last round can output read to disk if (asm_opt.roundID == asm_opt.number_of_round - 1) { start_time = Get_T(); - + Output_corrected_reads(); - + fprintf(stderr, "%-30s%18.2f\n\n", "Output time:", Get_T() - start_time); } } @@ -1357,7 +929,7 @@ UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact) { overlap_list->list[j].is_match = 3; } - + j++; inner_j++; } @@ -1400,7 +972,7 @@ void statistic(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long readNum) long long forward, reverse, strong, weak, exact, no_l_indel; no_l_indel = forward = reverse = exact = strong = weak = 0; long long i, j; - + for (i = 0; i < readNum; i++) { forward += paf[i].length; @@ -1495,7 +1067,7 @@ long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n) // if(bandLen == 0) bandLen = MIN(xRegionLen, yRegionLen); } - + ///do alignment forward kv_resize(uint8_t, x_num, (uint64_t)xRegionLen); kv_resize(uint8_t, y_num, (uint64_t)yRegionLen); @@ -1519,9 +1091,9 @@ UC_Read* g_read, UC_Read* overlap_read, uint8_t* c2n) char* y_string; Cigar_record* cigar; resize_Cigar_record_alloc(cigarline, overlap_list->length); - - + + for (i = 0; i < overlap_list->length; i++) { if(overlap_list->list[i].is_match == 1 || @@ -1554,10 +1126,10 @@ UC_Read* g_read, UC_Read* overlap_read, uint8_t* c2n) Get_READ_LENGTH(R_INF, overlap_list->list[i].x_id), Get_READ_LENGTH(R_INF, overlap_list->list[i].y_id), cigar, c2n); } - + } } - + } void* Final_overlap_calculate_heap_merge(void* arg) @@ -1617,7 +1189,7 @@ void* Final_overlap_calculate_heap_merge(void* arg) push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &overlap_list, 2); - + } finish_output_buffer(); @@ -1628,7 +1200,7 @@ void* Final_overlap_calculate_heap_merge(void* arg) destory_UC_Read(&g_read); destory_UC_Read(&overlap_read); destory_Cigar_record_alloc(&cigarline); - + pthread_mutex_lock(&statistics); asm_opt.complete_threads++; @@ -1685,7 +1257,7 @@ void Output_PAF() fprintf(output_file, "%d\t", sources[i].buffer[j].ml); fprintf(output_file, "%d\t", sources[i].buffer[j].bl); fprintf(output_file, "255\n"); - + } } @@ -1712,7 +1284,7 @@ int check_cluster(uint64_t* list, long long listLen, ma_hit_t_alloc* paf, float { A_edges++; } - + if(get_specific_overlap(&(paf[tn]), tn, qn) != -1) { A_edges++; @@ -1720,7 +1292,7 @@ int check_cluster(uint64_t* list, long long listLen, ma_hit_t_alloc* paf, float T_edges = T_edges + 2; } - + } if(A_edges >= (T_edges*threshold)) @@ -1747,7 +1319,7 @@ long long readNum, long long rescue_threshold, float cluster_threshold) kv_init(edge_vector_index); uint64_t flag; int index; - + for (i = 0; i < readNum; i++) { edge_vector.n = 0; @@ -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(); - Build_hash_table_multiple_thr(); + ha_idx = ha_gen_mzidx(&asm_opt, ha_flt_tab, asm_opt.roundID == 0? 0 : 1, &R_INF); pthread_t *_r_threads; @@ -1835,14 +1405,15 @@ void generate_overlaps(int last_round) *arg = i; pthread_create(_r_threads + i, NULL, Final_overlap_calculate_heap_merge, (void*)arg); } - + for (i = 0; i < asm_opt.thread_num; i++) pthread_join(_r_threads[i], NULL); 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); @@ -1850,7 +1421,7 @@ void generate_overlaps(int last_round) Output_PAF(); trio_partition(); - + 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, @@ -1860,20 +1431,18 @@ 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 { ///fprintf(stderr, "Cannot find overlap file. Please run the whole hifiasm.\n"); } - + clear_opt(&asm_opt, last_round); if(last_round == 0) @@ -1881,12 +1450,12 @@ void Correct_Reads(int last_round) generate_overlaps(last_round); return; } - + fprintf(stderr, "Error correction: Start the %d-th round ...\n", asm_opt.roundID); - - Counting_multiple_thr(); - Build_hash_table_multiple_thr(); + + ha_idx = ha_gen_mzidx(&asm_opt, ha_flt_tab, asm_opt.roundID == 0? 0 : 1, &R_INF); Overlap_calculate_multipe_thr(); + ha_idx_destroy(ha_idx); fprintf(stderr, "Error correction: The %d-th round has been completed.\n", asm_opt.roundID); diff --git a/CommandLines.h b/CommandLines.h index 2dfbbd5..5812255 100644 --- a/CommandLines.h +++ b/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 diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 9aa7f66..a8e3923 100644 --- a/Hash_Table.cpp +++ b/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<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<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; diff --git a/Hash_Table.h b/Hash_Table.h index 7090fdb..1af5fe6 100644 --- a/Hash_Table.h +++ b/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 diff --git a/htab.h b/htab.h index f710d18..4f68d96 100644 --- a/htab.h +++ b/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];