From b7e5d1c4d3f9afe334ecbffce475cfcab7b4f31a Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 1 Apr 2020 22:06:14 -0400 Subject: [PATCH] r176: enabled the minimap2 chaining heuristic --- CommandLines.h | 2 +- Hash_Table.cpp | 74 ++++++++++++++++++++------------------------------ Hash_Table.h | 19 ++++++------- anchor.cpp | 10 +------ htab.h | 4 +++ 5 files changed, 44 insertions(+), 65 deletions(-) diff --git a/CommandLines.h b/CommandLines.h index 29e02d3..f3e7161 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.3.0-dirty-r175" +#define HA_VERSION "0.3.0-dirty-r176" #define VERBOSE 0 #define VERBOSE_GFA 0 diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 397d2fd..5873b32 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -356,14 +356,12 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp) if(indels != dp->indels[i]) { - fprintf(stderr, "indels: %lld, dp->indels[i]: %lld\n", - indels, dp->indels[i]); + fprintf(stderr, "indels: %lld, dp->indels[i]: %ld\n", indels, (long)dp->indels[i]); } if(selfLen != dp->self_length[i]) { - fprintf(stderr, "selfLen: %lld, dp->self_length[i]: %lld\n", - selfLen, dp->self_length[i]); + fprintf(stderr, "selfLen: %lld, dp->self_length[i]: %ld\n", selfLen, (long)dp->self_length[i]); } } @@ -418,6 +416,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul resize_Chain_Data(dp, a_n); // fill the score and backtrack arrays + for (i = 0; i < a_n; ++i) dp->tmp[i] = -1; for (i = 0; i < a_n; ++i) { pos = a[i].offset; @@ -461,26 +460,18 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul score += dp->score[j]; - ///find a new max score - if(score > max_score) - { - max_score = score; - max_j = j; - max_indels = total_indels; - max_self_length = total_self_length; - /****************************may have bugs********************************/ - n_skip = 0; - /****************************may have bugs********************************/ - }/****************************may have bugs********************************/ - else - { - n_skip++; - if(n_skip > max_skip) - { - break; - } - } - /****************************may have bugs********************************/ + ///find a new max score + if (score > max_score) { + max_score = score; + max_j = j; + max_indels = total_indels; + max_self_length = total_self_length; + if (n_skip > 0) --n_skip; + } else if (dp->tmp[j] == i) { + if (++n_skip > max_skip) + break; + } + if (dp->pre[j] >= 0) dp->tmp[dp->pre[j]] = i; } dp->score[i] = max_score; @@ -636,14 +627,10 @@ void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_r } - chain_DP(candidates->list + sub_region_beg, - sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold, - 50, Get_READ_LENGTH((*R_INF), tmp_region.x_id), Get_READ_LENGTH((*R_INF), tmp_region.y_id)); + chain_DP(candidates->list + sub_region_beg, + sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold, + 25, Get_READ_LENGTH((*R_INF), tmp_region.x_id), Get_READ_LENGTH((*R_INF), tmp_region.y_id)); - // chain_DP_back(candidates->list + sub_region_beg, - // sub_region_end - sub_region_beg + 1, &(candidates->chainDP), &tmp_region, band_width_threshold); - - ///自己和自己重叠的要排除 ///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1) if (tmp_region.x_id != tmp_region.y_id) { @@ -720,12 +707,7 @@ void test_single_list(Candidates_list* candidates, k_mer_pos* n_list, uint64_t n void init_Chain_Data(Chain_Data* x) { - x->length = 0; - x->size = 0; - x->score = NULL; - x->pre = NULL; - x->indels = NULL; - x->self_length = NULL; + memset(x, 0, sizeof(Chain_Data)); } void clear_Chain_Data(Chain_Data* x) @@ -739,18 +721,20 @@ void destory_Chain_Data(Chain_Data* x) free(x->pre); free(x->indels); free(x->self_length); + free(x->tmp); } void resize_Chain_Data(Chain_Data* x, long long size) { - if(size > x->size) - { - x->size = size; - x->score = (long long*)realloc(x->score, x->size*sizeof(long long)); - x->pre = (long long*)realloc(x->pre, x->size*sizeof(long long)); - x->indels = (long long*)realloc(x->indels, x->size*sizeof(long long)); - x->self_length = (long long*)realloc(x->self_length, x->size*sizeof(long long)); - } + if (size > x->size) { + x->size = size; + kroundup64(x->size); + REALLOC(x->score, x->size); + REALLOC(x->pre, x->size); + REALLOC(x->indels, x->size); + REALLOC(x->self_length, x->size); + REALLOC(x->tmp, x->size); + } } void init_Candidates_list(Candidates_list* l) diff --git a/Hash_Table.h b/Hash_Table.h index 8e570ee..e890c9a 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -113,8 +113,7 @@ typedef struct overlap_region* list; uint64_t size; uint64_t length; - ///uint64_t mapped_overlaps_length; - long long mapped_overlaps_length; + int64_t mapped_overlaps_length; } overlap_region_alloc; typedef struct @@ -123,14 +122,14 @@ typedef struct uint32_t offset, self_offset; } k_mer_hit; -typedef struct -{ - long long* score; - long long* pre; - long long* indels; - long long* self_length; - long long length; - long long size; +typedef struct { + int32_t *score; + int64_t *pre; + int32_t *indels; + int32_t *self_length; + int64_t *tmp; + int64_t length; + int64_t size; } Chain_Data; typedef struct diff --git a/anchor.cpp b/anchor.cpp index 3b1cbca..95dfbda 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -62,7 +62,7 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg } if (ab->n_a > ab->m_a) { ab->m_a = ab->n_a; - ab->m_a = ab->m_a > 16? ab->m_a + (ab->m_a>>1) : 16; + kroundup64(ab->m_a); REALLOC(ab->a, ab->m_a); } for (i = 0, k = 0; i < ab->mz.n; ++i) { @@ -104,12 +104,4 @@ 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); - #if 0 - fprintf(stderr, "B\t%ld\t%ld\n", (long)overlap_list->length, (long)overlap_list->mapped_overlaps_length); - for (int i = 0; i < (int)overlap_list->length; ++i) { - overlap_region *r = &overlap_list->list[i]; - fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%d\t%d\t%c\t%d\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand], - (int)r->y_id, (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed, r->is_match); - } - #endif } diff --git a/htab.h b/htab.h index 65bbd96..332b35a 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 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 + #ifndef klib_unused #if (defined __clang__ && __clang_major__ >= 3) || (defined __GNUC__ && __GNUC__ >= 3) #define klib_unused __attribute__ ((__unused__))