xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-09-15 12:47:57 +08:00
Code Issues Packages Projects Releases Wiki Activity

r257: dropped contained unitigs

Browse Source
This commit is contained in:
Heng Li
2020-06-08 15:43:21 -04:00
parent afa93e4bbb
commit 36bfb02298
4 changed files with 109 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
CommandLines.cpp
View File

@@ -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);

3
CommandLines.h
View File

@@ -3,7 +3,7 @@
#include <pthread.h>
#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

5
Overlaps.h
View File

@@ -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 *

114
overlap2.cpp
View File

@@ -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;
}