mirror of
https://github.com/chhylp123/hifiasm.git
synced 2026-10-03 03:08:12 +08:00
update polishing
This commit is contained in:
chhylp123
+304
-456
@@ -9045,12 +9045,172 @@ ma_sub_t *coverage_cut, int max_hang, int min_ovlp)
|
||||
return 1;
|
||||
}
|
||||
|
||||
void reduce_ma_utg_t(ma_utg_t* collection, asg_t* read_g, ma_hit_t_alloc* sources,
|
||||
ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp)
|
||||
{
|
||||
asg_arc_t t_f;
|
||||
asg_arc_t* av = NULL;
|
||||
uint32_t m = 0, k, i, nv;
|
||||
for (i = 0; i < collection->n; i++)
|
||||
{
|
||||
if(collection->a[i] == (uint64_t)-1) continue;
|
||||
collection->a[m] = collection->a[i];
|
||||
m++;
|
||||
}
|
||||
collection->n = m;
|
||||
|
||||
|
||||
uint32_t totalLen = 0, v, w, l;
|
||||
for (i = 0; i < collection->n - 1; i++)
|
||||
{
|
||||
v = (uint64_t)(collection->a[i])>>32;
|
||||
w = (uint64_t)(collection->a[i + 1])>>32;
|
||||
|
||||
|
||||
|
||||
|
||||
/*******************************for debug************************************/
|
||||
l = (uint32_t)-1;
|
||||
av = asg_arc_a(read_g, v);
|
||||
nv = asg_arc_n(read_g, v);
|
||||
for (k = 0; k < nv; k++)
|
||||
{
|
||||
if(av[k].del) continue;
|
||||
if(av[k].v == w)
|
||||
{
|
||||
l = asg_arc_len(av[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(k == nv)
|
||||
{
|
||||
for (k = 0; k < edge->a.n; k++)
|
||||
{
|
||||
if(edge->a.a[k].del) continue;
|
||||
if((edge->a.a[k].ul>>32) == v && edge->a.a[k].v == w)
|
||||
{
|
||||
l = asg_arc_len(edge->a.a[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(k == edge->a.n)
|
||||
{
|
||||
if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, v, w, &t_f)==0)
|
||||
{
|
||||
fprintf(stderr, "####ERROR1: v>>1: %u, v&1: %u, w>>1: %u, w&1: %u, r_seq: %u\n",
|
||||
v>>1, v&1, w>>1, w&1, read_g->r_seq);
|
||||
}
|
||||
l = asg_arc_len(t_f);
|
||||
}
|
||||
}
|
||||
if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n");
|
||||
/*******************************for debug************************************/
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
collection->a[i] = v; collection->a[i] = collection->a[i]<<32;
|
||||
collection->a[i] = collection->a[i] | (uint64_t)(l);
|
||||
totalLen += l;
|
||||
}
|
||||
|
||||
if(i < collection->n)
|
||||
{
|
||||
if(collection->circ)
|
||||
{
|
||||
v = (uint64_t)(collection->a[i])>>32;
|
||||
w = (uint64_t)(collection->a[0])>>32;
|
||||
|
||||
|
||||
|
||||
|
||||
/*******************************for debug************************************/
|
||||
l = (uint32_t)-1;
|
||||
av = asg_arc_a(read_g, v);
|
||||
nv = asg_arc_n(read_g, v);
|
||||
for (k = 0; k < nv; k++)
|
||||
{
|
||||
if(av[k].del) continue;
|
||||
if(av[k].v == w)
|
||||
{
|
||||
l = asg_arc_len(av[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(k == nv)
|
||||
{
|
||||
for (k = 0; k < edge->a.n; k++)
|
||||
{
|
||||
if(edge->a.a[k].del) continue;
|
||||
if((edge->a.a[k].ul>>32) == v && edge->a.a[k].v == w)
|
||||
{
|
||||
l = asg_arc_len(edge->a.a[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(k == edge->a.n)
|
||||
{
|
||||
if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, v, w, &t_f)==0)
|
||||
{
|
||||
fprintf(stderr, "####ERROR1: v>>1: %u, v&1: %u, w>>1: %u, w&1: %u, r_seq: %u\n",
|
||||
v>>1, v&1, w>>1, w&1, read_g->r_seq);
|
||||
}
|
||||
l = asg_arc_len(t_f);
|
||||
}
|
||||
}
|
||||
if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n");
|
||||
/*******************************for debug************************************/
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
collection->a[i] = v; collection->a[i] = collection->a[i]<<32;
|
||||
collection->a[i] = collection->a[i] | (uint64_t)(l);
|
||||
|
||||
totalLen += l;
|
||||
}
|
||||
else
|
||||
{
|
||||
v = (uint64_t)(collection->a[i])>>32;
|
||||
l = read_g->seq[v>>1].len;
|
||||
collection->a[i] = v;
|
||||
collection->a[i] = collection->a[i]<<32;
|
||||
collection->a[i] = collection->a[i] | (uint64_t)(l);
|
||||
totalLen += l;
|
||||
}
|
||||
}
|
||||
|
||||
collection->len = totalLen;
|
||||
if(!collection->circ)
|
||||
{
|
||||
collection->start = collection->a[0]>>32;
|
||||
collection->end = (collection->a[collection->n-1]>>32)^1;
|
||||
}
|
||||
}
|
||||
|
||||
uint32_t polish_unitig(ma_utg_t* collection, asg_t* read_g, ma_hit_t_alloc* sources,
|
||||
ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp)
|
||||
{
|
||||
if(collection->m == 0) return 0;
|
||||
if(collection->n < 3) return 0;
|
||||
uint32_t i, k, v, pre, afte, nv, exactLen, inexactLen, tmp_exactLen, tmp_inexactLen, m = 0, skip = 0;
|
||||
uint32_t i, k, v, pre, afte, nv, exactLen, inexactLen, tmp_exactLen, tmp_inexactLen, skip = 0;
|
||||
asg_arc_t* av = NULL;
|
||||
asg_arc_t *pE = NULL, *aE = NULL;
|
||||
asg_arc_t t_f, t_b;
|
||||
@@ -9064,14 +9224,6 @@ ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp)
|
||||
if(pre == (uint32_t)-1) continue;
|
||||
if(afte == (uint32_t)-1) continue;
|
||||
|
||||
/****************************may have bugs********************************/
|
||||
///if(debug_purge_dup == 1 && ((v>>1)==1239235 || (v>>1)==4576917 || (v>>1)==4479766))
|
||||
if(debug_purge_dup == 1 && (v>>1)==4576917)
|
||||
{
|
||||
fprintf(stderr, "v>>1: %u, pre>>1: %u, afte>>1: %u, len: %u\n", v>>1, pre>>1, afte>>1, collection->len);
|
||||
}
|
||||
/****************************may have bugs********************************/
|
||||
|
||||
av = asg_arc_a(read_g, v^1);
|
||||
nv = asg_arc_n(read_g, v^1);
|
||||
for (k = 0; k < nv; k++)
|
||||
@@ -9124,15 +9276,6 @@ ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp)
|
||||
if(k == edge->a.n) fprintf(stderr, "ERROR\n");
|
||||
}
|
||||
|
||||
/****************************may have bugs********************************/
|
||||
///if(debug_purge_dup == 1 && ((v>>1)==1239235 || (v>>1)==4576917 || (v>>1)==4479766))
|
||||
if(debug_purge_dup == 1 && (v>>1)==4576917)
|
||||
{
|
||||
fprintf(stderr, "v>>1: %u, pE->el: %u, aE->el: %u\n", v>>1, pE->el, aE->el);
|
||||
}
|
||||
/****************************may have bugs********************************/
|
||||
|
||||
|
||||
if(pE->el == 1 && aE->el == 1) continue;
|
||||
|
||||
if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, pre,
|
||||
@@ -9147,31 +9290,12 @@ ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp)
|
||||
continue;
|
||||
}
|
||||
|
||||
|
||||
/****************************may have bugs********************************/
|
||||
///if(debug_purge_dup == 1 && ((v>>1)==1239235 || (v>>1)==4576917 || (v>>1)==4479766))
|
||||
if(debug_purge_dup == 1 && (v>>1)==4576917)
|
||||
{
|
||||
fprintf(stderr, "v>>1: %u, t_f.el: %u, t_b.el: %u\n", v>>1, t_f.el, t_b.el);
|
||||
}
|
||||
/****************************may have bugs********************************/
|
||||
|
||||
if(t_f.el == 0 || t_b.el == 0) continue;
|
||||
|
||||
get_overlapLen(v>>1, sources, &exactLen, &inexactLen);
|
||||
if(pE->el == 0)
|
||||
{
|
||||
get_overlapLen(pre>>1, sources, &tmp_exactLen, &tmp_inexactLen);
|
||||
|
||||
|
||||
/****************************may have bugs********************************/
|
||||
///if(debug_purge_dup == 1 && ((v>>1)==1239235 || (v>>1)==4576917 || (v>>1)==4479766))
|
||||
if(debug_purge_dup == 1 && (v>>1)==4576917)
|
||||
{
|
||||
fprintf(stderr, "v>>1: %u, exactLen: %u, inexactLen: %u, tmp_exactLen: %u, tmp_inexactLen: %u\n",
|
||||
v>>1, exactLen, inexactLen, tmp_exactLen, tmp_inexactLen);
|
||||
}
|
||||
/****************************may have bugs********************************/
|
||||
if(inexactLen < tmp_inexactLen) continue;
|
||||
if(inexactLen == tmp_inexactLen && exactLen > tmp_exactLen) continue;
|
||||
}
|
||||
@@ -9183,450 +9307,179 @@ ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp)
|
||||
if(inexactLen == tmp_inexactLen && exactLen > tmp_exactLen) continue;
|
||||
}
|
||||
|
||||
|
||||
/****************************may have bugs********************************/
|
||||
///if(debug_purge_dup == 1 && ((v>>1)==1239235 || (v>>1)==4576917 || (v>>1)==4479766))
|
||||
if(debug_purge_dup == 1 && (v>>1)==4576917)
|
||||
{
|
||||
fprintf(stderr, "v>>1: %u, i: %u\n", v>>1, i);
|
||||
}
|
||||
/****************************may have bugs********************************/
|
||||
|
||||
collection->a[i] = (uint64_t)-1;
|
||||
skip++;
|
||||
}
|
||||
|
||||
if(skip == 0) return 0;
|
||||
|
||||
m = 0;
|
||||
for (i = 0; i < collection->n; i++)
|
||||
{
|
||||
if(collection->a[i] == (uint64_t)-1) continue;
|
||||
collection->a[m] = collection->a[i];
|
||||
m++;
|
||||
}
|
||||
collection->n = m;
|
||||
|
||||
|
||||
uint32_t totalLen = 0, w, l;
|
||||
for (i = 0; i < collection->n - 1; i++)
|
||||
{
|
||||
v = (uint64_t)(collection->a[i])>>32;
|
||||
w = (uint64_t)(collection->a[i + 1])>>32;
|
||||
|
||||
|
||||
|
||||
|
||||
/*******************************for debug************************************/
|
||||
l = (uint32_t)-1;
|
||||
av = asg_arc_a(read_g, v);
|
||||
nv = asg_arc_n(read_g, v);
|
||||
for (k = 0; k < nv; k++)
|
||||
{
|
||||
if(av[k].del) continue;
|
||||
if(av[k].v == w)
|
||||
{
|
||||
l = asg_arc_len(av[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(k == nv)
|
||||
{
|
||||
for (k = 0; k < edge->a.n; k++)
|
||||
{
|
||||
if(edge->a.a[k].del) continue;
|
||||
if((edge->a.a[k].ul>>32) == v && edge->a.a[k].v == w)
|
||||
{
|
||||
l = asg_arc_len(edge->a.a[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(k == edge->a.n)
|
||||
{
|
||||
if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, v, w, &t_f)==0)
|
||||
{
|
||||
fprintf(stderr, "####ERROR1: v>>1: %u, v&1: %u, w>>1: %u, w&1: %u, r_seq: %u\n",
|
||||
v>>1, v&1, w>>1, w&1, read_g->r_seq);
|
||||
}
|
||||
l = asg_arc_len(t_f);
|
||||
}
|
||||
}
|
||||
if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n");
|
||||
/*******************************for debug************************************/
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
collection->a[i] = v; collection->a[i] = collection->a[i]<<32;
|
||||
collection->a[i] = collection->a[i] | (uint64_t)(l);
|
||||
totalLen += l;
|
||||
}
|
||||
|
||||
if(i < collection->n)
|
||||
{
|
||||
if(collection->circ)
|
||||
{
|
||||
v = (uint64_t)(collection->a[i])>>32;
|
||||
w = (uint64_t)(collection->a[0])>>32;
|
||||
|
||||
|
||||
|
||||
|
||||
/*******************************for debug************************************/
|
||||
l = (uint32_t)-1;
|
||||
av = asg_arc_a(read_g, v);
|
||||
nv = asg_arc_n(read_g, v);
|
||||
for (k = 0; k < nv; k++)
|
||||
{
|
||||
if(av[k].del) continue;
|
||||
if(av[k].v == w)
|
||||
{
|
||||
l = asg_arc_len(av[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(k == nv)
|
||||
{
|
||||
for (k = 0; k < edge->a.n; k++)
|
||||
{
|
||||
if(edge->a.a[k].del) continue;
|
||||
if((edge->a.a[k].ul>>32) == v && edge->a.a[k].v == w)
|
||||
{
|
||||
l = asg_arc_len(edge->a.a[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(k == edge->a.n)
|
||||
{
|
||||
if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, v, w, &t_f)==0)
|
||||
{
|
||||
fprintf(stderr, "####ERROR1: v>>1: %u, v&1: %u, w>>1: %u, w&1: %u, r_seq: %u\n",
|
||||
v>>1, v&1, w>>1, w&1, read_g->r_seq);
|
||||
}
|
||||
l = asg_arc_len(t_f);
|
||||
}
|
||||
}
|
||||
if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n");
|
||||
/*******************************for debug************************************/
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
collection->a[i] = v; collection->a[i] = collection->a[i]<<32;
|
||||
collection->a[i] = collection->a[i] | (uint64_t)(l);
|
||||
|
||||
totalLen += l;
|
||||
}
|
||||
else
|
||||
{
|
||||
v = (uint64_t)(collection->a[i])>>32;
|
||||
l = read_g->seq[v>>1].len;
|
||||
collection->a[i] = v;
|
||||
collection->a[i] = collection->a[i]<<32;
|
||||
collection->a[i] = collection->a[i] | (uint64_t)(l);
|
||||
totalLen += l;
|
||||
}
|
||||
}
|
||||
|
||||
collection->len = totalLen;
|
||||
if(!collection->circ)
|
||||
{
|
||||
collection->start = collection->a[0]>>32;
|
||||
collection->end = (collection->a[collection->n-1]>>32)^1;
|
||||
}
|
||||
reduce_ma_utg_t(collection, read_g, sources, coverage_cut, edge, max_hang, min_ovlp);
|
||||
|
||||
return 1;
|
||||
}
|
||||
|
||||
|
||||
int get_consensus_rate(ma_utg_t* collection, uint32_t cur_i, uint32_t next_i,
|
||||
asg_t* read_g, All_reads *RNF, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
|
||||
kvec_asg_arc_t_warp* edge, UC_Read* r_read, UC_Read* q_read, int max_hang, int min_ovlp,
|
||||
int* r_match, int* r_total)
|
||||
{
|
||||
(*r_match) = (*r_total) = 0;
|
||||
if(cur_i < 1) return -1;
|
||||
asg_arc_t *t = NULL, i_t;
|
||||
uint32_t v, w, k, v_beg, v_end, w_beg, w_end;
|
||||
int j;
|
||||
if(collection->a[cur_i] == (uint64_t)-1 || collection->a[next_i] == (uint64_t)-1) return 0;
|
||||
|
||||
v = (uint64_t)(collection->a[cur_i])>>32;
|
||||
w = (uint64_t)(collection->a[next_i])>>32;
|
||||
|
||||
if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, v, w, &i_t) == 0)
|
||||
{
|
||||
for (k = 0; k < edge->a.n; k++)
|
||||
{
|
||||
if(edge->a.a[k].del) continue;
|
||||
if((edge->a.a[k].ul>>32) == v && edge->a.a[k].v == w)
|
||||
{
|
||||
t = &(edge->a.a[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
t = &i_t;
|
||||
}
|
||||
|
||||
if(t == NULL) return -1;
|
||||
v_beg = 0; v_end = asg_arc_len(*t) - 1; r_read->length = 0;
|
||||
///cur_i must >= 1
|
||||
for (j = cur_i-1; j >= 0; j--)
|
||||
{
|
||||
if(collection->a[j] == (uint64_t)-1) continue;
|
||||
|
||||
w = (uint64_t)(collection->a[j])>>32;
|
||||
t = NULL;
|
||||
|
||||
if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, w, v, &i_t) == 0)
|
||||
{
|
||||
for (k = 0; k < edge->a.n; k++)
|
||||
{
|
||||
if(edge->a.a[k].del) continue;
|
||||
if((edge->a.a[k].ul>>32) == w && edge->a.a[k].v == v)
|
||||
{
|
||||
t = &(edge->a.a[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
t = &i_t;
|
||||
}
|
||||
|
||||
if(t == NULL) break;
|
||||
|
||||
w_beg = asg_arc_len(*t);
|
||||
w_end = MIN(((int)(w_beg + v_end)), ((int)(read_g->seq[w>>1].len-1)));
|
||||
// if(t->el == 1)
|
||||
// {
|
||||
// if(r_read->length == 0) recover_UC_sub_Read(r_read, v_beg, v_end - v_beg + 1, v&1, RNF, v>>1);
|
||||
// ///for debug: check if v[v_beg, v_end] == w[w_beg, w_end]
|
||||
// recover_UC_sub_Read(q_read, w_beg, w_end - w_beg +1, w&1, RNF, w>>1);
|
||||
// ///q_read->length<=r_read->length
|
||||
// if(if_exact_match(r_read->seq, r_read->length, q_read->seq, q_read->length,
|
||||
// 0, q_read->length-1, 0, q_read->length-1) == 0)
|
||||
// {
|
||||
// fprintf(stderr, "ERROR: cur_i: %u, j: %d, r_len: %lld, q_len: %lld\n", cur_i, j, r_read->length, q_read->length);
|
||||
// }
|
||||
// }
|
||||
|
||||
//filter overlaps which cannot cover the whole [v_beg, v_end]
|
||||
if(w_end - w_beg != v_end - v_beg) break;
|
||||
|
||||
(*r_total)++;
|
||||
|
||||
if(t->el == 1)
|
||||
{
|
||||
(*r_match)++;
|
||||
}
|
||||
else
|
||||
{
|
||||
if(r_read->length == 0) recover_UC_sub_Read(r_read, v_beg, v_end - v_beg + 1, v&1, RNF, v>>1);
|
||||
recover_UC_sub_Read(q_read, w_beg, w_end - w_beg +1, w&1, RNF, w>>1);
|
||||
///q_read->length<=r_read->length
|
||||
if(if_exact_match(r_read->seq, r_read->length, q_read->seq, q_read->length,
|
||||
0, q_read->length-1, 0, q_read->length-1) == 1)
|
||||
{
|
||||
(*r_match)++;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
uint32_t polish_unitig_advance(ma_utg_t* collection, asg_t* read_g, ma_hit_t_alloc* sources,
|
||||
ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp)
|
||||
return 1;
|
||||
}
|
||||
|
||||
uint32_t polish_unitig_advance(ma_utg_t* collection, asg_t* read_g, All_reads *RNF,
|
||||
ma_hit_t_alloc* sources, ma_sub_t *coverage_cut, kvec_asg_arc_t_warp* edge,
|
||||
UC_Read* r_read, UC_Read* q_read, int max_hang, int min_ovlp)
|
||||
{
|
||||
if(collection->m == 0) return 0;
|
||||
if(collection->n < 2) return 0;
|
||||
uint32_t i, k, v, pre, afte, nv, exactLen, inexactLen, tmp_exactLen, tmp_inexactLen, m = 0, skip = 0;
|
||||
asg_arc_t* av = NULL;
|
||||
asg_arc_t *pE = NULL, *aE = NULL;
|
||||
asg_arc_t t_f, t_b;
|
||||
|
||||
for (i = 1; i < collection->n; i++)
|
||||
{
|
||||
v = (uint64_t)(collection->a[i])>>32;
|
||||
pre = (uint64_t)(collection->a[i-1])>>32;
|
||||
afte = (uint64_t)(collection->a[i+1])>>32;
|
||||
if(v == (uint32_t)-1) continue;
|
||||
if(pre == (uint32_t)-1) continue;
|
||||
if(afte == (uint32_t)-1) continue;
|
||||
|
||||
av = asg_arc_a(read_g, v^1);
|
||||
nv = asg_arc_n(read_g, v^1);
|
||||
for (k = 0; k < nv; k++)
|
||||
{
|
||||
if(av[k].del) continue;
|
||||
if(av[k].v == (pre^1))
|
||||
{
|
||||
pE = &(av[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
if(k == nv)
|
||||
{
|
||||
for (k = 0; k < edge->a.n; k++)
|
||||
{
|
||||
if(edge->a.a[k].del) continue;
|
||||
if((edge->a.a[k].ul>>32) == (v^1) && edge->a.a[k].v == (pre^1))
|
||||
{
|
||||
pE = &(edge->a.a[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(k == edge->a.n) fprintf(stderr, "ERROR\n");
|
||||
}
|
||||
|
||||
av = asg_arc_a(read_g, v);
|
||||
nv = asg_arc_n(read_g, v);
|
||||
for (k = 0; k < nv; k++)
|
||||
{
|
||||
if(av[k].del) continue;
|
||||
if(av[k].v == afte)
|
||||
{
|
||||
aE = &(av[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(k == nv)
|
||||
{
|
||||
for (k = 0; k < edge->a.n; k++)
|
||||
{
|
||||
if(edge->a.a[k].del) continue;
|
||||
if((edge->a.a[k].ul>>32) == v && edge->a.a[k].v == afte)
|
||||
{
|
||||
aE = &(edge->a.a[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
if(k == edge->a.n) fprintf(stderr, "ERROR\n");
|
||||
}
|
||||
|
||||
|
||||
if(pE->el == 1 && aE->el == 1) continue;
|
||||
if(collection->n < 3) return 0;
|
||||
uint32_t i, skip = 0;
|
||||
int match_v, total_v, max_i, match_max, k;
|
||||
double match_rate, match_rate_max;
|
||||
|
||||
if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, pre,
|
||||
afte, &t_f) == 0)
|
||||
///we should be able to handle i = collection->n-1
|
||||
for (i = 1; i < collection->n-1; i++)
|
||||
{
|
||||
///in practice, collection->a[i] and collection->a[i+1] must be available
|
||||
///collection->a[index] might be unavailable only if index < i
|
||||
if(get_consensus_rate(collection, i, i+1, read_g, RNF, sources, coverage_cut,
|
||||
edge, r_read, q_read, max_hang, min_ovlp, &match_v, &total_v) != 1)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
match_rate = (total_v == 0)? 0:((double)(match_v)/(double)(total_v));
|
||||
///most reads support collection[i], so it is right
|
||||
if(match_v >= total_v * 0.5 && total_v > 0 && match_v > 0) continue;
|
||||
|
||||
if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, afte^1,
|
||||
pre^1, &t_b) == 0)
|
||||
max_i = i; match_max = match_v; match_rate_max = match_rate;
|
||||
for (k = i - 1; k >= 0; k--)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
|
||||
|
||||
if(t_f.el == 0 || t_b.el == 0) continue;
|
||||
|
||||
get_overlapLen(v>>1, sources, &exactLen, &inexactLen);
|
||||
if(pE->el == 0)
|
||||
{
|
||||
get_overlapLen(pre>>1, sources, &tmp_exactLen, &tmp_inexactLen);
|
||||
if(inexactLen < tmp_inexactLen) continue;
|
||||
if(inexactLen == tmp_inexactLen && exactLen > tmp_exactLen) continue;
|
||||
}
|
||||
|
||||
if(aE->el == 0)
|
||||
{
|
||||
get_overlapLen(afte>>1, sources, &tmp_exactLen, &tmp_inexactLen);
|
||||
if(inexactLen < tmp_inexactLen) continue;
|
||||
if(inexactLen == tmp_inexactLen && exactLen > tmp_exactLen) continue;
|
||||
}
|
||||
|
||||
|
||||
collection->a[i] = (uint64_t)-1;
|
||||
skip++;
|
||||
}
|
||||
|
||||
if(skip == 0) return 0;
|
||||
|
||||
m = 0;
|
||||
for (i = 0; i < collection->n; i++)
|
||||
{
|
||||
if(collection->a[i] == (uint64_t)-1) continue;
|
||||
collection->a[m] = collection->a[i];
|
||||
m++;
|
||||
}
|
||||
collection->n = m;
|
||||
|
||||
|
||||
uint32_t totalLen = 0, w, l;
|
||||
for (i = 0; i < collection->n - 1; i++)
|
||||
{
|
||||
v = (uint64_t)(collection->a[i])>>32;
|
||||
w = (uint64_t)(collection->a[i + 1])>>32;
|
||||
|
||||
|
||||
|
||||
|
||||
/*******************************for debug************************************/
|
||||
l = (uint32_t)-1;
|
||||
av = asg_arc_a(read_g, v);
|
||||
nv = asg_arc_n(read_g, v);
|
||||
for (k = 0; k < nv; k++)
|
||||
{
|
||||
if(av[k].del) continue;
|
||||
if(av[k].v == w)
|
||||
if(collection->a[k] == (uint64_t)-1) continue;
|
||||
///collection->a[k] might be unavailable, while collection->a[i+1] must be available
|
||||
///return -1 means there is no overlap from k to i+1
|
||||
if(get_consensus_rate(collection, k, i+1, read_g, RNF, sources, coverage_cut,
|
||||
edge, r_read, q_read, max_hang, min_ovlp, &match_v, &total_v) < 0)
|
||||
{
|
||||
l = asg_arc_len(av[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
///no read support k to i+1
|
||||
if(total_v == 0) break;
|
||||
|
||||
if(k == nv)
|
||||
{
|
||||
for (k = 0; k < edge->a.n; k++)
|
||||
match_rate = (total_v == 0)? 0:((double)(match_v)/(double)(total_v));
|
||||
if(match_rate > match_rate_max || (match_rate == match_rate_max && match_v > match_max))
|
||||
{
|
||||
if(edge->a.a[k].del) continue;
|
||||
if((edge->a.a[k].ul>>32) == v && edge->a.a[k].v == w)
|
||||
{
|
||||
l = asg_arc_len(edge->a.a[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(k == edge->a.n)
|
||||
{
|
||||
if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, v, w, &t_f)==0)
|
||||
{
|
||||
fprintf(stderr, "####ERROR1: v>>1: %u, v&1: %u, w>>1: %u, w&1: %u, r_seq: %u\n",
|
||||
v>>1, v&1, w>>1, w&1, read_g->r_seq);
|
||||
}
|
||||
l = asg_arc_len(t_f);
|
||||
max_i = k; match_max = match_v; match_rate_max = match_rate;
|
||||
}
|
||||
}
|
||||
if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n");
|
||||
/*******************************for debug************************************/
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
collection->a[i] = v; collection->a[i] = collection->a[i]<<32;
|
||||
collection->a[i] = collection->a[i] | (uint64_t)(l);
|
||||
totalLen += l;
|
||||
}
|
||||
|
||||
if(i < collection->n)
|
||||
{
|
||||
if(collection->circ)
|
||||
|
||||
///set [max_i+1, i] to be unavailable
|
||||
for (k = max_i+1; k <= (int)i; k++)
|
||||
{
|
||||
v = (uint64_t)(collection->a[i])>>32;
|
||||
w = (uint64_t)(collection->a[0])>>32;
|
||||
|
||||
|
||||
|
||||
|
||||
/*******************************for debug************************************/
|
||||
l = (uint32_t)-1;
|
||||
av = asg_arc_a(read_g, v);
|
||||
nv = asg_arc_n(read_g, v);
|
||||
for (k = 0; k < nv; k++)
|
||||
{
|
||||
if(av[k].del) continue;
|
||||
if(av[k].v == w)
|
||||
{
|
||||
l = asg_arc_len(av[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(k == nv)
|
||||
{
|
||||
for (k = 0; k < edge->a.n; k++)
|
||||
{
|
||||
if(edge->a.a[k].del) continue;
|
||||
if((edge->a.a[k].ul>>32) == v && edge->a.a[k].v == w)
|
||||
{
|
||||
l = asg_arc_len(edge->a.a[k]);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if(k == edge->a.n)
|
||||
{
|
||||
if(get_edge_from_source(sources, coverage_cut, NULL, max_hang, min_ovlp, v, w, &t_f)==0)
|
||||
{
|
||||
fprintf(stderr, "####ERROR1: v>>1: %u, v&1: %u, w>>1: %u, w&1: %u, r_seq: %u\n",
|
||||
v>>1, v&1, w>>1, w&1, read_g->r_seq);
|
||||
}
|
||||
l = asg_arc_len(t_f);
|
||||
}
|
||||
}
|
||||
if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n");
|
||||
/*******************************for debug************************************/
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
collection->a[i] = v; collection->a[i] = collection->a[i]<<32;
|
||||
collection->a[i] = collection->a[i] | (uint64_t)(l);
|
||||
|
||||
totalLen += l;
|
||||
}
|
||||
else
|
||||
{
|
||||
v = (uint64_t)(collection->a[i])>>32;
|
||||
l = read_g->seq[v>>1].len;
|
||||
collection->a[i] = v;
|
||||
collection->a[i] = collection->a[i]<<32;
|
||||
collection->a[i] = collection->a[i] | (uint64_t)(l);
|
||||
totalLen += l;
|
||||
if(collection->a[k] == (uint64_t)-1) continue;
|
||||
collection->a[k] = (uint64_t)-1;
|
||||
skip++;
|
||||
}
|
||||
}
|
||||
|
||||
collection->len = totalLen;
|
||||
if(!collection->circ)
|
||||
{
|
||||
collection->start = collection->a[0]>>32;
|
||||
collection->end = (collection->a[collection->n-1]>>32)^1;
|
||||
}
|
||||
if(skip == 0) return 1;
|
||||
|
||||
reduce_ma_utg_t(collection, read_g, sources, coverage_cut, edge, max_hang, min_ovlp);
|
||||
|
||||
return 1;
|
||||
}
|
||||
// generate unitig sequences
|
||||
|
||||
|
||||
|
||||
// generate unitig sequences
|
||||
@@ -9653,6 +9506,8 @@ ma_hit_t_alloc* sources, kvec_asg_arc_t_warp* edge, int max_hang, int min_ovlp)
|
||||
|
||||
polish_unitig(u, read_g, sources, coverage_cut, edge, max_hang, min_ovlp);
|
||||
|
||||
polish_unitig_advance(u, read_g, RNF, sources, coverage_cut, edge, &g_read, &tmp, max_hang, min_ovlp);
|
||||
|
||||
uint32_t l = 0;
|
||||
u->s = (char*)calloc(1, u->len + 1);
|
||||
memset(u->s, 'N', u->len);
|
||||
@@ -9734,8 +9589,8 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, uint8_t* r_flag)
|
||||
if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue;
|
||||
}
|
||||
if(read_g->seq[tn].del == 1) continue;
|
||||
C_bases += (Get_qe((*h)) - Get_qs((*h)));
|
||||
if(r_flag[tn] != 1) continue;
|
||||
C_bases += (Get_qe((*h)) - Get_qs((*h)));
|
||||
}
|
||||
}
|
||||
|
||||
@@ -23196,13 +23051,7 @@ R_to_U* ruIndex, float chimeric_rate, float drop_ratio, int max_hang, int min_ov
|
||||
adjust_utg_by_primary(&ug, sg, TRIO_THRES, sources, reverse_sources, coverage_cut,
|
||||
bubble_dist, tipsLen, tip_drop_ratio, stops_threshold, ruIndex, chimeric_rate, drop_ratio,
|
||||
max_hang, min_ovlp, &new_rtg_edges);
|
||||
/****************************may have bugs********************************/
|
||||
debug_purge_dup = 1;
|
||||
/****************************may have bugs********************************/
|
||||
ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp);
|
||||
/****************************may have bugs********************************/
|
||||
debug_purge_dup = 0;
|
||||
/****************************may have bugs********************************/
|
||||
|
||||
|
||||
fprintf(stderr, "Writing primary contig GFA to disk... \n");
|
||||
@@ -27354,9 +27203,8 @@ long long bubble_dist, int read_graph, int write)
|
||||
&R_INF, output_file_name);
|
||||
}
|
||||
|
||||
debug_info_of_specfic_read("m64062_190803_042216/177341795/ccs", sources, reverse_sources, -1, "beg");
|
||||
///debug_info_of_specfic_read("m64062_190803_042216/177341795/ccs", sources, reverse_sources, -1, "beg");
|
||||
// debug_info_of_specfic_read("m64062_190807_194840/126682874/ccs", sources, reverse_sources, -1, "beg");
|
||||
exit(1);
|
||||
|
||||
if (!(asm_opt.flag & HA_F_BAN_ASSEMBLY))
|
||||
{
|
||||
|
||||
Reference in New Issue
Block a user