diff --git a/CommandLines.cpp b/CommandLines.cpp index a5be845..6ae58c5 100644 --- a/CommandLines.cpp +++ b/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); diff --git a/CommandLines.h b/CommandLines.h index bf59be7..c327588 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#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; diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 96961de..4b105c4 100644 --- a/Hash_Table.cpp +++ b/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 diff --git a/Hash_Table.h b/Hash_Table.h index d7820ae..e9bebe2 100644 --- a/Hash_Table.h +++ b/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; diff --git a/anchor.cpp b/anchor.cpp index a46a39b..a3c7dc8 100644 --- a/anchor.cpp +++ b/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;