From 968b4caef95e3c6daa668eafe7a67eb98b1c3095 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 2 Apr 2020 21:55:00 -0400 Subject: [PATCH] r181: limit the max number of chains --- Assembly.cpp | 8 ++++---- CommandLines.cpp | 4 +++- CommandLines.h | 3 ++- Hash_Table.cpp | 6 ++---- Hash_Table.h | 16 ++++++++-------- Overlaps.cpp | 16 ++++++++-------- Trio.cpp | 4 ---- anchor.cpp | 38 ++++++++++++++++++++++++++++++++++---- htab.h | 4 ++++ 9 files changed, 65 insertions(+), 34 deletions(-) diff --git a/Assembly.cpp b/Assembly.cpp index a32204f..ef0d952 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -10,7 +10,7 @@ #include "Output.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 band_width_threshold, int keep_whole_chain); +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); All_reads R_INF; pthread_mutex_t statistics; @@ -423,7 +423,7 @@ void* Overlap_calculate_heap_merge(void* arg) { //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, &l, 0.02, 1); - ha_get_new_candidates(ab, i, &g_read, &overlap_list, &l, 0.02, 1); + ha_get_new_candidates(ab, i, &g_read, &overlap_list, &l, 0.02, asm_opt.max_n_chain, 1); clear_Cigar_record(¤t_cigar); clear_Round2_alignment(&second_round); @@ -543,7 +543,7 @@ void* Output_related_reads(void* arg) { //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, &l, 0.02, 1); - ha_get_new_candidates(ab, i, &g_read, &overlap_list, &l, 0.02, 1); + ha_get_new_candidates(ab, i, &g_read, &overlap_list, &l, 0.02, asm_opt.max_n_chain, 1); fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i), Get_NAME((R_INF), i)); @@ -1109,7 +1109,7 @@ void* Final_overlap_calculate_heap_merge(void* arg) { //get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0); - ha_get_new_candidates(ab, i, &g_read, &overlap_list, &l, 0.001, 0); + ha_get_new_candidates(ab, i, &g_read, &overlap_list, &l, 0.001, asm_opt.max_n_chain, 0); /** correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon, diff --git a/CommandLines.cpp b/CommandLines.cpp index 2ed649a..c6260a4 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -69,6 +69,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->high_factor = 5.0f; asm_opt->no_HPC = 0; asm_opt->no_kmer_flt = 0; + asm_opt->max_n_chain = 1000; asm_opt->k_mer_min_freq = 3; asm_opt->k_mer_max_freq = 66; asm_opt->load_index_from_disk = 1; @@ -300,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:lw:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:F", 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); @@ -321,6 +322,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) 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 a1c7424..ad150d0 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.3.0-dirty-r179" +#define HA_VERSION "0.3.0-dirty-r181" #define VERBOSE 0 #define VERBOSE_GFA 0 @@ -22,6 +22,7 @@ typedef struct { 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; diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 1c4a73e..e8bfd77 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -816,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; } diff --git a/Hash_Table.h b/Hash_Table.h index fff5945..16b14c8 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -67,15 +67,15 @@ typedef struct typedef struct { window_list* buffer; - long long length; - long long size; + int32_t length; + int32_t size; } window_list_alloc; typedef struct { uint64_t* buffer; - uint64_t length; - uint64_t size; + uint32_t length; + uint32_t size; } Fake_Cigar; typedef struct @@ -96,12 +96,12 @@ typedef struct uint32_t align_length; uint8_t is_match; uint8_t without_large_indel; + 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; @@ -180,7 +180,7 @@ static 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); } diff --git a/Overlaps.cpp b/Overlaps.cpp index fbd01d6..e74f1b3 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -4333,7 +4333,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)) @@ -4346,7 +4346,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]) @@ -9028,7 +9028,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); @@ -9042,7 +9042,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; @@ -9123,7 +9123,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************************************/ @@ -9172,7 +9172,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************************************/ @@ -19547,7 +19547,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); @@ -19578,7 +19578,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); diff --git a/Trio.cpp b/Trio.cpp index c198759..7e9d974 100644 --- a/Trio.cpp +++ b/Trio.cpp @@ -20,10 +20,6 @@ #define YAK_MAGIC "YAK\2" -#ifndef kroundup32 -#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x)) -#endif - #define yak_ch_eq(a, b) ((a)>>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(static klib_unused, yak_ht_t, yak_ht, uint64_t, yak_ch_hash, yak_ch_eq) diff --git a/anchor.cpp b/anchor.cpp index f313ff5..b84a685 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -17,6 +17,9 @@ 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; @@ -40,7 +43,7 @@ void ha_abuf_destroy(ha_abuf_t *ab) free(ab->seed); free(ab->a); free(ab->mz.a); free(ab); } -void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double band_width_threshold, int keep_whole_chain) +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; @@ -106,9 +109,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg } cl->length = ab->n_a; - calculate_overlap_region_by_chaining(cl, overlap_list, rid, ucr->length, &R_INF, band_width_threshold, keep_whole_chain); - - ks_introsort_or_xs(overlap_list->length, overlap_list->list); + 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) { @@ -120,4 +121,33 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg } } #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/htab.h b/htab.h index 332b35a..ba155f8 100644 --- a/htab.h +++ b/htab.h @@ -81,6 +81,10 @@ static inline uint64_t yak_hash_long(uint64_t x[4]) #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