From d3b9de41200651bdd72463ebe3cadfaa371f0b9f Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Wed, 16 Mar 2022 01:13:27 -0400 Subject: [PATCH] new gchain --- Assembly.cpp | 26 + Makefile | 5 +- Overlaps.cpp | 138 +++--- Overlaps.h | 16 +- gfa_ut.cpp | 1296 ++++++++++++++++++++++++++++++++++++++++++++++++++ gfa_ut.h | 18 + inter.cpp | 387 +++++++++++++-- inter.h | 2 + 8 files changed, 1771 insertions(+), 117 deletions(-) create mode 100644 gfa_ut.cpp create mode 100644 gfa_ut.h diff --git a/Assembly.cpp b/Assembly.cpp index d20bf6f..f789b47 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1476,6 +1476,31 @@ void Output_PAF() fprintf(stderr, "PAF has been written.\n"); } +void Output_yak_binning() +{ + fprintf(stderr, "Writing binning to disk ...... \n"); + char* paf_name = (char*)malloc(strlen(asm_opt.output_file_name)+50); + sprintf(paf_name, "%s.hap1.bin.log", asm_opt.output_file_name); + FILE* oh1 = fopen(paf_name, "w"); + sprintf(paf_name, "%s.hap2.bin.log", asm_opt.output_file_name); + FILE* oh2 = fopen(paf_name, "w"); + uint64_t i; + + for (i = 0; i < R_INF.total_reads; i++) { + if(R_INF.trio_flag[i]==FATHER) { + fprintf(oh1, "%.*s\n", (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + } + if(R_INF.trio_flag[i]==MOTHER) { + fprintf(oh2, "%.*s\n", (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); + } + } + + free(paf_name); + fclose(oh1); fclose(oh2); + fprintf(stderr, "Binning has been written.\n"); +} + + int check_cluster(uint64_t* list, long long listLen, ma_hit_t_alloc* paf, float threshold) { long long i, k; @@ -1723,6 +1748,7 @@ int ha_assemble(void) if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads(); if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF(); if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2; + if (asm_opt.fn_bin_yak[0] && asm_opt.fn_bin_yak[1]) Output_yak_binning(); } if (!ovlp_loaded) { ha_flt_tab = ha_idx = NULL; diff --git a/Makefile b/Makefile index 8371724..30c843c 100644 --- a/Makefile +++ b/Makefile @@ -7,7 +7,7 @@ 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 rcut.o horder.o \ - tovlp.o inter.o kalloc.o + tovlp.o inter.o kalloc.o gfa_ut.o EXE= hifiasm LIBS= -lz -lpthread -lm @@ -77,4 +77,5 @@ rcut.o: rcut.h horder.o: horder.h tovlp.o: tovlp.h inter.o: inter.h Process_Read.h -kalloc.o: kalloc.h \ No newline at end of file +kalloc.o: kalloc.h +gfa_ut.o: Overlaps.h \ No newline at end of file diff --git a/Overlaps.cpp b/Overlaps.cpp index 4bb95fa..66a5617 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -16,6 +16,7 @@ #include "rcut.h" #include "horder.h" #include "inter.h" +#include "gfa_ut.h" uint32_t debug_purge_dup = 0; @@ -1880,10 +1881,11 @@ ma_sub_t* max_left, ma_sub_t* max_right, float overlap_rate, uint32_t trio_flag) } -void collect_sides(ma_hit_t_alloc* paf, uint64_t rLen, ma_sub_t* max_left, ma_sub_t* max_right) +void collect_sides(uint32_t rid, ma_hit_t_alloc* pafs, all_ul_t *x, uint64_t rLen, ma_sub_t* max_left, ma_sub_t* max_right) { long long j; uint32_t qs, qe; + ma_hit_t_alloc *paf = (&pafs[rid]); for (j = 0; j < paf->length; j++) { if(paf->buffer[j].del) continue; @@ -1910,10 +1912,36 @@ void collect_sides(ma_hit_t_alloc* paf, uint64_t rLen, ma_sub_t* max_left, ma_su ///this overlap would be added to both b_left and b_right ///that is what we want } + + if(x) { + uint64_t *a = NULL, a_n, k; + uc_block_t *p = NULL; + a = get_hifi2ul_list(x, rid, &a_n); + for (k = 0; k < a_n; k++) { + p = &(x->a[a[k]>>32].bb.a[(uint32_t)(a[k])]); + if(p->hid&x->mm) continue;///should not happen + qs = p->ts; qe = p->te;///note here is ts && te, instead of qs && qe + + ///overlaps from left side + if(qs == 0){ + if(qs < max_left->s) max_left->s = qs; + if(qe > max_left->e) max_left->e = qe; + } + + ///overlaps from right side + if(qe == rLen){ + if(qs < max_right->s) max_right->s = qs; + if(qe > max_right->e) max_right->e = qe; + } + ///note: if (qs == 0 && qe == rLen) + ///this overlap would be added to both b_left and b_right + ///that is what we want + } + } } void collect_contain(ma_hit_t_alloc* paf1, ma_hit_t_alloc* paf2, uint64_t rLen, -ma_sub_t* max_left, ma_sub_t* max_right, float overlap_rate) +ma_sub_t* max_left, ma_sub_t* max_right, float overlap_rate, all_ul_t *x, uint64_t xid) { long long j, new_left_e, new_right_s; new_left_e = max_left->e; @@ -1978,6 +2006,34 @@ ma_sub_t* max_left, ma_sub_t* max_right, float overlap_rate) } } } + + if(x) { + uint64_t *a = NULL, a_n, k; + uc_block_t *p = NULL; + a = get_hifi2ul_list(x, xid, &a_n); + for (k = 0; k < a_n; k++) { + p = &(x->a[a[k]>>32].bb.a[(uint32_t)(a[k])]); + if(p->hid&x->mm) continue;///should not happen + qs = p->ts; qe = p->te;///note here is ts && te, instead of qs && qe + ///check contained overlaps + if(qs != 0 && qe != rLen) + { + ///[qs, qe), [max_left.s, max_left.e) + if(qs < max_left->e && qe > max_left->e && max_left->e - qs > (overlap_rate * (qe -qs))) + { + ///if(qe > max_left->e) max_left->e = qe; + if(qe > max_left->e && qe > new_left_e) new_left_e = qe; + } + + ///[qs, qe), [max_right.s, max_right.e) + if(qs < max_right->s && qe > max_right->s && qe - max_right->s > (overlap_rate * (qe -qs))) + { + ///if(qs < max_right->s) max_right->s = qs; + if(qs < max_right->s && qs < new_right_s) new_right_s = qs; + } + } + } + } max_left->e = new_left_e; @@ -2091,7 +2147,7 @@ void print_overlaps(ma_hit_t_alloc* paf, long long rLen, long long interval_s, l void detect_chimeric_reads(ma_hit_t_alloc* paf, long long n_read, uint64_t* readLen, -ma_sub_t* coverage_cut, float shift_rate) +ma_sub_t* coverage_cut, float shift_rate, all_ul_t *x) { double startTime = Get_T(); init_aux_table(); @@ -2109,7 +2165,7 @@ ma_sub_t* coverage_cut, float shift_rate) max_left.s = max_right.s = rLen; max_left.e = max_right.e = 0; - collect_sides(&(paf[i]), rLen, &max_left, &max_right); + collect_sides(i, paf, x, rLen, &max_left, &max_right); ///collect_sides(&(rev_paf[i]), rLen, &max_left, &max_right); ///that means this read is an end node if(max_left.s == rLen || max_right.s == rLen) @@ -2117,7 +2173,7 @@ ma_sub_t* coverage_cut, float shift_rate) continue; } - collect_contain(&(paf[i]), NULL, rLen, &max_left, &max_right, 0.1); + collect_contain(&(paf[i]), NULL, rLen, &max_left, &max_right, 0.1, x, i); ///collect_contain(&(paf[i]), &(rev_paf[i]), rLen, &max_left, &max_right, 0.1); ////shift_rate should be (asm_opt.max_ov_diff_final*2) @@ -4110,21 +4166,21 @@ int check_if_cross(asg_t *g, uint32_t v) f1 = detect_bubble_end_with_bubbles(g, N_list[0]^1, N_list[3]^1, &convex1, &l1, NULL); f2 = detect_bubble_end_with_bubbles(g, N_list[1]^1, N_list[2]^1, &convex2, &l2, NULL); - if(f1 && f2) + if(f1 && f2)///full bubble { if(l1 > min_thres && l2 > min_thres) { todel = 1; } } - else if(f1) + else if(f1)//semi bubble { if(l1 > min_thres) { todel = 1; } } - else if(f2) + else if(f2)//semi bubble { if(l2 > min_thres) { @@ -5071,10 +5127,6 @@ int asg_arc_del_trans(asg_t *g, int fuzz) return n_reduced; } - - - - ///max_ext is 4 int asg_cut_tip(asg_t *g, int max_ext) { @@ -5980,7 +6032,8 @@ ma_hit_t_alloc* reverse_sources, long long min_edge_length, R_to_U* ruIndex) else if(av[i].ol < drop_ratio_Len && check_if_diploid(v_max, av[i].v, g, reverse_sources, min_edge_length, ruIndex) != 1) { - av[i].ol = 1; + // av[i].ol = 1;///should be a bug + av[i].del = 1; asg_arc_del(g, av[i].v^1, av[i].ul>>32^1, 1); ++n_short; } @@ -6052,7 +6105,7 @@ static uint32_t asg_check_unambi1(asg_t *g, uint32_t v) return av[k].v; } ///to see if it is a long tip -static int asg_topocut_aux(asg_t *g, uint32_t v, int max_ext) +int asg_topocut_aux(asg_t *g, uint32_t v, int max_ext) { int32_t n_ext; for (n_ext = 1; n_ext < max_ext && v != (uint32_t)-1; ++n_ext) { @@ -6632,7 +6685,7 @@ ma_hit_t_alloc* reverse_sources, long long miniedgeLen, R_to_U* ruIndex) kv_push(uint32_t, b_r, w); aw = asg_arc_a(g, w); - min_edge = (u_int32_t)-1; + min_edge = (uint32_t)-1; for (t = 0; t < nw; t++) { if(aw[t].del) continue; @@ -19768,7 +19821,8 @@ long long miniedgeLen, R_to_U* ruIndex) { ///fprintf(stderr, "v: %u, v_max: %u, av[%d].v: %u\n", v>>1, v_max>>1, i, av[i].v>>1); - av[i].ol = 1; + // av[i].ol = 1;///should be a bug + av[i].del = 1; asg_arc_del(g, av[i].v^1, av[i].ul>>32^1, 1); n_reduced++; } @@ -25666,7 +25720,7 @@ void pre_clean(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *sg, uint3 ///remove isoloated single read if(pop_s_node) { - tri_flag += asg_arc_del_single_node_directly(sg, asm_opt.max_short_tip, sources); + tri_flag += asg_arc_del_single_node_directly(sg, asm_opt.max_short_tip, sources);///remove very small bubbles } // if ((!ha_opt_triobin(&asm_opt))&&(!ha_opt_hic(&asm_opt))) @@ -31101,27 +31155,6 @@ char *get_outfile_name(char* output_file_name) return buf; } -asg_t *build_init_sg(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int64_t n_read, -int64_t min_dp, uint64_t* readLen, int64_t mini_overlap_length, int64_t max_hang_length, -ma_sub_t *coverage_cut, R_to_U* ruIndex) -{ - asg_t *sg = NULL; - clean_weak_ma_hit_t(sources, reverse_sources, n_read); - ///ma_hit_sub is just use to init coverage_cut, - ///it seems we do not need ma_hit_cut & ma_hit_flt - ma_hit_sub(min_dp, sources, n_read, readLen, mini_overlap_length, &coverage_cut); - detect_chimeric_reads(sources, n_read, readLen, coverage_cut, asm_opt.max_ov_diff_final * 2.0); - - ma_hit_cut(sources, n_read, readLen, mini_overlap_length, &coverage_cut); - ///print_binned_reads(sources, n_read, coverage_cut); - ma_hit_flt(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length); - ///fix_binned_reads(sources, n_read, coverage_cut); - ///just need to deal with trio here - ma_hit_contained_advance(sources, n_read, coverage_cut, ruIndex, max_hang_length, mini_overlap_length); - sg = ma_sg_gen(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length); - return sg; -} - void create_ul_info(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int64_t max_hang, int64_t min_ovlp, int64_t gap_fuzz, int64_t min_dp, uint64_t* readLen, ma_sub_t *coverage_cut, R_to_U* ruIndex) { @@ -31138,6 +31171,7 @@ int64_t min_dp, uint64_t* readLen, ma_sub_t *coverage_cut, R_to_U* ruIndex) ul_load(&opt); } + void clean_graph( int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long n_read, uint64_t* readLen, long long mini_overlap_length, @@ -31155,10 +31189,8 @@ ma_sub_t **coverage_cut_ptr, int debug_g) init_bub_label_t(&b_mask_t, MIN(10, asm_opt.thread_num), sg->n_seq); goto debug_gfa; } - ///just for debug renew_graph_init(sources, reverse_sources, sg, coverage_cut, ruIndex, n_read); - ///it's hard to say which function is better ///normalize_ma_hit_t_single_side(sources, n_read); normalize_ma_hit_t_single_side_advance(sources, n_read); @@ -31182,12 +31214,10 @@ ma_sub_t **coverage_cut_ptr, int debug_g) } ///print_binned_reads(sources, n_read, coverage_cut); - ///ma_hit_sub is just use to init coverage_cut, ///it seems we do not need ma_hit_cut & ma_hit_flt ma_hit_sub(min_dp, sources, n_read, readLen, mini_overlap_length, &coverage_cut); - detect_chimeric_reads(sources, n_read, readLen, coverage_cut, asm_opt.max_ov_diff_final * 2.0); - + detect_chimeric_reads(sources, n_read, readLen, coverage_cut, asm_opt.max_ov_diff_final * 2.0, asm_opt.ar?&UL_INF:NULL); ma_hit_cut(sources, n_read, readLen, mini_overlap_length, &coverage_cut); ///print_binned_reads(sources, n_read, coverage_cut); ma_hit_flt(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length); @@ -31195,12 +31225,9 @@ ma_sub_t **coverage_cut_ptr, int debug_g) ///just need to deal with trio here ma_hit_contained_advance(sources, n_read, coverage_cut, ruIndex, max_hang_length, mini_overlap_length); sg = ma_sg_gen(sources, n_read, coverage_cut, max_hang_length, mini_overlap_length); - ///debug_info_of_specfic_node((char*)"m64043_200504_050026/93784180/ccs", sg, ruIndex, (char*)"sbsbsb"); init_bub_label_t(&b_mask_t, MIN(10, asm_opt.thread_num), sg->n_seq); - asg_arc_del_trans(sg, gap_fuzz); - asm_opt.coverage = get_coverage(sources, coverage_cut, n_read); if(VERBOSE >= 1) @@ -31210,7 +31237,9 @@ ma_sub_t **coverage_cut_ptr, int debug_g) output_read_graph(sg, coverage_cut, unlean_name, n_read); free(unlean_name); } - + ul_clean_gfa(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)); + /** asg_cut_tip(sg, asm_opt.max_short_tip); ///debug_info_of_specfic_node("m64043_200505_112554/8849050/ccs", sg, "inner_1"); ///drop_inexact_edegs_at_bubbles(sg, bubble_dist); @@ -31246,15 +31275,12 @@ ma_sub_t **coverage_cut_ptr, int debug_g) ///asg_arc_del_orthology(sg, reverse_sources, drop_ratio, asm_opt.max_short_tip); // asg_arc_del_orthology_multiple_way(sg, reverse_sources, drop_ratio, asm_opt.max_short_tip); // asg_cut_tip(sg, asm_opt.max_short_tip); - /****************************may have bugs********************************/ asg_arc_identify_simple_bubbles_multi(sg, &b_mask_t, 1); //reomve edge between two chromesomes //this node must be a single read asg_arc_del_false_node(sg, sources, asm_opt.max_short_tip); asg_cut_tip(sg, asm_opt.max_short_tip); - /****************************may have bugs********************************/ - /****************************may have bugs********************************/ ///asg_arc_identify_simple_bubbles_multi(sg, 1); asg_arc_identify_simple_bubbles_multi(sg, &b_mask_t, 0); @@ -31268,7 +31294,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g) asg_arc_del_short_diploid_by_exact(sg, asm_opt.max_short_tip, sources); } asg_cut_tip(sg, asm_opt.max_short_tip); - /****************************may have bugs********************************/ asg_arc_identify_simple_bubbles_multi(sg, &b_mask_t, 1); if (ha_opt_triobin(&asm_opt)) @@ -31285,7 +31310,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g) asg_arc_identify_simple_bubbles_multi(sg, &b_mask_t, 1); asg_arc_del_short_false_link(sg, 0.6, 0.85, bubble_dist, reverse_sources, asm_opt.max_short_tip, ruIndex); - + asg_arc_identify_simple_bubbles_multi(sg, &b_mask_t, 1); asg_arc_del_complex_false_link(sg, 0.6, 0.85, bubble_dist, reverse_sources, asm_opt.max_short_tip); @@ -31314,11 +31339,11 @@ ma_sub_t **coverage_cut_ptr, int debug_g) asg_arc_identify_simple_bubbles_multi(sg, &b_mask_t, 0); - asg_arc_del_too_short_overlaps(sg, 2000, min_ovlp_drop_ratio, reverse_sources, - asm_opt.max_short_tip, ruIndex); + asg_arc_del_too_short_overlaps(sg, 2000, min_ovlp_drop_ratio, reverse_sources, asm_opt.max_short_tip, ruIndex); asg_cut_tip(sg, asm_opt.max_short_tip); asg_arc_del_simple_circle_untig(sources, coverage_cut, sg, 100, 0); + **/ ///note: don't apply asg_arc_del_too_short_overlaps() after this function!!!! rescue_contained_reads_aggressive(NULL, sg, sources, coverage_cut, ruIndex, max_hang_length, mini_overlap_length, 10, 1, 0, NULL, NULL, &b_mask_t); @@ -31415,10 +31440,8 @@ long long bubble_dist, int read_graph, int write) asg_t *sg = NULL; ma_sub_t* coverage_cut = NULL; init_aux_table(); - ///actually min_thres = asm_opt.max_short_tip + 1 there are asm_opt.max_short_tip reads min_thres = asm_opt.max_short_tip + 1; - if (asm_opt.flag & HA_F_VERBOSE_GFA) { if(load_debug_graph(&sg, &sources, &coverage_cut, output_file_name, &reverse_sources, &ruIndex)) @@ -31434,19 +31457,16 @@ long long bubble_dist, int read_graph, int write) return; } } - if (asm_opt.write_index_to_disk && write) { write_all_data_to_disk(sources, reverse_sources, &R_INF, output_file_name); } - ///debug_info_of_specfic_read("m64011_190830_220126/31720629/ccs", sources, reverse_sources, -1, "beg"); if (!(asm_opt.flag & HA_F_BAN_ASSEMBLY)) { try_rescue_overlaps(sources, reverse_sources, n_read, 4); - clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length, max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio, output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 0); diff --git a/Overlaps.h b/Overlaps.h index 92a465e..2067309 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -141,14 +141,17 @@ typedef struct { #define u_trans_a(x, id) ((x).a + ((x).idx.a[(id)]>>32)) #define u_trans_n(x, id) ((uint32_t)((x).idx.a[(id)])) +#define OU_MASK (0x3fffU) typedef struct { uint64_t ul; uint32_t v; uint32_t ol:31, del:1; - uint8_t strong; + uint16_t ou:14, strong:1, no_l_indel:1; uint8_t el; - uint8_t no_l_indel; + // uint8_t strong; + // uint8_t el; + // uint8_t no_l_indel; } asg_arc_t; typedef struct { @@ -1048,9 +1051,12 @@ uint32_t tn, kv_u_trans_hit_t* ktb, uint32_t bn); void clean_u_trans_t_idx(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g); uint32_t test_dbug(ma_ug_t* ug, FILE* fp); void write_dbug(ma_ug_t* ug, FILE* fp); -asg_t *build_init_sg(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, int64_t n_read, -int64_t min_dp, uint64_t* readLen, int64_t mini_overlap_length, int64_t max_hang_length, -ma_sub_t *coverage_cut, R_to_U* ruIndex); +int asg_arc_identify_simple_bubbles_multi(asg_t *g, bub_label_t* x, int check_cross); +uint8_t get_tip_trio_infor(asg_t *sg, uint32_t begNode); +int asg_topocut_aux(asg_t *g, uint32_t v, int max_ext); +int asg_arc_del_triangular_directly(asg_t *g, long long min_edge_length, +ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex); +int asg_arc_del_short_diploid_by_exact(asg_t *g, int max_ext, ma_hit_t_alloc* sources); #define JUNK_COV 5 #define DISCARD_RATE 0.8 diff --git a/gfa_ut.cpp b/gfa_ut.cpp new file mode 100644 index 0000000..ebdb935 --- /dev/null +++ b/gfa_ut.cpp @@ -0,0 +1,1296 @@ +#include +#include +#include +#include +#include +#include "kthread.h" +#include "gfa_ut.h" +#include "CommandLines.h" +#include "Correct.h" + +#define generic_key(x) (x) +KRADIX_SORT_INIT(srt64, uint64_t, generic_key, 8) +#define OU_NOISY 2 +#define ASG_ET_MERGEABLE 0 +#define ASG_ET_TIP 1 +#define ASG_ET_MULTI_OUT 2 +#define ASG_ET_MULTI_NEI 3 + +typedef struct { + asg_t *g; + ma_hit_t_alloc *src; +} sset_aux; + +int32_t if_sup_chimeric(ma_hit_t_alloc* src, uint64_t rLen, asg64_v *b, int if_exact); + +void print_edge(asg_arc_t *t, const char *cmd) +{ + uint32_t v = t->ul>>32, w = t->v; + fprintf(stderr, "%s: v->%u(%c)[%u], w->%u(%c)[%u], el->%u, del->%u\n", cmd, v>>1, "+-"[v&1], v, w>>1, "+-"[w&1], w, t->el, t->del); +} + +void stats_chimeric(asg_t *g, ma_hit_t_alloc* src, asg64_v *in) +{ + asg64_v tx = {0,0,0}, *b = NULL; + uint32_t v, s[2] = {0}; + if(in) b = in; + else b = &tx; + b->n = 0; + + for (v = 0; v < g->n_seq; ++v) { + if (g->seq[v].del) continue; + s[if_sup_chimeric(&(src[v]), g->seq[v].len, b, 1)]++; + } + + fprintf(stderr, "[M::%s::] ==> # non-chimeric:%u, # chimeric:%u\n", __func__, s[0], s[1]); + if(!in) free(tx.a); +} + +static void stats_sysm_worker(void *_data, long eid, int tid) +{ + asg_t *g = (asg_t*)_data; + asg_arc_t *p = &(g->arc[eid]); + if(p->del) return; + uint32_t k, v = p->v^1, w = (p->ul>>32)^1, nv; asg_arc_t *av; + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + for (k = 0; k < nv; k++) { + if (av[k].del || av[k].v!=w) continue; + break; + } + assert(k < nv); + + v = p->ul>>32; w = p->v; + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + for (k = 0; k < nv; k++) { + if (av[k].del || av[k].v!=w) continue; + assert((uint32_t)eid == av-g->arc+k); + } +} + +void stats_sysm(asg_t *g) { + kt_for(asm_opt.thread_num, stats_sysm_worker, g, g->n_arc); + fprintf(stderr, "[M::%s::]", __func__); +} + +uint32_t get_arcs(asg_t *g, uint32_t v, uint32_t* idx, uint32_t idx_n) +{ + uint32_t i, kv = 0, an = asg_arc_n(g, v), beg = g->idx[v]>>32; + for (i = 0, kv = 0; i < an; i++) { + if(g->arc[beg+i].del) continue; + if(idx && kvb, v>>1); + if(b) kv_push(uint64_t, *b, v); + + if(kv == 0) return END_TIPS; + if(kv == 2) return TWO_OUTPUT; + if(kv > 2) return MUL_OUTPUT; + if((*occ) > lim) return LONG_TIPS; + w = g->arc[w].v; + ///up to here, kv=1 + ///kw must >= 1 + kw = get_arcs(g, w^1, NULL, 0); + v = w; + + if(kw == 2) return TWO_INPUT; + if(kw > 2) return MUL_INPUT; + if(v == s) return LOOP; + } + + return LONG_TIPS; +} + +static inline int asg_end(const asg_t *g, uint32_t v, uint64_t *lw, uint32_t *ou) +{ + ///v^1 is the another direction of v + uint32_t w, nv, nw, nw0, nv0 = asg_arc_n(g, v^1); + int i, i0 = -1; + asg_arc_t *aw, *av = asg_arc_a(g, v^1); + + ///if this arc has not been deleted + for (i = nv = 0; i < (int)nv0; ++i) + if (!av[i].del) i0 = i, ++nv; + + ///end without any out-degree + if (nv == 0) return ASG_ET_TIP; // tip + if (nv > 1) return ASG_ET_MULTI_OUT; // multiple outgoing arcs + ///until here, nv == 1 + if (lw) *lw = av[i0].ul<<32 | av[i0].v; + if (ou) *ou = av[i0].ou; + w = av[i0].v ^ 1; + nw0 = asg_arc_n(g, w); + aw = asg_arc_a(g, w); + for (i = nw = 0; i < (int)nw0; ++i) + if (!aw[i].del) ++nw; + + if (nw != 1) return ASG_ET_MULTI_NEI; + return ASG_ET_MERGEABLE; +} + +uint32_t asg_arc_cut_tips(asg_t *g, uint32_t max_ext, asg64_v *in, uint32_t is_ou) +{ + asg64_v tx = {0,0,0}, *b = NULL; + uint32_t n_vtx = g->n_seq<<1, v, w, i, k, cnt = 0, nv, kv, pb, ou, mm_ou; + asg_arc_t *av = NULL; uint64_t lw; + if(in) b = in; + else b = &tx; + b->n = 0; + for (v = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + + av = asg_arc_a(g, v^1); nv = asg_arc_n(g, v^1); + for (i = kv = 0; i < nv; i++) { + if (av[i].del) continue; + kv++; break; + } + + if(kv) continue; + kv = 1; mm_ou = (uint32_t)-1; ou = 0; + for (i = 0, w = v; i < max_ext; i++) { + if(asg_end(g, w^1, &lw, is_ou?&ou:NULL)!=0) break; + w = (uint32_t)lw; kv++; mm_ou = MIN(mm_ou, ou); + } + if(mm_ou == (uint32_t)-1) mm_ou = 0; + kv += mm_ou; i += mm_ou; + if(i < max_ext + (!!is_ou)) kv_push(uint64_t, *b, (((uint64_t)kv)<<32)|v); + } + + radix_sort_srt64(b->a, b->a + b->n); + + for (k = 0; k < b->n; k++) { + v = (uint32_t)(b->a[k]); + + if (g->seq[v>>1].del) continue; + + av = asg_arc_a(g, v^1); nv = asg_arc_n(g, v^1); + for (i = kv = 0; i < nv; i++) { + if (av[i].del) continue; + kv++; break; + } + + if(kv) continue; + pb = b->n; kv_push(uint64_t, *b, v); mm_ou = (uint32_t)-1; ou = 0; + for (i = 0, w = v; i < max_ext; i++) { + if(asg_end(g, w^1, &lw, is_ou?&ou:NULL)!=0) break; + w = (uint32_t)lw; kv_push(uint64_t, *b, lw); mm_ou = MIN(mm_ou, ou); + } + if(mm_ou == (uint32_t)-1) mm_ou = 0; + i += mm_ou; + + if(i < max_ext + (!!is_ou)) { + for (i = pb; i < b->n; i++) asg_seq_del(g, ((uint32_t)b->a[i])>>1); + cnt++; + } + b->n = pb; + } + + + /** + for (v = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + + av = asg_arc_a(g, v^1); nv = asg_arc_n(g, v^1); + for (i = kv = 0; i < nv; i++) { + if (av[i].del) continue; + kv++; + } + + if(kv) continue; + pb = b->n; kv_push(uint64_t, *b, v); + for (i = 0, w = v; i < max_ext; i++) { + if(asg_is_utg_end(g, w^1, &lw)!=0) break; + w = (uint32_t)lw; kv_push(uint64_t, *b, lw); + } + + if(i < max_ext) { + for (i = pb; i < b->n; i++) asg_seq_del(g, ((uint32_t)b->a[i])>>1); + cnt++; + } + b->n = pb; + } + **/ + // stats_sysm(g); + if(!in) free(tx.a); + if (cnt > 0) asg_cleanup(g); + + return cnt; +} + +static void update_sg_uo_t(void *data, long i, int tid) +{ + sset_aux *sl = (sset_aux *)data; + ma_hit_t_alloc *src = sl->src; asg_t *g = sl->g; + asg_arc_t *e = &(g->arc[i]); uint32_t k, qn, tn; + + for (k = 0; k < src[i].length; k++) { + qn = Get_qn(src[i].buffer[k]); + tn = Get_tn(src[i].buffer[k]); + if(qn == (e->ul>>33) && tn == (e->v>>1)) { + e->ou = (src[i].buffer[k].bl&OU_MASK); + break; + } + } + assert(k < src[i].length); +} + +void update_sg_uo(asg_t *g, ma_hit_t_alloc *src) +{ + sset_aux s; s.g = g; s.src = src; + kt_for(asm_opt.thread_num, update_sg_uo_t, &s, g->n_arc); +} + +int32_t if_sup_chimeric(ma_hit_t_alloc* src, uint64_t rLen, asg64_v *b, int if_exact) +{ + uint32_t k, qs, qe, l[2], r[2], st, bn; + int32_t dp, op; + l[0] = r[0] = rLen; l[1] = r[1] = 0; + for (k = 0; k < src->length; k++){ + if(src->buffer[k].del) continue; + if(if_exact && !(src->buffer[k].el)) continue; + + qs = Get_qs(src->buffer[k]); qe = Get_qe(src->buffer[k]); + + ///overlaps from left side + if(qs == 0){ + if(qs < l[0]) l[0] = qs; + if(qe > l[1]) l[1] = qe; + } + + ///overlaps from right side + if(qe == rLen){ + if(qs < r[0]) r[0] = qs; + if(qe > r[1]) r[1] = qe; + } + + ///note: if (qs == 0 && qe == rLen) + ///this overlap would be added to both b_left and b_right + ///that is what we want + } + if (l[1] > r[0]) return 0; + if (l[1] <= l[0] || r[1] <= r[0]) return 1; + + bn = b->n; + if(l[1] > l[0]) { + kv_push(uint64_t, *b, (l[0]<<1)); kv_push(uint64_t, *b, (l[1]<<1)|1); + } + if(r[1] > r[0]) { + kv_push(uint64_t, *b, (r[0]<<1)); kv_push(uint64_t, *b, (r[1]<<1)|1); + } + + ///check contained overlaps + for (k = 0; k < src->length; k++) { + if(src->buffer[k].del) continue; + if(if_exact && !(src->buffer[k].el)) continue; + + qs = Get_qs(src->buffer[k]); qe = Get_qe(src->buffer[k]); + if(qs == 0 || qe == rLen) continue; + + kv_push(uint64_t, *b, (qs<<1)); kv_push(uint64_t, *b, (qe<<1)|1); + } + radix_sort_srt64(b->a + bn, b->a + b->n); + l[0] = r[0] = rLen; l[1] = r[1] = 0; + + for (k = bn, dp = st = 0; k < b->n; k++) { + op = dp; + ///if a[j] is qe + if (b->a[k]&1) --dp; + else ++dp; + + if(op < 1 && dp >= 1) { + st = b->a[k]>>1; + } else if(op >= 1 && dp < 1) { + if(st == 0) l[0] = st, l[1] = b->a[k]>>1; + if((b->a[k]>>1) == rLen) r[0] = st, r[1] = b->a[k]>>1; + } + } + + b->n = bn; + if (l[1] > r[0]) return 0; + + return 1; +} + +///remove single node +void asg_arc_cut_chimeric(asg_t *g, ma_hit_t_alloc* src, asg64_v *in) +{ + asg64_v tx = {0,0,0}, *b = NULL; + uint32_t v, w, ei[2] = {0}, k, i, n_vtx = g->n_seq<<1; + uint32_t nw, el_n, cnt = 0; asg_arc_t *aw; + if(in) b = in; + else b = &tx; + b->n = 0; + + for (v = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + if(g->seq_vis[v] == 0) { + if((get_arcs(g, v, &(ei[0]), 1)!=1) || (get_arcs(g, v^1, &(ei[1]), 1)!=1)) continue; + assert((g->arc[ei[0]].ul>>32) == v && (g->arc[ei[1]].ul>>32) == (v^1)); + if((get_arcs(g, g->arc[ei[0]].v^1, NULL, 0)<2) || (get_arcs(g, g->arc[ei[1]].v^1, NULL, 0)<2)) continue; + if(g->arc[ei[0]].el) continue; + if(!if_sup_chimeric(&(src[v>>1]), g->seq[v>>1].len, b, 1)) continue; + kv_push(uint64_t, *b, (((uint64_t)(g->arc[ei[0]].ol))<<32)|((uint64_t)(ei[0]))); + } + } + + radix_sort_srt64(b->a, b->a + b->n); + ///here all edges are inexact matches + for (k = 0; k < b->n; k++) { + if(g->arc[(uint32_t)b->a[k]].del) continue; + v = g->arc[(uint32_t)b->a[k]].ul>>32; w = g->arc[(uint32_t)b->a[k]].v^1; + if(g->seq[v>>1].del || g->seq[w>>1].del) continue; + aw = asg_arc_a(g, w); nw = asg_arc_n(g, w); + if((get_arcs(g, v, &(ei[0]), 1)!=1) || (get_arcs(g, v^1, &(ei[1]), 1)!=1)) continue; + if((get_arcs(g, g->arc[ei[0]].v^1, NULL, 0)<2) || (get_arcs(g, g->arc[ei[1]].v^1, NULL, 0)<2)) continue; + + for (i = el_n = 0; i < nw; i++) { + if ((aw[i].del) || (aw[i].v==(v^1)) || (!aw[i].el)) continue; + el_n++; break; + } + + if(!el_n) continue; + asg_seq_del(g, v>>1); + cnt++; + } + // stats_sysm(g); + if(!in) free(tx.a); + if (cnt > 0) asg_cleanup(g); +} + +void asg_arc_cut_inexact(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, int32_t max_ext, uint32_t is_ou, uint32_t is_trio/**, asg64_v *dbg**/) +{ + asg64_v tx = {0,0,0}, *b = NULL; + uint32_t v, w, i, k, n_vtx = g->n_seq<<1; + asg_arc_t *av, *aw, *ve, *vmax, *we; uint32_t nv, nw, kv, kw, ol_max, ou_max, to_del, cnt = 0, mm_ol, mm_ou; + uint32_t trioF = (uint32_t)-1, ntrioF = (uint32_t)-1; + if(in) b = in; + else b = &tx; + b->n = 0; + + for (v = 0; v < n_vtx; ++v) { + if(g->seq[v>>1].del) continue; + if(g->seq_vis[v] == 0) { + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + if (nv < 2) continue; + + for (i = kv = 0; i < nv; ++i) { + if(av[i].del) continue; + kv++; + } + if(kv < 2) continue; + + for (i = 0; i < nv; ++i) { + if(av[i].del || av[i].el) continue; + kv_push(uint64_t, *b, (uint64_t)((((uint64_t)av[i].ol)<<32)|((uint64_t)(av-g->arc+i)))); + } + } + } + + radix_sort_srt64(b->a, b->a + b->n); + for (k = 0; k < b->n; k++) { + if(g->arc[(uint32_t)b->a[k]].del) continue; + assert((!g->arc[(uint32_t)b->a[k]].el)); + + v = g->arc[(uint32_t)b->a[k]].ul>>32; w = g->arc[(uint32_t)b->a[k]].v^1; + if(g->seq[v>>1].del || g->seq[w>>1].del) continue; + nv = asg_arc_n(g, v); nw = asg_arc_n(g, w); + av = asg_arc_a(g, v); aw = asg_arc_a(g, w); + // if(((v>>1) == 50356 && (w>>1) == 1276292)||((v>>1) == 1276292 && (w>>1) == 50356)) { + // fprintf(stderr, "[0]v->%u, w->%u, nv->%u, nw->%u\n", v, w, nv, nw); + // } + if(nv<=1 && nw <= 1) continue; + if(is_trio) { + if(get_arcs(g, v, NULL, 0)<=1 && get_arcs(g, w, NULL, 0)<=1) continue;///speedup + trioF = get_tip_trio_infor(g, v^1); + ntrioF = (trioF==FATHER? MOTHER : (trioF==MOTHER? FATHER : (uint32_t)-1)); + } + ve = &(g->arc[(uint32_t)b->a[k]]); + for (i = 0; i < nw; ++i) { + if (aw[i].v == (v^1)) { + we = &(aw[i]); + break; + } + } + ///mm_ol and mm_ou are used to make edge with long indel more easy to be cutted + mm_ol = MIN(ve->ol, we->ol); mm_ou = MIN(ve->ou, we->ou); + for (i = kv = ol_max = ou_max = 0, /**ve =**/ vmax = NULL; i < nv; ++i) { + if(av[i].del) continue; + // if(av[i].v == (w^1)) ve = &(av[i]); + kv++; + if(is_trio && get_tip_trio_infor(g, av[i].v) == ntrioF) continue; + if(ol_max < av[i].ol) ol_max = av[i].ol, vmax = &(av[i]); + if(ou_max < av[i].ou) ou_max = av[i].ou; + } + // if(((v>>1) == 50356 && (w>>1) == 1276292)||((v>>1) == 1276292 && (w>>1) == 50356)) { + // fprintf(stderr, "[1]v->%u, w->%u, kv->%u, ve->ol->%u, ol_max->%u\n", v, w, kv, ve->ol, ol_max); + // } + if (kv < 1) continue; + if (kv >= 2) { + if (/**ve->ol**/mm_ol >= ol_max) continue; + if (is_ou && /**ve->ou**/mm_ou >= ou_max) continue; + } + + for (i = kw = ol_max = ou_max = 0/**, we = NULL**/; i < nw; ++i) { + if(aw[i].del) continue; + // if(aw[i].v == (v^1)) we = &(aw[i]); + kw++; + if(is_trio && get_tip_trio_infor(g, aw[i].v) == ntrioF) continue; + if(ol_max < aw[i].ol) ol_max = aw[i].ol; + if(ou_max < aw[i].ou) ou_max = aw[i].ou; + } + // if(((v>>1) == 50356 && (w>>1) == 1276292)||((v>>1) == 1276292 && (w>>1) == 50356)) { + // fprintf(stderr, "[1]v->%u, w->%u, kw->%u, we->ol->%u, ol_max->%u\n", v, w, kw, we->ol, ol_max); + // } + if (kw < 1) continue; + if (kw >= 2) { + if (/**we->ol**/mm_ol >= ol_max) continue; + if (is_ou && /**we->ou**/mm_ou >= ou_max) continue; + } + if (kv <= 1 && kw <= 1) continue; + + to_del = 0; + ///if there is an inexact edge between two good reads + if(src[v>>1].is_fully_corrected == 1 && src[w>>1].is_fully_corrected == 1) { + if (kv > 1 && kw > 1) { + to_del = 1; + } else if (kw == 1) { + if (asg_topocut_aux(g, w^1, max_ext) < max_ext) to_del = 1; + } else if (kv == 1) { + if (asg_topocut_aux(g, v^1, max_ext) < max_ext) to_del = 1; + } + } + ///TODO: should check if the edge wmax also works + if(src[v>>1].is_fully_corrected == 1 && src[w>>1].is_fully_corrected == 0) { + if(vmax && vmax->v != ve->v && vmax->el == 1 && src[vmax->v>>1].is_fully_corrected == 1) { + if (kv > 1 && kw > 1) { + to_del = 1; + } else if (kw == 1) { + if (asg_topocut_aux(g, w^1, max_ext) < max_ext) to_del = 1; + } else if (kv == 1) { + if (asg_topocut_aux(g, v^1, max_ext) < max_ext) to_del = 1; + } + } + } + + // if(((v>>1) == 50356 && (w>>1) == 1276292)||((v>>1) == 1276292 && (w>>1) == 50356)) { + // fprintf(stderr, "[3]v->%u, w->%u, to_del->%u, el->%u, src[v>>1].is_fully_corrected->%u, src[w>>1].is_fully_corrected->%u\n", + // v, w, to_del, g->arc[(uint32_t)b->a[k]].el, src[v>>1].is_fully_corrected, src[w>>1].is_fully_corrected); + // fprintf(stderr, "[4]v->%u, w->%u, to_del->%u, vmax->v->%u, vmax->el->%u, src[vmax->v>>1].is_fully_corrected->%u\n", + // v, w, to_del, vmax->v, vmax->el, src[vmax->v>>1].is_fully_corrected); + // } + + if (to_del) { + ve->del = we->del = 1, ++cnt; + /** + if(dbg) { + kv_push(uint64_t, *dbg, ve - g->arc); + kv_push(uint64_t, *dbg, we - g->arc); + // if(((v>>1) == 50356 && (w>>1) == 1276292)||((v>>1) == 1276292 && (w>>1) == 50356)) { + // fprintf(stderr, "[5]v->%u, w->%u, (ve - g->arc)->%u, (we - g->arc)->%u\n", + // v, w, (uint32_t)(ve - g->arc), (uint32_t)(we - g->arc)); + // } + } + **/ + } + + } + // stats_sysm(g); + if(!in) free(tx.a); + if (cnt > 0) asg_cleanup(g); +} + +void asg_arc_cut_inexact_debug(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, int32_t max_ext, uint32_t is_ou, uint32_t is_trio, asg64_v *dbg) +{ + asg64_v tx = {0,0,0}, *b = NULL; + uint32_t v, w, i, k, nv, nw, kv, kw, iv, iw, n_vtx = g->n_seq<<1, to_del, cnt = 0, ov_max = 0, ow_max = 0, ov_max_i = 0; + asg_arc_t *av = NULL, *aw = NULL, *a = NULL; + if(in) b = in; + else b = &tx; + b->n = 0; + + for (v = 0; v < n_vtx; ++v) { + if(g->seq[v>>1].del) continue; + if(g->seq_vis[v] == 0) { + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + if (nv < 2) continue; + for (i = kv = 0; i < nv; ++i) { + if(av[i].del) continue; + kv++; + } + if(kv < 2) continue; + + for (i = 0; i < nv; ++i) { + if(av[i].del || av[i].el) continue; + kv_push(uint64_t, *b, (uint64_t)((uint64_t)av[i].ol << 32 | (av - g->arc + i))); + } + } + } + radix_sort_srt64(b->a, b->a + b->n); + + for (k = 0; k < b->n; k++) { + a = &g->arc[(uint32_t)b->a[k]]; + if(a->del) continue; + v = (a->ul)>>32, w = a->v^1; to_del = 0; + nv = asg_arc_n(g, v), nw = asg_arc_n(g, w); + if(((v>>1) == 50356 && (w>>1) == 1276292)||((v>>1) == 1276292 && (w>>1) == 50356)) { + fprintf(stderr, "[0]v->%u, w->%u, nv->%u, nw->%u\n", v, w, nv, nw); + } + if (nv == 1 && nw == 1) continue; + av = asg_arc_a(g, v), aw = asg_arc_a(g, w); + ov_max = ow_max = ov_max_i = 0; + + for (i = 0, kv = 0; i < nv; ++i) { + if (av[i].del) continue; + if (ov_max < av[i].ol) { + ov_max = av[i].ol; + ov_max_i = i; + } + ++kv; + } + if(((v>>1) == 50356 && (w>>1) == 1276292)||((v>>1) == 1276292 && (w>>1) == 50356)) { + fprintf(stderr, "[1]v->%u, w->%u, kv->%u, a->ol->%u, ov_max->%u\n", v, w, kv, a->ol, ov_max); + } + if (kv >= 2 && a->ol == ov_max) continue; + + for (i = 0, kw = 0; i < nw; ++i) { + if (aw[i].del) continue; + if (ow_max < aw[i].ol) { + ow_max = aw[i].ol; + } + ++kw; + } + if(((v>>1) == 50356 && (w>>1) == 1276292)||((v>>1) == 1276292 && (w>>1) == 50356)) { + fprintf(stderr, "[2]v->%u, w->%u, kw->%u, a->ol->%u, ow_max->%u\n", v, w, kw, a->ol, ow_max); + } + if (kw >= 2 && a->ol == ow_max) continue; + + if (kv <= 1 && kw <= 1) continue; + + ///to see which one is the current edge (from v and w) + for (iv = 0; iv < nv; ++iv) + if (av[iv].v == (w^1)) break; + for (iw = 0; iw < nw; ++iw) + if (aw[iw].v == (v^1)) break; + ///if one edge has been deleted, it should be deleted in both direction + if (av[iv].del && aw[iw].del) continue; + + ///if this edge is an inexact edge + if(a->el == 0 && src[v>>1].is_fully_corrected == 1 && src[w>>1].is_fully_corrected == 1) { + if (kv > 1 && kw > 1) { + to_del = 1; + } else if (kw == 1) { + if (asg_topocut_aux(g, w^1, max_ext) < max_ext) to_del = 1; + } else if (kv == 1) { + if (asg_topocut_aux(g, v^1, max_ext) < max_ext) to_del = 1; + } + } + + if(a->el == 0 && src[v>>1].is_fully_corrected == 1 && src[w>>1].is_fully_corrected == 0) { + if(av[ov_max_i].el == 1 && src[av[ov_max_i].v>>1].is_fully_corrected == 1) { + if (kv > 1 && kw > 1) { + to_del = 1; + } else if (kw == 1) { + if (asg_topocut_aux(g, w^1, max_ext) < max_ext) to_del = 1; + } else if (kv == 1) { + if (asg_topocut_aux(g, v^1, max_ext) < max_ext) to_del = 1; + } + } + } + + if(((v>>1) == 50356 && (w>>1) == 1276292)||((v>>1) == 1276292 && (w>>1) == 50356)) { + fprintf(stderr, "[3]v->%u, w->%u, to_del->%u, el->%u, src[v>>1].is_fully_corrected->%u, src[w>>1].is_fully_corrected->%u\n", + v, w, to_del, a->el, src[v>>1].is_fully_corrected, src[w>>1].is_fully_corrected); + fprintf(stderr, "[4]v->%u, w->%u, to_del->%u, av[ov_max_i].v->%u, av[ov_max_i].el->%u, src[av[ov_max_i].v>>1].is_fully_corrected->%u\n", + v, w, to_del, av[ov_max_i].v, av[ov_max_i].el, src[av[ov_max_i].v>>1].is_fully_corrected); + } + + if (to_del) { + av[iv].del = aw[iw].del = 1, ++cnt; + if(dbg) { + kv_push(uint64_t, *dbg, (av - g->arc + iv)); + kv_push(uint64_t, *dbg, (aw - g->arc + iw)); + if(((v>>1) == 50356 && (w>>1) == 1276292)||((v>>1) == 1276292 && (w>>1) == 50356)) { + fprintf(stderr, "[5]v->%u, w->%u, (av-g->arc+iv)->%u, (aw-g->arc+iw)->%u\n", + v, w, (uint32_t)(av - g->arc + iv), (uint32_t)(aw - g->arc + iw)); + } + } + } + } + + if(!in) free(tx.a); + if(dbg) { + for (i = 0; i < dbg->n; i++) g->arc[dbg->a[i]].del = 0; + } + // if (cnt > 0) asg_cleanup(g); +} + +uint32_t trans_path_check(uint32_t a, uint32_t b, asg_t *g, ma_hit_t_alloc *rev, R_to_U* rI, +uint32_t minLen, asg64_v *t) +{ + if(a == b) return -1; + uint32_t tn = t->n, m, e, l[2] = {0}; + uint64_t *x[2]; + m = follow_limit_path(g, a, &e, &(l[0]), t, (uint32_t)-1); + if(m == LOOP || l[0] <= minLen) { + t->n = tn; return -1; + } + + m = follow_limit_path(g, b, &e, &(l[1]), t, (uint32_t)-1); + if(m == LOOP || l[1] <= minLen) { + t->n = tn; return -1; + } + + x[0] = t->a + tn; x[1] = t->a + tn + l[0]; + if(l[0] > l[1]) { + x[0] = t->a + tn + l[0]; x[1] = t->a + tn; + m = l[0]; l[0] = l[1]; l[1] = m; + } + assert(l[0]+l[1]+tn==t->n); + + uint32_t i, k, qi, ti, isU; double max_count = 0, min_count = 0; + for (i = 0; i < l[1]; i++) g->seq_vis[x[1][i]>>1] = 1; + for (i = 0; i < l[0]; i++) { + qi = x[0][i]>>1; + for (k = 0; k < rev[qi].length; k++) { + ti = Get_tn(rev[qi].buffer[k]); + if(g->seq[ti].del == 1) { + get_R_to_U(rI, ti, &ti, &isU); + if(ti == (uint32_t)-1 || isU == 1 || g->seq[ti].del == 1) continue; + } + min_count++; max_count += g->seq_vis[ti]; + } + } + + for (i = 0; i < l[1]; i++) g->seq_vis[x[1][i]>>1] = 0; + t->n = tn; + + if(min_count == 0) return -1; + if(max_count == 0) return 0; + if((max_count/min_count)>0.3) return 1; + return 0; +} + +void asg_arc_cut_length(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_trio, +uint32_t is_topo, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len) +{ + asg64_v tx = {0,0,0}, *b = NULL; + uint32_t i, k, v, w, n_vtx = g->n_seq<<1, nv, nw, kv, kw, trioF = (uint32_t)-1, ntrioF = (uint32_t)-1, ol_max, ou_max, to_del, cnt = 0, mm_ol, mm_ou; + asg_arc_t *av, *aw, *ve, *we, *vl_max, *wl_max; + + if(in) b = in; + else b = &tx; + b->n = 0; + + for (v = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + if(g->seq_vis[v] == 0) { + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + if (nv < 2) continue; + + for (i = kv = 0; i < nv; ++i) { + if(av[i].del) continue; + kv++; + } + if(kv < 2) continue; + + for (i = 0; i < nv; ++i) { + if(av[i].del) continue; + if(max_drop_len && av[i].ol >= (*max_drop_len)) continue; + kv_push(uint64_t, *b, (((uint64_t)av[i].ol)<<32) | ((uint64_t)(av-g->arc+i))); + } + } + } + + if(rev && rI) memset(g->seq_vis, 0, g->n_seq*2*sizeof(uint8_t)); + radix_sort_srt64(b->a, b->a + b->n); + for (k = 0; k < b->n; k++) { + if(g->arc[(uint32_t)b->a[k]].del) continue; + + v = g->arc[(uint32_t)b->a[k]].ul>>32; w = g->arc[(uint32_t)b->a[k]].v^1; + if(g->seq[v>>1].del || g->seq[w>>1].del) continue; + nv = asg_arc_n(g, v); nw = asg_arc_n(g, w); + av = asg_arc_a(g, v); aw = asg_arc_a(g, w); + if(nv<=1 && nw <= 1) continue; + + if(is_trio) { + if(get_arcs(g, v, NULL, 0)<=1 && get_arcs(g, w, NULL, 0)<=1) continue;///speedup + trioF = get_tip_trio_infor(g, v^1); + ntrioF = (trioF==FATHER? MOTHER : (trioF==MOTHER? FATHER : (uint32_t)-1)); + } + + ve = &(g->arc[(uint32_t)b->a[k]]); + for (i = 0; i < nw; ++i) { + if (aw[i].v == (v^1)) { + we = &(aw[i]); + break; + } + } + ///mm_ol and mm_ou are used to make edge with long indel more easy to be cutted + mm_ol = MIN(ve->ol, we->ol); mm_ou = MIN(ve->ou, we->ou); + + for (i = kv = ol_max = ou_max = 0, /**ve =**/ vl_max = NULL; i < nv; ++i) { + if(av[i].del) continue; + // if(av[i].v == (w^1)) ve = &(av[i]); + kv++; + if(is_trio && get_tip_trio_infor(g, av[i].v) == ntrioF) continue; + if(ol_max < av[i].ol) ol_max = av[i].ol, vl_max = &(av[i]); + if(ou_max < av[i].ou) ou_max = av[i].ou; + } + if (kv < 1) continue; + if (kv >= 2) { + if (/**ve->ol**/mm_ol > ol_max*len_rat) continue; + if (is_ou && /**ve->ou**/mm_ou > ou_max*ou_rat) continue; + } + + + for (i = kw = ol_max = ou_max = 0, /**we =**/ wl_max = NULL; i < nw; ++i) { + if(aw[i].del) continue; + // if(aw[i].v == (v^1)) we = &(aw[i]); + kw++; + if(is_trio && get_tip_trio_infor(g, aw[i].v) == ntrioF) continue; + if(ol_max < aw[i].ol) ol_max = aw[i].ol, wl_max = &(aw[i]); + if(ou_max < aw[i].ou) ou_max = aw[i].ou; + } + if (kw < 1) continue; + if (kw >= 2) { + if (/**we->ol**/mm_ol > ol_max*len_rat) continue; + if (is_ou && /**we->ou**/mm_ou > ou_max*ou_rat) continue; + } + + if (kv <= 1 && kw <= 1) continue; + + to_del = 0; + if(is_topo) { + if (kv > 1 && kw > 1) { + to_del = 1; + } else if (kw == 1) { + if (asg_topocut_aux(g, w^1, max_ext) < max_ext) to_del = 1; + } else if (kv == 1) { + if (asg_topocut_aux(g, v^1, max_ext) < max_ext) to_del = 1; + } + } + + if(rev && rI) { + if((to_del == 0) && vl_max && (ve->v!=vl_max->v) && (trans_path_check(ve->v, vl_max->v, g, rev, rI, max_ext, b)==0)) { + to_del = 1; + } + if((to_del == 0) && wl_max && (we->v!=wl_max->v) && (trans_path_check(we->v, wl_max->v, g, rev, rI, max_ext, b)==0)) { + to_del = 1; + } + if(vl_max && wl_max) assert(ve->v!=vl_max->v||we->v!=wl_max->v); + } + + + if (to_del) { + ve->del = we->del = 1, ++cnt; + } + } + // stats_sysm(g); + if(!in) free(tx.a); + if (cnt > 0) asg_cleanup(g); +} + +uint32_t if_false_bub_links(uint32_t v, asg_t *g, buf_t *x, asg64_v *b, uint32_t bs, int32_t check_dist) +{ + uint32_t i, mm = 1; + if (g->seq[v>>1].del) return 0; + for (i = bs; i < b->n; i++) { + g->arc[b->a[i]].del = 1; + asg_arc_del(g, g->arc[b->a[i]].v^1, (g->arc[b->a[i]].ul>>32)^1, 1); + } + if (asg_arc_n(g, v) < 2 || get_arcs(g, v, NULL, 0) < 2) mm = 0; + + if(mm) { + mm = 0; + if(asg_bub_pop1_primary_trio(g, NULL, v, check_dist, x, (uint32_t)-1, (uint32_t)-1, 0, + NULL, NULL, NULL, 0, 0, NULL)) { + + for (i = bs; i < b->n; i++) { + g->arc[b->a[i]].del = 0; + asg_arc_del(g, g->arc[b->a[i]].v^1, (g->arc[b->a[i]].ul>>32)^1, 0); + } + + asg_arc_t *av = asg_arc_a(g, v); uint32_t nv = asg_arc_n(g, v); + for (i = 0, b->n = bs; i < nv; i++) { + if (av[i].del) continue; + av[i].del = 1; asg_arc_del(g, av[i].v^1, (av[i].ul>>32)^1, 1); + kv_push(uint64_t, *b, ((uint64_t)(av-g->arc+i))); + } + + if(asg_bub_pop1_primary_trio(g, NULL, x->S.a[0]^1, check_dist, x, (uint32_t)-1, (uint32_t)-1, 0, NULL, NULL, NULL, 0, 0, NULL)) { + mm = 1; + } + } + } + + for (i = bs; i < b->n; i++) { + g->arc[b->a[i]].del = 0; + asg_arc_del(g, g->arc[b->a[i]].v^1, (g->arc[b->a[i]].ul>>32)^1, 0); + } + + return mm; +} + +void asg_arc_cut_bub_links(asg_t *g, asg64_v *in, float len_rat, float sec_len_rat, float ou_rat, uint32_t is_ou, uint64_t check_dist, ma_hit_t_alloc *rev, R_to_U* rI, int32_t max_ext) +{ + asg64_v tx = {0,0,0}, *b = NULL; + uint32_t v, w, t, k, i, n_vtx = g->n_seq<<1, nv, nw, kv, kw, kol, bn, me, mu, cnt = 0, sec_check; + asg_arc_t *av, *aw, *ref; + buf_t x; memset(&x, 0, sizeof(x)); x.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); + + if(in) b = in; + else b = &tx; + b->n = 0; + + for (v = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + if(g->seq_vis[v] == 0) { + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + if (nv < 2) continue; + + for (i = kv = kol = 0; i < nv; ++i) { + if(av[i].del) continue; + kv++; kol += av[i].ol; + } + if(kv < 2) continue;//must have at least one exact and one inexact + + kv_push(uint64_t, *b, ((((uint64_t)(kol))<<32) | v)); + } + } + + if(rev && rI) memset(g->seq_vis, 0, g->n_seq*2*sizeof(uint8_t)); + radix_sort_srt64(b->a, b->a + b->n); bn = b->n; + for (k = 0; k < bn; k++) { + v = (uint32_t)b->a[k]; + if (g->seq[v>>1].del) continue; + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + if (nv < 2 || get_arcs(g, v, NULL, 0) < 2) continue; + + for (i = 0, b->n = bn, sec_check = 0; i < nv; i++) { + if (av[i].del) continue; + + w = av[i].v^1; nw = asg_arc_n(g, w); aw = asg_arc_a(g, w); + if(nw < 2) break; + + for (t = kw = 0, me = mu = (uint32_t)-1; t < nw; t++) { + if(aw[t].del) continue; + kw++; + if(aw[t].v == (v^1)) continue;//note: me is the shortest edge except aw[t], so here is continue + if(aw[t].ol < me) me = aw[t].ol; + if(aw[t].ou < mu) mu = aw[t].ou; + kv_push(uint64_t, *b, ((uint64_t)(aw-g->arc+t))); + } + if(kw < 2) break; + + if(av[i].ol > me*len_rat && av[i].ol > me*sec_len_rat) break; + if(av[i].ol > me*len_rat) sec_check++; + if(is_ou && av[i].ou > mu*ou_rat) break; + } + + if(i < nv) continue; + + if(sec_check) { + for (i = 0, ref = NULL; i < nv; i++) {//forward + if (av[i].del) continue; + if(!ref) { + ref = &(av[i]); + } else { + if(trans_path_check(ref->v, av[i].v, g, rev, rI, max_ext, b)!=1) break; + } + } + + if(i < nv) { + if (b->n < bn + 2) continue;///less than two edges + for (i = bn, ref = NULL; i < b->n; i++) { + if(g->arc[b->a[i]].del) continue; + if(get_arcs(g, g->arc[b->a[i]].ul>>32, NULL, 0)!=2) break; + if(!ref) { + ref = &(g->arc[b->a[i]]); + } else { + if(trans_path_check(ref->v, g->arc[b->a[i]].v, g, rev, rI, max_ext, b)!=1) break; + } + } + if(i < b->n) continue; + } + } + + if(if_false_bub_links(v, g, &x, b, bn, check_dist)) { + for (i = 0; i < nv; ++i) { + if (av[i].del) continue; + av[i].del = 1; asg_arc_del(g, av[i].v^1, (av[i].ul>>32)^1, 1); + } + cnt++; + } + } + + // stats_sysm(g); + if(!in) free(tx.a); + free(x.a); free(x.S.a); free(x.T.a); free(x.b.a); free(x.e.a); + if(cnt > 0) asg_cleanup(g); +} + +void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float ou_rat, uint32_t is_ou, bub_label_t *b_mask_t) +{ + asg64_v tx = {0,0,0}, *b = NULL; + uint32_t v, w, t, k, i, n_vtx = g->n_seq<<1, nv, nw, kv, kw, kol, me, mu, cnt = 0, bn; + asg_arc_t *av, *aw; + + if(in) b = in; + else b = &tx; + b->n = 0; + + for (v = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + if(g->seq_vis[v] == 0) { + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + if (nv < 2) continue; + + for (i = kv = kol = 0; i < nv; ++i) { + if(av[i].del) continue; + kv++; kol += av[i].ol; + } + if(kv < 2) continue;//must have at least one exact and one inexact + + kv_push(uint64_t, *b, ((((uint64_t)(kol))<<32) | v)); + } + } + + radix_sort_srt64(b->a, b->a + b->n); bn = b->n; + for (k = 0; k < bn; k++) { + v = (uint32_t)b->a[k]; + if (g->seq[v>>1].del) continue; + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + if (nv < 2 || get_arcs(g, v, NULL, 0) < 2) continue; + + for (i = 0; i < nv; i++) { + if (av[i].del) continue; + + w = av[i].v^1; nw = asg_arc_n(g, w); aw = asg_arc_a(g, w); + if(nw < 2) break; + + for (t = kw = 0, me = mu = (uint32_t)-1; t < nw; t++) { + if(aw[t].del) continue; + kw++; + if(aw[t].v == (v^1)) continue;//note: me is the shortest edge except aw[t], so here is continue + if(aw[t].ol < me) me = aw[t].ol; + if(aw[t].ou < mu) mu = aw[t].ou; + } + if(kw < 2) break; + + if(av[i].ol > me*len_rat) break; + if(is_ou && av[i].ou > mu*ou_rat) break; + } + + if(i < nv) continue; + + for (i = 0; i < nv; ++i) { + if (av[i].del) continue; + av[i].del = 1; asg_arc_del(g, av[i].v^1, (av[i].ul>>32)^1, 1); + kv_push(uint64_t, *b, ((uint64_t)(av-g->arc+i))); + } + // b->a[cnt++] = v; + } + + if(b->n > bn) { + asg_arc_identify_simple_bubbles_multi(g, b_mask_t, 0); + for (k = bn, cnt = 0; k < b->n; k++) { + // if(g->arc[b->a[k]].del) continue; + v = g->arc[b->a[k]].ul>>32; w = g->arc[b->a[k]].v; + if(g->seq_vis[v] || g->seq_vis[v^1] || g->seq_vis[w] || g->seq_vis[w^1]) { + cnt++; continue; + } + g->arc[b->a[k]].del = 0; asg_arc_del(g, g->arc[b->a[k]].v^1, (g->arc[b->a[k]].ul>>32)^1, 0); + } + } + // stats_sysm(g); + if(!in) free(tx.a); + if(cnt > 0) asg_cleanup(g); +} + +#define LIM_LEN 100 + +uint32_t asg_cut_semi_circ(asg_t *g, uint32_t lim_len, uint32_t is_clean) +{ + uint32_t v, t, k, e, ss, i, n_vtx = g->n_seq<<1, nv, kv, nw, cnt = 0; + asg_arc_t *av, *aw; + + for (v = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + + av = asg_arc_a(g, v^1); nv = asg_arc_n(g, v^1); + if(nv <= 1) continue; + for (i = kv = 0; i < nv; ++i) { + if(av[i].del) continue; + kv++; if(kv > 1) break; + } + if(kv <= 1) continue; + + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + if(nv < 1) continue; + for (i = kv = 0; i < nv; ++i) { + if(av[i].del) continue; + kv++; if(kv > 1) break; + } + if(kv != 1) continue; + + for (i = 0; i < nv; ++i) { + if(av[i].del) continue; + t = follow_limit_path(g, v, &e, &ss, NULL, lim_len); + if(ss > lim_len || t == LONG_TIPS || t == LOOP || t == END_TIPS) break;//as kv == 1 + aw = asg_arc_a(g, v^1); nw = asg_arc_n(g, v^1); + for (k = 0; k < nw; k++) { + if (aw[k].del) continue; + if (aw[k].v == (e^1)) { + aw[k].del = 1; + asg_arc_del(g, aw[k].v^1, (aw[k].ul>>32)^1, 1); + cnt++; + } + } + break; //as kv == 1 + } + } + + if(cnt > 0 && is_clean) asg_cleanup(g); + return cnt; +} + +uint32_t asg_cut_chimeric_bub(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t normal_len, uint32_t is_clean) +{ + asg64_v tx = {0,0,0}, *b = NULL; + uint32_t v, w, nw, k, n_vtx = g->n_seq<<1, ei[2] = {0}, e, ss, cnt = 0; + asg_arc_t *aw; + if(in) b = in; + else b = &tx; + b->n = 0; + // fprintf(stderr, "[M::%s]\n", __func__); + for (v = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + ///note: ei[0] and ei[1] are the edge idx + if((get_arcs(g, v, &(ei[0]), 1)!=1) || (get_arcs(g, v^1, &(ei[1]), 1)!=1)) continue; + assert((g->arc[ei[0]].ul>>32) == v && (g->arc[ei[1]].ul>>32) == (v^1)); + if((get_arcs(g, g->arc[ei[0]].v^1, NULL, 0)!=2) || (get_arcs(g, g->arc[ei[1]].v^1, NULL, 0)!=2)) continue; + if(!if_sup_chimeric(&(src[v>>1]), g->seq[v>>1].len, b, 1)) continue; + w = g->arc[ei[0]].v^1; + aw = asg_arc_a(g, w); nw = asg_arc_n(g, w); + for (k = 0; k < nw; k++) { + if (aw[k].del) continue; + if (aw[k].v == (v^1)) { + ss = k; continue; + } + break; + } + // assert(aw[ss].v == (v^1));//this assert does not work, just ignore + if(follow_limit_path(g, aw[k].v, &e, &ss, NULL, (uint32_t)-1) != TWO_INPUT) continue; + if(ss > normal_len) { + w = e; + aw = asg_arc_a(g, w); nw = asg_arc_n(g, w); + for (k = ss = 0; k < nw; k++) { + if (aw[k].del) continue; + ss++; e = aw[k].v; + if(ss > 1) break; + } + if(ss == 1 && e == g->arc[ei[1]].v) asg_seq_del(g, v>>1), cnt++; + } + } + + if(!in) free(tx.a); + if (is_clean && cnt > 0) asg_cleanup(g); + return cnt; +} + +void asg_iterative_semi_circ(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t normal_len, uint32_t pop_chimer) +{ + uint64_t occ = 0, s = 1; + while (s) { + s = asg_cut_semi_circ(g, LIM_LEN, 0); + if(pop_chimer) s = s + asg_cut_chimeric_bub(g, src, in, normal_len, 0); + occ += s; + } + + // stats_sysm(g); + if(occ) asg_cleanup(g); +} + +uint32_t asg_cut_large_indel(asg_t *g, asg64_v *in, int32_t max_ext, float ou_rat, uint32_t is_ou) +{ + asg64_v tx = {0,0,0}, *b = NULL; + uint32_t v, w, n_vtx = g->n_seq<<1, i, k, kv, kw, nv, nw, ou_max, to_del, cnt = 0; + asg_arc_t *av, *aw, *ve, *we; + if(in) b = in; + else b = &tx; + b->n = 0; + + for (v = 0; v < n_vtx; ++v) { + if (g->seq[v>>1].del) continue; + if(g->seq_vis[v] == 0) { + av = asg_arc_a(g, v); nv = asg_arc_n(g, v); + if (nv < 2) continue; + + for (i = kv = 0; i < nv; ++i) { + if(av[i].del) continue; + kv++; + } + if(kv < 2) continue; + + for (i = 0; i < nv; ++i) { + if(av[i].del || av[i].no_l_indel) continue; + ///means there is a large indel at this edge + kv_push(uint64_t, *b, (((uint64_t)av[i].ol)<<32) | ((uint64_t)(av-g->arc+i))); + } + } + } + + radix_sort_srt64(b->a, b->a + b->n); + for (k = 0; k < b->n; k++) { + if(g->arc[(uint32_t)b->a[k]].del) continue; + + v = g->arc[(uint32_t)b->a[k]].ul>>32; w = g->arc[(uint32_t)b->a[k]].v^1; + if(g->seq[v>>1].del || g->seq[w>>1].del) continue; + nv = asg_arc_n(g, v); nw = asg_arc_n(g, w); + av = asg_arc_a(g, v); aw = asg_arc_a(g, w); + if(nv<=1 && nw <= 1) continue; + + for (i = kv = ou_max = 0, ve = NULL; i < nv; ++i) { + if(av[i].del) continue; + if(av[i].v == (w^1)) ve = &(av[i]); + kv++; + if(ou_max < av[i].ou) ou_max = av[i].ou; + } + if (kv < 1) continue; + if (kv >= 2) { + if (is_ou && ve->ou > ou_max*ou_rat) continue; + } + + + for (i = kw = ou_max = 0, we = NULL; i < nw; ++i) { + if(aw[i].del) continue; + if(aw[i].v == (v^1)) we = &(aw[i]); + kw++; + if(ou_max < aw[i].ou) ou_max = aw[i].ou; + } + if (kw < 1) continue; + if (kw >= 2) { + if (is_ou && we->ou > ou_max*ou_rat) continue; + } + + if (kv <= 1 && kw <= 1) continue; + + to_del = 0; + if (kv > 1 && kw > 1) { + to_del = 1; + } else if (kw == 1) { + if (asg_topocut_aux(g, w^1, max_ext) < max_ext) to_del = 1; + } else if (kv == 1) { + if (asg_topocut_aux(g, v^1, max_ext) < max_ext) to_del = 1; + } + + if (to_del) { + ve->del = we->del = 1, ++cnt; + } + } + // stats_sysm(g); + if(!in) free(tx.a); + if (cnt > 0) asg_cleanup(g); + return cnt; +} + +void debug_edges(asg64_v *dbg, uint32_t *l, uint32_t l_n) { + uint32_t k, k_n, i, m; + for (i = k = 0; i < l_n; i++) { + fprintf(stderr, "# gid-%u: %u\n", i, l[i]); + for (k_n = k + l[i]; k < k_n; k++) { + dbg->a[k] <<= 32; dbg->a[k] += i; + } + } + fprintf(stderr, "# dbg->n: %u\n", (uint32_t)dbg->n); + + radix_sort_srt64(dbg->a, dbg->a + dbg->n); + for (k = 1, i = 0; k <= dbg->n; k++) { + if(k == dbg->n || (dbg->a[k]>>32) != (dbg->a[i]>>32)) { + if(k - i < l_n) { + + for (m = i; m < k; m++) { + fprintf(stderr, "eid->%lu, gid->%u\n", dbg->a[m]>>32, (uint32_t)dbg->a[m]); + } + } + i = k; + } + } +} + +void ul_clean_gfa(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) +{ + #define HARD_OU_DROP 0.75 + #define HARD_OL_DROP 0.6 + #define HARD_OL_SEC_DROP 0.85 + #define HARD_ORTHOLOGY_DROP 0.4 + double step = + (clean_round==1?max_ovlp_drop_ratio:((max_ovlp_drop_ratio-min_ovlp_drop_ratio)/(clean_round-1))); + double drop = min_ovlp_drop_ratio; + int64_t i; asg64_v bu = {0,0,0}; uint32_t l_drop = 2000; + if(is_ou) update_sg_uo(sg, src); + + asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + for (i = 0; i < clean_round; i++, drop += step) { + if(drop > max_ovlp_drop_ratio) drop = max_ovlp_drop_ratio; + // stats_chimeric(sg, src, &bu); + if(!is_ou) asg_iterative_semi_circ(sg, src, &bu, max_tip, 1); + // fprintf(stderr, "(0):i->%ld, drop->%f\n", i, drop); + asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1); + asg_arc_cut_chimeric(sg, src, &bu); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + + asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); + + // print_edge(sg->arc+45471, "a"); + asg_arc_cut_inexact(sg, src, &bu, max_tip, is_ou, is_trio/**, NULL**//**&dbg**/); + // debug_edges(&dbg, d, 2); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + + asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1); + asg_arc_cut_length(sg, &bu, max_tip, drop, ou_drop_rate, is_ou, is_trio, 1, NULL, NULL, NULL); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + + asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1); + asg_arc_cut_bub_links(sg, &bu, HARD_OL_DROP, HARD_OL_SEC_DROP, HARD_OU_DROP, is_ou, asm_opt.large_pop_bubble_size, rev, rI, max_tip); + + asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1); + asg_arc_cut_complex_bub_links(sg, &bu, HARD_OL_DROP, HARD_OU_DROP, is_ou, b_mask_t); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + } + if(!is_ou) asg_iterative_semi_circ(sg, src, &bu, max_tip, 1); + + asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); + asg_cut_large_indel(sg, &bu, max_tip, HARD_OU_DROP, is_ou);///shoule we ignore ou here? + asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + + ///asg_arc_del_triangular_directly might be unnecessary + asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); + asg_arc_cut_length(sg, &bu, max_tip, HARD_ORTHOLOGY_DROP/**min_ovlp_drop_ratio**/, ou_drop_rate, is_ou, 0/**is_trio**/, 0, rev, rI, NULL); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + + asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0); + asg_arc_cut_length(sg, &bu, max_tip, min_ovlp_drop_ratio, ou_drop_rate, is_ou, 0/**is_trio**/, 0, rev, rI, &l_drop); + asg_arc_cut_tips(sg, max_tip, &bu, is_ou); + + if(!is_ou) asg_cut_semi_circ(sg, LIM_LEN, 1); + + free(bu.a); +} \ No newline at end of file diff --git a/gfa_ut.h b/gfa_ut.h new file mode 100644 index 0000000..17f2745 --- /dev/null +++ b/gfa_ut.h @@ -0,0 +1,18 @@ +#ifndef __GFA_UT__ +#define __GFA_UT__ +#include "Overlaps.h" + +void ul_clean_gfa(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 asg_arc_cut_tips(asg_t *g, uint32_t max_ext, asg64_v *in, uint32_t is_ou); +void asg_iterative_semi_circ(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t normal_len, uint32_t pop_chimer, asg64_v *dbg); +void asg_arc_cut_chimeric(asg_t *g, ma_hit_t_alloc* src, asg64_v *in); +void asg_arc_cut_inexact(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, int32_t max_ext, uint32_t is_ou, uint32_t is_trio); +void asg_arc_cut_length(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_trio, +uint32_t is_topo, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len); +void asg_arc_cut_bub_links(asg_t *g, asg64_v *in, float len_rat, float sec_len_rat, float ou_rat, uint32_t is_ou, uint64_t check_dist, ma_hit_t_alloc *rev, R_to_U* rI, int32_t max_ext); +void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float ou_rat, uint32_t is_ou, bub_label_t *b_mask_t); +uint32_t asg_cut_large_indel(asg_t *g, asg64_v *in, int32_t max_ext, float ou_rat, uint32_t is_ou); +uint32_t asg_cut_semi_circ(asg_t *g, uint32_t lim_len, uint32_t is_clean); + +#endif diff --git a/inter.cpp b/inter.cpp index caa05a7..ac31933 100644 --- a/inter.cpp +++ b/inter.cpp @@ -549,7 +549,7 @@ static inline int32_t comput_sc(const mg128_t *ai, const mg128_t *aj, int32_t ma ///p[]: id of last ///f[]: the score ending at i, not always the peak ///v[]: keeps the peak score up to i; -///t[]: id of next +///t[]: used for buffer ///min_cnt = 2; min_sc = 30; extra_u = 0 ///u = mg_chain_backtrack(n, f, p, v, t, min_cnt, min_sc, 0, &n_u, &n_v); uint64_t *mg_chain_backtrack(void *km, int64_t n, const int32_t *f, const int64_t *p, int32_t *v, int32_t *t, int32_t min_cnt, int32_t min_sc, int32_t extra_u, int32_t *n_u_, int32_t *n_v_) @@ -669,7 +669,7 @@ mg128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int // fill the score and backtrack arrays for (i = st = 0, max_ii = -1; i < n; ++i) { int64_t max_j = -1, end_j; - ///max_f -> minimizer span length in query, which is the initial score + ///max_f -> score of minimizer int32_t max_f = normal_sc(a[i].y>>MG_SEED_WT_SHIFT, a[i].y>>32&0xff), n_skip = 0; ///until we are at the same rid, same direction, and the coordinates are close enough while (st < i && (a[i].x>>32 != a[st].x>>32 || a[i].x > a[st].x + max_dist_x)) ++st; @@ -687,10 +687,11 @@ mg128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int if (++n_skip > max_skip) break; } - if (p[j] >= 0) t[p[j]] = i; + if (p[j] >= 0) t[p[j]] = i;//p[]: prefix idx; means there is a chain longer than 2 } - end_j = j;///end_j might be > 0 + end_j = j;///end_j might be > 0; just the end idx of backwards ///if not close enough, select a new max + ///max_ii is just used to rescue best-score in case best-score appears before end_j if (max_ii < 0 || (int64_t)(a[i].x - a[max_ii].x) > (int64_t)max_dist_x) {///select a new max int32_t max = INT32_MIN; max_ii = -1; @@ -699,14 +700,15 @@ mg128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int } ///note: it will happen when `max_ii` < `end_j`; ///iteration is terminated at `end_j` mostly because of `max_skip` and `max_iter` + ///max_ii is just used to rescue best-score in case best-score appears before end_j if (max_ii >= 0 && max_ii < end_j) { int32_t tmp; tmp = comput_sc(&a[i], &a[max_ii], max_dist_x, max_dist_y, bw, chn_pen_gap); if (tmp != INT32_MIN && max_f < tmp + f[max_ii]) max_f = tmp + f[max_ii], max_j = max_ii; } - // v[] keeps the peak score up to i; f[] is the score ending at i, not always the peak - f[i] = max_f, p[i] = max_j; + // v[] keeps the peak score up to i (as score might decerase); f[] is the score ending at i, not always the peak + f[i] = max_f, p[i] = max_j;//p[]: prefix idx v[i] = max_j >= 0 && v[max_j] > max_f? v[max_j] : max_f; if (max_ii < 0 || ((int64_t)(a[i].x - a[max_ii].x) <= (int64_t)max_dist_x && f[max_ii] < f[i])) max_ii = i; @@ -719,7 +721,7 @@ mg128_t *mg_lchain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int kfree(km, a); kfree(km, v); return 0; } - //u[]: sc|occ of chains + //u[]: sc|occ of chains; chain is mostly sorted by the score; at least the first chain has the largest score //v[]: idx of each element return compact_a(km, n_u, u, n_v, v, a); } @@ -933,7 +935,7 @@ mg_pathv_t *mg_shortest_k(void *km0, const asg_t *g, uint32_t src, int32_t n_dst ///h_seeds+seeds+n_seeds ----> hash index of qs[0, ql) } **/ - + ///dst is how many candidates KCALLOC(km, dst_done, n_dst); KMALLOC(km, dst_group, n_dst); // multiple dst[] may have the same dst[].v. We need to group them first. @@ -978,12 +980,14 @@ mg_pathv_t *mg_shortest_k(void *km0, const asg_t *g, uint32_t src, int32_t n_dst ///each src corresponds to one node in the hash table , but corresponds to node in the AVL tree k = kh_put(sp, h, src, &absent);///here is a hash table q = &kh_val(h, k); + ///for normal graph traversal, one node just has one parental node; here each node has at most 16 parental nodes q->k = 1, q->p[0] = p, q->mlen = 0, q->qs = q->qe = -1; n_done = 0; ///the key of avl tree: #define sp_node_cmp(a, b) (((a)->di > (b)->di) - ((a)->di < (b)->di)) ///the higher bits of (*)->di is distance to src node ///so the key of avl tree is distance + ///in avl tree , one node might be saved multipe times while (kavl_size(head, root) > 0) {///thr first root is src int32_t i, nv; asg_arc_t *av; @@ -1013,10 +1017,10 @@ mg_pathv_t *mg_shortest_k(void *km0, const asg_t *g, uint32_t src, int32_t n_dst //src can reach ref id r->v; there might be not only one alignment chain in r->v //so we need to scan all of them for (j = 0; j < cnt; ++j) { - mg_path_dst_t *t = &dst[(int32_t)dst_group[off + j]]; + mg_path_dst_t *t = &dst[(int32_t)dst_group[off + j]];///t is a linear alignment at r->v int32_t done = 0; ///the src and dest are at the same ref id, say we directly find the shortest path - if (t->inner) { + if (t->inner) {//usually the first node, which is same to src done = 1; } else { int32_t mlen = 0, copy = 0; @@ -1028,17 +1032,20 @@ mg_pathv_t *mg_shortest_k(void *km0, const asg_t *g, uint32_t src, int32_t n_dst //if (mg_dbg_flag & MG_DBG_GC1) fprintf(stderr, " src=%c%s[%d],qlen=%d\tdst=%c%s[%d]\ttarget_distx=%d,target_hash=%x\tdistx=%d,mlen=%d,hash=%x\n", "><"[src&1], g->seg[src>>1].name, src, ql, "><"[t->v&1], g->seg[t->v>>1].name, t->v, t->target_dist - g->seg[src>>1].len, t->target_hash, dist - g->seg[src>>1].len, mlen, r->hash); // note: t indicates a linear alignmnet, instead of a node in graph ///target_dist should be the distance on query - if (t->n_path == 0) { // keep the shortest path + if (t->n_path == 0) { // means this alignment has never been visited before; keep the shortest path anyway copy = 1; - } else if (t->target_dist >= 0) { // we have a target distance; choose the closest - if (dist == t->target_dist && t->check_hash && r->hash == t->target_hash) { // we found the target path + // we have a target distance; choose the closest; + // there is already several paths reaching the linear alignment + } else if (t->target_dist >= 0) { + // we found the target path; hash is the path hash including multiple nodes, instead of node hash + if (dist == t->target_dist && t->check_hash && r->hash == t->target_hash) { copy = 1, done = 1; } else { int32_t d0 = t->dist, d1 = dist; d0 = d0 > t->target_dist? d0 - t->target_dist : t->target_dist - d0; d1 = d1 > t->target_dist? d1 - t->target_dist : t->target_dist - d1; ///if the new distance (d1) is smaller than the old distance (d0), update the results - ///in other words, the length of new path should be closer to t->target_dist + ///the length of new path should be closer to t->target_dist if (d1 - mlen/2 < d0 - t->mlen/2) copy = 1; } } @@ -2108,7 +2115,7 @@ st_mt_t *sp, mg_tbuf_t *b, int32_t w, int32_t k, int32_t hpc, int32_t mz_sd, int max_chain_gap_qry = max_chain_gap_ref = opt->max_gap; **/ max_chain_gap_qry = max_chain_gap_ref = qlen*2; - if (n_a == 0) { + if (n_a == 0) {//no matched minimizer if(a) kfree(b->km, a); a = 0, n_lc = 0, u = 0; } else { @@ -2116,8 +2123,8 @@ st_mt_t *sp, mg_tbuf_t *b, int32_t w, int32_t k, int32_t hpc, int32_t mz_sd, int opt->min_lc_cnt, opt->min_lc_score, opt->chn_pen_gap, n_a, a, &n_lc, &u, b->km); } - if (n_lc) {///n_lc is how many chain we found - lc = mg_lchain_gen(b->km, qlen, n_lc, u, a, ug); + if (n_lc) {///n_lc is how many linear chain we found + lc = mg_lchain_gen(b->km, qlen, n_lc, u, a, ug);//lc->the status of each chain; u->idx of each chain; for (i = 0; i < n_lc; ++i)///update a[] since ref_id|rev has already been saved to lc[].v mg_update_anchors(lc[i].cnt, &a[lc[i].off], n_mini_pos, mini_pos);///update a[].x } else lc = 0; @@ -2526,11 +2533,11 @@ double es_win_err(overlap_region* o, int64_t winLen, int64_t s, int64_t e) fprintf(stderr, "WARNNING-1, o->w_list_length->%u, o->x_id->%u, s->%ld, e->%ld, w_list_s->%lu, w_list_e->%lu, winLen->%ld, o->x_pos_s->%u, o->x_pos_e->%u, si->%ld, flag->%d\n", o->w_list_length, o->x_id, s, e, o->w_list[k].x_start, o->w_list[k].x_end, winLen, o->x_pos_s, o->x_pos_e, si, o->w_list[k].y_end); } - tLen += o->w_list[k].x_end+1-o->w_list[k].x_start; + tLen += ov/**o->w_list[k].x_end+1-o->w_list[k].x_start**/; if(o->w_list[k].y_end != -1) { tErr += (ov*o->w_list[k].error)/(o->w_list[k].x_end+1-o->w_list[k].x_start); } else { - tErr += o->w_list[k].x_end+1-o->w_list[k].x_start; + tErr += ov/**o->w_list[k].x_end+1-o->w_list[k].x_start**/; } k = ei; @@ -2540,11 +2547,11 @@ double es_win_err(overlap_region* o, int64_t winLen, int64_t s, int64_t e) fprintf(stderr, "WARNNING-2, o->w_list_length->%u, o->x_id->%u, s->%ld, e->%ld, w_list_s->%lu, w_list_e->%lu, winLen->%ld, o->x_pos_s->%u, o->x_pos_e->%u, ei->%ld, flag->%d\n", o->w_list_length, o->x_id, s, e, o->w_list[k].x_start, o->w_list[k].x_end, winLen, o->x_pos_s, o->x_pos_e, ei, o->w_list[k].y_end); } - tLen += o->w_list[k].x_end+1-o->w_list[k].x_start; + tLen += ov/**o->w_list[k].x_end+1-o->w_list[k].x_start**/; if(o->w_list[k].y_end != -1) { tErr += (ov*o->w_list[k].error)/(o->w_list[k].x_end+1-o->w_list[k].x_start); } else { - tErr += o->w_list[k].x_end+1-o->w_list[k].x_start; + tErr += ov/**o->w_list[k].x_end+1-o->w_list[k].x_start**/; } return ((double)tErr)/((double)tLen); @@ -2603,15 +2610,36 @@ int64_t gen_contain_chain(const ul_idx_t *uref, utg_ct_t *p, overlap_region* o, return 1; } -int64_t debug_utg_ct_t(const ul_idx_t *uref, overlap_region* o, utg_ct_t *ct_a, int64_t ct_n,haplotype_evdience *he_a, int64_t he_n) +int64_t debug_utg_ct_t(const ul_idx_t *uref, overlap_region* o, utg_ct_t *ct_a, int64_t ct_n, ma_utg_t *u, utg_ct_t *z, haplotype_evdience *he_a, int64_t he_n) { - int64_t k, i, ss, m = 0; + int64_t k, i, l, rs, re, ss, m = 0; utg_ct_t *p = NULL; - for (i = 0; i < ct_n; i++) { - p = &(ct_a[i]); + if(ct_a && ct_n) { + for (i = 0; i < ct_n; i++) { + p = &(ct_a[i]); + for (k = 0; k < he_n; k++) { + ss = o->y_pos_strand?uref->ug->u.a[o->y_id].len - he_a[k].cov - 1:he_a[k].cov; + if(ss >= p->s && ss < p->e) break; + } + if(k < he_n) m++; + } + } + if(u) { + for (i = l = 0; i < u->n; i++) { + rs = l; re = l + Get_READ_LENGTH(R_INF, (u->a[i]>>33)); + l += (uint32_t)u->a[i]; + for (k = 0; k < he_n; k++) { + ss = o->y_pos_strand?uref->ug->u.a[o->y_id].len - he_a[k].cov - 1:he_a[k].cov; + if(ss >= rs && ss < re) break; + } + if(k < he_n) m++; + } + } + + if(z) { for (k = 0; k < he_n; k++) { ss = o->y_pos_strand?uref->ug->u.a[o->y_id].len - he_a[k].cov - 1:he_a[k].cov; - if(ss >= p->s && ss < p->e) break;; + if(ss >= z->s && ss < z->e) break; } if(k < he_n) m++; } @@ -2680,6 +2708,63 @@ kv_ul_ov_t *chains, double diff_ec_ul, int64_t winLen, void *km) return t0; } +int64_t rescue_trans_ul_chains(const ul_idx_t *uref, overlap_region* o, haplotype_evdience *he_a, int64_t he_n, ma_utg_t *u, +kv_ul_ov_t *chains, double diff_ec_ul, int64_t winLen, void *km) +{ + uint64_t ys, ye, i, l; + int64_t k, ff, ss, t0 = 0; + utg_ct_t p; + + if(o->y_pos_strand == 0) { + ys = o->y_pos_s; ye = o->y_pos_e + 1; + for (i = k = l = 0; i < u->n; i++) { + p.x = u->a[i]>>32; p.s = l; p.e = l + Get_READ_LENGTH(R_INF, (u->a[i]>>33)); + l += (uint32_t)u->a[i]; + if(p.e <= ys) continue; + if(p.s >= ye) break; + for (ff = 1; k < he_n; k++) { + if(he_a[k].cov >= p.s && he_a[k].cov < p.e) { + ff = 0; + break; + } + if(he_a[k].cov >= p.e) break; + } + // if(ff == debug_utg_ct_t(uref, o, 0, 0, 0, &p, he_a, he_n)) fprintf(stderr, "ERROR\n"); + if(ff) { + ///push ovlp + t0 += gen_contain_chain(uref, &p, o, chains, diff_ec_ul, winLen, km); + } + // if(!ff) t0++; + } + } else { + ys = uref->ug->u.a[o->y_id].len - (o->y_pos_e+1); + ye = uref->ug->u.a[o->y_id].len - o->y_pos_s; + for (i = l = 0, k = he_n - 1; i < u->n; i++) { + p.x = u->a[i]>>32; p.s = l; p.e = l + Get_READ_LENGTH(R_INF, (u->a[i]>>33)); + l += (uint32_t)u->a[i]; + if(p.e <= ys) continue; + if(p.s >= ye) break; + for (ff = 1; k >= 0; k--) { + ss = uref->ug->u.a[o->y_id].len - he_a[k].cov - 1; + if(ss >= p.s && ss < p.e) { + ff = 0; + break; + } + if(ss >= p.e) break; + } + // if(ff == debug_utg_ct_t(uref, o, 0, 0, 0, &p, he_a, he_n)) fprintf(stderr, "ERROR\n"); + if(ff) { + ///push ovlp + t0 += gen_contain_chain(uref, &p, o, chains, diff_ec_ul, winLen, km); + } + // if(!ff) t0++; + } + } + // if(debug_utg_ct_t(uref, o, NULL, 0, u, he_a, he_n)!=t0) fprintf(stderr, "ERROR\n"); + // fprintf(stderr, "t0->%ld\n", t0); + return t0; +} + int64_t dedup_sort_ul_ov_t(ul_ov_t *a, int64_t a_n) { int64_t k, l, z, r, i; @@ -3069,6 +3154,132 @@ void fill_edge_weight(ul_ov_t *a, int64_t a_n, const ug_opt_t *uopt, int64_t bw, } **/ +int64_t get_ecov_adv(const ul_idx_t *uref, const ug_opt_t *uopt, uint32_t v, uint32_t w, int64_t qlen, int64_t bw, double diff_ec_ul, int64_t dq) +{ + int64_t dt = -1, dif, mm; + const asg_t *g = uref?uref->ug->g:NULL; + uint32_t nv, i; asg_arc_t *av = NULL; + if(g) { + nv = asg_arc_n(g, v); av = asg_arc_a(g, v); + for (i = 0; i < nv; i++) { + if(av[i].del || av[i].v != w) continue; + dt = av[i].ol; + break; + } + } + + if(dt < 0 && uopt) { + ma_hit_t_alloc* src = uopt->sources; + int64_t min_ovlp = uopt->min_ovlp; + int64_t max_hang = uopt->max_hang; + uint64_t z, qn, tn, x = v>>1; int32_t r = 1; asg_arc_t e; + for (z = 0; z < src[x].length; z++) { + qn = Get_qn(src[x].buffer[z]); tn = Get_tn(src[x].buffer[z]); + if(tn != (w>>1)) continue; + r = ma_hit2arc(&(src[x].buffer[z]), Get_READ_LENGTH(R_INF, qn), Get_READ_LENGTH(R_INF, tn), max_hang, asm_opt.max_hang_rate, min_ovlp, &e); + if(r < 0) { + if(r == MA_HT_QCONT || r == MA_HT_TCONT) { + if(src[x].buffer[z].rev == ((uint32_t)(v^w))) { + dt = Get_qe(src[x].buffer[z]) - Get_qs(src[x].buffer[z]); + if(dt < Get_te(src[x].buffer[z]) - Get_ts(src[x].buffer[z])) { + dt = Get_te(src[x].buffer[z]) - Get_ts(src[x].buffer[z]); + } + break; + } + } + continue; + } + if((e.ul>>32) != v || e.v != w) continue; + dt = e.ol; + break; + } + } + if(dt < 0) return 0; + dif = (dq>dt? dq-dt:dt-dq); + mm = MAX(dq, dt); mm *= diff_ec_ul; if(mm < bw) mm = bw; + // if((v>>1) == 1163 && (w>>1) == 1168) fprintf(stderr, ">>>>>>dis_q:%ld, dis_t:%ld, dif:%ld, mm:%ld\n", dis_q, dis_t, dif, mm); + if(dif <= mm) return 1; + return 0; +} + +int64_t gl_chain_advance(kv_ul_ov_t *res, kv_ul_ov_t *ex, const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, +double diff_ec_ul, int64_t qlen, uint64_t *srt, uint64_t *idx, uint64_t *track, uint64_t extend_check, void *km) +{ + uint32_t li_v, lj_v, rev_n; + int64_t mm_ovlp, x, i, j, k, sc, csc, mm_sc, mm_idx, qo; + ul_ov_t *li = NULL, *lj = NULL, rev_t; + radix_sort_ul_ov_srt_qe(res->a, res->a + res->n); + for (i = 0; i < (int64_t)res->n; ++i) { + li = &(res->a[i]); li_v = (li->tn<<1)|li->rev; + mm_ovlp = uref?max_ovlp(uref->ug->g, li_v^1):max_ovlp_src(uopt, li_v^1); + x = (li->qs + mm_ovlp)*diff_ec_ul; + if(x < bw) x = bw; + x += li->qs + mm_ovlp; + if (x > qlen+1) x = qlen+1; + x = find_ul_ov_max(i, res->a, x); + csc = uref?retrieve_u_cov_region(uref, li->tn, 0, li->ts, li->te, NULL):li->te-li->ts; + mm_sc = csc; mm_idx = -1; + for (j = x; j >= 0; --j) { // collect potential destination vertices + lj = &(res->a[j]); lj_v = (lj->tn<<1)|lj->rev; + if((!extend_check) && (lj->qe <= li->qs)) break; //evan this pair has a overlap, its length will be very small; just ignore + // if(lj->qs >= li->qs) continue; // lj is contained in li on the query coordinate + qo = infer_rovlp(li, lj, NULL, NULL); ///overlap length in query (UL read) + if(li_v != lj_v && get_ecov_adv(uref, uopt, li_v^1, lj_v^1, qlen, bw, diff_ec_ul, qo)) { + sc = csc + (track[j]>>32); + if(sc > mm_sc) mm_sc = sc, mm_idx = j; + } + } + // 4294967295L + track[i] = mm_sc; track[i] <<= 32; + track[i] |= (mm_idx>=0?mm_idx:((uint64_t)0x7FFFFFFF)); + srt[i] = mm_sc; srt[i] <<= 32; srt[i] |= i; + // fprintf(stderr, "+++i:%ld, mm_idx:%ld, mm_sc:%ld\n", i, mm_idx, mm_sc); + // fprintf(stderr, "[M::utg%.6d%c] qs->%u; qe->%u\n\n", li->tn+1, "lc"[uref->ug->u.a[li->tn].circ], li->qs, li->qe); + } + + + int64_t n_v, n_u, n_v0; + radix_sort_gfa64(srt, srt+res->n); ex->n = res->n; + for (k = (int64_t)res->n-1, n_v = n_u = 0; k >= 0; --k) { + // fprintf(stderr, "\nk:%ld\n", k); + n_v0 = n_v; + for (i = (uint32_t)srt[k]; i >= 0 && (track[i]&((uint64_t)0x80000000)) == 0;) { + ex->a[n_v++] = res->a[i]; track[i] |= ((uint64_t)0x80000000); + // fprintf(stderr, "+i:%ld, ", i); + // fprintf(stderr, "[M::utg%.6d%c] qs->%u; qe->%u\n", res->a[i].tn+1, "lc"[uref->ug->u.a[res->a[i].tn].circ], res->a[i].qs, res->a[i].qe); + if((track[i]&((uint64_t)0x7FFFFFFF)) == ((uint64_t)0x7FFFFFFF)) i = -1; + else i = track[i]&((uint64_t)0x7FFFFFFF); + } + if(n_v0 == n_v) continue; + ///keep the whole score; do not cut score like minigraph + sc = (i<0?(srt[k]>>32):((srt[k]>>32)-(track[i]>>32))); + // sc = srt[k]>>32; + idx[n_u++] = ((uint64_t)sc<<32)|(n_v-n_v0); + } + + for (k = 0, n_v = n_v0 = 0; k < n_u; k++) { + n_v0 = n_v; n_v += (uint32_t)idx[k]; + res->a[k].qn = idx[k]>>32; + res->a[k].ts = n_v0; res->a[k].te = n_v; + + rev_n = ((uint32_t)idx[k])>>1; + ///we need to consider contained reads; so determining qs is not such easy + res->a[k].qs = (uint32_t)-1; res->a[k].qe = ex->a[n_v0].qe; + for (i = 0; i < rev_n; i++) { + rev_t = ex->a[n_v0+i]; + ex->a[n_v0+i] = ex->a[n_v0+rev_n-i-1]; + ex->a[n_v0+rev_n-i-1] = rev_t; + if(res->a[k].qs > ex->a[n_v0+i].qs) res->a[k].qs = ex->a[n_v0+i].qs; + if(res->a[k].qs > ex->a[n_v0+rev_n-i-1].qs) res->a[k].qs = ex->a[n_v0+rev_n-i-1].qs; + } + if(i < ((uint32_t)idx[k]) && res->a[k].qs < ex->a[n_v0+i].qs) { + res->a[k].qs = ex->a[n_v0+i].qs; + } + } + res->n = n_u; + return res->n; +} + int64_t gl_chain_refine_advance(overlap_region_alloc* olist, Correct_dumy* dumy, haplotype_evdience_alloc *hap, glchain_t *ll, const ul_idx_t *uref, double diff_ec_ul, int64_t winLen, int64_t qlen, const ug_opt_t *uopt, void *km) { @@ -3078,20 +3289,29 @@ void *km) gl_chain_gen(olist, uref, idx, km); if(idx->n == 0) return 0; - uint64_t k, an, cn, si = 0, ei = 0, resc = 0, idx_pl = idx->n; + uint64_t k, an, cn, si = 0, ei = 0, resc = 0, resc_tk = 0, idx_pl = idx->n; + ma_utg_t *u = NULL; for (k = 0; k < olist->length; k++) { if(olist->list[k].is_match!=2) continue; - cn = ((uint32_t)(ct->idx.a[olist->list[k].y_id])); - if(cn==0) continue; - an = update_ava_het_site(hap, k, &si, &ei, cn); + an = update_ava_het_site(hap, k, &si, &ei, 1); // if(an != get_het_site(hap, k)) fprintf(stderr, "an->%lu, get_het_site->%lu\n", an, get_het_site(hap, k)); - - if(cn > 0 && an > 0) { + if(an == 0) { + fprintf(stderr, "ERROR\n"); + continue; + } + cn = ((uint32_t)(ct->idx.a[olist->list[k].y_id])); + if(cn > 0) { resc += rescue_contain_ul_chains(uref, &(olist->list[k]), hap->list+si, an, ct->rids.a + ((ct->idx.a[olist->list[k].y_id])>>32), cn, idx, diff_ec_ul, winLen, km); } + u = &(uref->ug->u.a[olist->list[k].y_id]); + if(u->n > 1) {///no redundant items here, + resc_tk += rescue_trans_ul_chains(uref, &(olist->list[k]), hap->list+si, an, u, + &(ll->tk), diff_ec_ul, winLen, km); + } + si = ei; } @@ -3106,6 +3326,14 @@ void *km) // } // } } + + // fprintf(stderr, "resc_tk->%ld\n", resc_tk); + if(resc_tk > 0) { + for (k = ll->tk.n - resc_tk; k < ll->tk.n; k++) { + kv_push_km(km, ul_ov_t, *idx, ll->tk.a[k]); + } + ll->tk.n -= resc_tk; + } if(idx->n > 0) { an = infer_read_ovlp(uref, olist, idx, &(ll->tk), diff_ec_ul, winLen, uopt, ct, km); @@ -3587,7 +3815,10 @@ void determine_connective(all_ul_t *m, const ug_opt_t *uopt, int64_t bw, double lk = &(p->bb.a[k]); lk_v = (((uint32_t)(lk->hid))<<1)|((uint32_t)(lk->rev)); if(lk->qe <= li->qs) break;//evan this pair has a overlap, its length will be very small; just ignore if((li_v == lk_v) || (lk->hid&m->mm)) continue; - if(li->qs <= 0) continue;///means the UL read does not longer than the overlap between li and lk + // if(li->qs <= 0) continue;///means the UL read does not longer than the overlap between li and lk + // if(lk->qs <= 0) continue;//the UL read should be cover the whole HiFi reads li and lk + if(((li->te - li->ts)*1.05) < Get_READ_LENGTH(R_INF, li->hid)) continue; + if(((lk->te - lk->ts)*1.05) < Get_READ_LENGTH(R_INF, lk->hid)) continue; x = /**((int64_t)(lk->qe))-((int64_t)(li->qs))**/infer_rovlp(NULL, NULL, li, lk); t = query_ovlp_src(uopt, li_v^1, lk_v^1, x, diff_ec_ul, &ol); if(t) { @@ -3612,6 +3843,13 @@ static void update_ovlp_src(void *data, long i, int tid) // callback for kt_for( } } +uint64_t* get_hifi2ul_list(all_ul_t *x, uint64_t hid, uint64_t* a_n) +{ + (*a_n) = x->ridx.idx.a[hid+1] - x->ridx.idx.a[hid]; + return x->ridx.occ.a + x->ridx.idx.a[hid];; +} + + static void update_ovlp_src_bl(void *data, long i, int tid) { uldat_t *sl = (uldat_t *)data; @@ -4364,27 +4602,74 @@ ul_contain *ul_contain_gen(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc* src, int64_t return p; } -cvert_t *cvert_t_gen(const ug_opt_t *uopt) + +void append_inexact_edges(ma_ug_t *ug, const ug_opt_t *uopt, asg_t *rg) { - uint64_t i, k; ma_utg_t *u = NULL; - cvert_t *p = NULL; CALLOC(p, 1); - p->rg = build_init_sg(uopt->sources, uopt->reverse_sources, R_INF.total_reads, uopt->min_dp, uopt->readLen, - uopt->min_ovlp, uopt->max_hang, uopt->coverage_cut, uopt->ruIndex); - p->ug = ma_ug_gen(p->rg); - MALLOC(p->idx, p->rg->n_seq); memset(p->idx, -1, sizeof((*(p->idx)))*p->rg->n_seq); - for (i = 0; i < p->ug->u.n; i++) { - u = &(p->ug->u.a[i]); - for (k = 0; k < u->n; k++) { - if(p->idx[u->a[k]>>33] == ((uint64_t)-1)) { - p->idx[u->a[k]>>33] = i; - p->idx[u->a[k]>>33] <<= 32; - p->idx[u->a[k]>>33] |= k; - } else { - p->idx[u->a[k]>>33] = ((uint32_t)-1); - } + uint32_t *idx = NULL, n_read = R_INF.total_reads, z, v, k, qn, tn, tu, ut_v, ut_w; + ma_utg_t *u = NULL; ma_hit_t_alloc *src = uopt->sources, *s = NULL; + int32_t r; asg_arc_t t, *p = NULL; + int64_t min_ovlp = uopt->min_ovlp, max_hang = uopt->max_hang; + + MALLOC(idx, n_read); memset(idx, -1, n_read*sizeof(*(idx))); + for (z = 0; z < ug->u.n; z++) { + u = &(ug->u.a[z]); + if(u->circ) continue; + idx[u->start>>1] = idx[u->end>>1] = z; + } + + for (z = 0; z < ug->u.n; z++) { + u = &(ug->u.a[z]); + if(u->circ) continue; + + v = u->start^1; s = &(src[v>>1]); ut_v = (z<<1); + for (k = 0; k < s->length; k++) { + if(s->buffer[k].el) continue;///we just need inexact edges + qn = Get_qn(s->buffer[k]); tn = Get_tn(s->buffer[k]); tu = idx[tn]; ut_w = (uint32_t)-1; + if(tu == (uint32_t)-1 || ug->g->seq[tu].del) continue; + if((Get_qe(s->buffer[k]) - Get_qs(s->buffer[k])) < min_ovlp) continue; + if((Get_te(s->buffer[k]) - Get_ts(s->buffer[k])) < min_ovlp) continue; + r = ma_hit2arc(&(s->buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + if(r < 0 || (t.ul>>32) != v) continue; + if(t.v == ug->u.a[tu].start) ut_w = tu<<1; + if(t.v == ug->u.a[tu].end) ut_w = (tu<<1)+1; + p = asg_arc_pushp(ug->g); + *p = t; p->ul = ut_v; p->ul <<= 32; p->ul += ((uint32_t)(t.ul)); p->v = ut_w; + } + + v = u->end^1; s = &(src[v>>1]); ut_v = (z<<1) + 1; + for (k = 0; k < s->length; k++) { + if(s->buffer[k].el) continue;///we just need inexact edges + qn = Get_qn(s->buffer[k]); tn = Get_tn(s->buffer[k]); tu = idx[tn]; ut_w = (uint32_t)-1; + if(tu == (uint32_t)-1 || ug->g->seq[tu].del) continue; + if((Get_qe(s->buffer[k]) - Get_qs(s->buffer[k])) < min_ovlp) continue; + if((Get_te(s->buffer[k]) - Get_ts(s->buffer[k])) < min_ovlp) continue; + r = ma_hit2arc(&(s->buffer[k]), rg->seq[qn].len, rg->seq[tn].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); + if(r < 0 || (t.ul>>32) != v) continue; + if(t.v == ug->u.a[tu].start) ut_w = tu<<1; + if(t.v == ug->u.a[tu].end) ut_w = (tu<<1)+1; + p = asg_arc_pushp(ug->g); + *p = t; p->ul = ut_v; p->ul <<= 32; p->ul += ((uint32_t)(t.ul)); p->v = ut_w; } } - return p; + + asg_cleanup(ug->g); + + uint32_t w, nv, n_asymm = 0; asg_arc_t *av = NULL; + for (z = 0; z < ug->g->n_arc; ++z) { + if(ug->g->arc[z].del) continue; + v = ug->g->arc[z].v^1; w = ug->g->arc[z].ul>>32^1; + nv = asg_arc_n(ug->g, v); av = asg_arc_a(ug->g, v); + for (k = 0; k < nv; ++k) + if ((!av[k].del) && av[k].v == w) break; + if (k == nv) ug->g->arc[z].del = 1, ++n_asymm; + } + + if(n_asymm) { + asg_cleanup(ug->g); + fprintf(stderr, "[M::%s::] # asymm edges: %u\n", __func__, n_asymm); + } + + free(idx); } ul_idx_t *dedup_HiFis(const ug_opt_t *uopt) diff --git a/inter.h b/inter.h index d15cce8..3e4e344 100644 --- a/inter.h +++ b/inter.h @@ -1,8 +1,10 @@ #ifndef __INTER__ #define __INTER__ #include "Overlaps.h" +#include "Process_Read.h" void ul_resolve(ma_ug_t *ug, const asg_t *rg, const ug_opt_t *uopt, int hap_n); void ul_load(const ug_opt_t *uopt); +uint64_t* get_hifi2ul_list(all_ul_t *x, uint64_t hid, uint64_t* a_n); #endif