From 46e899bbae649da0475405958516fe84730fbb20 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Thu, 1 Apr 2021 21:02:57 -0400 Subject: [PATCH] backup r321 --- CommandLines.cpp | 17 +- CommandLines.h | 6 +- Makefile | 3 +- Overlaps.cpp | 325 +++++++++++++++------ Overlaps.h | 3 +- Process_Read.h | 1 + Purge_Dups.cpp | 490 +++++++++++++------------------ Purge_Dups.h | 29 +- hic.cpp | 238 ++++++++++++++- partig.cpp | 732 +++++++++++++++++++++++++++++++++++++++++++++++ partig.h | 51 ++++ 11 files changed, 1494 insertions(+), 401 deletions(-) create mode 100644 partig.cpp create mode 100644 partig.h diff --git a/CommandLines.cpp b/CommandLines.cpp index 5586fa8..4d9105c 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -32,6 +32,8 @@ static ko_longopt_t long_options[] = { { "b-cov", ko_required_argument, 317 }, { "h-cov", ko_required_argument, 318 }, { "m-rate", ko_required_argument, 319 }, + { "b-partition", ko_no_argument, 320 }, + { "t-occ", ko_required_argument, 321 }, { 0, 0, 0 } }; @@ -87,6 +89,9 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " -d INT upper bound of the binned k-mer's frequency [%d]\n", asm_opt->mid_cnt); fprintf(stderr, " -3 FILE list of hap1/paternal read names []\n"); fprintf(stderr, " -4 FILE list of hap2/maternal read names []\n"); + fprintf(stderr, " --t-occ INT\n"); + fprintf(stderr, " force remove unitigs with >INT unexpected haplotype-specific reads;\n"); + fprintf(stderr, " ignore graph topology; [%d]\n", asm_opt->trio_flag_occ_thres); fprintf(stderr, " Purge-dups:\n"); fprintf(stderr, " -l INT purge level. 0: no purging; 1: light; 2/3: aggressive [0 for trio; 2 for unzip]\n"); @@ -96,6 +101,7 @@ void Print_H(hifiasm_opt_t* asm_opt) asm_opt->purge_overlap_len); fprintf(stderr, " --purge-cov INT\n"); fprintf(stderr, " coverage upper bound of Purge-dups [auto]\n"); + ///fprintf(stderr, " --high-het enable this mode for high heterozygosity sample [experimental, not stable]\n"); fprintf(stderr, " Hi-C-partition [experimental, not stable]:\n"); @@ -157,6 +163,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->recover_atg_cov_min = -1024; asm_opt->recover_atg_cov_max = INT_MAX; asm_opt->hom_global_coverage = -1; + asm_opt->hom_global_coverage_set = 0; asm_opt->bed_inconsist_rate = 70; asm_opt->hic_inconsist_rate = 30; ///asm_opt->bub_mer_length = 3; @@ -165,6 +172,8 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->b_high_cov = -1; asm_opt->m_rate = 0.75; asm_opt->hap_occ = 1; + asm_opt->polyploidy = 2; + asm_opt->trio_flag_occ_thres = 60; } void destory_enzyme(enzyme* f) @@ -616,7 +625,11 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 306) asm_opt->max_ov_diff_final = atof(opt.arg); else if (c == 307) asm_opt->extract_list = opt.arg; else if (c == 308) asm_opt->extract_iter = atoi(opt.arg); - else if (c == 309) asm_opt->hom_global_coverage = atoi(opt.arg); + else if (c == 309) + { + asm_opt->hom_global_coverage = atoi(opt.arg); + asm_opt->hom_global_coverage_set = 1; + } else if (c == 310) { char* s = NULL; @@ -636,6 +649,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) else if (c == 317) asm_opt->b_low_cov = atoi(opt.arg); else if (c == 318) asm_opt->b_high_cov = atoi(opt.arg); else if (c == 319) asm_opt->m_rate = atof(opt.arg); + else if (c == 320) asm_opt->flag |= HA_F_PARTITION; + else if (c == 321) asm_opt->trio_flag_occ_thres = atoi(opt.arg); else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg); diff --git a/CommandLines.h b/CommandLines.h index ece883d..1cb9909 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.14.2-r317" +#define HA_VERSION "0.14.2-r320" #define VERBOSE 0 @@ -18,6 +18,7 @@ #define HA_F_BAN_POST_JOIN 0x100 #define HA_F_BAN_ASSEMBLY 0x200 #define HA_F_HIGH_HET 0x400 +#define HA_F_PARTITION 0x800 #define HA_MIN_OV_DIFF 0.02 // min sequence divergence in an overlap @@ -74,6 +75,7 @@ typedef struct { int recover_atg_cov_min; int recover_atg_cov_max; int hom_global_coverage; + int hom_global_coverage_set; int bed_inconsist_rate; int hic_inconsist_rate; @@ -93,6 +95,8 @@ typedef struct { long long mem_buf; long long coverage; int hap_occ; + int polyploidy; + int trio_flag_occ_thres; } hifiasm_opt_t; diff --git a/Makefile b/Makefile index 2349e50..3f25f03 100644 --- a/Makefile +++ b/Makefile @@ -6,7 +6,7 @@ CPPFLAGS= INCLUDES= OBJS= CommandLines.o Process_Read.o Assembly.o Hash_Table.o \ POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \ - htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o hic.o + htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o hic.o partig.o EXE= hifiasm LIBS= -lz -lpthread -lm @@ -72,3 +72,4 @@ main.o: Levenshtein_distance.h htab.h sketch.o: kvec.h htab.h Process_Read.h Overlaps.h kdq.h CommandLines.h sys.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h hic.o: hic.h +partig.o: partig.h diff --git a/Overlaps.cpp b/Overlaps.cpp index 64b6ab5..7bcca74 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -45,7 +45,8 @@ KSORT_INIT_GENERIC(uint32_t) ///this value has been updated at the first line of build_string_graph_without_clean long long min_thres; - +uint32_t print_untig_by_read(ma_ug_t *g, const char* name, uint32_t in, ma_hit_t_alloc* sources, +ma_hit_t_alloc* reverse_sources, const char* info); void init_bub_label_t(bub_label_t* x, uint32_t n_thres, uint32_t n_reads) { @@ -11711,61 +11712,6 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) kv_destroy(new_rtg_edges.a); } -void classify_untigs_debug(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut, -ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, -kvec_asg_arc_t_warp* new_rtg_edges, int max_hang, int min_ovlp, trans_chain* t_ch) -{ - uint64_t i, dip_thre_max, dip_thres, n_utg; - uint8_t* primary_flag = (uint8_t*)calloc(sg->n_seq, sizeof(uint8_t)); - hap_cov_t *cov = init_hap_cov_t(ug, sg, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp, 0); - - int is_set = ((asm_opt.hom_global_coverage != -1)? 1 : 0); - - if(is_set == 0) - { - purge_dups(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, - asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, 0, 0, 1, cov, 0); - dip_thre_max = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE); - asm_opt.hom_global_coverage = -1; - } - else - { - dip_thre_max = asm_opt.hom_global_coverage; - } - dip_thre_max *= 0.75; - - for (i = 0; i < ug->g->n_seq; i++) - { - dip_thres = dip_thre_max; - if(ug->u.a[i].n <= dip_thre_max) dip_thres = dip_thre_max * 0.9; - - if(get_ug_coverage(&ug->u.a[i], sg, coverage_cut, sources, ruIndex, primary_flag)g->seq[i].c = 1; - if(get_ug_coverage_aggressive(ug, i, sg, coverage_cut, sources, ruIndex, primary_flag, &n_utg) dip_thres*2) ug->g->seq[i].c = 2; - } - } - else - { - ug->g->seq[i].c = 0; - } - - uint32_t flag = get_unitig_het_arb(&ug->u.a[i], t_ch->is_r_het, 2, 1, 0); - if(flag != ug->g->seq[i].c) - { - if(flag == 0 || ug->g->seq[i].c == 0) fprintf(stderr, "Attention\n"); - fprintf(stderr, "utg%.6lul, flag: %u, ug->g->seq[i].c: %u\n", i+1, flag, ug->g->seq[i].c); - } - } - free(primary_flag); - destory_hap_cov_t(&cov); - fprintf(stderr, "[M::%s] diploid coverage threshold: %lu\n", __func__, dip_thre_max); - exit(0); -} - - void set_ug_coverage_aggressive(ma_ug_t *ug, uint32_t uID, asg_t* read_g, const ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag, trans_chain* t_ch, long long het_cov_thres) @@ -12203,6 +12149,180 @@ bub_label_t* b_mask_t) 0.05, 0.9, max_hang, min_ovlp, 0, b_mask_t); } + +void set_trio_flag_by_cov(ma_ug_t *ug, hap_cov_t *cov) +{ + kvec_t(uint64_t) idx; kv_init(idx); + uint32_t i, k, j, qn, tn, s[2], flag; + ma_utg_t *u = NULL, *w = NULL; + hc_links link; + init_hc_links(&link, ug->g->n_seq, cov->t_ch); + for (i = 0, idx.n = 0; i < link.a.n; i++) + { + qn = i; + u = &(ug->u.a[qn]); + for (k = 0; k < u->n; k++) + { + if((R_INF.trio_flag[u->a[k]>>33]&SET_TRIO)==0) break; + } + + if(k >= u->n) continue; ///whole unitig is primary + + for (k = 0, s[0] = s[1] = 0; k < link.a.a[qn].f.n; k++) + { + tn = link.a.a[qn].f.a[k].uID; + w = &(ug->u.a[tn]); + for (j = 0; j < w->n; j++) + { + if((R_INF.trio_flag[w->a[j]>>33]&FATHER)||(R_INF.trio_flag[w->a[j]>>33]&MOTHER)) + { + s[0]++; + } + + if(R_INF.trio_flag[w->a[j]>>33]&SET_TRIO) + { + s[1]++; + } + } + } + + if(s[1] > 0) + { + kv_push(uint64_t, idx, (uint64_t)((uint32_t)-1 - s[0]) << 32 | (qn)); + } + } + + for (i = 0; i < idx.n; i++) + { + qn = (uint32_t)idx.a[i]; + u = &(ug->u.a[qn]); + for (k = 0, s[0] = s[1] = 0; k < link.a.a[qn].f.n; k++) + { + tn = link.a.a[qn].f.a[k].uID; + w = &(ug->u.a[tn]); + for (j = 0; j < w->n; j++) + { + if(R_INF.trio_flag[w->a[j]>>33]&FATHER) + { + s[0]++; + } + + if(R_INF.trio_flag[w->a[j]>>33]&MOTHER) + { + s[1]++; + } + } + } + + if(s[0] >= s[1]) + { + flag = MOTHER; + } + else + { + flag = FATHER; + } + for (k = 0; k < u->n; k++) + { + if(R_INF.trio_flag[u->a[k]>>33]&SET_TRIO) continue; + if(cov->t_ch->is_r_het[u->a[k]>>33] == N_HET) continue; + R_INF.trio_flag[u->a[k]>>33] |= flag; + } + } + + + for (i = 0, idx.n = 0; i < link.a.n; i++) + { + qn = i; + u = &(ug->u.a[qn]); + for (k = 0; k < u->n; k++) + { + if(R_INF.trio_flag[u->a[k]>>33]&FATHER) + { + R_INF.trio_flag[u->a[k]>>33] = FATHER; + } + else if(R_INF.trio_flag[u->a[k]>>33]&MOTHER) + { + R_INF.trio_flag[u->a[k]>>33] = MOTHER; + } + else + { + R_INF.trio_flag[u->a[k]>>33] = AMBIGU; + } + } + } + kv_destroy(idx); + destory_hc_links(&link); +} + +void print_r_het(hap_cov_t *cov, uint8_t* trio_flag, const char* cmd) +{ + if(cov && cov->t_ch) + { + fprintf(stderr, "\n+%s-is_r_het[1369536]=%u\n", cmd, cov->t_ch->is_r_het[1369536]); + fprintf(stderr, "+%s-is_r_het[5097804]=%u\n", cmd, cov->t_ch->is_r_het[5097804]); + fprintf(stderr, "+%s-is_r_het[603738]=%u\n", cmd, cov->t_ch->is_r_het[603738]); + } + if(trio_flag) + { + fprintf(stderr, "+%s-trio_flag[1369536]=%u\n", cmd, trio_flag[1369536]); + fprintf(stderr, "+%s-trio_flag[5097804]=%u\n", cmd, trio_flag[5097804]); + fprintf(stderr, "+%s-trio_flag[603738]=%u\n", cmd, trio_flag[603738]); + } +} + +void output_bp_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, +ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, +long long tipsLen, float tip_drop_ratio, long long stops_threshold, +R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, +bub_label_t* b_mask_t) +{ + hic_clean(sg); + + kvec_asg_arc_t_warp new_rtg_edges; + kv_init(new_rtg_edges.a); + ma_ug_t *ug = NULL; + ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); + + + hap_cov_t *cov = NULL; + asg_t *copy_sg = copy_read_graph(sg); + ma_ug_t *copy_ug = copy_untig_graph(ug); + ///asm_opt.purge_overlap_len = asm_opt.purge_overlap_len_hic; + ///asm_opt.purge_simi_thres = asm_opt.purge_simi_rate_hic; + adjust_utg_by_primary(©_ug, copy_sg, TRIO_THRES, sources, reverse_sources, coverage_cut, + tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio, + max_hang, min_ovlp, &new_rtg_edges, &cov, b_mask_t, 1); + + // print_untig_by_read(copy_ug, "m64011_190830_220126/175638789/ccs", 1369536, NULL, NULL, "sa"); + // print_untig_by_read(copy_ug, "m64012_190921_234837/21039588/ccs", 5097804, NULL, NULL, "sa"); + // print_untig_by_read(copy_ug, "m64011_190830_220126/88867583/ccs", 603738, NULL, NULL, "sa"); + + ma_ug_destroy(copy_ug); + asg_destroy(copy_sg); + + // print_untig_by_read(copy_ug, "m64011_190830_220126/175638789/ccs", 1369536, NULL, NULL, "sb"); + // print_untig_by_read(copy_ug, "m64012_190921_234837/21039588/ccs", 5097804, NULL, NULL, "sb"); + // print_untig_by_read(copy_ug, "m64011_190830_220126/88867583/ccs", 603738, NULL, NULL, "sb"); + + // print_r_het(cov, R_INF.trio_flag, "out-0"); + + set_trio_flag_by_cov(ug, cov); + + // print_r_het(cov, R_INF.trio_flag, "out-1"); + + + destory_hap_cov_t(&cov); + ma_ug_destroy(ug); + kv_destroy(new_rtg_edges.a); + + output_unitig_graph(sg, coverage_cut, output_file_name, sources, ruIndex, max_hang, min_ovlp); + output_trio_unitig_graph(sg, coverage_cut, output_file_name, FATHER, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, + 0.05, 0.9, max_hang, min_ovlp, 0, b_mask_t); + output_trio_unitig_graph(sg, coverage_cut, output_file_name, MOTHER, sources, reverse_sources, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, + 0.05, 0.9, max_hang, min_ovlp, 0, b_mask_t); +} + ma_ug_t* merge_utg(ma_ug_t **dest, ma_ug_t **src) { asg_t *g_d = (*dest)->g, *g_s = (*src)->g; @@ -13039,7 +13159,8 @@ R_to_U* ruIndex, uint32_t min_edge_length, float drop_ratio, uint32_t stops_thre } } - if(cov && operation != CUT) collect_trans_cov(__func__, &b_0, &b_1, ug, read_sg, cov); + ///note: we need to remove b_0, insetad of b_1 here + if(cov && operation != CUT) collect_trans_cov(__func__, &b_1, &b_0, ug, read_sg, cov); break; } @@ -13388,8 +13509,6 @@ R_to_U* ruIndex, uint32_t positive_flag, float drop_rate) } } - - del_node++; } @@ -13517,7 +13636,8 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex, uint32_ asg_seq_drop(g, b.b.a[k]>>1); } - if(cov) collect_trans_cov(__func__, &b_0, &b_1, ug, read_sg, cov); + ///note: we need to remove b_0, insetad of b_1 here + if(cov) collect_trans_cov(__func__, &b_1, &b_0, ug, read_sg, cov); ///lable the primary one b_0.b.n = 0; @@ -13960,7 +14080,7 @@ float drop_ratio, uint32_t trio_flag, float trio_drop_rate, hap_cov_t *cov) resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, trio_flag, drop_ratio); drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex); - all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, DOUBLE_CHECK_THRES); + all_to_all_deduplicate(ug, read_g, coverage_cut, sources, trio_flag, trio_drop_rate, reverse_sources, ruIndex, DOUBLE_CHECK_THRES, asm_opt.trio_flag_occ_thres); if(is_first) { is_first = 0; @@ -13986,6 +14106,8 @@ void print_graph_statistic(asg_t *g, const char* cmd) fprintf(stderr, "%s->n_node: %lu, n_arc: %lu, size: %lu\n", cmd, n_node, n_arc, size); } + + void clean_primary_untig_graph(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t* coverage_cut, long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, buf_t* b_0, uint8_t* visit, float density, @@ -13997,18 +14119,16 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov) int round = T_ROUND; ///uint32_t is_first = 1; - redo: + redo: ///print_graph_statistic(g, "beg"); ///print_debug_gfa(read_g, ug, coverage_cut, "debug_trans_ovlp_hg002", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); asg_pop_bubble_primary_trio(ug, NULL, (uint32_t)-1, DROP, cov, /**is_first?0:1**/1); ///untig_asg_arc_simple_large_bubbles_trio(ug, read_g, reverse_sources, 2, ruIndex, (uint32_t)-1, DROP, cov); ///is_first = 0; - if(just_bubble_pop == 0) { cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 2); } - long long pre_cons = get_graph_statistic(g); long long cur_cons = 0; while(pre_cons != cur_cons) @@ -14038,7 +14158,6 @@ int just_bubble_pop, float drop_ratio, hap_cov_t *cov) cut_trio_tip_primary(g, ug, tipsLen, (uint32_t)-1, 0, read_g, reverse_sources, ruIndex, 2); } - resolve_tangles(ug, read_g, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, (uint32_t)-1, drop_ratio); drop_semi_circle(ug, g, read_g, reverse_sources, ruIndex); unitig_arc_del_short_diploid_by_length_topo(g, ug, drop_ratio, asm_opt.max_short_tip, reverse_sources, 0, 1); @@ -14075,8 +14194,6 @@ void set_drop_trio_flag(ma_ug_t *ug) } } - - void update_unitig_graph(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate) @@ -14109,6 +14226,8 @@ uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate) if(R_INF.trio_flag[rId] == flag) flag_occ++; if(R_INF.trio_flag[rId] != flag) non_flag_occ++; } + + ///if(is_final_check && v == 0) fprintf(stderr, "flag: %u, flag_occ: %u, non_flag_occ: %u\n", flag, flag_occ, non_flag_occ); if(is_final_check == 0 && hap_label_occ == u->n) continue; ///if(is_double_check && non_flag_occ < u->n*DOUBLE_CHECK_THRES) continue; @@ -14169,7 +14288,6 @@ uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate) } - void get_candidate_uids(asg_t* nsg, ma_utg_t* nsu, kvec_t_u64_warp* u_vecs, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex) { @@ -14290,13 +14408,13 @@ uint32_t* non_require, uint32_t* ambigious) ///note: to use this function, don't renew unitig graph!!!!!!!!! void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, -ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate) +ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate, int non_tig_occ) { kvec_t_u64_warp u_vecs; kv_init(u_vecs.a); - uint32_t n_vtx, v, k, is_Unitig, uId, rId, convex, flag_occ, non_flag_occ, ambigious, /**is_ambigious,**/ flag, n_reduce = 1; + uint32_t n_vtx, v, k, is_Unitig, is_tig, uId, rId, convex, flag_occ, non_flag_occ, ambigious, /**is_ambigious,**/ flag, n_reduce = 1; asg_t* nsg = NULL; ma_utg_t* nsu = NULL; nsg = ug->g; @@ -14332,11 +14450,16 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate) if(nsu->m == 0) continue; if(nsg->seq[uId].del) continue; if(nsg->seq[uId].c == ALTER_LABLE) continue; + is_tig = 1; if(get_real_length(nsg, v^1, NULL) == 1) { get_real_length(nsg, v^1, &convex); - if(get_real_length(nsg, convex^1, NULL) == 1) continue; + if(get_real_length(nsg, convex^1, NULL) == 1) + { + is_tig = 0; + ///continue; + } } get_unitig_trio_flag(nsu, flag, &flag_occ, &non_flag_occ, &ambigious); @@ -14349,6 +14472,7 @@ ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate) else**/ { ///we may need it or not + if(is_tig == 0 && (int)flag_occ <= non_tig_occ) continue; if(flag_occ <= ((non_flag_occ+flag_occ)*drop_rate)) continue; if((flag_occ+non_flag_occ) == 0) continue; if(nsu->n >= 100) @@ -14698,7 +14822,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) uint32_t v, n_vtx = nsg->n_seq; hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp, asm_opt.purge_level_trio>0?1:0); - if(asm_opt.purge_level_primary > 0) + if(cov->t_ch) { set_r_het_flag(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, max_hang, min_ovlp, cov->t_ch); @@ -14742,6 +14866,7 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) clean_trio_untig_graph(*ug, read_g, coverage_cut, sources, reverse_sources, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, chimeric_rate, 0, 0, drop_ratio, flag, drop_rate, cov); + ///delete_useless_nodes(ug); delete_useless_trio_nodes(ug, read_g, coverage_cut, sources, ruIndex); @@ -14750,6 +14875,8 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 0, DOUBLE_CHECK_THRES, flag, drop_rate); + ///if(flag == MOTHER) print_debug_gfa(read_g, *ug, coverage_cut, "debug_trio_1", sources, ruIndex, asm_opt.max_hang_Len, asm_opt.min_overlap_Len); + renew_utg(ug, read_g, new_rtg_edges); @@ -14767,7 +14894,8 @@ kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t) } ///if(flag == MOTHER) print_untig_by_read(*ug, "m64043_200627_000137/124716590/ccs", 2789716, NULL, NULL, "beg"); - + + update_unitig_graph((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, 1, FINAL_DOUBLE_CHECK_THRES, flag, drop_rate); @@ -22502,8 +22630,8 @@ uint32_t collect_p_trans) uint32_t v, n_vtx = nsg->n_seq, k, rId, just_contain; ma_utg_t* u = NULL; hap_cov_t *cov = init_hap_cov_t(*ug, read_g, sources, ruIndex, reverse_sources, - coverage_cut, max_hang, min_ovlp, asm_opt.purge_level_primary>0?1:0); - if(asm_opt.purge_level_primary > 0) + coverage_cut, max_hang, min_ovlp, (asm_opt.purge_level_primary>0||i_cov)?1:0); + if(cov->t_ch) { set_r_het_flag(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, max_hang, min_ovlp, cov->t_ch); @@ -22519,13 +22647,13 @@ uint32_t collect_p_trans) nsg->seq[v].c = PRIMARY_LABLE; EvaluateLen((*ug)->u, v) = (*ug)->u.a[v].n; } + clean_primary_untig_graph(*ug, read_g, sources, reverse_sources, coverage_cut, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, NULL, NULL, 0, 0, 0, chimeric_rate, 0, 0, drop_ratio, cov); delete_useless_nodes(ug); renew_utg(ug, read_g, new_rtg_edges); if(i_cov) goto skip_purge; - ///classify_untigs_debug(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, max_hang, min_ovlp, cov->t_ch); if(asm_opt.purge_level_primary > 0) { @@ -22620,16 +22748,13 @@ uint32_t collect_p_trans) recover_utg_by_coverage(ug, read_g, coverage_cut, sources, ruIndex, cov->t_ch); if(i_cov) { - if(asm_opt.purge_level_primary > 0) + if(cov->t_ch && collect_p_trans) { - if(collect_p_trans) - { - // purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, - // asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, - // 0, 0, cov, collect_p_trans); - } + purge_dups(*ug, read_g, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, + asm_opt.purge_simi_thres, asm_opt.purge_overlap_len, max_hang, min_ovlp, drop_ratio, + 0, 0, cov, collect_p_trans); } - + (*i_cov) = cov; } else @@ -26624,7 +26749,6 @@ bub_label_t* b_mask_t) asg_destroy(copy_sg); copy_sg = NULL; - uint32_t beg_idx, occ; bubble_type bub; memset(&bub, 0, sizeof(bubble_type)); @@ -26634,15 +26758,12 @@ bub_label_t* b_mask_t) rescue_bubbles_by_contained_reads(ug, sg, sources, coverage_cut, ruIndex, max_hang, min_ovlp, chainLenThres, beg_idx, occ, &bub, b_mask_t); ///output_unitig_graph(sg, coverage_cut, (char*)"debug_1.rescue", sources, ruIndex, max_hang, min_ovlp); - - ma_ug_destroy(ug); ug = NULL; ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); reset_bub(&bub, ug, cov->t_ch, &new_rtg_edges); beg_idx = bub.f_bub; occ = bub.b_bub + bub.b_end_bub + bub.tangle_bub; rescue_bubbles_by_missing_ovlp(ug, sg, sources, coverage_cut, ruIndex, max_hang, min_ovlp, chainLenThres, beg_idx, occ, &bub, b_mask_t); ///output_unitig_graph(sg, coverage_cut, (char*)"debug_2.hic", sources, ruIndex, max_hang, min_ovlp); - ma_ug_destroy(ug); ug = NULL; ug = ma_ug_gen_primary(sg, PRIMARY_LABLE); reset_bub(&bub, ug, cov->t_ch, &new_rtg_edges); beg_idx = bub.f_bub; occ = bub.b_bub + bub.b_end_bub + bub.tangle_bub; @@ -26655,7 +26776,6 @@ bub_label_t* b_mask_t) reset_bub(&bub, ug, cov->t_ch, &new_rtg_edges); rescue_missing_hap_ovlp(ug, sg, sources, coverage_cut, max_hang, min_ovlp, &bub, gap_fuzz); } - destory_bubbles(&bub); @@ -28490,6 +28610,22 @@ int max_hang, int min_ovlp) } } +void set_hom_global_coverage(hifiasm_opt_t *opt, asg_t *sg, ma_sub_t* coverage_cut, +ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, int max_hang, int min_ovlp) +{ + if(opt->hom_global_coverage_set == 0) + { + ma_ug_t *ug = NULL; + ug = ma_ug_gen(sg); + + hap_cov_t *cov = init_hap_cov_t(ug, sg, sources, ruIndex, reverse_sources, coverage_cut, max_hang, min_ovlp, 0); + purge_dups(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, NULL, + opt->purge_simi_thres, opt->purge_overlap_len, max_hang, min_ovlp, 0, 0, 1, cov, 0); + destory_hap_cov_t(&cov); + + ma_ug_destroy(ug); + } +} void clean_graph( int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, @@ -28694,6 +28830,9 @@ ma_sub_t **coverage_cut_ptr, int debug_g) /*******************************for debug***************************************/ } + // set_hom_global_coverage(&asm_opt, sg, coverage_cut, sources, reverse_sources, ruIndex, + // max_hang_length, mini_overlap_length); + if (ha_opt_triobin(&asm_opt) && ha_opt_hic(&asm_opt)) { @@ -28725,6 +28864,14 @@ ma_sub_t **coverage_cut_ptr, int debug_g) 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, &b_mask_t); free(buf); } + else if(asm_opt.flag & HA_F_PARTITION) + { + char *buf = (char*)calloc(strlen(output_file_name) + 25, 1); + sprintf(buf, "%s.bp", output_file_name); + output_bp_graph(sg, coverage_cut, buf, sources, reverse_sources, (asm_opt.max_short_tip*2), + 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length, &b_mask_t); + free(buf); + } else { output_unitig_graph(sg, coverage_cut, output_file_name, sources, ruIndex, max_hang_length, mini_overlap_length); diff --git a/Overlaps.h b/Overlaps.h index b4d2c64..b65fbc0 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -1013,7 +1013,7 @@ uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, kvec_t_u32_warp* void rescue_missing_overlaps_aggressive(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t is_bubble_check, uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, bub_label_t* b_mask_t); void all_to_all_deduplicate(ma_ug_t* ug, asg_t* read_g, ma_sub_t* coverage_cut, -ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate); +ma_hit_t_alloc* sources, uint8_t postive_flag, float drop_rate, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, float double_check_rate, int non_tig_occ); void drop_semi_circle(ma_ug_t *ug, asg_t* nsg, asg_t* read_g, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex); void rescue_wrong_overlaps_to_unitigs(ma_ug_t *i_ug, asg_t *r_g, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, R_to_U* ruIndex, int max_hang, int min_ovlp, long long bubble_dist, kvec_asg_arc_t_warp* keep_edges, bub_label_t* b_mask_t); @@ -1056,6 +1056,7 @@ typedef struct{ #define N_HET 0 #define C_HET 1 #define P_HET 2 +#define S_HET 4 typedef struct{ kvec_t(uint32_t) uIDs; diff --git a/Process_Read.h b/Process_Read.h index 4b27d81..eeb5641 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -101,6 +101,7 @@ typedef struct #define MIX_TRIO 3 #define NON_TRIO 4 #define DROP 5 +#define SET_TRIO 8 typedef struct { diff --git a/Purge_Dups.cpp b/Purge_Dups.cpp index eca0e76..3371da8 100644 --- a/Purge_Dups.cpp +++ b/Purge_Dups.cpp @@ -8,6 +8,7 @@ #include "kthread.h" #include "kdq.h" #include "hic.h" +#include "partig.h" KDQ_INIT(uint64_t) @@ -49,33 +50,6 @@ typedef struct { #define MIXED 3 #define FLIP 4 -typedef struct { - uint8_t rev; - uint8_t type; - uint8_t status; - uint32_t x_beg_pos; - uint32_t x_end_pos; - uint32_t y_beg_pos; - uint32_t y_end_pos; - uint32_t x_beg_id; - uint32_t x_end_id; - uint32_t y_beg_id; - uint32_t y_end_id; - uint32_t xUid; - uint32_t yUid; - uint32_t weight; - long long score; -}hap_overlaps; - -typedef struct { - kvec_t(hap_overlaps) a; -}kvec_hap_overlaps; - -typedef struct { - kvec_hap_overlaps* x; - uint32_t num; -}hap_overlaps_list; - typedef struct { uint64_t* vote_counting; @@ -2377,7 +2351,7 @@ ma_hit_t_alloc* reverse_sources, long long xBegPos, long long xEndPos) // fprintf(stderr, "tailIndex->a.n: %u, xReads->n: %u, total_match: %lld, hap_match: %lld, inp_match: %lld\n", // tailIndex->a.n, xReads->n, (xEndPos - xBegPos + 1), hap_match, inp_match); - return (inp_match*3) - ((hap_match-inp_match)*1); + return ((double)(inp_match)*1) - ((double)(hap_match-inp_match)*0.334); } @@ -3139,7 +3113,7 @@ static void hap_alignment_advance_worker(void *_data, long eid, int tid) } -int inline get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn) +int get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn) { uint32_t i; for (i = 0; i < x->a.n; i++) @@ -3324,6 +3298,10 @@ double filter_rate) { homLen += ovlp; } + else if(asm_opt.polyploidy <= 2 && (a[k].h_status&S_HET)) + { + homLen += ovlp; + } else { hetLen += ovlp; @@ -3354,10 +3332,14 @@ double filter_rate) { homLen += ovlp; } + else if(asm_opt.polyploidy <= 2 && (a[k].h_status&S_HET)) + { + homLen += ovlp; + } else { hetLen += ovlp; - } + } } if(hetLen <= ((hetLen + homLen) * filter_rate)) @@ -4774,242 +4756,6 @@ void remove_contained_haplotig(hap_overlaps_list* all_ovlp, ma_ug_t *ug, asg_t* // } } -#define HAP1_LABLE 7 -#define HAP2_LABLE 8 -void init_contig_phase(hap_overlaps_list* all_ovlp, asg_t *purge_g, float drop_ratio) -{ - - uint32_t v, i, n_vtx = purge_g->n_seq * 2, beg, end, uId, yId, need_update; - long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen, hap1_weight, hap2_weight; - buf_t b_0; - memset(&b_0, 0, sizeof(buf_t)); - clean_purge_graph(purge_g, drop_ratio, 1); - - for (v = 0; v < n_vtx; ++v) - { - if(purge_g->seq[v>>1].c == ALTER_LABLE || purge_g->seq[v>>1].del) continue; - if(get_real_length(purge_g, v, NULL) != 1) continue; - if(get_real_length(purge_g, v^1, NULL) != 0) continue; - - beg = v; - b_0.b.n = 0; - if(get_unitig(purge_g, NULL, beg, &end, &nodeLen, &baseLen, &max_stop_nodeLen, - &max_stop_baseLen, 1, &b_0) == LOOP) - { - continue; - } - - for (i = 0; i < b_0.b.n; i++) - { - uId = b_0.b.a[i]; - if((i>>1) == 0) purge_g->seq[uId].c = HAP1_LABLE; - else purge_g->seq[uId].c = HAP2_LABLE; - } - } - - need_update = 1; - - while (need_update) - { - need_update = 0; - for (v = 0; v < all_ovlp->num; v++) - { - uId = v; - if(purge_g->seq[uId].c == HAP1_LABLE || purge_g->seq[uId].c == HAP2_LABLE) - { - continue; - } - if(all_ovlp->x[uId].a.n == 0) - { - purge_g->seq[uId].c = HAP1_LABLE; - continue; - } - - - hap1_weight = hap2_weight = 0; - for (i = 0; i < all_ovlp->x[uId].a.n; i++) - { - yId = all_ovlp->x[uId].a.a[i].yUid; - if(purge_g->seq[yId].c == HAP1_LABLE) hap1_weight += all_ovlp->x[uId].a.a[i].weight; - if(purge_g->seq[yId].c == HAP2_LABLE) hap2_weight += all_ovlp->x[uId].a.a[i].weight; - } - - if(hap1_weight == 0 && hap2_weight == 0) - { - need_update = 1; - continue; - } - - if(hap1_weight >= hap2_weight) - { - purge_g->seq[uId].c = HAP1_LABLE; - } - else - { - purge_g->seq[uId].c = HAP2_LABLE; - } - } - } - - free(b_0.b.a); -} - -void partition_contigs(hap_overlaps_list* all_ovlp, ma_ug_t *ug, asg_t *purge_g, hap_cov_t *cov, double keep_rate, -int max_hang, int min_ovlp, float drop_ratio) -{ - int r, index; - long long max_score; - uint32_t v, i, uId, xUid, is_contain = 0, m; - hap_overlaps *p = NULL; - asg_arc_t t, *p_t; - for (v = 0; v < all_ovlp->num; v++) - { - uId = v; - for (i = 0; i < all_ovlp->x[uId].a.n; i++) - { - if(all_ovlp->x[uId].a.a[i].status == DELETE) continue; - all_ovlp->x[uId].a.a[i].status = MIXED; - } - } - - for (v = 0; v < all_ovlp->num; v++) - { - uId = v; p = NULL; is_contain = 0; - if(all_ovlp->x[uId].a.n == 0) continue; - for (i = 0; i < all_ovlp->x[uId].a.n; i++) - { - if(p == NULL || p->score < all_ovlp->x[uId].a.a[i].score) - { - p = &(all_ovlp->x[uId].a.a[i]); - } - } - - max_score = p->score; - for (i = 0; i < all_ovlp->x[uId].a.n; i++) - { - if(all_ovlp->x[uId].a.a[i].type == YCX) - { - if(!filter_secondary_chain(p->score, all_ovlp->x[uId].a.a[i].score, MAX(0.95, keep_rate))) - { - continue; - } - - xUid = all_ovlp->x[uId].a.a[i].xUid; - purge_g->seq[xUid].c = ALTER_LABLE; - purge_g->seq[xUid].del = 1; - ///collect_trans_purge_cov(cov, ug, &(all_ovlp->x[uId].a.a[i]), 0); - - if(is_contain == 0 || all_ovlp->x[uId].a.a[i].score > max_score) - { - max_score = all_ovlp->x[uId].a.a[i].score; - } - is_contain = 1; - } - } - - if(is_contain) - { - for (i = 0; i < all_ovlp->x[uId].a.n; i++) - { - if(filter_secondary_chain(max_score, all_ovlp->x[uId].a.a[i].score, keep_rate)) - { - all_ovlp->x[uId].a.a[i].status = FLIP; - index = get_specific_hap_overlap(&(all_ovlp->x[all_ovlp->x[uId].a.a[i].yUid]), - all_ovlp->x[uId].a.a[i].yUid, all_ovlp->x[uId].a.a[i].xUid); - if(index == -1) fprintf(stderr, "ERROR 5\n"); - all_ovlp->x[all_ovlp->x[uId].a.a[i].yUid].a.a[index].status = FLIP; - - } - } - } - } - - for (v = 0; v < all_ovlp->num; v++) - { - uId = v; - ///has been removed as contained - if(purge_g->seq[uId].del || purge_g->seq[uId].c == ALTER_LABLE) continue; - for (i = 0; i < all_ovlp->x[uId].a.n; i++) - { - if(all_ovlp->x[uId].a.a[i].status == DELETE) continue; - - if(purge_g->seq[all_ovlp->x[uId].a.a[i].xUid].c == ALTER_LABLE|| - purge_g->seq[all_ovlp->x[uId].a.a[i].xUid].del|| - purge_g->seq[all_ovlp->x[uId].a.a[i].yUid].c == ALTER_LABLE|| - purge_g->seq[all_ovlp->x[uId].a.a[i].yUid].del) - { - continue; - } - - - ///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i])); - - r = get_hap_arch(&(all_ovlp->x[uId].a.a[i]), ug->u.a[all_ovlp->x[uId].a.a[i].xUid].len, - ug->u.a[all_ovlp->x[uId].a.a[i].yUid].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); - - if(r < 0) continue; - p_t = asg_arc_pushp(purge_g); - *p_t = t; - } - } - - asg_cleanup(purge_g); - asg_symm(purge_g); - - clean_purge_graph(purge_g, keep_rate, 0); - - asg_arc_t *av = NULL; - uint32_t n_vtx = purge_g->n_seq*2, nv, w; - for (v = 0; v < n_vtx; ++v) - { - if(purge_g->seq[v>>1].c == ALTER_LABLE || purge_g->seq[v>>1].del) continue; - nv = asg_arc_n(purge_g, v); - av = asg_arc_a(purge_g, v); - if (nv == 0) continue; - for (i = 0; i < nv; ++i) - { - if (av[i].del) continue; - w = av[i].v; - index = get_specific_hap_overlap(&(all_ovlp->x[v>>1]), v>>1, w>>1); - all_ovlp->x[v>>1].a.a[index].status = FLIP; - } - } - - p = NULL; - for (v = 0; v < all_ovlp->num; v++) - { - uId = v; - for (i = m = 0; i < all_ovlp->x[uId].a.n; i++) - { - if(all_ovlp->x[uId].a.a[i].status != FLIP) continue; - all_ovlp->x[uId].a.a[m] = all_ovlp->x[uId].a.a[i]; - if(p == NULL || p->score > all_ovlp->x[uId].a.a[m].score) - { - p = &(all_ovlp->x[uId].a.a[m]); - } - m++; - } - all_ovlp->x[uId].a.n = m; - } - - max_score = 0; - if(p && p->score <= 0) - { - max_score = ((p->score)*-1) + 1; - for (v = 0; v < all_ovlp->num; v++) - { - uId = v; - for (i = 0; i < all_ovlp->x[uId].a.n; i++) - { - all_ovlp->x[uId].a.a[i].score += max_score; - } - } - } - - init_contig_phase(all_ovlp, purge_g, drop_ratio); - -} - void purge_dups_back(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain, @@ -5243,7 +4989,7 @@ void debug_p_g_t(p_g_t* pg, hap_cov_t *cov, asg_t *read_g) offset += (uint32_t)u->a[i]; if(i >= sid && i <= eid) { - if((!!cov->t_ch->is_r_het[u->a[i]>>33]) != t->h_status) + if(cov->t_ch->is_r_het[u->a[i]>>33] != t->h_status) { fprintf(stderr, "ERROR-(-3): is_r_het: %u, h_status: %u\n", cov->t_ch->is_r_het[u->a[i]>>33], t->h_status); } @@ -5390,15 +5136,15 @@ p_g_t *init_p_g_t_back(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g) p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g) { - uint32_t v, uId, k, l, offset, l_pos, g_beg_idx, occ, ovlp, tLen, zLen; + uint32_t v, uId, k, l, offset, l_pos, g_beg_idx, occ/**, ovlp, tLen, zLen**/; p_g_t *pg = NULL; CALLOC(pg, 1); pg->ug = ug; asg_t* nsg = pg->ug->g; ma_utg_t *u = NULL; - p_node_t *t = NULL, *z = NULL; - asg_arc_t *e = NULL; + p_node_t *t = NULL/**, *z = NULL**/; + ///asg_arc_t *e = NULL; p_g_in_t *x = NULL; - pg->pg_het = asg_init(); + ///pg->pg_het = asg_init(); pg->pg_h_lev = asg_init(); kv_init(pg->pg_het_node); kv_init(pg->pg_h_lev_idx); @@ -5417,6 +5163,14 @@ p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g) pg->pg_h_lev->seq[uId].c = PRIMARY_LABLE; } + // if(asm_opt.polyploidy <= 2) + // { + // for (v = 0; v < cov->t_ch->r_num; v++) + // { + // if(cov->t_ch->is_r_het[v]&P_HET) cov->t_ch->is_r_het[v] |= S_HET; + // } + // } + for (v = 0; v < nsg->n_seq; v++) { uId = v; @@ -5427,18 +5181,18 @@ p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g) ///fprintf(stderr, "\n+v: %u, pg->pg_het_node.n: %u\n", v, (uint32_t)pg->pg_het_node.n); for (k = 1, l = 0, offset = 0, l_pos = 0; k <= u->n; ++k) { - ///if(k < u->n && cov->t_ch->is_r_het[u->a[k]>>33] == 0) fprintf(stderr, "sbsbs-uId=%u\n", uId); - if (k == u->n || (!!cov->t_ch->is_r_het[u->a[k]>>33]) != (!!cov->t_ch->is_r_het[u->a[l]>>33])) + ///if (k == u->n || (!!cov->t_ch->is_r_het[u->a[k]>>33]) != (!!cov->t_ch->is_r_het[u->a[l]>>33])) + if (k == u->n || cov->t_ch->is_r_het[u->a[k]>>33] != cov->t_ch->is_r_het[u->a[l]>>33]) { kv_pushp(p_node_t, pg->pg_het_node, &t); t->c_ug_id = uId; - t->h_status = (!!cov->t_ch->is_r_het[u->a[l]>>33]); + t->h_status = cov->t_ch->is_r_het[u->a[l]>>33]; t->baseBeg = l_pos; t->baseEnd = offset + read_g->seq[u->a[k-1]>>33].len - 1; t->nodeBeg = l; t->nodeEnd = k - 1; ///if(t->b_ug_id == (uint32_t)-1) fprintf(stderr, "xxxx\n"); - asg_seq_set(pg->pg_het, pg->pg_het_node.n-1, t->baseEnd+1-t->baseBeg, 0); + ///asg_seq_set(pg->pg_het, pg->pg_het_node.n-1, t->baseEnd+1-t->baseBeg, 0); ///fprintf(stderr, "l: %u, k: %u, u->n: %u, t->h_status: %u\n", l, k, u->n, t->h_status); l = k; l_pos = offset + (uint32_t)u->a[k-1]; @@ -5451,29 +5205,29 @@ p_g_t *init_p_g_t(ma_ug_t *ug, hap_cov_t *cov, asg_t *read_g) x->beg = g_beg_idx; x->occ = occ; ///fprintf(stderr, "-v: %u, pg->pg_het_node.n: %u\n", v, (uint32_t)pg->pg_het_node.n); - if(occ > 1) - { - for (k = g_beg_idx; (k + 1) < pg->pg_het_node.n; ++k) - { - t = &(pg->pg_het_node.a[k]); tLen = t->baseEnd + 1 - t->baseBeg; - z = &(pg->pg_het_node.a[k+1]); zLen = z->baseEnd + 1 - z->baseBeg; + // if(occ > 1) + // { + // for (k = g_beg_idx; (k + 1) < pg->pg_het_node.n; ++k) + // { + // t = &(pg->pg_het_node.a[k]); tLen = t->baseEnd + 1 - t->baseBeg; + // z = &(pg->pg_het_node.a[k+1]); zLen = z->baseEnd + 1 - z->baseBeg; - ovlp = ((MIN(t->baseEnd, z->baseEnd) >= MAX(t->baseBeg, z->baseBeg))? - MIN(t->baseEnd, z->baseEnd) - MAX(t->baseBeg, z->baseBeg) + 1 : 0); + // ovlp = ((MIN(t->baseEnd, z->baseEnd) >= MAX(t->baseBeg, z->baseBeg))? + // MIN(t->baseEnd, z->baseEnd) - MAX(t->baseBeg, z->baseBeg) + 1 : 0); - e = asg_arc_pushp(pg->pg_het); - e->ol = ovlp; - e->ul = (k<<1); e->ul <<= 32; e->ul += (tLen - ovlp); - e->v = ((k+1)<<1); e->del = 0; e->el = e->no_l_indel = e->strong = 1; + // e = asg_arc_pushp(pg->pg_het); + // e->ol = ovlp; + // e->ul = (k<<1); e->ul <<= 32; e->ul += (tLen - ovlp); + // e->v = ((k+1)<<1); e->del = 0; e->el = e->no_l_indel = e->strong = 1; - e = asg_arc_pushp(pg->pg_het); - e->ol = ovlp; - e->ul = ((k+1)<<1)+1; e->ul <<= 32; e->ul += (zLen - ovlp); - e->v = (k<<1)+1; e->del = 0; e->el = e->no_l_indel = e->strong = 1; - } - } + // e = asg_arc_pushp(pg->pg_het); + // e->ol = ovlp; + // e->ul = ((k+1)<<1)+1; e->ul <<= 32; e->ul += (zLen - ovlp); + // e->v = (k<<1)+1; e->del = 0; e->el = e->no_l_indel = e->strong = 1; + // } + // } } - asg_cleanup(pg->pg_het); + ///asg_cleanup(pg->pg_het); ///debug_p_g_t(pg, cov, read_g); ///print_p_g_t_interval(pg, cov); @@ -5493,6 +5247,147 @@ void destory_p_g_t(p_g_t **pg) } } + +void partition_contigs(hap_overlaps_list* all_ovlp, ma_ug_t *ug, hap_cov_t *cov, double keep_rate, +int max_hang, int min_ovlp, float drop_ratio, p_g_t *pg) +{ + int r, index; + uint32_t v, i, uId, m; + hap_overlaps *p = NULL; + asg_arc_t t, *p_t = NULL; + + for (v = 0; v < all_ovlp->num; v++) + { + uId = v; + ///has been removed as contained + if(pg->pg_h_lev->seq[uId].del || pg->pg_h_lev->seq[uId].c == ALTER_LABLE) continue; + for (i = 0; i < all_ovlp->x[uId].a.n; i++) + { + if(all_ovlp->x[uId].a.a[i].type == YCX) continue; + if(all_ovlp->x[uId].a.a[i].type == XCY) continue; + /****************************may have bugs********************************/ + if(all_ovlp->x[uId].a.a[i].score <= 0) continue; + /****************************may have bugs********************************/ + + if(pg->pg_h_lev->seq[all_ovlp->x[uId].a.a[i].xUid].c == ALTER_LABLE|| + pg->pg_h_lev->seq[all_ovlp->x[uId].a.a[i].xUid].del|| + pg->pg_h_lev->seq[all_ovlp->x[uId].a.a[i].yUid].c == ALTER_LABLE|| + pg->pg_h_lev->seq[all_ovlp->x[uId].a.a[i].yUid].del) + { + continue; + } + + + ///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i])); + + r = get_hap_arch(&(all_ovlp->x[uId].a.a[i]), ug->u.a[all_ovlp->x[uId].a.a[i].xUid].len, + ug->u.a[all_ovlp->x[uId].a.a[i].yUid].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + + if(r < 0) continue; + p_t = asg_arc_pushp(pg->pg_h_lev); + *p_t = t; + } + } + + asg_cleanup(pg->pg_h_lev); + asg_symm(pg->pg_h_lev); + + clean_purge_graph(pg->pg_h_lev, keep_rate, 0); + + asg_arc_t *av = NULL; + uint32_t n_vtx = pg->pg_h_lev->n_seq<<1, nv, a, b; + for (v = 0; v < n_vtx; v++) + { + av = asg_arc_a(pg->pg_h_lev, v); + nv = asg_arc_n(pg->pg_h_lev, v); + + for (i = 0; i < nv; ++i) + { + if(av[i].del) continue; + a = av[i].ul>>33; + b = av[i].v>>1; + + index = get_specific_hap_overlap(&(all_ovlp->x[a]), a, b); + p = &(all_ovlp->x[a].a.a[index]); + p->status = MIXED; + + index = get_specific_hap_overlap(&(all_ovlp->x[b]), b, a); + p = &(all_ovlp->x[b].a.a[index]); + p->status = MIXED; + } + } + + for (v = 0; v < all_ovlp->num; v++) + { + uId = v; + for (i = m = 0; i < all_ovlp->x[uId].a.n; i++) + { + p = (&all_ovlp->x[uId].a.a[i]); + /****************************may have bugs********************************/ + if(p->score <= 0) continue; + /****************************may have bugs********************************/ + if(p->type == YCX || p->type == XCY || p->status == MIXED) + { + all_ovlp->x[uId].a.a[m] = (*p); + all_ovlp->x[uId].a.a[m].status = SELF_EXIST; + m++; + } + } + all_ovlp->x[uId].a.n = m; + } + + for (v = 0; v < all_ovlp->num; v++) + { + uId = v; p = NULL; + if(all_ovlp->x[uId].a.n == 0) continue; + for (i = 0; i < all_ovlp->x[uId].a.n; i++) + { + if(p == NULL || p->score < all_ovlp->x[uId].a.a[i].score) + { + p = &(all_ovlp->x[uId].a.a[i]); + } + } + if(!p) continue; + + for (i = m = 0; i < all_ovlp->x[uId].a.n; i++) + { + if(all_ovlp->x[uId].a.a[i].type == YCX || all_ovlp->x[uId].a.a[i].type == XCY) + { + if(filter_secondary_chain(p->score, all_ovlp->x[uId].a.a[i].score, keep_rate)) + { + all_ovlp->x[uId].a.a[m] = all_ovlp->x[uId].a.a[i]; + m++; + } + } + else + { + all_ovlp->x[uId].a.a[m] = all_ovlp->x[uId].a.a[i]; + m++; + } + + } + all_ovlp->x[uId].a.n = m; + } + + + for (v = 0; v < all_ovlp->num; v++) + { + uId = v; + for (i = m = 0; i < all_ovlp->x[uId].a.n; i++) + { + p = (&all_ovlp->x[uId].a.a[i]); + index = get_specific_hap_overlap(&(all_ovlp->x[p->yUid]), p->yUid, p->xUid); + if(index != -1) + { + all_ovlp->x[uId].a.a[m] = (*p); + m++; + } + } + all_ovlp->x[uId].a.n = m; + } + +} + void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain, @@ -5573,6 +5468,13 @@ uint32_t just_coverage, hap_cov_t *cov, uint32_t collect_p_trans) ///if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp); filter_hap_overlaps_by_length(&all_ovlp, purege_minLen); + if(collect_p_trans) + { + pt_solve(&all_ovlp, cov->t_ch, ug, read_g, 0.8, R_INF.trio_flag); + goto end_coverage; + } + + pg = init_p_g_t(ug, cov, read_g); ///normalize_hap_overlaps(&all_ovlp, &back_all_ovlp); diff --git a/Purge_Dups.h b/Purge_Dups.h index 102376f..b3f15b2 100644 --- a/Purge_Dups.h +++ b/Purge_Dups.h @@ -13,6 +13,33 @@ #define REAL_ALTER_THRES 0.1 #define CHAIN_FILTER_RATE 0.7 +typedef struct { + uint8_t rev; + uint8_t type; + uint8_t status; + uint32_t x_beg_pos; + uint32_t x_end_pos; + uint32_t y_beg_pos; + uint32_t y_end_pos; + uint32_t x_beg_id; + uint32_t x_end_id; + uint32_t y_beg_id; + uint32_t y_end_id; + uint32_t xUid; + uint32_t yUid; + uint32_t weight; + long long score; +}hap_overlaps; + +typedef struct { + kvec_t(hap_overlaps) a; +}kvec_hap_overlaps; + +typedef struct { + kvec_hap_overlaps* x; + uint32_t num; +}hap_overlaps_list; + void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t purege_minLen, int max_hang, int min_ovlp, float drop_ratio, uint32_t just_contain, @@ -25,6 +52,6 @@ hap_cov_t* init_hap_cov_t(ma_ug_t *ug, asg_t* read_g, ma_hit_t_alloc* sources, R ma_sub_t *coverage_cut, int max_hang, int min_ovlp, uint32_t is_collect_trans); void destory_hap_cov_t(hap_cov_t **x); void chain_trans_ovlp(hap_cov_t *cov, ma_ug_t *ug, asg_t *read_sg, buf_t* xReads, uint32_t targetBaseLen, uint32_t* xEnd); - +int get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn); #endif \ No newline at end of file diff --git a/hic.cpp b/hic.cpp index 3c62887..18c1afb 100644 --- a/hic.cpp +++ b/hic.cpp @@ -19,7 +19,7 @@ KDQ_INIT(uint64_t) #define OFFSET_RATE 0.000000001 #define OFFSET_SECOND_RATE 0.0000000001 #define SCALL 10000 -#define OFFSET_RATE_MAX_W 20.8286263517*SCALL +#define OFFSET_RATE_MAX_W 6.90675477865*SCALL #define OFFSET_RATE_MIN_W 4.0000003e-10*SCALL #define HIC_COUNTER_BITS 12 @@ -4143,8 +4143,8 @@ void print_hc_links(hc_links* link, int dir, H_partition* hap) } } - fprintf(stderr, "self-utg%.6dl\tFW:%f\tRW:%f\tRT:%f\n**************************************************\n", - (int)(i+1), f_w, r_w, (f_w+r_w) != 0? r_w/(f_w+r_w):0); + // fprintf(stderr, "self-utg%.6dl\tFW:%f\tRW:%f\tRT:%f\n**************************************************\n", + // (int)(i+1), f_w, r_w, (f_w+r_w) != 0? r_w/(f_w+r_w):0); } } @@ -7437,6 +7437,7 @@ void update_bsg(asg_t *bsg, kvec_asg_arc_t_warp* edges) } void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub) { + double index_time = yak_realtime(); ma_ug_t *bub_ug = bub->b_ug; asg_t *bsg = bub->b_g; uint32_t k, i, v, n_vx, new_bub; @@ -7605,11 +7606,13 @@ void resolve_bubble_chain_tangle(ma_ug_t* ug, bubble_type* bub) kv_destroy(res_btg.a); kv_destroy(res_utg.a); kv_destroy(edges.a); ///print_debug_bubble_graph(bub, ug, asm_opt.output_file_name); + fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time); } void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint32_t is_end) { + double index_time = yak_realtime(); if(bub->b_ug) ma_ug_destroy(bub->b_ug); if(bub->chain_weight.a) kv_destroy(bub->chain_weight); kvec_t_u32_warp broken; @@ -7828,7 +7831,7 @@ void update_bubble_chain(ma_ug_t* ug, bubble_type* bub, uint32_t is_middle, uint } } **/ - + fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time); } void set_b_utg_weight_flag(bubble_type* bub, buf_t* b, uint32_t v, uint8_t* vis_flag, uint32_t flag, uint32_t* occ) @@ -8187,6 +8190,7 @@ int double_check_bub_branch(asg_arc_t *t, ma_ug_t *bs_ug, double *e_w, uint32_t void clean_bubble_chain_by_HiC(ma_ug_t* ug, hc_links* link, bubble_type* bub) { + double index_time = yak_realtime(); ma_ug_t *bs_ug = bub->b_ug; uint32_t v, u, i, m, max_i, nv, rv, n_vx, root, flag_pri = 1, flag_aux = 2, flag_ava = 4, occ; double w, cutoff = 2/**, max_w_cutoff = MAX(MIN(100*OFFSET_RATE_MIN_W, OFFSET_RATE_MAX_W/100), OFFSET_RATE_MIN_W)**/; @@ -8383,11 +8387,13 @@ void clean_bubble_chain_by_HiC(ma_ug_t* ug, hc_links* link, bubble_type* bub) free(vis); free(is_vis); free(is_used); free(dedup); free(b.b.a); free(e_w); free(e_occ); kv_destroy(stack.a); kv_destroy(result.a); kv_destroy(res_utg.a); kv_destroy(edges.a); ma_ug_destroy(back_bs_ug); + fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time); } void append_boundary_chain(ma_ug_t* ug, hc_links* link, bubble_type* bub) { + double index_time = yak_realtime(); ma_ug_t *bs_ug = bub->b_ug; uint32_t v, u, i, k, beg_idx, m, nv, n_vx, flag_pri = 1, flag_aux = 2, flag_ava = 4; uint32_t root, rId_0, ori_0, root_0, new_bub; @@ -8538,6 +8544,7 @@ void append_boundary_chain(ma_ug_t* ug, hc_links* link, bubble_type* bub) // } // } /*******************************for debug************************************/ + fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time); } int cmp_chain_hic_w_weight(const void * a, const void * b) @@ -9408,12 +9415,12 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, return 0; } - // for (i = 0; i < dis->n; i++) - // { - // if(i > 0 && dis->a[i].beg != dis->a[i-1].end) fprintf(stderr, "ERROR: dis->a[i].beg: %lu, dis->a[i-1].end: %lu\n", dis->a[i].beg, dis->a[i-1].end); - // fprintf(stderr, "beg: %lu, end: %lu, cnt_0: %lu, cnt_1: %lu, error_rate: %f\n", - // dis->a[i].beg, dis->a[i].end, dis->a[i].cnt_0, dis->a[i].cnt_1, (double)(dis->a[i].cnt_1)/(double)(dis->a[i].cnt_1 + dis->a[i].cnt_0)); - // } + for (i = 0; i < dis->n; i++) + { + if(i > 0 && dis->a[i].beg != dis->a[i-1].end) fprintf(stderr, "ERROR: dis->a[i].beg: %lu, dis->a[i-1].end: %lu\n", dis->a[i].beg, dis->a[i-1].end); + fprintf(stderr, "beg: %lu, end: %lu, cnt_0: %lu, cnt_1: %lu, error_rate: %f\n", + dis->a[i].beg, dis->a[i].end, dis->a[i].cnt_0, dis->a[i].cnt_1, (double)(dis->a[i].cnt_1)/(double)(dis->a[i].cnt_1 + dis->a[i].cnt_0)); + } LeastSquare_advance(dis, idx, med); // fprintf(stderr, "idx->a: %f, idx->b: %f, idx->frac: %f, med: %lu\n", @@ -9441,6 +9448,196 @@ int get_trans_rate_function(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, } +int get_trans_rate_function_advance(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, MT* M, H_partition* hap, trans_idx* dis) +{ + kvec_t(uint64_t) buf; + kv_init(buf); + uint64_t beg, end, cnt[2]; + uint64_t k, i, t_d, r_idx, f_idx, med = (uint64_t)-1; + int beg_status, end_status; + for (i = 0; i < link->a.n; i++) + { + for (k = 0; k < link->a.a[i].e.n; k++) + { + link->a.a[i].e.a[k].dis = (uint64_t)-1; + } + } + fill_utg_distance_multi(idx, link, M, bub); + + + buf.n = 0; + for (k = 0; k < hits->a.n; ++k) + { + beg = ((hits->a.a[k].s<<1)>>(64 - idx->uID_bits)); + end = ((hits->a.a[k].e<<1)>>(64 - idx->uID_bits)); + + if(IF_HOM(beg, *bub)) continue; + if(IF_HOM(end, *bub)) continue; + + + t_d = get_hic_distance(&(hits->a.a[k]), link, idx); + if(t_d == (uint64_t)-1) continue; + if(beg == end) + { + t_d = (t_d << 1); + } + else + { + beg_status = get_phase_status(hap, beg); + if(beg_status != 1 && beg_status != -1) continue; + end_status = get_phase_status(hap, end); + if(end_status != 1 && end_status != -1) continue; + if(beg_status != end_status) + { + t_d = (t_d << 1) + 1; + } + else + { + t_d = (t_d << 1); + } + } + + kv_push(uint64_t, buf, t_d); + } + + ///might have bias, we may not use right linkage larger than trans rc linkage + radix_sort_hc64(buf.a, buf.a+buf.n); + + for (k = 0, r_idx = f_idx = (uint64_t)-1; k < buf.n; k++) + { + if((buf.a[k]&1) == 0) r_idx = k; + if((buf.a[k]&1) == 1) f_idx = k; + } + buf.n = MIN(r_idx, f_idx); + + trans_p_t* p = NULL; + dis->n = 0; + uint64_t bin_size = MIN(2250, buf.n>>8), m; + i = 0; + while (i < buf.n) + { + kv_pushp(trans_p_t, *dis, &p); + p->beg = i; + + cnt[0] = cnt[1] = 0; + k = MIN(i+bin_size, buf.n); + for (; i < k; i++) + { + cnt[buf.a[i]&1]++; + } + p->end = i; + p->cnt_0 = cnt[0]; + p->cnt_1 = cnt[1]; + } + + i = m = 0; med = (uint64_t)-1; + while(i < dis->n) + { + if(dis->a[i].cnt_0 > 0 && dis->a[i].cnt_1 > 0) + { + dis->a[m] = dis->a[i]; + m++; + i++; + continue; + } + + if(med == (uint64_t)-1) med = buf.a[dis->a[i].beg]>>1; + k = i; cnt[0] = cnt[1] = 0; + for (; i < dis->n; i++) + { + cnt[0] += dis->a[i].cnt_0; + cnt[1] += dis->a[i].cnt_1; + if(cnt[0] > 0 && cnt[1] > 0) break; + } + + if(i < dis->n) + { + dis->a[m].cnt_0 = cnt[0]; + dis->a[m].cnt_1 = cnt[1]; + dis->a[m].beg = dis->a[k].beg; + dis->a[m].end = dis->a[i].end; + m++; + i++; + continue; + } + + cnt[0] -= dis->a[k].cnt_0; + cnt[1] -= dis->a[k].cnt_1; + while (1) + { + cnt[0] += dis->a[k].cnt_0; + cnt[1] += dis->a[k].cnt_1; + if(cnt[0] > 0 && cnt[1] > 0) break; + + if(k == 0) + { + k = (uint64_t)-1; + break; + } + k--; + } + + if(k != (uint64_t)-1) + { + dis->a[m].cnt_0 = cnt[0]; + dis->a[m].cnt_1 = cnt[1]; + dis->a[m].end = dis->a[i-1].end; + continue; + } + m = 0; + break; + } + dis->n = m; + if(dis->n == 0) + { + kv_destroy(buf); + return 0; + } + + for (i = 0; i < dis->n; i++) + { + dis->a[i].beg = buf.a[dis->a[i].beg]>>1; + dis->a[i].end = (buf.a[dis->a[i].end-1]>>1); + } + + for (i = 0; i < dis->n - 1; i++) + { + dis->a[i].end += ((dis->a[i+1].beg - dis->a[i].end)/2); + dis->a[i+1].beg = dis->a[i].end; + } + + for (i = 0; i < dis->n; i++) + { + if(i > 0 && dis->a[i].beg != dis->a[i-1].end) fprintf(stderr, "ERROR: dis->a[i].beg: %lu, dis->a[i-1].end: %lu\n", dis->a[i].beg, dis->a[i-1].end); + fprintf(stderr, "beg: %lu, end: %lu, cnt_0: %lu, cnt_1: %lu, error_rate: %f\n", + dis->a[i].beg, dis->a[i].end, dis->a[i].cnt_0, dis->a[i].cnt_1, (double)(dis->a[i].cnt_1)/(double)(dis->a[i].cnt_1 + dis->a[i].cnt_0)); + } + + LeastSquare_advance(dis, idx, med); + // fprintf(stderr, "idx->a: %f, idx->b: %f, idx->frac: %f, med: %lu\n", + // (double)idx->a, (double)idx->b, (double)idx->frac, med); + + dis->max = dis->a[dis->n-1].end; + + + kv_destroy(buf); + if(idx->a < 0) idx->a = 0; + if(idx->a == 0) + { + idx->b = MAX((((double)(dis->a[dis->n-1].cnt_1))/((double)(dis->a[dis->n-1].cnt_0 + dis->a[dis->n-1].cnt_1))), idx->b); + } + if(idx->b < 0 && get_trans(idx, dis->max) < 0) + { + idx->b = ((double)(dis->a[dis->n-1].cnt_1))/((double)(dis->a[dis->n-1].cnt_0 + dis->a[dis->n-1].cnt_1)); + } + + // fprintf(stderr, "idx->a: %f, idx->b: %f, idx->frac: %f, med: %lu\n", + // (double)idx->a, (double)idx->b, (double)idx->frac, med); + + return 1; +} + + void init_hic_p(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, kvec_hc_edge* back_hc_edge, MT* M, H_partition* hap, uint32_t ignore_dis) { @@ -9450,7 +9647,7 @@ kvec_hc_edge* back_hc_edge, MT* M, H_partition* hap, uint32_t ignore_dis) if(bub->round_id > 0 && ignore_dis == 0) { - is_comples_weight = get_trans_rate_function(idx, hits, link, bub, M, hap, &dis); + is_comples_weight = get_trans_rate_function_advance(idx, hits, link, bub, M, hap, &dis); } @@ -10304,6 +10501,7 @@ void print_contig_partition(H_partition* hap, const char* debug) void adjust_contig_partition(H_partition* hap, hc_links* link) { + double index_time = yak_realtime(); uint32_t i, k, *h0, h0_n, *h1, h1_n; uint32_t h0_status[4], h1_status[4], h0_status_max; int h0_h, h1_h; @@ -10404,6 +10602,7 @@ void adjust_contig_partition(H_partition* hap, hc_links* link) hap->g_p.a[i].weight[1] = get_cluster_weight(hap, link, h1, h1_n); hap->g_p.a[i].weight_convex = get_cluster_inner_weight(hap, link, h0, h0_n, h1, h1_n); } + fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time); } uint32_t get_weightest_uid(uint32_t* a, uint32_t n, H_partition* hap, uint8_t* hap_label_flag, @@ -11733,6 +11932,7 @@ int graph_bipartiteness(uint32_t* b_a, uint32_t b_a_n, uint8_t *color, hc_links* void assign_per_unitig_G_partition(G_partition* g_p, uint64_t hap_n, hc_links* link, bubble_type* bub, uint32_t bubble_first) { + double index_time = yak_realtime(); reset_G_partition(g_p, hap_n); partition_warp* res = NULL; @@ -11872,7 +12072,7 @@ uint32_t bubble_first) g_p->index[res->a.a[k]] = (g_p->index[res->a.a[k]] << 1) + 1; } } - + fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time); } typedef struct { @@ -12093,6 +12293,7 @@ bub_sort_vec* w_stack, uint8_t* vis, uint32_t* set_hap, uint32_t n_utg, uint32_t void sort_bubble_ele(G_partition* g_p, hc_links* link, bubble_type* bub, uint32_t n_utg) { + double index_time = yak_realtime(); kvec_t_u64_warp stack; kv_init(stack.a); bub_sort_vec w_stack; kv_init(w_stack); uint8_t* vis = NULL; MALLOC(vis, n_utg); @@ -12106,6 +12307,7 @@ void sort_bubble_ele(G_partition* g_p, hc_links* link, bubble_type* bub, uint32_ } kv_destroy(stack.a); kv_destroy(w_stack); free(vis); free(set_hap); + fprintf(stderr, "[M::%s::%.3f]\n", __func__, yak_realtime()-index_time); } uint32_t init_contig_partition(H_partition* hap, ha_ug_index* idx, bubble_type* bub, hc_links* link) @@ -13214,22 +13416,32 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx) bub.round_id = 0; bub.n_round = 2; for (bub.round_id = 0; bub.round_id < bub.n_round; bub.round_id++) { + fprintf(stderr, "0********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); identify_bubbles(idx->ug, &bub, idx->cov->t_ch->is_r_het); + fprintf(stderr, "1********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); if(bub.round_id == 0) { + fprintf(stderr, "2********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); collect_hc_links(sl.idx, &sl.hits, &link, &bub, &M); + fprintf(stderr, "3********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); collect_hc_reverse_links(&link, idx->ug, &bub); + fprintf(stderr, "4********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); } init_hic_p((ha_ug_index*)sl.idx, &sl.hits, &link, &bub, &back_hc_edge, &M, &hap, 0); + fprintf(stderr, "5********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); ///init_hic_p_new((ha_ug_index*)sl.idx, &sl.hits, idx->link, &bub, &back_hc_edge, &M); reset_H_partition(&hap, (bub.round_id == 0? 1 : 0)); + fprintf(stderr, "6********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); init_contig_partition(&hap, idx, &bub, &link); + fprintf(stderr, "7********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); phasing_improvement(&hap, &(hap.g_p), idx, &bub, &link); + fprintf(stderr, "8********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); label_unitigs(&(hap.g_p), idx->ug); + fprintf(stderr, "9********bub.round_id: %u(::%.3f)********\n", bub.round_id, yak_realtime()); ///print_hc_links(idx->link, 0, &hap); } - ///print_hc_links(idx->link, 0, &hap); + print_hc_links(&link, 0, &hap); cluster_contigs(&bub, idx, &sl.hits, &M, &hap, &link); diff --git a/partig.cpp b/partig.cpp new file mode 100644 index 0000000..9c90474 --- /dev/null +++ b/partig.cpp @@ -0,0 +1,732 @@ +#define __STDC_LIMIT_MACROS +#include +#include +#include "partig.h" +#include "Purge_Dups.h" +#include "Correct.h" +#include "ksort.h" + +#define generic_key(x) (x) +KRADIX_SORT_INIT(gfa64, uint64_t, generic_key, 8) + +#define pt_a(x, id) ((x).ma.a + ((x).idx.a[(id)]>>32)) +#define pt_n(x, id) ((uint32_t)((x).idx.a[(id)])) + +typedef struct { + int32_t topn; + int32_t n_perturb; + uint64_t seed; + double f_perturb; +} pt_svopt_t; + +typedef struct { + ///uint32_t m, n, *shuffled; + uint32_t *shuffled; + uint32_t off, size; // offset in pt_match_t::cc; size of the component + uint64_t *buf; + kvec_t(uint64_t) a; + int8_t *s, *s_tmp;///s is the status (haplotype) of each unitig: for backup +} solve_aux_t; + +static inline uint64_t kr_splitmix64(uint64_t x) +{ + uint64_t z = (x += 0x9E3779B97F4A7C15ULL); + z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL; + z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL; + return z ^ (z >> 31); +} + +static inline double kr_drand_r(uint64_t *x) +{ + union { uint64_t i; double d; } u; + *x = kr_splitmix64(*x); + u.i = 0x3FFULL << 52 | (*x) >> 12; + return u.d - 1.0; +} + +void pt_svopt_init(pt_svopt_t *opt) +{ + memset(opt, 0, sizeof(pt_svopt_t)); + opt->seed = 11; + opt->topn = 1<<30; + opt->n_perturb = 1000; + opt->f_perturb = 0.1; +} + +static void pt_pdist_idx(pt_match_t *ma) +{ + uint32_t st, i; + kv_resize(uint64_t, ma->idx, ma->n_seq); + ma->idx.n = ma->n_seq; + memset(ma->idx.a, 0, ma->idx.n*sizeof(uint64_t)); + for (st = 0, i = 1; i <= ma->ma.n; ++i) + if (i == ma->ma.n || ma->ma.a[i].sid[0] != ma->ma.a[st].sid[0]) + ma->idx.a[ma->ma.a[st].sid[0]] = (uint64_t)st << 32 | (i - st), st = i; +} + +static pt_match1_t *pt_pdist(const pt_match_t *ma, uint32_t sid1, uint32_t sid2) +{ + pt_match1_t *o = pt_a(*ma, sid1); + uint32_t n = pt_n(*ma, sid1), k; + for (k = 0; k < n; ++k) + if (o[k].sid[1] == sid2) + return &(o[k]); + return NULL; +} + +static void normalize_pdist(pt_match1_t *a, pt_match1_t *b) +{ + if(a->w >= b->w) + { + b->sid[0] = a->sid[1]; + b->sid[1] = a->sid[0]; + b->w = a->w; + } + else + { + a->sid[0] = b->sid[1]; + a->sid[1] = b->sid[0]; + a->w = b->w; + } +} + +uint32_t pt_pdist_symm(pt_match_t *ma) +{ + uint8_t *del = NULL; + uint32_t i, k, n = 0; + pt_match1_t *t = NULL, *m = NULL; + CALLOC(del, ma->ma.n); + + for (i = 0; i < ma->ma.n; ++i) { + m = &ma->ma.a[i]; + if (m->sid[0] == m->sid[1]) + { + del[i] = 1, ++n;///self overlap + continue; + } + t = pt_pdist(ma, m->sid[1], m->sid[0]); + if(!t) + { + del[i] = 1, ++n;///self overlap + continue; + } + normalize_pdist(m, t); + } + + if (n > 0) { + for (i = k = 0; i < ma->ma.n; ++i) + if (!del[i]) ma->ma.a[k++] = ma->ma.a[i]; + ma->ma.n = k; + pt_pdist_idx(ma); + } + + free(del); + return n; +} + + +static void pt_pdist_flt(pt_match_t *ma, uint32_t min_cnt, double drop_thres) +{ + uint32_t i, k, j, n, max, o; + uint8_t *mark = NULL; + CALLOC(mark, ma->ma.n); + for (i = 0; i < ma->n_seq; ++i) { + o = ma->idx.a[i] >> 32; + n = (uint32_t)ma->idx.a[i]; + max = 0; + if (n == 0) continue; + for (j = o; j < o + n; ++j) + max = max > ma->ma.a[j].w? max : ma->ma.a[j].w; + for (j = o; j < o + n; ++j) + if (ma->ma.a[j].w >= (max*drop_thres) || ma->ma.a[j].w + min_cnt >= max) + mark[j] = 1; + } + for (i = 0; i < ma->ma.n; ++i) + { + if (mark[i] == 0) continue; + o = ma->idx.a[ma->ma.a[i].sid[1]]>>32; + n = (uint32_t)ma->idx.a[ma->ma.a[i].sid[1]]; + for (j = o; j < o + n; ++j) + { + if (ma->ma.a[j].sid[1] == ma->ma.a[i].sid[0]) mark[j] = 1; + } + } + + for (i = k = 0; i < ma->ma.n; ++i) + if (mark[i]) ma->ma.a[k++] = ma->ma.a[i]; + ma->ma.n = k; + free(mark); + pt_pdist_idx(ma); + pt_pdist_symm(ma); +} + +pt_match_t *init_pt_match_t(hap_overlaps_list* ha, pt_g_t *x, double f_rate) +{ + pt_match_t *p = NULL; CALLOC(p, 1); p->n_seq = x->ug->g->n_seq; + kv_init(p->idx); kv_init(p->ma); + uint32_t v, i, k, qn, tn, qs, qe, ts, te, occ, as, ae; + uint64_t hetLen, homLen, oLen; + pt_node_t *a = NULL; + pt_match1_t *ma = NULL; + + for (v = 0; v < ha->num; v++) + { + for (i = 0; i < ha->x[v].a.n; i++) + { + if(ha->x[v].a.a[i].score <= 0) continue; + if(ha->x[v].a.a[i].xUid == ha->x[v].a.a[i].yUid) continue; + /*****************qn*****************/ + qn = ha->x[v].a.a[i].xUid; + qs = ha->x[v].a.a[i].x_beg_pos; + qe = ha->x[v].a.a[i].x_end_pos - 1; + + a = x->p.a + x->p_idx.a[qn]; + occ = x->p_idx.a[qn+1] - x->p_idx.a[qn]; + for (k = 0, hetLen = 0, homLen = 0; k < occ; k++) + { + as = a[k].bS; + ae = a[k].bE; + oLen = ((MIN(qe, ae) >= MAX(qs, as))? MIN(qe, ae) - MAX(qs, as) + 1 : 0); + if(homLen + hetLen > 0 && oLen == 0) break; + if(oLen == 0) continue; + if(a[k].hs == N_HET) + { + homLen += oLen; + } + else if(asm_opt.polyploidy <= 2 && (a[k].hs&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].hs&S_HET)) + { + homLen += oLen; + } + else + { + hetLen += oLen; + } + } + + if(hetLen <= ((hetLen + homLen)*f_rate)) continue; + /*****************qn*****************/ + + + /*****************tn*****************/ + tn = ha->x[v].a.a[i].yUid; + ts = ha->x[v].a.a[i].y_beg_pos; + te = ha->x[v].a.a[i].y_end_pos - 1; + + a = x->p.a + x->p_idx.a[tn]; + occ = x->p_idx.a[tn+1] - x->p_idx.a[tn]; + for (k = 0, hetLen = 0, homLen = 0; k < occ; k++) + { + as = a[k].bS; + ae = a[k].bE; + oLen = ((MIN(te, ae) >= MAX(ts, as))? MIN(te, ae) - MAX(ts, as) + 1 : 0); + if(homLen + hetLen > 0 && oLen == 0) break; + if(oLen == 0) continue; + if(a[k].hs == N_HET) + { + homLen += oLen; + } + else if(asm_opt.polyploidy <= 2 && (a[k].hs&P_HET))///if(asm_opt.polyploidy <= 2 && (a[k].hs&S_HET)) + { + homLen += oLen; + } + else + { + hetLen += oLen; + } + } + + if(hetLen <= ((hetLen + homLen)*f_rate)) continue; + /*****************tn*****************/ + + kv_pushp(pt_match1_t, p->ma, &ma); + ma->sid[0] = ha->x[v].a.a[i].xUid; + ma->sid[1] = ha->x[v].a.a[i].yUid; + ma->w = ha->x[v].a.a[i].score; + } + } + + pt_pdist_idx(p); + pt_pdist_symm(p); + pt_pdist_flt(p, 5, 0.5); + return p; +} + + +void debug_pt_g_t(pt_g_t *pg) +{ + fprintf(stderr, "0----------[M::%s]----------\n", __func__); + uint32_t i, offset, v, sid, eid, spos, epos, p_status, p_uid, occ; + ma_utg_t *u = NULL; + pt_node_t *a = NULL, *t = NULL; + + for (v = 0; v < pg->ug->u.n; v++) + { + a = pg->p.a + pg->p_idx.a[v]; + occ = pg->p_idx.a[v+1] - pg->p_idx.a[v]; + for (i = 0; i < occ; i++) + { + if(a[i].uID != v) fprintf(stderr, "ERROR-s\n"); + } + } + + for (v = 0, p_status = (uint32_t)-1, p_uid = (uint32_t)-1; v < pg->p.n; v++) + { + t = &(pg->p.a[v]); + sid = t->nS; + eid = t->nE; + spos = t->bS; + epos = t->bE; + if(p_uid == t->uID && p_status == t->hs) + { + fprintf(stderr, "ERROR-a\n"); + } + p_status = t->hs; + p_uid = t->uID; + u = &(pg->ug->u.a[t->uID]); + for (i = offset = 0; i < u->n; i++) + { + if(i == sid) + { + if(spos != offset) + { + fprintf(stderr, "ERROR-b\n"); + } + } + + if(i == eid) + { + if(epos != (offset+pg->rg->seq[u->a[i]>>33].len - 1)) + { + fprintf(stderr, "ERROR-c, real end: %u\n", + (uint32_t)(offset+pg->rg->seq[u->a[i]>>33].len - 1)); + } + } + + offset += (uint32_t)u->a[i]; + if(i >= sid && i <= eid) + { + if(pg->t_ch->is_r_het[u->a[i]>>33] != t->hs) + { + fprintf(stderr, "ERROR-d: is_r_het: %u, h_status: %u\n", pg->t_ch->is_r_het[u->a[i]>>33], t->hs); + } + } + } + } + + + pt_match1_t *o = NULL, *s = NULL; + uint32_t k, n, cnt; + for (i = 0; i < pg->e->n_seq; ++i) + { + o = pt_a(*(pg->e), i); n = pt_n(*(pg->e), i); + for (k = 0; k < n; ++k) + { + if(o[k].sid[0] != i) fprintf(stderr, "ERROR-g\n"); + s = pt_pdist(pg->e, o[k].sid[1], o[k].sid[0]); + if(!s) fprintf(stderr, "ERROR-e\n"); + if(s) + { + if(!(s->sid[0] == o[k].sid[1] && s->sid[1] == o[k].sid[0] && s->w == o[k].w)) + { + fprintf(stderr, "ERROR-f\n"); + } + } + } + + for (k = cnt = 0; k < pg->e->ma.n; ++k) + { + if(pg->e->ma.a[k].sid[0] == i) cnt++; + } + + if(cnt != n) fprintf(stderr, "ERROR-h\n"); + } +} + +pt_g_t *init_pt_g_t(hap_overlaps_list* ovlp, trans_chain* t_ch, ma_ug_t *ug, asg_t *read_g, double f_rate) +{ + uint32_t v, l, k, offset, l_pos; + pt_g_t *p = NULL; CALLOC(p, 1); + asg_t* nsg = ug->g; + ma_utg_t *u = NULL; + pt_node_t *t = NULL; + p->ug = ug; + p->rg = read_g; + p->t_ch = t_ch; + kv_init(p->info); p->info.n = p->info.m = p->ug->g->n_seq; CALLOC(p->info.a, p->info.n); + kv_init(p->p); + kv_init(p->p_idx); kv_push(uint32_t, p->p_idx, 0); + + for (v = 0; v < nsg->n_seq; v++) + { + if(nsg->seq[v].del || nsg->seq[v].c == ALTER_LABLE) continue; + + u = &(ug->u.a[v]); + for (k = 1, l = 0, offset = 0, l_pos = 0; k <= u->n; ++k) + { + if (k == u->n || t_ch->is_r_het[u->a[k]>>33] != t_ch->is_r_het[u->a[l]>>33]) + { + kv_pushp(pt_node_t, p->p, &t); + t->uID = v; + t->hs = t_ch->is_r_het[u->a[l]>>33]; + + t->bS = l_pos; + t->bE = offset + read_g->seq[u->a[k-1]>>33].len - 1; + + t->nS = l; + t->nE = k - 1; + l = k; + l_pos = offset + (uint32_t)u->a[k-1]; + } + offset += (uint32_t)u->a[k-1]; + } + kv_push(uint32_t, p->p_idx, p->p.n); + } + + p->e = init_pt_match_t(ovlp, p, f_rate); + return p; +} + +void destory_pt_g_t(pt_g_t **p) +{ + if(!p || !(*p)) return; + kv_destroy((*p)->p); + kv_destroy((*p)->info); + kv_destroy((*p)->p_idx); + kv_destroy((*p)->e->idx); + kv_destroy((*p)->e->ma); + free((*p)->e->cc); + free((*p)->e); + free((*p)); +} + +uint64_t *pt_cc_core(const pt_match_t *ma) +{ + uint32_t i, x, y, *flag; + uint64_t *group; + pt_match1_t *o = NULL; + kvec_t(uint32_t) stack; kv_init(stack); + + MALLOC(flag, ma->n_seq); + for (i = 0; i < ma->n_seq; ++i) + flag[i] = (uint32_t)-1; + + // connected componets + for (i = 0; i < ma->n_seq; ++i) { + if (flag[i] != (uint32_t)-1) continue; + stack.n = 0; + kv_push(uint32_t, stack, i); + while (stack.n > 0) { + uint32_t k, j, n; + stack.n--; + k = stack.a[stack.n]; + flag[k] = i;///group id + // n = (uint32_t)ma->idx[k]; + // s = ma->idx[k] >> 32; + o = pt_a(*ma, k); + n = pt_n(*ma, k); + for (j = 0; j < n; ++j) { + uint32_t t = o[j].sid[1]; + if (flag[t] != (uint32_t)-1) continue; + // if (ns == ms) PT_EXPAND(stack, ms); + // stack[ns++] = t; + kv_push(uint32_t, stack, t); + } + } + } + kv_destroy(stack); + + // precalculate the size of each group + CALLOC(group, ma->n_seq); + for (i = 0; i < ma->n_seq; ++i) + group[i] = (uint64_t)flag[i] << 32 | i; + radix_sort_gfa64(group, group + ma->n_seq); + for (i = 1, x = y = 0; i <= ma->n_seq; ++i) { + if (i == ma->n_seq || group[i]>>32 != group[x]>>32) { + uint32_t j; + for (j = x; j < i; ++j) + group[j] = (uint64_t)y << 32 | (uint32_t)group[j];///(group id)|first element in this group + ++y, x = i; + } + } + free(flag); + return group; +} + +void pt_cc(pt_match_t *ma) +{ + ma->cc = pt_cc_core(ma); +} + +///inspect top INT edges +static int64_t pt_score(const pt_match_t *ma, uint32_t topn, solve_aux_t *aux) +{ + uint32_t i; + int64_t z = 0; + for (i = 0; i < aux->size; ++i) {///aux->size: how many unitigs in this group + uint32_t k = (uint32_t)ma->cc[aux->off + i]; + uint32_t o = ma->idx.a[k] >> 32; + uint32_t n = (uint32_t)ma->idx.a[k], j; + for (j = 0; j < n; ++j) + aux->buf[j] = (uint64_t)((uint32_t)-1 - ma->ma.a[o + j].w) << 32 | (o + j); + radix_sort_gfa64(aux->buf, aux->buf + n); + for (j = 0; j < n && j < topn; ++j) { + const pt_match1_t *m = &ma->ma.a[(uint32_t)aux->buf[j]]; + z += -(int64_t)m->w * aux->s[m->sid[0]] * aux->s[m->sid[1]]; + } + } + return z; +} + +static int64_t pt_solve1_init_phase(const pt_match_t *ma, int32_t topn, uint64_t *x, solve_aux_t *aux) +{ + uint32_t i; + aux->a.n = 0; + for (i = 0; i < aux->size; ++i) { + uint32_t k = (uint32_t)ma->cc[aux->off + i];///unitig id + uint32_t o = ma->idx.a[k] >> 32;///group id + uint32_t n = (uint32_t)ma->idx.a[k], j; + aux->shuffled[i] = k;///init shuffled status + for (j = 0; j < n; ++j) { + const pt_match1_t *m = &ma->ma.a[o + j]; + ///first is weight, second part is edge id + kv_push(uint64_t, aux->a, (uint64_t)((uint32_t)-1 - m->w) << 32 | (o + j)); + } + } + radix_sort_gfa64(aux->a.a, aux->a.a + aux->a.n);///sort all edges in this group + ///randomly assign each unitig + for (i = 0; i < aux->a.n; ++i) { // from the strongest edge to the weakest + const pt_match1_t *m = &ma->ma.a[(uint32_t)aux->a.a[i]]; + ///sid[0]: query id, sid[1]: target id + ///the initial results of aux->s is 0 + if (aux->s[m->sid[0]] == 0 && aux->s[m->sid[1]] == 0) { + *x = kr_splitmix64(*x);// random number + aux->s[m->sid[0]] = *x&1? 1 : -1; + aux->s[m->sid[1]] = -aux->s[m->sid[0]]; + } else if (aux->s[m->sid[0]] == 0) { + aux->s[m->sid[0]] = -aux->s[m->sid[1]]; + } else if (aux->s[m->sid[1]] == 0) { + aux->s[m->sid[1]] = -aux->s[m->sid[0]]; + } + } + return pt_score(ma, topn, aux); +} + +///size is how many unitigs in this group +static void ks_shuffle_uint32_t(size_t n, uint32_t a[], uint64_t *x) +{ + size_t i, j; + for (i = n; i > 1; --i) { + uint32_t tmp; + j = (size_t)(kr_drand_r(x) * i);///semms + tmp = a[j]; a[j] = a[i-1]; a[i-1] = tmp; + } +} + +static void pt_solve1_perturb(const pt_svopt_t *opt, const pt_match_t *ma, uint64_t *x, solve_aux_t *aux) +{ + uint32_t i; + double y; + for (i = 0; i < aux->size; ++i) { + uint32_t k = (uint32_t)ma->cc[aux->off + i]; + y = kr_drand_r(x); + if (y < opt->f_perturb) + aux->s[k] = -aux->s[k]; + } +} + +static int64_t pt_solve1_optimize(const pt_match_t *ma, uint32_t topn, uint64_t *x, solve_aux_t *aux, uint32_t *n_iter) +{ + uint32_t i; + while (1) { + uint32_t n_flip = 0; + ++(*n_iter); + ks_shuffle_uint32_t(aux->size, aux->shuffled, x); + for (i = 0; i < aux->size; ++i) { + uint32_t k = aux->shuffled[i]; + uint32_t o = ma->idx.a[k] >> 32; + uint32_t n = (uint32_t)ma->idx.a[k], j; + uint64_t z[2]; + int8_t s; + for (j = 0; j < n; ++j) { + const pt_match1_t *m = &ma->ma.a[o + j]; + ///assert(m->sid[0] == k); + aux->buf[j] = (uint64_t)((uint32_t)-1 - m->w) << 32 | (o + j); + } + radix_sort_gfa64(aux->buf, aux->buf + n);///still sort by edge weight + for (j = 0, z[0] = z[1] = 0; j < n && j < topn; ++j) { + const pt_match1_t *m = &ma->ma.a[(uint32_t)aux->buf[j]]; + if (aux->s[m->sid[1]] > 0) z[0] += m->w; + else if (aux->s[m->sid[1]] < 0) z[1] += m->w; + } + if (z[0] == z[1]) continue; + s = z[0] > z[1]? -1 : 1; + if (aux->s[k] != s) + aux->s[k] = s, ++n_flip; + } + if (n_flip == 0) break; + } + return pt_score(ma, topn, aux); +} + +uint32_t pt_solve1(const pt_svopt_t *opt, const pt_match_t *ma, uint64_t *x, solve_aux_t *aux) +{ + uint32_t j, k, n_iter = 0; + int64_t sc_ori, sc_opt = -(1<<30), sc; + if (aux->size < 2) return 0;///how many unitigs + + // first guess + ///randomly assign haplotype status, and get a score + sc_ori = pt_solve1_init_phase(ma, opt->topn, x, aux); + if (aux->size == 2) return 0; + + // optimize + sc_opt = pt_solve1_optimize(ma, opt->topn, x, aux, &n_iter); + for (j = 0; j < aux->size; ++j) + aux->s_tmp[aux->shuffled[j]] = aux->s[aux->shuffled[j]]; + for (k = 0; k < (uint32_t)opt->n_perturb; ++k) { + pt_solve1_perturb(opt, ma, x, aux); + sc = pt_solve1_optimize(ma, opt->topn, x, aux, &n_iter); + if (sc > sc_opt) { + for (j = 0; j < aux->size; ++j) + aux->s_tmp[aux->shuffled[j]] = aux->s[aux->shuffled[j]]; + sc_opt = sc; + } else { + for (j = 0; j < aux->size; ++j) + aux->s[aux->shuffled[j]] = aux->s_tmp[aux->shuffled[j]]; + } + } + for (j = 0; j < aux->size; ++j) + aux->s[aux->shuffled[j]] = aux->s_tmp[aux->shuffled[j]]; + fprintf(stderr, "[%s] group:%d, size:%d, #edges:%u, #iter:%d, sc_ori:%ld, sc_opt:%ld\n", __func__, + (uint32_t)(ma->cc[aux->off]>>32), (uint32_t)(aux->size), (uint32_t)(aux->a.n), n_iter, (long)sc_ori, (long)sc_opt); + return n_iter; +} + +int8_t *pt_solve_core(const pt_svopt_t *opt, const pt_match_t *ma) +{ + int8_t *s; + uint32_t st, i, max = 0; + uint64_t x = opt->seed; + solve_aux_t *aux; + CALLOC(aux, 1); CALLOC(aux->s, ma->n_seq); CALLOC(aux->s_tmp, ma->n_seq); + kv_init(aux->a); + for (i = 0; i < ma->n_seq; ++i) {///count how many links for each unitig + uint32_t n = pt_n(*ma, i); + max = max > n? max : n; + } + MALLOC(aux->buf, max); + MALLOC(aux->shuffled, ma->n_seq); // FIXME: this is over-allocation for convenience + for (st = 0, i = 1; i <= ma->n_seq; ++i) { + if (i == ma->n_seq || ma->cc[st]>>32 != ma->cc[i]>>32) {///at same group + if (i - st >= 2) {///all unitigs in the same group + aux->off = st, aux->size = i - st; + pt_solve1(opt, ma, &x, aux); + } + st = i; + } + } + s = aux->s; + kv_destroy(aux->a); free(aux->buf); free(aux->shuffled); free(aux->s_tmp); + free(aux); + return s; +} + +void set_trio_flag(ma_ug_t *ug, asg_t *read_g, uint32_t uID, uint8_t* trio_flag, trans_chain* t_ch, + hap_overlaps_list* ha, pt_match_t *ma, int8_t s) +{ + uint32_t i; + ma_utg_t *u = &(ug->u.a[uID]); + for (i = 0; i < u->n; i++) + { + trio_flag[u->a[i]>>33] |= SET_TRIO; + if(t_ch->is_r_het[u->a[i]>>33] == N_HET) continue; + if(s == 0) + { + if(t_ch->is_r_het[u->a[i]>>33]&P_HET)//special case + { + trio_flag[u->a[i]>>33] |= FATHER; + } + continue; + } + trio_flag[u->a[i]>>33] |= (s > 0? FATHER:MOTHER); + } + + // pt_match1_t *o = pt_a(*ma, uID); + // uint32_t n = pt_n(*ma, uID), k, qn, tn, qs, qe, r_beg, r_end, offset, oLen, found; + // int index; + // for (k = 0; k < n; ++k) + // { + // qn = o[k].sid[0]; tn = o[k].sid[1]; qs = qe = (uint32_t)-1; + // index = get_specific_hap_overlap(&(ha->x[qn]), qn, tn); + // if(index != -1 && ha->x[qn].a.a[index].score == (long long)o[k].w) + // { + // qs = ha->x[qn].a.a[index].x_beg_pos; + // qe = ha->x[qn].a.a[index].x_end_pos - 1; + // } + // else + // { + // index = get_specific_hap_overlap(&(ha->x[tn]), tn, qn); + // if(index != -1 && ha->x[tn].a.a[index].score == (long long)o[k].w) + // { + // qs = ha->x[qn].a.a[index].y_beg_pos; + // qe = ha->x[qn].a.a[index].y_end_pos - 1; + // } + // } + + // if(qs == (uint32_t)-1 || qe == (uint32_t)-1) fprintf(stderr, "ERROR\n"); + // for (i = 0, offset = 0, found = 0; i < u->n; i++) + // { + // r_beg = offset; r_end = offset + (long long)(read_g->seq[u->a[i]>>33].len) - 1; + // offset += (uint32_t)u->a[i]; + // oLen = ((MIN(qe, r_end) >= MAX(qs, r_beg))? MIN(qe, r_end) - MAX(qs, r_beg) + 1 : 0); + // if(oLen > 0) found = 1; + // if(found == 1 && oLen == 0) break; + + // if(oLen > 0 && t_ch->is_r_het[u->a[i]>>33] != N_HET) + // { + // trio_flag[u->a[i]>>33] |= (s > 0? FATHER:MOTHER); + // } + + // } + // } +} + +void pt_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag) +{ + pt_svopt_t opt; + int8_t *s = NULL; + uint64_t *buf = NULL, i; + pt_svopt_init(&opt); + pt_g_t *pg = init_pt_g_t(ovlp, t_ch, ug, read_g, f_rate); + ///debug_pt_g_t(pg); + + pt_cc(pg->e); + s = pt_solve_core(&opt, pg->e); + + + MALLOC(buf, pg->e->ma.n); // FIXME: this is over-allocation for convenience + for (i = 0; i < pg->e->n_seq; ++i) { + uint64_t z[2]; + uint32_t o = pg->e->idx.a[i] >> 32; + uint32_t n = (uint32_t)pg->e->idx.a[i], j; + + set_trio_flag(ug, read_g, i, trio_flag, t_ch, ovlp, pg->e, s[i]); + + pg->info.a[i].s = s[i]; + for (j = 0; j < n; ++j) { + const pt_match1_t *m = &pg->e->ma.a[o + j]; + buf[j] = (uint64_t)((uint32_t)-1 - m->w) << 32 | (o + j); + } + radix_sort_gfa64(buf, buf + n); + for (j = 0, z[0] = z[1] = 0; j < n; ++j) { + const pt_match1_t *m = &pg->e->ma.a[(uint32_t)buf[j]]; + if (s[m->sid[1]] > 0) z[0] += m->w; + else if (s[m->sid[1]] < 0) z[1] += m->w; + } + pg->info.a[i].m[0] = z[0], pg->info.a[i].m[1] = z[1]; + } + + free(buf); + free(s); + destory_pt_g_t(&pg); +} \ No newline at end of file diff --git a/partig.h b/partig.h new file mode 100644 index 0000000..a5d259f --- /dev/null +++ b/partig.h @@ -0,0 +1,51 @@ +#ifndef __PURTIG__ +#define __PURTIG__ +#include +#include +#include "kvec.h" +#include "Overlaps.h" +#include "Purge_Dups.h" + +typedef struct { + uint32_t bS, bE; + uint32_t nS, nE; + uint32_t uID; + uint8_t hs; +}pt_node_t; + +typedef struct { + ///sid[0]: query id + ///sid[1]: target id + uint32_t sid[2]; + uint32_t w; +} pt_match1_t; + +typedef struct { + kvec_t(uint64_t) idx; + kvec_t(pt_match1_t) ma; + uint64_t* cc; + uint32_t n_seq; +} pt_match_t; + +typedef struct { + ///cnt1: how many unique minimizers + ///cnt2: how many non-unique minimizers + ///uint32_t cnt2, cnt1; + uint64_t m[2]; + int8_t s; +} pt_uinfo_t; + +typedef struct { + kvec_t(pt_node_t) p; + ma_ug_t *ug; + asg_t *rg; + trans_chain* t_ch; + ///kvec_t(int8_t) s; ///status + kvec_t(pt_uinfo_t) info; ///status + kvec_t(uint32_t) p_idx; + pt_match_t* e; +}pt_g_t; + +void pt_solve(hap_overlaps_list* ovlp, trans_chain* t_ch, ma_ug_t *ug, asg_t *read_g, double f_rate, uint8_t* trio_flag); + +#endif \ No newline at end of file