diff --git a/Assembly.cpp b/Assembly.cpp index 0259cf3..cf8c522 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1,437 +1,22 @@ -#include "Assembly.h" #include #include +#include #include +#include "Assembly.h" #include "Process_Read.h" #include "CommandLines.h" -#include "kmer.h" #include "Hash_Table.h" #include "POA.h" #include "Correct.h" #include "Output.h" -#include "Trio.h" +#include "htab.h" + +void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain); -Total_Count_Table TCB; -Total_Pos_Table PCB; 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 - 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; @@ -453,7 +38,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]); @@ -500,7 +85,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]); @@ -539,8 +124,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; @@ -557,13 +141,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]); @@ -590,24 +171,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; @@ -625,12 +203,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); @@ -640,7 +215,7 @@ char* new_read, int new_length, int correct_base) { fprintf(stderr, "error pre string\n"); } - + free(tmp_seq); @@ -653,7 +228,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; @@ -674,8 +248,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; @@ -690,7 +263,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); @@ -705,8 +278,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; @@ -718,12 +291,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; @@ -744,12 +314,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; @@ -791,7 +360,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; @@ -802,316 +371,185 @@ overlap_region_alloc* overlap_list, int flag) } } - return available_overlaps; - } -void get_new_candidates(long long readID, UC_Read* g_read, overlap_region_alloc* overlap_list, k_mer_pos_list_alloc* array_list, -HeapSq* heap, Candidates_list* l, double band_width_threshold, int keep_whole_chain) +typedef struct { + int is_final; + // chaining and overlapping related buffers + UC_Read self_read, ovlp_read; + Candidates_list clist; + overlap_region_alloc olist; + ha_abuf_t *ab; + // error correction related buffers + Cigar_record cigar1; + Graph POA_Graph; + Graph DAGCon; + Correct_dumy correct; + haplotype_evdience_alloc hap; + Round2_alignment round2; +} ha_ovec_buf_t; + +ha_ovec_buf_t *ha_ovec_init(int is_final) { - HPC_seq HPC_read; - Hash_code k_code; - long long avalible_k; - uint64_t code; - uint64_t end_pos; - k_mer_pos* list; - uint64_t list_length; - uint64_t sub_ID; - - clear_Heap(heap); - clear_Candidates_list(l); - - clear_k_mer_pos_list_alloc(array_list); - clear_overlap_region_alloc(overlap_list); - - recover_UC_Read(g_read, &R_INF, readID); - - ///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) - { - list_length = locate_Total_Pos_Table(&PCB, &k_code, &list, asm_opt.k_mer_length, &sub_ID); - - if (list_length != 0) - { - append_k_mer_pos_list_alloc(array_list, list, list_length, end_pos, 0); - } - } - } - else - { - avalible_k = 0; - init_Hash_code(&k_code); - } - } - - - - ///reverse complement strand - reverse_complement(g_read->seq, g_read->length); - init_HPC_seq(&HPC_read, g_read->seq, g_read->length); - init_Hash_code(&k_code); - avalible_k = 0; - - while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6) - { - if(code < 4) - { - k_mer_append(&k_code,code, asm_opt.k_mer_length); - 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); - if (list_length != 0) - { - append_k_mer_pos_list_alloc(array_list, list, list_length, end_pos, 1); - } - } - } - else - { - avalible_k = 0; - init_Hash_code(&k_code); - } - } - - merge_k_mer_pos_list_alloc_heap_sort(array_list, l, heap); - - calculate_overlap_region_by_chaining(l, overlap_list, readID, g_read->length, &R_INF, - band_width_threshold, keep_whole_chain); + ha_ovec_buf_t *b; + CALLOC(b, 1); + b->is_final = !!is_final; + init_UC_Read(&b->self_read); + init_UC_Read(&b->ovlp_read); + init_Candidates_list(&b->clist); + init_overlap_region_alloc(&b->olist); + b->ab = ha_abuf_init(); + if (!b->is_final) { + init_Cigar_record(&b->cigar1); + init_Graph(&b->POA_Graph); + init_Graph(&b->DAGCon); + init_Correct_dumy(&b->correct); + InitHaplotypeEvdience(&b->hap); + init_Round2_alignment(&b->round2); + } + return b; } +void ha_ovec_destroy(ha_ovec_buf_t *b) +{ + destory_UC_Read(&b->self_read); + destory_UC_Read(&b->ovlp_read); + destory_Candidates_list(&b->clist); + destory_overlap_region_alloc(&b->olist); + ha_abuf_destroy(b->ab); + if (!b->is_final) { + destory_Cigar_record(&b->cigar1); + destory_Graph(&b->POA_Graph); + destory_Graph(&b->DAGCon); + destory_Correct_dumy(&b->correct); + destoryHaplotypeEvdience(&b->hap); + destory_Round2_alignment(&b->round2); + } + free(b); +} + +static int64_t ha_Graph_mem(const Graph *g) +{ + int64_t i, mem = 0; + mem = sizeof(Graph) + g->node_q.size * 8 + g->g_nodes.size * sizeof(Node); + for (i = 0; i < (int64_t)g->g_nodes.size; ++i) { + Node *n = &g->g_nodes.list[i]; + mem += n->mismatch_edges.size * sizeof(Edge); + mem += n->deletion_edges.size * sizeof(Edge); + mem += n->insertion_edges.size * sizeof(Edge); + } + mem += g->g_nodes.sort.size * 9; + return mem; +} + +int64_t ha_ovec_mem(const ha_ovec_buf_t *b) +{ + int64_t i, mem = 0, mem_clist, mem_olist; + mem_clist = b->clist.size * sizeof(k_mer_hit) + b->clist.chainDP.size * 7 * 4; + mem_olist = b->olist.size * sizeof(overlap_region); + for (i = 0; i < (int64_t)b->olist.size; ++i) { + const overlap_region *r = &b->olist.list[i]; + mem_olist += r->w_list_size * sizeof(window_list); + mem_olist += r->f_cigar.size * 8; + mem_olist += r->boundary_cigars.size * sizeof(window_list); + } + mem = ha_abuf_mem(b->ab) + mem_clist + mem_olist; + if (!b->is_final) { + mem += sizeof(Cigar_record) + b->cigar1.lost_base_size + b->cigar1.size * 4; + mem += sizeof(Correct_dumy) + b->correct.size * 8; + mem += sizeof(Round2_alignment) + b->round2.cigar.size * 4 + b->round2.tmp_cigar.size * 4; + mem += sizeof(haplotype_evdience_alloc) + b->hap.size * sizeof(haplotype_evdience) + b->hap.snp_matrix_size + b->hap.snp_stat_size * sizeof(SnpStats); + mem += ha_Graph_mem(&b->POA_Graph); + mem += ha_Graph_mem(&b->DAGCon); + } + return mem; +} void* Overlap_calculate_heap_merge(void* arg) { - long long num_read_base = 0; - long long num_correct_base = 0; - long long num_recorrect_base = 0; - int fully_cov, abnormal; + long long num_read_base = 0; + long long num_correct_base = 0; + long long num_recorrect_base = 0; + long long mem_buf; + int fully_cov, abnormal; - int thr_ID = *((int*)arg); - long long i = 0; + int thr_ID = *((int*)arg); + long long i = 0; + ha_ovec_buf_t *b; - UC_Read g_read; - init_UC_Read(&g_read); + b = ha_ovec_init(0); + for (i = thr_ID; i < (long long)R_INF.total_reads; i = i + asm_opt.thread_num) { + //get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.02, 1); + ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1); - UC_Read overlap_read; - init_UC_Read(&overlap_read); + clear_Cigar_record(&b->cigar1); + clear_Round2_alignment(&b->round2); - Candidates_list l; - Graph POA_Graph; - Graph DAGCon; - init_Graph(&DAGCon); - init_Graph(&POA_Graph); - init_Candidates_list(&l); - k_mer_pos_list_alloc array_list; - init_k_mer_pos_list_alloc(&array_list); + correct_overlap(&b->olist, &R_INF, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon, + &b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal); - overlap_region_alloc overlap_list; - init_overlap_region_alloc(&overlap_list); + num_read_base += b->self_read.length; + num_correct_base += b->correct.corrected_base; + num_recorrect_base += b->round2.dumy.corrected_base; - HeapSq heap; + push_cigar(R_INF.cigars, i, &b->cigar1); + push_cigar(R_INF.second_round_cigar, i, &b->round2.cigar); - Init_Heap(&heap); + R_INF.paf[i].is_fully_corrected = 0; + if (fully_cov) { + if (get_cigar_errors(&b->cigar1) == 0 && get_cigar_errors(&b->round2.cigar) == 0) + R_INF.paf[i].is_fully_corrected = 1; + } + R_INF.paf[i].is_abnormal = abnormal; - Correct_dumy correct; - init_Correct_dumy(&correct); + push_overlaps(&(R_INF.paf[i]), &b->olist, 1, &R_INF, asm_opt.roundID%2); + push_overlaps(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, asm_opt.roundID%2); + } + finish_output_buffer(); + mem_buf = ha_ovec_mem(b); + ha_ovec_destroy(b); - - Output_buffer_sub_block current_sub_buffer; - - init_buffer_sub_block(¤t_sub_buffer); - - Cigar_record current_cigar; - init_Cigar_record(¤t_cigar); - - haplotype_evdience_alloc hap; - InitHaplotypeEvdience(&hap); - - Round2_alignment second_round; - init_Round2_alignment(&second_round); - - for (i = thr_ID; i < (long long)R_INF.total_reads; i = i + asm_opt.thread_num) - { - ///get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); - get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02, 1); - - clear_Cigar_record(¤t_cigar); - clear_Round2_alignment(&second_round); - - correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon, - ¤t_cigar, &hap, &second_round, 0, 1, &fully_cov, &abnormal); - - num_read_base += g_read.length; - num_correct_base += correct.corrected_base; - num_recorrect_base += second_round.dumy.corrected_base; - - push_cigar(R_INF.cigars, i, ¤t_cigar); - push_cigar(R_INF.second_round_cigar, i, &(second_round.cigar)); - - - R_INF.paf[i].is_fully_corrected = 0; - if(fully_cov) - { - if(get_cigar_errors(¤t_cigar) == 0 && - get_cigar_errors(&second_round.cigar) == 0) - { - R_INF.paf[i].is_fully_corrected = 1; - } - } - R_INF.paf[i].is_abnormal = abnormal; - - push_overlaps(&(R_INF.paf[i]), &overlap_list, 1, &R_INF, asm_opt.roundID%2); - push_overlaps(&(R_INF.reverse_paf[i]), &overlap_list, 2, &R_INF, asm_opt.roundID%2); - } - - - finish_output_buffer(); - destory_buffer_sub_block(¤t_sub_buffer); - destory_Candidates_list(&l); - destory_overlap_region_alloc(&overlap_list); - destory_Heap(&heap); - destory_k_mer_pos_list_alloc(&array_list); - destory_Graph(&POA_Graph); - destory_Graph(&DAGCon); - destory_UC_Read(&g_read); - destory_UC_Read(&overlap_read); - destory_Cigar_record(¤t_cigar); - destory_Correct_dumy(&correct); - destoryHaplotypeEvdience(&hap); - destory_Round2_alignment(&second_round); - - - - pthread_mutex_lock(&statistics); - asm_opt.num_bases += num_read_base; - asm_opt.num_corrected_bases += num_correct_base; - asm_opt.num_recorrected_bases += num_recorrect_base; - - asm_opt.complete_threads++; - if(asm_opt.complete_threads == asm_opt.thread_num) - { - fprintf(stderr, "total bases #: %lld\n", asm_opt.num_bases); - fprintf(stderr, "total corrected bases: %lld\n", asm_opt.num_corrected_bases); - fprintf(stderr, "total recorrected bases: %lld\n", asm_opt.num_recorrected_bases); - } + pthread_mutex_lock(&statistics); + asm_opt.num_bases += num_read_base; + asm_opt.num_corrected_bases += num_correct_base; + asm_opt.num_recorrected_bases += num_recorrect_base; + asm_opt.mem_buf += mem_buf; pthread_mutex_unlock(&statistics); - free(arg); - - return NULL; + return NULL; } void* Output_related_reads(void* arg) { - int thr_ID = *((int*)arg); - long long i = 0; + int thr_ID = *((int*)arg); + long long i = 0; + ha_ovec_buf_t *b; - UC_Read g_read; - init_UC_Read(&g_read); + long long required_read_name_length = strlen(asm_opt.required_read_name); + b = ha_ovec_init(0); + 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) + { + //get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.02, 1); + ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1); - UC_Read overlap_read; - init_UC_Read(&overlap_read); + fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i)); + recover_UC_Read(&b->self_read, &R_INF, i); + fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq); - Candidates_list l; - Graph POA_Graph; - Graph DAGCon; - init_Graph(&DAGCon); - init_Graph(&POA_Graph); - - - init_Candidates_list(&l); - //init_Candidates_list(&debug_l); - - k_mer_pos_list_alloc array_list; - init_k_mer_pos_list_alloc(&array_list); - - overlap_region_alloc overlap_list; - init_overlap_region_alloc(&overlap_list); - - HeapSq heap; - - Init_Heap(&heap); - - Correct_dumy correct; - init_Correct_dumy(&correct); - - - Output_buffer_sub_block current_sub_buffer; - - init_buffer_sub_block(¤t_sub_buffer); - - Cigar_record current_cigar; - init_Cigar_record(¤t_cigar); - - haplotype_evdience_alloc hap; - InitHaplotypeEvdience(&hap); - - - Round2_alignment second_round; - init_Round2_alignment(&second_round); - - - 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) - { - ////get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5); - get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.02, 1); - - fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), - Get_NAME((R_INF), i)); - recover_UC_Read(&g_read, &R_INF, i); - fprintf(stderr, "%.*s\n", (int)g_read.length, g_read.seq); - - - uint64_t k; - for (k = 0; k < overlap_list.length; k++) - { - fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF),overlap_list.list[k].y_id), - Get_NAME((R_INF),overlap_list.list[k].y_id)); - recover_UC_Read(&g_read, &R_INF, overlap_list.list[k].y_id); - fprintf(stderr, "%.*s\n", (int)g_read.length, g_read.seq); - } - - } - } - - finish_output_buffer(); - - destory_buffer_sub_block(¤t_sub_buffer); - destory_Candidates_list(&l); - destory_overlap_region_alloc(&overlap_list); - destory_Heap(&heap); - destory_k_mer_pos_list_alloc(&array_list); - destory_Graph(&POA_Graph); - destory_Graph(&DAGCon); - destory_UC_Read(&g_read); - destory_UC_Read(&overlap_read); - destory_Cigar_record(¤t_cigar); - destory_Correct_dumy(&correct); - destoryHaplotypeEvdience(&hap); - destory_Round2_alignment(&second_round); - - - free(arg); - - return NULL; + uint64_t k; + for (k = 0; k < b->olist.length; k++) { + fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), b->olist.list[k].y_id), Get_NAME((R_INF), b->olist.list[k].y_id)); + recover_UC_Read(&b->self_read, &R_INF, b->olist.list[k].y_id); + fprintf(stderr, "%.*s\n", (int)b->self_read.length, b->self_read.seq); + } + } + } + finish_output_buffer(); + ha_ovec_destroy(b); + return NULL; } inline long long get_N_occ(char* seq, long long length) @@ -1149,7 +587,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); @@ -1167,7 +605,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******************************/ @@ -1182,16 +620,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 @@ -1203,10 +641,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; @@ -1214,9 +652,9 @@ void* Save_corrected_reads(void* arg) } R_INF.read_length[i] = new_read_length; - - compress_base(Get_READ(R_INF, i), + + ha_compress_base(Get_READ(R_INF, i), new_read, new_read_length, &R_INF.N_site[i], N_occ); } @@ -1224,7 +662,6 @@ void* Save_corrected_reads(void* arg) destory_UC_Read(&g_read); free(first_round_read); free(second_round_read); - free(arg); return NULL; } @@ -1252,83 +689,37 @@ void Output_corrected_reads() fclose(output_file); } - -void Overlap_calculate_multipe_thr() +void ha_overlap_and_correct(int round) { - 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); - - int i = 0; + int i, *args; + pthread_t *_r_threads; + MALLOC(_r_threads, asm_opt.thread_num); + args = (int*)alloca(sizeof(int) * asm_opt.thread_num); + ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, round == 0? 0 : 1, &R_INF); // build the index + for (i = 0; i < asm_opt.thread_num; i++) { + args[i] = i; + if (!asm_opt.required_read_name) + pthread_create(_r_threads + i, NULL, Overlap_calculate_heap_merge, (void*)&args[i]); + else + pthread_create(_r_threads + i, NULL, Output_related_reads, (void*)&args[i]); + } for (i = 0; i < asm_opt.thread_num; i++) - { - int *arg = (int*)malloc(sizeof(*arg)); - *arg = i; - if(!asm_opt.required_read_name) - { - pthread_create(_r_threads + i, NULL, Overlap_calculate_heap_merge, (void*)arg); - } - else - { - 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); + ha_pt_destroy(ha_idx); + ha_idx = 0; + if (asm_opt.required_read_name) exit(0); // for debugging only - free(_r_threads); - - if(asm_opt.required_read_name) - { - exit(1); - } - - - destory_Total_Pos_Table(&PCB); - - fprintf(stderr, "All overlaps have been calculated.\n"); - - fprintf(stderr, "%-30s%18.2f\n\n", "Overlap calculation time:", Get_T() - start_time); - - start_time = Get_T(); - - _r_threads = (pthread_t *)malloc(sizeof(pthread_t)*asm_opt.thread_num); - for (i = 0; i < asm_opt.thread_num; i++) - { - int *arg = (int*)malloc(sizeof(*arg)); - *arg = i; - pthread_create(_r_threads + i, NULL, Save_corrected_reads, (void*)arg); + for (i = 0; i < asm_opt.thread_num; i++) { + args[i] = i; + pthread_create(_r_threads + i, NULL, Save_corrected_reads, (void*)&args[i]); } - - - for (i = 0; i < asm_opt.thread_num; i++) + for (i = 0; i < asm_opt.thread_num; i++) pthread_join(_r_threads[i], NULL); - free(_r_threads); - - 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); - } + free(_r_threads); } - - void update_overlaps(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf, UC_Read* g_read, UC_Read* overlap_read, int is_match, int is_exact) { @@ -1403,7 +794,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++; } @@ -1439,40 +830,31 @@ void update_exact_overlaps(overlap_region_alloc* overlap_list, UC_Read* g_read, } } - -void statistic(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long readNum) +void ha_print_ovlp_stat(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, long long readNum) { + long long forward, reverse, strong, weak, exact, no_l_indel; + long long i, j; - 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; - reverse += rev_paf[i].length; - for (j = 0; j < paf[i].length; j++) - { - if(paf[i].buffer[j].el == 1) exact++; - if(paf[i].buffer[j].ml == 1) strong++; - if(paf[i].buffer[j].ml == 0) weak++; - if(paf[i].buffer[j].no_l_indel == 1) no_l_indel++; - } - } - - - fprintf(stderr, "****************statistic for overlaps****************\n"); - fprintf(stderr, "overlaps #: %lld\n", forward); - fprintf(stderr, "strong overlaps #: %lld\n", strong); - fprintf(stderr, "weak overlaps #: %lld\n", weak); - fprintf(stderr, "exact overlaps #: %lld\n", exact); - fprintf(stderr, "inexact overlaps #: %lld\n", forward - exact); - fprintf(stderr, "overlaps without large indels#: %lld\n", no_l_indel); - fprintf(stderr, "reverse overlaps #: %lld\n", reverse); - fprintf(stderr, "****************statistic for overlaps****************\n"); + no_l_indel = forward = reverse = exact = strong = weak = 0; + for (i = 0; i < readNum; i++) { + forward += paf[i].length; + reverse += rev_paf[i].length; + for (j = 0; j < paf[i].length; j++) { + if (paf[i].buffer[j].el == 1) exact++; + if (paf[i].buffer[j].ml == 1) strong++; + if (paf[i].buffer[j].ml == 0) weak++; + if (paf[i].buffer[j].no_l_indel == 1) no_l_indel++; + } + } + fprintf(stderr, "[M::%s] # overlaps: %lld\n", __func__, forward); + fprintf(stderr, "[M::%s] # strong overlaps: %lld\n", __func__, strong); + fprintf(stderr, "[M::%s] # weak overlaps: %lld\n", __func__, weak); + fprintf(stderr, "[M::%s] # exact overlaps: %lld\n", __func__, exact); // this seems not right + fprintf(stderr, "[M::%s] # inexact overlaps: %lld\n", __func__, forward - exact); + fprintf(stderr, "[M::%s] # overlaps without large indels: %lld\n", __func__, no_l_indel); + fprintf(stderr, "[M::%s] # reverse overlaps: %lld\n", __func__, reverse); } - void fill_chain(Fake_Cigar* chain, char* x_string, char* y_string, long long xBeg, long long yBeg, long long x_readLen, long long y_readLen, Cigar_record* cigar, uint8_t* c2n) { @@ -1541,7 +923,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); @@ -1565,9 +947,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 || @@ -1600,51 +982,29 @@ 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) { int thr_ID = *((int*)arg); uint64_t i = 0; - - UC_Read g_read; - init_UC_Read(&g_read); - - UC_Read overlap_read; - init_UC_Read(&overlap_read); - - Candidates_list l; - init_Candidates_list(&l); - - k_mer_pos_list_alloc array_list; - init_k_mer_pos_list_alloc(&array_list); - - overlap_region_alloc overlap_list; - init_overlap_region_alloc(&overlap_list); - - HeapSq heap; - Init_Heap(&heap); - - Cigar_record_alloc cigarline; - init_Cigar_record_alloc(&cigarline); + ha_ovec_buf_t *b; uint8_t c2n[256]; memset(c2n, 4, 256); c2n[(uint8_t)'A'] = c2n[(uint8_t)'a'] = 0; c2n[(uint8_t)'C'] = c2n[(uint8_t)'c'] = 1; c2n[(uint8_t)'G'] = c2n[(uint8_t)'g'] = 2; c2n[(uint8_t)'T'] = c2n[(uint8_t)'t'] = 3; // build the encoding table - - - - + b = ha_ovec_init(1); for (i = thr_ID; i < R_INF.total_reads; i = i + asm_opt.thread_num) { + //get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0); + ha_get_new_candidates(b->ab, i, &b->self_read, &b->olist, &b->clist, 0.001, asm_opt.max_n_chain, 0); - get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, 0.001, 0); /** correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon, &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1, @@ -1652,57 +1012,27 @@ void* Final_overlap_calculate_heap_merge(void* arg) push_final_overlaps(&(R_INF.paf[i]), &overlap_list); **/ - overlap_region_sort_y_id(overlap_list.list, overlap_list.length); + overlap_region_sort_y_id(b->olist.list, b->olist.length); ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length); ma_hit_sort_tn(R_INF.reverse_paf[i].buffer, R_INF.reverse_paf[i].length); - reverse_complement(g_read.seq, g_read.length); - - update_overlaps(&overlap_list, &(R_INF.paf[i]), &g_read, &overlap_read, 1, 1); - update_overlaps(&overlap_list, &(R_INF.reverse_paf[i]), &g_read, &overlap_read, 2, 0); + update_overlaps(&b->olist, &(R_INF.paf[i]), &b->self_read, &b->ovlp_read, 1, 1); + update_overlaps(&b->olist, &(R_INF.reverse_paf[i]), &b->self_read, &b->ovlp_read, 2, 0); ///recover missing exact overlaps - update_exact_overlaps(&overlap_list, &g_read, &overlap_read); + update_exact_overlaps(&b->olist, &b->self_read, &b->ovlp_read); ///Final_phasing(&overlap_list, &cigarline, &g_read, &overlap_read, c2n); - push_final_overlaps(&(R_INF.paf[i]), R_INF.reverse_paf, - &overlap_list, 1); - push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, - &overlap_list, 2); - - + push_final_overlaps(&(R_INF.paf[i]), R_INF.reverse_paf, &b->olist, 1); + push_final_overlaps(&(R_INF.reverse_paf[i]), R_INF.reverse_paf, &b->olist, 2); } - finish_output_buffer(); - - destory_Candidates_list(&l); - destory_overlap_region_alloc(&overlap_list); - destory_Heap(&heap); - destory_k_mer_pos_list_alloc(&array_list); - destory_UC_Read(&g_read); - destory_UC_Read(&overlap_read); - destory_Cigar_record_alloc(&cigarline); - - - pthread_mutex_lock(&statistics); - asm_opt.complete_threads++; - if(asm_opt.complete_threads == asm_opt.thread_num) - { - if(VERBOSE >= 1) - { - statistic(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads); - } - } - pthread_mutex_unlock(&statistics); - free(arg); - + ha_ovec_destroy(b); return NULL; } - void Output_PAF() { - fprintf(stderr, "Writing PAF to disk ...... \n"); char* paf_name = (char*)malloc(strlen(asm_opt.output_file_name)+50); sprintf(paf_name, "%s.ovlp.paf", asm_opt.output_file_name); @@ -1710,8 +1040,6 @@ void Output_PAF() uint64_t i, j; ma_hit_t_alloc* sources = R_INF.paf; - - for (i = 0; i < R_INF.total_reads; i++) { for (j = 0; j < sources[i].length; j++) @@ -1739,7 +1067,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"); - + } } @@ -1749,7 +1077,6 @@ void Output_PAF() fprintf(stderr, "PAF has been written.\n"); } - int check_cluster(uint64_t* list, long long listLen, ma_hit_t_alloc* paf, float threshold) { long long i, k; @@ -1766,7 +1093,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++; @@ -1774,7 +1101,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)) @@ -1801,7 +1128,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; @@ -1866,95 +1193,64 @@ long long readNum, long long rescue_threshold, float cluster_threshold) fprintf(stderr, "[M::%s] took %0.2fs, revise edges #: %lld\n\n", __func__, Get_T()-startTime, revises); } - - -void generate_overlaps(int last_round) +void ha_overlap_final(void) { - double start_time = Get_T(); - asm_opt.roundID = asm_opt.number_of_round - last_round; - fprintf(stderr, "Begin calculting final overlaps ...\n"); + int i, *args; + pthread_t *_r_threads; - Counting_multiple_thr(); - Build_hash_table_multiple_thr(); + MALLOC(_r_threads, asm_opt.thread_num); + args = (int*)alloca(sizeof(int) * asm_opt.thread_num); - 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; - 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); - - fprintf(stderr, "Final overlaps have been calculated.\n"); - fprintf(stderr, "%-30s%18.2f\n\n", "Final overlaps calculation time:", Get_T() - start_time); - - 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, - asm_opt.large_pop_bubble_size, 0, 1); + ha_idx = ha_pt_gen(&asm_opt, ha_flt_tab, 1, &R_INF); + for (i = 0; i < asm_opt.thread_num; i++) { + args[i] = i; + pthread_create(_r_threads + i, NULL, Final_overlap_calculate_heap_merge, (void*)&args[i]); + } + for (i = 0; i < asm_opt.thread_num; i++) + pthread_join(_r_threads[i], NULL); + free(_r_threads); + ha_pt_destroy(ha_idx); + ha_idx = 0; + ///rescue_edges(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads, 4, 0.985); } - -void Correct_Reads(int last_round) +int ha_assemble(void) { - - 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); - } - else - { - ///fprintf(stderr, "Cannot find overlap file. Please run the whole hifiasm.\n"); - } - - clear_opt(&asm_opt, last_round); - - if(last_round == 0) - { - 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(); - Overlap_calculate_multipe_thr(); - - fprintf(stderr, "Error correction: The %d-th round has been completed.\n", asm_opt.roundID); - - Correct_Reads(last_round - 1); + int r, ovlp_loaded = 0; + if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) { + ovlp_loaded = 1; + fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage()); + } + if (!ovlp_loaded) { + // construct hash table for high occurrence k-mers + if (!asm_opt.no_kmer_flt) + ha_flt_tab = ha_ft_gen(&asm_opt, &R_INF); + // error correction + assert(asm_opt.number_of_round > 0); + for (r = 0; r < asm_opt.number_of_round; ++r) { + clear_opt(&asm_opt, r); // this update asm_opt.roundID and a few other fields + ha_overlap_and_correct(r); + fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> corrected reads for round %d\n", __func__, yak_realtime(), + yak_cpu_usage(), yak_peakrss_in_gb(), r + 1); + fprintf(stderr, "[M::%s] # bases: %lld; # corrected bases: %lld; # recorrected bases: %lld\n", __func__, + asm_opt.num_bases, asm_opt.num_corrected_bases, asm_opt.num_recorrected_bases); + fprintf(stderr, "[M::%s] size of buffer: %.3fGB\n", __func__, asm_opt.mem_buf / 1073741824.0); + } + Output_corrected_reads(); + fprintf(stderr, "[M::%s::%.3f*%.2f] ==> written corrected reads to disk\n", __func__, yak_realtime(), yak_cpu_usage()); + // overlap between corrected reads + clear_opt(&asm_opt, asm_opt.number_of_round); + ha_overlap_final(); + fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> found overlaps for the final round\n", __func__, yak_realtime(), + yak_cpu_usage(), yak_peakrss_in_gb()); + ha_print_ovlp_stat(R_INF.paf, R_INF.reverse_paf, R_INF.total_reads); + ha_ft_destroy(ha_flt_tab); + 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, asm_opt.large_pop_bubble_size, 0, !ovlp_loaded); + destory_All_reads(&R_INF); + return 0; } - - - - - - - - - - - - diff --git a/Assembly.h b/Assembly.h index 557aa8a..339d5e1 100644 --- a/Assembly.h +++ b/Assembly.h @@ -8,8 +8,6 @@ #define Get_Cigar_Type(RECORD) (RECORD&3) #define Get_Cigar_Length(RECORD) (RECORD>>2) -void Counting_multiple_thr(); -void Build_hash_table_multiple_thr(); -void Overlap_calculate_multipe_thr(); -void Correct_Reads(int last_round); +int ha_assemble(void); + #endif diff --git a/CommandLines.cpp b/CommandLines.cpp index 8734c47..8a1eb93 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -5,7 +5,6 @@ #include "ketopt.h" #include -#define VERSION "0.3.0" #define DEFAULT_OUTPUT "hifiasm.asm" hifiasm_opt_t asm_opt; @@ -23,12 +22,14 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, "Options:\n"); fprintf(stderr, " Assembly:\n"); fprintf(stderr, " -o FILE prefix of output files [%s]\n", asm_opt->output_file_name); - ///fprintf(stderr, " -c FILE file including trio information\n"); fprintf(stderr, " -t INT number of threads [%d]\n", asm_opt->thread_num); fprintf(stderr, " -r INT round of correction [%d]\n", asm_opt->number_of_round); fprintf(stderr, " -a INT round of assembly cleaning [%d]\n", asm_opt->clean_round); - fprintf(stderr, " -k INT k-mer length [%d] (must be < 64)\n", asm_opt->k_mer_length); - ///fprintf(stderr, " -w write all overlaps to disk, can accelerate assembly next time [%d]\n", asm_opt->write_index_to_disk); + fprintf(stderr, " -k INT k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length); + fprintf(stderr, " -w INT minimizer window size [%d]\n", asm_opt->mz_win); + fprintf(stderr, " -f INT number of bits for bloom filter [%d]\n", asm_opt->bf_shift); + fprintf(stderr, " -D FLOAT drop k-mers occuring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor); + fprintf(stderr, " -N INT consider up to INT overlaps for each oriented read [%d]\n", asm_opt->max_n_chain); ///fprintf(stderr, " -l load all overlaps from disk, can avoid overlap calculation [%d]\n", asm_opt->load_index_from_disk); ///fprintf(stderr, " -i ignore saved overlaps in *.ovlp*.bin files\n"); fprintf(stderr, " -i ignore saved overlaps in *.ovlp* files\n"); @@ -63,7 +64,13 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->pat_index = NULL; asm_opt->mat_index = NULL; asm_opt->thread_num = 1; - asm_opt->k_mer_length = 40; + asm_opt->k_mer_length = 51; + asm_opt->mz_win = 51; + asm_opt->bf_shift = 37; + asm_opt->high_factor = 5.0f; + asm_opt->no_HPC = 0; + asm_opt->no_kmer_flt = 0; + asm_opt->max_n_chain = 400; asm_opt->k_mer_min_freq = 3; asm_opt->k_mer_max_freq = 66; asm_opt->load_index_from_disk = 1; @@ -71,7 +78,6 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->number_of_round = 2; asm_opt->adapterLen = 0; asm_opt->clean_round = 4; - asm_opt->complete_threads = 0; asm_opt->small_pop_bubble_size = 100000; asm_opt->large_pop_bubble_size = 10000000; asm_opt->min_drop_rate = 0.2; @@ -94,13 +100,13 @@ void destory_opt(hifiasm_opt_t* asm_opt) } } -void clear_opt(hifiasm_opt_t* asm_opt, int last_round) +void clear_opt(hifiasm_opt_t* asm_opt, int round) { - asm_opt->complete_threads = 0; asm_opt->num_bases = 0; asm_opt->num_corrected_bases = 0; asm_opt->num_recorrected_bases = 0; - asm_opt->roundID = asm_opt->number_of_round - last_round; + asm_opt->mem_buf = 0; + asm_opt->roundID = round; } int check_file(char* name, const char* opt) @@ -295,7 +301,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) int c; - while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:lwm:n:r:a:b:z:x:y:p:c:d:M:P:i", 0)) >= 0) { + while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:lw:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:", 0)) >= 0) { if (c == 'h') { Print_H(asm_opt); @@ -303,16 +309,20 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) } else if (c == 'v') { - fprintf(stderr, "[Version] %s\n", VERSION); + puts(HA_VERSION); return 0; - } + } + else if (c == 'f') asm_opt->bf_shift = atoi(opt.arg); else if (c == 't') asm_opt->thread_num = atoi(opt.arg); else if (c == 'o') asm_opt->output_file_name = opt.arg; else if (c == 'r') asm_opt->number_of_round = atoi(opt.arg); else if (c == 'k') asm_opt->k_mer_length = atoi(opt.arg); else if (c == 'i') asm_opt->load_index_from_disk = 0; else if (c == 'l') asm_opt->load_index_from_disk = 1; - else if (c == 'w') asm_opt->write_index_to_disk = 1; + else if (c == 'w') asm_opt->mz_win = atoi(opt.arg); + else if (c == 'D') asm_opt->high_factor = atof(opt.arg); + else if (c == 'F') asm_opt->no_kmer_flt = 1; + else if (c == 'N') asm_opt->max_n_chain = atoi(opt.arg); else if (c == 'a') asm_opt->clean_round = atoi(opt.arg); else if (c == 'z') asm_opt->adapterLen = atoi(opt.arg); else if (c == 'b') asm_opt->required_read_name = opt.arg; diff --git a/CommandLines.h b/CommandLines.h index fa46850..7043dc7 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,6 +3,8 @@ #include +#define HA_VERSION "0.3.0-dirty-r192" + #define VERBOSE 0 #define VERBOSE_GFA 1 @@ -15,6 +17,12 @@ typedef struct { char* mat_index; int thread_num; int k_mer_length; + int mz_win; + int bf_shift; + float high_factor; + int no_HPC; + int no_kmer_flt; + int max_n_chain; int k_mer_min_freq; int k_mer_max_freq; int load_index_from_disk; @@ -22,7 +30,6 @@ typedef struct { int number_of_round; int adapterLen; int clean_round; - int complete_threads; int roundID; int max_hang_Len; int gap_fuzz; @@ -41,6 +48,7 @@ typedef struct { long long num_bases; long long num_corrected_bases; long long num_recorrected_bases; + long long mem_buf; long long coverage; } hifiasm_opt_t; @@ -48,8 +56,8 @@ 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); +void clear_opt(hifiasm_opt_t* asm_opt, int round); +int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt); double Get_T(void); -#endif \ No newline at end of file +#endif diff --git a/Correct.cpp b/Correct.cpp index 4ca6da9..1733e94 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -7213,8 +7213,6 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, haplotype_evdience_alloc* hap, Round2_alignment* second_round, int force_repeat, int is_consensus, int* fully_cov, int* abnormal) { - reverse_complement(g_read->seq, g_read->length); - clear_Correct_dumy(dumy, overlap_list); long long window_start, window_end; @@ -7261,7 +7259,6 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, (*fully_cov) = check_if_fully_covered(overlap_list, R_INF, g_read, dumy, g, abnormal); - } diff --git a/Correct.h b/Correct.h index 791106a..01f3c1c 100644 --- a/Correct.h +++ b/Correct.h @@ -1173,4 +1173,4 @@ long long* max_t_pos, long long* max_q_pos, long long* score, long long* droped) #define GAP_EXT_KSW 2 #define Z_DROP_KSW 400 #define BAND_KSW 50 -#endif \ No newline at end of file +#endif diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 19e4cd4..1c83cc0 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -1,12 +1,11 @@ #include #include #include +#include #include "Hash_Table.h" #include "Process_Read.h" #include "Correct.h" #include "CommandLines.h" -#include "kmer.h" -#include #include "ksort.h" pthread_mutex_t output_mutex; @@ -18,141 +17,6 @@ void overlap_region_sort_y_id(overlap_region *a, long long n) radix_sort_overlap_region_sort(a, a + n); } -void Init_Heap(HeapSq* HBT) -{ - HBT->MaxSize = 1000; - HBT->heap = (ElemType*)malloc(HBT->MaxSize*sizeof(ElemType)); - HBT->index_i = (uint64_t*)malloc(HBT->MaxSize*sizeof(uint64_t)); - HBT->len = 0; -} - -void destory_Heap(HeapSq* HBT) -{ - free(HBT->heap); - free(HBT->index_i); -} - -void clear_Heap(HeapSq* HBT) -{ - HBT->len = 0; -} - - - -inline int cmp_ElemType(ElemType* x, ElemType* y) -{ - - if (x->node.readID < y->node.readID) - { - return 1; - } - else if (x->node.readID > y->node.readID) - { - return 2; - } - else - { if (x->node.strand < y->node.strand) - { - return 1; - } - else if (x->node.strand > y->node.strand) - { - return 2; - } - else - { - - if(x->node.offset < y->node.offset) - { - return 1; - } - else if(x->node.offset > y->node.offset) - { - return 2; - } - else - { - if (x->node.self_offset < y->node.self_offset) - { - return 1; - } - else if (x->node.self_offset > y->node.self_offset) - { - return 2; - } - else ///if both r_pos_x and self_offset are equal, offset must be equal - { - return 0; - } - - } - } - } -} - - -inline void Insert_Heap(HeapSq* HBT, ElemType* x) -{ - long long i, j; - if (HBT->len == HBT->MaxSize) - { - HBT->MaxSize = 2*HBT->MaxSize; - HBT->heap = (ElemType*)realloc(HBT->heap, HBT->MaxSize*sizeof(ElemType)); - HBT->index_i = (uint64_t*)realloc(HBT->index_i,HBT->MaxSize*sizeof(uint64_t)); - } - HBT->heap[HBT->len] = *x; //add element to tail - HBT->len++; - i = HBT->len - 1; - while (i != 0) - { - j = (i - 1) / 2; - ///if (x >= HBT->heap[j]) - ///1: xheap[j])!=1) - break; - HBT->heap[i] = HBT->heap[j]; - i = j; - } - HBT->heap[i] = *x; -} - -inline int DeleteHeap(HeapSq* HBT, ElemType* get) -{ - ElemType temp, x; - int i, j; - if (HBT->len == 0) - { - return 0; - } - temp = HBT->heap[0]; - HBT->len--; - if (HBT->len == 0) - { - *get = temp; - return 2; - } - - x = HBT->heap[HBT->len]; - i = 0; - j = 2 * i + 1; - while (j <= HBT->len - 1) - { - ///if (j < HBT->len - 1 && HBT->heap[j] > HBT->heap[j+1]) - if (j < HBT->len - 1 && cmp_ElemType(&HBT->heap[j], &HBT->heap[j + 1]) == 2) - j++; - ///if (x <= HBT->heap[j]) - if (cmp_ElemType(&x, &HBT->heap[j])!=2) - break; - HBT->heap[i] = HBT->heap[j]; - i = j; - j = 2 * i + 1; - } - HBT->heap[i] = x; - - *get = temp; - return 1; -} - void init_overlap_region_alloc(overlap_region_alloc* list) { list->size = 1000; @@ -166,6 +30,7 @@ void init_overlap_region_alloc(overlap_region_alloc* list) init_window_list_alloc(&(list->list[i].boundary_cigars)); } } + void clear_overlap_region_alloc(overlap_region_alloc* list) { list->length = 0; @@ -194,7 +59,6 @@ void destory_overlap_region_alloc(overlap_region_alloc* list) free(list->list); } - int get_fake_gap_pos(Fake_Cigar* x, int index) { return (x->buffer[index]>>32); @@ -219,10 +83,8 @@ int get_fake_gap_shift(Fake_Cigar* x, int index) return result; } - - int append_inexact_overlap_region_alloc(overlap_region_alloc* list, overlap_region* tmp, -All_reads* R_INF, int add_beg_end) + All_reads* R_INF, int add_beg_end) { if (list->length + 1 > list->size) @@ -264,8 +126,6 @@ All_reads* R_INF, int add_beg_end) long long x_right_length = Get_READ_LENGTH((*R_INF), tmp->x_id) - tmp->x_pos_e - 1; long long y_right_length = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_e - 1; - - if(x_right_length <= y_right_length) { tmp->x_pos_e = Get_READ_LENGTH((*R_INF), tmp->x_id) - 1; @@ -276,7 +136,6 @@ All_reads* R_INF, int add_beg_end) tmp->x_pos_e = tmp->x_pos_e + y_right_length; tmp->y_pos_e = Get_READ_LENGTH((*R_INF), tmp->y_id) - 1; } - if (tmp->x_pos_strand == 1) { @@ -290,9 +149,6 @@ All_reads* R_INF, int add_beg_end) list->list[list->length].y_pos_s = Get_READ_LENGTH((*R_INF), tmp->y_id) - tmp->y_pos_e - 1; list->list[list->length].y_pos_strand = 1; - - - resize_fake_cigar(&(list->list[list->length].f_cigar), (tmp->f_cigar.length + 2)); if(add_beg_end == 1) { @@ -388,7 +244,6 @@ All_reads* R_INF, int add_beg_end) return 1; } - void append_overlap_region_alloc_debug(overlap_region_alloc* list, overlap_region* tmp) { @@ -414,67 +269,6 @@ void append_overlap_region_alloc_debug(overlap_region_alloc* list, overlap_regio list->length++; } - -int cmp_by_x_pos_s(const void * a, const void * b) -{ - if ((*(overlap_region*)a).x_pos_s > (*(overlap_region*)b).x_pos_s) - { - return 1; - } - else if ((*(overlap_region*)a).x_pos_s < (*(overlap_region*)b).x_pos_s) - { - return -1; - } - else - { - - if ((*(overlap_region*)a).x_pos_e > (*(overlap_region*)b).x_pos_e) - { - return 1; - } - else if ((*(overlap_region*)a).x_pos_e < (*(overlap_region*)b).x_pos_e) - { - return -1; - } - else - { - return 0; - } - - } -} - - -int cmp_by_x_pos_e(const void * a, const void * b) -{ - if ((*(overlap_region*)a).x_pos_e > (*(overlap_region*)b).x_pos_e) - { - return 1; - } - else if ((*(overlap_region*)a).x_pos_e < (*(overlap_region*)b).x_pos_e) - { - return -1; - } - else - { - - if ((*(overlap_region*)a).x_pos_s > (*(overlap_region*)b).x_pos_s) - { - return 1; - } - else if ((*(overlap_region*)a).x_pos_s < (*(overlap_region*)b).x_pos_s) - { - return -1; - } - else - { - return 0; - } - - } -} - - void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp) { long long i, j, current_j; @@ -503,21 +297,19 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp) if(indels != dp->indels[i]) { - fprintf(stderr, "indels: %lld, dp->indels[i]: %lld\n", - indels, dp->indels[i]); + fprintf(stderr, "indels: %lld, dp->indels[i]: %ld\n", indels, (long)dp->indels[i]); } if(selfLen != dp->self_length[i]) { - fprintf(stderr, "selfLen: %lld, dp->self_length[i]: %lld\n", - selfLen, dp->self_length[i]); + fprintf(stderr, "selfLen: %lld, dp->self_length[i]: %ld\n", selfLen, (long)dp->self_length[i]); } } } long long get_chainLen(long long x_beg, long long x_end, long long xLen, -long long y_beg, long long y_end, long long yLen) + long long y_beg, long long y_end, long long yLen) { if(x_beg <= y_beg) { @@ -548,13 +340,70 @@ long long y_beg, long long y_end, long long yLen) return x_end - x_beg + 1; } +static int32_t ha_kmer_hit_lis(int32_t n, const k_mer_hit *a, int32_t *b, int32_t *M) +{ + int32_t i, k, L = 0, *P = b; + for (i = 0; i < n; ++i) { + int32_t lo = 1, hi = L, newL; + while (lo <= hi) { + int32_t mid = (lo + hi + 1) >> 1; + if (a[M[mid]].offset < a[i].offset) lo = mid + 1; + else hi = mid - 1; + } + newL = lo, P[i] = M[newL - 1], M[newL] = i; + if (newL > L) L = newL; + } + k = M[L]; + memcpy(M, P, n * sizeof(int32_t)); + for (i = L - 1; i >= 0; --i) b[i] = k, k = M[k]; + return L; +} + +int32_t ha_chain_lis_core(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc, double bw_thres) +{ + int32_t *tmp = (int32_t*)dp->tmp; + int32_t i, m, *b = tmp, *M = tmp + n_a; + int32_t tot_indel = 0, tot_len = 0; + if (n_a < 2) return -1; + for (i = 1; i < n_a; ++i) + if (a[i-1].offset >= a[i].offset) + break; + if (i == n_a) { + for (i = 0; i < n_a; ++i) + b[i] = i; + m = n_a; + } else m = ha_kmer_hit_lis(n_a, a, b, M); + dp->score[0] = 0, dp->pre[0] = -1, dp->indels[0] = 0, dp->self_length[0] = 0; + for (i = 1; i < m; ++i) { + int32_t j0 = b[i-1], j1 = b[i], score, dg; + int32_t dx = (int32_t)a[j1].offset - (int32_t)a[j0].offset; + int32_t dy = (int32_t)a[j1].self_offset - (int32_t)a[j0].self_offset; + int32_t dd = dx > dy? dx - dy : dy - dx; + double gap_rate; + tot_indel += dd; + tot_len += dy; + if (tot_indel > tot_len * bw_thres) + break; + dg = dx < dy? dx : dy; + score = dg < min_sc? dg : min_sc; + gap_rate = (double)tot_indel / tot_len; + score -= (int)(gap_rate * score * bw_thres); + dp->score[i] = dp->score[i-1] + score; + dp->pre[i] = i - 1; + dp->indels[i] = tot_indel; + dp->self_length[i] = tot_len; + } + if (i < m) return -1; + for (i = 0; i < m; ++i) a[i] = a[b[i]]; + return m; +} ///double band_width_threshold = 0.05; void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, -double band_width_threshold, int max_skip, int x_readLen, int y_readLen) + double band_width_threshold, int max_skip, int x_readLen, int y_readLen) { long long i, j; - long long self_pos, pos, max_j, max_i, max_score, score, n_skip; + long long self_pos, pos, max_j, max_i, max_score, score; long long distance_pos, distance_self_pos, distance_gap, distance_min; ///double band_width_threshold = 0.05; double band_width_penalty = 1 / band_width_threshold; @@ -562,19 +411,29 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) long long max_indels, max_self_length; double gap_rate; long long total_indels, total_self_length; + int32_t ret; resize_Chain_Data(dp, a_n); + + ret = ha_chain_lis_core(a, a_n, dp, min_score, band_width_threshold); + if (ret > 0) { + a_n = ret; + goto skip_dp; + } + // fill the score and backtrack arrays + for (i = 0; i < a_n; ++i) dp->tmp[i] = -1; for (i = 0; i < a_n; ++i) { + int n_chn_skip = 0; + int n_max_skip = 0; + pos = a[i].offset; self_pos = a[i].self_offset; max_j = -1; max_score = min_score; - n_skip = 0; max_indels = 0; max_self_length = 0; - ///may have a pre-cut condition for j for (j = i - 1; j >= 0; --j) @@ -608,26 +467,23 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) score += dp->score[j]; - ///find a new max score - if(score > max_score) - { - max_score = score; - max_j = j; - max_indels = total_indels; - max_self_length = total_self_length; - /****************************may have bugs********************************/ - n_skip = 0; - /****************************may have bugs********************************/ - }/****************************may have bugs********************************/ - else - { - n_skip++; - if(n_skip > max_skip) - { - break; - } - } - /****************************may have bugs********************************/ + ///find a new max score + if (score > max_score) { + max_score = score; + max_j = j; + max_indels = total_indels; + max_self_length = total_self_length; + n_max_skip = 0; + if (n_chn_skip > 0) --n_chn_skip; + } else { + if (++n_max_skip > max_skip) + break; + if (dp->tmp[j] == i) { + if (++n_chn_skip > max_skip) + break; + } + } + if (dp->pre[j] >= 0) dp->tmp[dp->pre[j]] = i; } dp->score[i] = max_score; @@ -636,11 +492,9 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) dp->self_length[i] = max_self_length; } - ///debug_chain(a, a_n, dp); - - +skip_dp: max_score = -1; max_i = -1; @@ -652,12 +506,12 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) max_score = dp->score[i]; max_i = i; mini_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen, - a[i].offset, a[i].offset, y_readLen); + a[i].offset, a[i].offset, y_readLen); } else if(dp->score[i] == max_score) { tmp_xLen = get_chainLen(a[i].self_offset, a[i].self_offset, x_readLen, - a[i].offset, a[i].offset, y_readLen); + a[i].offset, a[i].offset, y_readLen); if(tmp_xLen < mini_xLen) { @@ -670,7 +524,6 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) } - clear_fake_cigar(&(result->f_cigar)); ///note a has been sorted by offset, that means has been sorted by query offset i = max_i; @@ -709,7 +562,6 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) } else { - while (i >= 0) { distance_self_pos = result->x_pos_e - a[i].self_offset; @@ -734,10 +586,8 @@ double band_width_threshold, int max_skip, int x_readLen, int y_readLen) } } - - 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) + uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end) { overlap_region tmp_region; long long i = 0; @@ -768,8 +618,6 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh ///here the strand of query is always 0 tmp_region.y_pos_strand = 0; - - sub_region_beg = i; sub_region_end = i; i++; @@ -789,32 +637,22 @@ uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_thresh continue; } - - chain_DP(candidates->list + sub_region_beg, - sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold, - 50, Get_READ_LENGTH((*R_INF), tmp_region.x_id), Get_READ_LENGTH((*R_INF), tmp_region.y_id)); + chain_DP(candidates->list + sub_region_beg, + sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold, + 25, Get_READ_LENGTH((*R_INF), tmp_region.x_id), Get_READ_LENGTH((*R_INF), tmp_region.y_id)); - // chain_DP_back(candidates->list + sub_region_beg, - // sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold); - - ///自己和自己重叠的要排除 ///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1) if (tmp_region.x_id != tmp_region.y_id) { append_inexact_overlap_region_alloc(overlap_list, &tmp_region, R_INF, add_beg_end); - ///append_inexact_overlap_region_alloc_back(overlap_list, &tmp_region, R_INF); } } destory_fake_cigar(&(tmp_region.f_cigar)); - - qsort(overlap_list->list, overlap_list->length, sizeof(overlap_region), cmp_by_x_pos_s); } - - 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) + int extra_begin, int extra_end, int error_threshold) { long long length = region->x_pos_e - region->x_pos_s + 1; @@ -843,73 +681,6 @@ int extra_begin, int extra_end, int error_threshold) region->w_list_length++; } - - -void init_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) -{ - list->size = 1000; - list->length = 0; - //list->list = (k_mer_pos_list*)malloc(sizeof(k_mer_pos_list)*list->size); - list->list = (k_mer_pos_list*)calloc(list->size, sizeof(k_mer_pos_list)); -} - -void clear_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) -{ - list->length = 0; -} - -void destory_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list) -{ - free(list->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) -{ - - if (list->length + 1 > list->size) - { - list->size = list->size * 2; - list->list = (k_mer_pos_list*)realloc(list->list, sizeof(k_mer_pos_list)*list->size); - } - - list->list[list->length].list = n_list; - list->list[list->length].length = n_length; - list->list[list->length].direction = n_direction; - list->list[list->length].end_pos = n_end_pos; - - list->length++; -} - - - - -int cmp_k_mer_pos_list(const void * a, const void * b) -{ - if ((*(k_mer_pos_list*)a).length > (*(k_mer_pos_list*)b).length) - { - return 1; - } - else if ((*(k_mer_pos_list*)a).length < (*(k_mer_pos_list*)b).length) - { - return -1; - } - else - { - return 0; - } - -} - -inline void append_pos_to_Candidates_list(Candidates_list* candidates, ElemType* x) -{ - candidates->list[candidates->length] = x->node; - candidates->length++; -} - - - - void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n_lengh, uint64_t end_pos, uint64_t strand) { uint64_t i; @@ -938,717 +709,12 @@ void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n fprintf(stderr, "ERROR 4\n"); } } - -} - - -void merge_k_mer_pos_list_alloc_heap_sort(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT) -{ - clear_Heap(HBT); - - uint64_t total_length = 0; - uint64_t i; - ElemType x, y; - //add the first element of each list to stack - for (i = 0; i < list->length; i++) - { - x.ID = i; - x.node.offset = list->list[i].list[0].offset; - x.node.readID = list->list[i].list[0].readID; - x.node.self_offset = list->list[i].end_pos; - x.node.strand = list->list[i].direction; - - Insert_Heap(HBT, &x); - - HBT->index_i[i] = 1; - - total_length = total_length + list->list[i].length; - } - - - candidates->length = 0; - if(total_length > (uint64_t)candidates->size) - { - candidates->size = total_length; - candidates->list = (k_mer_hit*)realloc(candidates->list, sizeof(k_mer_hit)*candidates->size); - candidates->tmp = (k_mer_hit*)realloc(candidates->tmp, sizeof(k_mer_hit)*candidates->size); - } - - uint64_t ID; - int flag; - ///while (flag = DeleteHeap(HBT, &x)) - while ((flag = DeleteHeap(HBT, &x))) - { - append_pos_to_Candidates_list(candidates, &x); - - i = HBT->index_i[x.ID]; - ID = x.ID; - - if (flag == 2) - { - for (; i < list->list[x.ID].length; i++) - { - y.ID = ID; - y.node.offset = list->list[x.ID].list[i].offset; - y.node.readID = list->list[x.ID].list[i].readID; - y.node.self_offset = list->list[x.ID].end_pos; - y.node.strand = list->list[x.ID].direction; - append_pos_to_Candidates_list(candidates, &y); - } - - break; - } - - if (i < list->list[x.ID].length) - { - y.ID = ID; - y.node.offset = list->list[x.ID].list[i].offset; - y.node.readID = list->list[x.ID].list[i].readID; - y.node.self_offset = list->list[x.ID].end_pos; - y.node.strand = list->list[x.ID].direction; - Insert_Heap(HBT, &y); - HBT->index_i[ID]++; - } - } -} - - - -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) - { - 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++) - { - kh_destroy(COUNT64, 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++) - { - kh_destroy(POS64, 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 -{ - long long* list; - uint64_t length; -} H_peaks; - -void insert_H_peaks(H_peaks* h, long long index, long long value) -{ - if(h->length <= (uint64_t)index) - { - long long newLen = index + 1; - h->list = (long long*)realloc(h->list, newLen*sizeof(long long)); - memset(h->list + h->length, 0, sizeof(long long) * (newLen - h->length)); - h->length = newLen; - } - - h->list[index] += value; -} - -inline void RC_Hash_code(Hash_code* code, Hash_code* rc_code, int k) -{ - rc_code->x[0] = 0; - rc_code->x[1] = 0; - int i; - for (i = 0; i < k; i++) - { - rc_code->x[0] = rc_code->x[0] << 1; - rc_code->x[1] = rc_code->x[1] << 1; - rc_code->x[0] |= (((uint64_t)((code->x[0] >> i) & 1))^((uint64_t)1)); - rc_code->x[1] |= (((uint64_t)((code->x[1] >> i) & 1))^((uint64_t)1)); - } -} - - - -void get_peak_debug(Total_Count_Table* TCB, long long* min, long long* max) -{ - int i; - Count_Table* h; - khint_t k; - long long c_count; - H_peaks LH; - LH.list = NULL; - LH.length = 0; - uint64_t sub_ID; - uint64_t sub_key; - Hash_code code, rc_code, debug_code; - char str[100]; - char rc_str[100]; - - - - 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 - { - sub_ID = i; - sub_key = kh_key(h, k); - - recover_hash_code(sub_ID, sub_key, &code, TCB->suffix_mode, - TCB->suffix_bits, asm_opt.k_mer_length); - RC_Hash_code(&code, &rc_code, asm_opt.k_mer_length); - RC_Hash_code(&rc_code, &debug_code, asm_opt.k_mer_length); - if(code.x[0] != debug_code.x[0] || code.x[1] != debug_code.x[1]) - { - fprintf(stderr, "error\n"); - } - - Hashcode_to_string(&code, str, asm_opt.k_mer_length); - Hashcode_to_string(&rc_code, rc_str, asm_opt.k_mer_length); - reverse_complement(str, asm_opt.k_mer_length); - if(memcmp(str, rc_str, asm_opt.k_mer_length) != 0) - { - fprintf(stderr, "error\n"); - int j; - for (j = 0; j < asm_opt.k_mer_length; j++) - { - fprintf(stderr, "%c",str[j]); - } - fprintf(stderr, "\n"); - - for (j = 0; j < asm_opt.k_mer_length; j++) - { - fprintf(stderr, "%c",rc_str[j]); - } - fprintf(stderr, "\n"); - - } - - - ///get_Total_Count_Table(&TCB, &k_code, k_mer_length); - - c_count = kh_value(h, k); - - if(get_Total_Count_Table(TCB, &code, asm_opt.k_mer_length) != c_count) - { - fprintf(stderr, "error\n"); - } - - insert_H_peaks(&LH, c_count, c_count); - } - } - } - - (*max) = -1; - (*min) = -1; - long long max_value = -1; - for (i = 0; i < (long long)LH.length; i++) - { - if(LH.list[i] >= max_value) - { - max_value = LH.list[i]; - (*max) = i; - } - } - - long long min_value = max_value; - for (i = 0; i < (long long)LH.length; i++) - { - if(LH.list[i] < min_value && LH.list[i] != 0) - { - min_value = LH.list[i]; - (*min) = i; - } - } - - for (i = 0; i < (long long)LH.length; i++) - { - ///fprintf(stderr, "%d, %d\n", i, LH.list[i]); - fprintf(stderr, "%lld\n", LH.list[i]); - } - - - - free(LH.list); -} - -///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; -} - -int get_total_freq(Total_Count_Table* TCB, uint64_t sub_ID, uint64_t sub_key, long long* T_count) -{ - Hash_code code, rc_code; - long long count, rc_count; - - recover_hash_code(sub_ID, sub_key, &code, TCB->suffix_mode, - TCB->suffix_bits, asm_opt.k_mer_length); - RC_Hash_code(&code, &rc_code, asm_opt.k_mer_length); - - count = get_Total_Count_Table(TCB, &code, asm_opt.k_mer_length); - rc_count = get_Total_Count_Table(TCB, &rc_code, asm_opt.k_mer_length); - (*T_count) = count + rc_count; - - if(count == 0) - { - return 0; - }///count > 0 && rc_count == 0 - else if(rc_count == 0) - { - return 1; - }///count > 0 && rc_count > 0 - else - { - int flag = cmp_Hash_code(&code, &rc_code); - - ///code > rc_code - if(flag > 0) - { - return 1; - }///code < rc_code - else if(flag < 0) - { - return 0; - } - else - { - (*T_count) = (*T_count)/2; - return 1; - } - } -} - -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 = kh_begin(h); 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); - - if(get_total_freq(TCB, sub_ID, sub_key, &count)==1) - { - 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 = kh_begin(h); 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 = 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 - { - sub_ID = i; - sub_key = kh_key(h, k); - get_total_freq(TCB, sub_ID, sub_key, &count); - - ///if (kh_value(h, k)>=k_mer_min_freq && kh_value(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_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; - } - } } void init_Chain_Data(Chain_Data* x) { - x->length = 0; - x->size = 0; - x->score = NULL; - x->pre = NULL; - x->indels = NULL; - x->self_length = NULL; + memset(x, 0, sizeof(Chain_Data)); } void clear_Chain_Data(Chain_Data* x) @@ -1656,101 +722,48 @@ void clear_Chain_Data(Chain_Data* x) x->length = 0; } - void destory_Chain_Data(Chain_Data* x) { free(x->score); free(x->pre); free(x->indels); free(x->self_length); + free(x->tmp); } - void resize_Chain_Data(Chain_Data* x, long long size) { - if(size > x->size) - { - x->size = size; - x->score = (long long*)realloc(x->score, x->size*sizeof(long long)); - x->pre = (long long*)realloc(x->pre, x->size*sizeof(long long)); - x->indels = (long long*)realloc(x->indels, x->size*sizeof(long long)); - x->self_length = (long long*)realloc(x->self_length, x->size*sizeof(long long)); - } + if (size + 1 > x->size) { + x->size = size + 1; + kroundup64(x->size); + REALLOC(x->score, x->size); + REALLOC(x->pre, x->size); + REALLOC(x->indels, x->size); + REALLOC(x->self_length, x->size); + REALLOC(x->tmp, x->size); + } } - void init_Candidates_list(Candidates_list* l) { l->length = 0; l->size = 0; l->list = NULL; - l->tmp = NULL; - l->foward_pos = 0; - l->rc_pos = 0; init_Chain_Data(&(l->chainDP)); } void clear_Candidates_list(Candidates_list* l) { l->length = 0; - l->foward_pos = 0; - l->rc_pos = 0; clear_Chain_Data(&(l->chainDP)); } void destory_Candidates_list(Candidates_list* l) { free(l->list); - free(l->tmp); destory_Chain_Data(&(l->chainDP)); } - - -int cmp_candidates_list(const void * a, const void * b) -{ - long long r_pos_a, r_pos_b; - - if((*((k_mer_hit*)a)).strand != (*((k_mer_hit*)b)).strand) - { - return (*((k_mer_hit*)a)).strand > (*((k_mer_hit*)b)).strand ? 1: -1; - } - else - { - if((*((k_mer_hit*)a)).readID != (*((k_mer_hit*)b)).readID) - { - return (*((k_mer_hit*)a)).readID > (*((k_mer_hit*)b)).readID ? 1: -1; - } - else - { - r_pos_a = (*((k_mer_hit*)a)).offset - (*((k_mer_hit*)a)).self_offset; - r_pos_b = (*((k_mer_hit*)b)).offset - (*((k_mer_hit*)b)).self_offset; - - if(r_pos_a != r_pos_b) - { - return r_pos_a > r_pos_b ? 1: -1; - } - else - { - if((*((k_mer_hit*)a)).self_offset != (*((k_mer_hit*)b)).self_offset) - { - return (*((k_mer_hit*)a)).self_offset > (*((k_mer_hit*)b)).self_offset ? 1: -1; - } - else - { - return 0; - } - - } - - } - - } -} - - - - void init_fake_cigar(Fake_Cigar* x) { x->buffer = NULL; @@ -1803,12 +816,10 @@ 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) { - if(size > x->size) - { + if (size > x->size) { x->size = size; - x->buffer = (uint64_t*)realloc(x->buffer, sizeof(uint64_t) * x->size); + REALLOC(x->buffer, x->size); } - x->length = 0; } @@ -1846,7 +857,5 @@ void resize_window_list_alloc(window_list_alloc* x, long long size) { x->buffer[i].error = -1; } - - x->length = 0; -} \ No newline at end of file +} diff --git a/Hash_Table.h b/Hash_Table.h index 4acd73a..d7820ae 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -1,13 +1,6 @@ #ifndef __HASHTABLE__ #define __HASHTABLE__ -#include "khash.h" -#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; +#include "htab.h" #define PREFIX_BITS 16 #define MAX_SUFFIX_BITS 64 @@ -37,26 +30,8 @@ typedef khash_t(POS64) Pos_Table; 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 -{ - uint64_t offset; - uint64_t readID; + uint32_t offset; + uint32_t readID:31, rev:1; } k_mer_pos; typedef struct @@ -70,17 +45,9 @@ typedef struct typedef struct { - k_mer_pos_list* list; - uint64_t size; - uint64_t length; -} k_mer_pos_list_alloc; - - -typedef struct -{ - int C_L[CIGAR_MAX_LENGTH]; - char C_C[CIGAR_MAX_LENGTH]; - int length; + int C_L[CIGAR_MAX_LENGTH]; + char C_C[CIGAR_MAX_LENGTH]; + int length; } CIGAR; typedef struct @@ -97,407 +64,96 @@ typedef struct CIGAR cigar; } window_list; - typedef struct { window_list* buffer; - long long length; - long long size; -}window_list_alloc; - + int32_t length; + int32_t size; +} window_list_alloc; typedef struct { uint64_t* buffer; - uint64_t length; - uint64_t size; -}Fake_Cigar; + uint32_t length; + uint32_t size; +} Fake_Cigar; typedef struct { - uint64_t x_id; + uint32_t x_id; ///the begining and end of the whole overlap - uint64_t x_pos_s; - uint64_t x_pos_e; - uint64_t x_pos_strand; + uint32_t x_pos_s; + uint32_t x_pos_e; + uint32_t x_pos_strand; - uint64_t y_id; - uint64_t y_pos_s; - uint64_t y_pos_e; - uint64_t y_pos_strand; + uint32_t y_id; + uint32_t y_pos_s; + uint32_t y_pos_e; + uint32_t y_pos_strand; - uint64_t overlapLen; - uint64_t shared_seed; - uint64_t align_length; - ///uint64_t total_errors; + uint32_t overlapLen; + uint32_t shared_seed; + uint32_t align_length; uint8_t is_match; uint8_t without_large_indel; - uint64_t non_homopolymer_errors; + int8_t strong; + uint32_t non_homopolymer_errors; window_list* w_list; - uint64_t w_list_size; - uint64_t w_list_length; - int8_t strong; + uint32_t w_list_size; + uint32_t w_list_length; Fake_Cigar f_cigar; window_list_alloc boundary_cigars; } overlap_region; - typedef struct { overlap_region* list; uint64_t size; uint64_t length; - ///uint64_t mapped_overlaps_length; - long long mapped_overlaps_length; + int64_t mapped_overlaps_length; } overlap_region_alloc; typedef struct { - ///uint64_t offset; - long long offset; - ///uint64_t self_offset; - long long self_offset; - uint64_t readID; - uint8_t strand; + uint32_t readID:31, strand:1; + uint32_t offset, self_offset; } k_mer_hit; - -typedef struct -{ - k_mer_hit node; - uint64_t ID; -} ElemType; - - -typedef struct -{ - ElemType* heap; - uint64_t* index_i; - int len; - int MaxSize; -} HeapSq; - -typedef struct -{ - long long* score; - long long* pre; - long long* indels; - long long* self_length; - long long length; - long long size; +typedef struct { + int32_t *score; + int64_t *pre; + int32_t *indels; + int32_t *self_length; + int64_t *tmp; // MUST BE 64-bit integer + int64_t length; + int64_t size; } Chain_Data; typedef struct { k_mer_hit* list; - k_mer_hit* tmp; long long length; long long size; - uint64_t foward_pos; - uint64_t rc_pos; 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; - - - - -inline uint64_t mod_d(uint64_t h_key, uint64_t low_key, uint64_t d) -{ - uint64_t result = (h_key >> 32) % d; - result = ((result << 32) + (h_key & (uint64_t)0xffffffff)) % d; - result = ((result << 32) + (low_key >> 32)) % d; - result = ((result << 32) + (low_key & (uint64_t)0xffffffff)) % d; - - return result; -} - - - -////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) -{ - uint64_t h_key, low_key; - h_key = low_key = 0; - - low_key = sub_ID << SAFE_SHIFT(suffix_bits); - low_key = low_key | sub_key; - - h_key = sub_ID >> (64 - suffix_bits); - - code->x[0] = code->x[1] = 0; - uint64_t mask = ALL >> (64 - k); - code->x[0] = low_key & mask; - - code->x[1] = h_key << (64 - k); - code->x[1] = code->x[1] | (low_key >> SAFE_SHIFT(k)); - - return 1; -} - -///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 might be 64，so it is unsafe - ///low_key = code->x[0] | (code->x[1] << k); - low_key = code->x[0] | (code->x[1] << SAFE_SHIFT(k)); - //k cannot be 0, so this shift is safe - h_key = code->x[1] >> (64 - k); - - if(mod_d(h_key, low_key, MODE_VALUE) > 3) - { - return 0; - } - - 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; - - return 1; -} - - -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 = kh_put(COUNT64, TCB->sub_h[sub_ID], sub_key, &absent); - if (absent) - { - kh_value(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_value(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; - - ///query hash table，key is k - t = kh_get(COUNT64, TCB->sub_h[sub_ID], sub_key); - - if (t != kh_end(TCB->sub_h[sub_ID])) - { - return kh_value(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; - - ///query hash table，key is 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].readID = readID|direction; - list[list[0].offset].offset = pos; - } - else - { - list[0].readID = readID; - ///list[0].readID = readID|direction; - list[0].offset = pos; - 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 merge_k_mer_pos_list_alloc_heap_sort(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT); - -void Init_Heap(HeapSq* HBT); -void destory_Heap(HeapSq* HBT); -void clear_Heap(HeapSq* HBT); - 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); @@ -505,14 +161,14 @@ 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) + +static inline long long y_start_offset(long long x_start, Fake_Cigar* o) { if(x_start == get_fake_gap_pos(o, o->length - 1)) { return get_fake_gap_shift(o, o->length - 1); } - long long i; for (i = 0; i < (long long)o->length; i++) { @@ -524,7 +180,7 @@ inline long long y_start_offset(long long x_start, Fake_Cigar* o) if(i == 0 || i == (long long)o->length) { - fprintf(stderr, "ERROR\n"); + fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__); exit(0); } @@ -532,24 +188,11 @@ inline long long y_start_offset(long long x_start, Fake_Cigar* o) return get_fake_gap_shift(o, i - 1); } -inline void print_fake_gap(Fake_Cigar* o) -{ - long long i; - for (i = 0; i < (long long)o->length; i++) - { - fprintf(stderr, "**i: %lld, gap_pos_in_x: %d, gap_shift: %d\n", - i, get_fake_gap_pos(o, i), - get_fake_gap_shift(o, i)); - } - -} - void resize_Chain_Data(Chain_Data* x, long long size); void init_window_list_alloc(window_list_alloc* x); 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); -void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, -double band_width_threshold, int max_skip, int x_readLen, int y_readLen); +void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, double band_width_threshold, int max_skip, int x_readLen, int y_readLen); -#endif \ No newline at end of file +#endif diff --git a/Makefile b/Makefile index e1ee15d..a0cff56 100644 --- a/Makefile +++ b/Makefile @@ -1,11 +1,12 @@ CXX= g++ -CXXFLAGS= -g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall #-fsanitize=address -fno-omit-frame-pointer#-Winline +CXXFLAGS= -g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall CPPFLAGS= INCLUDES= -OBJS= Output.o CommandLines.o Process_Read.o Assembly.o kmer.o Hash_Table.o \ - POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o #ksw2_extz2_sse.o +OBJS= Output.o CommandLines.o Process_Read.o Assembly.o Hash_Table.o \ + POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \ + htab.o hist.o sketch.o anchor.o sys.o EXE= hifiasm -LIBS= -lz -lpthread -lm #-fsanitize=address -fno-omit-frame-pointer +LIBS= -lz -lpthread -lm ifneq ($(asan),) CXXFLAGS+=-fsanitize=address @@ -30,28 +31,3 @@ depend: (LC_ALL=C; export LC_ALL; makedepend -Y -- $(CPPFLAGS) $(DFLAGS) -- *.cpp) # DO NOT DELETE - -Assembly.o: Assembly.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h -Assembly.o: CommandLines.h kmer.h Hash_Table.h khash.h POA.h Correct.h -Assembly.o: Levenshtein_distance.h Output.h Trio.h -CommandLines.o: CommandLines.h ketopt.h -Correct.o: Correct.h Hash_Table.h khash.h kmer.h Process_Read.h kseq.h -Correct.o: Overlaps.h kvec.h kdq.h CommandLines.h Levenshtein_distance.h -Correct.o: POA.h Assembly.h #ksw2.h -Hash_Table.o: Hash_Table.h khash.h kmer.h Process_Read.h kseq.h Overlaps.h -Hash_Table.o: kvec.h kdq.h CommandLines.h Correct.h Levenshtein_distance.h -Hash_Table.o: POA.h ksort.h -Levenshtein_distance.o: Levenshtein_distance.h -Output.o: Output.h CommandLines.h -Overlaps.o: Overlaps.h kvec.h kdq.h ksort.h Process_Read.h kseq.h -Overlaps.o: CommandLines.h Purge_Dups.h -POA.o: POA.h Hash_Table.h khash.h kmer.h Process_Read.h kseq.h Overlaps.h -POA.o: kvec.h kdq.h CommandLines.h Correct.h Levenshtein_distance.h -Process_Read.o: Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h -kmer.o: kmer.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h -main.o: CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h -main.o: Assembly.h Levenshtein_distance.h -Trio.o: Trio.h khashl.h kthread.h Process_Read.h CommandLines.h -kthread.o: kthread.h -Purge_Dups.o: Purge_Dups.h Overlaps.h Hash_Table.h Correct.h ksort.h -#ksw2_extz2_sse.o: ksw2.h \ No newline at end of file diff --git a/Overlaps.cpp b/Overlaps.cpp index 52edc3d..3aaaa1b 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -52,7 +52,6 @@ void ma_hit_sort_qns(ma_hit_t *a, long long n) radix_sort_hit_qns(a, a + n); } - void sort_kvec_t_u64_warp(kvec_t_u64_warp* u_vecs, uint32_t is_descend) { radix_sort_arch64(u_vecs->a.a, u_vecs->a.a + u_vecs->a.n); @@ -68,7 +67,6 @@ void sort_kvec_t_u64_warp(kvec_t_u64_warp* u_vecs, uint32_t is_descend) } } - asg_t *asg_init(void) { return (asg_t*)calloc(1, sizeof(asg_t)); @@ -4356,7 +4354,7 @@ uint32_t startNode, uint32_t endNode) } else { - fprintf(stderr, "ERROR\n"); + fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__); } if(asg_arc_a(g, N_list[2])[0].v == (N_list[0]^1)) @@ -4369,7 +4367,7 @@ uint32_t startNode, uint32_t endNode) } else { - fprintf(stderr, "ERROR\n"); + fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__); } if(N_list[3] != N_list[4]) @@ -4920,12 +4918,10 @@ int asg_arc_del_short_diploid_unclean(asg_t *g, float drop_ratio, ma_hit_t_alloc uint32_t detect_single_path_with_dels(asg_t *g, uint32_t begNode, uint32_t* endNode, long long* Len, buf_t* b) { - uint32_t v = begNode, w; uint32_t kv, kw; (*Len) = 0; - while (1) { (*Len)++; @@ -5594,6 +5590,7 @@ ma_hit_t_alloc* reverse_sources, long long min_edge_length, R_to_U* ruIndex) if(flag1 == LOOP || flag2 == LOOP) { + free(b_0.b.a); free(b_1.b.a); return -1; } @@ -5612,6 +5609,7 @@ ma_hit_t_alloc* reverse_sources, long long min_edge_length, R_to_U* ruIndex) if(l1 <= min_edge_length || l2 <= min_edge_length) { + free(b_0.b.a); free(b_1.b.a); return -1; } @@ -9051,7 +9049,7 @@ ma_sub_t *coverage_cut, int max_hang, int min_ovlp) break; } } - if(k == nv) fprintf(stderr, "ERROR\n"); + if(k == nv) fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__); av = asg_arc_a(read_g, v); nv = asg_arc_n(read_g, v); @@ -9065,7 +9063,7 @@ ma_sub_t *coverage_cut, int max_hang, int min_ovlp) } } - if(k == nv) fprintf(stderr, "ERROR\n"); + if(k == nv) fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__); if(pE->el == 1 && aE->el == 1) continue; @@ -9146,7 +9144,7 @@ ma_sub_t *coverage_cut, int max_hang, int min_ovlp) } l = asg_arc_len(t_f); } - if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n"); + if(l == (uint32_t)-1) fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__); /*******************************for debug************************************/ @@ -9195,7 +9193,7 @@ ma_sub_t *coverage_cut, int max_hang, int min_ovlp) } l = asg_arc_len(t_f); } - if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n"); + if(l == (uint32_t)-1) fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__); /*******************************for debug************************************/ @@ -11671,7 +11669,7 @@ void print_untig(ma_ug_t *g, uint32_t uId, const char* info, uint32_t is_print_r for (k = 0; k < u->n; k++) { fprintf(stderr, "%s: rId>>1: %lu, dir: %lu, name: %.*s\n", - info, u->a[k]>>33, (u->a[k]>>32)&1, + info, (unsigned long)(u->a[k]>>33), (unsigned long)((u->a[k]>>32)&1), (int)Get_NAME_LENGTH((R_INF), (u->a[k]>>33)), Get_NAME((R_INF), (u->a[k]>>33))); } } @@ -13640,8 +13638,6 @@ int load_ma_hit_ts(ma_hit_t_alloc** x, char* read_file_name) } - - void write_ma(ma_hit_t* x, FILE* fp) { fwrite(&(x->qns), sizeof(x->qns), 1, fp); @@ -13708,35 +13704,28 @@ All_reads *RNF, char* output_file_name) free(gfa_name); } -int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, -char* output_file_name) +int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, char* output_file_name) { - char* gfa_name = (char*)malloc(strlen(output_file_name)+25); - sprintf(gfa_name, "%s.ovlp", output_file_name); - if(!load_All_reads(&R_INF, gfa_name)) - { - return 0; - } - - sprintf(gfa_name, "%s.ovlp.source", output_file_name); - if(!load_ma_hit_ts(sources, gfa_name)) - { - return 0; - } - - - sprintf(gfa_name, "%s.ovlp.reverse", output_file_name); - if(!load_ma_hit_ts(reverse_sources, gfa_name)) - { - return 0; - } - - return 1; + char* gfa_name = (char*)malloc(strlen(output_file_name)+25); + sprintf(gfa_name, "%s.ovlp", output_file_name); + if (!load_All_reads(&R_INF, gfa_name)) { + free(gfa_name); + return 0; + } + sprintf(gfa_name, "%s.ovlp.source", output_file_name); + if (!load_ma_hit_ts(sources, gfa_name)) { + free(gfa_name); + return 0; + } + sprintf(gfa_name, "%s.ovlp.reverse", output_file_name); + if (!load_ma_hit_ts(reverse_sources, gfa_name)) { + free(gfa_name); + return 0; + } + free(gfa_name); + return 1; } - - - // count the number of outgoing arcs, excluding reduced arcs static inline int count_out(const asg_t *g, uint32_t v) { @@ -21006,7 +20995,7 @@ void merge_ug_nodes(ma_ug_t *ug, asg_t* read_g, kvec_t_u64_warp* array) if(aw[i].del) continue; if(aw[i].v == (v^1)) break; } - if(i == nw) fprintf(stderr, "ERROR\n"); + if(i == nw) fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__); kmp = kmp | (uint64_t)(aw[i].ol); @@ -21037,7 +21026,7 @@ void merge_ug_nodes(ma_ug_t *ug, asg_t* read_g, kvec_t_u64_warp* array) if(aw[i].del) continue; if(aw[i].v == (v^1)) break; } - if(i == nw) fprintf(stderr, "ERROR\n"); + if(i == nw) fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__); kmp = kmp | (uint64_t)(aw[i].ol); @@ -26198,7 +26187,7 @@ void fix_binned_reads(ma_hit_t_alloc* paf, uint64_t n_read, ma_sub_t* coverage_c } fprintf(stderr, "n_read: %lu, binned_reads: %lu, binned_error_reads: %lu\n", - n_read, binned_reads, binned_error_reads); + (unsigned long)n_read, (unsigned long)binned_reads, (unsigned long)binned_error_reads); } @@ -26301,24 +26290,22 @@ int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long n_read, uint64_t* readLen, long long mini_overlap_length, long long max_hang_length, long long clean_round, long long gap_fuzz, float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, -long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t* sg, -ma_sub_t* coverage_cut, int debug_g) +long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr, +ma_sub_t **coverage_cut_ptr, int debug_g) { + ma_sub_t *coverage_cut = *coverage_cut_ptr; + asg_t *sg = *sg_ptr; if(debug_g) goto debug_gfa; - ///just for debug renew_graph_init(sources, reverse_sources, sg, coverage_cut, ruIndex, n_read); - - ///it's hard to say which function is better ///normalize_ma_hit_t_single_side(sources, n_read); normalize_ma_hit_t_single_side_advance(sources, n_read); normalize_ma_hit_t_single_side_advance(reverse_sources, n_read); - // debug_info_of_specfic_read(">m64011_190830_220126/117834372/ccs", sources, reverse_sources, // -1, "clean"); @@ -26532,13 +26519,13 @@ ma_sub_t* coverage_cut, int debug_g) bubble_dist, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length); - output_contig_graph_alternative(sg, coverage_cut, output_file_name, sources, max_hang_length, mini_overlap_length); } + + *coverage_cut_ptr = coverage_cut; + *sg_ptr = sg; } - - void build_string_graph_without_clean( int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long n_read, uint64_t* readLen, long long mini_overlap_length, @@ -26551,8 +26538,7 @@ long long bubble_dist, int read_graph, int write) asg_t *sg = NULL; ma_sub_t* coverage_cut = NULL; - // debug_info_of_specfic_read("m64011_190329_072846/80545633/ccs", - // sources, reverse_sources, -1, "clean"); + // debug_info_of_specfic_read("m64011_190329_072846/80545633/ccs", sources, reverse_sources, -1, "clean"); ///actually min_thres = asm_opt.max_short_tip + 1 there are asm_opt.max_short_tip reads min_thres = asm_opt.max_short_tip + 1; @@ -26565,7 +26551,7 @@ long long bubble_dist, int read_graph, int write) clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length, max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio, - output_file_name, bubble_dist, read_graph, &ruIndex, sg, coverage_cut, 1); + output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 1); asg_destroy(sg); free(coverage_cut); destory_R_to_U(&ruIndex); @@ -26573,20 +26559,17 @@ long long bubble_dist, int read_graph, int write) } } - - if (asm_opt.write_index_to_disk && write) { write_all_data_to_disk(sources, reverse_sources, &R_INF, output_file_name); } - try_rescue_overlaps(sources, reverse_sources, n_read, 4); clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length, max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio, - output_file_name, bubble_dist, read_graph, &ruIndex, sg, coverage_cut, 0); + output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 0); asg_destroy(sg); free(coverage_cut); diff --git a/Overlaps.h b/Overlaps.h index 4e6bfce..03c5d09 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -47,9 +47,6 @@ #define CUT_DIF_HAP 12 - - - ///query is the read itself typedef struct { uint64_t qns; @@ -1057,4 +1054,4 @@ int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio); #define JUNK_COV 5 #define DISCARD_RATE 0.8 -#endif \ No newline at end of file +#endif diff --git a/POA.cpp b/POA.cpp index c519693..d5eae20 100644 --- a/POA.cpp +++ b/POA.cpp @@ -1,29 +1,27 @@ -#include "POA.h" #include +#include +#include "POA.h" #include "Correct.h" #include "Process_Read.h" #define INIT_EDGE_SIZE 50 #define INCREASE_EDGE_SIZE 5 #define INIT_NODE_SIZE 16000 - - +/******** + * Edge * + ********/ void init_Edge_alloc(Edge_alloc* list) { - if (list->list == NULL) - { - list->size = INIT_EDGE_SIZE; - list->length = 0; - list->delete_length = 0; - list->list = (Edge*)malloc(sizeof(Edge)*list->size); - } - else - { - list->length = 0; - list->delete_length = 0; - } - + if (list->list == NULL) { + list->size = INIT_EDGE_SIZE; + list->length = 0; + list->delete_length = 0; + list->list = (Edge*)malloc(sizeof(Edge)*list->size); + } else { + list->length = 0; + list->delete_length = 0; + } } void clear_Edge_alloc(Edge_alloc* list) @@ -34,28 +32,29 @@ void clear_Edge_alloc(Edge_alloc* list) void destory_Edge_alloc(Edge_alloc* list) { - free(list->list); + if (list && list->list) + free(list->list); } void append_Edge_alloc(Edge_alloc* list, uint64_t in_node, uint64_t out_node, uint64_t weight, uint64_t length) { - if (list->length + 1 > list->size) - { - list->size = list->size + INCREASE_EDGE_SIZE; - list->list = (Edge*)realloc(list->list, sizeof(Edge)*list->size); - } + if (list->length + 1 > list->size) { + uint64_t old_size = list->size; + list->size = list->size + INCREASE_EDGE_SIZE; + list->list = (Edge*)realloc(list->list, sizeof(Edge)*list->size); + memset(&list->list[old_size], 0, (list->size - old_size) * sizeof(Edge)); + } - list->list[list->length].in_node = in_node; - list->list[list->length].out_node = out_node; - list->list[list->length].weight = weight; - list->list[list->length].length = length; - list->list[list->length].num_insertions = 0; - list->list[list->length].self_edge_ID = list->length; + list->list[list->length].in_node = in_node; + list->list[list->length].out_node = out_node; + list->list[list->length].weight = weight; + list->list[list->length].length = length; + list->list[list->length].num_insertions = 0; + list->list[list->length].self_edge_ID = list->length; - list->length++; + list->length++; } - int add_and_check_bi_direction_edge(Graph* graph, Node* in_node, Node* out_node, uint64_t weight, uint64_t flag) { Edge* e_forward; @@ -81,7 +80,6 @@ int add_and_check_bi_direction_edge(Graph* graph, Node* in_node, Node* out_node, } } - void add_bi_direction_edge(Graph* graph, Node* in_node, Node* out_node, uint64_t weight, uint64_t flag) { @@ -95,8 +93,6 @@ void add_bi_direction_edge(Graph* graph, Node* in_node, Node* out_node, uint64_t = Output_Edges((*in_node)).length - 1; } - - int remove_and_check_bi_direction_edge_from_nodes(Graph* graph, Node* in_node, Node* out_node) { Edge* e_forward; @@ -135,8 +131,6 @@ int remove_and_check_bi_direction_edge_from_nodes(Graph* graph, Node* in_node, N } } - - int remove_and_check_bi_direction_edge_from_edge(Graph* graph, Edge* e) { Edge* e_forward; @@ -171,109 +165,71 @@ int remove_and_check_bi_direction_edge_from_edge(Graph* graph, Edge* e) { return 0; } - - } - +/******** + * Node * + ********/ void init_Node_alloc(Node_alloc* list) { - list->size = INIT_NODE_SIZE; - list->length = 0; - list->delete_length = 0; - list->list = (Node*)malloc(sizeof(Node)*list->size); - list->sort.size = 0; - list->sort.list = NULL; - list->sort.visit = NULL; - - list->sort.iterative_buffer = NULL; - list->sort.iterative_buffer_visit = NULL; - - - uint64_t i; - for (i = 0; i < list->size; i++) - { - list->list[i].insertion_edges.list=NULL; - list->list[i].mismatch_edges.list=NULL; - list->list[i].deletion_edges.list=NULL; - } + memset(list, 0, sizeof(Node_alloc)); + list->size = INIT_NODE_SIZE; + list->list = (Node*)calloc(list->size, sizeof(Node)); } void destory_Node_alloc(Node_alloc* list) { - uint64_t i =0; - for (i = 0; i < list->length; i++) - { - destory_Edge_alloc(&list->list[i].deletion_edges); - destory_Edge_alloc(&list->list[i].insertion_edges); - destory_Edge_alloc(&list->list[i].mismatch_edges); - } - - free(list->list); - free(list->sort.list); - free(list->sort.visit); - free(list->sort.iterative_buffer); - free(list->sort.iterative_buffer_visit); - ///free(list->topo_order); + uint64_t i; + for (i = 0; i < list->size; i++) { + destory_Edge_alloc(&list->list[i].deletion_edges); + destory_Edge_alloc(&list->list[i].insertion_edges); + destory_Edge_alloc(&list->list[i].mismatch_edges); + } + free(list->list); + free(list->sort.list); + free(list->sort.visit); + free(list->sort.iterative_buffer); + free(list->sort.iterative_buffer_visit); } void clear_Node_alloc(Node_alloc* list) { - uint64_t i =0; - for (i = 0; i < list->length; i++) - { - clear_Edge_alloc(&list->list[i].insertion_edges); - clear_Edge_alloc(&list->list[i].mismatch_edges); - clear_Edge_alloc(&list->list[i].deletion_edges); - } - - list->length = 0; - list->delete_length = 0; + uint64_t i =0; + for (i = 0; i < list->length; i++) { // TODO: is this list->size or list->length? The original version is list->length. + clear_Edge_alloc(&list->list[i].insertion_edges); + clear_Edge_alloc(&list->list[i].mismatch_edges); + clear_Edge_alloc(&list->list[i].deletion_edges); + } + list->length = 0; + list->delete_length = 0; } - uint64_t append_Node_alloc(Node_alloc* list, char base) { - - if (list->length + 1 > list->size) - { - uint64_t i = list->size; + if (list->length + 1 > list->size) { + uint64_t old_size = list->size; + list->size = list->size * 2; + list->list = (Node*)realloc(list->list, sizeof(Node) * list->size); + memset(&list->list[old_size], 0, (list->size - old_size) * sizeof(Node)); + } - list->size = list->size * 2; - list->list = (Node*)realloc(list->list, sizeof(Node)*list->size); - ///list->topo_order = (uint64_t*)realloc(list->topo_order, sizeof(uint64_t)*list->size); - - for (; i < list->size; i++) - { - list->list[i].deletion_edges.list=NULL; - list->list[i].insertion_edges.list=NULL; - list->list[i].mismatch_edges.list=NULL; - } - } - - list->list[list->length].ID = list->length; - list->list[list->length].base = base; - list->list[list->length].weight = 1; - list->list[list->length].num_insertions = 0; - init_Edge_alloc(&list->list[list->length].deletion_edges); - init_Edge_alloc(&list->list[list->length].insertion_edges); - init_Edge_alloc(&list->list[list->length].mismatch_edges); - - list->length++; + list->list[list->length].ID = list->length; + list->list[list->length].base = base; + list->list[list->length].weight = 1; + list->list[list->length].num_insertions = 0; + init_Edge_alloc(&list->list[list->length].deletion_edges); + init_Edge_alloc(&list->list[list->length].insertion_edges); + init_Edge_alloc(&list->list[list->length].mismatch_edges); - return list->length - 1; + list->length++; + + return list->length - 1; } - - - - - - - - - +/********* + * Graph * + *********/ void init_Graph(Graph* g) { @@ -310,12 +266,6 @@ void clear_Graph(Graph* g) clear_Queue(&(g->node_q)); } - - - - - - void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length, long long* startID, long long* endID) { long long firstID, lastID, nodeID, i; @@ -356,8 +306,6 @@ void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length, } - - void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_string, long long x_length, char* y_string, long long y_length, CIGAR* cigar, long long backbone_start, long long backbone_end) { @@ -422,9 +370,3 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri cigar_i++; } } - - - - - - diff --git a/Process_Read.cpp b/Process_Read.cpp index 2b4d3f4..e8cc448 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -1,25 +1,8 @@ -#include "Process_Read.h" #include #include #include #include -#include - - -gz_files fps; - -R_buffer RDB; -static uint64_t total_reads; - -pthread_mutex_t i_readinputMutex; -pthread_mutex_t i_queueMutex; -pthread_mutex_t i_terminateMutex; -pthread_cond_t i_flushCond; -pthread_cond_t i_readinputflushCond; -pthread_cond_t i_stallCond; -pthread_cond_t i_readinputstallCond; -pthread_mutex_t i_doneMutex; - +#include "Process_Read.h" uint8_t seq_nt6_table[256] = { 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, @@ -45,34 +28,28 @@ char bit_t_seq_table_rc[256][4] = {{0}}; char s_H[5] = {'A', 'C', 'G', 'T', 'N'}; char rc_Table[5] = {'T', 'G', 'C', 'A', 'N'}; - void init_All_reads(All_reads* r) { + memset(r, 0, sizeof(All_reads)); r->index_size = READ_INIT_NUMBER; r->read_length = (uint64_t*)malloc(sizeof(uint64_t)*r->index_size); - r->read_sperate = NULL; - r->N_site = 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; - r->trio_flag = NULL; } void destory_All_reads(All_reads* r) { uint64_t i = 0; - for (i = 0; i < r->total_reads; i++) - { + for (i = 0; i < r->total_reads; i++) { if (r->N_site[i] != NULL) - { free(r->N_site[i]); - } free(r->read_sperate[i]); + if (r->paf) free(r->paf[i].buffer); + if (r->reverse_paf) free(r->reverse_paf[i].buffer); } + free(r->paf); + free(r->reverse_paf); free(r->N_site); free(r->read_sperate); free(r->name); @@ -81,7 +58,6 @@ void destory_All_reads(All_reads* r) free(r->trio_flag); } - void write_All_reads(All_reads* r, char* read_file_name) { fprintf(stderr, "Writing reads to disk... \n"); @@ -112,9 +88,6 @@ void write_All_reads(All_reads* r, char* read_file_name) { fwrite(&zero, sizeof(zero), 1, fp); } - - - } fwrite(r->read_length, sizeof(uint64_t), r->total_reads, fp); @@ -132,16 +105,13 @@ void write_All_reads(All_reads* r, char* read_file_name) fprintf(stderr, "Reads has been written.\n"); } - - int load_All_reads(All_reads* r, char* read_file_name) { - fprintf(stderr, "Loading reads from disk... \n"); char* index_name = (char*)malloc(strlen(read_file_name)+15); sprintf(index_name, "%s.bin", read_file_name); FILE* fp = fopen(index_name, "r"); - if (!fp) - { + if (!fp) { + free(index_name); return 0; } int local_adapterLen; @@ -164,12 +134,10 @@ int load_All_reads(All_reads* r, char* read_file_name) r->N_site = (uint64_t**)malloc(sizeof(uint64_t*)*r->total_reads); for (i = 0; i < r->total_reads; i++) { - f_flag += fread(&zero, sizeof(zero), 1, fp); if (zero) { - r->N_site[i] = (uint64_t*)malloc(sizeof(uint64_t)*(zero + 1)); r->N_site[i][0] = zero; if (r->N_site[i][0]) @@ -181,7 +149,6 @@ int load_All_reads(All_reads* r, char* read_file_name) { r->N_site[i] = NULL; } - } r->read_length = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads); @@ -211,8 +178,6 @@ int load_All_reads(All_reads* r, char* read_file_name) r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads); r->second_round_cigar = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads); - r->paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads); - r->reverse_paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads); for (i = 0; i < r->total_reads; i++) { r->second_round_cigar[i].size = r->cigars[i].size = 0; @@ -222,8 +187,6 @@ int load_All_reads(All_reads* r, char* read_file_name) r->second_round_cigar[i].lost_base_size = r->cigars[i].lost_base_size = 0; r->second_round_cigar[i].lost_base_length = r->cigars[i].lost_base_length = 0; r->second_round_cigar[i].lost_base = r->cigars[i].lost_base = NULL; - init_ma_hit_t_alloc(&(r->paf[i])); - init_ma_hit_t_alloc(&(r->reverse_paf[i])); } free(index_name); @@ -233,32 +196,28 @@ int load_All_reads(All_reads* r, char* read_file_name) return 1; } - - -inline void insert_read(All_reads* r, kstring_t* read, kstring_t* name) +void ha_insert_read_len(All_reads *r, int read_len, int name_len) { r->total_reads++; - r->total_reads_bases = r->total_reads_bases + read->l; - r->total_name_length = r->total_name_length + name->l; + r->total_reads_bases += (uint64_t)read_len; + r->total_name_length += (uint64_t)name_len; - ///must +1 - if (r->index_size < r->total_reads + 2) - { + // must +1 + if (r->index_size < r->total_reads + 2) { r->index_size = r->index_size * 2 + 2; - r->read_length = (uint64_t*)realloc(r->read_length,sizeof(uint64_t)*(r->index_size)); + r->read_length = (uint64_t*)realloc(r->read_length, 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->name_index = (uint64_t*)realloc(r->name_index, sizeof(uint64_t) * r->name_index_size); } - r->read_length[r->total_reads - 1] = read->l; - r->name_index[r->total_reads] = r->name_index[r->total_reads-1] + name->l; + r->read_length[r->total_reads - 1] = read_len; + r->name_index[r->total_reads] = r->name_index[r->total_reads - 1] + name_len; } void malloc_All_reads(All_reads* r) { - r->read_size = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads); - memcpy (r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads); + memcpy(r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads); r->read_sperate = (uint8_t**)malloc(sizeof(uint8_t*)*r->total_reads); long long i = 0; @@ -313,7 +272,6 @@ void init_aux_table() bit_t_seq_table_rc[i][2] = RC_CHAR(bit_t_seq_table[i][1]); bit_t_seq_table_rc[i][3] = RC_CHAR(bit_t_seq_table[i][0]); } - } } @@ -338,18 +296,12 @@ void init_UC_Read(UC_Read* r) bit_t_seq_table_rc[i][2] = RC_CHAR(bit_t_seq_table[i][1]); bit_t_seq_table_rc[i][3] = RC_CHAR(bit_t_seq_table[i][0]); } - } - } - -void recover_UC_Read_sub_region_begin_end -(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID, int extra_begin, int extra_end) +void recover_UC_Read_sub_region_begin_end(char* r, long long start_pos, long long length, uint8_t strand, + All_reads* R_INF, long long ID, int extra_begin, int extra_end) { - - - long long readLen = Get_READ_LENGTH((*R_INF), ID); uint8_t* src = Get_READ((*R_INF), ID); @@ -357,17 +309,11 @@ void recover_UC_Read_sub_region_begin_end long long copyLen; long long end_pos = start_pos + length - 1; - - - if (strand == 0) { - i = start_pos; copyLen = 0; - - long long initLen = start_pos % 4; if (initLen != 0) @@ -376,8 +322,6 @@ void recover_UC_Read_sub_region_begin_end copyLen = copyLen + 4 - initLen; i = i + copyLen; } - - while (copyLen < length) { memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4); @@ -385,7 +329,6 @@ void recover_UC_Read_sub_region_begin_end i = i + 4; } - if (R_INF->N_site[ID]) { for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) @@ -400,17 +343,12 @@ void recover_UC_Read_sub_region_begin_end } } } - - } else { - start_pos = readLen - start_pos - 1; end_pos = readLen - end_pos - 1; - - ///start_pos > end_pos i = start_pos; copyLen = 0; @@ -436,7 +374,6 @@ void recover_UC_Read_sub_region_begin_end for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) { - if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos) { r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N'; @@ -447,21 +384,11 @@ void recover_UC_Read_sub_region_begin_end } } } - - - } - } - - - void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID) { - - - long long readLen = Get_READ_LENGTH((*R_INF), ID); uint8_t* src = Get_READ((*R_INF), ID); @@ -471,7 +398,6 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, if (strand == 0) { - i = start_pos; copyLen = 0; @@ -483,8 +409,7 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, copyLen = copyLen + 4 - initLen; i = i + copyLen; } - - + while (copyLen < length) { memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4); @@ -492,7 +417,6 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, i = i + 4; } - if (R_INF->N_site[ID]) { for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) @@ -507,17 +431,12 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, } } } - - } else { - start_pos = readLen - start_pos - 1; end_pos = readLen - end_pos - 1; - - ///start_pos > end_pos i = start_pos; copyLen = 0; @@ -543,7 +462,6 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) { - if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos) { r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N'; @@ -554,15 +472,11 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, } } } - - - } - } -void recover_UC_Read(UC_Read* r, All_reads* R_INF, uint64_t ID) +void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID) { r->length = Get_READ_LENGTH((*R_INF), ID); uint8_t* src = Get_READ((*R_INF), ID); @@ -623,7 +537,6 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID) index = index + 4; } - if (R_INF->N_site[ID]) { for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++) @@ -631,11 +544,8 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID) r->seq[r->length - R_INF->N_site[ID][i] - 1] = 'N'; } } - } - - #define COMPRESS_BASE {c = seq_nt6_table[(uint8_t)src[i]];\ if (c >= 4)\ {\ @@ -645,9 +555,8 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID) }\ i++;}\ -void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ) +void ha_compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ) { - ///N_site_lis saves the pos of all Ns in this read ///N_site_lis[0] is the number of Ns if (N_site_occ) @@ -666,10 +575,8 @@ void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_l uint8_t tmp = 0; uint8_t c = 0; - while (i + 4 <= src_l) { - tmp = 0; COMPRESS_BASE; @@ -705,351 +612,8 @@ void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_l dest[dest_i] = tmp; dest_i++; } - } - - - -void open_file(gz_files* nfps, char* name) -{ - nfps->fp = gzopen(name, "r"); - if(nfps->fp == 0) - { - fprintf(stderr, "[ERROR] Cannot find the input file: %s\n", name); - exit(0); - } - nfps->seq = kseq_init(nfps->fp); -} - - -void close_file(gz_files* nfps) -{ - kseq_destroy(nfps->seq); - gzclose(nfps->fp); -} - -void init_gz_files(hifiasm_opt_t* asm_opt) -{ - fps.idx = 0; - fps.num_reads = asm_opt->num_reads; - fps.reads = asm_opt->read_file_names; - fps.seq = NULL; - fps.fp = NULL; - if(fps.num_reads > 0) - { - open_file(&fps, fps.reads[fps.idx]); - fps.idx++; - } -} - -void destory_gz_files() -{ - close_file(&fps); -} - -int read_item() -{ - int l = kseq_read(fps.seq); - if(l >= 0 || (l < 0 && fps.idx >= fps.num_reads)) - { - return l; - } - ///l < 0 && fps.idx < fps.num_reads - close_file(&fps); - open_file(&fps, fps.reads[fps.idx]); - fps.idx++; - return read_item(); -} - - -inline void exchage_kstring_t(kstring_t* a, kstring_t* b) -{ - kstring_t tmp; - tmp = *a; - *a = *b; - *b = tmp; -} - -int get_read(kseq_t *s, int adapterLen) -{ - int l; - - ///if ((l = kseq_read(seq)) >= 0) - if ((l = read_item()) >= 0) - { - - exchage_kstring_t(&(fps.seq->comment), &s->comment); - exchage_kstring_t(&(fps.seq->name), &s->name); - exchage_kstring_t(&(fps.seq->qual), &s->qual); - exchage_kstring_t(&(fps.seq->seq), &s->seq); - - if(adapterLen > 0) - { - if((int)s->seq.l <= adapterLen*2) - { - s->seq.l = 0; - } - else - { - long long i; - for (i = 0; i < ((int)s->seq.l - adapterLen*2); i++) - { - s->seq.s[i] = s->seq.s[i + adapterLen]; - } - s->seq.l -= adapterLen*2; - } - - } - - return 1; - } - else - { - return 0; - } - - -} - -void init_R_buffer_block(R_buffer_block* curr_sub_block) -{ - curr_sub_block->read = (kseq_t*)calloc(RDB.block_inner_size, sizeof(kseq_t)); - 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; - RDB.num = 0; - RDB.block_inner_size = READ_BLOCK_SIZE; - RDB.size = thread_num*READ_BLOCK_NUM_PRE_THR; - - RDB.sub_block = (R_buffer_block*)malloc(sizeof(R_buffer_block)*RDB.size); - - int i = 0; - - for (i = 0; i < RDB.size; i++) - { - init_R_buffer_block(&RDB.sub_block[i]); - } - -} - - -void destory_R_buffer_block(R_buffer_block* curr_sub_block) -{ - kseq_destroy(curr_sub_block->read); -} - - -void destory_R_buffer() -{ - int i = 0; - - for (i = 0; i < RDB.size; i++) - { - destory_R_buffer_block(&RDB.sub_block[i]); - } - - free(RDB.sub_block); - -} - - -inline void load_read_block(R_buffer_block* read_batch, int batch_read_size, - int* return_file_flag, int is_insert, int adapterLen) -{ - int inner_i = 0; - int file_flag = 1; - - - - - while (inner_iread[inner_i], adapterLen); - - 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) - { - break; - } - } - - if (inner_i || file_flag) - { - file_flag = 1; - } - - *return_file_flag = file_flag; - read_batch->num = inner_i; - -} - - -inline void push_R_block(R_buffer_block* tmp_sub_block) -{ - - - ///only exchange pointers - kseq_t *k1; - k1 = RDB.sub_block[RDB.num].read; - - RDB.sub_block[RDB.num].read = tmp_sub_block->read; - - tmp_sub_block->read = k1; - - RDB.sub_block[RDB.num].num = tmp_sub_block->num; - tmp_sub_block->num = 0; - - RDB.num++; -} - - -inline void pop_R_block(R_buffer_block* curr_sub_block) -{ - RDB.num--; - - ///only exchange pointers - kseq_t *k1; - k1 = RDB.sub_block[RDB.num].read; - - RDB.sub_block[RDB.num].read = curr_sub_block->read; - - curr_sub_block->read = k1; - - curr_sub_block->num = RDB.sub_block[RDB.num].num; - RDB.sub_block[RDB.num].num = 0; - - - -} - - - -void* input_reads_muti_threads(void* arg) -{ - int is_insert = *((int*)arg); - - - total_reads = 0; - - - int file_flag = 1; - - R_buffer_block tmp_buf; - - init_R_buffer_block(&tmp_buf); - - - - while (1) - { - load_read_block(&tmp_buf, RDB.block_inner_size, &file_flag, is_insert, asm_opt.adapterLen); - - if (file_flag == 0) - { - break; - } - - - pthread_mutex_lock(&i_readinputMutex); - while (IS_FULL(RDB)) - { - - pthread_cond_signal(&i_readinputstallCond); - pthread_cond_wait(&i_readinputflushCond, &i_readinputMutex); - } - - - push_R_block(&tmp_buf); - - pthread_cond_signal(&i_readinputstallCond); - pthread_mutex_unlock(&i_readinputMutex); - } - - - pthread_mutex_lock(&i_readinputMutex); - RDB.all_read_end = 1; - pthread_cond_signal(&i_readinputstallCond); //important - pthread_mutex_unlock(&i_readinputMutex); - - destory_R_buffer_block(&tmp_buf); - - fprintf(stderr, "Reads #: %lu\n", (unsigned long)total_reads); - fprintf(stderr, "Bases #: %lu\n", (unsigned long)R_INF.total_reads_bases); - - - return NULL; -} - - - -int get_reads_mul_thread(R_buffer_block* curr_sub_block) -{ - - - pthread_mutex_lock(&i_readinputMutex); - - - while (IS_EMPTY(RDB) && RDB.all_read_end == 0) - { - - pthread_cond_signal(&i_readinputflushCond); - pthread_cond_wait(&i_readinputstallCond, &i_readinputMutex); - } - - - if (!IS_EMPTY(RDB)) - { - pop_R_block(curr_sub_block); - pthread_cond_signal(&i_readinputflushCond); - pthread_mutex_unlock(&i_readinputMutex); - - - return 1; - } - else - { - curr_sub_block->num = 0; - - pthread_cond_signal(&i_readinputstallCond); //important - - pthread_mutex_unlock(&i_readinputMutex); - - return 0; - } - - -} - - - - - - - - void reverse_complement(char* pattern, uint64_t length) { uint64_t i = 0; @@ -1069,63 +633,4 @@ void reverse_complement(char* pattern, uint64_t length) { pattern[end] = RC_CHAR(pattern[end]); } - -} - - - -typedef struct { - char* tmp; - long long tmpSize; - char* dest; - long long destSize; - FILE* fp; -} LineReader; - -int get_single_line(LineReader* line) -{ - long long currentLen = 0, getLen = 0; - line->tmp[0] = '\0'; - while (fgets(line->tmp, line->tmpSize, line->fp) != NULL) - { - getLen = strlen(line->tmp); - if(getLen + currentLen >= line->destSize) - { - line->destSize = getLen + currentLen + 1; - line->dest = (char*)realloc(line->dest, line->destSize); - } - memcpy(line->dest + currentLen, line->tmp, getLen+1); - currentLen = currentLen + getLen; - if(currentLen > 0 && line->dest[currentLen - 1] == '\n') - { - return 1; - } - } - - if(currentLen > 0) - { - return 1; - } - else - { - return 0; - } -} - -void get_trio_info(char* input, uint8_t* pm) -{ - uint32_t i; - (*pm) = AMBIGU; - for (i = 0; input[i] != '\0'; i++) - { - if(input[i] == '\t') - { - break; - } - } - - input[i] = '\0'; - i++; - if(input[i] == 'p') (*pm) = FATHER; - if(input[i] == 'm') (*pm) = MOTHER; } diff --git a/Process_Read.h b/Process_Read.h index 23764e6..ec9abd1 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -5,7 +5,6 @@ #include #include #include -#include "kseq.h" #include "Overlaps.h" #include "CommandLines.h" ///#include "Hash_Table.h" @@ -18,16 +17,11 @@ #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_READ_LENGTH(R_INF, ID) R_INF.read_length[(ID)] -#define Get_NAME_LENGTH(R_INF, ID) (R_INF.name_index[(ID)+1] - R_INF.name_index[(ID)]) +#define Get_READ_LENGTH(R_INF, ID) (R_INF).read_length[(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]>>2) + ID -#define Get_READ(R_INF, ID) R_INF.read_sperate[(ID)] -#define Get_NAME(R_INF, ID) R_INF.name + R_INF.name_index[(ID)] - - - -KSEQ_INIT(gzFile, gzread) - +#define Get_READ(R_INF, ID) (R_INF).read_sperate[(ID)] +#define Get_NAME(R_INF, ID) ((R_INF).name + (R_INF).name_index[(ID)]) extern uint8_t seq_nt6_table[256]; @@ -37,12 +31,9 @@ extern char s_H[5]; extern char rc_Table[5]; - #define RC_CHAR(x) rc_Table[seq_nt6_table[(uint8_t)x]] void init_aux_table(); -int get_read(kseq_t *s, int adapterLen); - typedef struct { @@ -76,7 +67,6 @@ inline void init_PAF_alloc(PAF_alloc* list) list->list = (PAF*)malloc(sizeof(PAF)*list->size); } - inline void append_PAF_alloc(PAF_alloc* list, PAF* e) { if(list->length+1 > list->size) @@ -89,9 +79,6 @@ inline void append_PAF_alloc(PAF_alloc* list, PAF* e) list->length++; } - - - typedef struct { /**[0-1] bits are type:**/ @@ -104,7 +91,7 @@ typedef struct uint32_t lost_base_length; uint32_t lost_base_size; uint32_t new_length; -}Compressed_Cigar_record; +} Compressed_Cigar_record; #define AMBIGU 0 #define FATHER 1 @@ -112,13 +99,13 @@ typedef struct #define MIX_TRIO 3 #define NON_TRIO 4 #define DROP 5 + typedef struct { uint64_t** N_site; ///uint8_t* read; char* name; - uint8_t** read_sperate; uint64_t* read_length; uint64_t* read_size; @@ -129,7 +116,6 @@ typedef struct ///uint64_t* index; uint64_t index_size; - ///name start pos in char* name uint64_t* name_index; uint64_t name_index_size; @@ -142,32 +128,10 @@ typedef struct ma_hit_t_alloc* paf; ma_hit_t_alloc* reverse_paf; - ma_sub_t* coverage_cut; - } All_reads; extern All_reads R_INF; -void malloc_All_reads(All_reads* r); - -typedef struct -{ - kseq_t* read; - long long num; - -} R_buffer_block; - - -typedef struct -{ - R_buffer_block* sub_block; - long long block_inner_size; - long long size; - long long num; - int all_read_end; -} R_buffer; - - typedef struct { char* seq; @@ -176,23 +140,12 @@ typedef struct long long RID; } UC_Read; -typedef struct -{ - gzFile fp; - kseq_t *seq; - char** reads; - int num_reads; - int idx; -} gz_files; - -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); -void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ); +void malloc_All_reads(All_reads* r); +void ha_insert_read_len(All_reads *r, int read_len, int name_len); +void ha_compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ); void init_UC_Read(UC_Read* r); -void recover_UC_Read(UC_Read* r, All_reads* R_INF, uint64_t ID); +void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID); void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID); void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID); void destory_UC_Read(UC_Read* r); @@ -201,12 +154,4 @@ void write_All_reads(All_reads* r, char* read_file_name); int load_All_reads(All_reads* r, char* read_file_name); void destory_All_reads(All_reads* r); -void destory_R_buffer_block(R_buffer_block* curr_sub_block); -void destory_R_buffer(); -void clear_R_buffer(); - -void init_gz_files(hifiasm_opt_t* asm_opt); -void destory_gz_files(); - - #endif diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index 6a9e0cb..a1d8f71 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -675,7 +675,7 @@ long long* r_y_interval_beg, long long* r_y_interval_end) xUid, 0, xReads->n - 1, Hap_rate, position_index, reverse_sources, read_g, ruIndex, rev, &y_interval_beg, &y_interval_end); /********************y*********************/ - } + } else abort(); (*r_x_interval_beg) = x_interval_beg; (*r_x_interval_end) = x_interval_end; diff --git a/Trio.cpp b/Trio.cpp index d9d1657..7e9d974 100644 --- a/Trio.cpp +++ b/Trio.cpp @@ -6,68 +6,35 @@ #include "khashl.h" // hash table #include "kthread.h" #include "Process_Read.h" -#include "Trio.h" +#include "htab.h" #include "CommandLines.h" -#define CALLOC(ptr, len) ((ptr) = (__typeof__(ptr))calloc((len), sizeof(*(ptr)))) -#define MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr)))) -#define REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr)))) +#define YAK_MAX_KMER 31 +#define YAK_COUNTER_BITS 10 // yak uses 10, but hifiasm uses 12; we have to copy over some yak code here due to this +#define YAK_N_COUNTS (1<>YAK_COUNTER_BITS == (b)>>YAK_COUNTER_BITS) // lower 8 bits for counts; higher bits for k-mer #define yak_ch_hash(a) ((a)>>YAK_COUNTER_BITS) -KHASHL_SET_INIT(, yak_ht_t, yak_ht, uint64_t, yak_ch_hash, yak_ch_eq) +KHASHL_SET_INIT(static klib_unused, yak_ht_t, yak_ht, uint64_t, yak_ch_hash, yak_ch_eq) -///#define CHUNK_SIZE 200000 +typedef struct { + struct yak_ht_t *h; +} yak_ch1_t; -unsigned char seq_nt4_table[256] = { // translate ACGT to 0123 - 0, 1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 -}; +typedef struct { + int k, pre, n_hash, n_shift; + uint64_t tot; + yak_ch1_t *h; +} yak_ch_t; -static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function -{ - key = (~key + (key << 21)) & mask; // key = (key << 21) - key - 1; - key = key ^ key >> 24; - key = ((key + (key << 3)) + (key << 8)) & mask; // key * 265 - key = key ^ key >> 14; - key = ((key + (key << 2)) + (key << 4)) & mask; // key * 21 - key = key ^ key >> 28; - key = (key + (key << 31)) & mask; - return key; -} - -static inline uint64_t yak_hash64_64(uint64_t key) -{ - key = ~key + (key << 21); - key = key ^ key >> 24; - key = (key + (key << 3)) + (key << 8); - key = key ^ key >> 14; - key = (key + (key << 2)) + (key << 4); - key = key ^ key >> 28; - key = key + (key << 31); - return key; -} - -static inline uint64_t yak_hash_long(uint64_t x[4]) -{ - int j = x[1] < x[3]? 0 : 1; - return yak_hash64_64(x[j<<1|0]) + yak_hash64_64(x[j<<1|1]); -} - -int yak_ch_get(const yak_ch_t *h, uint64_t x) +static int yak_ch_get(const yak_ch_t *h, uint64_t x) { int mask = (1<pre) - 1; yak_ht_t *g = h->h[x&mask].h; @@ -76,21 +43,7 @@ int yak_ch_get(const yak_ch_t *h, uint64_t x) return k == kh_end(g)? -1 : kh_key(g, k)&YAK_MAX_COUNT; } -yak_bf_t *yak_bf_init(int n_shift, int n_hashes) -{ - yak_bf_t *b; - void *ptr = 0; - if (n_shift + YAK_BLK_SHIFT > 64 || n_shift < YAK_BLK_SHIFT) return 0; - b = (yak_bf_t*)calloc(1, sizeof(yak_bf_t)); - b->n_shift = n_shift; - b->n_hashes = n_hashes; - posix_memalign(&ptr, 1<<(YAK_BLK_SHIFT-3), 1ULL<<(n_shift-3)); - b->b = (uint8_t*)ptr; - bzero(b->b, 1ULL<<(n_shift-3)); - return b; -} - -yak_ch_t *yak_ch_init(int k, int pre, int n_hash, int n_shift) +static yak_ch_t *yak_ch_init(int k, int pre) { yak_ch_t *h; int i; @@ -100,32 +53,10 @@ yak_ch_t *yak_ch_init(int k, int pre, int n_hash, int n_shift) CALLOC(h->h, 1<pre); for (i = 0; i < 1<pre; ++i) h->h[i].h = yak_ht_init(); - if (n_hash > 0 && n_shift > h->pre) { - h->n_hash = n_hash, h->n_shift = n_shift; - for (i = 0; i < 1<pre; ++i) - h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash); - } return h; } -void yak_bf_destroy(yak_bf_t *b) -{ - if (b == 0) return; - free(b->b); free(b); -} - -void yak_ch_destroy_bf(yak_ch_t *h) -{ - int i; - for (i = 0; i < 1<pre; ++i) { - if (h->h[i].b) - yak_bf_destroy(h->h[i].b); - h->h[i].b = 0; - } -} - - -yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ...) +static yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ...) { va_list ap; FILE *fp; @@ -161,13 +92,13 @@ yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ...) return 0; } - ch = ch0 == 0? yak_ch_init(t[0], t[1], 0, 0) : ch0; + ch = ch0 == 0? yak_ch_init(t[0], t[1]) : ch0; assert((int)t[0] == ch->k && (int)t[1] == ch->pre); for (i = 0; i < 1<pre; ++i) { yak_ht_t *h = ch->h[i].h; fread(t, 4, 2, fp); if (ch0 == 0) yak_ht_resize(h, t[0]); - for (j = 0; j < t[1]; ++j) { + for (j = 0; j < (int)t[1]; ++j) { uint64_t key; fread(&key, 8, 1, fp); if (mode == YAK_LOAD_ALL) { @@ -195,17 +126,39 @@ yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ...) return ch; } - -void yak_ch_destroy(yak_ch_t *h) +static void yak_ch_destroy(yak_ch_t *h) { int i; if (h == 0) return; - yak_ch_destroy_bf(h); for (i = 0; i < 1<pre; ++i) yak_ht_destroy(h->h[i].h); free(h->h); free(h); } +typedef struct { + int max; + uint32_t *s; +} tb_buf_t; + +typedef struct { + int k, n_threads, print_diff; + double ratio_thres; + const yak_ch_t *ch; + tb_buf_t *buf; + UC_Read *bseq; + All_reads* seq; +} tb_shared_t; + +typedef struct { + int c[16]; + int sc[2]; + int nk; +} tb_cnt_t; + +typedef struct { + int n_seq; + tb_shared_t *aux; +} tb_step_t; static char tb_classify(const int sc[2], const int *c, int k, double ratio_thres) { @@ -229,8 +182,6 @@ static char tb_classify(const int sc[2], const int *c, int k, double ratio_thres static void tb_worker(void *_data, long k, int tid) { - ///tb_step_t *t = (tb_step_t*)_data; - ///tb_shared_t *aux = t->aux; tb_shared_t *aux = (tb_shared_t*)_data; UC_Read *s = &aux->bseq[tid]; recover_UC_Read(s, aux->seq, k); @@ -266,12 +217,7 @@ static void tb_worker(void *_data, long k, int tid) if (++l >= aux->k) { int type = 0, c1, c2; uint64_t y; - - - //++t->cnt[k].nk; ++cnt.nk; - - if (aux->ch->k < 32) y = yak_hash64(x[0] < x[1]? x[0] : x[1], mask); else @@ -282,28 +228,18 @@ static void tb_worker(void *_data, long k, int tid) if (c1 == 2 && c2 == 0) type = 1; else if (c2 == 2 && c1 == 0) type = 2; b->s[i] = type; - - - ///++t->cnt[k].c[flag]; ++cnt.c[flag]; - - // if (aux->print_diff && (flag>>2&3) != (flag&3)) - // printf("D\t%s\t%d\t%d\t%d\n", s->name, i, flag&3, flag>>2&3); } } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; } for (l = 0, i = 1; i <= s->length; ++i) { if (i == s->length || b->s[i] != b->s[l]) { if (b->s[l] > 0 && i - l >= aux->k - 4) - { - ///t->cnt[k].sc[b->s[l] - 1] += i - l; cnt.sc[b->s[l] - 1] += i - l; - } l = i; } } - int *c = cnt.c; char type; type = tb_classify(cnt.sc, c, aux->k, aux->ratio_thres); @@ -312,12 +248,9 @@ static void tb_worker(void *_data, long k, int tid) if(type == 'm') aux->seq->trio_flag[k] = MOTHER; } - - -void trio_partition() +void trio_partition(void) { - if(asm_opt.pat_index == NULL || asm_opt.mat_index == NULL) - { + if (asm_opt.pat_index == NULL || asm_opt.mat_index == NULL) { memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t)); return; } @@ -326,32 +259,26 @@ void trio_partition() fprintf(stderr, "Start trio binning ...... \n"); yak_ch_t *ch; - int i/**, min_cnt = 2, mid_cnt = 5**/; + int i /**, min_cnt = 2, mid_cnt = 5**/; tb_shared_t aux; memset(&aux, 0, sizeof(tb_shared_t)); aux.n_threads = asm_opt.thread_num, aux.print_diff = 0; aux.ratio_thres = 0.33; aux.seq = &R_INF; - ch = yak_ch_restore_core(0, asm_opt.pat_index, YAK_LOAD_TRIOBIN1, asm_opt.min_cnt, asm_opt.mid_cnt); ch = yak_ch_restore_core(ch, asm_opt.mat_index, YAK_LOAD_TRIOBIN2, asm_opt.min_cnt, asm_opt.mid_cnt); - - aux.k = ch->k; aux.ch = ch; aux.buf = (tb_buf_t*)calloc(aux.n_threads, sizeof(tb_buf_t)); aux.bseq = (UC_Read*)calloc(aux.n_threads, sizeof(UC_Read)); for (i = 0; i < aux.n_threads; ++i) - { init_UC_Read(&aux.bseq[i]); - } kt_for(aux.n_threads, tb_worker, &aux, aux.seq->total_reads); - for (i = 0; i < aux.n_threads; ++i) - { + for (i = 0; i < aux.n_threads; ++i) { free(aux.buf[i].s); destory_UC_Read(&aux.bseq[i]); } @@ -361,4 +288,4 @@ void trio_partition() fprintf(stderr, "Trio binning has been done.\n"); fprintf(stderr, "%-30s%18.2f\n\n", "Trio binning time:", Get_T() - start_time); -} \ No newline at end of file +} diff --git a/Trio.h b/Trio.h deleted file mode 100644 index d5319f5..0000000 --- a/Trio.h +++ /dev/null @@ -1,71 +0,0 @@ -#ifndef __TRIO__ -#define __TRIO__ -#include - -#define YAK_MAX_KMER 31 -#define YAK_COUNTER_BITS 10 -#define YAK_N_COUNTS (1< +#include "htab.h" +#include "ksort.h" +#include "Hash_Table.h" + +typedef struct { // this struct is not strictly necessary; we can use k_mer_pos instead, with modifications + uint64_t srt; + uint32_t self_off; + uint32_t other_off; +} anchor1_t; + +#define an_key1(a) ((a).srt) +#define an_key2(a) ((a).self_off) +KRADIX_SORT_INIT(ha_an1, anchor1_t, an_key1, 8) +KRADIX_SORT_INIT(ha_an2, anchor1_t, an_key2, 4) + +#define oreg_xs_lt(a, b) (((uint64_t)(a).x_pos_s<<32|(a).x_pos_e) < ((uint64_t)(b).x_pos_s<<32|(b).x_pos_e)) +KSORT_INIT(or_xs, overlap_region, oreg_xs_lt) + +#define oreg_ss_lt(a, b) ((a).shared_seed > (b).shared_seed) // in the decending order +KSORT_INIT(or_ss, overlap_region, oreg_ss_lt) + +typedef struct { + int n; + const ha_idxpos_t *a; +} seed1_t; + +struct ha_abuf_s { + uint64_t n_a, m_a; + uint32_t old_mz_m; + ha_mz1_v mz; + seed1_t *seed; + anchor1_t *a; +}; + +ha_abuf_t *ha_abuf_init(void) +{ + return (ha_abuf_t*)calloc(1, sizeof(ha_abuf_t)); +} + +void ha_abuf_destroy(ha_abuf_t *ab) +{ + free(ab->seed); free(ab->a); free(ab->mz.a); free(ab); +} + +uint64_t ha_abuf_mem(const ha_abuf_t *ab) +{ + return ab->m_a * sizeof(anchor1_t) + ab->mz.m * (sizeof(ha_mz1_t) + sizeof(seed1_t)) + sizeof(ha_abuf_t); +} + +void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain) +{ + extern void *ha_flt_tab; + extern ha_pt_t *ha_idx; + uint32_t i; + uint64_t k, l; + + // prepare + clear_Candidates_list(cl); + clear_overlap_region_alloc(overlap_list); + recover_UC_Read(ucr, &R_INF, rid); + ab->mz.n = 0, ab->n_a = 0; + + // get the list of anchors + ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !asm_opt.no_HPC, &ab->mz, ha_flt_tab); + if (ab->mz.m > ab->old_mz_m) { + ab->old_mz_m = ab->mz.m; + REALLOC(ab->seed, ab->old_mz_m); + } + for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) { + ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &ab->seed[i].n); + ab->n_a += ab->seed[i].n; + } + if (ab->n_a > ab->m_a) { + ab->m_a = ab->n_a; + kroundup64(ab->m_a); + REALLOC(ab->a, ab->m_a); + } + for (i = 0, k = 0; i < ab->mz.n; ++i) { + int j; + ha_mz1_t *z = &ab->mz.a[i]; + seed1_t *s = &ab->seed[i]; + for (j = 0; j < s->n; ++j) { + const ha_idxpos_t *y = &s->a[j]; + anchor1_t *an = &ab->a[k++]; + uint8_t rev = z->rev == y->rev? 0 : 1; + an->other_off = y->pos; + an->self_off = rev? ucr->length - 1 - (z->pos + 1 - z->span) : z->pos; + an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off; + } + } + + // sort anchors + radix_sort_ha_an1(ab->a, ab->a + ab->n_a); + for (k = 1, l = 0; k <= ab->n_a; ++k) { + if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) { + if (k - l > 1) + radix_sort_ha_an2(ab->a + l, ab->a + k); + l = k; + } + } + + // copy over to _cl_ + if (ab->m_a >= (uint64_t)cl->size) { + cl->size = ab->m_a; + REALLOC(cl->list, cl->size); + } + for (k = 0; k < ab->n_a; ++k) { + k_mer_hit *p = &cl->list[k]; + p->readID = ab->a[k].srt >> 33; + p->strand = ab->a[k].srt >> 32 & 1; + p->offset = ab->a[k].other_off; + p->self_offset = ab->a[k].self_off; + } + cl->length = ab->n_a; + + calculate_overlap_region_by_chaining(cl, overlap_list, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain); + + #if 0 + if (overlap_list->length > 2000) { + fprintf(stderr, "B\t%ld\t%ld\t%ld\n", (long)rid, (long)overlap_list->length, (long)Get_READ_LENGTH(R_INF, rid)); + for (int i = 0; i < (int)overlap_list->length; ++i) { + overlap_region *r = &overlap_list->list[i]; + fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand], + (int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed); + } + } + #endif + + if ((int)overlap_list->length > max_n_chain) { + uint32_t n[2], s[2]; + n[0] = n[1] = 0, s[0] = s[1] = 0; + for (i = 0; i < (uint32_t)overlap_list->length; ++i) { + const overlap_region *r = &overlap_list->list[i]; + int dir = r->x_pos_s == 0? 0 : 1; + ++n[dir]; + if ((int)n[dir] == max_n_chain) s[dir] = r->shared_seed; + } + if (s[0] > 0 || s[1] > 0) { + for (i = 0, k = 0; i < (uint32_t)overlap_list->length; ++i) { + overlap_region *r = &overlap_list->list[i]; + int dir = r->x_pos_s == 0? 0 : 1; + if (r->shared_seed > s[dir]) { + if ((uint32_t)k != i) { + overlap_region t; + t = overlap_list->list[k]; + overlap_list->list[k] = overlap_list->list[i]; + overlap_list->list[i] = t; + } + ++k; + } + } + overlap_list->length = k; + } + } + + ks_introsort_or_xs(overlap_list->length, overlap_list->list); +} diff --git a/hist.cpp b/hist.cpp new file mode 100644 index 0000000..9a1e31e --- /dev/null +++ b/hist.cpp @@ -0,0 +1,92 @@ +#include +#include "htab.h" + +static void ha_hist_line(int c, int x, int exceed, int64_t cnt) +{ + int j; + if (c >= 0) fprintf(stderr, "[M::%s] %5d: ", __func__, c); + else fprintf(stderr, "[M::%s] %5s: ", __func__, "rest"); + for (j = 0; j < x; ++j) fputc('*', stderr); + if (exceed) fputc('>', stderr); + fprintf(stderr, " %lld\n", (long long)cnt); +} + +int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het) +{ + const int hist_max = 100; + int i, start, low_i, max_i, max2_i, max3_i; + int64_t max, max2, max3, min; + + // find the low point from the left + *peak_het = -1; + start = cnt[1] > 0? 1 : 2; + low_i = start; + for (i = low_i + 1; i < n_cnt; ++i) + if (cnt[i] > cnt[i-1]) break; + low_i = i - 1; + fprintf(stderr, "[M::%s] lowest: count[%d] = %ld\n", __func__, low_i, (long)cnt[low_i]); + if (low_i == n_cnt - 1) return -1; // low coverage + + // find the highest peak + max_i = low_i + 1, max = cnt[max_i]; + for (i = low_i + 1; i < n_cnt; ++i) + if (cnt[i] > max) + max = cnt[i], max_i = i; + fprintf(stderr, "[M::%s] highest: count[%d] = %ld\n", __func__, max_i, (long)cnt[max_i]); + + // print histogram + for (i = start; i < n_cnt; ++i) { + int x, exceed = 0; + x = (int)((double)hist_max * cnt[i] / cnt[max_i] + .499); + if (x > hist_max) exceed = 1, x = hist_max; // may happen if cnt[2] is higher + if (i > max_i && x == 0) break; + ha_hist_line(i, x, exceed, cnt[i]); + } + { + int x, exceed = 0; + int64_t rest = 0; + for (; i < n_cnt; ++i) rest += cnt[i]; + x = (int)((double)hist_max * rest / cnt[max_i] + .499); + if (x > hist_max) exceed = 1, x = hist_max; + ha_hist_line(-1, x, exceed, rest); + } + + // look for smaller peak on the low end + max2 = -1; max2_i = -1; + for (i = max_i - 1; i > low_i; --i) { + if (cnt[i] >= cnt[i-1] && cnt[i] >= cnt[i+1]) { + if (cnt[i] > max2) max2 = cnt[i], max2_i = i; + } + } + if (max2_i > low_i && max2_i < max_i) { + for (i = max2_i + 1, min = max; i < max_i; ++i) + if (cnt[i] < min) min = cnt[i]; + if (max2 < max * 0.05 || min > max2 * 0.95) + max2 = -1, max2_i = -1; + } + if (max2 > 0) fprintf(stderr, "[M::%s] left: count[%d] = %ld\n", __func__, max2_i, (long)cnt[max2_i]); + else fprintf(stderr, "[M::%s] left: none\n", __func__); + + // look for smaller peak on the high end + max3 = -1; max3_i = -1; + for (i = max_i + 1; i < n_cnt - 1; ++i) { + if (cnt[i] >= cnt[i-1] && cnt[i] >= cnt[i+1]) { + if (cnt[i] > max3) max3 = cnt[i], max3_i = i; + } + } + if (max3_i > max_i) { + for (i = max_i + 1, min = max; i < max3_i; ++i) + if (cnt[i] < min) min = cnt[i]; + if (max3 < max * 0.05 || min > max3 * 0.95 || max3_i > max_i * 2.5) + max3 = -1, max3_i = -1; + } + if (max3 > 0) fprintf(stderr, "[M::%s] right: count[%d] = %ld\n", __func__, max3_i, (long)cnt[max3_i]); + else fprintf(stderr, "[M::%s] right: none\n", __func__); + if (max3_i > 0) { + *peak_het = max_i; + return max3_i; + } else { + if (max2_i > 0) *peak_het = max2_i; + return max_i; + } +} diff --git a/htab.cpp b/htab.cpp new file mode 100644 index 0000000..d96b21d --- /dev/null +++ b/htab.cpp @@ -0,0 +1,832 @@ +#include +#include +#include +#include +#include +#include +#include "kthread.h" +#include "khashl.h" +#include "kseq.h" +#include "ksort.h" +#include "htab.h" + +#define YAK_COUNTER_BITS 12 +#define YAK_N_COUNTS (1<bf_shift = 0; + o->bf_n_hash = 4; + o->k = 31; + o->w = 1; + o->pre = YAK_COUNTER_BITS; + o->n_thread = 4; + o->chunk_size = 20000000; +} + +/************************ + * Blocked bloom filter * + ************************/ + +#define YAK_BLK_SHIFT 9 // 64 bytes, the size of a cache line +#define YAK_BLK_MASK ((1<<(YAK_BLK_SHIFT)) - 1) + +typedef struct { + int n_shift, n_hashes; + uint8_t *b; +} yak_bf_t; + +yak_bf_t *yak_bf_init(int n_shift, int n_hashes) +{ + yak_bf_t *b; + void *ptr = 0; + if (n_shift + YAK_BLK_SHIFT > 64 || n_shift < YAK_BLK_SHIFT) return 0; + CALLOC(b, 1); + b->n_shift = n_shift; + b->n_hashes = n_hashes; + posix_memalign(&ptr, 1<<(YAK_BLK_SHIFT-3), 1ULL<<(n_shift-3)); + b->b = (uint8_t*)ptr; + bzero(b->b, 1ULL<<(n_shift-3)); + return b; +} + +void yak_bf_destroy(yak_bf_t *b) +{ + if (b == 0) return; + free(b->b); free(b); +} + +int yak_bf_insert(yak_bf_t *b, uint64_t hash) +{ + int x = b->n_shift - YAK_BLK_SHIFT; + uint64_t y = hash & ((1ULL<> x & YAK_BLK_MASK; + int h2 = hash >> b->n_shift & YAK_BLK_MASK; + uint8_t *p = &b->b[y<<(YAK_BLK_SHIFT-3)]; + int i, z = h1, cnt = 0; + if ((h2&31) == 0) h2 = (h2 + 1) & YAK_BLK_MASK; // otherwise we may repeatedly use a few bits + for (i = 0; i < b->n_hashes; z = (z + h2) & YAK_BLK_MASK) { + uint8_t *q = &p[z>>3], u; + u = 1<<(z&7); + cnt += !!(*q & u); + *q |= u; + ++i; + } + return cnt; +} + +/******************** + * Count hash table * + ********************/ + +#define yak_ct_eq(a, b) ((a)>>YAK_COUNTER_BITS == (b)>>YAK_COUNTER_BITS) // lower 8 bits for counts; higher bits for k-mer +#define yak_ct_hash(a) ((a)>>YAK_COUNTER_BITS) +KHASHL_SET_INIT(static klib_unused, yak_ct_t, yak_ct, uint64_t, yak_ct_hash, yak_ct_eq) + +typedef struct { + yak_ct_t *h; + yak_bf_t *b; +} ha_ct1_t; + +typedef struct { + int k, pre, n_hash, n_shift; + uint64_t tot; + ha_ct1_t *h; +} ha_ct_t; + +static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift) +{ + ha_ct_t *h; + int i; + if (pre < YAK_COUNTER_BITS) return 0; + CALLOC(h, 1); + h->k = k, h->pre = pre; + CALLOC(h->h, 1<pre); + for (i = 0; i < 1<pre; ++i) + h->h[i].h = yak_ct_init(); + if (n_hash > 0 && n_shift > h->pre) { + h->n_hash = n_hash, h->n_shift = n_shift; + for (i = 0; i < 1<pre; ++i) + h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash); + } + return h; +} + +static void ha_ct_destroy_bf(ha_ct_t *h) +{ + int i; + for (i = 0; i < 1<pre; ++i) { + if (h->h[i].b) + yak_bf_destroy(h->h[i].b); + h->h[i].b = 0; + } +} + +static void ha_ct_destroy(ha_ct_t *h) +{ + int i; + if (h == 0) return; + ha_ct_destroy_bf(h); + for (i = 0; i < 1<pre; ++i) + yak_ct_destroy(h->h[i].h); + free(h->h); free(h); +} + +static int ha_ct_insert_list(ha_ct_t *h, int create_new, int n, const uint64_t *a) +{ + int j, mask = (1<pre) - 1, n_ins = 0; + ha_ct1_t *g; + if (n == 0) return 0; + g = &h->h[a[0]&mask]; + for (j = 0; j < n; ++j) { + int ins = 1, absent; + uint64_t x = a[j] >> h->pre; + khint_t k; + if ((a[j]&mask) != (a[0]&mask)) continue; + if (create_new) { + if (g->b) + ins = (yak_bf_insert(g->b, x) == h->n_hash); + if (ins) { + k = yak_ct_put(g->h, x << YAK_COUNTER_BITS | (g->b? 1 : 0), &absent); + if (absent) ++n_ins; + if ((kh_key(g->h, k)&YAK_MAX_COUNT) < YAK_MAX_COUNT) + ++kh_key(g->h, k); + } + } else { + k = yak_ct_get(g->h, x<h) && (kh_key(g->h, k)&YAK_MAX_COUNT) < YAK_MAX_COUNT) + ++kh_key(g->h, k); + } + } + return n_ins; +} + +/*** generate histogram ***/ + +typedef struct { + uint64_t c[YAK_N_COUNTS]; +} buf_cnt_t; + +typedef struct { + const ha_ct_t *h; + buf_cnt_t *cnt; +} hist_aux_t; + +static void worker_ct_hist(void *data, long i, int tid) // callback for kt_for() +{ + hist_aux_t *a = (hist_aux_t*)data; + uint64_t *cnt = a->cnt[tid].c; + yak_ct_t *g = a->h->h[i].h; + khint_t k; + for (k = 0; k < kh_end(g); ++k) + if (kh_exist(g, k)) + ++cnt[kh_key(g, k)&YAK_MAX_COUNT]; +} + +static void ha_ct_hist(const ha_ct_t *h, int64_t cnt[YAK_N_COUNTS], int n_thread) +{ + hist_aux_t a; + int i, j; + a.h = h; + memset(cnt, 0, YAK_N_COUNTS * sizeof(uint64_t)); + CALLOC(a.cnt, n_thread); + kt_for(n_thread, worker_ct_hist, &a, 1<pre); + for (i = 0; i < YAK_N_COUNTS; ++i) cnt[i] = 0; + for (j = 0; j < n_thread; ++j) + for (i = 0; i < YAK_N_COUNTS; ++i) + cnt[i] += a.cnt[j].c[i]; + free(a.cnt); +} + +/*** shrink a hash table ***/ + +typedef struct { + int min, max; + ha_ct_t *h; +} shrink_aux_t; + +static void worker_ct_shrink(void *data, long i, int tid) // callback for kt_for() +{ + shrink_aux_t *a = (shrink_aux_t*)data; + ha_ct_t *h = a->h; + yak_ct_t *g = h->h[i].h, *f; + khint_t k; + f = yak_ct_init(); + yak_ct_resize(f, kh_size(g)); + for (k = 0; k < kh_end(g); ++k) { + if (kh_exist(g, k)) { + int absent, c = kh_key(g, k) & YAK_MAX_COUNT; + if (c >= a->min && c <= a->max) + yak_ct_put(f, kh_key(g, k), &absent); + } + } + yak_ct_destroy(g); + h->h[i].h = f; +} + +static void ha_ct_shrink(ha_ct_t *h, int min, int max, int n_thread) +{ + int i; + shrink_aux_t a; + a.h = h, a.min = min, a.max = max; + kt_for(n_thread, worker_ct_shrink, &a, 1<pre); + for (i = 0, h->tot = 0; i < 1<pre; ++i) + h->tot += kh_size(h->h[i].h); +} + +/*********************** + * Position hash table * + ***********************/ + +KHASHL_MAP_INIT(static klib_unused, yak_pt_t, yak_pt, uint64_t, uint64_t, yak_ct_hash, yak_ct_eq) +#define generic_key(x) (x) +KRADIX_SORT_INIT(ha64, uint64_t, generic_key, 8) + +typedef struct { + yak_pt_t *h; + uint64_t n; + ha_idxpos_t *a; +} ha_pt1_t; + +struct ha_pt_s { + int k, pre; + uint64_t tot, tot_pos; + ha_pt1_t *h; +}; + +typedef struct { + const ha_ct_t *ct; + ha_pt_t *pt; +} pt_gen_aux_t; + +static void worker_pt_gen(void *data, long i, int tid) // callback for kt_for() +{ + pt_gen_aux_t *a = (pt_gen_aux_t*)data; + ha_pt1_t *b = &a->pt->h[i]; + yak_ct_t *g = a->ct->h[i].h; + khint_t k; + for (k = 0, b->n = 0; k != kh_end(g); ++k) { + if (kh_exist(g, k)) { + int absent; + khint_t l; + l = yak_pt_put(b->h, kh_key(g, k) >> a->ct->pre << YAK_COUNTER_BITS, &absent); + kh_val(b->h, l) = b->n; + b->n += kh_key(g, k) & YAK_MAX_COUNT; + } + } + yak_ct_destroy(g); + a->ct->h[i].h = 0; + CALLOC(b->a, b->n); +} + +ha_pt_t *ha_pt_gen(ha_ct_t *ct, int n_thread) +{ + pt_gen_aux_t a; + int i; + ha_pt_t *pt; + ha_ct_destroy_bf(ct); + CALLOC(pt, 1); + pt->k = ct->k, pt->pre = ct->pre, pt->tot = ct->tot; + CALLOC(pt->h, 1<pre); + for (i = 0; i < 1<pre; ++i) { + pt->h[i].h = yak_pt_init(); + yak_pt_resize(pt->h[i].h, kh_size(ct->h[i].h)); + } + a.ct = ct, a.pt = pt; + kt_for(n_thread, worker_pt_gen, &a, 1<pre); + free(ct->h); free(ct); + return pt; +} + +int ha_pt_insert_list(ha_pt_t *h, int n, const ha_mz1_t *a) +{ + int j, mask = (1<pre) - 1, n_ins = 0; + ha_pt1_t *g; + if (n == 0) return 0; + g = &h->h[a[0].x&mask]; + for (j = 0; j < n; ++j) { + uint64_t x = a[j].x >> h->pre; + khint_t k; + int n; + ha_idxpos_t *p; + if ((a[j].x&mask) != (a[0].x&mask)) continue; + k = yak_pt_get(g->h, x<h)) continue; + n = kh_key(g->h, k) & YAK_MAX_COUNT; + assert(n < YAK_MAX_COUNT); + p = &g->a[kh_val(g->h, k) + n]; + p->rid = a[j].rid, p->rev = a[j].rev, p->pos = a[j].pos, p->span = a[j].span; + //(uint64_t)a[j].rid<<36 | (uint64_t)a[j].rev<<35 | (uint64_t)a[j].pos<<8 | (uint64_t)a[j].span; + ++kh_key(g->h, k); + ++n_ins; + } + return n_ins; +} +/* +static void worker_pt_sort(void *data, long i, int tid) +{ + ha_pt_t *h = (ha_pt_t*)data; + ha_pt1_t *g = &h->h[i]; + khint_t k; + for (k = 0; k < kh_end(g->h); ++k) { + int n; + uint64_t *p; + if (!kh_exist(g->h, k)) continue; + n = kh_key(g->h, k) & YAK_MAX_COUNT; + p = &g->a[kh_val(g->h, k)]; + radix_sort_ha64(p, p + n); + } +} + +void ha_pt_sort(ha_pt_t *h, int n_thread) +{ + kt_for(n_thread, worker_pt_sort, h, 1<pre); +} +*/ +void ha_pt_destroy(ha_pt_t *h) +{ + int i; + if (h == 0) return; + for (i = 0; i < 1<pre; ++i) { + yak_pt_destroy(h->h[i].h); + free(h->h[i].a); + } + free(h->h); free(h); +} + +const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n) +{ + khint_t k; + const ha_pt1_t *g = &h->h[hash & ((1ULL<pre) - 1)]; + *n = 0; + k = yak_pt_get(g->h, hash >> h->pre << YAK_COUNTER_BITS); + if (k == kh_end(g->h)) return 0; + *n = kh_key(g->h, k) & YAK_MAX_COUNT; + return &g->a[kh_val(g->h, k)]; +} + +/********************************** + * Buffer for counting all k-mers * + **********************************/ + +typedef struct { + int n, m; + uint64_t n_ins; + uint64_t *a; + ha_mz1_t *b; +} ch_buf_t; + +static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k-mer $y to a linear buffer +{ + int pre = y & ((1<n == b->m) { + b->m = b->m < 8? 8 : b->m + (b->m>>1); + REALLOC(b->a, b->m); + } + b->a[b->n++] = y; +} + +static inline void pt_insert_buf(ch_buf_t *buf, int p, const ha_mz1_t *y) +{ + int pre = y->x & ((1<n == b->m) { + b->m = b->m < 8? 8 : b->m + (b->m>>1); + REALLOC(b->b, b->m); + } + b->b[b->n++] = *y; +} + +static void count_seq_buf(ch_buf_t *buf, int k, int p, int len, const char *seq) // insert k-mers in $seq to linear buffer $buf +{ + int i, l; + uint64_t x[4], mask = (1ULL<>1)) & mask; + x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift; + x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; + if (++l >= k) + ct_insert_buf(buf, p, yak_hash_long(x)); + } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart + } +} + +static void count_seq_buf_HPC(ch_buf_t *buf, int k, int p, int len, const char *seq) // insert k-mers in $seq to linear buffer $buf +{ + int i, l, last = -1; + uint64_t x[4], mask = (1ULL<>1)) & mask; + x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift; + x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; + if (++l >= k) + ct_insert_buf(buf, p, yak_hash_long(x)); + last = c; + } + } else l = 0, last = -1, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart + } +} + +/****************** + * K-mer counting * + ******************/ + +KSEQ_INIT(gzFile, gzread) + +#define HAF_COUNT_EXACT 0x1 +#define HAF_COUNT_ALL 0x2 +#define HAF_RS_WRITE_LEN 0x4 +#define HAF_RS_WRITE_SEQ 0x8 +#define HAF_RS_READ 0x10 + +typedef struct { // global data structure for kt_pipeline() + const yak_copt_t *opt; + const void *flt_tab; + int flag, create_new, is_store; + uint64_t n_seq; + kseq_t *ks; + UC_Read ucr; + ha_ct_t *ct; + ha_pt_t *pt; + const All_reads *rs_in; + All_reads *rs_out; +} pl_data_t; + +typedef struct { // data structure for each step in kt_pipeline() + pl_data_t *p; + uint64_t n_seq0; + int n_seq, m_seq, sum_len, nk; + int *len; + char **seq; + ha_mz1_v *mz_buf; + ha_mz1_v *mz; + ch_buf_t *buf; +} st_data_t; + +static void worker_for_insert(void *data, long i, int tid) // callback for kt_for() +{ + st_data_t *s = (st_data_t*)data; + ch_buf_t *b = &s->buf[i]; + if (s->p->pt) + b->n_ins += ha_pt_insert_list(s->p->pt, b->n, b->b); + else + b->n_ins += ha_ct_insert_list(s->p->ct, s->p->create_new, b->n, b->a); +} + +static void worker_for_mz(void *data, long i, int tid) +{ + st_data_t *s = (st_data_t*)data; + ha_mz1_v *b = &s->mz_buf[tid]; + s->mz_buf[tid].n = 0; + ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab); + s->mz[i].n = s->mz[i].m = b->n; + MALLOC(s->mz[i].a, b->n); + memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t)); +} + +static void *worker_count(void *data, int step, void *in) // callback for kt_pipeline() +{ + pl_data_t *p = (pl_data_t*)data; + if (step == 0) { // step 1: read a block of sequences + int ret; + st_data_t *s; + CALLOC(s, 1); + s->p = p; + s->n_seq0 = p->n_seq; + if (p->rs_in && (p->flag & HAF_RS_READ)) { + while (p->n_seq < p->rs_in->total_reads) { + int l; + recover_UC_Read(&p->ucr, p->rs_in, p->n_seq); + l = p->ucr.length; + if (s->n_seq == s->m_seq) { + s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1); + REALLOC(s->len, s->m_seq); + REALLOC(s->seq, s->m_seq); + } + MALLOC(s->seq[s->n_seq], l); + memcpy(s->seq[s->n_seq], p->ucr.seq, l); + s->len[s->n_seq++] = l; + ++p->n_seq; + s->sum_len += l; + s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0; + if (s->sum_len >= p->opt->chunk_size) + break; + } + } else { + while ((ret = kseq_read(p->ks)) >= 0) { + int l = p->ks->seq.l; + if (p->n_seq >= 1<<28) { + fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28); + exit(1); + } + if (p->rs_out) { + if (p->flag & HAF_RS_WRITE_LEN) { + assert(p->n_seq == p->rs_out->total_reads); + ha_insert_read_len(p->rs_out, l, p->ks->name.l); + } else if (p->flag & HAF_RS_WRITE_SEQ) { + int i, n_N; + assert(l == (int)p->rs_out->read_length[p->n_seq]); + for (i = n_N = 0; i < l; ++i) // count number of ambiguous bases + if (seq_nt4_table[(uint8_t)p->ks->seq.s[i]] >= 4) + ++n_N; + ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s, l, &p->rs_out->N_site[p->n_seq], n_N); + memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l); + } + } + if (s->n_seq == s->m_seq) { + s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1); + REALLOC(s->len, s->m_seq); + REALLOC(s->seq, s->m_seq); + } + MALLOC(s->seq[s->n_seq], l); + memcpy(s->seq[s->n_seq], p->ks->seq.s, l); + s->len[s->n_seq++] = l; + ++p->n_seq; + s->sum_len += l; + s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0; + if (s->sum_len >= p->opt->chunk_size) + break; + } + } + if (s->sum_len == 0) free(s); + else return s; + } else if (step == 1) { // step 2: extract k-mers + st_data_t *s = (st_data_t*)in; + int i, n_pre = 1<opt->pre, m; + // allocate the k-mer buffer + CALLOC(s->buf, n_pre); + m = (int)(s->nk * 1.2 / n_pre) + 1; + for (i = 0; i < n_pre; ++i) { + s->buf[i].m = m; + if (p->pt) MALLOC(s->buf[i].b, m); + else MALLOC(s->buf[i].a, m); + } + // fill the buffer + if (p->opt->w == 1) { // enumerate all k-mers + for (i = 0; i < s->n_seq; ++i) { + if (p->opt->is_HPC) + count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]); + else + count_seq_buf(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]); + if (!p->is_store) free(s->seq[i]); + } + } else { // minimizers only + uint32_t j; + // compute minimizers + CALLOC(s->mz, s->n_seq); + CALLOC(s->mz_buf, p->opt->n_thread); + kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq); + for (i = 0; i < p->opt->n_thread; ++i) + free(s->mz_buf[i].a); + free(s->mz_buf); + // insert minimizers + if (p->pt) { + for (i = 0; i < s->n_seq; ++i) + for (j = 0; j < s->mz[i].n; ++j) + pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]); + } else { + for (i = 0; i < s->n_seq; ++i) + for (j = 0; j < s->mz[i].n; ++j) + ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x); + } + for (i = 0; i < s->n_seq; ++i) { + free(s->mz[i].a); + if (!p->is_store) free(s->seq[i]); + } + free(s->mz); + } + free(s->seq); free(s->len); + s->seq = 0, s->len = 0; + return s; + } else if (step == 2) { // step 3: insert k-mers to hash table + st_data_t *s = (st_data_t*)in; + int i, n = 1<opt->pre; + uint64_t n_ins = 0; + kt_for(p->opt->n_thread, worker_for_insert, s, n); + for (i = 0; i < n; ++i) { + n_ins += s->buf[i].n_ins; + if (p->pt) free(s->buf[i].b); + else free(s->buf[i].a); + } + if (p->ct) p->ct->tot += n_ins; + if (p->pt) p->pt->tot_pos += n_ins; + free(s->buf); + #if 0 + fprintf(stderr, "[M::%s::%.3f*%.2f] processed %ld sequences; %ld %s in the hash table\n", __func__, + yak_realtime(), yak_cpu_usage(), (long)s->n_seq0 + s->n_seq, + (long)(p->pt? p->pt->tot_pos : p->ct->tot), p->pt? "positions" : "distinct k-mers"); + #endif + free(s); + } + return 0; +} + +static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt_t *p0, ha_ct_t *c0, const void *flt_tab, All_reads *rs) +{ + int read_rs = (rs && (flag & HAF_RS_READ)); + pl_data_t pl; + gzFile fp = 0; + memset(&pl, 0, sizeof(pl_data_t)); + if (read_rs) { + pl.rs_in = rs; + init_UC_Read(&pl.ucr); + } else { + if ((fp = gzopen(fn, "r")) == 0) return 0; + pl.ks = kseq_init(fp); + } + if (rs && (flag & (HAF_RS_WRITE_LEN|HAF_RS_WRITE_SEQ))) + pl.rs_out = rs; + pl.flt_tab = flt_tab; + pl.opt = opt; + pl.flag = flag; + if (p0) { + pl.pt = p0, pl.create_new = 0; + assert(p0->k == opt->k && p0->pre == opt->pre); + } else if (c0) { + pl.ct = c0, pl.create_new = 0; + assert(c0->k == opt->k && c0->pre == opt->pre); + } else { + pl.create_new = 1; + pl.ct = ha_ct_init(opt->k, opt->pre, opt->bf_n_hash, opt->bf_shift); + } + kt_pipeline(3, worker_count, &pl, 3); + if (read_rs) { + destory_UC_Read(&pl.ucr); + } else { + kseq_destroy(pl.ks); + gzclose(fp); + } + return pl.ct; +} + +ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const void *flt_tab, All_reads *rs) +{ + int i; + yak_copt_t opt; + ha_ct_t *h = 0; + assert(!(flag & HAF_RS_WRITE_LEN) || !(flag & HAF_RS_WRITE_SEQ)); // not both + if (rs) { + if (flag & HAF_RS_WRITE_LEN) + init_All_reads(rs); + else if (flag & HAF_RS_WRITE_SEQ) + malloc_All_reads(rs); + } + yak_copt_init(&opt); + opt.k = asm_opt->k_mer_length; + opt.is_HPC = !asm_opt->no_HPC; + opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win; + opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift; + opt.n_thread = asm_opt->thread_num; + for (i = 0; i < asm_opt->num_reads; ++i) + h = yak_count(&opt, asm_opt->read_file_names[i], flag, p0, h, flt_tab, rs); + if (h && opt.bf_shift > 0) + ha_ct_destroy_bf(h); + return h; +} + +/*************************** + * High count filter table * + ***************************/ + +KHASHL_SET_INIT(static klib_unused, yak_ft_t, yak_ft, uint64_t, kh_hash_dummy, kh_eq_generic) + +static yak_ft_t *gen_hh(const ha_ct_t *h) +{ + int i; + yak_ft_t *hh; + hh = yak_ft_init(); + yak_ft_resize(hh, h->tot * 2); + for (i = 0; i < 1<pre; ++i) { + yak_ct_t *ht = h->h[i].h; + khint_t k; + for (k = 0; k < kh_end(ht); ++k) { + if (kh_exist(ht, k)) { + uint64_t y = kh_key(ht, k) >> h->pre << YAK_COUNTER_BITS | i; + int absent; + yak_ft_put(hh, y, &absent); + } + } + } + return hh; +} + +int ha_ft_isflt(const void *hh, uint64_t y) +{ + yak_ft_t *h = (yak_ft_t*)hh; + khint_t k; + k = yak_ft_get(h, y); + return k == kh_end(h)? 0 : 1; +} + +void ha_ft_destroy(void *h) +{ + if (h) yak_ft_destroy((yak_ft_t*)h); +} + +/************************* + * High-level interfaces * + *************************/ + +void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs) +{ + yak_ft_t *flt_tab; + int64_t cnt[YAK_N_COUNTS]; + int peak_hom, peak_het, cutoff; + ha_ct_t *h; + h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN, NULL, NULL, rs); + ha_ct_hist(h, cnt, asm_opt->thread_num); + peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het); + if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het); + cutoff = (int)(peak_hom * asm_opt->high_factor); + if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1; + ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num); + flt_tab = gen_hh(h); + ha_ct_destroy(h); + fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> filtered out %ld k-mers occurring %d or more times\n", __func__, + yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb(), (long)kh_size(flt_tab), cutoff); + return (void*)flt_tab; +} + +ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs) +{ + int64_t cnt[YAK_N_COUNTS], tot_cnt; + int peak_hom, peak_het, i, extra_flag1, extra_flag2; + ha_ct_t *ct; + ha_pt_t *pt; + if (read_from_store) { + extra_flag1 = extra_flag2 = HAF_RS_READ; + } else if (rs->total_reads == 0) { + extra_flag1 = HAF_RS_WRITE_LEN; + extra_flag2 = HAF_RS_WRITE_SEQ; + } else { + extra_flag1 = HAF_RS_WRITE_SEQ; + extra_flag2 = HAF_RS_READ; + } + ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs); + fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__, + yak_realtime(), yak_cpu_usage(), (long)ct->tot); + ha_ct_hist(ct, cnt, asm_opt->thread_num); + fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]); + peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het); + if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het); + if (flt_tab == 0) { + int cutoff = (int)(peak_hom * asm_opt->high_factor); + if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1; + ha_ct_shrink(ct, 2, cutoff, asm_opt->thread_num); + for (i = 2, tot_cnt = 0; i <= cutoff; ++i) tot_cnt += cnt[i] * i; + } else { + ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num); + for (i = 2, tot_cnt = 0; i <= YAK_MAX_COUNT - 1; ++i) tot_cnt += cnt[i] * i; + } + pt = ha_pt_gen(ct, asm_opt->thread_num); + ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, pt, flt_tab, rs); + assert((uint64_t)tot_cnt == pt->tot_pos); + //ha_pt_sort(pt, asm_opt->thread_num); + fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__, + yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos); + return pt; +} diff --git a/htab.h b/htab.h new file mode 100644 index 0000000..6298923 --- /dev/null +++ b/htab.h @@ -0,0 +1,103 @@ +#ifndef __HA_HTAB_H__ +#define __HA_HTAB_H__ +#define __STDC_LIMIT_MACROS +#include +#include "Process_Read.h" +#include "CommandLines.h" + +typedef struct { + uint64_t x; + 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; // actually it is not necessary to keep span in the index +} ha_idxpos_t; + +typedef struct { uint32_t n, m; ha_mz1_t *a; } ha_mz1_v; + +struct ha_pt_s; +typedef struct ha_pt_s ha_pt_t; + +struct ha_abuf_s; +typedef struct ha_abuf_s ha_abuf_t; + +extern const unsigned char seq_nt4_table[256]; +extern void *ha_flt_tab; +extern ha_pt_t *ha_idx; + +void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs); +int ha_ft_isflt(const void *hh, uint64_t y); +void ha_ft_destroy(void *h); + +ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs); +void ha_pt_destroy(ha_pt_t *h); +const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n); + +ha_abuf_t *ha_abuf_init(void); +void ha_abuf_destroy(ha_abuf_t *ab); +uint64_t ha_abuf_mem(const ha_abuf_t *ab); + +double yak_cputime(void); +void yak_reset_realtime(void); +double yak_realtime(void); +long yak_peakrss(void); +double yak_peakrss_in_gb(void); +double yak_cpu_usage(void); + +void trio_partition(void); + +void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf); +int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het); + +static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function +{ + key = (~key + (key << 21)) & mask; // key = (key << 21) - key - 1; + key = key ^ key >> 24; + key = ((key + (key << 3)) + (key << 8)) & mask; // key * 265 + key = key ^ key >> 14; + key = ((key + (key << 2)) + (key << 4)) & mask; // key * 21 + key = key ^ key >> 28; + key = (key + (key << 31)) & mask; + return key; +} + +static inline uint64_t yak_hash64_64(uint64_t key) +{ + key = ~key + (key << 21); + key = key ^ key >> 24; + key = (key + (key << 3)) + (key << 8); + key = key ^ key >> 14; + key = (key + (key << 2)) + (key << 4); + key = key ^ key >> 28; + key = key + (key << 31); + return key; +} + +static inline uint64_t yak_hash_long(uint64_t x[4]) +{ + int j = x[1] < x[3]? 0 : 1; + return yak_hash64_64(x[j<<1|0]) + yak_hash64_64(x[j<<1|1]); +} + +#define CALLOC(ptr, len) ((ptr) = (__typeof__(ptr))calloc((len), sizeof(*(ptr)))) +#define MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr)))) +#define REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr)))) + +#ifndef kroundup32 +#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x)) +#endif + +#ifndef kroundup64 +#define kroundup64(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, x|=(x)>>32, ++(x)) +#endif + +#ifndef klib_unused +#if (defined __clang__ && __clang_major__ >= 3) || (defined __GNUC__ && __GNUC__ >= 3) +#define klib_unused __attribute__ ((__unused__)) +#else +#define klib_unused +#endif +#endif /* klib_unused */ + +#endif // __YAK_H__ diff --git a/khash.h b/khash.h deleted file mode 100644 index 71b1e3e..0000000 --- a/khash.h +++ /dev/null @@ -1,669 +0,0 @@ -/* The MIT License - - Copyright (c) 2008, 2009, 2011 by Attractive Chaos - - Permission is hereby granted, free of charge, to any person obtaining - a copy of this software and associated documentation files (the - "Software"), to deal in the Software without restriction, including - without limitation the rights to use, copy, modify, merge, publish, - distribute, sublicense, and/or sell copies of the Software, and to - permit persons to whom the Software is furnished to do so, subject to - the following conditions: - - The above copyright notice and this permission notice shall be - included in all copies or substantial portions of the Software. - - THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF - MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS - BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN - ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN - CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE - SOFTWARE. -*/ - -/* - An example: - -#include "khash.h" -KHASH_MAP_INIT_INT(32, char) -int main() { - int ret, is_missing; - khiter_t k; - khash_t(32) *h = kh_init(32); - k = kh_put(32, h, 5, &ret); - kh_value(h, k) = 10; - k = kh_get(32, h, 10); - is_missing = (k == kh_end(h)); - k = kh_get(32, h, 5); - kh_del(32, h, k); - for (k = kh_begin(h); k != kh_end(h); ++k) - if (kh_exist(h, k)) kh_value(h, k) = 1; - kh_destroy(32, h); - return 0; -} -*/ - -/* - 2013-05-02 (0.2.8): - - * Use quadratic probing. When the capacity is power of 2, stepping function - i*(i+1)/2 guarantees to traverse each bucket. It is better than double - hashing on cache performance and is more robust than linear probing. - - In theory, double hashing should be more robust than quadratic probing. - However, my implementation is probably not for large hash tables, because - the second hash function is closely tied to the first hash function, - which reduce the effectiveness of double hashing. - - Reference: http://research.cs.vt.edu/AVresearch/hashing/quadratic.php - - 2011-12-29 (0.2.7): - - * Minor code clean up; no actual effect. - - 2011-09-16 (0.2.6): - - * The capacity is a power of 2. This seems to dramatically improve the - speed for simple keys. Thank Zilong Tan for the suggestion. Reference: - - - http://code.google.com/p/ulib/ - - http://nothings.org/computer/judy/ - - * Allow to optionally use linear probing which usually has better - performance for random input. Double hashing is still the default as it - is more robust to certain non-random input. - - * Added Wang's integer hash function (not used by default). This hash - function is more robust to certain non-random input. - - 2011-02-14 (0.2.5): - - * Allow to declare global functions. - - 2009-09-26 (0.2.4): - - * Improve portability - - 2008-09-19 (0.2.3): - - * Corrected the example - * Improved interfaces - - 2008-09-11 (0.2.2): - - * Improved speed a little in kh_put() - - 2008-09-10 (0.2.1): - - * Added kh_clear() - * Fixed a compiling error - - 2008-09-02 (0.2.0): - - * Changed to token concatenation which increases flexibility. - - 2008-08-31 (0.1.2): - - * Fixed a bug in kh_get(), which has not been tested previously. - - 2008-08-31 (0.1.1): - - * Added destructor -*/ - - -#ifndef __AC_KHASH_H -#define __AC_KHASH_H - -/*! - @header - - Generic hash table library. - */ - -#define AC_VERSION_KHASH_H "0.2.8" - -#include -#include -#include -#include - - -/* compiler specific configuration */ - -#if UINT_MAX == 0xffffffffu -typedef unsigned int khint32_t; -#elif ULONG_MAX == 0xffffffffu -typedef unsigned long khint32_t; -#endif - -#if ULONG_MAX == ULLONG_MAX -typedef unsigned long khint64_t; -#else -typedef unsigned long long khint64_t; -#endif - -#ifndef kh_inline -#ifdef _MSC_VER -#define kh_inline __inline -#else -#define kh_inline inline -#endif -#endif /* kh_inline */ - -#ifndef klib_unused -#if (defined __clang__ && __clang_major__ >= 3) || (defined __GNUC__ && __GNUC__ >= 3) -#define klib_unused __attribute__ ((__unused__)) -#else -#define klib_unused -#endif -#endif /* klib_unused */ - -typedef khint32_t khint_t; -typedef khint_t khiter_t; - -#define __ac_isempty(flag, i) ((flag[i>>4]>>((i&0xfU)<<1))&2) -#define __ac_isdel(flag, i) ((flag[i>>4]>>((i&0xfU)<<1))&1) -#define __ac_iseither(flag, i) ((flag[i>>4]>>((i&0xfU)<<1))&3) -#define __ac_set_isdel_false(flag, i) (flag[i>>4]&=~(1ul<<((i&0xfU)<<1))) -#define __ac_set_isempty_false(flag, i) (flag[i>>4]&=~(2ul<<((i&0xfU)<<1))) -#define __ac_set_isboth_false(flag, i) (flag[i>>4]&=~(3ul<<((i&0xfU)<<1))) -#define __ac_set_isdel_true(flag, i) (flag[i>>4]|=1ul<<((i&0xfU)<<1)) - -#define __ac_fsize(m) ((m) < 16? 1 : (m)>>4) - -#ifndef kroundup32 -#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x)) -#endif - -#ifndef kcalloc -#define kcalloc(N,Z) calloc(N,Z) -#endif -#ifndef kmalloc -#define kmalloc(Z) malloc(Z) -#endif -#ifndef krealloc -#define krealloc(P,Z) realloc(P,Z) -#endif -#ifndef kfree -#define kfree(P) free(P) -#endif - -static const double __ac_HASH_UPPER = 0.77; - -#define __KHASH_TYPE(name, khkey_t, khval_t) \ - typedef struct kh_##name##_s { \ - khint_t n_buckets, size, n_occupied, upper_bound; \ - khint32_t *flags; \ - khkey_t *keys; \ - khval_t *vals; \ - } kh_##name##_t; - -#define __KHASH_PROTOTYPES(name, khkey_t, khval_t) \ - extern kh_##name##_t *kh_init_##name(void); \ - extern void kh_destroy_##name(kh_##name##_t *h); \ - extern void kh_clear_##name(kh_##name##_t *h); \ - extern khint_t kh_get_##name(const kh_##name##_t *h, khkey_t key); \ - extern int kh_resize_##name(kh_##name##_t *h, khint_t new_n_buckets); \ - extern khint_t kh_put_##name(kh_##name##_t *h, khkey_t key, int *ret); \ - extern void kh_del_##name(kh_##name##_t *h, khint_t x);\ - extern void kh_write_##name(kh_##name##_t *h, FILE* fp);\ - extern void kh_load_##name(kh_##name##_t *h, FILE* fp); - -#define __KHASH_IMPL(name, SCOPE, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) \ - SCOPE kh_##name##_t *kh_init_##name(void) { \ - return (kh_##name##_t*)kcalloc(1, sizeof(kh_##name##_t)); \ - } \ - SCOPE void kh_destroy_##name(kh_##name##_t *h) \ - { \ - if (h) { \ - kfree((void *)h->keys); kfree(h->flags); \ - kfree((void *)h->vals); \ - kfree(h); \ - } \ - } \ - SCOPE void kh_clear_##name(kh_##name##_t *h) \ - { \ - if (h && h->flags) { \ - memset(h->flags, 0xaa, __ac_fsize(h->n_buckets) * sizeof(khint32_t)); \ - h->size = h->n_occupied = 0; \ - } \ - } \ - SCOPE khint_t kh_get_##name(const kh_##name##_t *h, khkey_t key) \ - { \ - if (h->n_buckets) { \ - khint_t k, i, last, mask, step = 0; \ - mask = h->n_buckets - 1; \ - k = __hash_func(key); i = k & mask; \ - last = i; \ - while (!__ac_isempty(h->flags, i) && (__ac_isdel(h->flags, i) || !__hash_equal(h->keys[i], key))) { \ - i = (i + (++step)) & mask; \ - if (i == last) return h->n_buckets; \ - } \ - return __ac_iseither(h->flags, i)? h->n_buckets : i; \ - } else return 0; \ - } \ - SCOPE int kh_resize_##name(kh_##name##_t *h, khint_t new_n_buckets) \ - { /* This function uses 0.25*n_buckets bytes of working space instead of [sizeof(key_t+val_t)+.25]*n_buckets. */ \ - khint32_t *new_flags = 0; \ - khint_t j = 1; \ - { \ - kroundup32(new_n_buckets); \ - if (new_n_buckets < 4) new_n_buckets = 4; \ - if (h->size >= (khint_t)(new_n_buckets * __ac_HASH_UPPER + 0.5)) j = 0; /* requested size is too small */ \ - else { /* hash table size to be changed (shrink or expand); rehash */ \ - new_flags = (khint32_t*)kmalloc(__ac_fsize(new_n_buckets) * sizeof(khint32_t)); \ - if (!new_flags) return -1; \ - memset(new_flags, 0xaa, __ac_fsize(new_n_buckets) * sizeof(khint32_t)); \ - if (h->n_buckets < new_n_buckets) { /* expand */ \ - khkey_t *new_keys = (khkey_t*)krealloc((void *)h->keys, new_n_buckets * sizeof(khkey_t)); \ - if (!new_keys) { kfree(new_flags); return -1; } \ - h->keys = new_keys; \ - if (kh_is_map) { \ - khval_t *new_vals = (khval_t*)krealloc((void *)h->vals, new_n_buckets * sizeof(khval_t)); \ - if (!new_vals) { kfree(new_flags); return -1; } \ - h->vals = new_vals; \ - } \ - } /* otherwise shrink */ \ - } \ - } \ - if (j) { /* rehashing is needed */ \ - for (j = 0; j != h->n_buckets; ++j) { \ - if (__ac_iseither(h->flags, j) == 0) { \ - khkey_t key = h->keys[j]; \ - khval_t val; \ - khint_t new_mask; \ - new_mask = new_n_buckets - 1; \ - if (kh_is_map) val = h->vals[j]; \ - __ac_set_isdel_true(h->flags, j); \ - while (1) { /* kick-out process; sort of like in Cuckoo hashing */ \ - khint_t k, i, step = 0; \ - k = __hash_func(key); \ - i = k & new_mask; \ - while (!__ac_isempty(new_flags, i)) i = (i + (++step)) & new_mask; \ - __ac_set_isempty_false(new_flags, i); \ - if (i < h->n_buckets && __ac_iseither(h->flags, i) == 0) { /* kick out the existing element */ \ - { khkey_t tmp = h->keys[i]; h->keys[i] = key; key = tmp; } \ - if (kh_is_map) { khval_t tmp = h->vals[i]; h->vals[i] = val; val = tmp; } \ - __ac_set_isdel_true(h->flags, i); /* mark it as deleted in the old hash table */ \ - } else { /* write the element and jump out of the loop */ \ - h->keys[i] = key; \ - if (kh_is_map) h->vals[i] = val; \ - break; \ - } \ - } \ - } \ - } \ - if (h->n_buckets > new_n_buckets) { /* shrink the hash table */ \ - h->keys = (khkey_t*)krealloc((void *)h->keys, new_n_buckets * sizeof(khkey_t)); \ - if (kh_is_map) h->vals = (khval_t*)krealloc((void *)h->vals, new_n_buckets * sizeof(khval_t)); \ - } \ - kfree(h->flags); /* free the working space */ \ - h->flags = new_flags; \ - h->n_buckets = new_n_buckets; \ - h->n_occupied = h->size; \ - h->upper_bound = (khint_t)(h->n_buckets * __ac_HASH_UPPER + 0.5); \ - } \ - return 0; \ - } \ - SCOPE khint_t kh_put_##name(kh_##name##_t *h, khkey_t key, int *ret) \ - { \ - khint_t x; \ - if (h->n_occupied >= h->upper_bound) { /* update the hash table */ \ - if (h->n_buckets > (h->size<<1)) { \ - if (kh_resize_##name(h, h->n_buckets - 1) < 0) { /* clear "deleted" elements */ \ - *ret = -1; return h->n_buckets; \ - } \ - } else if (kh_resize_##name(h, h->n_buckets + 1) < 0) { /* expand the hash table */ \ - *ret = -1; return h->n_buckets; \ - } \ - } /* TODO: to implement automatically shrinking; resize() already support shrinking */ \ - { \ - khint_t k, i, site, last, mask = h->n_buckets - 1, step = 0; \ - x = site = h->n_buckets; k = __hash_func(key); i = k & mask; \ - if (__ac_isempty(h->flags, i)) x = i; /* for speed up */ \ - else { \ - last = i; \ - while (!__ac_isempty(h->flags, i) && (__ac_isdel(h->flags, i) || !__hash_equal(h->keys[i], key))) { \ - if (__ac_isdel(h->flags, i)) site = i; \ - i = (i + (++step)) & mask; \ - if (i == last) { x = site; break; } \ - } \ - if (x == h->n_buckets) { \ - if (__ac_isempty(h->flags, i) && site != h->n_buckets) x = site; \ - else x = i; \ - } \ - } \ - } \ - if (__ac_isempty(h->flags, x)) { /* not present at all */ \ - h->keys[x] = key; \ - __ac_set_isboth_false(h->flags, x); \ - ++h->size; ++h->n_occupied; \ - *ret = 1; \ - } else if (__ac_isdel(h->flags, x)) { /* deleted */ \ - h->keys[x] = key; \ - __ac_set_isboth_false(h->flags, x); \ - ++h->size; \ - *ret = 2; \ - } else *ret = 0; /* Don't touch h->keys[x] if present and not deleted */ \ - return x; \ - } \ - SCOPE void kh_del_##name(kh_##name##_t *h, khint_t x) \ - { \ - if (x != h->n_buckets && !__ac_iseither(h->flags, x)) { \ - __ac_set_isdel_true(h->flags, x); \ - --h->size; \ - } \ - } \ - SCOPE void kh_write_##name(kh_##name##_t *h, FILE* fp)\ - {\ - fwrite(&(h->n_buckets), sizeof(khint_t), 1, fp);\ - fwrite(&(h->size), sizeof(khint_t), 1, fp);\ - fwrite(&(h->n_occupied), sizeof(khint_t), 1, fp);\ - fwrite(&(h->upper_bound), sizeof(khint_t), 1, fp);\ - if (h->n_buckets)\ - {\ - fwrite(h->flags, sizeof(khint32_t), __ac_fsize(h->n_buckets), fp);\ - fwrite(h->keys, sizeof(khkey_t), h->n_buckets, fp);\ - fwrite(h->vals, sizeof(khval_t), h->n_buckets, fp);\ - }\ - } \ - SCOPE void kh_load_##name(kh_##name##_t *h, FILE* fp)\ - {\ - int f_flag;\ - f_flag = fread(&(h->n_buckets), sizeof(khint_t), 1, fp);\ - f_flag += fread(&(h->size), sizeof(khint_t), 1, fp);\ - f_flag += fread(&(h->n_occupied), sizeof(khint_t), 1, fp);\ - f_flag += fread(&(h->upper_bound), sizeof(khint_t), 1, fp);\ - if (h->n_buckets)\ - {\ - h->flags = (khint32_t*)kmalloc(__ac_fsize(h->n_buckets) * sizeof(khint32_t));\ - f_flag += fread(h->flags, sizeof(khint32_t), __ac_fsize(h->n_buckets), fp);\ - h->keys = (khkey_t*)kmalloc(sizeof(khkey_t)*h->n_buckets);\ - f_flag += fread(h->keys, sizeof(khkey_t), h->n_buckets, fp);\ - h->vals = (khval_t*)kmalloc(sizeof(khval_t)*h->n_buckets);\ - f_flag += fread(h->vals, sizeof(khval_t), h->n_buckets, fp);\ - }\ - } - -#define KHASH_DECLARE(name, khkey_t, khval_t) \ - __KHASH_TYPE(name, khkey_t, khval_t) \ - __KHASH_PROTOTYPES(name, khkey_t, khval_t) - -#define KHASH_INIT2(name, SCOPE, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) \ - __KHASH_TYPE(name, khkey_t, khval_t) \ - __KHASH_IMPL(name, SCOPE, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) - -#define KHASH_INIT(name, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) \ - KHASH_INIT2(name, static kh_inline klib_unused, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) - -/* --- BEGIN OF HASH FUNCTIONS --- */ - -/*! @function - @abstract Integer hash function - @param key The integer [khint32_t] - @return The hash value [khint_t] - */ -#define kh_int_hash_func(key) (khint32_t)(key) -/*! @function - @abstract Integer comparison function - */ -#define kh_int_hash_equal(a, b) ((a) == (b)) -/*! @function - @abstract 64-bit integer hash function - @param key The integer [khint64_t] - @return The hash value [khint_t] - */ -#define kh_int64_hash_func(key) (khint32_t)((key)>>33^(key)^(key)<<11) -/*! @function - @abstract 64-bit integer comparison function - */ -#define kh_int64_hash_equal(a, b) ((a) == (b)) -/*! @function - @abstract const char* hash function - @param s Pointer to a null terminated string - @return The hash value - */ -static kh_inline khint_t __ac_X31_hash_string(const char *s) -{ - khint_t h = (khint_t)*s; - if (h) for (++s ; *s; ++s) h = (h << 5) - h + (khint_t)*s; - return h; -} -/*! @function - @abstract Another interface to const char* hash function - @param key Pointer to a null terminated string [const char*] - @return The hash value [khint_t] - */ -#define kh_str_hash_func(key) __ac_X31_hash_string(key) -/*! @function - @abstract Const char* comparison function - */ -#define kh_str_hash_equal(a, b) (strcmp(a, b) == 0) - -static kh_inline khint_t __ac_Wang_hash(khint_t key) -{ - key += ~(key << 15); - key ^= (key >> 10); - key += (key << 3); - key ^= (key >> 6); - key += ~(key << 11); - key ^= (key >> 16); - return key; -} -#define kh_int_hash_func2(key) __ac_Wang_hash((khint_t)key) - -/* --- END OF HASH FUNCTIONS --- */ - -/* Other convenient macros... */ - -/*! - @abstract Type of the hash table. - @param name Name of the hash table [symbol] - */ -#define khash_t(name) kh_##name##_t - -/*! @function - @abstract Initiate a hash table. - @param name Name of the hash table [symbol] - @return Pointer to the hash table [khash_t(name)*] - */ -#define kh_init(name) kh_init_##name() - -/*! @function - @abstract Destroy a hash table. - @param name Name of the hash table [symbol] - @param h Pointer to the hash table [khash_t(name)*] - */ -#define kh_destroy(name, h) kh_destroy_##name(h) - -/*! @function - @abstract Reset a hash table without deallocating memory. - @param name Name of the hash table [symbol] - @param h Pointer to the hash table [khash_t(name)*] - */ -#define kh_clear(name, h) kh_clear_##name(h) - -/*! @function - @abstract Resize a hash table. - @param name Name of the hash table [symbol] - @param h Pointer to the hash table [khash_t(name)*] - @param s New size [khint_t] - */ -#define kh_resize(name, h, s) kh_resize_##name(h, s) - -/*! @function - @abstract Insert a key to the hash table. - @param name Name of the hash table [symbol] - @param h Pointer to the hash table [khash_t(name)*] - @param k Key [type of keys] - @param r Extra return code: -1 if the operation failed; - 0 if the key is present in the hash table; - 1 if the bucket is empty (never used); 2 if the element in - the bucket has been deleted [int*] - @return Iterator to the inserted element [khint_t] - */ -#define kh_put(name, h, k, r) kh_put_##name(h, k, r) - -/*! @function - @abstract Retrieve a key from the hash table. - @param name Name of the hash table [symbol] - @param h Pointer to the hash table [khash_t(name)*] - @param k Key [type of keys] - @return Iterator to the found element, or kh_end(h) if the element is absent [khint_t] - */ -#define kh_get(name, h, k) kh_get_##name(h, k) - -/*! @function - @abstract Remove a key from the hash table. - @param name Name of the hash table [symbol] - @param h Pointer to the hash table [khash_t(name)*] - @param k Iterator to the element to be deleted [khint_t] - */ -#define kh_del(name, h, k) kh_del_##name(h, k) - -/*! @function - @abstract Test whether a bucket contains data. - @param h Pointer to the hash table [khash_t(name)*] - @param x Iterator to the bucket [khint_t] - @return 1 if containing data; 0 otherwise [int] - */ -#define kh_exist(h, x) (!__ac_iseither((h)->flags, (x))) - -/*! @function - @abstract Get key given an iterator - @param h Pointer to the hash table [khash_t(name)*] - @param x Iterator to the bucket [khint_t] - @return Key [type of keys] - */ -#define kh_key(h, x) ((h)->keys[x]) - -/*! @function - @abstract Get value given an iterator - @param h Pointer to the hash table [khash_t(name)*] - @param x Iterator to the bucket [khint_t] - @return Value [type of values] - @discussion For hash sets, calling this results in segfault. - */ -#define kh_val(h, x) ((h)->vals[x]) - -/*! @function - @abstract Alias of kh_val() - */ -#define kh_value(h, x) ((h)->vals[x]) - -/*! @function - @abstract Get the start iterator - @param h Pointer to the hash table [khash_t(name)*] - @return The start iterator [khint_t] - */ -#define kh_begin(h) (khint_t)(0) - -/*! @function - @abstract Get the end iterator - @param h Pointer to the hash table [khash_t(name)*] - @return The end iterator [khint_t] - */ -#define kh_end(h) ((h)->n_buckets) - -/*! @function - @abstract Get the number of elements in the hash table - @param h Pointer to the hash table [khash_t(name)*] - @return Number of elements in the hash table [khint_t] - */ -#define kh_size(h) ((h)->size) - -/*! @function - @abstract Get the number of buckets in the hash table - @param h Pointer to the hash table [khash_t(name)*] - @return Number of buckets in the hash table [khint_t] - */ -#define kh_n_buckets(h) ((h)->n_buckets) - -/*! @function - @abstract Iterate over the entries in the hash table - @param h Pointer to the hash table [khash_t(name)*] - @param kvar Variable to which key will be assigned - @param vvar Variable to which value will be assigned - @param code Block of code to execute - */ -#define kh_foreach(h, kvar, vvar, code) { khint_t __i; \ - for (__i = kh_begin(h); __i != kh_end(h); ++__i) { \ - if (!kh_exist(h,__i)) continue; \ - (kvar) = kh_key(h,__i); \ - (vvar) = kh_val(h,__i); \ - code; \ - } } - -/*! @function - @abstract Iterate over the values in the hash table - @param h Pointer to the hash table [khash_t(name)*] - @param vvar Variable to which value will be assigned - @param code Block of code to execute - */ -#define kh_foreach_value(h, vvar, code) { khint_t __i; \ - for (__i = kh_begin(h); __i != kh_end(h); ++__i) { \ - if (!kh_exist(h,__i)) continue; \ - (vvar) = kh_val(h,__i); \ - code; \ - } } - -/* More convenient interfaces */ - -/*! @function - @abstract Instantiate a hash set containing integer keys - @param name Name of the hash table [symbol] - */ -#define KHASH_SET_INIT_INT(name) \ - KHASH_INIT(name, khint32_t, char, 0, kh_int_hash_func, kh_int_hash_equal) - -/*! @function - @abstract Instantiate a hash map containing integer keys - @param name Name of the hash table [symbol] - @param khval_t Type of values [type] - */ -#define KHASH_MAP_INIT_INT(name, khval_t) \ - KHASH_INIT(name, khint32_t, khval_t, 1, kh_int_hash_func, kh_int_hash_equal) - -/*! @function - @abstract Instantiate a hash set containing 64-bit integer keys - @param name Name of the hash table [symbol] - */ -#define KHASH_SET_INIT_INT64(name) \ - KHASH_INIT(name, khint64_t, char, 0, kh_int64_hash_func, kh_int64_hash_equal) - -/*! @function - @abstract Instantiate a hash map containing 64-bit integer keys - @param name Name of the hash table [symbol] - @param khval_t Type of values [type] - */ -#define KHASH_MAP_INIT_INT64(name, khval_t) \ - KHASH_INIT(name, khint64_t, khval_t, 1, kh_int64_hash_func, kh_int64_hash_equal) - -typedef const char *kh_cstr_t; -/*! @function - @abstract Instantiate a hash map containing const char* keys - @param name Name of the hash table [symbol] - */ -#define KHASH_SET_INIT_STR(name) \ - KHASH_INIT(name, kh_cstr_t, char, 0, kh_str_hash_func, kh_str_hash_equal) - -/*! @function - @abstract Instantiate a hash map containing const char* keys - @param name Name of the hash table [symbol] - @param khval_t Type of values [type] - */ -#define KHASH_MAP_INIT_STR(name, khval_t) \ - KHASH_INIT(name, kh_cstr_t, khval_t, 1, kh_str_hash_func, kh_str_hash_equal) - - - - -#define kh_write(name, h, fp) kh_write_##name(h, fp) - -#define kh_load(name, h, fp) kh_load_##name(h, fp) - - -#endif /* __AC_KHASH_H */ diff --git a/khashl.h b/khashl.h index 2d0f413..9381850 100644 --- a/khashl.h +++ b/khashl.h @@ -276,6 +276,7 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26 KHASHL_INIT(KH_LOCAL, HType, prefix##_m, HType##_m_bucket_t, prefix##_m_hash, prefix##_m_eq) \ SCOPE HType *prefix##_init(void) { return prefix##_m_init(); } \ SCOPE void prefix##_destroy(HType *h) { prefix##_m_destroy(h); } \ + SCOPE void prefix##_resize(HType *h, khint_t new_n_buckets) { prefix##_m_resize(h, new_n_buckets); } \ SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_m_bucket_t t; t.key = key; return prefix##_m_getp(h, &t); } \ SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_m_del(h, k); } \ SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_m_bucket_t t; t.key = key; return prefix##_m_putp(h, &t, absent); } diff --git a/kmer.cpp b/kmer.cpp deleted file mode 100644 index 0d234cf..0000000 --- a/kmer.cpp +++ /dev/null @@ -1,195 +0,0 @@ -#include -#include -#include "kmer.h" - -void init_HPC_seq(HPC_seq* seq, char* str, long long l) -{ - seq->i = 0; - seq->l = l; - seq->N_occ = 0; - seq->str = str; -} - -void init_Hash_code(Hash_code* code) -{ - code->x[0] = 0; - code->x[1] = 0; -} - - - -void init_small_hash_table(small_hash_table* x) -{ - x->size = 0; - x->buffer = NULL; - x->length = 0; -} - -void clear_small_hash_table(small_hash_table* x) -{ - x->length = 0; -} - -void resize_small_hash_table(small_hash_table* x, uint64_t size) -{ - if(size > x->size) - { - x->size = size; - x->buffer = (k_v*)realloc(x->buffer, x->size*sizeof(k_v)); - } -} - -void destory_small_hash_table(small_hash_table* x) -{ - free(x->buffer); -} - - -void add_small_hash_table(small_hash_table* x, k_v* element) -{ - if(x->length + 1 > x->size) - { - x->size = (x->length + 1) * 2; - x->buffer = (k_v*)realloc(x->buffer, x->size*sizeof(k_v)); - } - - x->buffer[x->length] = (*element); - x->length++; -} - -//x > y, return 1; x < y, return -1, x == y, return 0 -int compare_k_mer(k_v* x, k_v* y) -{ - if(x->key.x[1] != y->key.x[1]) - { - return x->key.x[1] > y->key.x[1] ? 1: -1; - } - else - { - if(x->key.x[0] != y->key.x[0]) - { - return x->key.x[0] > y->key.x[0] ? 1: -1; - } - else - { - return 0; - } - } - -} - -int cmp_k_mer_kv(const void * a, const void * b) -{ - int flag = compare_k_mer((k_v*)a, (k_v*)b); - - if(flag == 0) - { - if ((*(k_v*)a).value != (*(k_v*)b).value) - { - return (*(k_v*)a).value > (*(k_v*)b).value ? 1: -1; - } - else - { - return 0; - } - - } - else - { - return flag; - } -} - -void sort_small_hash_table(small_hash_table* x) -{ - qsort(x->buffer, x->length, sizeof(k_v), cmp_k_mer_kv); -} - - -inline long long firstEqual(k_v* arr, long long arrLen, k_v* key) -{ - long long L = 0, R = arrLen - 1; //[L, R] - long long mid; - int flag; - while( L <= R) - { - mid = L + (R - L)/2; - - flag = compare_k_mer(&(arr[mid]), key); - - ///arr[mid] >= key - if(flag >= 0) - { - R = mid - 1; - } - else - { - L = mid + 1; - } - } - - - if(L < arrLen && (flag = compare_k_mer(&(arr[L]), key) == 0)) - { - return L; - } - - return -1; -} - -inline long long lastEqual(k_v* arr, long long arrLen, k_v* key) -{ - long long L = 0, R = arrLen - 1; //[L, R] - long long mid; - int flag; - while( L <= R) - { - mid = L + (R - L)/2; - flag = compare_k_mer(&(arr[mid]), key); - ///arr[mid] <= key - if(flag <= 0) - { - L = mid + 1; - } - else - { - R = mid - 1; - } - } - - if(R >= 0 && ((flag = compare_k_mer(&(arr[R]), key)) == 0)) - { - return R; - } - - return -1; -} - -int query_small_hash_table(small_hash_table* target, k_v* query, long long* l_end, long long* r_end) -{ - (*l_end) = -1; - (*r_end) = -1; - long long left_end; - long long right_end; - - left_end = firstEqual(target->buffer, target->length, query); - - if(left_end != -1) - { - right_end = lastEqual(target->buffer + left_end, target->length - left_end, query) + left_end; - - (*l_end) = left_end; - (*r_end) = right_end; - - if(right_end == -1) - { - fprintf(stderr, "error\n"); - } - - - return right_end - left_end + 1; - } - - return 0; -} - diff --git a/kmer.h b/kmer.h deleted file mode 100644 index 016a531..0000000 --- a/kmer.h +++ /dev/null @@ -1,114 +0,0 @@ -#ifndef __KMER__ -#define __KMER__ -#include "Process_Read.h" - -///#define ALL (0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff) -#define ALL (0xffffffffffffffff) -/****************************may have bugs********************************/ -#define SAFE_SHIFT(k) k & ((k < 64)?ALL:0) -/****************************may have bugs********************************/ - - - - - - -typedef struct -{ - //can represent at most 64-mer - uint64_t x[2]; -} Hash_code; - -typedef struct { - Hash_code key; ///k-mer itself - uint64_t value; ///offset -} k_v; - -typedef struct { - k_v* buffer; - uint32_t size; - uint32_t length; -} small_hash_table; - -void init_small_hash_table(small_hash_table* x); -void clear_small_hash_table(small_hash_table* x); -void resize_small_hash_table(small_hash_table* x, uint64_t size); -void destory_small_hash_table(small_hash_table* x); -void add_small_hash_table(small_hash_table* x, k_v* element); -void sort_small_hash_table(small_hash_table* x); -int compare_k_mer(k_v* x, k_v* y); -int query_small_hash_table(small_hash_table* target, k_v* query, long long* l_end, long long* r_end); - - -typedef struct -{ - char* str; - long long l; - long long i; - long long N_occ; - -} HPC_seq; - - -inline uint64_t get_HPC_code(HPC_seq* seq, uint64_t* end_pos) -{ - - if(seq->i < seq ->l) - { - uint8_t code = seq_nt6_table[(uint8_t)seq->str[seq->i]]; - - (*end_pos) = seq->i; - - for (; seq->i < seq->l; seq->i++) - { - ///number of Ns - if (seq_nt6_table[(uint8_t)seq->str[seq->i]] >= 4) - { - seq->N_occ++; - } - - if (seq_nt6_table[(uint8_t)seq->str[seq->i]] != code) - { - break; - } - } - - return (uint64_t)code; - } - else - { - ///end - return 6; - } - -} - -inline void k_mer_append(Hash_code* code, uint64_t c, int k) -{ - - uint64_t mask = ALL >> (64 -k); - - code->x[0] = ((code->x[0]<<1) | (c&1)) & mask; - code->x[1] = ((code->x[1]<<1) | (c>>1)) & mask; -} - -inline void Hashcode_to_string(Hash_code* code, char* str, int k) -{ - uint8_t c; - int i; - for (i = 0; i < k; i++) - { - c = (code->x[1] >> (k - i - 1)) & ((uint64_t)1); - c = c << 1; - c = c | ((code->x[0] >> (k - i - 1)) & ((uint64_t)1)); - - str[i] = s_H[c]; - } - -} - -void init_HPC_seq(HPC_seq* seq, char* str, long long l); -void init_Hash_code(Hash_code* code); - - -#endif \ No newline at end of file diff --git a/main.cpp b/main.cpp index 76fcfd8..ccc4ed1 100644 --- a/main.cpp +++ b/main.cpp @@ -4,16 +4,20 @@ #include "Process_Read.h" #include "Assembly.h" #include "Levenshtein_distance.h" +#include "htab.h" int main(int argc, char *argv[]) { + int i, ret; + yak_reset_realtime(); init_opt(&asm_opt); - if (!CommandLine_process(argc, argv, &asm_opt)) return 1; - - Correct_Reads(asm_opt.number_of_round); - + ret = ha_assemble(); destory_opt(&asm_opt); - - return 0; + fprintf(stderr, "[M::%s] Version: %s\n", __func__, HA_VERSION); + fprintf(stderr, "[M::%s] CMD:", __func__); + for (i = 0; i < argc; ++i) + fprintf(stderr, " %s", argv[i]); + fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, yak_realtime(), yak_cputime(), yak_peakrss_in_gb()); + return ret; } diff --git a/sketch.cpp b/sketch.cpp new file mode 100644 index 0000000..c3546ed --- /dev/null +++ b/sketch.cpp @@ -0,0 +1,110 @@ +#include +#include +#include +#include +#include "kvec.h" +#include "htab.h" + +typedef struct { // a simplified version of kdq + int front, count; + int a[64]; +} tiny_queue_t; + +static inline void tq_push(tiny_queue_t *q, int x) +{ + q->a[((q->count++) + q->front) & 0x3f] = x; +} + +static inline int tq_shift(tiny_queue_t *q) +{ + int x; + if (q->count == 0) return -1; + x = q->a[q->front++]; + q->front &= 0x3f; + --q->count; + return x; +} + +/** + * Find symmetric (w,k)-minimizers on a DNA sequence + * + * @param str DNA sequence + * @param len length of $str + * @param w find a minimizer for every $w consecutive k-mers + * @param k k-mer size + * @param rid reference ID; will be copied to the output $p array + * @param is_hpc homopolymer-compressed or not + * @param p minimizers + */ +void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf) +{ + static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 }; + uint64_t shift1 = k - 1, mask = (1ULL< 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63)); + memset(buf, 0xff, w * 16); + memset(&tq, 0, sizeof(tiny_queue_t)); + kv_resize(ha_mz1_t, *p, p->n + len/w); + + for (i = l = buf_pos = min_pos = 0; i < len; ++i) { + int c = seq_nt4_table[(uint8_t)str[i]]; + ha_mz1_t info = dummy; + if (c < 4) { // not an ambiguous base + int z; + if (is_hpc) { + int skip_len = 1; + if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) { + for (skip_len = 2; i + skip_len < len; ++skip_len) + if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c) + break; + i += skip_len - 1; // put $i at the end of the current homopolymer run + } + tq_push(&tq, skip_len); + kmer_span += skip_len; + if (tq.count > k) kmer_span -= tq_shift(&tq); + } else kmer_span = l + 1 < k? l + 1 : k; + kmer[0] = (kmer[0] << 1 | (c&1)) & mask; // forward k-mer + kmer[1] = (kmer[1] << 1 | (c>>1)) & mask; + kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; // reverse k-mer + kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1; + if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand + z = kmer[1] < kmer[3]? 0 : 1; // strand + ++l; + if (l >= k && kmer_span < 256) { + uint64_t y; + y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]); + if (hf == 0 || ha_ft_isflt(hf, y) == 0) + info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span; + } + } else l = 0, tq.count = tq.front = 0, kmer_span = 0; + buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below + if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet + for (j = buf_pos + 1; j < w; ++j) + if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]); + for (j = 0; j < buf_pos; ++j) + if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]); + } + if (info.x <= min.x) { // a new minimum; then write the old min + if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); + min = info, min_pos = buf_pos; + } else if (buf_pos == min_pos) { // old min has moved outside the window + if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min); + for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers + if (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer + for (j = 0; j <= buf_pos; ++j) + if (min.x >= buf[j].x) min = buf[j], min_pos = j; + if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers + for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted + if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]); + for (j = 0; j <= buf_pos; ++j) + if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]); + } + } + if (++buf_pos == w) buf_pos = 0; + } + if (min.x != UINT64_MAX) + kv_push(ha_mz1_t, *p, min); +} diff --git a/sys.cpp b/sys.cpp new file mode 100644 index 0000000..87b5ac2 --- /dev/null +++ b/sys.cpp @@ -0,0 +1,53 @@ +#include +#include +#include "htab.h" + +int yak_verbose = 3; + +static double yak_realtime0; + +double yak_cputime(void) +{ + struct rusage r; + getrusage(RUSAGE_SELF, &r); + return r.ru_utime.tv_sec + r.ru_stime.tv_sec + 1e-6 * (r.ru_utime.tv_usec + r.ru_stime.tv_usec); +} + +static inline double yak_realtime_core(void) +{ + struct timeval tp; + struct timezone tzp; + gettimeofday(&tp, &tzp); + return tp.tv_sec + tp.tv_usec * 1e-6; +} + +void yak_reset_realtime(void) +{ + yak_realtime0 = yak_realtime_core(); +} + +double yak_realtime(void) +{ + return yak_realtime_core() - yak_realtime0; +} + +long yak_peakrss(void) +{ + struct rusage r; + getrusage(RUSAGE_SELF, &r); +#ifdef __linux__ + return r.ru_maxrss * 1024; +#else + return r.ru_maxrss; +#endif +} + +double yak_peakrss_in_gb(void) +{ + return yak_peakrss() / 1073741824.0; +} + +double yak_cpu_usage(void) +{ + return (yak_cputime() + 1e-9) / (yak_realtime() + 1e-9); +}