r226b: prefer unique anchors

CommandLines.cpp

@@ -75,6 +75,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
 	asm_opt->mz_win = 51;
 	asm_opt->bf_shift = 37;
 	asm_opt->high_factor = 5.0f;
+	asm_opt->hom_cov = 20;
 	asm_opt->max_n_chain = 100;
     asm_opt->k_mer_min_freq = 3;
     asm_opt->k_mer_max_freq = 66;
@@ -117,6 +118,7 @@ void ha_opt_reset_to_round(hifiasm_opt_t* asm_opt, int round)
 void ha_opt_update_cov(hifiasm_opt_t *opt, int hom_cov)
 {
 	int max_n_chain = (int)(hom_cov * opt->high_factor + .499);
+	opt->hom_cov = hom_cov;
 	if (opt->max_n_chain < max_n_chain)
 		opt->max_n_chain = max_n_chain;
 	fprintf(stderr, "[M::%s] updated max_n_chain to %d\n", __func__, opt->max_n_chain);

CommandLines.h

@@ -3,7 +3,7 @@
 
 #include <pthread.h>
 
-#define HA_VERSION "0.3.0-dirty-r226"
+#define HA_VERSION "0.3.0-dirty-r226b"
 
 #define VERBOSE 0
 
@@ -27,6 +27,7 @@ typedef struct {
 	int mz_win;
 	int bf_shift;
 	float high_factor; // coverage cutoff set to high_factor*hom_cov
+	int hom_cov;
 	int max_n_chain; // fall-back max number of chains to consider
     int k_mer_min_freq;
     int k_mer_max_freq;

Hash_Table.cpp

@@ -385,6 +385,7 @@ int32_t ha_chain_lis_core(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min
 		dg = dx < dy? dx : dy;
 		if (dd > THRESHOLD_MAX_SIZE && dd > dg * bw_thres) break;
 		score = dg < min_sc? dg : min_sc;
+		if (!a[j1].good) score >>= 1;
 		gap_rate = (double)tot_indel / tot_len;
 		score -= (int)(gap_rate * score * bw_pen);
 		dp->score[i] = dp->score[i-1] + score;
@@ -459,6 +460,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul
             ///min distance
             distance_min = distance_pos < distance_self_pos? distance_pos:distance_self_pos;
             score = distance_min < min_score? distance_min : min_score;
+			if (!a[j].good) score >>= 1;
 
             gap_rate = (double)((double)(total_indels)/(double)(total_self_length));
             ///if the gap rate > 0.06, score will be negative

Hash_Table.h

@@ -92,7 +92,7 @@ typedef struct
     uint32_t y_pos_strand;
 
     uint32_t overlapLen;
-    uint32_t shared_seed;
+    int32_t shared_seed;
     uint32_t align_length;
     uint8_t is_match;
     uint8_t without_large_indel;
@@ -117,7 +117,7 @@ typedef struct
 
 typedef struct
 {
-	uint32_t readID:31, strand:1;
+	uint32_t readID:30, strand:1, good:1;
 	uint32_t offset, self_offset;
 } k_mer_hit;
 

anchor.cpp

@@ -3,9 +3,11 @@
 #include "ksort.h"
 #include "Hash_Table.h"
 
+#define HA_KMER_GOOD_RATIO 0.333
+
 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 self_off:31, good:1;
 	uint32_t other_off;
 } anchor1_t;
 
@@ -21,7 +23,7 @@ KSORT_INIT(or_xs, overlap_region, oreg_xs_lt)
 KSORT_INIT(or_ss, overlap_region, oreg_ss_lt)
 
 typedef struct {
-	int n;
+	int n, good;
 	const ha_idxpos_t *a;
 } seed1_t;
 
@@ -54,6 +56,8 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	extern ha_pt_t *ha_idx;
 	uint32_t i;
 	uint64_t k, l;
+	double low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO;
+	double high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO);
 
 	// prepare
     clear_Candidates_list(cl);
@@ -68,8 +72,11 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 		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;
+		int n;
+		ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &n);
+		ab->seed[i].n = n;
+		ab->seed[i].good = (n > low_occ && n < high_occ);
+		ab->n_a += n;
 	}
 	if (ab->n_a > ab->m_a) {
 		ab->m_a = ab->n_a;
@@ -86,6 +93,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 			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->good = s->good;
 			an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off;
 		}
 	}
@@ -111,6 +119,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 		p->strand = ab->a[k].srt >> 32 & 1;
 		p->offset = ab->a[k].other_off;
 		p->self_offset = ab->a[k].self_off;
+		p->good = ab->a[k].good;
 	}
 	cl->length = ab->n_a;