diff --git a/CommandLines.h b/CommandLines.h index 0cb1a4d..6099955 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -4,7 +4,7 @@ #include #include -#define HA_VERSION "0.15.1-r336" +#define HA_VERSION "0.15.3-r339" #define VERBOSE 0 diff --git a/Overlaps.cpp b/Overlaps.cpp index 2c4ebbf..79fddc5 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -13295,9 +13295,9 @@ long long gap_fuzz, bub_label_t* b_mask_t) if(t_ch) destory_trans_chain(&t_ch); - // gfa_name = (char*)malloc(strlen(output_file_name)+25); + // char* gfa_name = (char*)malloc(strlen(output_file_name)+50); // sprintf(gfa_name, "%s.after.clean_d_utg.noseq.gfa", output_file_name); - // output_file = fopen(gfa_name, "w"); + // FILE* output_file = fopen(gfa_name, "w"); // ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file); // fclose(output_file); // free(gfa_name); @@ -16308,7 +16308,7 @@ uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate) if(is_final_check) { - if(non_flag_occ < u->n*double_check_rate) continue; + /**if(non_flag_occ < u->n*double_check_rate) continues**/; } else { @@ -16361,6 +16361,88 @@ uint8_t is_final_check, float double_check_rate, uint8_t flag, float drop_rate) } +void force_trio_clean(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 flag, float self_drop_rate, float contig_drop_rate, uint32_t min_occ) +{ + asg_t* nsg = ug->g; + uint32_t beg, end, n_vtx = nsg->n_seq<<1, k, i, rId, tf_occ, tnf_occ, flag_occ, non_flag_occ, n_reduce = 1; + long long nodeLen, baseLen, max_stop_nodeLen, max_stop_baseLen; + buf_t b; memset(&b, 0, sizeof(buf_t)); + ma_utg_t *u = NULL; + uint8_t* primary_flag = (uint8_t*)calloc(read_g->n_seq, sizeof(uint8_t)); + + while (n_reduce) + { + n_reduce = 0; + n_vtx = nsg->n_seq; + for (beg = 0; beg < n_vtx; ++beg) + { + if(nsg->seq[beg>>1].del || asg_arc_n(nsg, beg) <= 0 || get_real_length(nsg, beg, NULL)<=0) + { + continue; + } + if(get_real_length(nsg, beg^1, NULL) == 1)///check if beg is the tig end + { + get_real_length(nsg, beg^1, &end); + if(get_real_length(nsg, end^1, NULL) == 1) continue; + } + b.b.n = 0; tf_occ = tnf_occ = 0; + get_unitig(nsg, ug, beg, &end, &nodeLen, &baseLen, &max_stop_nodeLen, &max_stop_baseLen, 1, &b); + + for (i = 0; i < b.b.n; i++) + { + u = &((ug)->u.a[b.b.a[i]>>1]); + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + if(R_INF.trio_flag[rId] == AMBIGU) continue; + if(R_INF.trio_flag[rId] == DROP) continue; + if(R_INF.trio_flag[rId] == flag) tf_occ++; + if(R_INF.trio_flag[rId] != flag) tnf_occ++; + } + } + if(tnf_occ <= ((tnf_occ+tf_occ)*contig_drop_rate)) continue; + + for (i = 0; i < b.b.n; i++) + { + flag_occ = non_flag_occ = 0; + u = &((ug)->u.a[b.b.a[i]>>1]); + for (k = 0; k < u->n; k++) + { + rId = u->a[k]>>33; + if(R_INF.trio_flag[rId] == AMBIGU) continue; + if(R_INF.trio_flag[rId] == DROP) continue; + if(R_INF.trio_flag[rId] == flag) flag_occ++; + if(R_INF.trio_flag[rId] != flag) non_flag_occ++; + } + + if(non_flag_occ <= min_occ) continue; + if(non_flag_occ <= ((non_flag_occ+flag_occ)*self_drop_rate)) continue; + if(non_flag_occ <= ((tf_occ+tnf_occ)*contig_drop_rate)) continue; + if(if_primary_unitig(u, read_g, coverage_cut, sources, ruIndex, primary_flag)) + { + continue; + } + + if(u->m != 0) + { + u->circ = u->end = u->len = u->m = u->n = u->start = 0; + free(u->a); + u->a = NULL; + } + asg_seq_del(nsg, b.b.a[i]>>1); + n_reduce++; + } + + } + } + + free(b.b.a); + asg_cleanup(nsg); + free(primary_flag); +} + 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) { @@ -16943,6 +17025,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); + force_trio_clean((*ug), read_g, coverage_cut, sources, reverse_sources, ruIndex, flag, 0.55, 0.01, 5); + ///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); diff --git a/hic.cpp b/hic.cpp index 4284a23..8029243 100644 --- a/hic.cpp +++ b/hic.cpp @@ -6173,7 +6173,8 @@ inline double get_trans_weight_advance(const ha_ug_index* idx, uint64_t x, trans } else { - rate = 0.2; + // rate = 0.2; + rate = 0.4; } @@ -13633,11 +13634,11 @@ void resolve_bubble_chain_by_hic(ha_ug_index *idx, kv_u_trans_t *ta, bubble_type // } // } - for (i = 0; i < edges.a.n; i++) - { - fprintf(stderr, "s-utg%.6lul<------>d-utg%.6ul\n", (edges.a.a[i].ul>>33) + 1, (edges.a.a[i].v>>1) + 1); + // for (i = 0; i < edges.a.n; i++) + // { + // fprintf(stderr, "s-utg%.6lul<------>d-utg%.6ul\n", (edges.a.a[i].ul>>33) + 1, (edges.a.a[i].v>>1) + 1); - } + // } /*******************************for debug************************************/ if(edges.a.n > 0) { @@ -15250,11 +15251,11 @@ void print_kv_weight(kv_u_trans_t *ta) fprintf(stderr, "*********ta->n: %u\n", (uint32_t)ta->n); for (i = 0; i < ta->n; i++) { - fprintf(stderr, "+qn(%u)->tn(%u): %f\n", ta->a[i].qn, ta->a[i].tn, ta->a[i].nw); + fprintf(stderr, "+s-utg%.6ul->d-utg%.6ul: %f\n", ta->a[i].qn+1, ta->a[i].tn+1, ta->a[i].nw); get_u_trans_spec(ta, ta->a[i].tn, ta->a[i].qn, &e, NULL); if(e) { - fprintf(stderr, "-tn(%u)->qn(%u): %f\n", e->qn, e->tn, e->nw); + fprintf(stderr, "-d-utg%.6ul->s-utg%.6ul: %f\n", e->qn+1, e->tn+1, e->nw); } else { @@ -15693,12 +15694,12 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o resolve_tangles_hic(idx, &bub, &sl.hits, &k_trans); measure_distance(idx, idx->ug, &sl.hits, &link, &bub, &(idx->t_ch->k_trans)); - if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_ps_t(&s, asm_opt.output_file_name)) - { - bub.round_id = bub.n_round; - label_unitigs_sm(s->s, NULL, idx->ug); - goto skip_flipping; - } + // if((asm_opt.flag & HA_F_VERBOSE_GFA) && load_ps_t(&s, asm_opt.output_file_name)) + // { + // bub.round_id = bub.n_round; + // label_unitigs_sm(s->s, NULL, idx->ug); + // goto skip_flipping; + // } s = init_ps_t(11, idx->ug->g->n_seq); for (bub.round_id = 0; bub.round_id < bub.n_round; bub.round_id++) { @@ -15726,9 +15727,9 @@ int hic_short_align(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_o label_unitigs(&(hap.g_p), idx->ug); **/ } - write_ps_t(s, asm_opt.output_file_name); + // write_ps_t(s, asm_opt.output_file_name); - skip_flipping: + // skip_flipping: verbose_het_stat(&bub); // print_kv_weight(&k_trans); @@ -15886,7 +15887,7 @@ int hic_short_align_poy(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, // } /*******************************for debug************************************/ } - // write_psg_t(s, asm_opt.output_file_name); + write_psg_t(s, asm_opt.output_file_name); skip_flipping: verbose_het_stat(&bub);