From 36bfb02298715a55a8cad149606c721561a03fec Mon Sep 17 00:00:00 2001 From: Heng Li Date: Mon, 8 Jun 2020 15:43:21 -0400 Subject: [PATCH] r257: dropped contained unitigs --- CommandLines.cpp | 2 + CommandLines.h | 3 +- Overlaps.h | 5 ++- overlap2.cpp | 114 +++++++++++++++++++++++++++++++++++++++++------ 4 files changed, 109 insertions(+), 15 deletions(-) diff --git a/CommandLines.cpp b/CommandLines.cpp index b165697..d4881ca 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -19,6 +19,7 @@ static ko_longopt_t long_options[] = { { "max-od-final", ko_no_argument, 306 }, { "ex-list", ko_required_argument, 307 }, { "ex-iter", ko_required_argument, 308 }, + { "keep-contained",ko_no_argument, 309 }, { 0, 0, 0 } }; @@ -404,6 +405,7 @@ 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->flag |= HA_F_KEEP_CONTAINED; 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 d6ff1b5..9b492ad 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -3,7 +3,7 @@ #include -#define HA_VERSION "0.7-dirty-r255" +#define HA_VERSION "0.7-dirty-r257" #define VERBOSE 0 @@ -15,6 +15,7 @@ #define HA_F_SKIP_TRIOBIN 0x20 #define HA_F_PURGE_CONTAIN 0x40 #define HA_F_PURGE_JOIN 0x80 +#define HA_F_KEEP_CONTAINED 0x100 #define HA_MIN_OV_DIFF 0.02 // min sequence divergence in an overlap diff --git a/Overlaps.h b/Overlaps.h index 05cd11b..4543540 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -83,7 +83,6 @@ char* output_file_name); typedef struct { uint32_t s:31, del:1, e; uint8_t c; - uint8_t weak_contain; } ma_sub_t; void ma_hit_sub(int min_dp, ma_hit_t_alloc* sources, long long n_read, uint64_t* readLen, @@ -519,6 +518,10 @@ void asg_con_sort(asg_t *g); void asg_con_index(asg_t *g); asg_t *ma_sg_gen(const ma_hit_t_alloc* sources, long long n_read, const ma_sub_t *coverage_cut, int max_hang, int min_ovlp); int asg_arc_del_trans(asg_t *g, int fuzz); +int asg_drop_contained_utg(asg_t *g); + +#define asg_con_n(g, v) ((g)->contain? (uint32_t)(g)->con_idx[(v)] : 0) +#define asg_con_a(g, v) ((g)->contain? &(g)->contain[(g)->con_idx[(v)]>>32] : 0) /******************* * Other rountines * diff --git a/overlap2.cpp b/overlap2.cpp index c2a2157..ee1a8d8 100644 --- a/overlap2.cpp +++ b/overlap2.cpp @@ -2,6 +2,7 @@ #include "utils.h" #include "CommandLines.h" #include "Overlaps.h" +#include "Process_Read.h" /******************************* * Dropping strong containment * @@ -62,9 +63,8 @@ void ma_hit_contained_advance(ma_hit_t_alloc *sources, long long n_read, ma_sub_ r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); //assert(r != MA_HT_INT && r != MA_HT_SHORT_OVLP); if (r == MA_HT_QCONT) { - if (h->ml || 1) { + if (h->ml || (asm_opt.flag & HA_F_KEEP_CONTAINED) == 0) { h->del = 1; - ++n_strong_contain; delete_single_edge(sources, coverage_cut, Get_tn(*h), Get_qn(*h)); delete_all_edges(sources, coverage_cut, Get_qn(*h)); set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0); @@ -73,14 +73,12 @@ void ma_hit_contained_advance(ma_hit_t_alloc *sources, long long n_read, ma_sub_ // set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0); // sq->del = 1; // set_R_to_U(ruIndex, Get_qn(*h), Get_tn(*h), 0); - } else { - sq->weak_contain = 1; - ++n_weak_contain; } + if (h->ml) ++n_strong_contain; + else ++n_weak_contain; } else if (r == MA_HT_TCONT) { - if (h->ml || 1) { + if (h->ml || (asm_opt.flag & HA_F_KEEP_CONTAINED) == 0) { h->del = 1; - ++n_strong_contain; delete_single_edge(sources, coverage_cut, Get_tn(*h), Get_qn(*h)); delete_all_edges(sources, coverage_cut, Get_tn(*h)); set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0); @@ -89,10 +87,9 @@ void ma_hit_contained_advance(ma_hit_t_alloc *sources, long long n_read, ma_sub_ // set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0); // st->del = 1; // set_R_to_U(ruIndex, Get_tn(*h), Get_qn(*h), 0); - } else { - st->weak_contain = 1; - ++n_weak_contain; } + if (h->ml) ++n_strong_contain; + else ++n_weak_contain; } } } @@ -177,7 +174,6 @@ asg_t *ma_sg_gen(const ma_hit_t_alloc* sources, long long n_read, const ma_sub_t ql = coverage_cut[qn].e - coverage_cut[qn].s; tl = coverage_cut[tn].e - coverage_cut[tn].s; r = ma_hit2arc(h, ql, tl, max_hang, asm_opt.max_hang_rate, min_ovlp, &t); - assert(r >= 0); if (r >= 0) { p = asg_arc_pushp(g); *p = t; @@ -265,9 +261,28 @@ int asg_arc_del_trans(asg_t *g, int fuzz) for (i = 0; i < nv; ++i) { uint32_t w = av[i].v; info[w].mark = g->seq[w>>1].del? 2 : 1; + //if (asg_con_n(g, w>>1) > 0) info[w].mark = 2; info[w].len = asg_arc_len(av[i]); } + // remove contained reads + for (i = 0; i < nv; ++i) { + uint32_t j, nw, w = av[i].v; + uint64_t *aw; + if (info[w].mark != 1) continue; + nw = asg_con_n(g, w>>1); + if (nw == 0) continue; + aw = asg_con_a(g, w>>1); + for (j = 0; j < nw; ++j) { + uint32_t x = (uint32_t)aw[j]; + if (w&1) x ^= 1; + if (info[x].mark == 1 && info[x].len <= info[w].len) + break; + } + if (j < nw) info[w].mark = 2; + //if (nw > 0) fprintf(stderr, "X\t%.*s\t%.*s\n", (int)Get_NAME_LENGTH(R_INF, w>>1), Get_NAME(R_INF, w>>1), (int)Get_NAME_LENGTH(R_INF, (uint32_t)aw[j]>>1), Get_NAME(R_INF, (uint32_t)aw[j]>>1)); + } + // length of node (not overlap length) // av[nv-1] is longest out-dege /** @@ -284,14 +299,17 @@ int asg_arc_del_trans(asg_t *g, int fuzz) for (i = 0; i < nv; ++i) { uint32_t w = av[i].v; uint32_t j, nw = asg_arc_n(g, w); + uint32_t is_con = (asg_con_n(g, w>>1) > 0); asg_arc_t *aw = asg_arc_a(g, w); if (info[w].mark != 1) continue; for (j = 0; j < nw; ++j) { uint32_t x, sum = asg_arc_len(aw[j]) + asg_arc_len(av[i]); if (sum > L) break; x = aw[j].v; - if (info[x].mark == 1 && sum < info[x].len + fuzz && sum + fuzz > info[x].len) - info[x].mark = 2; + if (info[x].mark == 1 && sum < info[x].len + fuzz && sum + fuzz > info[x].len) { + if (!is_con || asg_con_n(g, x>>1) > 0) + info[x].mark = 2; + } } } #if 0 @@ -303,6 +321,7 @@ int asg_arc_del_trans(asg_t *g, int fuzz) if (info[aw[j].v].mark) info[aw[j].v].mark = 2; } #endif + // remove edges for (i = 0; i < nv; ++i) { if (info[av[i].v].mark == 2) av[i].del = 1, ++n_reduced; @@ -314,6 +333,75 @@ int asg_arc_del_trans(asg_t *g, int fuzz) if (n_reduced) { asg_cleanup(g); asg_symm(g); + asg_drop_contained_utg(g); } + fprintf(stderr, "[M::%s] transitively reduced %d arcs\n", __func__, n_reduced); return n_reduced; } + +#define GFA_VT_MERGEABLE 0 +#define GFA_VT_TIP 1 +#define GFA_VT_MULTI_OUT 2 +#define GFA_VT_MULTI_IN 3 + +static inline int32_t gfa_deg(const asg_t *g, uint32_t v, uint32_t *w) +{ + uint32_t i, nv, nv0, k; + const asg_arc_t *av; + if (w) *w = (uint32_t)-1; + if (g->seq[v>>1].del) return 0; + nv0 = k = asg_arc_n(g, v); + av = asg_arc_a(g, v); + for (i = nv = 0; i < nv0; ++i) + if (!av[i].del) + ++nv, k = i; + if (w) *w = nv == 1? av[k].v : (uint32_t)-1; + return nv; +} + +static inline int32_t gfa_vtype(const asg_t *g, uint32_t v, uint32_t *w_) +{ + int32_t nv, nw; + uint32_t w; + nv = gfa_deg(g, v, &w); + if (w_) *w_ = w; + if (nv == 0) return GFA_VT_TIP; + if (nv > 1) return GFA_VT_MULTI_OUT; + nw = gfa_deg(g, w^1, 0); + return nw == 1? GFA_VT_MERGEABLE : GFA_VT_MULTI_IN; +} + +int asg_drop_contained_utg(asg_t *g) +{ + uint32_t n_vtx = g->n_seq * 2, v, cnt = 0; + if (g->contain == 0) return 0; + for (v = 0; v < n_vtx; ++v) { + int32_t vt, is_contained; + uint32_t w; + if (g->seq[v>>1].del) continue; + if (asg_con_n(g, v>>1) == 0) continue; + vt = gfa_vtype(g, v^1, &w); + if (vt == GFA_VT_MERGEABLE) continue; + w = v, is_contained = 1; + while (1) { + if (asg_con_n(g, w>>1) == 0) { + is_contained = 0; + break; + } + vt = gfa_vtype(g, w, &w); + if (vt != GFA_VT_MERGEABLE) break; + } + if (is_contained) { + w = v; + while (1) { + ++cnt; + asg_seq_del(g, w>>1); + vt = gfa_vtype(g, w, &w); + if (vt != GFA_VT_MERGEABLE) break; + } + } + } + if (cnt > 0) asg_cleanup(g); + fprintf(stderr, "[M::%s] drop %d reads in contained unitigs\n", __func__, cnt); + return cnt; +}