xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-10-10 07:40:55 +08:00
Code Issues Packages Projects Releases Wiki Activity

version with wrong realignment

Browse Source
This commit is contained in:
chhylp123
2022-11-07 02:35:26 -05:00
parent 612c1ab712
commit 46acc8a2f6
6 changed files with 555 additions and 110 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 <pthread.h>
#include <stdint.h> #include <stdint.h>
#define HA_VERSION "0.17.1-r428" #define HA_VERSION "0.17.2-r432"
#define VERBOSE 0 #define VERBOSE 0
Correct.cpp
+151 -39
View File
@@ -12989,19 +12989,24 @@ int64_t ql, uint64_t rid, ul_ov_t *des)
int64_t push_adp_k_hits(Candidates_list *cl, int64_t cln, uint64_t qs, uint64_t qe, uint64_t ts, uint64_t te, uint64_t readID, int64_t ci) int64_t push_adp_k_hits(Candidates_list *cl, int64_t cln, uint64_t qs, uint64_t qe, uint64_t ts, uint64_t te, uint64_t readID, int64_t ci, int64_t dbgid)
{ {
int64_t k = ci; k_mer_hit *p = NULL, *m = NULL; int64_t k = ci, m = -1; k_mer_hit *p = NULL;
if(cl->length > cln) { if(cl->length > cln) {
m = &(cl->list[cl->length-1]); assert(is_alnw(*m)); m = cl->length-1; assert(is_alnw(cl->list[m]));
} }
for (; k < cln && cl->list[k].readID == readID && cl->list[k].self_offset < qs; k++) { for (; k < cln && cl->list[k].readID == readID && cl->list[k].self_offset < qs; k++) {
p = &(cl->list[k]); p = &(cl->list[k]);
// fprintf(stderr, "[M::%s::] p->q::%ld, p->t::%ld\n", __func__, p->self_offset, p->offset); // fprintf(stderr, "[M::%s::] p->q::%ld, p->t::%ld\n", __func__, p->self_offset, p->offset);
if(p->offset >= ts) continue; if(p->offset >= ts) continue;
if((m) && ((m->offset >= p->offset) || (m->self_offset >= p->self_offset))) continue; if((m>=0) && ((cl->list[m].offset >= p->offset) || (cl->list[m].self_offset >= p->self_offset))) {
continue;
}
kv_pushp_cl(k_mer_hit, (*cl), &p); *p = cl->list[k]; p->readID = p->cnt; kv_pushp_cl(k_mer_hit, (*cl), &p); *p = cl->list[k]; p->readID = p->cnt;
// fprintf(stderr, "[M::%s::+]\n", __func__); // if(dbgid == 109111 || dbgid == 75436) {
// fprintf(stderr, "[M::%s::k->%ld] self_offset::%u, offset::%u, cl->length::%lld, cln::%ld, p->readID::%u\n",
// __func__, k, cl->list[k].self_offset, cl->list[k].offset, cl->length, cln, p->readID);
// }
} }
if(qs == (uint64_t)-1 || qe == (uint64_t)-1) return k; if(qs == (uint64_t)-1 || qe == (uint64_t)-1) return k;
// if(!(k >= cln || cl->list[k].self_offset >= qe)) { // if(!(k >= cln || cl->list[k].self_offset >= qe)) {
@@ -13028,7 +13033,7 @@ int64_t push_adp_k_hits(Candidates_list *cl, int64_t cln, uint64_t qs, uint64_t
int64_t gen_weight_khits0(uint32_t qs, uint32_t qe, k_mer_hit *a, int64_t an, int64_t k, uint64_t dp)///[qs, qe) int64_t gen_weight_khits0(uint32_t qs, uint32_t qe, k_mer_hit *a, int64_t an, int64_t k, uint64_t dp)///[qs, qe)
{ {
for (; k >= 0 && a[k].self_offset >= qs; k--); for (; k >= 0 && a[k].self_offset >= qs; k--);
for (((k>=0)?k:0); k < an && a[k].self_offset < qe; k++) { for (k = ((k>=0)?k:0); k < an && a[k].self_offset < qe; k++) {
if((a[k].self_offset >= qs) && (a[k].self_offset < qe) && (!is_alnw(a[k]))) { if((a[k].self_offset >= qs) && (a[k].self_offset < qe) && (!is_alnw(a[k]))) {
a[k].readID = dp; a[k].readID = dp;
} }
@@ -13055,7 +13060,25 @@ void gen_weight_khits(asg64_v* idx, k_mer_hit *a, int64_t an)
} }
} }
int64_t gen_cns_chain(overlap_region *z, Candidates_list *cl, asg64_v* iidx, int64_t max_lgap, double sgap_rate, int64_t need_filter_khit) void prt_khit(Candidates_list *cl, overlap_region_alloc* ol, overlap_region *z, uint64_t utg_id, const char *cmd)
{
uint64_t k, cid; int64_t i;
if(!z) {
for (k = 0; k < ol->length; k++) {
z = &(ol->list[k]);
if(z->y_id == utg_id) break;
}
}
if((z) || (z->y_id == utg_id)) {
fprintf(stderr, "[M::%s::cmd->%s] ******\n", __func__, (char *)cmd);
for (i = z->shared_seed, cid = cl->list[i].readID; i < cl->length && cl->list[i].readID == cid; i++) {
fprintf(stderr, "[M::%s::i->%lu] qoff::%u, toff::%u, cid::%lu, cl->length::%lld\n", __func__, i,
cl->list[i].self_offset, cl->list[i].offset, cid, cl->length);
}
}
}
int64_t gen_cns_chain(overlap_region_alloc* ol, overlap_region *z, Candidates_list *cl, asg64_v* iidx, int64_t max_lgap, double sgap_rate, int64_t need_filter_khit)
{ {
int64_t k, wn = z->w_list.n, aln_n, qs, qe, ts, te, ci, id, rcn = cl->length, kn; k_mer_hit *ka; int64_t k, wn = z->w_list.n, aln_n, qs, qe, ts, te, ci, id, rcn = cl->length, kn; k_mer_hit *ka;
if(wn <= 0) return 0; if(wn <= 0) return 0;
@@ -13068,29 +13091,41 @@ int64_t gen_cns_chain(overlap_region *z, Candidates_list *cl, asg64_v* iidx, int
qe = z->w_list.a[k].x_end+1; te = z->w_list.a[k].y_end+1; qe = z->w_list.a[k].x_end+1; te = z->w_list.a[k].y_end+1;
} else { } else {
if(qs != -1) { if(qs != -1) {
// fprintf(stderr, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld, cn::%lld\n", // if(z->y_id == 109111 || z->y_id == 75436) {
// __func__, qs, qe, ts, te, ci, cl->length); // fprintf(stderr, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld, cn::%lld\n",
ci = push_adp_k_hits(cl, rcn, qs, qe, ts, te, id, ci); // __func__, qs, qe, ts, te, ci, cl->length);
// }
ci = push_adp_k_hits(cl, rcn, qs, qe, ts, te, id, ci, z->y_id);
aln_n++;//[qs, qe); [ts, te) aln_n++;//[qs, qe); [ts, te)
} }
qs = qe = ts = te = -1; qs = qe = ts = te = -1;
} }
} }
if(qs != -1) { if(qs != -1) {
// fprintf(stderr, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld, cn::%lld\n", // if(z->y_id == 109111 || z->y_id == 75436) {
// __func__, qs, qe, ts, te, ci, cl->length); // fprintf(stderr, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld, cn::%lld\n",
ci = push_adp_k_hits(cl, rcn, qs, qe, ts, te, id, ci); // __func__, qs, qe, ts, te, ci, cl->length);
// }
ci = push_adp_k_hits(cl, rcn, qs, qe, ts, te, id, ci, z->y_id);
aln_n++;//[qs, qe); [ts, te) aln_n++;//[qs, qe); [ts, te)
} }
// fprintf(stderr, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld, cn::%lld\n", // if(z->y_id == 109111 || z->y_id == 75436) {
// __func__, qs, qe, ts, te, ci, cl->length); // fprintf(stderr, "\n[M::%s::] q::[%ld, %ld), t::[%ld, %ld), ci::%ld, cn::%lld, rcn::%ld\n",
push_adp_k_hits(cl, rcn, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, id, ci); // __func__, qs, qe, ts, te, ci, cl->length, rcn);
// }
push_adp_k_hits(cl, rcn, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, id, ci, z->y_id);
ka = cl->list+rcn; kn = cl->length-rcn; ka = cl->list+rcn; kn = cl->length-rcn;
// if(z->y_id == 66) { // if(z->y_id == 66) {
// fprintf(stderr, "[M::%s::] kn::%ld, q_pos::%u, t_pos::%u\n", __func__, kn, ka[kn-1].self_offset, ka[kn-1].offset); // fprintf(stderr, "[M::%s::] kn::%ld, q_pos::%u, t_pos::%u\n", __func__, kn, ka[kn-1].self_offset, ka[kn-1].offset);
// } // }
// if(z->y_id == 109111 || z->y_id == 75436) {
// prt_khit(cl, ol, NULL, 109111, "sa");
// }
if(iidx) gen_weight_khits(iidx, ka, kn); if(iidx) gen_weight_khits(iidx, ka, kn);
// if(z->y_id == 109111 || z->y_id == 75436) {
// prt_khit(cl, ol, NULL, 109111, "sb");
// }
if(need_filter_khit) { if(need_filter_khit) {
kn = lchain_dp_trace(ka, kn, max_lgap, sgap_rate, SGAP); cl->length = rcn + kn; kn = lchain_dp_trace(ka, kn, max_lgap, sgap_rate, SGAP); cl->length = rcn + kn;
} }
@@ -13313,6 +13348,7 @@ void count_k_hits_filter(overlap_region_alloc* ol, Candidates_list *cl, asg64_v*
} }
void count_k_hits_adv(All_reads *rref, const ul_idx_t *uref, char* qstr, UC_Read *buf, void count_k_hits_adv(All_reads *rref, const ul_idx_t *uref, char* qstr, UC_Read *buf,
overlap_region_alloc* ol, Candidates_list *cl, asg64_v* ii, Chain_Data *dp, double sgap_rate, overlap_region_alloc* ol, Candidates_list *cl, asg64_v* ii, Chain_Data *dp, double sgap_rate,
uint64_t khit, uint64_t basec) uint64_t khit, uint64_t basec)
@@ -13326,7 +13362,7 @@ uint64_t khit, uint64_t basec)
refine_khits(&(ol->list[k]), cl, dp, sgap_rate); refine_khits(&(ol->list[k]), cl, dp, sgap_rate);
} }
radix_sort_bc64(ii->a, ii->a+ii->n); radix_sort_bc64(ii->a, ii->a+ii->n);
// prt_khit(cl, ol, NULL, 109111, "a");
for (k = m = 0; k < on; k++) { for (k = m = 0; k < on; k++) {
z = &(ol->list[(uint32_t)ii->a[k]]); i = z->shared_seed; z = &(ol->list[(uint32_t)ii->a[k]]); i = z->shared_seed;
pid = cl->list[i].readID; z->shared_seed = m; pid = cl->list[i].readID; z->shared_seed = m;
@@ -13336,15 +13372,20 @@ uint64_t khit, uint64_t basec)
cl->list[m].cnt = 0; cl->list[m].cnt = 0;
t = cl->list[m].self_offset; t <<= 32; t |= m; t = cl->list[m].self_offset; t <<= 32; t |= m;
kv_push(uint64_t, (*ii), t); kv_push(uint64_t, (*ii), t);
// if(z->y_id == 109111) {
// fprintf(stderr, "[M::%s::m->%ld] qoff::%u, toff::%u\n", __func__, m,
// cl->list[m].self_offset, cl->list[m].offset);
// }
m++; m++;
} }
} }
cl->length = cn = m; srt = ii->a+on; srt_n = ii->n-on; cl->length = cn = m; srt = ii->a+on; srt_n = ii->n-on;
// prt_khit(cl, ol, NULL, 109111, "b");
radix_sort_bc64(srt, srt+srt_n); radix_sort_bc64(srt, srt+srt_n);
if(basec) { if(basec) {
resize_UC_Read(buf, (khit<<1)); str0 = buf->seq; str1 = buf->seq + khit; resize_UC_Read(buf, (khit<<1)); str0 = buf->seq; str1 = buf->seq + khit;
} }
// prt_khit(cl, ol, NULL, 109111, "c");
for (k = 1, l = 0; k <= srt_n; k++) { for (k = 1, l = 0; k <= srt_n; k++) {
if(k == cn || (srt[l]>>32) != (srt[k]>>32)) { if(k == cn || (srt[l]>>32) != (srt[k]>>32)) {
ff = k - l; ff = k - l;
@@ -13360,6 +13401,7 @@ uint64_t khit, uint64_t basec)
l = k; l = k;
} }
} }
// prt_khit(cl, ol, NULL, 109111, "d");
} }
@@ -14436,6 +14478,15 @@ ul_ov_t *res, int64_t wl, int64_t ql, int64_t tl, int64_t ch_i)
// __func__, ch_a[iend-1].self_offset, ch_a[iend-1].offset); // __func__, ch_a[iend-1].self_offset, ch_a[iend-1].offset);
// } // }
// if(z->y_id == 109111) {
// int64_t di;
// fprintf(stderr, "[M::%s] ch_n::%ld, s::%lu, e::%lu, ibeg::%ld, iend::%ld\n", __func__, ch_n, s, e, ibeg, iend);
// for (di = 0; di < ch_n; di++) {
// fprintf(stderr, "[M::%s::di->%ld] qoff::%u, toff::%u\n", __func__, di,
// ch_a[di].self_offset, ch_a[di].offset);
// }
// }
if(ibeg >= 0) { if(ibeg >= 0) {
res->qs = ch_a[ibeg].self_offset; res->qs = ch_a[ibeg].self_offset;
res->ts = ch_a[ibeg].offset; res->ts = ch_a[ibeg].offset;
@@ -14537,10 +14588,15 @@ overlap_region *aux_o, double e_rate)
} else { } else {
is_done = hc_aln_exz_adv(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], mode, wl, exz, ql, e_rate, MAX_SIN_L, MAX_SIN_E, FORCE_SIN_L, -1, aux_o); is_done = hc_aln_exz_adv(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], mode, wl, exz, ql, e_rate, MAX_SIN_L, MAX_SIN_E, FORCE_SIN_L, -1, aux_o);
} }
if(!is_done) { if(!is_done) {
push_unmap_alnw(aux_o, q[0], q[1]-1, t[0], t[1]-1, mode); push_unmap_alnw(aux_o, q[0], q[1]-1, t[0], t[1]-1, mode);
// chain_win_aln(z, dp, cl, q[0], q[1], t[0], t[1], ql, tl, wl, exz); // chain_win_aln(z, dp, cl, q[0], q[1], t[0], t[1], ql, tl, wl, exz);
} }
// if(aux_o->y_id == 109111) {
// fprintf(stderr, "<[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld, is_done::%ld\n",
// __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode, is_done);
// }
} }
return 0; return 0;
} }
@@ -14584,20 +14640,26 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u
if((mode == 0) && is_alnw(ch_a[l]) && is_alnw(ch_a[i]) if((mode == 0) && is_alnw(ch_a[l]) && is_alnw(ch_a[i])
&& (ch_a[l].strand == 0) && (ch_a[i].strand == 1)) { && (ch_a[l].strand == 0) && (ch_a[i].strand == 1)) {
is_done = 1; is_done = 1;
// fprintf(stderr, "*[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n", // if(aux_o->y_id == 109111) {
// __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode); // fprintf(stderr, "*[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n",
// __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode);
// }
// push_replace_alnw(aux_o, q[0], q[1]-1, t[0], t[1]-1, mode);///no need this // push_replace_alnw(aux_o, q[0], q[1]-1, t[0], t[1]-1, mode);///no need this
push_replace_alnw_adv(z, wl, aux_o, q[0], q[1]-1, t[0], t[1]-1, mode); push_replace_alnw_adv(z, wl, aux_o, q[0], q[1]-1, t[0], t[1]-1, mode);
} else if(mode != 3) { } else if(mode != 3) {
if(mode == 1 || mode == 2) adjust_ext_offset(&(q[0]), &(q[1]), &(t[0]), &(t[1]), ql, tl, 0, mode); if(mode == 1 || mode == 2) adjust_ext_offset(&(q[0]), &(q[1]), &(t[0]), &(t[1]), ql, tl, 0, mode);
// fprintf(stderr, "#[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n", // if(aux_o->y_id == 109111) {
// __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode); // fprintf(stderr, "#[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n",
// __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode);
// }
is_done = hc_aln_exz_adv(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], mode, wl, exz, ql, e_rate, MAX_CNS_L, MAX_CNS_E, FORCE_CNS_L, -1, aux_o); is_done = hc_aln_exz_adv(z, uref, hpc_g, rref, qstr, tu, q[0], q[1], t[0], t[1], mode, wl, exz, ql, e_rate, MAX_CNS_L, MAX_CNS_E, FORCE_CNS_L, -1, aux_o);
} }
if(!is_done) {///postprocess if(!is_done) {///postprocess
// fprintf(stderr, ">[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n", // if(aux_o->y_id == 109111) {
// __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode); // fprintf(stderr, ">[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n",
// __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode);
// }
is_done = ovlp_base_aln_all(z, ch_a, ch_n, l, i, uref, hpc_g, rref, qstr, tu, ov, ql, tl, wl, exz, aux_o, e_rate); is_done = ovlp_base_aln_all(z, ch_a, ch_n, l, i, uref, hpc_g, rref, qstr, tu, ov, ql, tl, wl, exz, aux_o, e_rate);
} }
// if(rid == (uint64_t)-1) { // if(rid == (uint64_t)-1) {
@@ -15034,26 +15096,45 @@ UC_Read *tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql,
for (i = 0; i < on; i++) { for (i = 0; i < on; i++) {
// fprintf(stderr, "[M::%s::i->%ld] ovq::[%u, %u), ovt::[%u, %u), hits::[%d, %d)\n", __func__, i, // if(aux_o->y_id == 109111) {
// ov->qs, ov->qe, ov->ts, ov->te, // fprintf(stderr, "[M::%s::i->%ld] ovq::[%u, %u), ovt::[%u, %u), hits::[%d, %d)\n", __func__, i,
// (ov->qn!=((uint32_t)-1))?(int32_t)ov->qn:-1, (int32_t)ov->tn); // ov->qs, ov->qe, ov->ts, ov->te,
// (ov->qn!=((uint32_t)-1))?(int32_t)ov->qn:-1, (int32_t)ov->tn);
// }
assert((i<=0)||(ov[i].qs>ov[i-1].qe)); assert((i<=0)||(ov[i].qs>ov[i-1].qe));
ovlp_base_aln(z, ch_a, ch_n, &(ov[i]), wl, uref, hpc_g, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, rid); ovlp_base_aln(z, ch_a, ch_n, &(ov[i]), wl, uref, hpc_g, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, rid);
} }
int64_t aux_n = aux_o->w_list.n; int64_t aux_n = aux_o->w_list.n;
///for debug
// if(aux_o->y_id == 109111) {
// for (i = 0; i < ((int64_t)aux_o->w_list.n); i++) {
// fprintf(stderr, "-0-[aln::i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, clen::%u, extra_end::%d, mode::%d\n", i,
// aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start,
// aux_o->w_list.a[i].x_start, aux_o->w_list.a[i].x_end+1,
// aux_o->w_list.a[i].y_start, aux_o->w_list.a[i].y_end+1,
// aux_o->w_list.a[i].error, aux_o->w_list.a[i].clen,
// aux_o->w_list.a[i].extra_end, aux_o->w_list.a[i].error_threshold);
// }
// }
for (i = 0; i < aux_n; i++) { for (i = 0; i < aux_n; i++) {
if(!(is_ualn_win(aux_o->w_list.a[i]))) continue; if(!(is_ualn_win(aux_o->w_list.a[i]))) continue;
// if((aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start) <= FORCE_CNS_L) {
// fprintf(stderr, "[aln::i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, clen::%u, mode::%d\n", i,
// aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start,
// aux_o->w_list.a[i].x_start, aux_o->w_list.a[i].x_end+1,
// aux_o->w_list.a[i].y_start, aux_o->w_list.a[i].y_end+1,
// aux_o->w_list.a[i].error, aux_o->w_list.a[i].clen, aux_o->w_list.a[i].error_threshold);
// }
//will overwrite ch_a; does not matter //will overwrite ch_a; does not matter
rechain_aln(z, cl, aux_o, i, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, tl, h_khit, rid); rechain_aln(z, cl, aux_o, i, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, tl, h_khit, rid);
} }
///for debug
// if(aux_o->y_id == 109111) {
// for (i = 0; i < ((int64_t)aux_o->w_list.n); i++) {
// fprintf(stderr, "-2-[aln::i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, clen::%u, mode::%d\n", i,
// aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start,
// aux_o->w_list.a[i].x_start, aux_o->w_list.a[i].x_end+1,
// aux_o->w_list.a[i].y_start, aux_o->w_list.a[i].y_end+1,
// aux_o->w_list.a[i].error, aux_o->w_list.a[i].clen, aux_o->w_list.a[i].error_threshold);
// }
// }
if(((int64_t)aux_o->w_list.n) > aux_n) { if(((int64_t)aux_o->w_list.n) > aux_n) {
for (i = m = 0; i < ((int64_t)aux_o->w_list.n); i++) { for (i = m = 0; i < ((int64_t)aux_o->w_list.n); i++) {
if((i < aux_n) && (is_ualn_win(aux_o->w_list.a[i]))) continue; if((i < aux_n) && (is_ualn_win(aux_o->w_list.a[i]))) continue;
@@ -15779,7 +15860,7 @@ void prt_overlap_region_phase_stat(overlap_region *z)
} }
} }
void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t *uopt, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, asg64_v* buf1) void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t *uopt, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, asg64_v* buf1, int64_t ulid)
{ {
int64_t on = ol->length, k, i, zwn, q[2], t[2], w[2]; int64_t on = ol->length, k, i, zwn, q[2], t[2], w[2];
uint64_t m; overlap_region *z; ul_ov_t *cp; uint64_t m; overlap_region *z; ul_ov_t *cp;
@@ -15787,6 +15868,12 @@ void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t
kv_resize(ul_ov_t, *c_idx, ol->length); kv_resize(ul_ov_t, *c_idx, ol->length);
for (k = idx->n = c_idx->n = 0; k < on; k++) { for (k = idx->n = c_idx->n = 0; k < on; k++) {
z = &(ol->list[k]); zwn = z->w_list.n; z = &(ol->list[k]); zwn = z->w_list.n;
// if(ulid == 35437 && z->y_id == 109111) {
// fprintf(stderr, "+0+[M::%s::utg%.6dl::%c] ulid::%ld, all::%u, non-best::%ld, best::%u\n", __func__,
// (int32_t)z->y_id+1, "+-"[z->y_pos_strand], ulid, z->overlapLen, zwn, z->align_length);
// prt_overlap_region_stat(&(ol->list[k]));
// prt_overlap_region_phase_stat(&(ol->list[k]));
// }
// z->align_length = z->overlapLen = z->x_pos_e+1-z->x_pos_s; // z->align_length = z->overlapLen = z->x_pos_e+1-z->x_pos_s;
z->align_length = 0; z->overlapLen = (uint32_t)-1; z->align_length = 0; z->overlapLen = (uint32_t)-1;
z->non_homopolymer_errors = 0; z->non_homopolymer_errors = 0;
@@ -15832,6 +15919,12 @@ void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t
ovlp_cur_xoff(*cp) = z->w_list.a[w[0]].x_start; ///cur xpos ovlp_cur_xoff(*cp) = z->w_list.a[w[0]].x_start; ///cur xpos
ovlp_cur_coff(*cp) = 0; ///cur cigar off in cur window ovlp_cur_coff(*cp) = 0; ///cur cigar off in cur window
} }
// if(ulid == 35437 && z->y_id == 109111) {
// fprintf(stderr, "+1+[M::%s::utg%.6dl::%c] ulid::%ld, all::%u, non-best::%ld, best::%u\n", __func__,
// (int32_t)z->y_id+1, "+-"[z->y_pos_strand], ulid, z->overlapLen, zwn, z->align_length);
// prt_overlap_region_stat(&(ol->list[k]));
// prt_overlap_region_phase_stat(&(ol->list[k]));
// }
} }
radix_sort_bc64(idx->a, idx->a+idx->n); radix_sort_bc64(idx->a, idx->a+idx->n);
gen_gov_idx(ol, uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, buf1); gen_gov_idx(ol, uref, uopt, G_CHAIN_BW, N_GCHAIN_RATE, buf1);
@@ -15869,6 +15962,15 @@ void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t
for (k = 0; k < on; k++) { for (k = 0; k < on; k++) {
z = &(ol->list[k]); z = &(ol->list[k]);
// if(ulid == 35437 && z->y_id == 109111) {
// fprintf(stderr, "+2+[M::%s::utg%.6dl::%c] ulid::%ld, all::%u, non-best::%ld, best::%u\n", __func__,
// (int32_t)z->y_id+1, "+-"[z->y_pos_strand], ulid, z->overlapLen, zwn, z->align_length);
// prt_overlap_region_stat(&(ol->list[k]));
// prt_overlap_region_phase_stat(&(ol->list[k]));
// }
z->overlapLen = z->x_pos_e+1-z->x_pos_s; z->overlapLen = z->x_pos_e+1-z->x_pos_s;
z->non_homopolymer_errors = 0; zwn = 0; z->non_homopolymer_errors = 0; zwn = 0;
for (i = m = 0; i < z->align_length; i++) { for (i = m = 0; i < z->align_length; i++) {
@@ -15879,7 +15981,14 @@ void region_phase(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t
} }
m++; m++;
} }
z->w_list.n = m; assert(zwn <= z->overlapLen); z->w_list.n = m;
// if(ulid == 35437 && z->y_id == 109111/**zwn > z->overlapLen**/) {
// fprintf(stderr, "+3+[M::%s::utg%.6dl::%c] ulid::%ld, all::%u, non-best::%ld, best::%u\n", __func__,
// (int32_t)z->y_id+1, "+-"[z->y_pos_strand], ulid, z->overlapLen, zwn, z->align_length);
// prt_overlap_region_stat(&(ol->list[k]));
// prt_overlap_region_phase_stat(&(ol->list[k]));
// }
assert(zwn <= z->overlapLen);
z->align_length = z->overlapLen - zwn; z->align_length = z->overlapLen - zwn;
// fprintf(stderr, "[M::%s::utg%.6dl::%c] all::%u, non-best::%ld, best::%u\n", __func__, // fprintf(stderr, "[M::%s::utg%.6dl::%c] all::%u, non-best::%ld, best::%u\n", __func__,
// (int32_t)z->y_id+1, "+-"[z->y_pos_strand], z->overlapLen, zwn, z->align_length); // (int32_t)z->y_id+1, "+-"[z->y_pos_strand], z->overlapLen, zwn, z->align_length);
@@ -15899,14 +16008,17 @@ int64_t max_lgap)
int64_t k, l, ch_n, a_n = aln->n; uint64_t pqn, pk; overlap_region *z; //k_mer_hit *ch_a; int64_t k, l, ch_n, a_n = aln->n; uint64_t pqn, pk; overlap_region *z; //k_mer_hit *ch_a;
// count_k_hits(rref, uref, qstr, tu, ol, cl, buf, khit, base_chekc_k_hit); // count_k_hits(rref, uref, qstr, tu, ol, cl, buf, khit, base_chekc_k_hit);
count_k_hits_adv(rref, uref, qstr, tu, ol, cl, buf, &(cl->chainDP), e_rate, khit, base_chekc_k_hit); count_k_hits_adv(rref, uref, qstr, tu, ol, cl, buf, &(cl->chainDP), e_rate, khit, base_chekc_k_hit);
// prt_khit(cl, ol, NULL, 109111, "e");
for (k = 1, l = 0, pqn = 0; k <= a_n; k++) { for (k = 1, l = 0, pqn = 0; k <= a_n; k++) {
if(k == a_n || aln->a[l].qn != aln->a[k].qn) { if(k == a_n || aln->a[l].qn != aln->a[k].qn) {
z = &(ol->list[aln->a[l].qn]); assert(z->align_length == l); z = &(ol->list[aln->a[l].qn]); assert(z->align_length == l);
// fprintf(stderr, "[M::%s::] oid::[%lu, %u)\n", __func__, pqn, aln->a[l].qn); // fprintf(stderr, "[M::%s::] oid::[%lu, %u)\n", __func__, pqn, aln->a[l].qn);
for (pk = pqn; pk < aln->a[l].qn; pk++) ol->list[pk].w_list.n = 0; for (pk = pqn; pk < aln->a[l].qn; pk++) ol->list[pk].w_list.n = 0;
pqn = aln->a[l].qn+1; pqn = aln->a[l].qn+1;
// fprintf(stderr, "[M::%s::utg%.6dl::%c]\n", __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand]);
ch_n = gen_cns_chain(z, cl, iidx, max_lgap, e_rate, 0); // prt_khit(cl, ol, NULL, 109111, "f");
ch_n = gen_cns_chain(ol, z, cl, iidx, max_lgap, e_rate, 0);
// prt_khit(cl, ol, NULL, 109111, "g");
if(ch_n) { if(ch_n) {
cigar_gen_by_chain_adv(z, cl, cl->length, ch_n, aln->a+l, k-l, wl, uref, hpc_g, rref, qstr, tu, exz, aux_o, e_rate, ql, rid, khit); cigar_gen_by_chain_adv(z, cl, cl->length, ch_n, aln->a+l, k-l, wl, uref, hpc_g, rref, qstr, tu, exz, aux_o, e_rate, ql, rid, khit);
} }
@@ -16233,7 +16345,7 @@ void ul_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur
copy_asg_arr(hap->snp_srt, iidx); copy_asg_arr(v_idx->a, buf); copy_asg_arr((*stb), buf1); copy_asg_arr(hap->snp_srt, iidx); copy_asg_arr(v_idx->a, buf); copy_asg_arr((*stb), buf1);
copy_asg_arr(iidx, hap->snp_srt); copy_asg_arr(buf, v_idx->a); copy_asg_arr(buf1, (*stb)); copy_asg_arr(iidx, hap->snp_srt); copy_asg_arr(buf, v_idx->a); copy_asg_arr(buf1, (*stb));
region_phase(ol, uref, uopt, aln, &iidx, &buf, &buf1); region_phase(ol, uref, uopt, aln, &iidx, &buf, &buf1, sid);
copy_asg_arr(hap->snp_srt, iidx); copy_asg_arr(v_idx->a, buf); copy_asg_arr((*stb), buf1); copy_asg_arr(hap->snp_srt, iidx); copy_asg_arr(v_idx->a, buf); copy_asg_arr((*stb), buf1);
} }
} }
Overlaps.cpp
+1
View File
@@ -930,6 +930,7 @@ void normalize_ma_hit_t_single_side_advance(ma_hit_t_alloc* sources, long long n
} }
} }
typedef struct { typedef struct {
kvec_t_u64_warp *buf; kvec_t_u64_warp *buf;
ma_hit_t_alloc* src; ma_hit_t_alloc* src;
gfa_ut.cpp
+319 -32
View File
@@ -441,7 +441,7 @@ uint32_t get_arcs(asg_t *g, uint32_t v, uint32_t* idx, uint32_t idx_n)
return kv; return kv;
} }
#define flex_arcs0(res, fg, id) ((res) = ((!((id)&(((uint32_t)(0x80000000)))))?(&((fg).g->arc[(id)])):(&((fg).a[(id)])))); #define flex_arcs0(res, fg, id) ((res) = ((!((id)&(((uint32_t)(0x80000000)))))?(&((fg).g->arc[(id)])):(&((fg).a[(id)-(((uint32_t)(0x80000000)))]))));
uint32_t get_flex_arcs(flex_asg_t *fg, uint32_t v, uint32_t* idx, uint32_t idx_n) uint32_t get_flex_arcs(flex_asg_t *fg, uint32_t v, uint32_t* idx, uint32_t idx_n)
{ {
@@ -683,6 +683,38 @@ void recover_contain_g(asg_t *g, ma_hit_t_alloc *src, R_to_U* ruIndex, int64_t m
memset(g->seq_vis, 0, (sizeof(*(g->seq_vis))*(g->n_seq<<1))); memset(g->seq_vis, 0, (sizeof(*(g->seq_vis))*(g->n_seq<<1)));
} }
// static void update_norm_arc(void *data, long i, int tid)
// {
// sset_aux *sl = (sset_aux *)data; ma_hit_t *h, *z; asg_arc_t t;
// ma_hit_t_alloc *src = sl->src; int64_t r, idx;
// uint32_t k, rr, qn, tn, is_u;
// ma_hit_t_alloc *x = &(src[i]);
// for (k = 0; k < x->length; k++) {
// h = &(x->buffer[k]);
// qn = Get_qn((*h)); tn = Get_tn((*h));
// if(h->bl < sl->ul_occ) continue;
// if(g->seq[tn].del) {
// get_R_to_U(ridx, tn, &rr, &is_u);
// if(rr == (uint32_t)-1 || is_u == 1) continue;
// }
// r = ma_hit2arc(h, g->seq[qn].len, g->seq[tn].len, sl->max_hang, asm_opt.max_hang_rate, sl->min_ovlp, &t);
// if(r < 0) continue;
// idx = get_specific_overlap(&(src[tn]), tn, qn);
// z = &(src[tn].buffer[idx]);
// assert(z->bl == h->bl);
// r = ma_hit2arc(z, g->seq[tn].len, g->seq[qn].len, sl->max_hang, asm_opt.max_hang_rate, sl->min_ovlp, &t);
// if(r < 0) continue;
// h->del = 0; if(!(g->seq[tn].del)) z->del = 0;
// g->seq_vis[i] = 1;
// }
// }
// void normalize_ma_hit_t_mul(ma_hit_t_alloc *src, uint32_t n_src)
// {
// sset_aux s; s.src = src;
// kt_for(asm_opt.thread_num, update_norm_arc, &s, n_src);
// }
static void normalize_gou0(void *data, long i, int tid) static void normalize_gou0(void *data, long i, int tid)
{ {
sset_aux *sl = (sset_aux *)data; sset_aux *sl = (sset_aux *)data;
@@ -9982,6 +10014,7 @@ int usg_naive_topocut_aux_sec(usg_t *g, uint32_t v0, int max_ext)
kv++; kv++;
} }
if(kv != 1) break; if(kv != 1) break;
if(v == v0) return 0;///circle, it is ok to remove it
} }
v = v0^1; v = v0^1;
@@ -10007,6 +10040,7 @@ int usg_naive_topocut_aux_sec(usg_t *g, uint32_t v0, int max_ext)
kv++; kv++;
} }
if(kv != 1) break; if(kv != 1) break;
if(v == (v0^1)) return 0;///circle, it is ok to remove it
} }
n_ext -= g->a[v0>>1].occ; n_ext -= g->a[v0>>1].occ;
@@ -10021,12 +10055,14 @@ int usg_naive_topocut_aux_sec(usg_t *g, uint32_t v0, int max_ext)
void usg_arc_cut_length(usg_t *g, asg64_v *in_0, asg64_v *in_1, int32_t max_ext, float len_rat, uint32_t is_trio, void usg_arc_cut_length(usg_t *g, asg64_v *in_0, asg64_v *in_1, int32_t max_ext, float len_rat, uint32_t is_trio,
uint32_t is_topo, uint32_t *max_drop_len) uint32_t is_topo, uint32_t *max_drop_len)
{ {
// fprintf(stderr, "+[M::%s::] max_ext::%d, len_rat::%f\n", __func__, max_ext, len_rat); // if(len_rat > 0.7) {
// fprintf(stderr, "+[M::%s::] max_ext::%d, len_rat::%f\n", __func__, max_ext, len_rat);
// }
asg64_v tx = {0,0,0}, tz = {0,0,0}, *b = NULL, *ub = NULL; asg64_v tx = {0,0,0}, tz = {0,0,0}, *b = NULL, *ub = NULL;
uint32_t i, k, v, w, n_vtx = g->n<<1, nv, nw, kv, kw, /**trioF = (uint32_t)-1, ntrioF = (uint32_t)-1,**/ ol_max, ou_max, to_del, cnt = 0, mm_ol; uint32_t i, k, v, w, n_vtx = g->n<<1, nv, nw, kv, kw, /**trioF = (uint32_t)-1, ntrioF = (uint32_t)-1,**/ ol_max, ou_max, to_del, cnt = 0, mm_ol;
usg_arc_t *av, *aw, *ve, *we; uint64_t x, kocc[2], ou; uint8_t *f; CALLOC(f, g->n); usg_arc_t *av, *aw, *ve, *we; uint64_t x, kocc[2], ou; uint8_t *f; CALLOC(f, g->n);
b = ((in_0)?(in_0):(&tx)); ub = ((in_1)?(in_1):(&tz)); b = ((in_0)?(in_0):(&tx)); ub = ((in_1)?(in_1):(&tz));
for (v = 0, b->n = ub->n = 0; v < n_vtx; ++v) { for (v = 0, b->n = ub->n = 0; v < n_vtx; ++v) {
if (g->a[v>>1].del) continue; if (g->a[v>>1].del) continue;
av = usg_arc_a(g, v); nv = usg_arc_n(g, v); av = usg_arc_a(g, v); nv = usg_arc_n(g, v);
@@ -10050,6 +10086,10 @@ uint32_t is_topo, uint32_t *max_drop_len)
kv_push(uint64_t, *ub, ((((uint64_t)(v))<<32)|((uint64_t)(i)))); kv_push(uint64_t, *ub, ((((uint64_t)(v))<<32)|((uint64_t)(i))));
} }
} }
// if(len_rat > 0.7) {
// fprintf(stderr, "[M::%s::] max_ext::%d, len_rat::%f, b->n::%u, ub->n::%u\n",
// __func__, max_ext, len_rat, (uint32_t)b->n, (uint32_t)ub->n);
// }
radix_sort_srt64(b->a, b->a + b->n); radix_sort_srt64(b->a, b->a + b->n);
for (k = 0; k < b->n; k++) { for (k = 0; k < b->n; k++) {
@@ -10099,6 +10139,10 @@ uint32_t is_topo, uint32_t *max_drop_len)
} }
if (kv <= 1 && kw <= 1) continue; if (kv <= 1 && kw <= 1) continue;
// if(len_rat > 0.7) {
// fprintf(stderr, "0[M::%s::] v>>1::%u(%c), w>>1::%u(%c), kv::%u, kw::%u\n",
// __func__, v>>1, "+-"[v&1], w>>1, "+-"[w&1], kv, kw);
// }
to_del = 1; to_del = 1;
if(is_topo) { if(is_topo) {
@@ -10111,6 +10155,10 @@ uint32_t is_topo, uint32_t *max_drop_len)
if (usg_naive_topocut_aux(g, v^1, max_ext, f, b, ub) < max_ext) to_del = 1; if (usg_naive_topocut_aux(g, v^1, max_ext, f, b, ub) < max_ext) to_del = 1;
} }
} }
// if(len_rat > 0.7) {
// fprintf(stderr, "1[M::%s::] v>>1::%u(%c), w>>1::%u(%c), kv::%u, kw::%u, to_del::%u\n",
// __func__, v>>1, "+-"[v&1], w>>1, "+-"[w&1], kv, kw, to_del);
// }
if (to_del) { if (to_del) {
ve->del = we->del = 1; ve->del = we->del = 1;
@@ -10130,6 +10178,10 @@ uint32_t is_topo, uint32_t *max_drop_len)
} }
} }
// if(len_rat > 0.7) {
// fprintf(stderr, "2[M::%s::] v>>1::%u(%c), w>>1::%u(%c), kv::%u, kw::%u, to_del::%u\n",
// __func__, v>>1, "+-"[v&1], w>>1, "+-"[w&1], kv, kw, to_del);
// }
} }
if(in_0) free(tx.a); if(in_1) free(tz.a); if(in_0) free(tx.a); if(in_1) free(tz.a);
@@ -10836,7 +10888,7 @@ void integer_realign_g(ul_resolve_t *uidx, usg_t *ng, uinfo_srt_warp_t *seq, uin
for (t = 0; t < zt->n; t++) { for (t = 0; t < zt->n; t++) {
j = l + t - zt->n; j = l + t - zt->n;
vj = (zt->a[t]<<1)|(seq->a[z-1].v&1); vj = (zt->a[t]<<1)|(seq->a[z-1].v&1);
if(check_hybrid_connect(ng, seq_id, vj, z-1, vi, z)) { if(check_hybrid_connect(ng, seq_id, vj, z-1, vi, z)) {///seq_id:: integer contig id
sc = csc + (((uint32_t)-1) - (track[j]>>32)); sc = csc + (((uint32_t)-1) - (track[j]>>32));
if(sc > mm_sc) { if(sc > mm_sc) {
mm_sc = sc; mm_idx = j; mm_sc = sc; mm_idx = j;
@@ -10915,7 +10967,9 @@ void integer_realign_g(ul_resolve_t *uidx, usg_t *ng, uinfo_srt_warp_t *seq, uin
} }
buf->u.n = n_u; buf->res_dump.n = l; buf->u.n = n_u; buf->res_dump.n = l;
for (k = 0; k < buf->u.n; k++) { for (k = 0; k < buf->u.n; k++) {
t = (uint32_t)-1; t <<= 32; t |= (((uint32_t)buf->u.a[k])-(buf->u.a[k]>>32)); // t = (uint32_t)-1; t <<= 32;
t = seq_id; t <<= 32; t |= ((uint64_t)0x8000000000000000);
t |= (((uint32_t)buf->u.a[k])-(buf->u.a[k]>>32));
kv_push(uint64_t, buf->res_dump, t); kv_push(uint64_t, buf->res_dump, t);
n_v0 = buf->u.a[k]>>32; n_v = (uint32_t)buf->u.a[k]; n_v0 = buf->u.a[k]>>32; n_v = (uint32_t)buf->u.a[k];
for (z = n_v0; z < n_v; z++) { for (z = n_v0; z < n_v; z++) {
@@ -10951,7 +11005,7 @@ uint64_t get_ug_occ_v(uint32_t i_ug_occ)
return v; return v;
} }
///this function might be wrong
uint32_t usg_unique_arcs_cluster(asg64_v *b64, uint64_t a_n, uint64_t *idx, uint64_t *integ_seq) uint32_t usg_unique_arcs_cluster(asg64_v *b64, uint64_t a_n, uint64_t *idx, uint64_t *integ_seq)
{ {
uint64_t bn = b64->n, k, v; uint64_t bn = b64->n, k, v;
@@ -10959,17 +11013,22 @@ uint32_t usg_unique_arcs_cluster(asg64_v *b64, uint64_t a_n, uint64_t *idx, uint
memset(b64->a + bn, -1, sizeof((*(b64->a)))*a_n); memset(b64->a + bn, -1, sizeof((*(b64->a)))*a_n);
uint64_t *cidx = b64->a + bn, s, e, i, zs, ze, z; uint64_t *cidx = b64->a + bn, s, e, i, zs, ze, z;
///b64.a[0, a_n]:: all resolvable paths with unique beg && end
///b64.a[a_n, bn]:: (raw unitig/non-unqiue node id)|(resolvable path id)
///idx:: the idx for b64.a[a_n, bn]
for (i = 0; i < a_n; i++) {///available interval with beg/end with unique arcs for (i = 0; i < a_n; i++) {///available interval with beg/end with unique arcs
s = b64->a[i]>>32; e = (uint32_t)(b64->a[i]); assert(e > s); s = b64->a[i]>>32; e = (uint32_t)(b64->a[i]); assert(e > s);
if(cidx[i] != (uint64_t)-1) continue; if(cidx[i] != (uint64_t)-1) continue;
for (k = s + 1; k < e; k++) {///note: here is [s, e] for (k = s + 1; k < e; k++) {///note: here is [s, e]
v = integ_seq[k]; v = integ_seq[k];///v is the raw unitig id
zs = (idx[v]<<1)>>33; ze = (uint32_t)idx[v]; zs = (idx[v]<<1)>>33; ze = (uint32_t)idx[v];
for (z = zs; z < ze; z++) { for (z = zs; z < ze; z++) {
// if((b64->a[z]>>32)!=v) { // if((b64->a[z]>>32)!=v) {
// fprintf(stderr, "[M::%s::] v::%lu, b64->a[z]::%lu, zs::%lu, ze::%lu, a_n::%lu\n", // fprintf(stderr, "[M::%s::] v::%lu, b64->a[z]::%lu, zs::%lu, ze::%lu, a_n::%lu\n",
// __func__, v, b64->a[z]>>32, zs, ze, a_n); // __func__, v, b64->a[z]>>32, zs, ze, a_n);
// } // }
///(b64->a[z]>>32):: raw unitig id
///(uint32_t)b64->a[z]:: available interval id
assert((b64->a[z]>>32)==v); assert((b64->a[z]>>32)==v);
cidx[((uint32_t)b64->a[z])] = (i<<32)|((uint32_t)b64->a[z]); cidx[((uint32_t)b64->a[z])] = (i<<32)|((uint32_t)b64->a[z]);
} }
@@ -11001,7 +11060,8 @@ uint32_t usg_unique_arcs_cluster(asg64_v *b64, uint64_t a_n, uint64_t *idx, uint
if(k == a_n || (cidx[i]>>32) != (cidx[k]>>32)) { if(k == a_n || (cidx[i]>>32) != (cidx[k]>>32)) {
for (z = i; z < k; z++) { for (z = i; z < k; z++) {
assert(cidx[z] != (uint64_t)-1); assert(cidx[z] != (uint64_t)-1);
b64->a[b64->n++] = (v<<32)|((uint32_t)cidx[z]);///cluest integer seqs-> (cluster id)|(integer seq id) ///cluest integer seqs-> (cluster id)|(available interval id)
b64->a[b64->n++] = (v<<32)|((uint32_t)cidx[z]);
} }
i = k; v++; i = k; v++;
} }
@@ -11011,6 +11071,104 @@ uint32_t usg_unique_arcs_cluster(asg64_v *b64, uint64_t a_n, uint64_t *idx, uint
return v;///how many cluster return v;///how many cluster
} }
void iter_unique_arcs(asg64_v *buf, asg64_v *b64, uint64_t a_n, uint64_t *idx, uint64_t *integ_seq, uint64_t *cidx, uint64_t i0)
{
uint64_t s, e, x, k, v, zs, ze, z;
buf->n = 0;
kv_push(uint64_t, *buf, i0);
while(buf->n) {
x = kv_pop(*buf);
if(cidx[x] != (uint64_t)-1) continue;
cidx[x] = (i0<<32)|(x);
s = b64->a[x]>>32; e = (uint32_t)(b64->a[x]); assert(e > s);
for (k = s + 1; k < e; k++) {///note: here is [s, e]
v = integ_seq[k];///v is the raw unitig id
zs = (idx[v]<<1)>>33; ze = (uint32_t)idx[v];
for (z = zs; z < ze; z++) {
///(b64->a[z]>>32):: raw unitig id
///(uint32_t)b64->a[z]:: available interval id
assert((b64->a[z]>>32)==v);
if(cidx[((uint32_t)b64->a[z])] != (uint64_t)-1) {
assert((cidx[((uint32_t)b64->a[z])]>>32)==i0);
continue;
}
// cidx[((uint32_t)b64->a[z])] = (i0<<32)|((uint32_t)b64->a[z]);
kv_push(uint64_t, *buf, ((uint32_t)b64->a[z]));
}
v = integ_seq[k]^1;
zs = (idx[v]<<1)>>33; ze = (uint32_t)idx[v];
for (z = zs; z < ze; z++) {
///(b64->a[z]>>32):: raw unitig id
///(uint32_t)b64->a[z]:: available interval id
assert((b64->a[z]>>32)==v);
if(cidx[((uint32_t)b64->a[z])] != (uint64_t)-1) {
assert((cidx[((uint32_t)b64->a[z])]>>32)==i0);
continue;
}
// cidx[((uint32_t)b64->a[z])] = (i0<<32)|((uint32_t)b64->a[z]);
kv_push(uint64_t, *buf, ((uint32_t)b64->a[z]));
}
}
v = integ_seq[s];
zs = (idx[v]<<1)>>33; ze = (uint32_t)idx[v];
for (z = zs; z < ze; z++) {
assert((b64->a[z]>>32)==v);
if(cidx[((uint32_t)b64->a[z])] != (uint64_t)-1) {
assert((cidx[((uint32_t)b64->a[z])]>>32)==i0);
continue;
}
// cidx[((uint32_t)b64->a[z])] = (i0<<32)|((uint32_t)b64->a[z]);
kv_push(uint64_t, *buf, ((uint32_t)b64->a[z]));
}
v = integ_seq[e]^1;
zs = (idx[v]<<1)>>33; ze = (uint32_t)idx[v];
for (z = zs; z < ze; z++) {
assert((b64->a[z]>>32)==v);
if(cidx[((uint32_t)b64->a[z])] != (uint64_t)-1) {
assert((cidx[((uint32_t)b64->a[z])]>>32)==i0);
continue;
}
// cidx[((uint32_t)b64->a[z])] = (i0<<32)|((uint32_t)b64->a[z]);
kv_push(uint64_t, *buf, ((uint32_t)b64->a[z]));
}
}
}
uint32_t usg_unique_arcs_cluster_adv(asg64_v *b64, uint64_t a_n, uint64_t *idx, uint64_t *integ_seq, asg64_v *buf)
{
uint64_t bn = b64->n, k, z, v; buf->n = 0;
kv_resize(uint64_t, *b64, b64->n + a_n); b64->n += a_n;
memset(b64->a + bn, -1, sizeof((*(b64->a)))*a_n);
uint64_t *cidx = b64->a + bn, i;
///b64.a[0, a_n]:: all resolvable paths with unique beg && end
///b64.a[a_n, bn]:: (raw unitig/non-unqiue node id)|(resolvable path id)
///idx:: the idx for b64.a[a_n, bn]
for (i = 0; i < a_n; i++) {///available interval with beg/end with unique arcs
if(cidx[i] != (uint64_t)-1) continue;
iter_unique_arcs(buf, b64, a_n, idx, integ_seq, cidx, i);
}
radix_sort_srt64(cidx, cidx + a_n); b64->n = a_n;
for (i = 0, k = 1, v = 0; k <= a_n; k++) {
if(k == a_n || (cidx[i]>>32) != (cidx[k]>>32)) {
for (z = i; z < k; z++) {
assert(cidx[z] != (uint64_t)-1);
///cluest integer seqs-> (cluster id)|(available interval id)
b64->a[b64->n++] = (v<<32)|((uint32_t)cidx[z]);
}
i = k; v++;
}
}
buf->n = 0;
assert(b64->n == (a_n<<1));
return v;///how many cluster
}
uint32_t ava_pass_unique_bridge(uint64_t *idx, uint64_t *integer_seq, uint64_t s, uint64_t e) uint32_t ava_pass_unique_bridge(uint64_t *idx, uint64_t *integer_seq, uint64_t s, uint64_t e)
{ {
uint64_t k; uint64_t k;
@@ -11030,14 +11188,15 @@ uint32_t ava_pass_unique_bridge(uint64_t *idx, uint64_t *integer_seq, uint64_t s
uint32_t ava_pass_unique_bridge_tips(usg_t *g, asg64_v *b64, uint64_t g_s, uint64_t g_e, uint64_t *integer_seq, uint8_t *f, uint64_t max_ext) uint32_t ava_pass_unique_bridge_tips(usg_t *g, asg64_v *b64, uint64_t g_s, uint64_t g_e, uint64_t *integer_seq, uint8_t *f, uint64_t max_ext)
{ {
uint64_t i, k, z, s, e, v, nv, bn = b64->n, kv, n_ext = 0; usg_arc_t *av; uint64_t i, k, z, s, e, v, nv, bn = b64->n, kv, n_ext = 0; usg_arc_t *av;
for (i = g_s; i < g_e; i++) { for (i = g_s; i < g_e; i++) {///available intervals within the same cluster
s = b64->a[((uint32_t)b64->a[i])]>>32; e = ((uint32_t)b64->a[((uint32_t)b64->a[i])]); assert(s < e); s = b64->a[((uint32_t)b64->a[i])]>>32; e = ((uint32_t)b64->a[((uint32_t)b64->a[i])]); assert(s < e);
for (k = s + 1; k < e; k++) {///note: here is [s, e] for (k = s + 1; k < e; k++) {///note: here is [s, e]
f[integer_seq[k]] = f[integer_seq[k]^1] = 1; f[integer_seq[k]] = f[integer_seq[k]^1] = 1;
} }
f[integer_seq[s]] = f[integer_seq[e]^1] = 1; f[integer_seq[s]] = f[integer_seq[e]^1] = 1;
} }
///collect nodes within raw unitig graph that are linked by the clusters but not in the cluster
for (i = g_s; i < g_e; i++) { for (i = g_s; i < g_e; i++) {
s = b64->a[((uint32_t)b64->a[i])]>>32; e = ((uint32_t)b64->a[((uint32_t)b64->a[i])]); assert(s < e); s = b64->a[((uint32_t)b64->a[i])]>>32; e = ((uint32_t)b64->a[((uint32_t)b64->a[i])]); assert(s < e);
for (k = s + 1; k < e; k++) {///note: here is [s, e] for (k = s + 1; k < e; k++) {///note: here is [s, e]
@@ -11233,7 +11392,7 @@ void update_usg_t_threading_0(usg_t *ng, uint64_t *a, uint64_t a_n, uint32_t *oc
void update_usg_t_threading(ul_resolve_t *uidx, usg_t *ng, uint64_t *arcs, uint64_t *arcs_g, uint64_t arcs_gn, uint64_t *integ_seq, uint32_t *occ, asg64_v *b) void update_usg_t_threading(ul_resolve_t *uidx, usg_t *ng, uint64_t *arcs, uint64_t *arcs_g, uint64_t arcs_gn, uint64_t *integ_seq, uint32_t *occ, asg64_v *b)
{ {
uint64_t i, k, s, e, nvtx = ng->n<<1; memset(occ, 0, sizeof((*occ))*nvtx); uint64_t i, k, s, e, nvtx = ng->n<<1; memset(occ, 0, sizeof((*occ))*nvtx);
for (i = 0; i < arcs_gn; i++) { for (i = 0; i < arcs_gn; i++) {///set cluster
s = arcs[((uint32_t)arcs_g[i])]>>32; e = ((uint32_t)arcs[((uint32_t)arcs_g[i])]); assert(s < e); s = arcs[((uint32_t)arcs_g[i])]>>32; e = ((uint32_t)arcs[((uint32_t)arcs_g[i])]); assert(s < e);
for (k = s + 1; k < e; k++) {///note: here is [s, e] for (k = s + 1; k < e; k++) {///note: here is [s, e]
occ[integ_seq[k]]++; occ[integ_seq[k]^1]++; occ[integ_seq[k]]++; occ[integ_seq[k]^1]++;
@@ -11525,10 +11684,43 @@ ma_ug_t *ma_ug_hybrid_gen(usg_t *g)
void prt_thread_info(uint64_t *interval, uint64_t interval_n, uint64_t *cluster, uint64_t cluster_n,
const char *nn)
{
char* gfa_name = NULL; MALLOC(gfa_name, strlen(nn)+70);
sprintf(gfa_name, "%s.thread_info.log", nn);
FILE* fp = fopen(gfa_name, "w"); free(gfa_name);
if (!fp) return;
uint64_t k;
if(interval) {
for (k = 0; k < interval_n; k++) fprintf(fp,"it_val::[%lu, %u)\n", interval[k]>>32, (uint32_t)interval[k]);
}
if(cluster) {
for (k = 0; k < cluster_n; k++) fprintf(fp,"cluster::%lu\tit_id::%u)\n", cluster[k]>>32, (uint32_t)cluster[k]);
}
fclose(fp);
}
void prt_intg_info(uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, const char *nn)
{
char* gfa_name = NULL; MALLOC(gfa_name, strlen(nn)+70);
sprintf(gfa_name, "%s.intg_info.log", nn);
FILE* fp = fopen(gfa_name, "w"); free(gfa_name);
if (!fp) return;
uint64_t k, i, s, e;
for (i = 0; i < int_idx_n; i++) {///scan all integer contigs
s = int_idx[i]>>32; e = s + ((uint32_t)int_idx[i]);
assert(e > s + 1);//the length is at least 2
fprintf(fp,"idx::[%lu, %lu)\n", s, e);
for (k = s; k < e; k++) {
fprintf(fp,"%lu\n", int_a[k]);
}
}
fclose(fp);
}
void prt_usg_t(ul_resolve_t *uidx, usg_t *ng, const char *cmd);
uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, uint32_t max_ext) uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, uint32_t max_ext)
{ {
@@ -11548,6 +11740,8 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint
} }
ma_ug_t *un_g = ma_ug_hybrid_gen(ng); ma_utg_t *u; ma_ug_t *un_g = ma_ug_hybrid_gen(ng); ma_utg_t *u;
// print_debug_gfa(uidx->sg, un_g, uidx->uopt->coverage_cut, "iig0", uidx->uopt->sources,
// uidx->uopt->ruIndex, uidx->uopt->max_hang, uidx->uopt->min_ovlp, 0, 0, 0);
int32_t ui, un; int32_t ui, un;
for (k = 0; k < un_g->u.n; k++) {///all unitigs of raw utg for (k = 0; k < un_g->u.n; k++) {///all unitigs of raw utg
u = &(un_g->u.a[k]); u = &(un_g->u.a[k]);
@@ -11570,7 +11764,7 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint
ma_ug_destroy(un_g); ma_ug_destroy(un_g);
// fprintf(stderr, ">>>>>>[M::%s::] int_idx[0]::%lu, int_idx[1]::%lu\n", __func__, int_idx[0], int_idx[1]); // fprintf(stderr, ">>>>>>[M::%s::] int_idx[0]::%lu, int_idx[1]::%lu\n", __func__, int_idx[0], int_idx[1]);
// prt_intg_info(int_idx, int_idx_n, int_a, "intg");
for (i = b64.n = ua_n = a_n = 0; i < int_idx_n; i++) {///scan all integer contigs for (i = b64.n = ua_n = a_n = 0; i < int_idx_n; i++) {///scan all integer contigs
s = int_idx[i]>>32; e = s + ((uint32_t)int_idx[i]); s = int_idx[i]>>32; e = s + ((uint32_t)int_idx[i]);
assert(e > s + 1);//the length is at least 2 assert(e > s + 1);//the length is at least 2
@@ -11603,12 +11797,13 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint
} }
} }
assert(a_n + ua_n == b64.n); assert(a_n + ua_n == b64.n);
// prt_thread_info(b64.a, a_n+ua_n, NULL, 0, "tt_minus");
// fprintf(stderr, "**0**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); // fprintf(stderr, "**0**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n);
uint64_t *i_idx, n_clus; CALLOC(i_idx, ng->n<<1); uint64_t *i_idx, n_clus; CALLOC(i_idx, ng->n<<1);
radix_sort_srt64(b64.a, b64.a + b64.n);///keeps the coordinates within int_a[] radix_sort_srt64(b64.a, b64.a + b64.n);///keeps the coordinates within int_a[]
// prt_thread_info(b64.a, a_n, NULL, 0, "tt0");
/*********debugging*********/ /*********debugging*********/
// for (i = 0; i < a_n; i++) { ///available intervals // for (i = 0; i < a_n; i++) { ///available intervals
// s = b64.a[i]>>32; e = (uint32_t)b64.a[i]; // s = b64.a[i]>>32; e = (uint32_t)b64.a[i];
@@ -11641,6 +11836,7 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint
i_idx[int_a[s]] |= ((uint64_t)0x8000000000000000); i_idx[int_a[s]] |= ((uint64_t)0x8000000000000000);
i_idx[int_a[e]^1] |= ((uint64_t)0x8000000000000000); i_idx[int_a[e]^1] |= ((uint64_t)0x8000000000000000);
} }
// prt_thread_info(b64.a, a_n, NULL, 0, "tt1");
///unavailable intervals are useless ///unavailable intervals are useless
b64.n = a_n; assert(ua_n == 0); b64.n = a_n; assert(ua_n == 0);
for (i = 0; i < a_n; i++) { ///available intervals for (i = 0; i < a_n; i++) { ///available intervals
@@ -11648,7 +11844,7 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint
assert(!(b64.a[i]&((uint64_t)0x8000000000000000))); assert(!(b64.a[i]&((uint64_t)0x8000000000000000)));
for (k = s + 1; k < e; k++) {///note: here is [s, e]; s && e are unique, but [s+1, e-1] are not unique for (k = s + 1; k < e; k++) {///note: here is [s, e]; s && e are unique, but [s+1, e-1] are not unique
kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[k]]++; kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[k]]++;
(*pz) = int_a[k]; (*pz) <<= 32; (*pz) |= i;///(raw unitig node id)|(integer contig id) (*pz) = int_a[k]; (*pz) <<= 32; (*pz) |= i;///(raw unitig/non-unqiue node id)|(integer contig id)
kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[k]^1]++; kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[k]^1]++;
(*pz) = int_a[k]^1; (*pz) <<= 32; (*pz) |= i; (*pz) = int_a[k]^1; (*pz) <<= 32; (*pz) |= i;
@@ -11659,23 +11855,31 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint
kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[e]^1]++; kv_pushp(uint64_t, b64, &pz); //i_idx[int_a[e]^1]++;
(*pz) = int_a[e]^1; (*pz) <<= 32; (*pz) |= i; (*pz) = int_a[e]^1; (*pz) <<= 32; (*pz) |= i;
} }
// prt_thread_info(b64.a, a_n, NULL, 0, "tt2");
// fprintf(stderr, "**1**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); // fprintf(stderr, "**1**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n);
///index ///index
radix_sort_srt64(b64.a + a_n, b64.a + b64.n); radix_sort_srt64(b64.a + a_n, b64.a + b64.n);///(raw unitig node id)|(integer contig id)
for (k = a_n + 1, i = a_n; k <= b64.n; k++) { for (k = a_n + 1, i = a_n; k <= b64.n; k++) {
if(k == b64.n || (b64.a[k]>>32) != (b64.a[i]>>32)) { if(k == b64.n || (b64.a[k]>>32) != (b64.a[i]>>32)) {
i_idx[b64.a[i]>>32] |= (((uint64_t)i)<<32)|((uint64_t)k); i_idx[b64.a[i]>>32] |= (((uint64_t)i)<<32)|((uint64_t)k);///b64.a[i]>>32 appear once (unique ends)/multipe times
i = k; i = k;
} }
} }
// prt_thread_info(b64.a, a_n, NULL, 0, "tt3");
// fprintf(stderr, "**2**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n); // fprintf(stderr, "**2**[M::%s::] a_n::%lu, ua_n::%lu\n", __func__, a_n, ua_n);
n_clus = usg_unique_arcs_cluster(&b64, a_n, i_idx, int_a); ///b64.a[0, a_n]:: all resolvable paths with unique beg && end
///b64.a[a_n, b64.n]:: (raw unitig/non-unqiue node id)|(resolvable path id)
// n_clus = usg_unique_arcs_cluster(&b64, a_n, i_idx, int_a);///this function might be wrong
n_clus = usg_unique_arcs_cluster_adv(&b64, a_n, i_idx, int_a, &ub64);
// fprintf(stderr, "**3**[M::%s::] a_n::%lu, n_clus::%lu\n", __func__, a_n, n_clus); // fprintf(stderr, "**3**[M::%s::] a_n::%lu, n_clus::%lu\n", __func__, a_n, n_clus);
// prt_thread_info(b64.a, a_n, NULL, 0, "tt4");
assert(b64.n == (a_n<<1)); assert(b64.n == (a_n<<1));
for (k = a_n + 1, i = a_n, mm = a_n; k <= b64.n; k++) { for (k = a_n + 1, i = a_n, mm = a_n; k <= b64.n; k++) {
if(k == b64.n || (b64.a[k]>>32) != (b64.a[i]>>32)) { if(k == b64.n || (b64.a[k]>>32) != (b64.a[i]>>32)) {
for (z = i; z < k; z++) {///all intger seqs within the same cluster for (z = i; z < k; z++) {///all intger seqs within the same cluster
s = b64.a[((uint32_t)b64.a[z])]>>32; e = ((uint32_t)b64.a[((uint32_t)b64.a[z])]); assert(e > s); s = b64.a[((uint32_t)b64.a[z])]>>32; e = ((uint32_t)b64.a[((uint32_t)b64.a[z])]); assert(e > s);
///[s, e]:: available interval
if(!ava_pass_unique_bridge(i_idx, int_a, s, e)) break; if(!ava_pass_unique_bridge(i_idx, int_a, s, e)) break;
} }
if(z >= k) {///all arcs in this cluster is fine -> each of arch is reliable if(z >= k) {///all arcs in this cluster is fine -> each of arch is reliable
@@ -11685,6 +11889,7 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint
} }
} }
assert(n_clus == 0); assert(n_clus == 0);
// prt_thread_info(b64.a, a_n, NULL, 0, "tt5");
// fprintf(stderr, "**4**[M::%s::] a_n::%lu, n_clus::%lu\n", __func__, a_n, n_clus); // fprintf(stderr, "**4**[M::%s::] a_n::%lu, n_clus::%lu\n", __func__, a_n, n_clus);
b64.n = mm; n_clus = 0; b64.n = mm; n_clus = 0;
for (k = a_n + 1, i = a_n, mm = a_n; k <= b64.n; k++) { for (k = a_n + 1, i = a_n, mm = a_n; k <= b64.n; k++) {
@@ -11696,9 +11901,10 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint
i = k; i = k;
} }
} }
// prt_thread_info(b64.a, a_n, NULL, 0, "tt6");
b64.n = mm; b64.n = mm;
// fprintf(stderr, "**5**[M::%s::] a_n::%lu, n_clus::%lu, ng->n::%lu\n", __func__, a_n, n_clus, ng->n); // fprintf(stderr, "**5**[M::%s::] a_n::%lu, n_clus::%lu, ng->n::%lu\n", __func__, a_n, n_clus, ng->n);
// prt_thread_info(b64.a, a_n, b64.a + a_n, b64.n - a_n, "thred");
if(n_clus > 0) { if(n_clus > 0) {
update_usg_t_threading(uidx, ng, b64.a, b64.a + a_n, b64.n - a_n, int_a, ng_occ, &ub64); update_usg_t_threading(uidx, ng, b64.a, b64.a + a_n, b64.n - a_n, int_a, ng_occ, &ub64);
} }
@@ -11711,24 +11917,68 @@ uint64_t gen_unique_g_adv(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint
void u2g_hybrid_detan(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, asg64_v *in, asg64_v *ib) void u2g_hybrid_detan(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, asg64_v *in, asg64_v *ib)
{ {
uint64_t k, i, x; asg64_v tx = {0,0,0}, tb = {0,0,0}, *ob = NULL, *ub = NULL; uint64_t k, i, x, m, *tmp, sn; asg64_v tx = {0,0,0}, tb = {0,0,0}, *ob = NULL, *ub = NULL;
ob = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); ob->n = ub->n = 0; ob = (in?(in):(&tx)); ub = (ib?(ib):(&tb)); ob->n = ub->n = 0;
for (k = 0; k < uidx->str_b.n_thread; k++) { for (k = 0; k < uidx->str_b.n_thread; k++) {
uidx->str_b.buf[k].res_dump.n = uidx->str_b.buf[k].u.n = uidx->str_b.buf[k].o.n = 0; uidx->str_b.buf[k].res_dump.n = uidx->str_b.buf[k].u.n = uidx->str_b.buf[k].o.n = 0;
} }
kt_for(uidx->str_b.n_thread, worker_integer_realign_g, uidx, uidx->uovl.i_ug->u.n); kt_for(uidx->str_b.n_thread, worker_integer_realign_g, uidx, uidx->uovl.i_ug->u.n);
for (k = ob->n = ub->n = 0; k < uidx->str_b.n_thread; k++) {
for (k = ob->n = ub->n = m = 0; k < uidx->str_b.n_thread; k++) {
for (i = 0; i < uidx->str_b.buf[k].res_dump.n; i++) { for (i = 0; i < uidx->str_b.buf[k].res_dump.n; i++) {
if((uidx->str_b.buf[k].res_dump.a[i]>>32)==((uint32_t)-1)) { x = uidx->str_b.buf[k].res_dump.a[i];
x = ob->n; x <<= 32; x |= ((uint32_t)uidx->str_b.buf[k].res_dump.a[i]); kv_push(uint64_t, *ob, x);//aln details
kv_push(uint64_t, *ub, x); if(x&((uint64_t)0x8000000000000000)) {
x -= ((uint64_t)0x8000000000000000); x >>= 32; x <<= 32;//seq_id
x |= ob->n;//offset
kv_push(uint64_t, *ub, x);///idx:: seq_id|offset_in_ob
} else { } else {
kv_push(uint64_t, *ob, uidx->str_b.buf[k].res_dump.a[i]); m++;
} }
} }
} }
radix_sort_srt64(ub->a, ub->a + ub->n);
kv_resize(uint64_t, *ub, ub->n+m); tmp = ub->a + ub->n; m = 0;
for (k = 0; k < ub->n; k++) {
sn = ((uint32_t)(ob->a[((uint32_t)ub->a[k])-1]));
memcpy(tmp + m, ob->a + ((uint32_t)ub->a[k]), sn*sizeof((*tmp)));
ub->a[k] = m; ub->a[k] <<= 32; ub->a[k] |= sn;//offset_in_ob|occ
m += sn;
}
assert(m <= ob->n);
memcpy(ob->a, tmp, m*sizeof((*tmp))); ob->n = m;
// for (k = 0, p = NULL; k < ub->n; k++) {
// p = &(ob->a[((uint32_t)ub->a[k])-1]);
// x = ub->a[k]<<32;//offset
// x |= ((uint32_t)(*p));///occ
// ub->a[k] = x;
// (*p) >>= 32; (*p) <<= 32; (*p) |= k;
// }
// for (k = m = 0; k < ob->n; k++) {
// if(ob->a[k]&((uint64_t)0x8000000000000000)) {
// x = m; x <<= 32; x |= ((uint32_t)ub->a[(uint32_t)ob->a[k]]);
// ub->a[(uint32_t)ob->a[k]] = x;
// } else {
// ob->a[m++] = ob->a[k];
// }
// }
// ob->n = m;
// for (k = ob->n = ub->n = 0; k < uidx->str_b.n_thread; k++) {
// for (i = 0; i < uidx->str_b.buf[k].res_dump.n; i++) {
// if((uidx->str_b.buf[k].res_dump.a[i]>>32)==((uint32_t)-1)) {
// x = ob->n; x <<= 32; x |= ((uint32_t)uidx->str_b.buf[k].res_dump.a[i]);
// kv_push(uint64_t, *ub, x);///idx:: offset_in_ob|occ
// } else {
// kv_push(uint64_t, *ob, uidx->str_b.buf[k].res_dump.a[i]);//aln details
// }
// }
// }
/*********debugging*********/ /*********debugging*********/
fprintf(stderr, "[M::%s::] # iug::%u, # gchain::%u\n", __func__, (uint32_t)uidx->uovl.i_ug->u.n, (uint32_t)ub->n); fprintf(stderr, "[M::%s::] # iug::%u, # gchain::%u\n", __func__, (uint32_t)uidx->uovl.i_ug->u.n, (uint32_t)ub->n);
// for (k = 0; k < ub->n; k++) { // for (k = 0; k < ub->n; k++) {
@@ -11747,6 +11997,7 @@ void u2g_hybrid_detan(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, asg64_v *
// u_ug = gen_unique_g(uidx, ng, ng_occ, ub->a, ub->n, ob->a);//ma_ug_hybrid_gen(ng); // u_ug = gen_unique_g(uidx, ng, ng_occ, ub->a, ub->n, ob->a);//ma_ug_hybrid_gen(ng);
///debug ///debug
debug_sysm_usg_t(ng, __func__); debug_sysm_usg_t(ng, __func__);
// prt_usg_t(uidx, ng, "ng4");
if(gen_unique_g_adv(uidx, ng, ub->a, ub->n, ob->a, max_ext)) { if(gen_unique_g_adv(uidx, ng, ub->a, ub->n, ob->a, max_ext)) {
// usg_arc_t *z = get_usg_arc(ng, 2, 576), *q = get_usg_arc(ng, 577, 3); // usg_arc_t *z = get_usg_arc(ng, 2, 576), *q = get_usg_arc(ng, 577, 3);
// fprintf(stderr, "xxxx0xxx[M::%s::] p->del::%u, q->del::%u\n", // fprintf(stderr, "xxxx0xxx[M::%s::] p->del::%u, q->del::%u\n",
@@ -11762,6 +12013,7 @@ void u2g_hybrid_detan(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, asg64_v *
debug_sysm_usg_t(ng, __func__); debug_sysm_usg_t(ng, __func__);
// fprintf(stderr, "-[M::%s::] ng->n::%u\n", __func__, (uint32_t)ng->n); // fprintf(stderr, "-[M::%s::] ng->n::%u\n", __func__, (uint32_t)ng->n);
} }
// prt_usg_t(uidx, ng, "ng_dbg");
if(!in) free(tx.a); if(!ib) free(tb.a); if(!in) free(tx.a); if(!ib) free(tb.a);
} }
@@ -11825,31 +12077,65 @@ void u2g_hybrid_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, usg_t *ng, asg64_v *
((ulopt->max_ovlp_drop_ratio-ulopt->min_ovlp_drop_ratio)/(ulopt->clean_round-1))); ((ulopt->max_ovlp_drop_ratio-ulopt->min_ovlp_drop_ratio)/(ulopt->clean_round-1)));
double drop = ulopt->min_ovlp_drop_ratio; ///CALLOC(iug->g->seq_vis, iug->g->n_seq*2); double drop = ulopt->min_ovlp_drop_ratio; ///CALLOC(iug->g->seq_vis, iug->g->n_seq*2);
fprintf(stderr, "\n[M::%s::] Starting hybrid clean, mm_tip::%ld\n", __func__, mm_tip); fprintf(stderr, "\n[M::%s::] Starting hybrid clean, mm_tip::%ld\n", __func__, mm_tip);
// prt_usg_t(uidx, ng, "ng0");
usg_arc_cut_tips(ng, mm_tip, 0, b); usg_arc_cut_tips(ng, mm_tip, 0, b);
// prt_usg_t(uidx, ng, "ng1");
// char sb[1000];
for (ss = 1; ss <= 1/**6**/; ss++) { for (ss = 1; ss <= 1/**6**/; ss++) {
mm_tip = ulopt->max_tip_hifi*ss; mm_tip = ulopt->max_tip_hifi*ss;
fprintf(stderr, "\n[M::%s::] ss::%ld, mm_tip::%ld, ulopt->clean_round::%ld\n",
__func__, ss, mm_tip, ulopt->clean_round);
for (i = 0, drop = ulopt->min_ovlp_drop_ratio; i < ulopt->clean_round; i++, drop += step) { for (i = 0, drop = ulopt->min_ovlp_drop_ratio; i < ulopt->clean_round; i++, drop += step) {
if(drop > ulopt->max_ovlp_drop_ratio) drop = ulopt->max_ovlp_drop_ratio; if(drop > ulopt->max_ovlp_drop_ratio) drop = ulopt->max_ovlp_drop_ratio;
// fprintf(stderr, "-0-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop);
// sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_a", ss, i, drop);
// prt_usg_t(uidx, ng, sb);
usg_arc_cut_length(ng, b, ub, mm_tip>>1, drop, ulopt->is_trio, 1, NULL); usg_arc_cut_length(ng, b, ub, mm_tip>>1, drop, ulopt->is_trio, 1, NULL);
// fprintf(stderr, "-1-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop);
// sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_b", ss, i, drop);
// prt_usg_t(uidx, ng, sb);
usg_arc_cut_tips(ng, mm_tip, 0, b); usg_arc_cut_tips(ng, mm_tip, 0, b);
// fprintf(stderr, "-2-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop);
// sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_c", ss, i, drop);
// prt_usg_t(uidx, ng, sb);
usg_arc_cut_length(ng, b, ub, mm_tip, drop, ulopt->is_trio, 1, NULL); usg_arc_cut_length(ng, b, ub, mm_tip, drop, ulopt->is_trio, 1, NULL);
// fprintf(stderr, "-3-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop);
// sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_d", ss, i, drop);
// prt_usg_t(uidx, ng, sb);
usg_arc_cut_tips(ng, mm_tip, 1, b); usg_arc_cut_tips(ng, mm_tip, 1, b);
// sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_e", ss, i, drop);
// prt_usg_t(uidx, ng, sb);
// fprintf(stderr, "-4-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop);
} }
drop = 1; drop = 1;
// fprintf(stderr, "-0-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop);
// sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_a", ss, i, drop);
// prt_usg_t(uidx, ng, sb);
usg_arc_cut_length(ng, b, ub, mm_tip>>1, drop, ulopt->is_trio, 1, NULL); usg_arc_cut_length(ng, b, ub, mm_tip>>1, drop, ulopt->is_trio, 1, NULL);
// fprintf(stderr, "-1-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop);
// sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_b", ss, i, drop);
// prt_usg_t(uidx, ng, sb);
usg_arc_cut_tips(ng, mm_tip, 0, b); usg_arc_cut_tips(ng, mm_tip, 0, b);
// fprintf(stderr, "-2-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop);
// sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_c", ss, i, drop);
// prt_usg_t(uidx, ng, sb);
usg_arc_cut_length(ng, b, ub, mm_tip, drop, ulopt->is_trio, 1, NULL); usg_arc_cut_length(ng, b, ub, mm_tip, drop, ulopt->is_trio, 1, NULL);
// fprintf(stderr, "-3-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop);
// sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_d", ss, i, drop);
// prt_usg_t(uidx, ng, sb);
usg_arc_cut_tips(ng, mm_tip, 1, b); usg_arc_cut_tips(ng, mm_tip, 1, b);
// sprintf(sb, "ng_ss::%ld_i::%ld_drop::%f_e", ss, i, drop);
// prt_usg_t(uidx, ng, sb);
// fprintf(stderr, "-4-[M::%s::] i::%ld, drop::%f\n", __func__, i, drop);
} }
// prt_usg_t(uidx, ng, "ng2");
///debug ///debug
debug_sysm_usg_t(ng, __func__); debug_sysm_usg_t(ng, __func__);
/******for debug******/ /******for debug******/
// prt_usg_t(uidx, ng, "ng1"); prt_usg_t(uidx, ng, "ng_dbg");
/******for debug******/ /******for debug******/
// u2g_hybrid_extend(ng, NULL, b, ub); // u2g_hybrid_extend(ng, NULL, b, ub);
@@ -11928,9 +12214,9 @@ ma_ug_t *gen_hybrid_ug(ul_resolve_t *uidx, usg_t *ng)
ug->g->seq[i].c = PRIMARY_LABLE; ug->g->seq[i].c = PRIMARY_LABLE;
u = &(ug->u.a[i]); u = &(ug->u.a[i]);
if(u->m == 0) continue; if(u->m == 0) continue;
fprintf(stderr, "+[M::%s::] i::%u\n", __func__, i); // fprintf(stderr, "+[M::%s::] i::%u\n", __func__, i);
merge_hybrid_utg_content(u, uidx->l1_ug, uidx->sg, ng, &e); merge_hybrid_utg_content(u, uidx->l1_ug, uidx->sg, ng, &e);
fprintf(stderr, "-[M::%s::] i::%u\n", __func__, i); // fprintf(stderr, "-[M::%s::] i::%u\n", __func__, i);
ug->g->seq[i].len = u->len; ug->g->seq[i].len = u->len;
} }
kv_destroy(e.a); kv_destroy(e.a);
@@ -12690,6 +12976,7 @@ double max_ovlp_drop_ratio, int64_t max_tip, bub_label_t *b_mask_t, uint32_t is_
print_raw_uls_seq(uidx, asm_opt.output_file_name); print_raw_uls_seq(uidx, asm_opt.output_file_name);
ul_re_correct(uidx, 3); ul_re_correct(uidx, 3);
init_ulg_opt_t(&uu, uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.55, max_tip, max_tip<<1, b_mask_t, is_trio); init_ulg_opt_t(&uu, uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, 0.55, max_tip, max_tip<<1, b_mask_t, is_trio);
print_debug_gfa(sg, init_ug, uopt->coverage_cut, "UL.debug0", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1);
/**ul2ul_idx_t *u2o = **/gen_ul2ul(uidx, uopt, &uu, 0); /**ul2ul_idx_t *u2o = **/gen_ul2ul(uidx, uopt, &uu, 0);
// print_ul_alignment(init_ug, &UL_INF, 47072, "after-3"); // print_ul_alignment(init_ug, &UL_INF, 47072, "after-3");
inter.cpp
+82 -37
View File
@@ -5184,7 +5184,7 @@ const ul_idx_t *uref, const ug_opt_t *uopt, int64_t bw, double diff_ec_ul, uint6
if(lj->qs >= li->qs) return INT32_MIN; if(lj->qs >= li->qs) return INT32_MIN;
uint32_t li_v = (li->tn<<1)|li->rev, lj_v = (lj->tn<<1)|lj->rev; uint32_t li_v = (li->tn<<1)|li->rev, lj_v = (lj->tn<<1)|lj->rev;
int64_t qo = infer_rovlp(li, lj, NULL, NULL, ridx, ug), trans_l = 0, sec_err = 0, sc; ///overlap length in query (UL read) int64_t qo = infer_rovlp(li, lj, NULL, NULL, ridx, ug), trans_l = 0, sec_err = 0, sc; ///overlap length in query (UL read)
if(li_v != lj_v && get_ecov_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, mode, NULL)) { if(/**li_v != lj_v &&**/ get_ecov_adv(uref, uopt, li_v^1, lj_v^1, bw, diff_ec_ul, qo, mode, NULL)) {
trans_l = get_overlap_region_sub_err(&(ol[li->qn]), tc, lj->qe, &sec_err); trans_l = get_overlap_region_sub_err(&(ol[li->qn]), tc, lj->qe, &sec_err);
// int64_t trans_l_debug, sec_err_debug; // int64_t trans_l_debug, sec_err_debug;
@@ -5445,7 +5445,7 @@ int64_t *t, int64_t n_ext)
if(lj->qs >= li->qs) continue; if(lj->qs >= li->qs) continue;
set_ul_ov_t_by_mg_lchain_t(&uj, lj); set_ul_ov_t_by_mg_lchain_t(&uj, lj);
qo = infer_rovlp(&ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug); ///overlap length in query (UL read) qo = infer_rovlp(&ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug); ///overlap length in query (UL read)
if(li->v!=lj->v && get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, ng_diff_thre, qo, 0, NULL)) { if(/**li->v!=lj->v &&**/get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, ng_diff_thre, qo, 0, NULL)) {
sc = csc + f[j]; sc = csc + f[j];
if(sc > mm_sc) { if(sc > mm_sc) {
mm_sc = sc, mm_idx = j; mm_sc = sc, mm_idx = j;
@@ -5473,7 +5473,7 @@ int64_t *t, int64_t n_ext)
if(lj->qe+G_CHAIN_INDEL > li->qs && lj->qs < li->qs) { if(lj->qe+G_CHAIN_INDEL > li->qs && lj->qs < li->qs) {
set_ul_ov_t_by_mg_lchain_t(&uj, lj); set_ul_ov_t_by_mg_lchain_t(&uj, lj);
qo = infer_rovlp(&ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug);///overlap length in query (UL read) qo = infer_rovlp(&ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug);///overlap length in query (UL read)
if(li->v!=lj->v && get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, ng_diff_thre, qo, 0, NULL)) { if(/**li->v!=lj->v &&**/get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, ng_diff_thre, qo, 0, NULL)) {
sc = csc + f[max_ii]; sc = csc + f[max_ii];
if(sc > mm_sc) { if(sc > mm_sc) {
mm_sc = sc; mm_idx = max_ii; mm_sc = sc; mm_idx = max_ii;
@@ -5651,7 +5651,7 @@ st_mt_t *bf, Chain_Data* dp, int64_t max_skip, int64_t max_iter, int64_t max_dis
if((!is_f) && (lj->qe+G_CHAIN_INDEL > li->qs)) { if((!is_f) && (lj->qe+G_CHAIN_INDEL > li->qs)) {
set_ul_ov_t_by_mg_lchain_t(&uj, lj); set_ul_ov_t_by_mg_lchain_t(&uj, lj);
qo = infer_rovlp(&ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug); qo = infer_rovlp(&ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug);
if(li->v!=lj->v && get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, N_GCHAIN_RATE, qo, 0, NULL)) { if(/**li->v!=lj->v &&**/ get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, N_GCHAIN_RATE, qo, 0, NULL)) {
is_f = 1; if(n_skip > 0) n_skip--; is_f = 1; if(n_skip > 0) n_skip--;
if(n_skip < (max_skip>>1)) n_skip= (max_skip>>1); if(n_skip < (max_skip>>1)) n_skip= (max_skip>>1);
} }
@@ -5882,7 +5882,7 @@ int64_t need_srt)
if((!is_f) && (lj->qe+G_CHAIN_INDEL > li->qs)) { if((!is_f) && (lj->qe+G_CHAIN_INDEL > li->qs)) {
set_ul_ov_t_by_mg_lchain_t(&uj, lj); set_ul_ov_t_by_mg_lchain_t(&uj, lj);
qo = infer_rovlp(&ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug); qo = infer_rovlp(&ui, &uj, NULL, NULL, NULL, (ma_ug_t *)ug);
if(li->v!=lj->v && get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, N_GCHAIN_RATE, qo, 0, NULL)) { if(/**li->v!=lj->v &&**/ get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, N_GCHAIN_RATE, qo, 0, NULL)) {
is_f = 1; if(n_skip > 0) n_skip--; is_f = 1; if(n_skip > 0) n_skip--;
if(n_skip < (max_skip>>1)) n_skip= (max_skip>>1); if(n_skip < (max_skip>>1)) n_skip= (max_skip>>1);
} }
@@ -6095,9 +6095,12 @@ const asg_t *g, st_mt_t *dst_done, vec_sp_node_t *out, vec_mg_pathv_t *res, vec_
if(res->a[i].v == (uint32_t)-1) { if(res->a[i].v == (uint32_t)-1) {
gchains->a[i+gchains->n] = a[res->a[i].pre + g_item->rs]; gchains->a[i+gchains->n] = a[res->a[i].pre + g_item->rs];
gchains->a[i+gchains->n].dist_pre = res->a[i].d; gchains->a[i+gchains->n].dist_pre = res->a[i].d;
// if(ulid == 14714) {
// fprintf(stderr, "+[M::%s::]\tutg%.6dl(%c)\n", __func__, // fprintf(stderr, "+[M::%s::]\tutg%.6dl(%c)\tq::[%d, %d)\tt::[%d, %d)\n", __func__,
// (int32_t)(gchains->a[i+gchains->n].v>>1)+1, "+-"[gchains->a[i+gchains->n].v&1]); // (int32_t)(gchains->a[i+gchains->n].v>>1)+1, "+-"[gchains->a[i+gchains->n].v&1],
// gchains->a[i+gchains->n].qs, gchains->a[i+gchains->n].qe,
// gchains->a[i+gchains->n].rs, gchains->a[i+gchains->n].re);
// }
} else { } else {
gchains->a[i+gchains->n].v = res->a[i].v; gchains->a[i+gchains->n].v = res->a[i].v;
gchains->a[i+gchains->n].off = -1; gchains->a[i+gchains->n].off = -1;
@@ -6108,6 +6111,12 @@ const asg_t *g, st_mt_t *dst_done, vec_sp_node_t *out, vec_mg_pathv_t *res, vec_
// fprintf(stderr, "aaaaaaa, ulid->%ld\n", ulid); // fprintf(stderr, "aaaaaaa, ulid->%ld\n", ulid);
// fprintf(stderr, "-[M::%s::]\tutg%.6dl(%c)\n", __func__, // fprintf(stderr, "-[M::%s::]\tutg%.6dl(%c)\n", __func__,
// (int32_t)(gchains->a[i+gchains->n].v>>1)+1, "+-"[gchains->a[i+gchains->n].v&1]); // (int32_t)(gchains->a[i+gchains->n].v>>1)+1, "+-"[gchains->a[i+gchains->n].v&1]);
// if(ulid == 14714) {
// fprintf(stderr, "-[M::%s::]\tutg%.6dl(%c)\tq::[%d, %d)\tt::[%d, %d)\n", __func__,
// (int32_t)(gchains->a[i+gchains->n].v>>1)+1, "+-"[gchains->a[i+gchains->n].v&1],
// gchains->a[i+gchains->n].qs, gchains->a[i+gchains->n].qe,
// gchains->a[i+gchains->n].rs, gchains->a[i+gchains->n].re);
// }
} }
} }
g_item->cnt = res->n; g_item->cnt = res->n;
@@ -8803,6 +8812,7 @@ static void worker_for_ul_rescall_alignment(void *data, long i, int tid) // call
// if(s->id+i!=3046/** && s->id+i!=3111**/) return; // if(s->id+i!=3046/** && s->id+i!=3111**/) return;
// if((s->id+i!=871) && (s->id+i!=963) && (s->id+i!=980)) return; // if((s->id+i!=871) && (s->id+i!=963) && (s->id+i!=980)) return;
// if(s->id+i!=963) return; // if(s->id+i!=963) return;
// if(s->id+i != 35437) return;
// fprintf(stderr, "\n[M::%s] rid::%ld, len::%lu, name::%.*s\n", __func__, s->id+i, s->len[i], // fprintf(stderr, "\n[M::%s] rid::%ld, len::%lu, name::%.*s\n", __func__, s->id+i, s->len[i],
// (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a); // (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a);
@@ -9603,7 +9613,7 @@ utg_rid_dt *get_r_ug_region(utg_rid_t *idx, uint64_t *n, uint64_t rid)
return (*n)?idx->p.a + idx->idx[rid]:NULL; return (*n)?idx->p.a + idx->idx[rid]:NULL;
} }
void rov2uov(uint64_t rid, const ul_idx_t *uref, utg_rid_dt *ru_map, uc_block_t *rovlp, ul_ov_t *res, uint32_t adjust_rev) uint64_t rov2uov(uint64_t rid, const ul_idx_t *uref, utg_rid_dt *ru_map, uc_block_t *rovlp, ul_ov_t *res, uint32_t adjust_rev, int64_t ulid)
{ {
uint64_t ori = ru_map->u&1, ts, te; uint64_t ori = ru_map->u&1, ts, te;
if(!ori) { if(!ori) {
@@ -9612,15 +9622,27 @@ void rov2uov(uint64_t rid, const ul_idx_t *uref, utg_rid_dt *ru_map, uc_block_t
ts = uref->r_ug->rg->seq[rid].len - rovlp->te; ts = uref->r_ug->rg->seq[rid].len - rovlp->te;
te = uref->r_ug->rg->seq[rid].len - rovlp->ts; te = uref->r_ug->rg->seq[rid].len - rovlp->ts;
} }
// if(ulid == 14714) {
// fprintf(stderr, "[M::%s::]\tori::%lu\trovlp->rev::%u\tro_t::[%u, %u)\tt::[%lu, %lu)\toff::%u\n", __func__,
// ori, rovlp->rev, rovlp->ts, rovlp->te, ts, te, ru_map->off);
// }
ts += ru_map->off; te += ru_map->off; ts += ru_map->off; te += ru_map->off;
memset(res, 0, sizeof(*res)); if(ts >= 0 && te <= uref->ug->g->seq[ru_map->u>>1].len) {
res->qn = 0; res->qs = rovlp->qs; res->qe = rovlp->qe; memset(res, 0, sizeof(*res));
res->tn = ru_map->u>>1; res->ts = ts; res->te = te; res->qn = 0; res->qs = rovlp->qs; res->qe = rovlp->qe;
res->el = rovlp->el; res->rev = (rovlp->rev == ori?0:1); res->tn = ru_map->u>>1; res->ts = ts; res->te = te;
if(adjust_rev && res->rev) {///for linear chaining res->el = rovlp->el; res->rev = (rovlp->rev == ori?0:1);
res->ts = uref->ug->g->seq[res->tn].len - te; if(adjust_rev && res->rev) {///for linear chaining
res->te = uref->ug->g->seq[res->tn].len - ts; res->ts = uref->ug->g->seq[res->tn].len - te;
} res->te = uref->ug->g->seq[res->tn].len - ts;
}
return 1;
}
return 0;
// if(ulid == 14714) {
// fprintf(stderr, "[M::%s::]\tulen::%u\trlen::%u\tro_t::[%u, %u)\tt::[%lu, %lu)\toff::%u\n", __func__,
// uref->ug->g->seq[res->tn].len, uref->r_ug->rg->seq[rid].len, rovlp->ts, rovlp->te, ts, te, ru_map->off);
// }
} }
void print_ul_ov_t(ul_ov_t *xs, const char* cmd) void print_ul_ov_t(ul_ov_t *xs, const char* cmd)
@@ -9629,9 +9651,9 @@ void print_ul_ov_t(ul_ov_t *xs, const char* cmd)
"+-"[xs->rev], (int)Get_NAME_LENGTH(R_INF, ((xs->tn<<1)>>1)), Get_NAME(R_INF, ((xs->tn<<1)>>1)), xs->ts, xs->te); "+-"[xs->rev], (int)Get_NAME_LENGTH(R_INF, ((xs->tn<<1)>>1)), Get_NAME(R_INF, ((xs->tn<<1)>>1)), xs->ts, xs->te);
} }
void gl_rg2ug_gen(ul_vec_t *r_cl, kv_ul_ov_t *u_cl, const ul_idx_t *uref, uint64_t is_el, uint64_t n_pchain) void gl_rg2ug_gen(ul_vec_t *r_cl, kv_ul_ov_t *u_cl, const ul_idx_t *uref, uint64_t is_el, uint64_t n_pchain, int64_t ulid)
{ {
uint64_t k, a_k, a_n; uc_block_t *z; utg_rid_dt *a; ul_ov_t *p; uint64_t k, a_k, a_n; uc_block_t *z; utg_rid_dt *a; ul_ov_t p;
u_cl->n = 0; u_cl->n = 0;
for (k = 0; k < r_cl->bb.n; k++) { for (k = 0; k < r_cl->bb.n; k++) {
z = &(r_cl->bb.a[k]); z = &(r_cl->bb.a[k]);
@@ -9641,14 +9663,22 @@ void gl_rg2ug_gen(ul_vec_t *r_cl, kv_ul_ov_t *u_cl, const ul_idx_t *uref, uint64
a = get_r_ug_region(uref->r_ug, &a_n, z->hid); a = get_r_ug_region(uref->r_ug, &a_n, z->hid);
if(!a) continue; if(!a) continue;
for (a_k = 0; a_k < a_n; a_k++) { for (a_k = 0; a_k < a_n; a_k++) {
kv_pushp(ul_ov_t, *u_cl, &p); // if(ulid == 14714) {
// fprintf(stderr, "\n+[M::%s::] %u\t%u\t%c\t%.*s(%u)\t%u\t%u\n", __func__, z->qs, z->qe, "+-"[z->rev], // fprintf(stderr, "\n+[M::%s::]\tq::[%u, %u)\t%c\t%.*s(%u)\tt::[%u, %u)\n", __func__, z->qs, z->qe, "+-"[z->rev],
// (int)Get_NAME_LENGTH(R_INF, z->hid), Get_NAME(R_INF, z->hid), (uint32_t)Get_READ_LENGTH(R_INF, z->hid), z->ts, z->te); // (int)Get_NAME_LENGTH(R_INF, z->hid), Get_NAME(R_INF, z->hid), (uint32_t)Get_READ_LENGTH(R_INF, z->hid),
// fprintf(stderr, "*[M::%s::] utg%.6d%c(%u)\t%c\t%u\n", __func__, // z->ts, z->te);
// (int32_t)(a[a_k].u>>1)+1, "lc"[uref->ug->u.a[a[a_k].u>>1].circ], uref->ug->u.a[a[a_k].u>>1].len, // fprintf(stderr, "*[M::%s::] utg%.6d%c(%u)\t%c\toff::%u\tpos::%u\n", __func__,
// "+-"[a[a_k].u&1], a[a_k].off); // (int32_t)(a[a_k].u>>1)+1, "lc"[uref->ug->u.a[a[a_k].u>>1].circ], uref->ug->u.a[a[a_k].u>>1].len,
rov2uov(z->hid, uref, &(a[a_k]), z, p, 1); // "+-"[a[a_k].u&1], a[a_k].off, a[a_k].pos);
p->el = 1; p->tn <<= 1; p->tn |= p->rev; p->qn = k/**uref->r_ug->idx[z->hid] + a_k**/;//for linear chain // }
if(!rov2uov(z->hid, uref, &(a[a_k]), z, &p, 1, ulid)) continue;
p.el = 1; p.tn <<= 1; p.tn |= p.rev; p.qn = k/**uref->r_ug->idx[z->hid] + a_k**/;//for linear chain
kv_push(ul_ov_t, *u_cl, p);
// if(ulid == 14714) {
// fprintf(stderr, "[M::%s::]\tutg%.6dl(%c)\tq::[%d, %d)\tt::[%d, %d)\ttlen::%u\n", __func__,
// (int32_t)(p->tn>>1)+1, "+-"[p->tn&1], p->qs, p->qe,
// p->ts, p->te, uref->ug->g->seq[p->tn>>1].len);
// }
// if(k == 2) { // if(k == 2) {
// fprintf(stderr, "[M::%s::] p->ts:%u, p->te:%u, z->ts:%u, z->te:%u, a[a_k].off:%u\n", __func__, p->ts, p->te, z->ts, z->te, a[a_k].off); // fprintf(stderr, "[M::%s::] p->ts:%u, p->te:%u, z->ts:%u, z->te:%u, a[a_k].off:%u\n", __func__, p->ts, p->te, z->ts, z->te, a[a_k].off);
// } // }
@@ -10142,7 +10172,7 @@ void extend_end_coord(mg_lchain_t *li, ul_ov_t *ui, const int64_t qlen, const in
} }
} }
void dump_linear_chain(ma_ug_t *ug, kv_ul_ov_t *lidx, vec_mg_lchain_t *res, int64_t qlen) void dump_linear_chain(ma_ug_t *ug, kv_ul_ov_t *lidx, vec_mg_lchain_t *res, int64_t qlen, int64_t ulid)
{ {
uint64_t i; int64_t iqs, iqe, its, ite; mg_lchain_t *p; uint64_t i; int64_t iqs, iqe, its, ite; mg_lchain_t *p;
kv_resize(mg_lchain_t, *res, lidx->n); kv_resize(mg_lchain_t, *res, lidx->n);
@@ -10154,12 +10184,21 @@ void dump_linear_chain(ma_ug_t *ug, kv_ul_ov_t *lidx, vec_mg_lchain_t *res, int6
p->off = i; p->score = lidx->a[i].sec; p->off = i; p->score = lidx->a[i].sec;
p->qs = lidx->a[i].qs; p->qe = lidx->a[i].qe; p->qs = lidx->a[i].qs; p->qe = lidx->a[i].qe;
p->rs = lidx->a[i].ts; p->re = lidx->a[i].te; p->rs = lidx->a[i].ts; p->re = lidx->a[i].te;
// if(ulid == 14714) {
// fprintf(stderr, "+[M::%s::]\tutg%.6dl(%c)\tq::[%d, %d)\tqlen::%ld\tt::[%d, %d)\ttlen::%u\n", __func__,
// (int32_t)(p->v>>1)+1, "+-"[p->v&1], p->qs, p->qe, qlen,
// p->rs, p->re, ug->g->seq[p->v>>1].len);
// }
extend_end_coord(p, NULL, qlen, ug->g->seq[p->v>>1].len, &iqs, &iqe, &its, &ite); extend_end_coord(p, NULL, qlen, ug->g->seq[p->v>>1].len, &iqs, &iqe, &its, &ite);
p->qs = iqs; p->qe = iqe; p->rs = its; p->re = ite; p->qs = iqs; p->qe = iqe; p->rs = its; p->re = ite;
// if(!ugl_cover_check(p->rs, p->re, &(ug->u.a[p->v>>1]))) res->n--; // if(!ugl_cover_check(p->rs, p->re, &(ug->u.a[p->v>>1]))) res->n--;
// fprintf(stderr, "chain_id:%d\t%u\t%u\t%c\tutg%.6dl(%u)\t%u\t%u\n", // fprintf(stderr, "chain_id:%d\t%u\t%u\t%c\tutg%.6dl(%u)\t%u\t%u\n",
// res->a[k].off, res->a[k].qs, res->a[k].qe, "+-"[res->a[k].v&1], (int32_t)(res->a[k].v>>1)+1, // res->a[k].off, res->a[k].qs, res->a[k].qe, "+-"[res->a[k].v&1], (int32_t)(res->a[k].v>>1)+1,
// g->seq[res->a[k].v>>1].len, res->a[k].rs, res->a[k].re); // g->seq[res->a[k].v>>1].len, res->a[k].rs, res->a[k].re);
// if(ulid == 14714) {
// fprintf(stderr, "-[M::%s::]\tutg%.6dl(%c)\tq::[%d, %d)\tt::[%d, %d)\n", __func__,
// (int32_t)(p->v>>1)+1, "+-"[p->v&1], p->qs, p->qe, p->rs, p->re);
// }
} }
} }
@@ -11670,7 +11709,7 @@ mg_lchain_t *uo, kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn)
} }
void gl_ug2rg_gen(const asg_t *rg, ul_vec_t *rch, ma_ug_t *ug, mg_lchain_t *uo, vec_mg_lchain_t *res, int64_t tOff, void gl_ug2rg_gen(const asg_t *rg, ul_vec_t *rch, ma_ug_t *ug, mg_lchain_t *uo, vec_mg_lchain_t *res, int64_t tOff,
kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn) int64_t ulid, kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn)
{ {
// fprintf(stderr, "\n[M::%s::] uo->qs:%d, uo->qe:%d\n", __func__, uo->qs, uo->qe); // fprintf(stderr, "\n[M::%s::] uo->qs:%d, uo->qe:%d\n", __func__, uo->qs, uo->qe);
///uo is a unitig alignment ///uo is a unitig alignment
@@ -11684,6 +11723,10 @@ kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn)
if(p.s >= re) break; if(p.s >= re) break;
assert(extract_rovlp_by_ug(&p, uo, res, tOff)); assert(extract_rovlp_by_ug(&p, uo, res, tOff));
// if(!extract_rovlp_by_ug(&p, uo, res, tOff)) {
// fprintf(stderr, "[M::%s::ulid::%ld]u_rs->%lu, u_re->%lu, p.s->%u, p.e->%u\n", __func__, ulid, rs, re, p.s, p.e);
// exit(1);
// }
res->a[res->n-1].score = uo->v; res->a[res->n-1].cnt = i; res->a[res->n-1].score = uo->v; res->a[res->n-1].cnt = i;
// fprintf(stderr, "[M::%s::]u_rs->%lu, u_re->%lu, p.s->%u, p.e->%u\n", __func__, rs, re, p.s, p.e); // fprintf(stderr, "[M::%s::]u_rs->%lu, u_re->%lu, p.s->%u, p.e->%u\n", __func__, rs, re, p.s, p.e);
@@ -11972,7 +12015,7 @@ int64_t flat_rovlp_chain(ma_ug_t *ug, mg_lchain_t *x, int64_t x_n, Chain_Data* d
} }
void gen_rovlp_chain_by_ul(const asg_t *rg, ul_vec_t *rch, const ul_idx_t *uref, kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn, mg_lchain_t *a, int64_t a_n, vec_mg_lchain_t *res, Chain_Data* dp, void gen_rovlp_chain_by_ul(const asg_t *rg, ul_vec_t *rch, const ul_idx_t *uref, kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn, mg_lchain_t *a, int64_t a_n, vec_mg_lchain_t *res, Chain_Data* dp,
int64_t dp_max_skip, int64_t dp_max_iter, int64_t dp_max_dis) int64_t dp_max_skip, int64_t dp_max_iter, int64_t dp_max_dis, int64_t ulid)
{ {
if(a_n == 0) return; if(a_n == 0) return;
int64_t k, l, res_n0 = res->n, tt = 0; ma_ug_t *ug = uref->ug; int64_t k, l, res_n0 = res->n, tt = 0; ma_ug_t *ug = uref->ug;
@@ -11981,7 +12024,7 @@ int64_t dp_max_skip, int64_t dp_max_iter, int64_t dp_max_dis)
for (k = 0, l = ug->g->seq[a[0].v>>1].len; k < a_n; k++) { for (k = 0, l = ug->g->seq[a[0].v>>1].len; k < a_n; k++) {
// fprintf(stderr, ">k->%ld, ls->%ld, le->%ld, rev->%c\n", k, l - ug->g->seq[a[k].v>>1].len, l, "+-"[a[k].v&1]); // fprintf(stderr, ">k->%ld, ls->%ld, le->%ld, rev->%c\n", k, l - ug->g->seq[a[k].v>>1].len, l, "+-"[a[k].v&1]);
l -= ug->g->seq[a[k].v>>1].len; l -= ug->g->seq[a[k].v>>1].len;
gl_ug2rg_gen(rg, rch, ug, &(a[k]), res, l, raw_idx, raw_chn); gl_ug2rg_gen(rg, rch, ug, &(a[k]), res, l, ulid, raw_idx, raw_chn);
l += ug->g->seq[a[k].v>>1].len + a[k].dist_pre; l += ug->g->seq[a[k].v>>1].len + a[k].dist_pre;
} }
mg_lchain_t *x = res->a + res_n0; int64_t x_n = res->n - res_n0, fn; mg_lchain_t *x = res->a + res_n0; int64_t x_n = res->n - res_n0, fn;
@@ -12296,10 +12339,12 @@ int64_t dp_max_skip, int64_t dp_max_iter, int64_t dp_max_dis)
int64_t k, ucn = uc->n; mg_lchain_t *ix; int64_t k, ucn = uc->n; mg_lchain_t *ix;
for (k = 0, swap->n = 0; k < ucn; k += ix->cnt + 1) { for (k = 0, swap->n = 0; k < ucn; k += ix->cnt + 1) {
ix = &(uc->a[k]); assert(ix->v == (uint32_t)-1); ix = &(uc->a[k]); assert(ix->v == (uint32_t)-1);
// fprintf(stderr, "\n[M::%s::ucn->%ld, k->%ld, kcnt->%d]\n", __func__, ucn, k, ix->cnt); // if(ulid == 14714) {
// print_debug_gchain(uref, uc->a + k + 1, ix->cnt, rch); // fprintf(stderr, "\n[M::%s::ucn->%ld, k->%ld, kcnt->%d]\n", __func__, ucn, k, ix->cnt);
// print_debug_gchain(uref, uc->a + k + 1, ix->cnt, rch);
// }
gen_rovlp_chain_by_ul(rg, rch, uref, raw_idx, raw_chn, uc->a + k + 1, ix->cnt, swap, dp, gen_rovlp_chain_by_ul(rg, rch, uref, raw_idx, raw_chn, uc->a + k + 1, ix->cnt, swap, dp,
dp_max_skip, dp_max_iter, dp_max_dis); dp_max_skip, dp_max_iter, dp_max_dis, ulid);
} }
///up to now, given a <x> in swap ///up to now, given a <x> in swap
@@ -12345,7 +12390,7 @@ int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid, Chain_Data* dp, c
// if(ulid != 86660) return 0; // if(ulid != 86660) return 0;
kv_ul_ov_t *idx = &(ll->lo), *init = &(ll->tk); int64_t max_idx; kv_ul_ov_t *idx = &(ll->lo), *init = &(ll->tk); int64_t max_idx;
idx->n = init->n = 0; idx->n = init->n = 0;
gl_rg2ug_gen(rch, idx, uref, 1, 2); gl_rg2ug_gen(rch, idx, uref, 1, 2, ulid);
if(idx->n == 0) return 0; if(idx->n == 0) return 0;
///generate linear chains ///generate linear chains
gen_linear_chains(idx, init, uref, uopt, bw, diff_ec_ul, rch->rlen, dp); gen_linear_chains(idx, init, uref, uopt, bw, diff_ec_ul, rch->rlen, dp);
@@ -12358,7 +12403,7 @@ int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid, Chain_Data* dp, c
// fprintf(stderr, "\n++[M::%s::%.*s(id:%ld), len:%u] idx->n:%lu\n", __func__, UL_INF.nid.a[ulid].n, UL_INF.nid.a[ulid].a, // fprintf(stderr, "\n++[M::%s::%.*s(id:%ld), len:%u] idx->n:%lu\n", __func__, UL_INF.nid.a[ulid].n, UL_INF.nid.a[ulid].a,
// ulid, rch->rlen, (uint64_t)idx->n); // ulid, rch->rlen, (uint64_t)idx->n);
dump_linear_chain(uref->ug, idx, &(gdp->l), rch->rlen); dump_linear_chain(uref->ug, idx, &(gdp->l), rch->rlen, ulid);
if(gdp->l.n == 0) return 0; if(gdp->l.n == 0) return 0;
// fprintf(stderr, "\n+++[M::%s::id->%ld, len->%u] idx->n:%lu\n", __func__, ulid, rch->rlen, (uint64_t)idx->n); // fprintf(stderr, "\n+++[M::%s::id->%ld, len->%u] idx->n:%lu\n", __func__, ulid, rch->rlen, (uint64_t)idx->n);
// kv_resize(uint64_t, ll->srt.a, idx->n); kv_resize(uint64_t, hap->snp_srt, idx->n); kv_resize(uint64_t, gdp->v, idx->n); // kv_resize(uint64_t, ll->srt.a, idx->n); kv_resize(uint64_t, hap->snp_srt, idx->n); kv_resize(uint64_t, gdp->v, idx->n);
@@ -14645,7 +14690,7 @@ ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_c
///for debug interval ///for debug interval
if(!load_all_ul_t(&UL_INF, gfa_name, &R_INF, ug)) { if(!load_all_ul_t(&UL_INF, gfa_name, &R_INF, ug)) {
gen_UL_reovlps(&sl, ug, sg, gfa_name, cutoff); gen_UL_reovlps(&sl, ug, sg, gfa_name, cutoff);
// exit(1); exit(1);
write_all_ul_t(&UL_INF, gfa_name, ug); write_all_ul_t(&UL_INF, gfa_name, ug);
} else if(double_check_cache){ } else if(double_check_cache){
if(drenew_UL_reovlps(&sl, ug, sg, gfa_name, cutoff)) { if(drenew_UL_reovlps(&sl, ug, sg, gfa_name, cutoff)) {
inter.h
+1 -1
View File
@@ -13,7 +13,7 @@
#define G_CHAIN_GAP 0.1 #define G_CHAIN_GAP 0.1
#define UG_SKIP 5 #define UG_SKIP 5
#define RG_SKIP 25 #define RG_SKIP 25
#define UG_SKIP_GRAPH_N 50 #define UG_SKIP_GRAPH_N 72
#define UG_SKIP_N 100 #define UG_SKIP_N 100
#define UG_ITER_N 5000 #define UG_ITER_N 5000
#define UG_DIS_N 50000 #define UG_DIS_N 50000