diff --git a/Correct.cpp b/Correct.cpp index d84acbc..1f9b9a9 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -539,6 +539,236 @@ char* r_string) } } + +void verify_ul_window(long long window_start, long long window_end, overlap_region_alloc* overlap_list,Correct_dumy* dumy, const ma_ug_t *uref, +char* r_string) +{ + long long i; + long long currentID; + long long x_start, y_start, o_len; + long long Window_Len = WINDOW + (THRESHOLD << 1); + char* x_string = NULL; + char* y_string = NULL; + long long x_end, x_len; + int end_site; + unsigned int error; + int groupLen = 0; + int return_sites[GROUP_SIZE]; + unsigned int return_sites_error[GROUP_SIZE]; + uint64_t overlapID[GROUP_SIZE]; + uint64_t y_startGroup[GROUP_SIZE]; + int y_extra_begin[GROUP_SIZE]; + int y_extra_end[GROUP_SIZE]; + int error_threshold[GROUP_SIZE]; + int extra_begin; + int extra_end; + + ///here are overlaps fully covered by WINDOW + for (i = 0; i < (long long)dumy->length; i++) + { + extra_begin = extra_end = 0; + ///if the window has been fully covered, the interval at x is [window_start, window_end] + x_len = WINDOW; + currentID = dumy->overlapID[i]; + x_start = window_start; + ///offset of y + y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s; + /****************************may have bugs********************************/ + y_start += y_start_offset(x_start, &(overlap_list->list[currentID].f_cigar)); + /****************************may have bugs********************************/ + + + if(!determine_overlap_region(THRESHOLD, y_start, overlap_list->list[currentID].y_id, Window_Len, R_INF, + &extra_begin, &extra_end, &y_start, &o_len)) + { + append_window_list(&overlap_list->list[currentID], window_start, window_end, + -1, -1, -1, -1, -1, -1); + continue; + } + + fill_subregion(dumy->overlap_region_group[groupLen], y_start, o_len, overlap_list->list[currentID].y_pos_strand, + R_INF, overlap_list->list[currentID].y_id, extra_begin, extra_end); + + y_extra_begin[groupLen] = extra_begin; + y_extra_end[groupLen] = extra_end; + overlapID[groupLen] = currentID; + y_startGroup[groupLen] = y_start; + error_threshold[groupLen] = THRESHOLD; + x_string = r_string + x_start; + groupLen++; + + + if (groupLen == GROUP_SIZE) + { + Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1], + dumy->overlap_region_group[2], dumy->overlap_region_group[3], Window_Len, x_string, WINDOW, + return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE); + groupLen = 0; + + + if (return_sites_error[0]!=(unsigned int)-1) + { + overlap_list->list[overlapID[0]].align_length += x_len; + + append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, + y_startGroup[0], y_startGroup[0] + return_sites[0], (int)return_sites_error[0], + y_extra_begin[0], y_extra_end[0], error_threshold[0]); + } + else + { + append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], -1, -1, + y_extra_begin[0], y_extra_end[0], error_threshold[0]); + } + + + if (return_sites_error[1]!=(unsigned int)-1) + { + overlap_list->list[overlapID[1]].align_length += x_len; + + append_window_list(&overlap_list->list[overlapID[1]], window_start, window_end, + y_startGroup[1], y_startGroup[1] + return_sites[1], (int)return_sites_error[1], + y_extra_begin[1], y_extra_end[1], error_threshold[1]); + } + else + { + append_window_list(&overlap_list->list[overlapID[1]], window_start, window_end, y_startGroup[1], -1, -1, + y_extra_begin[1], y_extra_end[1], error_threshold[1]); + } + + + if (return_sites_error[2]!=(unsigned int)-1) + { + overlap_list->list[overlapID[2]].align_length += x_len; + + append_window_list(&overlap_list->list[overlapID[2]], window_start, window_end, + y_startGroup[2], y_startGroup[2] + return_sites[2], (int)return_sites_error[2], + y_extra_begin[2], y_extra_end[2], error_threshold[2]); + } + else + { + append_window_list(&overlap_list->list[overlapID[2]], window_start, window_end, y_startGroup[2], -1, -1, + y_extra_begin[2], y_extra_end[2], error_threshold[2]); + } + + + if (return_sites_error[3]!=(unsigned int)-1) + { + overlap_list->list[overlapID[3]].align_length += x_len; + + append_window_list(&overlap_list->list[overlapID[3]], window_start, window_end, + y_startGroup[3], y_startGroup[3] + return_sites[3], (int)return_sites_error[3], + y_extra_begin[3], y_extra_end[3], error_threshold[3]); + } + else + { + append_window_list(&overlap_list->list[overlapID[3]], window_start, window_end, y_startGroup[3], -1, -1, + y_extra_begin[3], y_extra_end[3], error_threshold[3]); + } + } + } + + if (groupLen == 1) + { + end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], Window_Len, x_string, WINDOW, THRESHOLD, &error); + + if (error!=(unsigned int)-1) + { + overlap_list->list[overlapID[0]].align_length += x_len; + + append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, + y_startGroup[0], y_startGroup[0] + end_site, (int)error, + y_extra_begin[0], y_extra_end[0], error_threshold[0]); + } + else + { + append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], -1, -1, + y_extra_begin[0], y_extra_end[0], error_threshold[0]); + } + } + else if (groupLen > 1) + { + Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1], + dumy->overlap_region_group[2], dumy->overlap_region_group[3], Window_Len, x_string, WINDOW, + return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE); + + for (i = 0; i < groupLen; i++) + { + if (return_sites_error[i]!=(unsigned int)-1) + { + overlap_list->list[overlapID[i]].align_length += x_len; + append_window_list(&overlap_list->list[overlapID[i]], window_start, window_end, + y_startGroup[i], y_startGroup[i] + return_sites[i], (int)return_sites_error[i], + y_extra_begin[i], y_extra_end[i], error_threshold[i]); + } + else + { + append_window_list(&overlap_list->list[overlapID[i]], window_start, window_end, y_startGroup[i], -1, -1, + y_extra_begin[i], y_extra_end[i], error_threshold[i]); + } + + } + + groupLen = 0; + } + + long long reverse_i = dumy->size - 1; + int threshold; + + ///here are overlaps partially covered by WINDOW + for (i = 0; i < (long long)dumy->lengthNT; i++) + { + extra_begin = extra_end = 0; + currentID = dumy->overlapID[reverse_i--]; + x_start = MAX(window_start, (long long)overlap_list->list[currentID].x_pos_s); + x_end = MIN(window_end, (long long)overlap_list->list[currentID].x_pos_e); + + ///overlap length between [window_start, window_end] + x_len = x_end - x_start + 1; + threshold = x_len * asm_opt.max_ov_diff_ec; + /****************************may have bugs********************************/ + threshold = Adjust_Threshold(threshold, x_len); + /****************************may have bugs********************************/ + + ///offset of y + y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s; + /****************************may have bugs********************************/ + y_start += y_start_offset(x_start, &(overlap_list->list[currentID].f_cigar)); + /****************************may have bugs********************************/ + + Window_Len = x_len + (threshold << 1); + + if(!determine_overlap_region(threshold, y_start, overlap_list->list[currentID].y_id, Window_Len, R_INF, + &extra_begin, &extra_end, &y_start, &o_len)) + { + append_window_list(&overlap_list->list[currentID], x_start, x_end, + -1, -1, -1, -1, -1, -1); + continue; + } + + fill_subregion(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, + R_INF, overlap_list->list[currentID].y_id, extra_begin, extra_end); + + x_string = r_string + x_start; + y_string = dumy->overlap_region; + + end_site = Reserve_Banded_BPM(y_string, Window_Len, x_string, x_len, threshold, &error); + + if (error!=(unsigned int)-1) + { + overlap_list->list[currentID].align_length += x_len; + append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, y_start + end_site, (int)error, + extra_begin, extra_end, threshold); + } + else + { + append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, -1, -1, + extra_begin, extra_end, threshold); + } + } +} + + + inline double trim_error_rate(overlap_region_alloc* overlap_list, long long ID) { long long tLen, tError,i, subWinLen, subWinNum; @@ -7745,6 +7975,63 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF, } +void correct_ul_overlap(overlap_region_alloc* overlap_list, const ma_ug_t *uref, + UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, + Graph* g, Graph* DAGCon, Cigar_record* current_cigar, + haplotype_evdience_alloc* hap, Round2_alignment* second_round, + int force_repeat, int is_consensus, int* fully_cov, int* abnormal, double max_ov_diff_ec) +{ + + clear_Correct_dumy(dumy, overlap_list); + + long long window_start, window_end; + + Window_Pool w_inf; + + init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/max_ov_diff_ec)); + + int flag = 0; + + while(get_Window(&w_inf, &window_start, &window_end) && flag != -2) + { + dumy->length = 0; + dumy->lengthNT = 0; + flag = get_interval(window_start, window_end, overlap_list, dumy); + + switch (flag) + { + case 1: ///no match here + break; + case 0: ///no match here + break; + case -2: ///if flag == -2, loop would be terminated + break; + } + + + ///dumy->lengthNT represent how many overlaps that the length of them is not equal to WINDOW; may larger or less than WINDOW + ///dumy->length represent how many overlaps that the length of them is WINDOW + ///now the windows which are larger than WINDOW are verified one-by-one, to improve it, we can do it group-bygroup + verify_ul_window(window_start, window_end, overlap_list, dumy, uref, g_read->seq); + } + + /** + // recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read); + // partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat); + recalcate_window_advance(overlap_list, R_INF, g_read, dumy, overlap_read); + partition_overlaps_advance(overlap_list, R_INF, g_read, overlap_read, dumy, hap, force_repeat); + + if(is_consensus) + { + generate_consensus(overlap_list, R_INF, g_read, dumy, g, DAGCon, current_cigar, second_round); + } + + + (*fully_cov) = check_if_fully_covered(overlap_list, R_INF, g_read, dumy, g, abnormal); + **/ +} + + void init_Cigar_record(Cigar_record* dummy) { dummy->length = 0; diff --git a/Correct.h b/Correct.h index f8efab9..f15c420 100644 --- a/Correct.h +++ b/Correct.h @@ -1106,6 +1106,12 @@ void add_new_cell_to_cigar_record_with_different_base(Cigar_record* dummy, uint3 void add_existing_cell_to_cigar_record_with_different_base(Cigar_record* dummy, uint32_t len, uint32_t type, char* seq); +void correct_ul_overlap(overlap_region_alloc* overlap_list, const ma_ug_t *uref, + UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, + Graph* g, Graph* DAGCon, Cigar_record* current_cigar, + haplotype_evdience_alloc* hap, Round2_alignment* second_round, + int force_repeat, int is_consensus, int* fully_cov, int* abnormal, double max_ov_diff_ec); + /*** type: 0. match diff --git a/Process_Read.cpp b/Process_Read.cpp index 38fb881..c3a325a 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -754,10 +754,12 @@ void init_all_ul_t(all_ul_t *x, All_reads *hR) { void destory_all_ul_t(all_ul_t *x) { uint64_t i; for (i = 0; i < x->n; i++) { - free(x->a[i].n_n); free(x->a[i].N_site.a); - free(x->a[i].r_base.a); free(x->a[i].bb.a); + free(x->a[i].N_site.a); free(x->a[i].r_base.a); free(x->a[i].bb.a); } free(x->a); + + for (i = 0; i < x->nid.n; i++) free(x->nid.a[i].a); + free(x->nid.a); } void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn, uint64_t nn_offset) @@ -843,23 +845,29 @@ void push_subblock_original_bases(char* str, all_ul_t *x, ul_vec_t *p, uint32_t void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, int64_t str_l, ul_ov_t *o, int64_t on) { int64_t i, mine, maxs, ovlp, end; ul_vec_t *p = NULL; + nid_t *np = NULL; ul_ov_t *z = NULL, *zp = NULL; uc_block_t *b = NULL; - if(rid) fprintf(stderr, "rid:%lu\n", *rid); - if(rid == NULL) { - kv_pushp(ul_vec_t, *x, &p); - memset(p, 0, sizeof(*p)); - fprintf(stderr, "x->n:%u\n", x->n); - } - else { - p = &(x->a[*rid]); + + if(id) { + kv_pushp(nid_t, x->nid, &np); + np->n = id_l; MALLOC(np->a, np->n+1); memcpy(np->a, id, id_l); np->a[id_l] = '\0'; } - if(id && id_l > 0) { - free(p->n_n); - p->n_l = id_l; MALLOC(p->n_n, p->n_l+1); memcpy(p->n_n, id, id_l); p->n_n[id_l] = '\0'; - } - if(str && str_l > 0) { + if(str) { + if(rid == NULL) { + kv_pushp(ul_vec_t, *x, &p); + memset(p, 0, sizeof(*p)); + } else { + if((*rid) >= x->m) kv_resize(ul_vec_t, *x, (*rid) + 1); + if((*rid) >= x->n) { + memset(x->a+x->n, 0, sizeof(*p)*((*rid) + 1 - x->n)); + x->n = (*rid) + 1; + } + p = &(x->a[(*rid)]); + } + + p->bb.n = p->N_site.n = p->r_base.n = 0; p->rlen = str_l; diff --git a/Process_Read.h b/Process_Read.h index d77017d..0c46d32 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -171,9 +171,11 @@ typedef struct typedef struct { - char *n_n; - uint32_t n_l; + char *a; uint32_t n; +} nid_t; +typedef struct +{ kvec_t(uint8_t) r_base; uint32_t rlen; @@ -183,6 +185,7 @@ typedef struct typedef struct { + kvec_t(nid_t) nid; ul_vec_t *a; size_t n, m; All_reads *hR; diff --git a/anchor.cpp b/anchor.cpp index bad4b55..d852e53 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -199,6 +199,133 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg ///ks_introsort_or_xs(overlap_list->length, overlap_list->list); } +void ha_get_new_ul_candidates(ha_abufl_t *ab, int64_t rid, char* rs, int64_t rl, uint64_t mz_w, uint64_t mz_k, 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, kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp) +{ + uint32_t i; + uint64_t k, l; + uint32_t high_occ = asm_opt.hom_cov >= 1?asm_opt.hom_cov:1; + + + // prepare + clear_Candidates_list(cl); + clear_overlap_region_alloc(overlap_list); + ab->mz.n = 0, ab->n_a = 0; + + // get the list of anchors + mz2_ha_sketch(rs, rl, mz_w, mz_k, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0); + + // minimizer of queried read + if (ab->mz.m > ab->old_mz_m) { + ab->old_mz_m = ab->mz.m; + REALLOC(ab->seed, ab->old_mz_m); + } + for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) { + int n; + ab->seed[i].a = ha_ptl_get(ha_idx, ab->mz.a[i].x, &n); + ab->seed[i].n = n; + ab->n_a += n; + } + if (ab->n_a > ab->m_a) { + ab->m_a = ab->n_a; + kroundup64(ab->m_a); + REALLOC(ab->a, ab->m_a); + } + for (i = 0, k = 0; i < ab->mz.n; ++i) { + int j; + ///z is one of the minimizer + ha_mzl_t *z = &ab->mz.a[i]; + seedl_t *s = &ab->seed[i]; + for (j = 0; j < s->n; ++j) { + const ha_idxposl_t *y = &s->a[j]; + anchor1_t *an = &ab->a[k++]; + uint8_t rev = z->rev == y->rev? 0 : 1; + an->other_off = y->pos; + an->self_off = rev? rl - 1 - (z->pos + 1 - z->span) : z->pos; + an->cnt = s->n; + an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off; + } + } + + // sort anchors + radix_sort_ha_an1(ab->a, ab->a + ab->n_a); + for (k = 1, l = 0; k <= ab->n_a; ++k) { + if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) { + if (k - l > 1) + radix_sort_ha_an2(ab->a + l, ab->a + k); + l = k; + } + } + + + // copy over to _cl_ + if (ab->m_a >= (uint64_t)cl->size) { + cl->size = ab->m_a; + REALLOC(cl->list, cl->size); + } + for (k = 0; k < ab->n_a; ++k) { + k_mer_hit *p = &cl->list[k]; + p->readID = ab->a[k].srt >> 33; + p->strand = ab->a[k].srt >> 32 & 1; + p->offset = ab->a[k].other_off; + p->self_offset = ab->a[k].self_off; + if(ab->a[k].cnt <= high_occ){ + p->cnt = 1; + } + else{ + p->cnt = 1 + ((ab->a[k].cnt + (high_occ<<1) - 1)/(high_occ<<1)); + p->cnt = pow(p->cnt, 1.1); + } + } + cl->length = ab->n_a; + + calculate_overlap_region_by_chaining(cl, overlap_list, chain_idx, rid, rl, &R_INF, bw_thres, keep_whole_chain, f_cigar); + + #if 0 + if (overlap_list->length > 0) { + fprintf(stderr, "B\t%ld\t%ld\t%d\n", (long)rid, (long)overlap_list->length, rlen); + 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%ld\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, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed, ha_ov_type(r, rlen)); + } + } + #endif + + if ((int)overlap_list->length > max_n_chain) { + int32_t w, n[4], s[4]; + n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0; + ks_introsort_or_ss(overlap_list->length, overlap_list->list); + for (i = 0; i < (uint32_t)overlap_list->length; ++i) { + const overlap_region *r = &overlap_list->list[i]; + w = ha_ov_type(r, rl); + ++n[w]; + if ((int)n[w] == max_n_chain) s[w] = r->shared_seed; + } + if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) { + // n[0] = n[1] = n[2] = n[3] = 0; + for (i = 0, k = 0; i < (uint32_t)overlap_list->length; ++i) { + overlap_region *r = &overlap_list->list[i]; + w = ha_ov_type(r, rl); + // ++n[w]; + // if (((int)n[w] <= max_n_chain) || (r->shared_seed >= s[w] && s[w] >= (asm_opt.k_mer_length<<1))) { + if (r->shared_seed >= s[w]) { + if ((uint32_t)k != i) { + overlap_region t; + t = overlap_list->list[k]; + overlap_list->list[k] = overlap_list->list[i]; + overlap_list->list[i] = t; + } + ++k; + } + } + overlap_list->length = k; + } + } + + ///ks_introsort_or_xs(overlap_list->length, overlap_list->list); +} + void calculate_ug_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, kvec_t_u64_warp* chain_idx, uint64_t readID, ma_utg_v *ua, double band_width_threshold, int add_beg_end, overlap_region* f_cigar, long long mz_occ, double mz_rate) @@ -626,6 +753,19 @@ void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overl ks_introsort_or_xs(overlap_list->length, overlap_list->list); } + +void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, 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, kvec_t_u64_warp* chain_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp) +{ + extern void *ha_flt_tab; + extern ha_pt_t *ha_idx; + + ha_get_new_ul_candidates(ab, rid, rs, rl, mz_w, mz_k, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar, dbg_ct, sp); + ks_introsort_or_xs(overlap_list->length, overlap_list->list); +} + + + void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list) { ks_introsort_or_xs(overlap_list->length, overlap_list->list); diff --git a/htab.cpp b/htab.cpp index 6d36cf5..a207807 100644 --- a/htab.cpp +++ b/htab.cpp @@ -1141,11 +1141,11 @@ ha_pt_t *ha_pt_ul_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_ ///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen if (flt_tab == 0) { if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1; - ha_ct_shrink(ct, 2, cutoff, asm_opt->thread_num); + ha_ct_shrink(ct, /**2**/1, cutoff, asm_opt->thread_num); } else { ///Note: here is just to remove minimizer appearing YAK_MAX_COUNT times ///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen - ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num); + ha_ct_shrink(ct, /**2**/1, YAK_MAX_COUNT - 1, asm_opt->thread_num); } pt = ha_pt_gen(ct, asm_opt->thread_num, 1); diff --git a/inter.cpp b/inter.cpp index 02616f1..d5b4b73 100644 --- a/inter.cpp +++ b/inter.cpp @@ -18,6 +18,9 @@ #include "Assembly.h" KSEQ_INIT(gzFile, gzread) +void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, 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, kvec_t_u64_warp* chain_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp); + #define MG_SEED_IGNORE (1ULL<<41) #define MG_SEED_TANDEM (1ULL<<42) #define MG_SEED_KEPT (1ULL<<43) @@ -68,6 +71,8 @@ typedef struct { int max_lc_skip, max_lc_iter, min_lc_cnt, min_lc_score, max_gc_skip, ref_bonus; int min_gc_cnt, min_gc_score, sub_diff, best_n; float chn_pen_gap, mask_level, pri_ratio; + ///base-alignment + double bw_thres, diff_ec_ul; int max_n_chain; } mg_idxopt_t; typedef struct { @@ -162,6 +167,7 @@ typedef struct { // global data structure for kt_pipeline() uint64_t total_pair; mg_gres_a hits; mg_dbn_t nn; + uint64_t num_bases, num_corrected_bases, num_recorrected_bases; } uldat_t; typedef struct { @@ -293,9 +299,10 @@ typedef struct { // data structure for each step in kt_pipeline() mg_gchains_t **gcs;///useless mg_tbuf_t **buf;///useless ha_ovec_buf_t **hab; + uint64_t num_bases, num_corrected_bases, num_recorrected_bases; } utepdat_t; -void init_mg_opt(mg_idxopt_t *opt, int is_HPC, int k, int w, int hap_n) +void init_mg_opt(mg_idxopt_t *opt, int is_HPC, int k, int w, int hap_n, int max_n_chain, double bw_thres, double diff_ec_ul) { opt->k = k; opt->w = w; @@ -320,6 +327,9 @@ void init_mg_opt(mg_idxopt_t *opt, int is_HPC, int k, int w, int hap_n) opt->sub_diff = 6; opt->best_n = 5; opt->pri_ratio = 0.8f; + opt->max_n_chain = max_n_chain; + opt->bw_thres = bw_thres; + opt->diff_ec_ul = diff_ec_ul; } void uidx_l_build(ma_ug_t *ug, mg_idxopt_t *opt, int cutoff) @@ -2116,6 +2126,30 @@ static void worker_for_ul_alignment(void *data, long i, int tid) // callback for s->opt->is_HPC, asm_opt.mz_sample_dist, asm_opt.mz_rewin, s->opt, s->uopt, &(s->gcs[i])); } +static void worker_for_ul_scall_alignment(void *data, long i, int tid) // callback for kt_for() +{ + utepdat_t *s = (utepdat_t*)data; + ha_ovec_buf_t *b = s->hab[tid]; + int64_t rid = s->id+i; + int fully_cov, abnormal; + + ha_get_ul_candidates_interface(b->abl, rid, s->seq[i], s->len[i], s->opt->w, s->opt->k, &b->olist, &b->olist_hp, &b->clist, s->opt->bw_thres, + s->opt->max_n_chain, 1, &(b->k_flag), &b->r_buf, &(b->tmp_region), NULL, &(b->sp)); + // mg_map_frag(s->ha_flt_tab, s->ha_idx, s->ug, s->rg, s->id+i, s->len[i], s->seq[i], &(s->mzs[tid]), &(s->sps[tid]), s->buf[tid], s->opt->w, s->opt->k, + // s->opt->is_HPC, asm_opt.mz_sample_dist, asm_opt.mz_rewin, s->opt, s->uopt, &(s->gcs[i])); + clear_Cigar_record(&b->cigar1); + clear_Round2_alignment(&b->round2); + + b->self_read.seq = s->seq[i]; b->self_read.length = s->len[i]; b->self_read.size = 0; + correct_ul_overlap(&b->olist, s->ug, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon, + &b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal, s->opt->diff_ec_ul); + + b->num_read_base += b->self_read.length; + b->num_correct_base += b->correct.corrected_base; + b->num_recorrect_base += b->round2.dumy.corrected_base; + memset(&b->self_read, 0, sizeof(b->self_read)); +} + void dump_gaf(mg_gres_a *hits, const mg_gchains_t *gs, uint32_t only_p) { if (gs == NULL || gs->n_gc == 0 || gs->n_lc == 0) return; @@ -2297,6 +2331,7 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac uint64_t i; CALLOC(s->hab, p->n_thread); for (i = 0; i < p->n_thread; ++i) s->hab[i] = ha_ovec_init(0, 0, 1); + kt_for(p->n_thread, worker_for_ul_scall_alignment, s, s->n); ///debug /** uint64_t i; @@ -2319,7 +2354,12 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac free(s->mzs[i].a); free(s->sps[i].a); } **/ - for (i = 0; i < p->n_thread; ++i) ha_ovec_destroy(s->hab[i]); + for (i = 0; i < p->n_thread; ++i) { + s->num_bases += s->hab[i]->num_read_base; + s->num_corrected_bases += s->hab[i]->num_correct_base; + s->num_recorrected_bases += s->hab[i]->num_recorrect_base; + ha_ovec_destroy(s->hab[i]); + } free(s->hab); // free(s->buf); free(s->mzs); free(s->sps); return s; @@ -2327,6 +2367,9 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac else if (step == 2) { // step 3: dump utepdat_t *s = (utepdat_t*)in; uint64_t i, rid; + p->num_bases += s->num_bases; + p->num_corrected_bases += s->num_corrected_bases; + p->num_recorrected_bases += s->num_recorrected_bases; for (i = 0; i < (uint64_t)s->n; ++i) { ///debug /** @@ -2367,8 +2410,9 @@ int scall_ul_pipeline(uldat_t* sl, const enzyme *fn) sl->hits.total_base = sl->total_base; sl->hits.total_pair = sl->total_pair; fprintf(stderr, "[M::%s::%.3f] ==> Qualification\n", __func__, yak_realtime()-index_time); - fprintf(stderr, "[M::%s::%.3f] ==> # reads: %lu, # bases: %lu\n", __func__, yak_realtime()-index_time, - UL_INF.n, sl->total_base); + fprintf(stderr, "[M::%s::] ==> # reads: %lu, # bases: %lu\n", __func__, UL_INF.n, sl->total_base); + fprintf(stderr, "[M::%s::] ==> # bases: %lu; # corrected bases: %lu; # recorrected bases: %lu\n", + __func__, sl->num_bases, sl->num_corrected_bases, sl->num_recorrected_bases); return 1; } @@ -2377,13 +2421,11 @@ int scall_ul_pipeline(uldat_t* sl, const enzyme *fn) int print_ul_rs(all_ul_t *U_INF) { uint32_t i; - ul_vec_t *p = NULL; UC_Read ur; init_UC_Read(&ur); for (i = 0; i < U_INF->n; i++) { - p = &(U_INF->a[i]); retrieve_ul_t(&ur, NULL, U_INF, i, 0, 0, -1); - fprintf(stderr, ">%s\n", p->n_n); + fprintf(stderr, ">%s\n", U_INF->nid.a[i].a); fprintf(stderr, "%.*s\n", (int)ur.length, ur.seq); } @@ -2864,7 +2906,7 @@ void ul_resolve(ma_ug_t *ug, const asg_t *rg, const ug_opt_t *uopt, int hap_n) { fprintf(stderr, "[M::%s::] ==> UL\n", __func__); mg_idxopt_t opt; - init_mg_opt(&opt, 0, 19, 10, hap_n); + init_mg_opt(&opt, 0, 19, 10, hap_n, 0, 0, 0.05); int exist = (asm_opt.load_index_from_disk? uidx_load(&ha_flt_tab, &ha_idx, asm_opt.output_file_name) : 0); if(exist == 0) uidx_build(ug, &opt); if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name); @@ -2884,7 +2926,7 @@ int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const enzyme *fn, void *ha sl.uopt = uopt; scall_ul_pipeline(&sl, fn); // print_ul_rs(&UL_INF); - debug_retrieve_rc_sub(&UL_INF, &R_INF, &(ug->u), 100); + // debug_retrieve_rc_sub(&UL_INF, &R_INF, &(ug->u), 100); // if(!load_ul_hits(&sl.hits, &sl.nn, asm_opt.output_file_name)) { // scall_ul_pipeline(&sl, fn); // write_ul_hits(&sl.hits, &sl.nn, asm_opt.output_file_name); @@ -2952,14 +2994,15 @@ void ul_load(const ug_opt_t *uopt) fprintf(stderr, "[M::%s::] ==> UL\n", __func__); mg_idxopt_t opt; ma_ug_t *ug = dedup_HiFis(uopt->sources, uopt->min_ovlp, uopt->max_hang, uopt->gap_fuzz); - int cutoff = asm_opt.hom_cov * asm_opt.high_factor; - init_aux_table(); - init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff); - /** + int cutoff; + init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov); + cutoff = asm_opt.max_n_chain; + init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, 0.05, 0.05); + int exist = (asm_opt.load_index_from_disk? uidx_load(&ha_flt_tab, &ha_idx, asm_opt.output_file_name) : 0); if(exist == 0) uidx_l_build(ug, &opt, cutoff); - if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name); - **/ + if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name); + ul_v_call(&opt, uopt, asm_opt.ar, ha_flt_tab, ha_idx, ug); ma_ug_destroy(ug); destory_all_ul_t(&UL_INF); } \ No newline at end of file