diff --git a/Assembly.cpp b/Assembly.cpp index 8edfe4a..48ae9c5 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -12,7 +12,7 @@ #include "kthread.h" void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, -int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf); +int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar); void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list); All_reads R_INF; @@ -457,6 +457,7 @@ typedef struct { kvec_t_u32_warp b_buf; kvec_t_u64_warp r_buf; kvec_t_u8_warp k_flag; + overlap_region tmp_region; } ha_ovec_buf_t; ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov) @@ -469,6 +470,7 @@ ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov) init_Candidates_list(&b->clist); init_overlap_region_alloc(&b->olist); init_overlap_region_alloc(&b->olist_hp); + init_fake_cigar(&(b->tmp_region.f_cigar)); kv_init(b->b_buf.a); kv_init(b->r_buf.a); kv_init(b->k_flag.a); @@ -492,6 +494,7 @@ void ha_ovec_destroy(ha_ovec_buf_t *b) destory_overlap_region_alloc(&b->olist); destory_overlap_region_alloc(&b->olist_hp); ha_abuf_destroy(b->ab); + destory_fake_cigar(&(b->tmp_region.f_cigar)); kv_destroy(b->b_buf.a); kv_destroy(b->r_buf.a); kv_destroy(b->k_flag.a); @@ -558,7 +561,7 @@ static void worker_ovec(void *data, long i, int tid) int fully_cov, abnormal; ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, - 0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &(R_INF.paf[i]), &(R_INF.reverse_paf[i])); + 0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region)); clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2); @@ -581,10 +584,10 @@ static void worker_ovec(void *data, long i, int tid) R_INF.paf[i].is_abnormal = abnormal; R_INF.trio_flag[i] = AMBIGU; - if(ha_idx_hp == NULL) - { - R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &b->b_buf, &b->r_buf, &(b->k_flag), RESEED_HP_RATE, NULL); - } + // if(ha_idx_hp == NULL) + // { + // R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &(b->k_flag), RESEED_HP_RATE, Get_READ_LENGTH(R_INF, i), NULL); + // } if (R_INF.trio_flag[i] != AMBIGU || b->save_ov) { int is_rev = (asm_opt.number_of_round % 2 == 0); @@ -615,7 +618,7 @@ static void worker_ovec_related_reads(void *data, long i, int tid) int fully_cov, abnormal; ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, - 0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &(R_INF.paf[i]), &(R_INF.reverse_paf[i])); + 0.02, asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region)); clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2); @@ -679,10 +682,10 @@ static void worker_ovec_related_reads(void *data, long i, int tid) } R_INF.trio_flag[i] = AMBIGU; - if(ha_idx_hp == NULL) - { - R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &b->b_buf, &b->r_buf, &(b->k_flag), RESEED_HP_RATE, R_INF_FLAG.fp); - } + // if(ha_idx_hp == NULL) + // { + // R_INF.trio_flag[i] += collect_hp_regions(&b->olist, &R_INF, &(b->k_flag), RESEED_HP_RATE, Get_READ_LENGTH(R_INF, i), R_INF_FLAG.fp); + // } fprintf(R_INF_FLAG.fp, "R_INF.trio_flag[%ld]: %u\n", i, R_INF.trio_flag[i]); @@ -891,7 +894,7 @@ void ha_overlap_and_correct(int round) free(e[i].second_round_read); } free(e); - debug_print_pob_regions(); + ///debug_print_pob_regions(); } @@ -1291,7 +1294,7 @@ static void worker_ov_final(void *data, long i, int tid) //get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0); ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, 0.001, - asm_opt.max_n_chain, 0, &(b->k_flag), &(R_INF.paf[i]), &(R_INF.reverse_paf[i])); + asm_opt.max_n_chain, 0, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region)); overlap_region_sort_y_id(b->olist.list, b->olist.length); ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length); @@ -1357,7 +1360,7 @@ static void worker_ov_final_high_het(void *data, long i, int tid) ha_ovec_buf_t *b = ((ha_ovec_buf_t**)data)[tid]; ha_get_candidates_interface(b->ab, i, &b->self_read, &b->olist, &b->olist_hp, &b->clist, HIGH_HET_ERROR_RATE, - asm_opt.max_n_chain, 1, &(b->k_flag), &(R_INF.paf[i]), &(R_INF.reverse_paf[i])); + asm_opt.max_n_chain, 1, &(b->k_flag), &b->r_buf, &(R_INF.paf[i]), &(R_INF.reverse_paf[i]), &(b->tmp_region)); overlap_region_sort_y_id(b->olist.list, b->olist.length); ma_hit_sort_tn(R_INF.paf[i].buffer, R_INF.paf[i].length); diff --git a/Correct.cpp b/Correct.cpp index 0389134..86d2285 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -12,6 +12,8 @@ #define generic_key(x) (x) KRADIX_SORT_INIT(b32, uint32_t, generic_key, 4) +int ha_ov_type(const overlap_region *r, uint32_t len); + void clear_Round2_alignment(Round2_alignment* h) { @@ -7268,12 +7270,15 @@ kvec_t_u32_warp* b, kvec_t_u64_warp* r, int min_dp, int min_len) } } -int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u32_warp* b, kvec_t_u64_warp* r, kvec_t_u8_warp* k_flag, float hp_rate, FILE* fp) +int collect_hp_regions_back(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u32_warp* b, kvec_t_u64_warp* r, kvec_t_u8_warp* k_flag, float hp_rate, FILE* fp) { - int i, k, qs, qe, ava_k_mer = 0, hp_k_mer = 0; - int min_dp = RESEED_DP; + int i, k, qs, qe, ava_k_mer = 0, hp_k_mer = 0, min_dp; + // min_dp = RESEED_DP; + // if(asm_opt.hom_cov > 0) min_dp = asm_opt.hom_cov * RESEED_PEAK_RATE; + // if(min_dp > RESEED_DP) min_dp = RESEED_DP; + min_dp = RESEED_DP; if(asm_opt.hom_cov > 0) min_dp = asm_opt.hom_cov * RESEED_PEAK_RATE; - if(min_dp > RESEED_DP) min_dp = RESEED_DP; + if(asm_opt.het_cov > 0) min_dp = asm_opt.het_cov * RESEED_PEAK_RATE; collect_no_cov_regions(olist, R_INF, b, r, min_dp, RESEED_LEN); for (i = 0; i < (int)r->a.n; i++) @@ -7291,21 +7296,86 @@ int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u32 if(fp) fprintf(fp, "qs: %d, qe: %d, ava_k_mer: %d, hp_k_mer: %d\n", qs, qe, ava_k_mer, hp_k_mer); } - if(fp) fprintf(fp, "ava_k_mer: %d, hp_k_mer: %d, hp_rate: %f\n", ava_k_mer, hp_k_mer, hp_rate); + if(fp) fprintf(fp, "ava_k_mer: %d, hp_k_mer: %d, hp_rate: %f, min_dp: %d, a.n: %d\n", ava_k_mer, hp_k_mer, hp_rate, min_dp, (int)r->a.n); - if(fp) - { - for (k = 0; k < (int)k_flag->a.n; k++) - { - if(k_flag->a.a[k] > 0) fprintf(fp, "(%d) %u\n", k, k_flag->a.a[k]); - } - } + // if(fp) + // { + // for (k = 0; k < (int)k_flag->a.n; k++) + // { + // if(k_flag->a.a[k] > 0) fprintf(fp, "(%d) %u\n", k, k_flag->a.a[k]); + // } + // } if(hp_k_mer > ava_k_mer*hp_rate) return 1; ///must use '>' instead of '>=' r->a.n = 0; return 0; } + +int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u8_warp* k_flag, +float hp_rate, int rlen, FILE* fp) +{ + int i, ava_k_mer = 0, hp_k_mer = 0, vLen, min_dp; + int32_t w, n[4]; + n[0] = n[1] = n[2] = n[3] = 0; + min_dp = RESEED_DP; + if(asm_opt.hom_cov > 0) min_dp = asm_opt.hom_cov * RESEED_PEAK_RATE; + if(min_dp > RESEED_DP) min_dp = RESEED_DP; + overlap_region* ov = NULL; + + for (i = 0; i < (long long)olist->length; i++) + { + ov = &(olist->list[i]); + if (ov->is_match != 1 && ov->is_match != 2) continue; + + w = ha_ov_type(ov, rlen); + ++n[w]; + } + + if(fp) fprintf(fp, "n[0]: %d, n[1]: %d, n[2]: %d, n[3]: %d\n", n[0], n[1], n[2], n[3]); + + // n[0] += n[2]; + // n[1] += n[2]; + + if(n[0] < min_dp) + { + ava_k_mer = hp_k_mer = 0; + vLen = MIN(k_flag->a.n, RESEED_LEN); + for (i = 0; i < vLen; i++) + { + if(k_flag->a.a[i] > 1) ava_k_mer++; + if(k_flag->a.a[i] > 2) hp_k_mer++; + } + if(hp_k_mer > ava_k_mer*hp_rate) return 1; + } + + + if(n[1] < min_dp) + { + ava_k_mer = hp_k_mer = 0; + vLen = MIN(k_flag->a.n, RESEED_LEN); + for (i = k_flag->a.n - vLen; i < (int)k_flag->a.n; i++) + { + if(k_flag->a.a[i] > 1) ava_k_mer++; + if(k_flag->a.a[i] > 2) hp_k_mer++; + } + if(hp_k_mer > ava_k_mer*hp_rate) return 1; + } + + + if(fp) fprintf(fp, "ava_k_mer: %d, hp_k_mer: %d, hp_rate: %f, min_dp: %d\n", ava_k_mer, hp_k_mer, hp_rate, min_dp); + + // if(fp) + // { + // for (k = 0; k < (int)k_flag->a.n; k++) + // { + // if(k_flag->a.a[k] > 0) fprintf(fp, "(%d) %u\n", k, k_flag->a.a[k]); + // } + // } + + return 0; +} + void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, Graph* g, Graph* DAGCon, Cigar_record* current_cigar, diff --git a/Correct.h b/Correct.h index 8bb4735..f51a5df 100644 --- a/Correct.h +++ b/Correct.h @@ -1173,7 +1173,7 @@ void correct_overlap_high_het(overlap_region_alloc* overlap_list, All_reads* R_I UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read); long long get_affine_gap_score(overlap_region* ovc, UC_Read* g_read, UC_Read* overlap_read, uint8_t* x_num, uint8_t* y_num, uint64_t EstimateXOlen, uint64_t EstimateYOlen); -int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u32_warp* b, kvec_t_u64_warp* r, kvec_t_u8_warp* k_flag, float hp_rate, FILE* fp); +int collect_hp_regions(overlap_region_alloc* olist, All_reads* R_INF, kvec_t_u8_warp* k_flag, float hp_rate, int rlen, FILE* fp); #define FORWARD_KSW 0 #define BACKWARD_KSW 1 diff --git a/Hash_Table.cpp b/Hash_Table.cpp index fdda1e0..f1eaeca 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -305,6 +305,9 @@ void debug_chain(k_mer_hit* a, long long a_n, Chain_Data* dp) } } + + + long long get_chainLen(long long x_beg, long long x_end, long long xLen, long long y_beg, long long y_end, long long yLen) { @@ -337,6 +340,50 @@ long long get_chainLen(long long x_beg, long long x_end, long long xLen, return x_end - x_beg + 1; } + +void debug_chain_single_site(k_mer_hit* a, long long a_n, Chain_Data* dp, int x_readLen, int y_readLen, int s_index) +{ + long long j, current_j = s_index; + long long selfLen = 0, indels = 0; + long long distance_self_pos, distance_pos, distance_gap; + + j = s_index; + while (j >= 0) + { + current_j = j; + + j = dp->pre[j]; + + if(j != -1) + { + distance_self_pos = a[current_j].self_offset - a[j].self_offset; + distance_pos = a[current_j].offset - a[j].offset; + distance_gap = distance_pos > distance_self_pos? distance_pos - distance_self_pos : distance_self_pos - distance_pos; + + indels += distance_gap; + selfLen += distance_self_pos; + } + fprintf(stderr, "j: %lld, score: %lld, occ: %d, pre_j: %lld\n", + current_j, (long long)dp->score[current_j], dp->occ[current_j], j); + } + + fprintf(stderr, "s_self_offset: %u, s_offset: %u, e_self_offset: %u, e_offset: %u, ovlp length: %lld, x_readLen: %d, y_readLen: %d\n", + a[s_index].self_offset, a[s_index].offset, a[current_j].self_offset, a[current_j].offset, + get_chainLen(a[s_index].self_offset, a[current_j].self_offset, x_readLen, + a[s_index].offset, a[current_j].offset, y_readLen), x_readLen, y_readLen); + + if(indels != dp->indels[s_index]) + { + fprintf(stderr, "indels: %lld, dp->indels[i]: %ld\n", indels, (long)dp->indels[s_index]); + } + + if(selfLen != dp->self_length[s_index]) + { + fprintf(stderr, "selfLen: %lld, dp->self_length[i]: %ld\n", selfLen, (long)dp->self_length[s_index]); + } + fprintf(stderr,"\n"); +} + int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc, double bw_thres) { int32_t i, tot_indel = 0, tot_len = 0; @@ -348,7 +395,7 @@ int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc if (i < n_a) return -1; bw_pen = 1.0 / bw_thres; dp->score[0] = a[0].good? min_sc : min_sc>>1; - dp->pre[0] = -1, dp->indels[0] = 0, dp->self_length[0] = 0; + dp->pre[0] = -1, dp->indels[0] = 0, dp->self_length[0] = 0, dp->occ[0] = 1; for (i = 1; i < n_a; ++i) { int32_t score, dg; int32_t dx = (int32_t)a[i].offset - (int32_t)a[i-1].offset; @@ -368,6 +415,7 @@ int32_t ha_chain_check(k_mer_hit *a, int32_t n_a, Chain_Data *dp, int32_t min_sc dp->pre[i] = i - 1; dp->indels[i] = tot_indel; dp->self_length[i] = tot_len; + dp->occ[i] = i + 1; } if (i < n_a) return -1; return n_a; @@ -391,16 +439,7 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul resize_Chain_Data(dp, a_n); ret = ha_chain_check(a, a_n, dp, min_score, band_width_threshold); - /***************************************debug**************************************/ - if(a_n > 0 && Get_NAME_LENGTH((R_INF),a[0].readID)==strlen("m64062_190803_042216/128778853/ccs")) - { - if (memcmp("m64062_190803_042216/128778853/ccs", Get_NAME((R_INF), a[0].readID), - Get_NAME_LENGTH((R_INF), a[0].readID)) == 0) - { - fprintf(stderr, "ret: %d\n", ret); - } - } - /***************************************debug**************************************/ + if (ret > 0) { a_n = ret; goto skip_dp; @@ -445,7 +484,8 @@ 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; + ///if (!a[j].good) score = (score >> 1) + (score & 1); + 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 @@ -454,7 +494,7 @@ 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) { + if (score > max_score) {///must use > instead of >= max_score = score; max_j = j; max_indels = total_indels; @@ -476,10 +516,12 @@ void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* resul dp->pre[i] = max_j; dp->indels[i] = max_indels; dp->self_length[i] = max_self_length; + dp->occ[i] = 1; + if(max_j != -1) dp->occ[i] = dp->occ[max_j] + 1; } ///debug_chain(a, a_n, dp); - + skip_dp: max_score = -1; @@ -510,17 +552,6 @@ skip_dp: } - /***************************************debug**************************************/ - if(a_n > 0 && Get_NAME_LENGTH((R_INF),a[0].readID)==strlen("m64062_190803_042216/128778853/ccs")) - { - if (memcmp("m64062_190803_042216/128778853/ccs", Get_NAME((R_INF), a[0].readID), - Get_NAME_LENGTH((R_INF), a[0].readID)) == 0) - { - fprintf(stderr, "max_i: %lld\n", max_i); - } - } - /***************************************debug**************************************/ - clear_fake_cigar(&(result->f_cigar)); ///note a has been sorted by offset, that means has been sorted by query offset i = max_i; @@ -583,7 +614,7 @@ skip_dp: } } -void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, +void calculate_overlap_region_by_chaining_back(Candidates_list* candidates, overlap_region_alloc* overlap_list, uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end) { overlap_region tmp_region; @@ -648,6 +679,71 @@ void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_r destory_fake_cigar(&(tmp_region.f_cigar)); } + +void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, kvec_t_u64_warp* chain_idx, + uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end, overlap_region* f_cigar) +{ + long long i = 0; + uint64_t current_ID; + uint64_t current_stand; + + if (candidates->length == 0) + { + return; + } + + long long sub_region_beg; + long long sub_region_end; + + clear_fake_cigar(&((*f_cigar).f_cigar)); + + i = 0; + while (i < candidates->length) + { + chain_idx->a.n = 0; + current_ID = candidates->list[i].readID; + current_stand = candidates->list[i].strand; + + ///reference read + (*f_cigar).x_id = readID; + (*f_cigar).x_pos_strand = current_stand; + ///query read + (*f_cigar).y_id = current_ID; + ///here the strand of query is always 0 + (*f_cigar).y_pos_strand = 0; + + sub_region_beg = i; + sub_region_end = i; + i++; + + while (i < candidates->length + && + current_ID == candidates->list[i].readID + && + current_stand == candidates->list[i].strand) + { + sub_region_end = i; + i++; + } + + if ((*f_cigar).x_id == (*f_cigar).y_id) + { + continue; + } + + chain_DP(candidates->list + sub_region_beg, + sub_region_end - sub_region_beg + 1, &(candidates->chainDP), f_cigar, band_width_threshold, + 25, Get_READ_LENGTH((*R_INF), (*f_cigar).x_id), Get_READ_LENGTH((*R_INF), (*f_cigar).y_id)); + + ///if (tmp_region.x_id != tmp_region.y_id && tmp_region.shared_seed > 1) + if ((*f_cigar).x_id != (*f_cigar).y_id) + { + append_inexact_overlap_region_alloc(overlap_list, f_cigar, R_INF, add_beg_end); + } + } +} + + void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error, int extra_begin, int extra_end, int error_threshold) { @@ -725,6 +821,7 @@ void destory_Chain_Data(Chain_Data* x) free(x->pre); free(x->indels); free(x->self_length); + free(x->occ); free(x->tmp); } @@ -737,6 +834,7 @@ void resize_Chain_Data(Chain_Data* x, long long size) REALLOC(x->pre, x->size); REALLOC(x->indels, x->size); REALLOC(x->self_length, x->size); + REALLOC(x->occ, x->size); REALLOC(x->tmp, x->size); } } diff --git a/Hash_Table.h b/Hash_Table.h index 157adc5..e0742a4 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -119,6 +119,7 @@ typedef struct { int64_t *pre; int32_t *indels; int32_t *self_length; + int32_t *occ; int64_t *tmp; // MUST BE 64-bit integer int64_t length; int64_t size; @@ -144,8 +145,8 @@ int extra_begin, int extra_end, int error_threshold); void overlap_region_sort_y_id(overlap_region *a, long long n); -void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, -uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end); +void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, kvec_t_u64_warp* chain_idx, +uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end, overlap_region* f_cigar); void init_fake_cigar(Fake_Cigar* x); void destory_fake_cigar(Fake_Cigar* x); diff --git a/Process_Read.h b/Process_Read.h index fc29101..d259bd4 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -22,7 +22,8 @@ ///#define Get_READ(R_INF, ID) R_INF.read + (R_INF.index[ID]>>2) + ID #define Get_READ(R_INF, ID) (R_INF).read_sperate[(ID)] #define Get_NAME(R_INF, ID) ((R_INF).name + (R_INF).name_index[(ID)]) - +#define CHECK_BY_NAME(R_INF, NAME, ID) (Get_NAME_LENGTH((R_INF),(ID))==strlen((NAME)) && \ + memcmp((NAME), Get_NAME((R_INF), (ID)), Get_NAME_LENGTH((R_INF),(ID))) == 0) extern uint8_t seq_nt6_table[256]; extern char bit_t_seq_table[256][4]; diff --git a/anchor.cpp b/anchor.cpp index a1f9ed3..7e0bd55 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -50,7 +50,7 @@ uint64_t ha_abuf_mem(const ha_abuf_t *ab) return ab->m_a * sizeof(anchor1_t) + ab->mz.m * (sizeof(ha_mz1_t) + sizeof(seed1_t)) + sizeof(ha_abuf_t); } -static int ha_ov_type(const overlap_region *r, uint32_t len) +int ha_ov_type(const overlap_region *r, uint32_t len) { if (r->x_pos_s == 0 && r->x_pos_e == len - 1) return 2; // contained in a longer read else if (r->x_pos_s > 0 && r->x_pos_e < len - 1) return 3; // containing a shorter read @@ -58,7 +58,7 @@ static int ha_ov_type(const overlap_region *r, uint32_t len) } 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, kvec_t_u8_warp* k_flag, -void *ha_flt_tab, ha_pt_t *ha_idx) +kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar) { uint32_t i, rlen; uint64_t k, l; @@ -130,22 +130,10 @@ void *ha_flt_tab, ha_pt_t *ha_idx) p->offset = ab->a[k].other_off; p->self_offset = ab->a[k].self_off; p->good = ab->a[k].good; - - /***************************************debug**************************************/ - if(Get_NAME_LENGTH((R_INF),p->readID)==strlen("m64062_190803_042216/128778853/ccs")) - { - if (memcmp("m64062_190803_042216/128778853/ccs", Get_NAME((R_INF), p->readID), - Get_NAME_LENGTH((R_INF), p->readID)) == 0) - { - fprintf(stderr, "(%lu) readID: %u, strand: %u, offset: %u, self_offset: %u\n", - k, p->readID, p->strand, p->offset, p->self_offset); - } - } - /***************************************debug**************************************/ } cl->length = ab->n_a; - calculate_overlap_region_by_chaining(cl, overlap_list, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain); + calculate_overlap_region_by_chaining(cl, overlap_list, chain_idx, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain, f_cigar); #if 0 if (overlap_list->length > 0) { @@ -217,14 +205,14 @@ void lable_matched_ovlp(overlap_region_alloc* overlap_list, ma_hit_t_alloc* paf) void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, -int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf) +int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, ma_hit_t_alloc* paf, ma_hit_t_alloc* rev_paf, overlap_region* f_cigar) { extern void *ha_flt_tab; extern ha_pt_t *ha_idx; extern void *ha_flt_tab_hp; extern ha_pt_t *ha_idx_hp; - ha_get_new_candidates(ab, rid, ucr, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, ha_flt_tab, ha_idx); + ha_get_new_candidates(ab, rid, ucr, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar); if(ha_idx_hp) { @@ -254,7 +242,7 @@ int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, ma_hit_t_alloc* p overlap_list->length = k; - ha_get_new_candidates(ab, rid, ucr, overlap_list_hp, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, ha_flt_tab_hp, ha_idx_hp); + ha_get_new_candidates(ab, rid, ucr, overlap_list_hp, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab_hp, ha_idx_hp, f_cigar); if(overlap_list->length + overlap_list_hp->length > overlap_list->size) {