diff --git a/Assembly.cpp b/Assembly.cpp index f789b47..e61427c 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -21,7 +21,7 @@ void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list); All_reads R_INF; Debug_reads R_INF_FLAG; -all_ul_t UL_INF; +all_ul_t UL_INF, ULG_INF; uint32_t *het_cnt = NULL; void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length, char* new_read, int* new_length) diff --git a/CommandLines.h b/CommandLines.h index 0c0b4db..469c110 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.16.4-r393" +#define HA_VERSION "0.16.4-r396" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index 81e91df..0e5153c 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -8245,6 +8245,16 @@ int asg_arc_del_false_node(asg_t *g, ma_hit_t_alloc* sources, int max_ext) return n_cut; } +int64_t count_edges_v_w(asg_t *g, uint32_t v, uint32_t w) +{ + asg_arc_t* av = asg_arc_a(g, v); uint32_t nv = asg_arc_n(g, v), i, occ = 0; + for (i = 0; i < nv; i++) { + if(av[i].del) continue; + if(av[i].v == w) occ++; + } + return occ; +} + void update_ug_ou(ma_ug_t *ug, asg_t *sg) { uint32_t k, i, uv, uw, rv, rw, nv; asg_arc_t *ue, *av; @@ -9451,7 +9461,7 @@ UC_Read* r_read, UC_Read* q_read, int max_hang, int min_ovlp, kvec_asg_arc_t_war ma_ug_t *gen_polished_ug(const ug_opt_t *uopt, asg_t *sg) { kvec_asg_arc_t_warp e, d; - uint32_t i, k; ma_utg_t *u; + uint32_t i; ma_utg_t *u; kv_init(e.a); kv_init(d.a); ma_ug_t *ug = ma_ug_gen(sg); UC_Read g_read, tmp; @@ -9468,7 +9478,7 @@ ma_ug_t *gen_polished_ug(const ug_opt_t *uopt, asg_t *sg) destory_UC_Read(&g_read); destory_UC_Read(&tmp); uint32_t n_vtx = ug->g->n_seq*2, v, nv, vLen = 0; - asg_arc_t* av = NULL, *p = NULL; + asg_arc_t* av = NULL; for (v = 0; v < n_vtx; ++v) { if (ug->g->seq[v>>1].del) continue; av = asg_arc_a(ug->g, v); nv = asg_arc_n(ug->g, v); @@ -9480,14 +9490,14 @@ ma_ug_t *gen_polished_ug(const ug_opt_t *uopt, asg_t *sg) } } - if(d.a.n > 0) { - for (k = 0; k < d.a.n; k++) { - p = asg_arc_pushp(sg); - *p = d.a.a[k]; - } - free(sg->idx); sg->idx = 0; sg->is_srt = 0; - asg_cleanup(sg); - } + // if(d.a.n > 0) { + // for (k = 0; k < d.a.n; k++) { + // p = asg_arc_pushp(sg); + // *p = d.a.a[k]; + // } + // free(sg->idx); sg->idx = 0; sg->is_srt = 0; + // asg_cleanup(sg); + // } kv_destroy(e.a); kv_destroy(d.a); return ug; @@ -17718,6 +17728,53 @@ int load_ma_hit_ts(ma_hit_t_alloc** x, char* read_file_name) { read_ma(&((*x)[i].buffer[k]), fp); } + // fread((*x)[i].buffer, sizeof((*((*x)[i].buffer))), (*x)[i].length, fp); + } + + free(index_name); + fclose(fp); + fprintf(stderr, "ma_hit_ts has been read.\n"); + + return 1; +} + + +int load_debug_ma_hit_ts(ma_hit_t_alloc** x, char* read_file_name) +{ + fprintf(stderr, "Loading ma_hit_ts from disk... \n"); + char* index_name = (char*)malloc(strlen(read_file_name)+15); + sprintf(index_name, "%s.bin", read_file_name); + FILE* fp = fopen(index_name, "r"); + if(!fp) + { + return 0; + } + + + long long n_read; + long long i/**, k**/; + int f_flag; + f_flag = fread(&n_read, sizeof(n_read), 1, fp); + (*x) = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*n_read); + + + for (i = 0; i < n_read; i++) + { + f_flag += fread(&((*x)[i].is_fully_corrected), sizeof((*x)[i].is_fully_corrected), 1, fp); + f_flag += fread(&((*x)[i].is_abnormal), sizeof((*x)[i].is_abnormal), 1, fp); + f_flag += fread(&((*x)[i].length), sizeof((*x)[i].length), 1, fp); + (*x)[i].size = (*x)[i].length; + + (*x)[i].buffer = NULL; + if((*x)[i].length == 0) continue; + + (*x)[i].buffer = (ma_hit_t*)malloc(sizeof(ma_hit_t)*(*x)[i].length); + + // for (k = 0; k < (*x)[i].length; k++) + // { + // read_ma(&((*x)[i].buffer[k]), fp); + // } + fread((*x)[i].buffer, sizeof((*((*x)[i].buffer))), (*x)[i].length, fp); } free(index_name); @@ -17769,6 +17826,37 @@ void write_ma_hit_ts(ma_hit_t_alloc* x, long long n_read, char* read_file_name) { write_ma(x[i].buffer + k, fp); } + // fwrite(x[i].buffer, sizeof((*(x[i].buffer))), x[i].length, fp); + } + + + free(index_name); + fflush(fp); + fclose(fp); + fprintf(stderr, "ma_hit_ts has been written.\n"); +} + + +void write_debug_ma_hit_ts(ma_hit_t_alloc* x, long long n_read, char* read_file_name) +{ + fprintf(stderr, "Writing ma_hit_ts to disk... \n"); + char* index_name = (char*)malloc(strlen(read_file_name)+15); + sprintf(index_name, "%s.bin", read_file_name); + FILE* fp = fopen(index_name, "w"); + long long i/**, k**/; + fwrite(&n_read, sizeof(n_read), 1, fp); + + + for (i = 0; i < n_read; i++) + { + fwrite(&(x[i].is_fully_corrected), sizeof(x[i].is_fully_corrected), 1, fp); + fwrite(&(x[i].is_abnormal), sizeof(x[i].is_abnormal), 1, fp); + fwrite(&(x[i].length), sizeof(x[i].length), 1, fp); + // for (k = 0; k < x[i].length; k++) + // { + // write_ma(x[i].buffer + k, fp); + // } + fwrite(x[i].buffer, sizeof((*(x[i].buffer))), x[i].length, fp); } @@ -24460,18 +24548,20 @@ int load_coverage_cut(ma_sub_t** coverage_cut, char* read_file_name) f_flag += fread(&n_read, sizeof(n_read), 1, fp); (*coverage_cut) = (ma_sub_t*)malloc(sizeof(ma_sub_t)*n_read); - uint64_t i = 0, tmp; - for (i = 0; i < n_read; i++) - { - f_flag += fread(&tmp, sizeof(tmp), 1, fp); - (*coverage_cut)[i].c = tmp; - f_flag += fread(&tmp, sizeof(tmp), 1, fp); - (*coverage_cut)[i].del = tmp; - f_flag += fread(&tmp, sizeof(tmp), 1, fp); - (*coverage_cut)[i].e = tmp; - f_flag += fread(&tmp, sizeof(tmp), 1, fp); - (*coverage_cut)[i].s = tmp; - } + // uint64_t i = 0, tmp; + // for (i = 0; i < n_read; i++) + // { + // f_flag += fread(&tmp, sizeof(tmp), 1, fp); + // (*coverage_cut)[i].c = tmp; + // f_flag += fread(&tmp, sizeof(tmp), 1, fp); + // (*coverage_cut)[i].del = tmp; + // f_flag += fread(&tmp, sizeof(tmp), 1, fp); + // (*coverage_cut)[i].e = tmp; + // f_flag += fread(&tmp, sizeof(tmp), 1, fp); + // (*coverage_cut)[i].s = tmp; + // } + fread((*coverage_cut), sizeof((*((*coverage_cut)))), n_read, fp); + free(index_name); fflush(fp); fclose(fp); @@ -24485,18 +24575,19 @@ int write_coverage_cut(ma_sub_t* coverage_cut, char* read_file_name, uint64_t n_ sprintf(index_name, "%s.bin", read_file_name); FILE* fp = fopen(index_name, "w"); fwrite(&n_read, sizeof(n_read), 1, fp); - uint64_t i = 0, tmp; - for (i = 0; i < n_read; i++) - { - tmp = coverage_cut[i].c; - fwrite(&tmp, sizeof(tmp), 1, fp); - tmp = coverage_cut[i].del; - fwrite(&tmp, sizeof(tmp), 1, fp); - tmp = coverage_cut[i].e; - fwrite(&tmp, sizeof(tmp), 1, fp); - tmp = coverage_cut[i].s; - fwrite(&tmp, sizeof(tmp), 1, fp); - } + // uint64_t i = 0, tmp; + // for (i = 0; i < n_read; i++) + // { + // tmp = coverage_cut[i].c; + // fwrite(&tmp, sizeof(tmp), 1, fp); + // tmp = coverage_cut[i].del; + // fwrite(&tmp, sizeof(tmp), 1, fp); + // tmp = coverage_cut[i].e; + // fwrite(&tmp, sizeof(tmp), 1, fp); + // tmp = coverage_cut[i].s; + // fwrite(&tmp, sizeof(tmp), 1, fp); + // } + fwrite(coverage_cut, sizeof((*(coverage_cut))), n_read, fp); free(index_name); fflush(fp); fclose(fp); @@ -24557,7 +24648,7 @@ int write_asg_t(asg_t *sg, char* read_file_name) char* index_name = (char*)malloc(strlen(read_file_name)+15); sprintf(index_name, "%s.bin", read_file_name); FILE* fp = fopen(index_name, "w"); - uint32_t tmp, i; + uint32_t tmp/**, i**/; tmp = sg->n_arc; fwrite(&tmp, sizeof(tmp), 1, fp); @@ -24588,37 +24679,39 @@ int write_asg_t(asg_t *sg, char* read_file_name) fwrite(sg->idx, sizeof(sg->idx[0]), Len, fp); - for (i = 0; i < sg->n_arc; i++) - { - tmp = sg->arc[i].del; - fwrite(&tmp, sizeof(tmp), 1, fp); - tmp = sg->arc[i].el; - fwrite(&tmp, sizeof(tmp), 1, fp); - tmp = sg->arc[i].no_l_indel; - fwrite(&tmp, sizeof(tmp), 1, fp); - tmp = sg->arc[i].ol; - fwrite(&tmp, sizeof(tmp), 1, fp); - tmp = sg->arc[i].strong; - fwrite(&tmp, sizeof(tmp), 1, fp); + // for (i = 0; i < sg->n_arc; i++) + // { + // tmp = sg->arc[i].del; + // fwrite(&tmp, sizeof(tmp), 1, fp); + // tmp = sg->arc[i].el; + // fwrite(&tmp, sizeof(tmp), 1, fp); + // tmp = sg->arc[i].no_l_indel; + // fwrite(&tmp, sizeof(tmp), 1, fp); + // tmp = sg->arc[i].ol; + // fwrite(&tmp, sizeof(tmp), 1, fp); + // tmp = sg->arc[i].strong; + // fwrite(&tmp, sizeof(tmp), 1, fp); - uint64_t tmp_64; - tmp_64 = sg->arc[i].ul; - fwrite(&tmp_64, sizeof(tmp_64), 1, fp); + // uint64_t tmp_64; + // tmp_64 = sg->arc[i].ul; + // fwrite(&tmp_64, sizeof(tmp_64), 1, fp); - tmp = sg->arc[i].v; - fwrite(&tmp, sizeof(tmp), 1, fp); - } + // tmp = sg->arc[i].v; + // fwrite(&tmp, sizeof(tmp), 1, fp); + // } + fwrite(sg->arc, sizeof((*(sg->arc))), sg->n_arc, fp); - for (i = 0; i < sg->n_seq; i++) - { - tmp = sg->seq[i].c; - fwrite(&tmp, sizeof(tmp), 1, fp); - tmp = sg->seq[i].del; - fwrite(&tmp, sizeof(tmp), 1, fp); - tmp = sg->seq[i].len; - fwrite(&tmp, sizeof(tmp), 1, fp); - } + // for (i = 0; i < sg->n_seq; i++) + // { + // tmp = sg->seq[i].c; + // fwrite(&tmp, sizeof(tmp), 1, fp); + // tmp = sg->seq[i].del; + // fwrite(&tmp, sizeof(tmp), 1, fp); + // tmp = sg->seq[i].len; + // fwrite(&tmp, sizeof(tmp), 1, fp); + // } + fwrite(sg->seq, sizeof((*(sg->seq))), sg->n_seq, fp); free(index_name); fflush(fp); @@ -24637,7 +24730,7 @@ int load_asg_t(asg_t **sg, char* read_file_name) { return 0; } - uint32_t tmp, i; + uint32_t tmp/**, i**/; (*sg) = (asg_t*)calloc(1, sizeof(asg_t)); int f_flag = 0; @@ -24676,41 +24769,43 @@ int load_asg_t(asg_t **sg, char* read_file_name) (*sg)->seq = (asg_seq_t*)malloc(sizeof(asg_seq_t)*(*sg)->m_seq); - for (i = 0; i < (*sg)->n_arc; i++) - { - f_flag += fread(&tmp, sizeof(tmp), 1, fp); - (*sg)->arc[i].del = tmp; + // for (i = 0; i < (*sg)->n_arc; i++) + // { + // f_flag += fread(&tmp, sizeof(tmp), 1, fp); + // (*sg)->arc[i].del = tmp; - f_flag += fread(&tmp, sizeof(tmp), 1, fp); - (*sg)->arc[i].el = tmp; + // f_flag += fread(&tmp, sizeof(tmp), 1, fp); + // (*sg)->arc[i].el = tmp; - f_flag += fread(&tmp, sizeof(tmp), 1, fp); - (*sg)->arc[i].no_l_indel = tmp; + // f_flag += fread(&tmp, sizeof(tmp), 1, fp); + // (*sg)->arc[i].no_l_indel = tmp; - f_flag += fread(&tmp, sizeof(tmp), 1, fp); - (*sg)->arc[i].ol = tmp; + // f_flag += fread(&tmp, sizeof(tmp), 1, fp); + // (*sg)->arc[i].ol = tmp; - f_flag += fread(&tmp, sizeof(tmp), 1, fp); - (*sg)->arc[i].strong = tmp; + // f_flag += fread(&tmp, sizeof(tmp), 1, fp); + // (*sg)->arc[i].strong = tmp; - uint64_t tmp_64; - f_flag += fread(&tmp_64, sizeof(tmp_64), 1, fp); - (*sg)->arc[i].ul = tmp_64; + // uint64_t tmp_64; + // f_flag += fread(&tmp_64, sizeof(tmp_64), 1, fp); + // (*sg)->arc[i].ul = tmp_64; - f_flag += fread(&tmp, sizeof(tmp), 1, fp); - (*sg)->arc[i].v = tmp; - } + // f_flag += fread(&tmp, sizeof(tmp), 1, fp); + // (*sg)->arc[i].v = tmp; + // } + fread((*sg)->arc, sizeof((*((*sg)->arc))), (*sg)->n_arc, fp); - for (i = 0; i < (*sg)->n_seq; i++) - { - f_flag += fread(&tmp, sizeof(tmp), 1, fp); - (*sg)->seq[i].c = tmp; - f_flag += fread(&tmp, sizeof(tmp), 1, fp); - (*sg)->seq[i].del = tmp; - f_flag += fread(&tmp, sizeof(tmp), 1, fp); - (*sg)->seq[i].len = tmp; - } + // for (i = 0; i < (*sg)->n_seq; i++) + // { + // f_flag += fread(&tmp, sizeof(tmp), 1, fp); + // (*sg)->seq[i].c = tmp; + // f_flag += fread(&tmp, sizeof(tmp), 1, fp); + // (*sg)->seq[i].del = tmp; + // f_flag += fread(&tmp, sizeof(tmp), 1, fp); + // (*sg)->seq[i].len = tmp; + // } + fread((*sg)->seq, sizeof((*((*sg)->seq))), (*sg)->n_seq, fp); free(index_name); fflush(fp); @@ -24720,18 +24815,18 @@ int load_asg_t(asg_t **sg, char* read_file_name) } int write_debug_graph(asg_t *sg, ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, -char* output_file_name, long long n_read, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex) +char* output_file_name, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex) { char* gfa_name = (char*)malloc(strlen(output_file_name)+55); ////write_All_reads(&R_INF, gfa_name); sprintf(gfa_name, "%s.all.debug.source", output_file_name); - write_ma_hit_ts(sources, R_INF.total_reads, gfa_name); + write_debug_ma_hit_ts(sources, R_INF.total_reads, gfa_name); sprintf(gfa_name, "%s.all.debug.reverse", output_file_name); - write_ma_hit_ts(reverse_sources, R_INF.total_reads, gfa_name); + write_debug_ma_hit_ts(reverse_sources, R_INF.total_reads, gfa_name); sprintf(gfa_name, "%s.all.debug.coverage_cut", output_file_name); - write_coverage_cut(coverage_cut, gfa_name, n_read); + write_coverage_cut(coverage_cut, gfa_name, R_INF.total_reads); sprintf(gfa_name, "%s.all.debug.ruIndex", output_file_name); write_ruIndex(ruIndex, gfa_name); sprintf(gfa_name, "%s.all.debug.asg_t", output_file_name); @@ -24791,13 +24886,13 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex) sprintf(gfa_name, "%s.all.debug.source", output_file_name); - if(!load_ma_hit_ts(sources, gfa_name)) + if(!load_debug_ma_hit_ts(sources, gfa_name)) { return 0; } sprintf(gfa_name, "%s.all.debug.reverse", output_file_name); - if(!load_ma_hit_ts(reverse_sources, gfa_name)) + if(!load_debug_ma_hit_ts(reverse_sources, gfa_name)) { return 0; } @@ -31392,7 +31487,14 @@ ma_sub_t **coverage_cut_ptr, int debug_g) output_read_graph(sg, coverage_cut, unlean_name, n_read); free(unlean_name); } - + + + if (asm_opt.flag & HA_F_VERBOSE_GFA) + { + write_debug_graph(sg, sources, coverage_cut, output_file_name, reverse_sources, ruIndex); + debug_gfa:; + if(asm_opt.ar) create_ul_info(sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex); + } gen_ug_opt_t(&uopt, sources, reverse_sources, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, readLen, coverage_cut, ruIndex); ul_clean_gfa(&uopt, sg, sources, reverse_sources, ruIndex, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.6, asm_opt.max_short_tip, &b_mask_t, !!asm_opt.ar, ha_opt_triobin(&asm_opt), UL_COV_THRES, o_file); @@ -31529,11 +31631,12 @@ ma_sub_t **coverage_cut_ptr, int debug_g) output_contig_graph_primary_pre(sg, coverage_cut, o_file, sources, reverse_sources, asm_opt.small_pop_bubble_size, asm_opt.max_short_tip, ruIndex, max_hang_length, mini_overlap_length); + /** if (asm_opt.flag & HA_F_VERBOSE_GFA) { write_debug_graph(sg, sources, coverage_cut, output_file_name, n_read, reverse_sources, ruIndex); debug_gfa:; - } + }**/ if(asm_opt.fn_bin_poy) { diff --git a/Overlaps.h b/Overlaps.h index 4a397c9..b750254 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -1086,6 +1086,7 @@ void output_unitig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_na ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp); void flat_soma_v(asg_t *sg, ma_hit_t_alloc* sources, R_to_U* ruIndex); void hic_clean(asg_t* read_g); +int64_t count_edges_v_w(asg_t *g, uint32_t v, uint32_t w); #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/Process_Read.h b/Process_Read.h index 7352be5..d7cd042 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -202,6 +202,7 @@ typedef struct } all_ul_t; extern all_ul_t UL_INF; +extern all_ul_t ULG_INF; // extern uint32_t *het_cnt; void init_All_reads(All_reads* r); diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 7303c05..bd74f8c 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -1328,6 +1328,8 @@ void fill_containment_by_ul(asg_t *g, ma_hit_t_alloc *src) } + + void ul_clean_gfa(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, int64_t clean_round, double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, double ou_drop_rate, int64_t max_tip, bub_label_t *b_mask_t, int32_t is_ou, int32_t is_trio, uint32_t ou_thres, char *o_file) { @@ -1345,6 +1347,8 @@ double ou_drop_rate, int64_t max_tip, bub_label_t *b_mask_t, int32_t is_ou, int3 // debug_info_of_specfic_node("m64011_190830_220126/95028102/ccs", sg, rI, "beg"); asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + // fprintf(stderr, "[M::%s] count_edges_v_w(sg, 49778, 49847)->%ld\n", __func__, count_edges_v_w(sg, 49778, 49847)); + for (i = 0; i < clean_round; i++, drop += step) { if(drop > max_ovlp_drop_ratio) drop = max_ovlp_drop_ratio; // fprintf(stderr, "(0):i->%ld, drop->%f\n", i, drop); @@ -1406,7 +1410,6 @@ double ou_drop_rate, int64_t max_tip, bub_label_t *b_mask_t, int32_t is_ou, int3 // mini_overlap_length, bubble_dist, 10); set_hom_global_coverage(&asm_opt, sg, uopt->coverage_cut, src, rev, rI, uopt->max_hang, uopt->min_ovlp); rescue_bubble_by_chain(sg, uopt->coverage_cut, src, rev, (asm_opt.max_short_tip*2), 0.15, 3, rI, 0.05, 0.9, uopt->max_hang, uopt->min_ovlp, 10, uopt->gap_fuzz, b_mask_t); - output_unitig_graph(sg, uopt->coverage_cut, o_file, src, rI, uopt->max_hang, uopt->min_ovlp); // flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL; flat_soma_v(sg, src, rI); @@ -1418,6 +1421,10 @@ double ou_drop_rate, int64_t max_tip, bub_label_t *b_mask_t, int32_t is_ou, int3 } if(is_ou) { + for (i = 0; i < (int64_t)sg->n_seq; ++i) { + if(sg->seq[i].del) continue; + sg->seq[i].c = PRIMARY_LABLE; + } hic_clean(sg); ul_realignment(uopt, sg); // if(ul_refine_alignment(uopt, sg)) update_sg_uo(sg, src); diff --git a/inter.cpp b/inter.cpp index bd44e60..7e8e083 100644 --- a/inter.cpp +++ b/inter.cpp @@ -6621,7 +6621,7 @@ int64_t *n_u_, int64_t *n_v_) *n_u_ = *n_v_ = 0; for (i = 0, k = 0; i < n; ++i) { if(f[i] >= 0) { - srt[k] = (uint64_t)f[i]; srt[k] <<= 32; srt[k] |= ((uint64_t)i)<<1; k++; + srt[k] = (uint64_t)f[i]; srt[k] <<= 32; srt[k] |= (((uint64_t)i)<<1); k++; } } n_srt = k; @@ -7346,9 +7346,9 @@ void update_exist_chain_adv(const ul_idx_t *uref, ul_ov_t *ch, uint64_t *idx, in void dedup_second_chain_adv(const ul_idx_t *uref, ul_ov_t *gb, int64_t gb_n, mg_lchain_t *chain_a, -kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn, uint64_t *b, int64_t qlen) +kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn, uint64_t *b, int64_t qlen, int64_t ulid) { - int64_t k, i, z, ss, ee, b_n, n_s; uint64_t qs, qe, ts, te; uint32_t mk = 0x80000000, pi; mg_lchain_t nn; + int64_t k, i, z, ss, ee, b_n, n_s; uint64_t qs, qe, ts, te; uint32_t mk = 0x80000000/**, pi**/; mg_lchain_t nn; int64_t iqs, iqe, its, ite, tsc; for (z = gb_n - 1; z >= 0; z--) {///start from the best chain @@ -7373,6 +7373,24 @@ kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn, uint64_t *b, int64_t qlen) assert(b_n > 0); assert(raw_idx->a[chain_a[k].off].qs == qs && raw_idx->a[chain_a[k].off].qe == qe && raw_idx->a[chain_a[k].off].ts == ts && raw_idx->a[chain_a[k].off].te == te); + // if(!(raw_idx->a[chain_a[k].off].qs == qs && raw_idx->a[chain_a[k].off].qe == qe && + // raw_idx->a[chain_a[k].off].ts == ts && raw_idx->a[chain_a[k].off].te == te)) { + // fprintf(stderr, "\n[M::%s::ulid->%ld******] raw_idx_offset:%d, qs:%lu, qe:%lu, ts:%lu, te:%lu, raw_idx->qs:%u, raw_idx->qe:%u, raw_idx->ts:%u, raw_idx->te:%u\n", + // __func__, ulid, chain_a[k].off, qs, qe, ts, te, + // raw_idx->a[chain_a[k].off].qs, raw_idx->a[chain_a[k].off].qe, + // raw_idx->a[chain_a[k].off].ts, raw_idx->a[chain_a[k].off].te); + + // for (i = raw_idx->a[chain_a[k].off].qn;i>=0;) { + // if(((raw_chn->a[i].qs<<1)>>1) < qs) qs = ((raw_chn->a[i].qs<<1)>>1); + // if(raw_chn->a[i].ts < ts) ts = raw_chn->a[i].ts; + // if(raw_chn->a[i].qe > qe) qe = raw_chn->a[i].qe; + // if(raw_chn->a[i].te > te) te = raw_chn->a[i].te; + // fprintf(stderr, "[M::%s->pieces (%ld)] qs->%u, qe->%u, ts->%u, te->%u\n", __func__, i, + // ((raw_chn->a[i].qs<<1)>>1), raw_chn->a[i].qe, raw_chn->a[i].ts, raw_chn->a[i].te); + // if(raw_chn->a[i].tn == (uint32_t)-1) i = -1; + // else i = raw_chn->a[i].tn; + // } + // } update_exist_chain_adv(uref, raw_chn->a, b, b_n, raw_idx->a[chain_a[k].off].tn, &nn); nn.v = (raw_idx->a[chain_a[k].off].tn<<1)|raw_idx->a[chain_a[k].off].rev; @@ -7382,10 +7400,10 @@ kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn, uint64_t *b, int64_t qlen) tsc += (chain_a[k].score - nn.score); if(n_s) { - raw_idx->a[chain_a[k].off].qn = b[b_n-1]; - for (i = 0, pi = (uint32_t)-1; i < b_n; i++) { - raw_chn->a[b[i]].tn = pi; pi = b[i]; - } + // raw_idx->a[chain_a[k].off].qn = b[b_n-1]; + // for (i = 0, pi = (uint32_t)-1; i < b_n; i++) { + // raw_chn->a[b[i]].tn = pi; pi = b[i]; + // } ///don't update chain_a[k] as it will be used for taceback in the next step chain_a[k].score = nn.score; chain_a[k].qs = nn.qs; chain_a[k].qe = nn.qe; @@ -7452,7 +7470,7 @@ const asg_t *g, st_mt_t *dst_done, vec_sp_node_t *out, vec_mg_pathv_t *res, uint } if(is_done == 0) { - for (i = gb_n - 2; i >= 0; i--) { + for (i = gb_n - 2; i >= 0; i--) {///from the second best chain p = &(gb[i]); ovlp = ((MIN(m->qe, p->qe) > MAX(m->qs, p->qs))? (MIN(m->qe, p->qe) - MAX(m->qs, p->qs)):0); novlp = (p->qe - p->qs) - ovlp; @@ -7466,7 +7484,7 @@ const asg_t *g, st_mt_t *dst_done, vec_sp_node_t *out, vec_mg_pathv_t *res, uint if(is_done == 0) { if(((m->qe - m->qs) > (qlen*primary_fragment_cov_rate)) || ((m->qe - m->qs) > mini_primary_fragment_len)) { // dedup_second_chain(uref, idx->a, idx_n, m.ts, m.te, a, raw_idx, raw_chn, b, bw, diff_ec_ul, qlen); - if(raw_chn && raw_idx) dedup_second_chain_adv(uref, gb, gb_n, a, raw_idx, raw_chn, b, qlen); + if(raw_chn && raw_idx) dedup_second_chain_adv(uref, gb, gb_n, a, raw_idx, raw_chn, b, qlen, ulid); for (k = gb_n-1, n_mchain = 0; k >= 0; k--) { m = &(gb[k]);///max chain // fprintf(stderr, "++[M::%s::k->%ld] score->%u, qs->%u, qe->%u\n", __func__, k, m->qn, m->qs, m->qe); @@ -8017,12 +8035,39 @@ vec_mg_lchain_t *uc, vec_mg_lchain_t *swap) dd_ul_vec_t(uref, swap->a, swap->n, rch); } +void print_ru_raw_chains(kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn, vec_mg_lchain_t *gch, ul_vec_t *rch, ma_ug_t *ug) +{ + int64_t k, i; uint64_t ts, te, qs, qe; + for (k = 0; k < (int64_t)gch->n; k++) { + i = raw_idx->a[gch->a[k].off].qn; + qs = ((raw_chn->a[i].qs<<1)>>1); qe = raw_chn->a[i].qe; + fprintf(stderr, "\n[M::%s->overall chain (%ld)] utg%.6d%c(%c), qs->%u, qe->%u, qlen->%u, ts->%u, te->%u, tlen->%u\n", __func__, k, + (int32_t)(gch->a[k].v>>1)+1, "lc"[ug->u.a[(gch->a[k].v>>1)].circ], "+-"[(gch->a[k].v&1)], + raw_idx->a[gch->a[k].off].qs, raw_idx->a[gch->a[k].off].qe, rch->rlen, + raw_idx->a[gch->a[k].off].ts, raw_idx->a[gch->a[k].off].te, ug->u.a[(gch->a[k].v>>1)].len); + ts = raw_chn->a[i].ts; te = raw_chn->a[i].te; + for (;i>=0;) { + if(((raw_chn->a[i].qs<<1)>>1) < qs) qs = ((raw_chn->a[i].qs<<1)>>1); + if(raw_chn->a[i].ts < ts) ts = raw_chn->a[i].ts; + if(raw_chn->a[i].qe > qe) qe = raw_chn->a[i].qe; + if(raw_chn->a[i].te > te) te = raw_chn->a[i].te; + + fprintf(stderr, "[M::%s->chain pieces (%ld)] qs->%u, qe->%u, ts->%u, te->%u\n", __func__, i, + ((raw_chn->a[i].qs<<1)>>1), raw_chn->a[i].qe, raw_chn->a[i].ts, raw_chn->a[i].te); + + if(raw_chn->a[i].tn == (uint32_t)-1) i = -1; + else i = raw_chn->a[i].tn; + } + + } + +} ///sps and hap are just vector for uint64_t; used for buffer uint32_t direct_gchain(mg_tbuf_t *b, ul_vec_t *rch, glchain_t *ll, gdpchain_t *gdp, st_mt_t *sps, haplotype_evdience_alloc *hap, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid) { - // if(ulid != 814) return 0; + // if(ulid != 7768/** && ulid != 44522**/) return 0; kv_ul_ov_t *idx = &(ll->lo), *init = &(ll->tk); int64_t max_idx; idx->n = init->n = 0; gl_rg2ug_gen(rch, idx, uref, 1, 2); @@ -8038,7 +8083,8 @@ int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid) // fprintf(stderr, "\n+++[M::%s::id->%ld, len->%u] idx->n:%lu\n", __func__, ulid, rch->rlen, (uint64_t)idx->n); // kv_resize(uint64_t, ll->srt.a, idx->n); kv_resize(uint64_t, hap->snp_srt, idx->n); kv_resize(uint64_t, gdp->v, idx->n); // occ = gl_chain_advance(&(gdp->l), &(gdp->swap), uref, uopt, G_CHAIN_BW, diff_ec_ul, qlen, UG_SKIP, dumy->overlapID, ll->srt.a.a, hap->snp_srt.a, G_CHAIN_TRANS_WEIGHT, 0, NULL, uref->ug, debug_i, km); - + // print_ru_raw_chains(idx, init, &(gdp->l), rch, uref->ug); + ///buffer kv_resize(uint64_t, ll->srt.a, gdp->l.n); kv_resize(uint64_t, hap->snp_srt, gdp->l.n); kv_resize(uint64_t, gdp->v, gdp->l.n); kv_resize(int64_t, gdp->f, gdp->l.n);