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

fix gap alignment extention

Browse Source
This commit is contained in:
chhylp123
2022-06-23 22:11:28 -04:00
parent 3934167770
commit a82e05c9f2
5 changed files with 703 additions and 98 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.h
+1 -1
View File
@@ -4,7 +4,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.16.4-r405"
#define HA_VERSION "0.16.5-r408"
#define VERBOSE 0
gfa_ut.cpp
+614 -77
View File
@@ -213,6 +213,7 @@ typedef struct {
kvec_t(ul_snp_t) snp;
poa_g_t pg;
kvec_t(uint64_t) res_dump;
uint64_t n_correct, n_circle;
}integer_t;
typedef struct {
@@ -221,6 +222,11 @@ typedef struct {
uint64_t n_thread;
}integer_ml_t;
typedef struct {
asg_t *g;
uint64_t n[2], tot;
} cul_g_t;
typedef struct{
ma_ug_t *init_ug;
ma_ug_t *l0_ug;
@@ -233,6 +239,7 @@ typedef struct{
ubuf_t buf;
ul_str_idx_t pstr;
integer_ml_t str_b;
cul_g_t *cg;
// ul_path_srt_t psrt;
}ul_resolve_t;
@@ -2917,7 +2924,10 @@ ma_ug_t *ug, ul_str_idx_t *str_idx, all_ul_t *ul_idx, ul_chain_t *res)
sc += a[i].sc;
if(a[i].tk <= a[i-1].tk) break;///== means there is a circle
if(((uint32_t)a[i].tn_rev_qk) <= ((uint32_t)a[i-1].tn_rev_qk)) break;
if(!dis_check_integer_aln_t(ul_idx, str_idx, ug, &(a[i]), &(a[i-1]), qid, tid, is_rev, 0.08, 2000)) break;
if(str_idx && ul_idx && (!dis_check_integer_aln_t(ul_idx, str_idx,
ug, &(a[i]), &(a[i-1]), qid, tid, is_rev, 0.08, 2000))) {
break;
}
}
// if(tid == 269 || tid == 276 || tid == 277 || tid == 278) fprintf(stderr, "tid->%u, i->%ld, an->%ld\n", tid, i, a_n);
if(i >= a_n) {
@@ -2941,7 +2951,10 @@ ma_ug_t *ug, ul_str_idx_t *str_idx, all_ul_t *ul_idx, ul_chain_t *res)
///qk of lk and li might be equal
if(lk->tk >= li->tk || ((uint32_t)lk->tn_rev_qk) >= ((uint32_t)li->tn_rev_qk)) continue;
// if(is_circle && (!dis_check_integer_aln_t(ul_idx, str_idx, ug, li, lk, qid, tid, is_rev, 0.08))) continue;
if(!dis_check_integer_aln_t(ul_idx, str_idx, ug, li, lk, qid, tid, is_rev, 0.08, 2000)) continue;
if(str_idx && ul_idx && (!dis_check_integer_aln_t(ul_idx, str_idx,
ug, li, lk, qid, tid, is_rev, 0.08, 2000))) {
continue;
}
sc = csc + f[k];
if(sc > max_f) {
max_f = sc; max_k = k;
@@ -3198,27 +3211,57 @@ void print_integer_seq(ma_ug_t *ug, ul_str_t *str, int64_t id, int64_t is_header
fprintf(stderr,"\n");
}
void print_integer_g(poa_g_t *g, ma_ug_t *ug)
#define poa_arc_n(g, v) ((uint32_t)(g)->idx.a[(v)])
#define poa_arc_a(g, v) (&(g)->arc.a[(g)->idx.a[(v)]>>32])
void print_integer_g(poa_g_t *g, ma_ug_t *ug, uint64_t is_gfa)
{
uint64_t i, v, w;
for (i = 0; i < g->seq.n; i++) {
fprintf(stderr, "Node::[M::%s::i->%lu] utg%.6d%c(%c)\n", __func__, i, (int32_t)(g->seq.a[i].nid>>1)+1,
"lc"[ug->u.a[g->seq.a[i].nid>>1].circ], "+-"[g->seq.a[i].nid&1]);
if(!is_gfa) {
for (i = 0; i < g->seq.n; i++) {
fprintf(stderr, "Node::[M::%s::i->%lu] utg%.6d%c(%c)\n", __func__, i, (int32_t)(g->seq.a[i].nid>>1)+1,
"lc"[ug->u.a[g->seq.a[i].nid>>1].circ], "+-"[g->seq.a[i].nid&1]);
}
for (i = 0; i < g->arc.n; i++) {
v = g->arc.a[i].ul>>32; w = g->arc.a[i].v;
if(v&1) continue;
v = g->seq.a[v>>1].nid; w = g->seq.a[w>>1].nid;
fprintf(stderr, "Arch::[M::%s::w->%u] utg%.6d%c(%c)<nid::%lu> -> utg%.6d%c(%c)<nid::%u>\n", __func__, (uint32_t)g->arc.a[i].ul,
(int32_t)(v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1], g->arc.a[i].ul>>33,
(int32_t)(w>>1)+1, "lc"[ug->u.a[w>>1].circ], "+-"[w&1], g->arc.a[i].v>>1);
}
for (i = 0; i < g->seq.n; i++) {
fprintf(stderr, "Srt::[M::%s::i->%lu] utg%.6d%c(%c)<nid::%u>\n", __func__, i,
(int32_t)(g->seq.a[g->srt_b.res.a[i]].nid>>1)+1,
"lc"[ug->u.a[g->seq.a[g->srt_b.res.a[i]].nid>>1].circ], "+-"[g->seq.a[g->srt_b.res.a[i]].nid&1], g->srt_b.res.a[i]);
}
} else {
char name[32];
for (i = 0; i < g->seq.n; i++) {
sprintf(name, "dtg%.6d", (int)i);
fprintf(stderr, "S\t%s\t*\tLN:i:%u\trd:i:utg%.6d%c(%c)\n", name, g->srt_b.res2nid.a[i],
(int32_t)(g->seq.a[i].nid>>1)+1, "lc"[ug->u.a[g->seq.a[i].nid>>1].circ], "+-"[g->seq.a[i].nid&1]);
}
uint32_t u, nu, j; poa_arc_t *au;
for (i = 0; i < g->seq.n; ++i) {
u = i<<1;
au = poa_arc_a(g, u); nu = poa_arc_n(g, u);
for (j = 0; j < nu; j++) {
v = au[j].v;
fprintf(stderr, "L\tdtg%.6d\t%c\tdtg%.6d\t%c\t%dM\n",
(int)(u>>1), "+-"[u&1], (int)(v>>1), "+-"[v&1], asg_arc_len(au[j]));
}
u = (i<<1) + 1;
au = poa_arc_a(g, u); nu = poa_arc_n(g, u);
for (j = 0; j < nu; j++) {
v = au[j].v;
fprintf(stderr, "L\tdtg%.6d\t%c\tdtg%.6d\t%c\t%dM\n",
(int)(u>>1), "+-"[u&1], (int)(v>>1), "+-"[v&1], asg_arc_len(au[j]));
}
}
}
for (i = 0; i < g->arc.n; i++) {
v = g->arc.a[i].ul>>32; w = g->arc.a[i].v;
if(v&1) continue;
v = g->seq.a[v>>1].nid; w = g->seq.a[w>>1].nid;
fprintf(stderr, "Arch::[M::%s::w->%u] utg%.6d%c(%c) -> utg%.6d%c(%c)\n", __func__, (uint32_t)g->arc.a[i].ul,
(int32_t)(v>>1)+1, "lc"[ug->u.a[v>>1].circ], "+-"[v&1],
(int32_t)(w>>1)+1, "lc"[ug->u.a[w>>1].circ], "+-"[w&1]);
}
for (i = 0; i < g->seq.n; i++) {
fprintf(stderr, "Srt::[M::%s::i->%lu] utg%.6d%c(%c)\n", __func__, i,
(int32_t)(g->seq.a[g->srt_b.res.a[i]].nid>>1)+1,
"lc"[ug->u.a[g->seq.a[g->srt_b.res.a[i]].nid>>1].circ], "+-"[g->seq.a[g->srt_b.res.a[i]].nid&1]);
}
}
void print_cns_seq(ma_ug_t *ug, ul_str_t *str, uint64_t *cns_seq, uint64_t cns_occ)
@@ -3581,9 +3624,6 @@ uint64_t *cns, uint64_t cns_occ)
return 1;
}
#define poa_arc_n(g, v) ((uint32_t)(g)->idx.a[(v)])
#define poa_arc_a(g, v) (&(g)->arc.a[(g)->idx.a[(v)]>>32])
void reset_poa_g_t(poa_g_t *g)
{
g->seq.n = g->arc.n = g->idx.n = 0; g->update_arc = g->update_seq = g->e_idx.n = 0;
@@ -3770,6 +3810,7 @@ void append_integer_seq_frag(ma_ug_t *ug, uc_block_t *raw, poa_g_t *g, int64_t g
update_poa_nid_occ(ug, raw, g, g_end, str, (is_rev?(0):(str_occ-1)), is_rev);
if(str_occ < 2) return;
insert_poa_nodes_0(ug, raw, g, (is_rev?(str+1):(str)), str_occ-1, is_rev, str_id, str_off+(is_rev?(1):(0)));
push_poa_arch_0(g, g->seq.n-1, g_end, str_id, str_off + (is_rev?(1):(str_occ-2)),
str_off + (is_rev?(0):(str_occ-1)));
return;
@@ -3847,8 +3888,11 @@ void append_aligned_integer_seq_by_aln_pair(ma_ug_t *ug, uc_block_t *raw, poa_g_
qstr_n = p_str+1-(poa_str_idx(((uint32_t)a[k].tn_rev_qk), str_occ, is_rev));
qoff = poa_str_idx(((uint32_t)a[k].tn_rev_qk), str_occ, is_rev) + str_off;
}
// fprintf(stderr, "+[M::%s::k->%ld] p_g::%ld, c_g::%u, p_str::%ld, c_str::%ld\n",
// __func__, k, p_g, gidx[a[k].tk], p_str, poa_str_idx(((uint32_t)a[k].tn_rev_qk), str_occ, is_rev));
if(str_id == 47072) {
fprintf(stderr, "+[M::%s::k->%ld] p_g::%ld, c_g::%u, p_str::%ld, c_str::%ld, str_off::%ld, qoff::%lu, qstr_n::%lu\n",
__func__, k, p_g, gidx[a[k].tk], p_str, poa_str_idx(((uint32_t)a[k].tn_rev_qk), str_occ, is_rev),
str_off, qoff, qstr_n);
}
append_integer_seq_frag(ug, raw, g, p_g, gidx[a[k].tk], qstr_a, qstr_n, is_rev, str_id, qoff);
@@ -3865,7 +3909,7 @@ void append_aligned_integer_seq_by_aln_pair(ma_ug_t *ug, uc_block_t *raw, poa_g_
}
void topo_srt_gen(poa_g_t *g)
uint32_t topo_srt_gen(poa_g_t *g, uint32_t debug_qid, ma_ug_t *ug)
{
uint32_t k, v, w, a_n; poa_arc_t *a;
kv_resize(uint32_t, g->srt_b.ind, g->seq.n);
@@ -3890,6 +3934,15 @@ void topo_srt_gen(poa_g_t *g)
if(g->srt_b.ind.a[w] == 0) kv_push(uint32_t, g->srt_b.stack, w);
}
}
// if(!(g->srt_b.res.n == g->seq.n)) {
// fprintf(stderr, "[M::%s::] debug_qid::%u, g->seq.n::%u, g->srt_b.res.n::%u\n",
// __func__, debug_qid, (uint32_t)g->seq.n, (uint32_t)g->srt_b.res.n);
// print_integer_g(g, ug, 1);
// for (k = 0; k < g->seq.n; k++) {
// if(g->srt_b.ind.a[k] > 0) fprintf(stderr, "circle->nid::%u\n", k);
// }
// }
if(g->srt_b.res.n != g->seq.n) return 0;///there is a circle
assert(g->srt_b.res.n == g->seq.n);
for (k = 0; k < g->seq.n; k++) {
g->srt_b.res2nid.a[g->srt_b.res.a[k]] = k;
@@ -3897,6 +3950,7 @@ void topo_srt_gen(poa_g_t *g)
// g->srt_b.aln.a[k] <<= 32; g->srt_b.aln.a[k] += k;
}
// radix_sort_srt64(g->srt_b.aln.a, g->srt_b.aln.a + g->srt_b.aln.n);
return 1;
}
void init_poa_dp(ma_ug_t *ug, poa_dp_t *dp, poa_g_t *g, uint64_t g_occ, uint64_t *str, uint64_t str_occ, uint64_t is_rev, uc_block_t *raw, integer_t *buf)
@@ -3936,9 +3990,9 @@ void init_poa_dp(ma_ug_t *ug, poa_dp_t *dp, poa_g_t *g, uint64_t g_occ, uint64_t
}
}
void update_poa_dp(poa_g_t *g)
uint32_t update_poa_dp(poa_g_t *g, uint32_t debug_qid, ma_ug_t *debug_ug)
{
uint32_t is_srt = 0/**, is_up_aln = 0**/, k;
uint32_t is_srt = 0/**, is_up_aln = 0**/, k, is_circle = 0;
if(g->seq.n > g->update_seq || g->arc.n > g->update_arc) {
if(g->seq.n > g->update_seq) {
kv_resize(uint32_t, g->srt_b.res, g->seq.n);
@@ -3960,15 +4014,17 @@ void update_poa_dp(poa_g_t *g)
}
if(k < g->arc.n) is_srt = 1;
}
// fprintf(stderr, "[M::%s::] is_srt::%u\n", __func__, is_srt);
clean_poa_g_t(g);
if(is_srt) {
topo_srt_gen(g);
is_circle = 1 - topo_srt_gen(g, debug_qid, debug_ug);
}
// else if(is_up_aln) {
// radix_sort_srt64(g->srt_b.aln.a, g->srt_b.aln.a + g->srt_b.aln.n);
// }
}
return is_circle;
}
/**
@@ -4101,7 +4157,7 @@ ul_chain_t *idx, int64_t idx_n, int64_t qid, integer_t *buf)
}
**/
int64_t suffix_gorder_check(poa_g_t *g, integer_t *buf, uint64_t gk_0, uint64_t gk_1, int64_t update_vis)
int64_t suffix_gorder_check(poa_g_t *g, integer_t *buf, uint64_t gk_0, uint64_t gk_1, int64_t update_vis, uint64_t set_flag)
{
uint32_t *g_idx = g->srt_b.res.a, *n2gidx = g->srt_b.res2nid.a, a_n, v, init_n, k;
poa_arc_t *a;
@@ -4115,17 +4171,17 @@ int64_t suffix_gorder_check(poa_g_t *g, integer_t *buf, uint64_t gk_0, uint64_t
kv_push(uint32_t, buf->vis, v);
while (buf->vis.n > init_n) {
v = buf->vis.a[--buf->vis.n];
buf->vis.a[n2gidx[v>>1]] = gk_0;
buf->vis.a[n2gidx[v>>1]] = set_flag;
a_n = poa_arc_n(g, v); a = poa_arc_a(g, v);
for (k = 0; k < a_n; k++) {
if(buf->vis.a[n2gidx[a[k].v>>1]] == gk_0) continue;
if(buf->vis.a[n2gidx[a[k].v>>1]] == set_flag) continue;
kv_push(uint32_t, buf->vis, a[k].v);
}
}
assert(buf->vis.n == init_n);
}
if(buf->vis.a[gk_1] == gk_0) return 1;
if(buf->vis.a[gk_1] == set_flag) return 1;
return 0;
}
@@ -4134,17 +4190,18 @@ int64_t integer_g_chain(poa_g_t *g, ma_ug_t *ug, integer_aln_t *a, int64_t a_n,
res->v = res->s = res->e = (uint32_t)-1; res->sc = (uint64_t)-1;
res->q_sidx = res->q_eidx = res->t_sidx = res->t_eidx = (uint32_t)-1;
if(a_n <= 0) return 0;
int64_t i, k, *p, *f, tf, ti, csc, sc, max_f, max_k, vis_i, pas; integer_aln_t *li, *lk;
int64_t i, k, *p, *f, tf, ti, csc, sc, max_f, max_k, vis_i, pas; integer_aln_t *li, *lk;
buf->vis.n = 0; vis_i = -1;
for (i = 1; i < a_n; ++i) {//already sorted by qk
if(((uint32_t)a[i].tn_rev_qk) <= ((uint32_t)a[i-1].tn_rev_qk)) break; ///== means there is a circle
if(a[i].tk == a[i-1].tk) break;
if(a[i].tk < a[i-1].tk) {
pas = suffix_gorder_check(g, buf, i, i-1, vis_i==i?0:1); vis_i = i;
pas = suffix_gorder_check(g, buf, a[i].tk, a[i-1].tk, vis_i==i?0:1, i); vis_i = i;
if(pas) break;
}
}
fprintf(stderr, "[M::%s::] i::%ld, a_n::%ld\n", __func__, i, a_n);
if(i >= a_n) {
res->s = 0; res->e = a_n;
return 1;
@@ -4154,17 +4211,22 @@ int64_t integer_g_chain(poa_g_t *g, ma_ug_t *ug, integer_aln_t *a, int64_t a_n,
kv_resize(int64_t, buf->p, (uint64_t)a_n); p = buf->p.a;
kv_resize(int64_t, buf->f, (uint64_t)a_n); f = buf->f.a;
tf = ti = -1;
tf = ti = -1; buf->vis.n = 0; vis_i = -1;
for (i = 0; i < a_n; ++i) {
li = &(a[i]); csc = li->sc;
max_f = csc; max_k = -1;
for (k = i-1; k >= 0; --k) {
lk = &(a[k]);
///qk of lk and li might be equal
if(((uint32_t)lk->tn_rev_qk) >= ((uint32_t)li->tn_rev_qk)) continue;
if(lk->tk == li->tk) continue;
if(lk->tk > li->tk) {
pas = suffix_gorder_check(g, buf, i, k, vis_i==i?0:1);
pas = suffix_gorder_check(g, buf, a[i].tk, a[k].tk, vis_i==i?0:1, i);
// if(i == 17 && a[i].tk == 6 && k == 16 && a[k].tk == 60) {
// fprintf(stderr, "******i::%ld, a[i].tk::%u, k::%ld, a[k].tk::%u, vis_i::%ld, pas::%ld\n",
// i, a[i].tk, k, a[k].tk, vis_i, pas);
// }
vis_i = i; if(pas) continue;
}
sc = csc + f[k];
@@ -4176,6 +4238,8 @@ int64_t integer_g_chain(poa_g_t *g, ma_ug_t *ug, integer_aln_t *a, int64_t a_n,
if(tf < max_f) {
tf = max_f; ti = i;
}
fprintf(stderr, "[M::%s::i->%ld] qk::%u, tk::%u, max_k::%ld, max_f::%ld, ti::%ld\n",
__func__, i, (uint32_t)li->tn_rev_qk, li->tk, max_k, max_f, ti);
}
if(ti < 0) return 0;
@@ -4189,8 +4253,9 @@ int64_t integer_g_chain(poa_g_t *g, ma_ug_t *ug, integer_aln_t *a, int64_t a_n,
}
void poa_chain_0(poa_g_t *g, ma_ug_t *ug, uc_block_t *raw, ul_str_t *str, int64_t s, int64_t e, int64_t is_rev, integer_t *buf, int64_t str_id)
void poa_chain_0(poa_g_t *g, ma_ug_t *ug, uc_block_t *raw, ul_str_t *str, int64_t s, int64_t e, int64_t is_rev, integer_t *buf, int64_t str_id, int64_t debug_qid, uint32_t *is_circle)
{
(*is_circle) = 0;
if(e <= s) return;
uint32_t *g_idx = g->srt_b.res.a; uc_block_t *z; integer_aln_t *b; ul_chain_t rr; int64_t i, k, n;
uint64_t *pat = (is_rev?(str->a + str->cn - e):(str->a + s)), pp; int64_t pat_n = e - s;
@@ -4236,7 +4301,7 @@ void poa_chain_0(poa_g_t *g, ma_ug_t *ug, uc_block_t *raw, ul_str_t *str, int64_
i = integer_g_chain(g, ug, buf->b.a, buf->b.n, buf, &rr); assert(i);
if(!i) return; buf->b.n = rr.e; assert(buf->b.n);
append_aligned_integer_seq_by_aln_pair(ug, raw, g, g->seq.n, pat, pat_n, is_rev, buf->b.a, buf->b.n, str_id, (is_rev?(str->cn - e):(s)));
update_poa_dp(g);
(*is_circle) = update_poa_dp(g, debug_qid, ug);
}
void gen_cns_by_poa(poa_g_t *g)
@@ -4304,20 +4369,23 @@ void gen_cns_by_poa(poa_g_t *g)
// }
}
void poa_cns_chain(poa_g_t *g, all_ul_t *ul_idx, ma_ug_t *ug, ul_str_t *str, ul_chain_t *idx, int64_t idx_n, int64_t qid, integer_t *buf)
void poa_cns_chain(poa_g_t *g, all_ul_t *ul_idx, ma_ug_t *ug, ul_str_t *str, ul_chain_t *idx, int64_t idx_n, int64_t qid, integer_t *buf, uint32_t *is_circle)
{
int64_t k, tid, is_rev;
(*is_circle) = 0;
int64_t k, tid, is_rev;
reset_poa_g_t(g);
append_unmatch_integer_seq(g, ug, ul_idx->a[qid].bb.a, &(str[qid]), 0, str[qid].cn, 0, qid);
clean_poa_g_t(g); topo_srt_gen(g);
clean_poa_g_t(g); (*is_circle) = 1 - topo_srt_gen(g, qid, ug);
if((*is_circle)) return;
for (k = 0; k < idx_n; k++) {
tid = idx[k].v>>1; is_rev = idx[k].v&1;
// fprintf(stderr, "\n[M::%s::] k::%ld, tid::%ld, is_rev::%ld\n", __func__, k, tid, is_rev);
// print_integer_seq(ug, str, tid, 1);
poa_chain_0(g, ug, ul_idx->a[tid].bb.a, &(str[tid]), idx[k].t_sidx, idx[k].t_eidx, is_rev, buf, tid);
// print_integer_g(g, ug);
fprintf(stderr, "\n[M::%s::] k::%ld, tid::%ld, is_rev::%ld\n", __func__, k, tid, is_rev);
print_integer_seq(ug, str, tid, 1);
poa_chain_0(g, ug, ul_idx->a[tid].bb.a, &(str[tid]), idx[k].t_sidx, idx[k].t_eidx, is_rev, buf, tid, qid, is_circle);
if((*is_circle)) return;
print_integer_g(g, ug, 1);
}
gen_cns_by_poa(g);
@@ -4325,8 +4393,10 @@ void poa_cns_chain(poa_g_t *g, all_ul_t *ul_idx, ma_ug_t *ug, ul_str_t *str, ul_
uint64_t cal_forward_dis(asg_t *g, uc_block_t *a, uint32_t s, uint32_t e)
{
uint32_t i, li = (uint32_t)-1, v, w, nv, z; int64_t l; asg_arc_t *av;
for(i = s, l = 0, v = w = (uint32_t)-1; i != (uint32_t)-1 && i <= e; i = a[i].aidx) {
uint32_t i, li, v, w, nv, z; int64_t l; asg_arc_t *av;
fprintf(stderr, "\n[M::%s::] s::%u, e::%u\n", __func__, s, e);
///TODO: a[i].aidx might be < i; if s == e, then i <= e might be wrong, cannot pass the assert(li == e);
for(i = s, l = 0, v = w = (uint32_t)-1, li = s; i != (uint32_t)-1 && i <= /**!=**/ e; i = a[i].aidx) {
w = (((uint32_t)(a[i].hid))<<1)|((uint32_t)(a[i].rev));
if(v != (uint32_t)-1) {
av = asg_arc_a(g, v); nv = asg_arc_n(g, v);
@@ -4341,6 +4411,8 @@ uint64_t cal_forward_dis(asg_t *g, uc_block_t *a, uint32_t s, uint32_t e)
}
}
v = w; li = i;
fprintf(stderr, "[M::%s::] i::%u, li::%u, a[i].aidx::%u, a[i].qs::%u, a[i].qe::%u\n", __func__, i, li, a[i].aidx,
a[i].qs, a[i].qe);
}
assert(li == e);
if(l < 0) l = 0;
@@ -4356,12 +4428,14 @@ uint64_t cal_integer_match_dis(ma_ug_t *ug, uc_block_t *a, int64_t k_0, int64_t
} else {
i = k_1; k = k_0;
}
fprintf(stderr, "k_0::%ld, k_1::%ld\n", k_0, k_1);
uint32_t li, lk, pk, bi = i; int64_t l;
for (li = i, l = 0; i != (uint32_t)-1 && i >= k; i = a[i].pidx) {
li = i; if(a[i].pidx != (uint32_t)-1 && i > k) l += a[i].pdis;
}
fprintf(stderr, "+bi::%u, li::%u, i::%u, l::%ld\n", bi, li, i, l);
if(is_rev) l = cal_forward_dis(ug->g, a, li, bi);
fprintf(stderr, "++bi::%u, li::%u, i::%u, l::%ld\n", bi, li, i, l);
if(li == k) {///direct path
(*is_g_connect) = 1;
return l;
@@ -4369,11 +4443,13 @@ uint64_t cal_integer_match_dis(ma_ug_t *ug, uc_block_t *a, int64_t k_0, int64_t
i = li;
assert(i > k); pk = k;
for (lk = k; k != (uint32_t)-1 && k <= i; k = a[k].aidx) lk = k;
fprintf(stderr, "-pk::%u, lk::%u, k::%u\n", pk, lk, k);
if(!is_rev) {
for (k = lk; k != (uint32_t)-1 && k != pk; k = a[k].pidx) l += a[k].pdis;
} else {
l += cal_forward_dis(ug->g, a, pk, lk);
}
fprintf(stderr, "--pk::%u, lk::%u, k::%u, l::%ld\n", pk, lk, k, l);
if(l < 0) l = 0;
k = lk;
assert(i > k);
@@ -4459,10 +4535,21 @@ uint64_t cal_integer_most_dis(uint64_t *a, uint64_t a_n, double cluster_rate)
return ((uint32_t)a[max_i])>>1;
}
uint32_t poa_g_arc_w(poa_g_t *pg, uint32_t v, uint32_t w)
{
poa_arc_t *av; uint32_t an, k;
an = poa_arc_n(pg, v); av = poa_arc_a(pg, v);
for (k = 0; k < an; k++) {
if(av[k].v == w) return (uint32_t)av[k].ul;
}
return 0;
}
void update_raw_integer_seq(poa_g_t *pg, ma_ug_t *ug, uint32_t *cns_seq, uint32_t cns_occ, all_ul_t *ul_idx, ul_str_t *str, uint32_t qid, integer_t *buf, ul_chain_t *idx_a, uint64_t idx_n)
{
if(cns_occ <= 0) return;
radix_sort_emap_t_srt(pg->e_idx.a, pg->e_idx.a + pg->e_idx.n);
///no need to check edges if there is just one node -> no dege
if(cns_occ > 1) radix_sort_emap_t_srt(pg->e_idx.a, pg->e_idx.a + pg->e_idx.n);
kv_resize(uint64_t, buf->u, cns_occ); buf->u.n = cns_occ - 1; memset(buf->u.a, 0, sizeof(*(buf->u.a)*buf->u.n));
uint64_t *arc_idx = buf->u.a, arc_idx_n = buf->u.n, x; uint64_t k, l, i, n = pg->e_idx.n, fe = cns_occ - 1;
for (k = 1, l = 0; k <= n && fe > 0; k++) {
@@ -4494,12 +4581,17 @@ void update_raw_integer_seq(poa_g_t *pg, ma_ug_t *ug, uint32_t *cns_seq, uint32_
t = qid; t <<= 32; t |= ((uint64_t)(0xffffffff));
kv_push(uint64_t, buf->res_dump, t);
t = pg->seq.a[cns_seq[0]>>1].nid; t |= ((uint64_t)(0xffffffff00000000));
t = pg->seq.a[cns_seq[0]>>1].nid;
// t = pg->seq.a[cns_seq[0]>>1].nid; t |= ((uint64_t)(0xffffffff00000000));
t |= ((uint64_t)(ug->g->seq[pg->seq.a[cns_seq[0]>>1].nid>>1].len))<<32;
kv_push(uint64_t, buf->res_dump, t);
buf->n_correct++;
}
for (k = 0; k < arc_idx_n; k++) {
// csn_v = pg->seq.a[cns_seq[k]>>1].nid; cns_w = pg->seq.a[cns_seq[k+1]>>1].nid;
e_s = arc_idx[k]>>32; e_e = (uint32_t)arc_idx[k]; assert(e_e > e_s);
e_s = arc_idx[k]>>32; e_e = (uint32_t)arc_idx[k];
fprintf(stderr, "[M::%s::] k::%lu, arc_idx_n::%lu, e_s::%u, e_e::%u\n", __func__, k, arc_idx_n, e_s, e_e);
assert(poa_g_arc_w(pg, cns_seq[k], cns_seq[k+1]) == (e_e - e_s)); assert(e_e > e_s);
buf->o.n = 0; kv_resize(uint64_t, buf->o, e_e - e_s); con_occ = 0;
for (i = e_s; i < e_e; i++) {
g_arc = &(pg->e_idx.a[i]);
@@ -4507,9 +4599,11 @@ void update_raw_integer_seq(poa_g_t *pg, ma_ug_t *ug, uint32_t *cns_seq, uint32_
is_rev = ((g_arc->ule>>32) > ((uint32_t)g_arc->ule)?1:0);
assert((pg->seq.a[g_arc->pge>>32].nid^(is_rev?1:0)) == ((uint32_t)str[g_arc->ulid].a[g_arc->ule>>32]));
assert((pg->seq.a[(uint32_t)g_arc->pge].nid^(is_rev?1:0)) == ((uint32_t)str[g_arc->ulid].a[(uint32_t)g_arc->ule]));
fprintf(stderr, "+++i::%lu, target_ulid::%u, str_sidx::%u, str_eidx::%u, str_cn::%u, ul_idx->a[g_arc->ulid].bb.n::%u\n",
i - e_s, g_arc->ulid, (uint32_t)(g_arc->ule>>32), (uint32_t)g_arc->ule, str[g_arc->ulid].cn, (uint32_t)ul_idx->a[g_arc->ulid].bb.n);
dd = cal_integer_match_dis(ug, ul_idx->a[g_arc->ulid].bb.a, str[g_arc->ulid].a[g_arc->ule>>32]>>32,
str[g_arc->ulid].a[(uint32_t)g_arc->ule]>>32, is_rev, &is_g_connect);
fprintf(stderr, "---i::%lu, dd::%lu, is_g_connect::%u\n", i - e_s, dd, is_g_connect);
dd <<= 1; if(is_rev) dd += 1;
if(is_g_connect) {
con_occ++; buf->o.a[buf->o.n] = dd;
@@ -4532,9 +4626,12 @@ void update_raw_integer_seq(poa_g_t *pg, ma_ug_t *ug, uint32_t *cns_seq, uint32_
void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, uint32_t is_hom)
{
// if(qid != 281) return;
uint64_t k, z, m_het, m_het_occ, ref_occ, b_n, m; uint32_t vk, vz; integer_aln_t *p; ul_chain_t sc;
ul_str_idx_t *str_idx = &(uidx->pstr); ma_ug_t *ug = uidx->l1_ug;
// if(qid != 3165) return;
// if(qid != 17165) return;
// if(qid != 24100) return;///circle
if(qid != 27512) return;
uint64_t k, z, m_het, m_het_occ, ref_occ, b_n, m; uint32_t vk, vz, is_circle = 0; integer_aln_t *p; ul_chain_t sc;
ul_str_idx_t *str_idx = &(uidx->pstr); ma_ug_t *ug = uidx->l1_ug;
uint64_t *hid_a, hid_n; uc_block_t *xi;
ul_str_t *str = &(str_idx->str.a[qid]);
if(str->cn < 2) return; ///directly filter out too short UL
@@ -4552,10 +4649,7 @@ void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, uint32_
}
// if((!is_hom) && (m_het < 2) && (m_het > 0)) return;///if all matched unitigs are hom, is ok
if(m_het == 0 || m_het_occ == 0) is_hom = 1;
// fprintf(stderr, "\n");
// print_integer_seq(ug, str_idx->str.a, qid, 1);
print_ul_alignment(ug, &UL_INF, 27512, "inner-0");
for (k = 0, buf->b.n = 0; k < str->cn; k++) {
vk = (uint32_t)str->a[k];
hid_a = str_idx->occ.a + str_idx->idx.a[vk>>1];
@@ -4577,7 +4671,7 @@ void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, uint32_
if(p->sc > m) p->sc = m;
}
}
print_ul_alignment(ug, &UL_INF, 27512, "inner-1");
radix_sort_integer_aln_t_srt(buf->b.a, buf->b.a + buf->b.n);
b_n = buf->b.n; buf->sc.n = 0;
for (k = 1, z = 0; k < b_n; k++) {
@@ -4594,13 +4688,13 @@ void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, uint32_
z = k;
}
}
print_ul_alignment(ug, &UL_INF, 27512, "inner-2");
uint64_t *o, o_n, cns_het, cns_het_occ, ref_cns_occ, corrected = 0;
o_n = integer_chain_dp(uidx->bub, buf, str_idx->str.a, buf->b.a, buf->sc.a, buf->sc.n, qid, is_hom, 2, &corrected);
assert(o_n <= str->cn);
if(o_n <= 0) return;
if(corrected) return;
print_ul_alignment(ug, &UL_INF, 27512, "inner-3");
for (k = cns_het = cns_het_occ = ref_cns_occ = 0, o = buf->o.a; k < o_n; k++) {
if((!is_hom) && (!IF_HOM((((uint32_t)str->a[o[k]])>>1), (*uidx->bub)))) {
@@ -4616,11 +4710,12 @@ void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, uint32_
} else {
if(ref_cns_occ <= (ref_occ*0.25)) return;
}
// print_cns_seq(ug, str, o, o_n);
print_ul_alignment(ug, &UL_INF, 27512, "inner-4");
fprintf(stderr, "\n");
print_integer_seq(ug, str_idx->str.a, qid, 1);
// print_aln_seq(ug, uidx->idx, qid, 1);
print_cns_seq(ug, str, o, o_n);
print_ul_alignment(ug, &UL_INF, 27512, "inner-5");
for (k = m = 0; k < buf->sc.n; k++) {
// fprintf(stderr, "[M::%s::k->%lu] m::%lu\n", __func__, k, m);
@@ -4639,9 +4734,15 @@ void integer_candidate(ul_resolve_t *uidx, integer_t *buf, uint32_t qid, uint32_
}
buf->sc.n = m;
if(m <= 0) return;
print_ul_alignment(ug, &UL_INF, 27512, "inner-6");
// if(m != str->cn) print_integer_ovlps(uidx->l1_ug, str_idx->str.a, buf->b.a, buf->b.n, buf->sc.a, buf->sc.n, qid, m);
poa_cns_chain(&(buf->pg), uidx->idx, ug, str_idx->str.a, buf->sc.a, buf->sc.n, qid, buf);
// print_res_seq(&(buf->pg), ug, buf->pg.srt_b.res.a, buf->pg.srt_b.res.n);
poa_cns_chain(&(buf->pg), uidx->idx, ug, str_idx->str.a, buf->sc.a, buf->sc.n, qid, buf, &is_circle);
if(is_circle) {
buf->n_circle++;
return;
}
print_ul_alignment(ug, &UL_INF, 27512, "inner-7");
print_res_seq(&(buf->pg), ug, buf->pg.srt_b.res.a, buf->pg.srt_b.res.n);
// radix_sort_ul_chain_t_srt(buf->sc.a, buf->sc.a + buf->sc.n);
// integer_phase(str_idx->str.a, buf, buf->sc.a, buf->sc.n, buf->b.a, qid);
@@ -4665,6 +4766,357 @@ static void worker_integer_correction(void *data, long i, int tid) // callback f
integer_candidate(uidx, buf, i, (asm_opt.purge_level_primary == 0?1:0));
}
static void worker_integer_postprecess(void *data, long i, int tid) // callback for kt_for()
{
ul_resolve_t *uidx = (ul_resolve_t *)data;
integer_ml_t *sl = &(uidx->str_b);
integer_t *buf = &(sl->buf[tid]);
// uc_block_t *uls; uint64_t uls_n; uint32_t v;
// uint64_t *srt_a, srt_n, is_circle = ((uidx->psrt.idx.a[i]&((uint64_t)(0x100000000)))?1:0);
// srt_a = uidx->psrt.srt.a + (uint32_t)uidx->psrt.idx.a[i]; srt_n = uidx->psrt.idx.a[i]>>33;
// if(srt_n == 0) return;
// integer_candidate(uidx, srt_a, srt_n, is_circle, buf);
integer_candidate(uidx, buf, i, (asm_opt.purge_level_primary == 0?1:0));
}
void print_primary_ul_chain(ul_str_t *p_str, ul_vec_t *ul, ma_ug_t *ug)
{
uc_block_t *z; uint64_t k, dd = 0;
if(p_str) {
for (k = 0; k < p_str->cn; k++) {
z = &(ul->bb.a[p_str->a[k]>>32]);
fprintf(stderr, "[M::%s::pk->%lu] utg%.6d%c(%c::len->%u), qs::%u, qe::%u, ts::%u, te::%u\n", __func__, k,
z->hid+1, "lc"[ug->u.a[z->hid].circ], "+-"[z->rev], ug->u.a[z->hid].len,
z->qs, z->qe, z->ts, z->te);
dd++;
}
} else {
for (k = 0; k < ul->bb.n; k++) {
z = &(ul->bb.a[k]);
if(!z->pchain) continue;
fprintf(stderr, "[M::%s::rk->%lu] utg%.6d%c(%c::len->%u), qs::%u, qe::%u, ts::%u, te::%u\n", __func__, k,
z->hid+1, "lc"[ug->u.a[z->hid].circ], "+-"[z->rev], ug->u.a[z->hid].len,
z->qs, z->qe, z->ts, z->te);
dd++;
}
}
if(dd) fprintf(stderr, "*******************************************\n");
}
void push_integer_seq_exact(ul_vec_t *res, ma_ug_t *ug, uint64_t *seq, uint64_t seq_n, uint64_t *off, uint32_t tid)
{
if(seq_n <= 0) return;
uint64_t k, v, ql, ul; uc_block_t *x;
// fprintf(stderr, "[M::%s::] old_len::%u, new_len::%u\n", __func__, res->rlen, (uint32_t)off[seq_n-1]);
res->rlen = (uint32_t)off[seq_n-1]; res->bb.n = 0; kv_resize(uc_block_t, res->bb, seq_n);
for (k = 0; k < seq_n; k++) {
v = (uint32_t)seq[k];
kv_pushp(uc_block_t, res->bb, &x);
x->hid = v>>1; x->rev = !!(v&1); x->pchain = 1; x->el = 1; x->base = 0;
x->qs = off[k]>>32; x->qe = (uint32_t)off[k];
ql = x->qe - x->qs; ul = ug->g->seq[x->hid].len;
if((ul == ql) || (k > 0 && k + 1 < seq_n)) {
x->ts = 0; x->te = ul;
} else {
// if(!(k == 0 || k + 1 == seq_n)){
// fprintf(stderr, "[M::%s::] k::%lu, ql::%lu, ul::%lu\n", __func__, k, ql, ul);
// }
// assert(k == 0 || k + 1 == seq_n);
if(k == 0) {
if(x->rev) {
x->ts = 0; x->te = MIN(ql, ul);
} else {
x->ts = ul - MIN(ql, ul); x->te = ul;
}
}
if(k + 1 == seq_n) {
if(!x->rev) {
x->ts = 0; x->te = MIN(ql, ul);
} else {
x->ts = ul - MIN(ql, ul); x->te = ul;
}
}
}
x->pidx = x->aidx = x->pdis = (uint32_t)-1;
if(k > 0 && (seq[k]&((uint64_t)(0x8000000000000000)))) {
x->pidx = k - 1; x->pdis = (seq[k]<<1)>>33;
res->bb.a[x->pidx].aidx = k;
}
}
assert(res->rlen == x->qe);
// res->rlen = ((uint32_t)-1) - tid;
// print_primary_ul_chain(NULL, res, ug);
}
uint32_t double_check_gconnect(asg_t *g, uint32_t v, uint32_t w, uint32_t v2w_d)
{
asg_arc_t *av = asg_arc_a(g, v); uint32_t nv = asg_arc_n(g, v), k;
for (k = 0; k < nv; k++) {
if(av[k].del) continue;
if(av[k].v == w) break;
}
if(k < nv) return 1;
return 0;
}
void update_integer_seq(ul_resolve_t *uidx, integer_t *buf, uint32_t id, uint64_t *seq, uint64_t seq_n, uint32_t tid)
{
// if((((uint32_t)-1) - uidx->idx->a[id].rlen) < asm_opt.thread_num) {
// fprintf(stderr, "id->%u, c_tid->%u, l_tid->%u\n", id, tid, (((uint32_t)-1) - uidx->idx->a[id].rlen));
// exit(0);
// }
// assert(uidx->idx->a[id].rlen!=(uint32_t)-1);
// fprintf(stderr, "dd->%u\n", uidx->idx->a[id].dd);
// if(id != 1487) return;
if(seq_n <= 0) return;
buf->n_correct++;
ul_str_t *str = &(uidx->pstr.str.a[id]); uc_block_t *z; ul_chain_t sc, msc;
uint64_t k, i, pp; integer_aln_t *b; ma_ug_t *ug = uidx->l1_ug;
assert((seq[0]>>32) == (ug->g->seq[((uint32_t)seq[0])>>1].len));
// fprintf(stderr, "\n[M::%s::id->%u] seq_n::%lu\n", __func__, id, seq_n);
// print_primary_ul_chain(str, &(uidx->idx->a[id]), ug);
if(seq_n == 1) {
pp = ug->g->seq[((uint32_t)seq[0])>>1].len;
push_integer_seq_exact(&(uidx->idx->a[id]), ug, seq, seq_n, &pp, tid);
return;
}
buf->u.n = 0;
for (k = 0; k < str->cn; k++) {///old seq
pp = ((uint32_t)str->a[k])>>1; pp <<= 33; pp += ((uint64_t)(0x100000000));
pp += (k<<1); pp += ((uint32_t)str->a[k])&1;
kv_push(uint64_t, buf->u, pp);
// fprintf(stderr, "ok::%lu, ov::%u\n", k, ((uint32_t)str->a[k]));
}
for (k = 0; k < seq_n; k++) {///new seq
pp = ((uint32_t)seq[k])>>1; pp <<= 33;
pp += (k<<1); pp += ((uint32_t)seq[k])&1;
kv_push(uint64_t, buf->u, pp);
// fprintf(stderr, "nk::%lu, nv::%u\n", k, ((uint32_t)seq[k]));
// fprintf(stderr, ">>>>>>[M::%s::c_k->%lu] utg%.6d%c(%c::len->%u)\n", __func__, k,
// (((uint32_t)seq[k])>>1)+1, "lc"[ug->u.a[(((uint32_t)seq[k])>>1)].circ],
// "+-"[((uint32_t)seq[k])&1], ug->u.a[(((uint32_t)seq[k])>>1)].len);
}
radix_sort_srt64(buf->u.a, buf->u.a + buf->u.n);
for (k = 0, buf->b.n = 0; k < buf->u.n; k++) {
if(buf->u.a[k]&((uint64_t)(0x100000000))) continue;///skip nodes in the old seq
for (i = k+1; (i < buf->u.n) && ((buf->u.a[k]>>33) == (buf->u.a[i]>>33)); i++) {
if((buf->u.a[i]&((uint64_t)(0x100000000))) == 0) continue;///skip nodes in the new seq
if((buf->u.a[k]&1) != (buf->u.a[i]&1)) continue;///ignore reverse alignment;
///a[k] is the new seq (q); a[i] is old seq (t)
kv_pushp(integer_aln_t, buf->b, &b);
b->vq = ((buf->u.a[k]>>33)<<1) + (buf->u.a[k]&1);
b->tn_rev_qk = ((uint32_t)buf->u.a[k])>>1;
b->tk = ((uint32_t)buf->u.a[i])>>1;
z = &(uidx->idx->a[id].bb.a[str->a[b->tk]>>32]);
assert(((z->hid<<1)+z->rev)==(((buf->u.a[i]>>33)<<1) + (buf->u.a[i]&1)));
if((buf->u.a[k]&1) != (buf->u.a[i]&1)) {
b->tk = str->cn - b->tk - 1;
b->tn_rev_qk |= ((uint64_t)(0x100000000));
}
b->sc = ug_occ_w(z->ts, z->te, &(ug->u.a[z->hid]));
assert(b->sc > 0);
}
}
msc.v = msc.s = msc.e = (uint32_t)-1; msc.sc = 0;
if(buf->b.n > 0) {
radix_sort_integer_aln_t_srt(buf->b.a, buf->b.a + buf->b.n); ///sorted by rev|qk
for (k = 1, i = 0; k <= buf->b.n; k++) {
// fprintf(stderr, "[M::%s::z->%lu] rev::%u, qk::%u, tk::%u\n",
// __func__, k-1, (uint32_t)(!!(buf->b.a[k-1].tn_rev_qk>>32)), (uint32_t)buf->b.a[k-1].tn_rev_qk, buf->b.a[k-1].tk);
if(k == buf->b.n || (buf->b.a[i].tn_rev_qk>>32) != (buf->b.a[k].tn_rev_qk>>32)) {
if(integer_chain(id, buf->b.a + i, k - i, i, buf, ug, NULL, NULL, &sc) && sc.v != (uint32_t)-1) {
if(msc.s == (uint32_t)-1 || msc.sc < sc.sc) msc = sc;
}
i = k;
}
}
}
// fprintf(stderr, "[M::%s::id->%u] seq_n::%lu, align_n::%u, buf->b.n::%u\n",
// __func__, id, seq_n, msc.e - msc.s, (uint32_t)buf->b.n);
uint64_t l, v, pd, uls, ule, p_ls, p_le, is_rev, qk, tk, bc; integer_aln_t *x; uc_block_t *tb;
buf->u.n = 0; kv_resize(uint64_t, buf->u, seq_n);
// assert((seq[0]>>32) == (ug->g->seq[((uint32_t)seq[0])>>1].len));
for (k = l = 0, p_ls = p_le = (uint64_t)-1; k < seq_n; k++) {
v = (uint32_t)seq[k]; pd = (seq[k]<<1)>>33;
ule = l + pd; uls = ((ule >= ug->g->seq[v>>1].len)?(ule - ug->g->seq[v>>1].len):(0));
bc = 0;
if(k > 0){
bc = (!!(seq[k]&((uint64_t)(0x8000000000000000))));///if connected in the graph
if(bc) bc = double_check_gconnect(ug->g, ((uint32_t)seq[k])^1, ((uint32_t)seq[k-1])^1, pd);
// fprintf(stderr, "bc->%lu\n", bc);
}
if(p_ls != (uint64_t)-1) {
///uls should uls>=p_ls && uls<=p_le
if(uls < p_ls) uls = p_ls;
if(bc && uls > p_le) uls = p_le;
///ule should ule > p_le
if(ule <= p_le) ule = p_le + 1;
}
buf->u.a[k] = uls; buf->u.a[k] <<= 32; buf->u.a[k] |= ule;
p_ls = uls; p_le = ule;
l = ule;
// fprintf(stderr, "[init_k->%lu] uls::%lu, ule::%lu, pd::%lu\n", k, uls, ule, pd);
}
buf->u.n = seq_n;
int64_t beg_nl, end_nl; ma_utg_t *u;
v = (uint32_t)seq[0]; u = &(ug->u.a[v>>1]); beg_nl = u->len;
v = (uint32_t)seq[seq_n-1]; u = &(ug->u.a[v>>1]); end_nl = u->len;
// fprintf(stderr, "+[M::%s] beg_nl::%lu, end_nl::%lu\n", __func__, beg_nl, end_nl);
///q -> new seq; t -> old seq
if(msc.s != (uint32_t)-1 && msc.e > msc.s) {
int64_t qoff, toff, beg_cut = 0, end_cut = 0; uint32_t r_end;
is_rev = ((buf->b.a[msc.s].tn_rev_qk>>32)&1);
///beg
x = &(buf->b.a[msc.s]);
qk = (uint32_t)(x->tn_rev_qk); tk = ((is_rev == 0)? (x->tk):(str->cn - x->tk - 1));
tb = &(uidx->idx->a[id].bb.a[str->a[tk]>>32]);
if((!!tb->rev) == (!!is_rev)) {
if(tb->te == ug->u.a[tb->hid].len) r_end = 1;
else r_end = 0;
} else {
if(tb->ts == 0) r_end = 1;
else r_end = 0;
}
if(r_end) {///start from right end
qoff = (uint32_t)buf->u.a[qk]; toff = (is_rev?(uidx->idx->a[id].rlen-tb->qs):(tb->qe));
} else {///start from left end
qoff = buf->u.a[qk]>>32; toff = (is_rev?(uidx->idx->a[id].rlen-tb->qe):(tb->qs));
}
if(qk == 0) toff = tb->te - tb->ts;
if(toff < qoff) beg_cut = qoff - toff;
// fprintf(stderr, "[M::%s::beg::is_rev->%lu] qoff::%ld, toff::%ld, r_end::%u\n", __func__, is_rev, qoff, toff, r_end);
///end
x = &(buf->b.a[msc.e-1]);
qk = (uint32_t)(x->tn_rev_qk); tk = ((is_rev == 0)? (x->tk):(str->cn - x->tk - 1));
tb = &(uidx->idx->a[id].bb.a[str->a[tk]>>32]);
if((!!tb->rev) == (!!is_rev)) {
if(tb->ts == 0) r_end = 0;
else r_end = 1;
} else {
if(tb->te == ug->u.a[tb->hid].len) r_end = 0;
else r_end = 1;
}
if(!r_end) {
qoff = l - (buf->u.a[qk]>>32); toff = (is_rev?(tb->qe):(uidx->idx->a[id].rlen-tb->qs));
} else {
qoff = l - (uint32_t)(buf->u.a[qk]); toff = (is_rev?(tb->qs):(uidx->idx->a[id].rlen-tb->qe));
}
if(qk+1==seq_n) toff = tb->te - tb->ts;
if(toff < qoff) end_cut = qoff - toff;
// fprintf(stderr, "[M::%s::end::is_rev->%lu] qoff::%ld, toff::%ld, r_end::%u\n", __func__, is_rev, qoff, toff, r_end);
if(beg_cut > 0) {
v = (uint32_t)seq[0]; u = &(ug->u.a[v>>1]);
if(beg_cut < u->len) beg_cut = u->len - beg_cut;
else beg_cut = 0;
if(v&1) {
beg_nl = (int64_t)(Get_READ_LENGTH(R_INF, (u->a[0]>>33)));
} else {
beg_nl = (int64_t)(Get_READ_LENGTH(R_INF, (u->a[u->n-1]>>33)));
}
if(beg_cut > beg_nl) beg_nl = beg_cut;
if(beg_nl == u->len) beg_cut = 0;
else beg_cut = 1;
}
if(end_cut > 0) {
v = (uint32_t)seq[seq_n-1]; u = &(ug->u.a[v>>1]);
if(end_cut < u->len) end_cut = u->len - end_cut;
else end_cut = 0;
if(v&1) {
end_nl = Get_READ_LENGTH(R_INF, (u->a[u->n-1]>>33));
} else {
end_nl = Get_READ_LENGTH(R_INF, (u->a[0]>>33));
}
if(end_cut > end_nl) end_nl = end_cut;
if(end_nl == u->len) end_cut = 0;
else end_cut = 1;
}
// fprintf(stderr, "++[M::%s] beg_nl::%lu, end_nl::%lu\n", __func__, beg_nl, end_nl);
if(beg_cut || end_cut) {
k = 0; uls = 0; ule = beg_nl;
buf->u.a[k] = uls; buf->u.a[k] <<= 32; buf->u.a[k] |= ule;
l = ule; p_ls = uls; p_le = ule;
for (k += 1; k + 1 < seq_n; k++) {
v = (uint32_t)seq[k]; pd = (seq[k]<<1)>>33;
ule = l + pd; uls = ((ule >= ug->g->seq[v>>1].len)?(ule - ug->g->seq[v>>1].len):(0));
bc = 0;
if(k > 0){
bc = (!!(seq[k]&((uint64_t)(0x8000000000000000))));///if connected in the graph
if(bc) bc = double_check_gconnect(ug->g, ((uint32_t)seq[k])^1, ((uint32_t)seq[k-1])^1, pd);
// fprintf(stderr, "bc->%lu\n", bc);
}
if(p_ls != (uint64_t)-1) {
///uls should uls>=p_ls && uls<=p_le
if(uls < p_ls) uls = p_ls;
if(bc && uls > p_le) uls = p_le;
///ule should ule > p_le
if(ule <= p_le) ule = p_le + 1;
}
buf->u.a[k] = uls; buf->u.a[k] <<= 32; buf->u.a[k] |= ule;
p_ls = uls; p_le = ule;
l = ule;
}
assert(k < seq_n);
v = (uint32_t)seq[k]; pd = (seq[k]<<1)>>33;
ule = l + pd; uls = ((ule >= ug->g->seq[v>>1].len)?(ule - ug->g->seq[v>>1].len):(0));
ule = uls + end_nl;
bc = 0;
if(k > 0){
bc = (!!(seq[k]&((uint64_t)(0x8000000000000000))));///if connected in the graph
if(bc) bc = double_check_gconnect(ug->g, ((uint32_t)seq[k])^1, ((uint32_t)seq[k-1])^1, pd);
// fprintf(stderr, "bc->%lu\n", bc);
}
if(p_ls != (uint64_t)-1) {
///uls should uls>=p_ls && uls<=p_le
if(uls < p_ls) uls = p_ls;
if(bc && uls > p_le) uls = p_le;
///ule should ule > p_le
if(ule <= p_le) ule = p_le + 1;
}
if(ule - uls < (uint64_t)end_nl) ule = uls + end_nl;
buf->u.a[k] = uls; buf->u.a[k] <<= 32; buf->u.a[k] |= ule;
p_ls = uls; p_le = ule;
l = ule;
}
}
push_integer_seq_exact(&(uidx->idx->a[id]), ug, seq, seq_n, buf->u.a, tid);
}
static void worker_integer_update(void *data, long i, int tid) // callback for kt_for()
{
ul_resolve_t *uidx = (ul_resolve_t *)data;
integer_ml_t *sl = &(uidx->str_b);
integer_t *buf = &(sl->buf[i]);///normally should be buf = &(sl->buf[tid])
uint64_t k, l;
for (k = 0, l = (uint64_t)-1; k <= buf->res_dump.n; k++) {
if(k == buf->res_dump.n || (buf->res_dump.a[k]&((uint64_t)(0xffffffff))) == ((uint64_t)(0xffffffff))) {
if(l != (uint64_t)-1 && k - l > 1) {
update_integer_seq(uidx, buf, buf->res_dump.a[l]>>32, buf->res_dump.a + l + 1, k - l - 1, tid);
}
l = k;
}
}
}
void integer_correction(ul_resolve_t *uidx)
{
@@ -4672,6 +5124,89 @@ void integer_correction(ul_resolve_t *uidx)
init_integer_ml_t(&sl, uidx, asm_opt.thread_num);
}
uint64_t clean_ul_re_correct_buf(ul_resolve_t *uidx, uint64_t clean_dump, uint64_t *tot_circle)
{
uint64_t k, occ, n_circle;
for (k = occ = n_circle = 0; k < uidx->str_b.n_thread; k++) {
occ += uidx->str_b.buf[k].n_correct;
n_circle += uidx->str_b.buf[k].n_circle;
uidx->str_b.buf[k].n_correct = 0;
uidx->str_b.buf[k].n_circle = 0;
if(clean_dump) uidx->str_b.buf[k].res_dump.n = 0;
uidx->str_b.buf[k].q.n = 0;
uidx->str_b.buf[k].t.n = 0;
uidx->str_b.buf[k].b.n = 0;
uidx->str_b.buf[k].f.n = 0;
uidx->str_b.buf[k].p.n = 0;
uidx->str_b.buf[k].o.n = 0;
uidx->str_b.buf[k].u.n = 0;
uidx->str_b.buf[k].vis.n = 0;
uidx->str_b.buf[k].sc.n = 0;
uidx->str_b.buf[k].snp.n = 0;
}
(*tot_circle) = n_circle;
return occ;
}
void rebuid_idx(ul_resolve_t *uidx)
{
free(uidx->idx->ridx.idx.a); free(uidx->idx->ridx.occ.a);
memset(&(uidx->idx->ridx), 0, sizeof((uidx->idx->ridx)));
filter_ul_ug(uidx->l1_ug);
gen_ul_vec_rid_t(uidx->idx, NULL, uidx->l1_ug);
update_ug_arch_ul_mul(uidx->l1_ug);
free(uidx->pstr.idx.a); free(uidx->pstr.occ.a); free(uidx->pstr.str.a);
memset(&(uidx->pstr), 0, sizeof(uidx->pstr));
init_ul_str_idx_t(uidx);
}
void ul_re_correct(ul_resolve_t *uidx, uint64_t n_r)
{
uint64_t k, occ, n_circle; uidx->str_b.n_thread = 1;
for (k = 0; k < n_r; k++) {
kt_for(uidx->str_b.n_thread, worker_integer_correction, uidx, uidx->idx->n);
occ = clean_ul_re_correct_buf(uidx, 0, &n_circle);
fprintf(stderr, "+[M::%s::round->%lu] # corrected UL reads::%lu, # circle UL reads::%lu\n",
__func__, k, occ, n_circle);
kt_for(uidx->str_b.n_thread, worker_integer_update, uidx, uidx->str_b.n_thread);
occ = clean_ul_re_correct_buf(uidx, 1, &n_circle);
fprintf(stderr, "-[M::%s::round->%lu] # corrected UL reads::%lu, # circle UL reads::%lu\n",
__func__, k, occ, n_circle);
rebuid_idx(uidx);
}
// kt_for(uidx->str_b.n_thread, worker_integer_postprecess, uidx, uidx->idx->n);
}
uint64_t str_occ_w(ul_str_t *str, ul_vec_t *raw, ma_ug_t *ug)
{
uint64_t occ, k; uc_block_t *z;
for (k = occ = 0; k < str->cn; k++) {
z = &(raw->bb.a[str->a[k]>>32]);
assert(((z->hid<<1)+z->rev)==((uint32_t)str->a[k]));
occ += ug_occ_w(z->ts, z->te, &(ug->u.a[z->hid]));
}
return occ;
}
void gen_cul_g_t(ul_resolve_t *uidx)
{
ma_ug_t *ug = uidx->l1_ug; uint64_t k; ul_str_idx_t *str_idx = &(uidx->pstr);
CALLOC(uidx->cg, 1);
uidx->cg->n[0] = uidx->idx->n; uidx->cg->n[1] = ug->g->n_seq;
uidx->cg->tot = uidx->cg->n[0] + uidx->cg->n[1];
uidx->cg->g = asg_init();
for (k = 0; k < uidx->cg->tot; k++) {
if(k < uidx->cg->n[0]) {
asg_seq_set(uidx->cg->g, k, str_occ_w(&(str_idx->str.a[k]), &(uidx->idx->a[k]), ug),
((str_idx->str.a[k].cn>1)?0:1));
} else {
asg_seq_set(uidx->cg->g, k, ug->u.a[k-uidx->cg->n[0]].n, 0);
}
}
}
void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg)
{
@@ -4683,15 +5218,17 @@ void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg)
hic_clean(sg);
ma_ug_t *init_ug = ul_realignment(uopt, sg);
filter_sg_by_ug(sg, init_ug, uopt);
print_ul_alignment(init_ug, &UL_INF, 47072, "after-0");
bub = gen_bubble_chain(sg, init_ug, uopt, &r_het);
print_ul_alignment(init_ug, &UL_INF, 47072, "after-1");
ul_resolve_t *uidx = init_ul_resolve_t(sg, init_ug, bub, &UL_INF, r_het);
kt_for(asm_opt.thread_num, worker_integer_correction, uidx, uidx->idx->n);
print_ul_alignment(init_ug, &UL_INF, 47072, "after-2");
ul_re_correct(uidx, 3);
print_ul_alignment(init_ug, &UL_INF, 47072, "after-3");
print_debug_ul("UL.debug", init_ug, sg, uopt->coverage_cut, uopt->sources, uopt->ruIndex, bub, &UL_INF);
// print_debug_ul("UL.debug", init_ug, sg, uopt->coverage_cut, uopt->sources, uopt->ruIndex, bub, &UL_INF);
resolve_dip_bub_chains(uidx);
// resolve_dip_bub_chains(uidx);
// free(r_het); destory_bubbles(bub); free(bub);
}
gfa_ut.h
-1
View File
@@ -15,5 +15,4 @@ void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float o
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);
void ul_realignment_gfa(ug_opt_t *uopt, asg_t *sg);
#endif
inter.cpp
+84 -19
View File
@@ -4750,8 +4750,9 @@ int64_t adjust_utg_chain_qoffset(uint64_t *r_srt, int64_t *r_pos, int64_t *q_pos
}
assert(k < 3);
rdis = r_pos[(uint32_t)r_srt[k+1]] - r_pos[(uint32_t)r_srt[k]];
qdis = q_pos[(uint32_t)r_srt[k+1]] - q_pos[(uint32_t)r_srt[k]];
rdis = r_pos[(uint32_t)r_srt[k+1]] - r_pos[(uint32_t)r_srt[k]];
if(qdis < 0) return -1;
return q_pos[(uint32_t)r_srt[k]] + get_offset_adjust(r_off-r_pos[(uint32_t)r_srt[k]], rdis, qdis);
}
@@ -4763,7 +4764,8 @@ void update_uovlp_chain_qse(ma_ug_t *ug, int64_t sidx, int64_t eidx, mg_lchain_t
assert(sidx>=0 && eidx<a_n);
int64_t r_pos[4], q_pos[4]; uint64_t r_srt[4];
int64_t i, rs, re;//qs or qe might be -1, while rs and re should >= 0
int64_t i, rs, re, prs, pre, ars, are;//qs or qe might be -1, while rs and re should >= 0
int64_t pqs, pqe, aqs, aqe, fail_s, fail_e;
if(sidx >= 0) {
get_u_offset(ug, &(a[sidx]), &r_pos[0], &r_pos[1], &q_pos[0], &q_pos[1]);
} else {
@@ -4775,6 +4777,11 @@ void update_uovlp_chain_qse(ma_ug_t *ug, int64_t sidx, int64_t eidx, mg_lchain_t
} else {
get_u_offset(ug, &(a[a_n-1]), &r_pos[2], &r_pos[3], &q_pos[2], &q_pos[3]);
}
prs = r_pos[0]; pre = r_pos[1]; ars = r_pos[2]; are = r_pos[3];
pqs = q_pos[0]; pqe = q_pos[1]; aqs = q_pos[2]; aqe = q_pos[3];
// fprintf(stderr, "\n[M::%s::] sidx->%ld<r[0]::%ld, r[1]::%ld, q[0]::%ld, q[1]::%ld>, eidx->%ld<r[2]::%ld, r[3]::%ld, q[2]::%ld, q[3]::%ld>\n", __func__,
// sidx, r_pos[0], r_pos[1], q_pos[0], q_pos[1],
// eidx, r_pos[2], r_pos[3], q_pos[2], q_pos[3]);
// assert((left_q[0] >= 0 && left_q[1] >= 0) || (right_q[0] >= 0 && right_q[1] >= 0)); ///assert(re >= rs);
///for ug chains, sidx >= 0 && eidx < a_n
assert(q_pos[0] >= 0 && q_pos[1] >= 0 && q_pos[2] >= 0 && q_pos[3] >= 0);
@@ -4789,10 +4796,26 @@ void update_uovlp_chain_qse(ma_ug_t *ug, int64_t sidx, int64_t eidx, mg_lchain_t
// assert(r_pos[(uint32_t)r_srt[3]] == (r_srt[3]>>32));
for (i = sidx+1; i < eidx; i++) {
get_u_offset(ug, &(a[i]), &rs, &re, NULL, NULL);
a[i].qs = adjust_utg_chain_qoffset(r_srt, r_pos, q_pos, rs);
a[i].qe = adjust_utg_chain_qoffset(r_srt, r_pos, q_pos, re);
assert(rs >= prs && rs <= ars && re >= pre && re <= are); assert(rs <= re);
a[i].qs = adjust_utg_chain_qoffset(r_srt, r_pos, q_pos, rs); fail_s = 1;
a[i].qe = adjust_utg_chain_qoffset(r_srt, r_pos, q_pos, re); fail_e = 1;
if(a[i].qs >= 0 && a[i].qs >= pqs && a[i].qs <= aqs) fail_s = 0;
if(a[i].qe >= 0 && a[i].qe >= pqe && a[i].qe <= aqe) fail_e = 0;
if(a[i].qs > a[i].qe) fail_s = fail_e = 1;
if(fail_s || fail_e) {
a[i].qe = pqe + get_offset_adjust(re-pre, are-pre, aqe-pqe);///first priority
if(aqs < a[i].qe) {
a[i].qs = pqs + get_offset_adjust(rs-prs, ars-prs, aqs-pqs);
} else {
a[i].qs = pqs + get_offset_adjust(rs-prs, re-prs, a[i].qe-pqs);
}
}
assert(a[i].qs >= pqs && a[i].qs <= aqs && a[i].qe >= pqe && a[i].qe <= aqe); assert(a[i].qs <= a[i].qe);
// fprintf(stderr, "[M::%s::i->%ld] rs::%ld, re::%ld, a[i].qs::%d, a[i].qe::%d\n",
// __func__, i, rs, re, a[i].qs, a[i].qe);
if(i + 1 < eidx) {
q_pos[0] = a[i].qs; q_pos[1] = a[i].qe; r_pos[0] = rs; r_pos[1] = re;
@@ -4805,7 +4828,11 @@ void update_uovlp_chain_qse(ma_ug_t *ug, int64_t sidx, int64_t eidx, mg_lchain_t
// assert(r_pos[(uint32_t)r_srt[1]] == (r_srt[1]>>32));
// assert(r_pos[(uint32_t)r_srt[2]] == (r_srt[2]>>32));
// assert(r_pos[(uint32_t)r_srt[3]] == (r_srt[3]>>32));
// fprintf(stderr, "[Srt::i->%ld] <r[0]::%ld, r[1]::%ld, q[0]::%ld, q[1]::%ld>, <r[2]::%ld, r[3]::%ld, q[2]::%ld, q[3]::%ld>\n", i,
// r_pos[0], r_pos[1], q_pos[0], q_pos[1],
// r_pos[2], r_pos[3], q_pos[2], q_pos[3]);
}
prs = rs; pre = re; pqs = a[i].qs; pqe = a[i].qe;
}
@@ -4862,20 +4889,18 @@ void debug_update_uovlp_chain_qse(ma_ug_t *ug, mg_lchain_t *a, int64_t a_n, int6
void fill_unaligned_alignments(ma_ug_t *ug, mg_lchain_t *a, int64_t a_n, int64_t offset, int64_t ulid)
{
// fprintf(stderr, "a_n->%ld, offset->%ld\n", a_n, offset);
// fprintf(stderr, "\n[M::%s::] a_n->%ld, offset->%ld\n", __func__, a_n, offset);
if(a_n == 0) return;
int64_t k, l;
for (k = 0, l = ug->g->seq[a[0].v>>1].len; k < a_n; k++) {
l -= ug->g->seq[a[k].v>>1].len;
// fprintf(stderr, "k->%ld, l->%ld [M::utg%.6u%c::%c]\n", k, l,
// (a[k].v>>1)+1, "lc"[ug->u.a[a[k].v>>1].circ], "+-"[a[k].v&1]);
if(a[k].off < 0) a[k].qs = a[k].qe = -1;
a[k].off = l; a[k].hash_pre = (uint32_t)-1;
if(k > 0) a[k].hash_pre = offset + k - 1;
// fprintf(stderr, "k->%ld, l->%ld, [M::utg%.6u%c::%c::len->%u], qs->%u, qe->%u, rs->%u, re->%u\n", k, l,
// (a[k].v>>1)+1, "lc"[ug->u.a[a[k].v>>1].circ], "+-"[a[k].v&1], ug->u.a[a[k].v>>1].len,
// a[k].qs, a[k].qe, a[k].rs, a[k].re);
l += ug->g->seq[a[k].v>>1].len + a[k].dist_pre;
}
// debug_update_uovlp_chain_qse(ug, a, a_n, ulid);
@@ -5231,12 +5256,16 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call
int64_t /**rid = s->id+i,**/ winLen = MIN((((double)THRESHOLD_MAX_SIZE)/s->opt->diff_ec_ul), WINDOW);
// uint64_t align = 0;
int fully_cov, abnormal;
if(UL_INF.a[s->id+i].rlen != s->len[i]) {
fprintf(stderr, "[M::%s] rid:%ld, s->len:%lu, UL_INF->rlen:%u\n", __func__, s->id+i, s->len[i], UL_INF.a[s->id+i].rlen);
}
// if(UL_INF.a[s->id+i].rlen != s->len[i]) {
// fprintf(stderr, "[M::%s] rid:%ld, s->len:%lu, UL_INF->rlen:%u\n", __func__, s->id+i, s->len[i], UL_INF.a[s->id+i].rlen);
// }
assert(UL_INF.a[s->id+i].rlen == s->len[i]);
// void *km = s->buf?(s->buf[tid]?s->buf[tid]->km:NULL):NULL;
// if(s->id+i!=3373) return;
// if(s->id+i!=41927 && s->id+i!=47072 && s->id+i!=67641 && s->id+i!=90305 && s->id+i!=698342 && s->id+i!=329421) {
// return;
// }
// if(s->id+i!=41927) return;
// fprintf(stderr, "\n[M::%s] rid:%ld, len:%lu\n", __func__, s->id+i, s->len[i]);
// if (memcmp(UL_INF.nid.a[s->id+i].a, "d0aab024-b3a7-40fb-83cc-22c3d6d951f8", UL_INF.nid.a[s->id+i].n-1)) return;
// fprintf(stderr, "[M::%s::] ==> len: %lu\n", __func__, s->len[i]);
@@ -5302,6 +5331,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call
// if(l1 == 0 && l2 > 0) fprintf(stderr, "[M::%s::%lu::no_match]\n", UL_INF.nid.a[s->id+i].a, s->len[i]);
// fprintf(stderr, "[M::%s::%lu::] l1->%u; l2->%u\n", UL_INF.nid.a[s->id+i].a, s->len[i], l1, l2);
// fprintf(stderr, "[M::%s::rid->%ld] done\n", __func__, s->id+i);
exit(1);
}
void dump_gaf(mg_gres_a *hits, const mg_gchains_t *gs, uint32_t only_p)
@@ -5683,6 +5713,7 @@ static void *worker_ul_rescall_pipeline(void *data, int step, void *in) // callb
if(UL_INF.a[rid].dd == 0 && p->ucr_s && p->ucr_s->flag == 1) {
assert(s->seq[i]);
// if(s->seq[i] == NULL) fprintf(stderr, "[M::%s::]rid->%ld, len->%lu\n", __func__, rid, s->len[i]);
///for debug interval
write_compress_base_disk(p->ucr_s->fp, rid, s->seq[i], s->len[i], &(p->ucr_s->u));
}
// if(UL_INF.a[rid].dd) fprintf(stderr, "rid->%ld\n", rid);
@@ -8553,6 +8584,11 @@ static void update_ug_arch_ul(void *data, long i, int tid) // callback for kt_fo
}
}
void update_ug_arch_ul_mul(ma_ug_t *ug)
{
kt_for(asm_opt.thread_num, update_ug_arch_ul, ug, ug->g->n_arc);
}
static void filter_short_ulalignments(void *data, long i, int tid) // callback for kt_for()
{
const ma_ug_t *ug = (ma_ug_t *)data;
@@ -8563,10 +8599,13 @@ static void filter_short_ulalignments(void *data, long i, int tid) // callback f
// }
for (k = a_n - 1; k >= 0; k--) {
p = &(a[k]);
// if(i == 4) {
// fprintf(stderr, "[M::%s::k->%ld] p->ts::%u, p->te::%u, p->pchain::%u, p->pidx::%u, p->aidx::%u\n",
// __func__, k, p->ts, p->te, p->pchain, p->pidx, p->aidx);
// if(i == 27512) {
// fprintf(stderr, "+[M::%s::k->%ld::a_n->%ld] p->ts::%u, p->te::%u, p->pchain::%u, p->pidx::%u, p->aidx::%u, p->qs::%u, p->qe::%u\n",
// __func__, k, a_n, p->ts, p->te, p->pchain, p->pidx, p->aidx, p->qs, p->qe);
// }
if(p->qs > p->qe) {
fprintf(stderr, "[M::%s] id->%ld, qs->%u, qe->%u\n", __func__, i, p->qs, p->qe);
}
if(p->base || (!p->el) || (!p->pchain)) continue;
if(p->pidx == (uint32_t)-1) {
if(!ugl_cover_check(p->ts, p->te, &(ug->u.a[p->hid]))) {
@@ -8598,6 +8637,10 @@ static void filter_short_ulalignments(void *data, long i, int tid) // callback f
for (k = a_n - 1; k >= 0; k--) {
p = &(a[k]);
// if(i == 27512) {
// fprintf(stderr, "-[M::%s::k->%ld::a_n->%ld] p->ts::%u, p->te::%u, p->pchain::%u, p->pidx::%u, p->aidx::%u, p->qs::%u, p->qe::%u\n",
// __func__, k, a_n, p->ts, p->te, p->pchain, p->pidx, p->aidx, p->qs, p->qe);
// }
if(p->base || (!p->el) || (!p->pchain)) continue;
if(p->pidx != (uint32_t)-1) {
assert(a[p->pidx].aidx == (uint32_t)k);
@@ -8611,6 +8654,21 @@ static void filter_short_ulalignments(void *data, long i, int tid) // callback f
}
void print_ul_alignment(ma_ug_t *ug, all_ul_t *aln, uint32_t id, const char* cmd)
{
uc_block_t *a = NULL; int64_t k, a_n;
a = aln->a[id].bb.a; a_n = aln->a[id].bb.n;
fprintf(stderr, "\n%s::[M::%s::ul_id->%u::a_n->%ld]\n", cmd, __func__, id, a_n);
for (k = 0; k < a_n; k++) {
fprintf(stderr, "[k->%ld::utg%.6d%c(len->%u)]\tts::%u\tte::%u\t%c\tqs::%u\tqe::%u\tpchain::%u\tpidx::%u\taidx::%u\tpdis::%u\n",
k, a[k].hid + 1, "lc"[ug->u.a[a[k].hid].circ], ug->u.a[a[k].hid].len,
a[k].ts, a[k].te, "+-"[a[k].rev], a[k].qs, a[k].qe, a[k].pchain, a[k].pidx, a[k].aidx, a[k].pdis);
}
}
void filter_ul_ug(ma_ug_t *ug)
{
kt_for(asm_opt.thread_num, filter_short_ulalignments, ug, UL_INF.n);
@@ -10250,13 +10308,20 @@ ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg)
// debug_sl_compress_base_disk_0(&sl, asm_opt.ar);
// detect_outlier_len("ul_realignment");
clear_all_ul_t(&UL_INF);
///for debug interval
if(!load_all_ul_t(&UL_INF, gfa_name, &R_INF, ug)) {
gen_UL_reovlps(&sl, ug, sg, gfa_name, cutoff);
write_all_ul_t(&UL_INF, gfa_name, ug);
}
print_ul_alignment(ug, &UL_INF, 41927, "init-0");
filter_ul_ug(ug);
print_ul_alignment(ug, &UL_INF, 41927, "init-1");
gen_ul_vec_rid_t(&UL_INF, NULL, ug);
kt_for(asm_opt.thread_num, update_ug_arch_ul, ug, ug->g->n_arc);
print_ul_alignment(ug, &UL_INF, 41927, "init-2");
update_ug_arch_ul_mul(ug);
print_ul_alignment(ug, &UL_INF, 41927, "init-3");
// kt_for(asm_opt.thread_num, update_ug_arch_ul, ug, ug->g->n_arc);
// print_all_ul_t_stat(&UL_INF);
// kt_for(sl.n_thread, update_ovlp_src, &sl, R_INF.total_reads);
// kt_for(sl.n_thread, update_ovlp_src_bl, &sl, R_INF.total_reads);
inter.h
+4
View File
@@ -11,5 +11,9 @@ ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg);
int32_t write_all_ul_t(all_ul_t *x, char* file_name, ma_ug_t *ug);
int32_t load_all_ul_t(all_ul_t *x, char* file_name, All_reads *hR, ma_ug_t *ug);
uint32_t ugl_cover_check(uint64_t is, uint64_t ie, ma_utg_t *u);
void filter_ul_ug(ma_ug_t *ug);
void gen_ul_vec_rid_t(all_ul_t *x, All_reads *rdb, ma_ug_t *ug);
void update_ug_arch_ul_mul(ma_ug_t *ug);
void print_ul_alignment(ma_ug_t *ug, all_ul_t *aln, uint32_t id, const char* cmd);
#endif