r326: select high-occ k-mers (experimental)

CommandLines.cpp

@@ -102,7 +102,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->k_mer_length = 51;
 	asm_opt->mz_win = 51;
 	asm_opt->bf_shift = 37;
-	asm_opt->max_kmer_cnt = 4000;
+	asm_opt->max_kmer_cnt = 2000;
 	asm_opt->high_factor = 5.0;
 	asm_opt->max_ov_diff_ec = 0.04;
 	asm_opt->max_ov_diff_final = 0.03;

CommandLines.h

@@ -3,7 +3,7 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.13-r308"
+#define HA_VERSION "0.13-r326-dirty"
 
 #define VERBOSE 0
 

htab.cpp

@@ -189,7 +189,7 @@ static int ha_ct_insert_list(ha_ct_t *h, int create_new, int n, const uint64_t *
 		///so low 12 bits are not useful
 		uint64_t x = a[j] >> h->pre;
 		khint_t k;
-		if ((a[j]&mask) != (a[0]&mask)) continue;
+		assert((a[j]&mask) == (a[0]&mask));
 		if (create_new) {
 			///for 0-th counting, g->b = NULL
 			if (g->b)
@@ -365,9 +365,9 @@ int ha_pt_insert_list(ha_pt_t *h, int n, const ha_mz1_t *a)
 		khint_t k;
 		int n;
 		ha_idxpos_t *p;
-		if ((a[j].x&mask) != (a[0].x&mask)) continue;
+		assert((a[j].x&mask) == (a[0].x&mask));
 		k = yak_pt_get(g->h, x<<YAK_COUNTER_BITS);
-		if (k == kh_end(g->h)) continue;
+		if (k == kh_end(g->h)) continue; // TODO: understand why we sometimes come here
 		n = kh_key(g->h, k) & YAK_MAX_COUNT;
 		assert(n < YAK_MAX_COUNT);
 		p = &g->a[kh_val(g->h, k) + n];
@@ -518,7 +518,7 @@ 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; ///number of total reads
+	uint64_t n_mz, n_seq; ///number of total reads
 	kseq_t *ks;
 	UC_Read ucr;
 	ha_ct_t *ct;
@@ -573,8 +573,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 		s->n_seq0 = p->n_seq;
 		if (p->rs_in && (p->flag & HAF_RS_READ)) {
 			while (p->n_seq < p->rs_in->total_reads) {
-				if((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU)
+				if ((p->flag & HAF_SKIP_READ) && p->rs_in->trio_flag[p->n_seq] != AMBIGU) {
-				{
 					++p->n_seq;
 					continue;
 				}
@@ -685,6 +684,7 @@ static void *worker_count(void *data, int step, void *in) // callback for kt_pip
 						ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);
 			}
 			for (i = 0; i < s->n_seq; ++i) {
+				p->n_mz += s->mz[i].n;
 				free(s->mz[i].a);
 				if (!p->is_store) free(s->seq[i]);
 			}
@@ -763,6 +763,7 @@ static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt
 		gzclose(fp);
 	}
 	*n_seq = pl.n_seq;
+	if (pl.opt->w > 1) fprintf(stderr, "[M::%s] collected %ld minimizers\n", __func__, (long)pl.n_mz);
 	return pl.ct;
 }
 

sketch.cpp

@@ -25,6 +25,44 @@ static inline int tq_shift(tiny_queue_t *q)
 	return x;
 }
 
+static inline int mzcmp(const ha_mz1_t *a, const ha_mz1_t *b)
+{
+	return a->rid < b->rid? -1 : a->rid > b->rid? 1 : ((a->x > b->x) - (a->x < b->x));
+}
+
+static void select_mz(ha_mz1_v *p, int len)
+{
+	static const ha_mz1_t dummy = { UINT64_MAX, (1<<28) - 1, 0, 0 };
+	int32_t i, last0 = -1, n = (int32_t)p->n, m = 0;
+	if (n == 0 || n == 1) return;
+	for (i = 0; i < n; ++i)
+		if (p->a[i].rid != 0) ++m;
+	if (m == 0) return; // no high-frequency k-mers; do nothing
+	for (i = 0; i <= n; ++i) {
+		if (i == n || p->a[i].rid == 0) {
+			if (i - last0 > 1) {
+				int32_t ps = last0 < 0? 0 : p->a[last0].pos;
+				int32_t pe = i == n? len : p->a[i].pos;
+				int32_t j, st = last0 + 1, en = i;
+				ha_mz1_t min1 = dummy, min2 = dummy;
+				int32_t min1_i = -1, min2_i = -1;
+				for (j = st; j < en; ++j) { // choose up to two minimum k-mers
+					if (mzcmp(&p->a[j], &min1) < 0) min2 = min1, min2_i = min1_i, min1 = p->a[j], min1_i = j;
+					else if (mzcmp(&p->a[j], &min2) < 0) min2 = p->a[j], min2_i = j;
+				}
+				if (min1_i >= 0 && p->a[min1_i].rid < pe - ps) p->a[min1_i].rid = 0;
+				if (min2_i >= 0 && p->a[min2_i].rid < pe - ps) p->a[min2_i].rid = 0;
+			}
+			last0 = i;
+		}
+	}
+	for (i = n = 0; i < (int32_t)p->n; ++i) // squeeze out filtered minimizers
+		if (p->a[i].rid == 0)
+			p->a[n++] = p->a[i];
+//	fprintf(stderr, "X\tn0=%d,n1=%d,m=%d\n", p->n, n, m);
+	p->n = n;
+}
+
 /**
  * Find symmetric (w,k)-minimizers on a DNA sequence
  *
@@ -43,7 +81,7 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 	 uint64_t rid:28, pos:27, rev:1, span:8;
 	 **/
 	extern void *ha_ct_table;
-	static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 };
+	static const ha_mz1_t dummy = { UINT64_MAX, (1<<28) - 1, 0, 0 };
 	uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0};
 	int i, j, l, buf_pos, min_pos, kmer_span = 0;
 	ha_mz1_t buf[256], min = dummy;
@@ -96,13 +134,12 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 			++l;
 			if (l >= k && kmer_span < 256) {
 				uint64_t y;
-				int32_t cnt = 0, filtered = 0;
+				int32_t cnt, filtered;
 				y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
-				if (hf != 0) cnt = ha_ft_cnt(hf, y);
+				cnt = hf? ha_ft_cnt(hf, y) : 0;
-				filtered = (cnt > 0);
+				filtered = (cnt >= 1<<28);
 				if (dbg_ct != NULL) kv_push(uint64_t, dbg_ct->a, ((((uint64_t)(query_ct_index(ha_ct_table, y))<<1)|filtered)<<32)|(uint64_t)(i));
-				if (cnt < 1<<28 && filtered == 0)
+				if (!filtered) info.x = y, info.rid = cnt, info.pos = i, info.rev = z, info.span = kmer_span;
-					info.x = y, info.rid = cnt, info.pos = i, info.rev = z, info.span = kmer_span;
 				if (k_flag != NULL) k_flag->a.a[i]++;
 				if (k_flag != NULL && filtered > 0) k_flag->a.a[i]++;
 			}
@@ -115,9 +152,9 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 		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]);
+				if (mzcmp(&min, &buf[j]) == 0 && 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 (mzcmp(&min, &buf[j]) == 0 && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
 		}
 		/**
 		 * There are three cases:
@@ -135,21 +172,22 @@ void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc,
 			///buf_pos == min_pos, means current minimizer has moved outside the window
 			///so for now we need to find a new minimizer at the current window (w k-mers)
 			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
+				if (mzcmp(&min, &buf[j]) >= 0) 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 (mzcmp(&min, &buf[j]) >= 0) 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]);
+					if (mzcmp(&min, &buf[j]) == 0 && 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 (mzcmp(&min, &buf[j]) == 0 && 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);
+	select_mz(p, len);
 	for (i = 0; i < (int)p->n; ++i) // populate .rid as this was keeping counts
 		p->a[i].rid = rid;
 }