r181: limit the max number of chains

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(&current_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,

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;

CommandLines.h

@@ -3,7 +3,7 @@
 
 #include <pthread.h>
 
-#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;

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;
 }
 

Hash_Table.h

@@ -67,15 +67,15 @@ typedef struct
 typedef struct
 {
     window_list* buffer;
-    long long length;
+    int32_t length;
-    long long size;
+    int32_t size;
 } window_list_alloc;
 
 typedef struct
 {
     uint64_t* buffer;
-    uint64_t length;
+	uint32_t length;
-    uint64_t size;
+	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;
+    uint32_t w_list_size;
-    uint64_t w_list_length;
+    uint32_t w_list_length;
-    int8_t strong;
     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);
     }
 

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

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)

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);
+	calculate_overlap_region_by_chaining(cl, overlap_list, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain);
-
-	ks_introsort_or_xs(overlap_list->length, overlap_list->list);
 
 	#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);
 }

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