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

update graph chain

Browse Source
This commit is contained in:
chhylp123
2022-04-16 22:17:52 -04:00
parent 9c79bdfe7f
commit d862711623
1 changed files with 685 additions and 105 deletions
Download Patch File Download Diff File Expand all files Collapse all files
inter.cpp
+648 -68
View File
@@ -50,6 +50,10 @@ void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint6
#define MG_SHORT_K_EXT 10000 ///1000 in minigraph #define MG_SHORT_K_EXT 10000 ///1000 in minigraph
#define GC_OFFSET_RATE 0.0001
#define GC_OFFSET_POS 8
#define SEC_LEN_DIF 0.03
#define generic_key(x) (x) #define generic_key(x) (x)
KRADIX_SORT_INIT(gfa64, uint64_t, generic_key, 8) KRADIX_SORT_INIT(gfa64, uint64_t, generic_key, 8)
@@ -144,6 +148,12 @@ typedef struct { uint64_t x, y; } mg128_t;
KRADIX_SORT_INIT(128x, mg128_t, sort_key_128x, 8) KRADIX_SORT_INIT(128x, mg128_t, sort_key_128x, 8)
void radix_sort_128x(mg128_t *beg, mg128_t *end); void radix_sort_128x(mg128_t *beg, mg128_t *end);
#define mg_pathv_t_v_srt_key(x) ((x).v)
KRADIX_SORT_INIT(mg_pathv_t_v_srt, mg_pathv_t, mg_pathv_t_v_srt_key, member_size(mg_pathv_t, v))
#define mg_pathv_t_d_srt_key(x) ((x).d)
KRADIX_SORT_INIT(mg_pathv_t_d_srt, mg_pathv_t, mg_pathv_t_d_srt_key, member_size(mg_pathv_t, d))
typedef struct { typedef struct {
int32_t off, cnt; int32_t off, cnt;
uint32_t v; uint32_t v;
@@ -4709,7 +4719,7 @@ void gl_rg2ug_gen(ul_vec_t *r_cl, kv_ul_ov_t *u_cl, const ul_idx_t *uref, uint64
// (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, // (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,
// "+-"[a[a_k].u&1], a[a_k].off); // "+-"[a[a_k].u&1], a[a_k].off);
rov2uov(z->hid, uref, &(a[a_k]), z, p, 1); rov2uov(z->hid, uref, &(a[a_k]), z, p, 1);
p->tn <<= 1; p->tn |= p->rev; p->qn = uref->r_ug->idx[z->hid] + a_k;//for linear chain p->el = 1; p->tn <<= 1; p->tn |= p->rev; p->qn = uref->r_ug->idx[z->hid] + a_k;//for linear chain
// fprintf(stderr, "[M::%s::id->%ld] idx->n:%lu\n", __func__, ulid, (uint64_t)idx->n); // fprintf(stderr, "[M::%s::id->%ld] idx->n:%lu\n", __func__, ulid, (uint64_t)idx->n);
// fprintf(stderr, "-[M::%s::] %u\t%u\t%c\tutg%.6d%c(%u)\t%u\t%u\n", __func__, p->qs, p->qe, "+-"[p->rev], // fprintf(stderr, "-[M::%s::] %u\t%u\t%c\tutg%.6d%c(%u)\t%u\t%u\n", __func__, p->qs, p->qe, "+-"[p->rev],
// (int32_t)(p->tn>>1)+1, "lc"[uref->ug->u.a[p->tn>>1].circ], uref->ug->u.a[p->tn>>1].len, p->ts, p->te); // (int32_t)(p->tn>>1)+1, "lc"[uref->ug->u.a[p->tn>>1].circ], uref->ug->u.a[p->tn>>1].len, p->ts, p->te);
@@ -4812,19 +4822,20 @@ double diff_ec_ul, int64_t qlen, int64_t max_skip, uint64_t *idx, uint64_t *trac
if(sv[i].ts < ch[n_u].ts) ch[n_u].ts = sv[i].ts; if(sv[i].ts < ch[n_u].ts) ch[n_u].ts = sv[i].ts;
if(sv[i].qe > ch[n_u].qe) ch[n_u].qe = sv[i].qe; if(sv[i].qe > ch[n_u].qe) ch[n_u].qe = sv[i].qe;
if(sv[i].te > ch[n_u].te) ch[n_u].te = sv[i].te; if(sv[i].te > ch[n_u].te) ch[n_u].te = sv[i].te;
ch[n_u].qn = i;//start idx of read alignment in chain // ch[n_u].qn = i;//start idx of read alignment in chain
i = pop_pre(track[i]); i = pop_pre(track[i]);
} }
adjust_rev_tse(&(ch[n_u]), ug->g->seq[ch[n_u].tn].len, &its, &ite); adjust_rev_tse(&(ch[n_u]), ug->g->seq[ch[n_u].tn].len, &its, &ite);
ch[n_u].ts = its; ch[n_u].te = ite; ch[n_u].sec = (sc>0x3FFFFFFF?0x3FFFFFFF:sc); ch[n_u].ts = its; ch[n_u].te = ite; ch[n_u].sec = (sc>0x3FFFFFFF?0x3FFFFFFF:sc);
ch[n_u].qn += chain_offset; //start idx of read alignment in chain // ch[n_u].qn += chain_offset; //start idx of read alignment in chain
ch[n_u].tn = k + chain_offset; //end idx of read alignment in chain // ch[n_u].tn = k + chain_offset; //end idx of read alignment in chain
ch[n_u].qn = k + chain_offset; //end idx of read alignment in chain
n_u++; n_u++;
} }
for (i = 0; i < ch_n; ++i) { for (i = 0; i < ch_n; ++i) {
adjust_rev_tse(&(sv[i]), ug->g->seq[sv[i].tn].len, &its, &ite); adjust_rev_tse(&(sv[i]), ug->g->seq[sv[i].tn].len, &its, &ite);
sv[i].ts = its; sv[i].te = ite; sv[i].ts = its; sv[i].te = ite;
k = pop_pre(track[i]); sv[i].qn = k>=0?k+chain_offset:(uint32_t)-1; k = pop_pre(track[i]); sv[i].tn = k>=0?k+chain_offset:(uint32_t)-1;
} }
return n_u; return n_u;
} }
@@ -5077,7 +5088,8 @@ void dump_linear_chain(asg_t *g, kv_ul_ov_t *lidx, kv_ul_ov_t *autom, vec_mg_lch
res->n = 0; kv_resize(mg_lchain_t, *res, lidx->n); res->n = lidx->n; res->n = 0; kv_resize(mg_lchain_t, *res, lidx->n); res->n = lidx->n;
for (k = 0; k < lidx->n; k++) { for (k = 0; k < lidx->n; k++) {
memset(&(res->a[k]), 0, sizeof(res->a[k])); memset(&(res->a[k]), 0, sizeof(res->a[k]));
res->a[k].v = (autom->a[lidx->a[k].tn].tn<<1)|lidx->a[k].rev; // res->a[k].v = (autom->a[lidx->a[k].tn].tn<<1)|lidx->a[k].rev;
res->a[k].v = (lidx->a[k].tn<<1)|(lidx->a[k].rev);
///.off -> idx of original chain; cnt -> score of the chain ///.off -> idx of original chain; cnt -> score of the chain
res->a[k].off = k; res->a[k].score = lidx->a[k].sec; res->a[k].off = k; res->a[k].score = lidx->a[k].sec;
res->a[k].qs = lidx->a[k].qs; res->a[k].qe = lidx->a[k].qe; res->a[k].qs = lidx->a[k].qs; res->a[k].qe = lidx->a[k].qe;
@@ -5145,10 +5157,182 @@ inline int32_t cal_gchain_sc(const mg_path_dst_t *dj, const mg_lchain_t *li, con
return sc; return sc;
} }
void set_trans_arr(uint64_t *trans, vec_sp_node_t *out, int64_t idx)
{
int64_t i;
for (i = idx; i >=0; ) {
// if(out->a[idx]->v == 65) {
// fprintf(stderr, "+[M::%s::] out->a[%ld]->v:%u\n", __func__, i, out->a[i]->v);
// }
trans[i]++;
i = out->a[i]->pre;
}
}
int64_t select_mul_way_nodes(mg_pathv_t *a, int64_t a_n, vec_sp_node_t *out, float len_dif, int32_t m_pathn, uint64_t *flag)
{
int64_t i, k, pd, kd, occ, tt = 0; uint32_t pp; mg_pathv_t *p = NULL;
if(a_n > 1) radix_sort_mg_pathv_t_d_srt(a, a + a_n);
for (i = 0; i < a_n; i++) {
p = &(a[i]); pp = (p->d<<1)>>1; pd = out->a[p->pre]->di>>32; occ = 1;
if(p->d&0x80000000) break;
for (k = 0; k < a_n; k++) {
if(k == i) continue; ///same alignment
if(((a[k].d<<1)>>1) == pp) continue; //same path
if(a[k].v!=p->v) continue;
kd = out->a[a[k].pre]->di>>32;
if(kd >= (pd*(1-len_dif)) && kd <= (pd*(1+len_dif))) occ++;
}
if(occ >= m_pathn) {
flag[p->pre] = 1; tt++;
}
}
return tt;
}
int32_t phase_mul_ways(vec_mg_pathv_t *res, st_mt_t *dst_done, vec_sp_node_t *out, int32_t n_dst, mg_path_dst_t *dst, int32_t max_k, float len_dif)
{
int64_t i, j, z, zl, n = 0, n_mpath, kk_p, od, res_n = res->n, pid; mg_pathv_t *h;
dst_done->n = 0; kv_resize(uint64_t, *dst_done, out->n);
uint64_t *trans = dst_done->a; memset(dst_done->a, 0, out->n*sizeof(*(dst_done->a)));
n_mpath = 0;
for (i = 0; i < n_dst; ++i) { // mark dst vertices with a target distance
mg_path_dst_t *t = &dst[i];
if (t->n_path > 0 && t->target_dist >= 0 && t->path_end >= 0){
assert((int32_t)(out->a[t->path_end]->di>>32) == t->target_dist);
if(t->n_path >= max_k) {
t->n_path = 0;
for (z = zl = t->path_end; z >= 0;) {
zl = z;
z = out->a[z]->pre;
}
n += 2; trans[t->path_end] = trans[zl] = 1;
} else {
kk_p = 0;
for (j = t->path_end; j < (int32_t)out->n; j++) {
od = out->a[j]->di>>32;
if(od >= (t->target_dist*(1-len_dif)) &&
od <= (t->target_dist*(1+len_dif))) {
if(out->a[j]->v == t->v) kk_p++;
} else {
break;
}
}
for (j = t->path_end-1; j >=0; j--) {
od = out->a[j]->di>>32;
if(od >= (t->target_dist*(1-len_dif)) &&
od <= (t->target_dist*(1+len_dif))) {
if(out->a[j]->v == t->v) kk_p++;
} else {
break;
}
}
assert(kk_p > 0 && kk_p <= t->n_path);
n_mpath += kk_p;
}
}
}
// if(detect_mul_way && src == 74) {
// fprintf(stderr, "+[M::%s::] src:%u, dst:%u, n_mpath:%d\n", __func__, src, dst[0].v, n_mpath);
// }
if(n_mpath > 1) {
for (i = 0, pid = 0; i < n_dst; ++i) { // mark dst vertices with a target distance
mg_path_dst_t *t = &dst[i];
if (t->n_path > 0 && t->target_dist >= 0 && t->path_end >= 0){
assert((int32_t)(out->a[t->path_end]->di>>32) == t->target_dist);
assert(t->n_path < max_k);
for (j = t->path_end; j < (int32_t)out->n; j++) {
od = out->a[j]->di>>32;
if(od >= (t->target_dist*(1-len_dif)) &&
od <= (t->target_dist*(1+len_dif))) {
if(out->a[j]->v == t->v) {
for (z = j; z >= 0;) {
kv_pushp(mg_pathv_t, *res, &h);
h->v = out->a[z]->v; h->pre = z;
h->d = pid; if(j!=t->path_end) h->d |= 0x80000000;
z = out->a[z]->pre;
}
pid++;
}
} else {
break;
}
}
for (j = t->path_end-1; j >=0; j--) {
od = out->a[j]->di>>32;
if(od >= (t->target_dist*(1-len_dif)) &&
od <= (t->target_dist*(1+len_dif))) {
if(out->a[j]->v == t->v) {
for (z = j; z >= 0;) {
kv_pushp(mg_pathv_t, *res, &h);
h->v = out->a[z]->v; h->pre = z;
h->d = pid; if(j!=t->path_end) h->d |= 0x80000000;
z = out->a[z]->pre;
}
pid++;
}
} else {
break;
}
}
}
}
radix_sort_mg_pathv_t_v_srt(res->a + res_n, res->a + res->n);
for (i = res_n+1, j = res_n/**, n = 0**/; i <= (int64_t)res->n; ++i) {
if (i == (int64_t)res->n || res->a[i].v != res->a[j].v) {
n += select_mul_way_nodes(res->a + j, i - j, out, len_dif, n_mpath, trans);
j = i;
}
}
res->n = res_n;
}
if(n > 0) {//found some nodes
for (i = n = 0; (uint32_t)i < out->n; ++i) { // generate coordinate translations
if (trans[i]) {
trans[i] = n++;
} else {
trans[i] = (uint32_t)-1;
}
}
kv_resize(mg_pathv_t, *res, res->n + n); //res->n += n;
for (i = 0; (uint32_t)i < out->n; ++i) { // generate the backtrack array
mg_pathv_t *p;
if (trans[i] == (uint32_t)-1) continue;
p = &res->a[trans[i]+res->n];
p->v = out->a[i]->v, p->d = out->a[i]->di >> 32;
if(out->a[i]->pre < 0) {
p->pre = out->a[i]->pre;
} else {
if(trans[out->a[i]->pre] == (uint32_t)-1) p->pre = -2;
else p->pre = trans[out->a[i]->pre];
}
}
res->n += n;
for (i = 0; i < n_dst; ++i) // translate "path_end"
if (dst[i].path_end >= 0)
dst[i].path_end = trans[dst[i].path_end];
}
return n_mpath;
}
///max_dist is like the overlap length in string graph ///max_dist is like the overlap length in string graph
///first_src_ban do not allow co-linear chain at the same node ///first_src_ban do not allow co-linear chain at the same node
void hc_shortest_k(void *km0, const asg_t *g, uint32_t src, int32_t n_dst, mg_path_dst_t *dst, int32_t max_dist, int32_t max_k, void hc_shortest_k(void *km0, const asg_t *g, uint32_t src, int32_t n_dst, mg_path_dst_t *dst, int32_t max_dist, int32_t max_k,
st_mt_t *dst_done, uint64_t *dst_group, vec_sp_node_t *out, vec_mg_pathv_t *res, uint64_t first_src_ban) st_mt_t *dst_done, uint64_t *dst_group, vec_sp_node_t *out, vec_mg_pathv_t *res, uint64_t first_src_ban,
uint64_t detect_mul_way, float len_dif)
{ {
sp_node_t *p, *root = 0; sp_node_t *p, *root = 0;
sp_topk_t *q; sp_topk_t *q;
@@ -5298,14 +5482,18 @@ st_mt_t *dst_done, uint64_t *dst_group, vec_sp_node_t *out, vec_mg_pathv_t *res,
if (t->target_dist >= 0) { if (t->target_dist >= 0) {
///src is from li from li to lj, so the dis is generally increased; dijkstra algorithm ///src is from li from li to lj, so the dis is generally increased; dijkstra algorithm
///target_dist should be the distance on query ///target_dist should be the distance on query
if (dist == t->target_dist && t->check_hash && r->hash == t->target_hash) done = 1; if (dist == t->target_dist && t->check_hash && r->hash == t->target_hash) {
else if ((dist > t->target_dist + MG_SHORT_K_EXT) && (dist > (t->target_dist>>4))) done = 1; done = 1;
} else if ((dist > t->target_dist + MG_SHORT_K_EXT) && (dist > (t->target_dist>>4)) &&
(dist > (t->target_dist*1.25))) {
done = 1;
}
} }
} }
++t->n_path;///we found a path to the alignment t ++t->n_path;///we found a path to the alignment t
if (t->n_path >= max_k) done = 1; if (t->n_path >= max_k) done = 1;
} }
if (dst_done->a[off + j] == 0 && done) if (detect_mul_way == 0 && dst_done->a[off + j] == 0 && done)
dst_done->a[off + j] = 1, ++n_done; dst_done->a[off + j] = 1, ++n_done;
} }
///if all alignments have been settle down ///if all alignments have been settle down
@@ -5376,8 +5564,14 @@ st_mt_t *dst_done, uint64_t *dst_group, vec_sp_node_t *out, vec_mg_pathv_t *res,
if (dst[i].n_path > 0) ++n_found;///n_path might be larger than 16 if (dst[i].n_path > 0) ++n_found;///n_path might be larger than 16
///we can assume n_pathv = NULL for now ///we can assume n_pathv = NULL for now
if (n_found > 0 && res) { // then generate the backtrack array if (n_found > 0 && res) { // then generate the backtrack array
int32_t n; dst_done->n = 0; kv_resize(uint64_t, *dst_done, out->n); int32_t n, n_mpath = 1; dst_done->n = 0; kv_resize(uint64_t, *dst_done, out->n);
uint64_t *trans = dst_done->a; memset(dst_done->a, 0, out->n*sizeof(*(dst_done->a))); uint64_t *trans = dst_done->a; memset(dst_done->a, 0, out->n*sizeof(*(dst_done->a)));
if(detect_mul_way) {
n_mpath = phase_mul_ways(res, dst_done, out, n_dst, dst, max_k, len_dif);
}
if(n_mpath == 1) {
// KCALLOC(km, trans, n_out); // used to squeeze unused elements in out[] // KCALLOC(km, trans, n_out); // used to squeeze unused elements in out[]
///n_out: how many times that nodes in graph have been visited ///n_out: how many times that nodes in graph have been visited
///note one node might be visited multiples times ///note one node might be visited multiples times
@@ -5387,15 +5581,15 @@ st_mt_t *dst_done, uint64_t *dst_group, vec_sp_node_t *out, vec_mg_pathv_t *res,
if (t->n_path > 0 && t->target_dist >= 0 && t->path_end >= 0) if (t->n_path > 0 && t->target_dist >= 0 && t->path_end >= 0)
trans[(uint32_t)out->a[t->path_end]->di] = 1;///(int32_t)out[]->di: traverse track corresponds to the alignment chain dst[] trans[(uint32_t)out->a[t->path_end]->di] = 1;///(int32_t)out[]->di: traverse track corresponds to the alignment chain dst[]
} }
for (i = 0; (uint32_t)i < out->n; ++i) { // mark dst vertices without a target distance // for (i = 0; (uint32_t)i < out->n; ++i) { // mark dst vertices without a target distance
k = kh_get(sp2, h2, out->a[i]->v); // k = kh_get(sp2, h2, out->a[i]->v);
if (k != kh_end(h2)) { // TODO: check if this is correct! // if (k != kh_end(h2)) { // TODO: check if this is correct!
int32_t off = kh_val(h2, k)>>32, cnt = (int32_t)kh_val(h2, k); // int32_t off = kh_val(h2, k)>>32, cnt = (int32_t)kh_val(h2, k);
for (j = off; j < off + cnt; ++j) // for (j = off; j < off + cnt; ++j)
if (dst[j].target_dist < 0) // if (dst[j].target_dist < 0)
trans[i] = 1; // trans[i] = 1;
} // }
} // }
for (i = (int32_t)(out->n) - 1; i >= 0; --i) // mark all predecessors for (i = (int32_t)(out->n) - 1; i >= 0; --i) // mark all predecessors
if (trans[i] && out->a[i]->pre >= 0) if (trans[i] && out->a[i]->pre >= 0)
trans[out->a[i]->pre] = 1; trans[out->a[i]->pre] = 1;
@@ -5416,6 +5610,7 @@ st_mt_t *dst_done, uint64_t *dst_group, vec_sp_node_t *out, vec_mg_pathv_t *res,
if (dst[i].path_end >= 0) if (dst[i].path_end >= 0)
dst[i].path_end = trans[dst[i].path_end]; dst[i].path_end = trans[dst[i].path_end];
} }
}
km_destroy(km); km_destroy(km);
} }
@@ -5457,13 +5652,57 @@ int64_t *n_u_, int64_t *n_v_)
return n_u; return n_u;
} }
int64_t hc_gchain1_dp(void *km, const ma_ug_t *ug, vec_mg_lchain_t *lc, vec_mg_lchain_t *sw, vec_mg_path_dst_t *dst, vec_sp_node_t *out, vec_mg_pathv_t *path, void set_ul_ov_t_by_mg_lchain_t(ul_ov_t *u, mg_lchain_t *l)
int64_t qlen, const ug_opt_t *uopt, int64_t bw, double diff_thre, uint64_t *srt, st_mt_t *bf, int64_t *f, uint64_t *p, uint64_t *v) {
u->tn = l->v>>1; u->rev = (l->v&1);
u->ts = l->rs; u->te = l->re;
u->qs = l->qs; u->qe = l->qe;
}
uint64_t primary_chain_check(uint64_t *idx, int64_t idx_n, mg_lchain_t *a)
{
if(idx_n <= 0) return 0;
ul_ov_t m; memset(&m, 0, sizeof(m)); int64_t m_sc = -1, i, a_n; uint64_t s_idx, e_idx, ovlp, novlp;
for (i = a_n = 0; i < idx_n; ++i) {
if(((int64_t)(idx[i]>>32)) > m_sc) {
m_sc = ((int64_t)(idx[i]>>32)); m.qn = i;
m.ts = a_n; m.te = a_n + ((uint32_t)idx[i]);
m.qs = a[m.ts].qs; m.qe = a[m.te-1].qe;
}
a_n += ((uint32_t)idx[i]);
}
assert(a[m.ts].qs<=a[m.te-1].qs && a[m.te-1].qe>=a[m.ts].qe);
for (i = a_n = 0; i < idx_n; ++i) {
s_idx = a[a_n].qs; a_n += ((uint32_t)idx[i]); e_idx = a[a_n-1].qe;
if(i == m.qn) continue;
ovlp = ((MIN(m.qe, e_idx) > MAX(m.qs, s_idx))? (MIN(m.qe, e_idx) - MAX(m.qs, s_idx)):0);
novlp = (e_idx - s_idx) - ovlp;
if(novlp > ((m.qe-m.qs)*GC_OFFSET_RATE) && novlp > GC_OFFSET_POS) break;
}
if(i >= idx_n) return 1;
return 0;
}
int64_t g_adjacent_dis(const asg_t *g, uint32_t v, uint32_t w)
{
uint32_t nv, i; asg_arc_t *av = NULL;
nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
for (i = 0; i < nv; i++) {
if(av[i].del || av[i].v != w) continue;
return (uint32_t)av[i].ul;
}
return -1;
}
int64_t hc_gchain1_dp(void *km, const ul_idx_t *uref, const ma_ug_t *ug, vec_mg_lchain_t *lc, vec_mg_lchain_t *sw, vec_mg_path_dst_t *dst, vec_sp_node_t *out, vec_mg_pathv_t *path,
int64_t qlen, const ug_opt_t *uopt, int64_t bw, double diff_thre, double ng_diff_thre, uint64_t *srt, st_mt_t *bf, int64_t *f, uint64_t *p, uint64_t *v)
{ {
bf->n = 0; bf->n = 0;
if(lc->n == 0) return 0; if(lc->n == 0) return 0;
int64_t i, j, lc_n = lc->n, n_ext, mm_ovlp, target_dist, max_target_dist, x, m_idx, m_sc; int64_t i, j, lc_n = lc->n, n_ext, mm_ovlp, target_dist, max_target_dist, x, m_idx, m_sc, qo, sc;
mg_lchain_t *r, *li, *lj; mg_path_dst_t *q; asg_t *g = ug->g; uint64_t isolated, *u; int64_t max_f, max_j = -1, max_d = -1, max_inner = 0, share; uint32_t max_hash = 0; int64_t k, k0, n_u, n_v, ni;
mg_lchain_t *r, *li, *lj; mg_path_dst_t *q; asg_t *g = ug->g; uint64_t isolated, *u; ul_ov_t ui, uj;
for (i = n_ext = 0; i < lc_n; i++) { for (i = n_ext = 0; i < lc_n; i++) {
r = &lc->a[i]; r->dist_pre = -1; isolated = 0;///dist_pre -> parent in graph chain r = &lc->a[i]; r->dist_pre = -1; isolated = 0;///dist_pre -> parent in graph chain
if((r->re < g->seq[r->v>>1].len) && (r->rs > 0)) isolated = 1;///UL contained in one vertice if((r->re < g->seq[r->v>>1].len) && (r->rs > 0)) isolated = 1;///UL contained in one vertice
@@ -5488,6 +5727,67 @@ int64_t qlen, const ug_opt_t *uopt, int64_t bw, double diff_thre, uint64_t *srt,
memcpy(lc->a, sw->a, lc_n *sizeof((*(lc->a)))); memcpy(lc->a, sw->a, lc_n *sizeof((*(lc->a))));
// fprintf(stderr, "[M::%s::] n_ext:%ld, lc_n:%ld\n", __func__, n_ext, lc_n); // fprintf(stderr, "[M::%s::] n_ext:%ld, lc_n:%ld\n", __func__, n_ext, lc_n);
//first non-gap chain
for (i = 0; i < n_ext; ++i) { // core loop
li = &lc->a[i]; set_ul_ov_t_by_mg_lchain_t(&ui, li);
mm_ovlp = max_ovlp(g, li->v^1);
x = (li->qs + mm_ovlp)*diff_thre; if(x < bw) x = bw;
x += li->qs + mm_ovlp;
if (x > qlen+1) x = qlen+1;
x = find_mg_lchain_max(i, lc->a, x+G_CHAIN_INDEL);
// fprintf(stderr, "\nli->(%ld)\tutg%.6d%c(%u)\tqs:%u\tqe:%u\t%c\trs:%u\tre:%u\tsrc:%u\tscore:%d, x:%ld\n",
// i, (int32_t)(li->v>>1)+1, "lc"[ug->u.a[li->v>>1].circ], ug->u.a[li->v>>1].len,
// li->qs, li->qe, "+-"[li->v&1], li->rs, li->re, li->v^1, li->score, x);
max_f = li->score, max_j = -1;
// collect potential destination vertices
for (j = x; j >= 0; --j) {
lj = &lc->a[j]; ///extend_end_coord(lj, qlen, g->seq[lj->v>>1].len, &jqs, &jqe, &jrs, &jre);
//even this pair has a overlap, its length will be very small; just ignore; only for non-gapped chains
if(lj->qe+G_CHAIN_INDEL <= li->qs) break;
if(lj->qs >= li->qs+G_CHAIN_INDEL) continue;
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)
if(li->v!=lj->v && get_ecov_adv(uref, uopt, li->v^1, lj->v^1, bw, ng_diff_thre, qo, 0, &share)) {
sc = li->score + f[j];
if(sc > max_f) {
max_f = sc; max_j = j;
}
}
}
f[i] = max_f, p[i] = max_j<0?(uint64_t)-1:max_j;
li->dist_pre = max_j<0?-1:g_adjacent_dis(g, li->v^1, lc->a[max_j].v^1); li->inner_pre = 0;
li->hash_pre = max_j<0?0:(__ac_Wang_hash((li->v^1))+__ac_Wang_hash((lc->a[max_j].v^1)));
// fprintf(stderr, "i->%ld, utg%.6d%c->utg%.6d%c, max_f:%ld\n", i, (int32_t)(li->v>>1)+1, "lc"[ug->u.a[li->v>>1].circ],
// max_j<0?0:(int32_t)(lc->a[max_j].v>>1)+1, max_j<0?'*':"lc"[ug->u.a[lc->a[max_j].v>>1].circ], max_f);
}
kv_resize(uint64_t, *bf, (uint64_t)lc_n); u = bf->a;
hc_chain_backtrack(n_ext, f, p, srt, u, v, &n_u, &n_v);
for (i = 0; i < lc_n - n_ext; ++i) {
u[n_u++] = (((uint64_t)lc->a[n_ext + i].score)<<32) | 1;
v[n_v++] = n_ext + i;
}
sw->n = 0; kv_resize(mg_lchain_t, *sw, (uint64_t)n_v); m_idx = m_sc = -1; bf->n = 0;
for (i = 0, k = 0; i < n_u; ++i) {
k0 = k, ni = (int32_t)u[i];
for (j = 0; j < ni; ++j) {
sw->a[k++] = lc->a[v[k0 + (ni - j - 1)]];
}
if(m_idx < 0 || m_sc < ((int64_t)(u[i]>>32))) {
m_idx = i; m_sc = ((int64_t)(u[i]>>32));
}
}
assert(k == n_v); bf->n = n_u;
if(primary_chain_check(u, n_u, sw->a)) {
memcpy(lc->a, sw->a, n_v*sizeof(mg_lchain_t));
return m_idx;
}
bf->n = 0;
///then gapped-chaining
for (i = 0; i < n_ext; ++i) { // core loop for (i = 0; i < n_ext; ++i) { // core loop
li = &lc->a[i]; li = &lc->a[i];
mm_ovlp = max_ovlp(g, li->v^1); mm_ovlp = max_ovlp(g, li->v^1);
@@ -5533,15 +5833,14 @@ int64_t qlen, const ug_opt_t *uopt, int64_t bw, double diff_thre, uint64_t *srt,
} }
// confirm reach-ability // confirm reach-ability
int64_t max_f = li->score, max_j = -1, max_d = -1, max_inner = 0; uint32_t max_hash = 0; max_f = li->score, max_j = -1, max_d = -1, max_inner = 0; max_hash = 0;
if(dst->n) { if(dst->n) {
max_target_dist *= (1+diff_thre); if(max_target_dist < bw) max_target_dist = bw; max_target_dist *= (1+diff_thre); if(max_target_dist < bw) max_target_dist = bw;
hc_shortest_k(km, g, li->v^1, dst->n, dst->a, max_target_dist, MG_MAX_SHORT_K, bf, srt, out, NULL, 1); hc_shortest_k(km, g, li->v^1, dst->n, dst->a, max_target_dist, MG_MAX_SHORT_K, bf, srt, out, NULL, 1, 0, 0);
// remove unreachable destinations // remove unreachable destinations
//TODO: check sequence identity //TODO: check sequence identity
for (j = 0; j < (int64_t)dst->n; ++j) { for (j = 0; j < (int64_t)dst->n; ++j) {
mg_path_dst_t *dj = &dst->a[j]; mg_path_dst_t *dj = &dst->a[j];
int32_t sc;
if (dj->n_path == 0) continue; // unreachable if (dj->n_path == 0) continue; // unreachable
sc = cal_gchain_sc(dj, li, lc->a, f, bw, diff_thre, W_CHN_PEN_GAP); sc = cal_gchain_sc(dj, li, lc->a, f, bw, diff_thre, W_CHN_PEN_GAP);
@@ -5566,7 +5865,6 @@ int64_t qlen, const ug_opt_t *uopt, int64_t bw, double diff_thre, uint64_t *srt,
// max_j<0?0:(int32_t)(lc->a[max_j].v>>1)+1, max_j<0?'*':"lc"[ug->u.a[lc->a[max_j].v>>1].circ], max_f); // max_j<0?0:(int32_t)(lc->a[max_j].v>>1)+1, max_j<0?'*':"lc"[ug->u.a[lc->a[max_j].v>>1].circ], max_f);
} }
int64_t k, k0, n_u, n_v, ni;
kv_resize(uint64_t, *bf, (uint64_t)lc_n); u = bf->a; kv_resize(uint64_t, *bf, (uint64_t)lc_n); u = bf->a;
hc_chain_backtrack(n_ext, f, p, srt, u, v, &n_u, &n_v); hc_chain_backtrack(n_ext, f, p, srt, u, v, &n_u, &n_v);
for (i = 0; i < lc_n - n_ext; ++i) { for (i = 0; i < lc_n - n_ext; ++i) {
@@ -5574,7 +5872,7 @@ int64_t qlen, const ug_opt_t *uopt, int64_t bw, double diff_thre, uint64_t *srt,
v[n_v++] = n_ext + i; v[n_v++] = n_ext + i;
} }
sw->n = 0; kv_resize(mg_lchain_t, *sw, (uint64_t)n_v); m_idx = m_sc = -1; sw->n = 0; kv_resize(mg_lchain_t, *sw, (uint64_t)n_v); m_idx = m_sc = -1; bf->n = 0;
for (i = 0, k = 0; i < n_u; ++i) { for (i = 0, k = 0; i < n_u; ++i) {
k0 = k, ni = (int32_t)u[i]; k0 = k, ni = (int32_t)u[i];
for (j = 0; j < ni; ++j) { for (j = 0; j < ni; ++j) {
@@ -5590,18 +5888,32 @@ int64_t qlen, const ug_opt_t *uopt, int64_t bw, double diff_thre, uint64_t *srt,
} }
void debug_gchain(const asg_t *g, mg_lchain_t *a, uint64_t n) void debug_gchain(void *km, const asg_t *g, mg_lchain_t *a, uint64_t n, st_mt_t *dst_done, vec_sp_node_t *out)
{ {
uint64_t k, i, v, w, nv; asg_arc_t *av; uint64_t k, i, v, w, nv; int64_t dd; asg_arc_t *av; mg_path_dst_t dst; uint64_t dst_group;
for (k = 1; k < n; k++) { for (k = 1; k < n; k++) {
v = a[k-1].v; w = a[k].v; v = a[k].v^1; w = a[k-1].v^1;
nv = asg_arc_n(g, v); av = asg_arc_a(g, v); nv = asg_arc_n(g, v); av = asg_arc_a(g, v);
for (i = 0; i < nv; i++) { for (i = 0; i < nv; i++) {
if(av[i].v == w) break; if(av[i].v == w) break;
} }
if(i >= nv) { // if(i >= nv) {
// fprintf(stderr, "[M::%s::]\n", __func__); // // fprintf(stderr, "[M::%s::]\n", __func__);
fprintf(stderr, "[M::%s::]\tutg%.6dl(%c)\t->\tutg%.6dl(%c)\n", __func__, (int32_t)(v>>1)+1, "+-"[v&1], (int32_t)(w>>1)+1, "+-"[w&1]); // fprintf(stderr, "[M::%s::]\tutg%.6dl(%c)\t->\tutg%.6dl(%c)\n", __func__, (int32_t)(v>>1)+1, "+-"[v&1], (int32_t)(w>>1)+1, "+-"[w&1]);
// }
if(i < nv) {
dd = (int64_t)((uint32_t)(av[i].ul));
} else {
memset(&dst, 0, sizeof(dst));
dst.v = w;
dst.target_dist = a[k-1].dist_pre;
dst.target_hash = 0; dst.check_hash = 0;
hc_shortest_k(km, g, v, 1, &dst, dst.target_dist, MG_MAX_SHORT_K, dst_done, &dst_group, out, NULL, 1, 0, 0);
dd = dst.dist;
}
if(a[k-1].dist_pre != dd) {
fprintf(stderr, "[M::%s::]\tutg%.6dl(%c)\t->\tutg%.6dl(%c)\tdist_pre:%d\td:%ld\n", __func__, (int32_t)(v>>1)+1, "+-"[v&1],
(int32_t)(w>>1)+1, "+-"[w&1], a[k-1].dist_pre, dd);
} }
} }
} }
@@ -5625,37 +5937,70 @@ void debug_gchain2(const asg_t *g, mg_pathv_t *a, uint64_t n)
void reverse_track(mg_pathv_t *a, uint64_t a_n) void reverse_track(mg_pathv_t *a, uint64_t a_n)
{ {
uint64_t k; mg_pathv_t z; int64_t k, hn = (a_n>>1); mg_pathv_t z;
for (k = 0; k < (a_n>>1); k++) { for (k = a_n-1; k >= 1; k--) a[k].d -= a[k-1].d;
for (k = 0; k < hn; k++) {
z = a[k]; a[k] = a[a_n - k - 1]; a[a_n - k - 1] = z; z = a[k]; a[k] = a[a_n - k - 1]; a[a_n - k - 1] = z;
a[k].v ^= 1; a[a_n - k - 1].v ^= 1; a[k].v ^= 1; a[a_n - k - 1].v ^= 1;
} }
if(a_n&1) a[k].v ^= 1; if(a_n&1) a[k].v ^= 1;
} }
void dbg_print(mg_pathv_t *a, int64_t a_n)
{
int64_t k;
for (k = 0; k < a_n; k++) {
if(a[k].pre == -2) break;
}
if(k < a_n) {
fprintf(stderr, "+[M::%s::] src:utg%.6dl(v:%u), dst:utg%.6dl(v:%u)\n", __func__,
(int32_t)(a[0].v>>1)+1, a[0].v, (int32_t)(a[a_n-1].v>>1)+1, a[a_n-1].v);
for (k = 0; k < a_n; k++) {
fprintf(stderr, "-[M::%s::] utg%.6dl(v:%u), pre:%d, d:%u\n", __func__, (int32_t)(a[k].v>>1)+1, a[k].v,
a[k].pre, a[k].d);
}
}
}
uint32_t gen_gchain_track(void *km, mg_lchain_t *a, int64_t a_n, const asg_t *g, uint32_t gen_gchain_track(void *km, mg_lchain_t *a, int64_t a_n, const asg_t *g,
st_mt_t *dst_done, vec_sp_node_t *out, vec_mg_pathv_t *res) st_mt_t *dst_done, vec_sp_node_t *out, vec_mg_pathv_t *res)
{ {
int64_t k, p_n; mg_lchain_t *l0, *l1; mg_path_dst_t dst; uint64_t dst_group; mg_pathv_t *p; int64_t k, p_n/**, trav_occ = 0**/; mg_lchain_t *l0, *l1; mg_path_dst_t dst; uint64_t dst_group; mg_pathv_t *p;
res->n = 0; kv_pushp(mg_pathv_t, *res, &p); p->v = p->d = (uint32_t)-1; p->d = 0; res->n = 0; kv_pushp(mg_pathv_t, *res, &p); p->v = (uint32_t)-1; p->d = 0; p->pre = 0;
for (k = 1; k < a_n; k++) { for (k = 1; k < a_n; k++) {
l0 = a + k - 1; l1 = a + k; l0 = a + k - 1; l1 = a + k;
assert(!l1->inner_pre); assert(!l1->inner_pre); assert(l1->dist_pre >= 0);
memset(&dst, 0, sizeof(dst)); memset(&dst, 0, sizeof(dst));
dst.v = l0->v^1; dst.v = l0->v^1;
assert(l1->dist_pre >= 0); assert(l1->dist_pre >= 0);
dst.target_dist = l1->dist_pre; dst.target_dist = l1->dist_pre;
dst.target_hash = l1->hash_pre; dst.target_hash = l1->hash_pre;
dst.check_hash = 1; p_n = res->n; dst.check_hash = 1; p_n = res->n;
hc_shortest_k(km, g, l1->v^1, 1, &dst, dst.target_dist, MG_MAX_SHORT_K, dst_done, &dst_group, out, res, 1); if((dst.target_hash != (__ac_Wang_hash((l1->v^1))+__ac_Wang_hash(dst.v))) ||
(g_adjacent_dis(g, l1->v^1, dst.v) != dst.target_dist)) {
hc_shortest_k(km, g, l1->v^1, 1, &dst, dst.target_dist*(1+SEC_LEN_DIF), MG_MAX_SHORT_K, dst_done, &dst_group, out, res, 1, 1, SEC_LEN_DIF);
// debug_gchain2(g, res->a + p_n, res->n - p_n); // debug_gchain2(g, res->a + p_n, res->n - p_n);
// fprintf(stderr, "[M::%s::n->%ld]\tutg%.6dl(%c)\t->\tutg%.6dl(%c)\n", __func__, res->n - p_n, // fprintf(stderr, "[M::%s::n->%ld]\tutg%.6dl(%c)\t->\tutg%.6dl(%c)\n", __func__, res->n - p_n,
// (int32_t)(l0->v>>1)+1, "+-"[l0->v&1], (int32_t)(l1->v>>1)+1, "+-"[l1->v&1]); // (int32_t)(l0->v>>1)+1, "+-"[l0->v&1], (int32_t)(l1->v>>1)+1, "+-"[l1->v&1]);
assert(res->n - p_n > 1); assert(dst.target_hash == dst.hash); res->n--;
// fprintf(stderr, "\n-[M::%s::res->n->%u::p_n->%ld] utg%.6dl(v:%u) -> utg%.6dl(v:%u)\n",
// __func__, (uint32_t)res->n, p_n, (int32_t)(l1->v>>1)+1, l1->v, (int32_t)(l0->v>>1)+1, l0->v);
// dbg_print(res->a + p_n, res->n - p_n);
assert(res->n - p_n > 1); assert(dst.target_hash == dst.hash);
res->a[p_n-1].d = res->a[res->n-1].d - res->a[res->n-2].d; res->n--;
reverse_track(res->a + p_n, res->n - p_n); res->n--;///reomve l1 from res reverse_track(res->a + p_n, res->n - p_n); res->n--;///reomve l1 from res
kv_pushp(mg_pathv_t, *res, &p); p->v = p->d = (uint32_t)-1; p->d = k; // trav_occ++;
} else {
res->a[p_n-1].d = dst.target_dist;
} }
kv_pushp(mg_pathv_t, *res, &p); p->v = (uint32_t)-1; p->pre = k; p->d = 0;
}
// fprintf(stderr, "[M::%s::]\ta_n:%ld\ttrav_occ:%ld\n", __func__, a_n, trav_occ);
return res->n; return res->n;
} }
@@ -5672,25 +6017,241 @@ void print_chain(mg_lchain_t *a, uint32_t a_n)
} }
} }
void dedup_second_chain(uint64_t *a, int64_t a_n, int64_t p_sidx, int64_t p_eidx, mg_lchain_t *chain_a, void update_exist_chain(const ul_idx_t *uref, ul_ov_t *ch, uint64_t *idx, int64_t idx_n, int64_t tid, int64_t bw,
kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn) double diff_ec_ul, mg_lchain_t *res)
{ {
int64_t k, i; int64_t i, j, cov_i, i_qs, i_qe, i_ts, i_te, j_qs, j_qe, j_ts, j_te, dq, dt, dd, mm, sc = 0;
int64_t tlen = uref->ug->g->seq[tid].len; memset(res, 0, sizeof(*res));
ul_ov_t *li, *lj;
if(idx_n <= 0) return;
i = idx_n - 1;
res->qs = (ch[idx[i]].qs<<1)>>1; res->qe = ch[idx[i]].qe;
res->rs = ch[idx[i]].ts; res->re = ch[idx[i]].te;
for (; i >= 0; i--) {
li = &(ch[idx[i]]);
i_qs = (li->qs<<1)>>1; i_qe = li->qe;
i_ts = (li->rev?(tlen-li->te):(li->ts));
i_te = (li->rev?(tlen-li->ts):(li->te));
if((int64_t)((li->qs<<1)>>1) < res->qs) res->qs = ((li->qs<<1)>>1);
if((int64_t)li->ts < res->rs) res->rs = li->ts;
if((int64_t)li->qe > res->qe) res->qe = li->qe;
if((int64_t)li->te > res->re) res->re = li->te;
cov_i = 0; j = i + 1;
if(j < idx_n) {
lj = &(ch[idx[j]]);
j_qs = (lj->qs<<1)>>1; j_qe = lj->qe;
j_ts = (lj->rev?(tlen-lj->te):(lj->ts));
j_te = (lj->rev?(tlen-lj->ts):(lj->te));
if(j_qs <= i_qs && j_qe <= i_qe && j_ts <= i_ts && j_te <= i_te) {///co-linear
assert(li->rev == lj->rev);
if(j_qs == i_qs && j_qe == i_qe && j_ts == i_ts && j_te == i_te) continue;
dq = i_qe - j_qs; dt = i_te - j_ts;
dd = (dq>dt? dq-dt:dt-dq);
mm = MAX(dq, dt); mm *= diff_ec_ul; if(mm < bw) mm = bw;
if(dd <= mm) {///pass distance checking
sc += get_add_cov_score(uref, lj->ts, lj->te, li->ts, li->te, tid, &cov_i);
}
}
} else {
sc += retrieve_u_cov_region(uref, tid, 0, li->ts, li->te, &cov_i);
}
}
res->score = (sc>0x3FFFFFFF?0x3FFFFFFF:sc);
}
void debug_ll_chains(const ul_idx_t *uref, uint64_t *ix, int64_t ix_n, int64_t p_sidx, int64_t p_eidx, mg_lchain_t *chain_a,
kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn, uint64_t *b, int64_t bw, double diff_ec_ul, int64_t qlen)
{
int64_t k, i, z, a_n, ss, ee, b_n; uint64_t qs, qe, ts, te; uint32_t mk = 0x80000000; mg_lchain_t nn;
int64_t iqs, iqe, its, ite, tsc;
for (k = p_sidx; k < p_eidx; k++) { for (k = p_sidx; k < p_eidx; k++) {
i = raw_idx->a[chain_a[k].off].tn; i = raw_idx->a[chain_a[k].off].qn;
qs = raw_chn->a[i].qs; qe = raw_chn->a[i].qe;
ts = raw_chn->a[i].ts; te = raw_chn->a[i].te;
// fprintf(stderr, "++++++++++++[M::%s::idx:%ld]\n", __func__, i);
for (;i>=0;) { for (;i>=0;) {
i = raw_chn->a[i].qn == (uint32_t)-1?-1:raw_chn->a[i].qn; // fprintf(stderr, "--[M::%s::i->%ld]\n", __func__, i);
if(raw_chn->a[i].qs < qs) qs = raw_chn->a[i].qs;
if(raw_chn->a[i].ts < ts) ts = raw_chn->a[i].ts;
if(raw_chn->a[i].qe > qe) qe = raw_chn->a[i].qe;
if(raw_chn->a[i].te > te) te = raw_chn->a[i].te;
raw_chn->a[i].qs |= mk;
if(raw_chn->a[i].tn == (uint32_t)-1) i = -1;
else i = raw_chn->a[i].tn;
}
assert(raw_idx->a[chain_a[k].off].qs == qs && raw_idx->a[chain_a[k].off].qe == qe &&
raw_idx->a[chain_a[k].off].ts == ts && raw_idx->a[chain_a[k].off].te == te);
}
///dedup
for (z = a_n = 0; z < ix_n; ++z) {
ss = a_n; ee = a_n + ((uint32_t)ix[z]); tsc = 0;
/**if(ss != p_sidx || ee != p_eidx)**/ {
for (k = ss; k < ee; k++) {
i = raw_idx->a[chain_a[k].off].qn; b_n = 0;
qs = ((raw_chn->a[i].qs<<1)>>1); qe = raw_chn->a[i].qe;
ts = raw_chn->a[i].ts; te = raw_chn->a[i].te;
// fprintf(stderr, "++++++++++++[M::%s::idx:%ld]\n", __func__, i);
for (;i>=0;) {
// fprintf(stderr, "--[M::%s::i->%ld]\n", __func__, i);
if(((raw_chn->a[i].qs<<1)>>1) < qs) qs = ((raw_chn->a[i].qs<<1)>>1);
if(raw_chn->a[i].ts < ts) ts = raw_chn->a[i].ts;
if(raw_chn->a[i].qe > qe) qe = raw_chn->a[i].qe;
if(raw_chn->a[i].te > te) te = raw_chn->a[i].te;
// raw_chn->a[i].qs |= mk;
//update here!!!!!!!
// if(!(raw_chn->a[i].qs&mk)) b[b_n++] = i;
b[b_n++] = i;
if(raw_chn->a[i].tn == (uint32_t)-1) i = -1;
else i = raw_chn->a[i].tn;
}
assert(raw_idx->a[chain_a[k].off].qs == qs && raw_idx->a[chain_a[k].off].qe == qe &&
raw_idx->a[chain_a[k].off].ts == ts && raw_idx->a[chain_a[k].off].te == te);
update_exist_chain(uref, raw_chn->a, b, b_n, raw_idx->a[chain_a[k].off].tn, bw, diff_ec_ul, &nn);
nn.v = (raw_idx->a[chain_a[k].off].tn<<1)|raw_idx->a[chain_a[k].off].rev;
extend_end_coord(&nn, NULL, qlen, uref->ug->g->seq[raw_idx->a[chain_a[k].off].tn].len,
&iqs, &iqe, &its, &ite);
nn.qs = iqs; nn.qe = iqe; nn.rs = its; nn.re = ite;
if(!(chain_a[k].score == nn.score && chain_a[k].qs == nn.qs && chain_a[k].qe == nn.qe
&& chain_a[k].rs == nn.rs && chain_a[k].re == nn.re)){
fprintf(stderr, "[M::%s::] chain_a[k].score->%d, nn.score->%d\n", __func__, chain_a[k].score, nn.score);
fprintf(stderr, "[M::%s::] chain_a[k].qs->%d, nn.qs->%d, chain_a[k].qe->%d, nn.qe->%d, chain_a[k].rs->%d, nn.rs->%d, chain_a[k].re->%d, nn.re->%d\n", __func__,
chain_a[k].qs, nn.qs, chain_a[k].qe, nn.qe, chain_a[k].rs, nn.rs, chain_a[k].re, nn.re);
}
assert(chain_a[k].score == nn.score && chain_a[k].qs == nn.qs && chain_a[k].qe == nn.qe
&& chain_a[k].rs == nn.rs && chain_a[k].re == nn.re);
tsc += nn.score;
}
assert(tsc >= ((int64_t)(ix[z]>>32)));
}
a_n += ((uint32_t)ix[z]);
}
for (k = p_sidx; k < p_eidx; k++) {
i = raw_idx->a[chain_a[k].off].qn;
qs = ((raw_chn->a[i].qs<<1)>>1); qe = raw_chn->a[i].qe;
ts = raw_chn->a[i].ts; te = raw_chn->a[i].te;
for (;i>=0;) {
// fprintf(stderr, "--[M::%s::i->%ld]\n", __func__, i);
if(raw_chn->a[i].qs&mk) raw_chn->a[i].qs -= mk;
if(raw_chn->a[i].qs < qs) qs = raw_chn->a[i].qs;
if(raw_chn->a[i].ts < ts) ts = raw_chn->a[i].ts;
if(raw_chn->a[i].qe > qe) qe = raw_chn->a[i].qe;
if(raw_chn->a[i].te > te) te = raw_chn->a[i].te;
if(raw_chn->a[i].tn == (uint32_t)-1) i = -1;
else i = raw_chn->a[i].tn;
}
assert(raw_idx->a[chain_a[k].off].qs == qs && raw_idx->a[chain_a[k].off].qe == qe &&
raw_idx->a[chain_a[k].off].ts == ts && raw_idx->a[chain_a[k].off].te == te);
} }
} }
void dedup_second_chain(const ul_idx_t *uref, uint64_t *ix, int64_t ix_n, int64_t p_sidx, int64_t p_eidx, mg_lchain_t *chain_a,
kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn, uint64_t *b, int64_t bw, double diff_ec_ul, int64_t qlen)
{
int64_t k, i, z, a_n, ss, ee, b_n; uint64_t qs, qe, ts, te; uint32_t mk = 0x80000000; mg_lchain_t nn;
int64_t iqs, iqe, its, ite, tsc;
for (k = p_sidx; k < p_eidx; k++) {
i = raw_idx->a[chain_a[k].off].qn;
qs = raw_chn->a[i].qs; qe = raw_chn->a[i].qe;
ts = raw_chn->a[i].ts; te = raw_chn->a[i].te;
// fprintf(stderr, "++++++++++++[M::%s::idx:%ld]\n", __func__, i);
for (;i>=0;) {
// fprintf(stderr, "--[M::%s::i->%ld]\n", __func__, i);
if(raw_chn->a[i].qs < qs) qs = raw_chn->a[i].qs;
if(raw_chn->a[i].ts < ts) ts = raw_chn->a[i].ts;
if(raw_chn->a[i].qe > qe) qe = raw_chn->a[i].qe;
if(raw_chn->a[i].te > te) te = raw_chn->a[i].te;
raw_chn->a[i].qs |= mk;
if(raw_chn->a[i].tn == (uint32_t)-1) i = -1;
else i = raw_chn->a[i].tn;
}
assert(raw_idx->a[chain_a[k].off].qs == qs && raw_idx->a[chain_a[k].off].qe == qe &&
raw_idx->a[chain_a[k].off].ts == ts && raw_idx->a[chain_a[k].off].te == te);
} }
uint32_t gen_max_gchain(void *km, int64_t ulid, st_mt_t *idx, vec_mg_lchain_t *e, int64_t qlen, float primary_cov_rate,
float primary_fragment_cov_rate, float primary_fragment_second_score_rate, const asg_t *g, st_mt_t *dst_done, ///dedup
vec_sp_node_t *out, vec_mg_pathv_t *res) for (z = a_n = 0; z < ix_n; ++z) {
ss = a_n; ee = a_n + ((uint32_t)ix[z]); tsc = 0;
if(ss != p_sidx || ee != p_eidx) {
for (k = ss; k < ee; k++) {
i = raw_idx->a[chain_a[k].off].qn; b_n = 0;
qs = ((raw_chn->a[i].qs<<1)>>1); qe = raw_chn->a[i].qe;
ts = raw_chn->a[i].ts; te = raw_chn->a[i].te;
// fprintf(stderr, "++++++++++++[M::%s::idx:%ld]\n", __func__, i);
for (;i>=0;) {
// fprintf(stderr, "--[M::%s::i->%ld]\n", __func__, i);
if(((raw_chn->a[i].qs<<1)>>1) < qs) qs = ((raw_chn->a[i].qs<<1)>>1);
if(raw_chn->a[i].ts < ts) ts = raw_chn->a[i].ts;
if(raw_chn->a[i].qe > qe) qe = raw_chn->a[i].qe;
if(raw_chn->a[i].te > te) te = raw_chn->a[i].te;
// raw_chn->a[i].qs |= mk;
//update here!!!!!!!
if(!(raw_chn->a[i].qs&mk)) b[b_n++] = i;
if(raw_chn->a[i].tn == (uint32_t)-1) i = -1;
else i = raw_chn->a[i].tn;
}
assert(raw_idx->a[chain_a[k].off].qs == qs && raw_idx->a[chain_a[k].off].qe == qe &&
raw_idx->a[chain_a[k].off].ts == ts && raw_idx->a[chain_a[k].off].te == te);
update_exist_chain(uref, raw_chn->a, b, b_n, raw_idx->a[chain_a[k].off].tn, bw, diff_ec_ul, &nn);
nn.v = (raw_idx->a[chain_a[k].off].tn<<1)|raw_idx->a[chain_a[k].off].rev;
extend_end_coord(&nn, NULL, qlen, uref->ug->g->seq[raw_idx->a[chain_a[k].off].tn].len, &iqs, &iqe, &its, &ite);
nn.qs = iqs; nn.qe = iqe; nn.rs = its; nn.re = ite;
assert(chain_a[k].score >= nn.score);
tsc += (chain_a[k].score - nn.score);
}
}
///TODO: also update qs, qe
tsc = ((int64_t)(ix[z]>>32)) - tsc; if(tsc < 0) tsc = 0;
ix[z] <<= 32; ix[z] >>= 32; ix[z] |= ((uint64_t)tsc)<<32;
a_n += ((uint32_t)ix[z]);
}
for (k = p_sidx; k < p_eidx; k++) {
i = raw_idx->a[chain_a[k].off].qn;
qs = ((raw_chn->a[i].qs<<1)>>1); qe = raw_chn->a[i].qe;
ts = raw_chn->a[i].ts; te = raw_chn->a[i].te;
for (;i>=0;) {
// fprintf(stderr, "--[M::%s::i->%ld]\n", __func__, i);
if(raw_chn->a[i].qs&mk) raw_chn->a[i].qs -= mk;
if(raw_chn->a[i].qs < qs) qs = raw_chn->a[i].qs;
if(raw_chn->a[i].ts < ts) ts = raw_chn->a[i].ts;
if(raw_chn->a[i].qe > qe) qe = raw_chn->a[i].qe;
if(raw_chn->a[i].te > te) te = raw_chn->a[i].te;
if(raw_chn->a[i].tn == (uint32_t)-1) i = -1;
else i = raw_chn->a[i].tn;
}
assert(raw_idx->a[chain_a[k].off].qs == qs && raw_idx->a[chain_a[k].off].qe == qe &&
raw_idx->a[chain_a[k].off].ts == ts && raw_idx->a[chain_a[k].off].te == te);
}
}
uint32_t gen_max_gchain(void *km, const ul_idx_t *uref, int64_t ulid, st_mt_t *idx, vec_mg_lchain_t *e, kv_ul_ov_t *raw_idx, kv_ul_ov_t *raw_chn,
int64_t qlen, float primary_cov_rate, float primary_fragment_cov_rate, float primary_fragment_second_score_rate, const asg_t *g, st_mt_t *dst_done,
vec_sp_node_t *out, vec_mg_pathv_t *res, uint64_t *b, int64_t bw, double diff_ec_ul)
{ {
if(idx->n <= 0) return 0; if(idx->n <= 0) return 0;
int64_t a_n, idx_n = idx->n, i, m_sc = 0, is_done = 0; uint64_t s_idx, e_idx, om, ok, ovlp; int64_t a_n, idx_n = idx->n, i, m_sc = 0, is_done = 0; uint64_t s_idx, e_idx, om, ok, ovlp, novlp;
ul_ov_t m; memset(&m, 0, sizeof(m)); m_sc = -1; mg_lchain_t *a = e->a; ul_ov_t m; memset(&m, 0, sizeof(m)); m_sc = -1; mg_lchain_t *a = e->a;
for (i = a_n = 0; i < idx_n; ++i) { for (i = a_n = 0; i < idx_n; ++i) {
if(((int64_t)(idx->a[i]>>32)) > m_sc) { if(((int64_t)(idx->a[i]>>32)) > m_sc) {
@@ -5698,23 +6259,37 @@ vec_sp_node_t *out, vec_mg_pathv_t *res)
m.ts = a_n; m.te = a_n + ((uint32_t)idx->a[i]); m.ts = a_n; m.te = a_n + ((uint32_t)idx->a[i]);
m.qs = a[m.ts].qs; m.qe = a[m.te-1].qe; m.qs = a[m.ts].qs; m.qe = a[m.te-1].qe;
} }
// dedup_second_chain(NULL, 0, a_n, a_n + ((uint32_t)idx->a[i]), a, raw_idx, raw_chn);
a_n += ((uint32_t)idx->a[i]); a_n += ((uint32_t)idx->a[i]);
} }
assert(a[m.ts].qs<=a[m.te-1].qs && a[m.te-1].qe>=a[m.ts].qe); assert(a[m.ts].qs<=a[m.te-1].qs && a[m.te-1].qe>=a[m.ts].qe);
// print_chain(a + m.ts, m.te - m.ts); // print_chain(a + m.ts, m.te - m.ts);
///for debug
// dedup_second_chain(uref, idx->a, idx_n, m.ts, m.te, a, raw_idx, raw_chn, b, bw, diff_ec_ul);
// debug_ll_chains(uref, idx->a, idx_n, m.ts, m.te, a, raw_idx, raw_chn, b, bw, diff_ec_ul, qlen);
if((m.qe - m.qs) > (qlen*primary_cov_rate)) is_done = 1; if((m.qe - m.qs) > (qlen*primary_cov_rate)) is_done = 1;
if(is_done == 0) { if(is_done == 0) {
// for (i = a_n = 0; i < idx_n; ++i) {
// s_idx = a[a_n].qs; a_n += ((uint32_t)idx->a[i]); e_idx = a[a_n-1].qe;
// if(i == m.qn) continue;
// if(s_idx < m.qs || e_idx < m.qs || s_idx > m.qe || e_idx > m.qe) break;
// }
// if(i >= idx_n) is_done = 2;///no alignment that is on the left or the right side of the primary chain
for (i = a_n = 0; i < idx_n; ++i) { for (i = a_n = 0; i < idx_n; ++i) {
s_idx = a[a_n].qs; a_n += ((uint32_t)idx->a[i]); e_idx = a[a_n-1].qe; s_idx = a[a_n].qs; a_n += ((uint32_t)idx->a[i]); e_idx = a[a_n-1].qe;
if(i == m.qn) continue; if(i == m.qn) continue;
if(s_idx < m.qs || e_idx < m.qs || s_idx > m.qe || e_idx > m.qe) break; ovlp = ((MIN(m.qe, e_idx) > MAX(m.qs, s_idx))? (MIN(m.qe, e_idx) - MAX(m.qs, s_idx)):0);
novlp = (e_idx - s_idx) - ovlp;
if(novlp > ((m.qe-m.qs)*GC_OFFSET_RATE) && novlp > GC_OFFSET_POS) break;
} }
if(i >= idx_n) is_done = 2;///no alignment that is on the left or the right side of the primary chain if(i >= idx_n) is_done = 2;///no alignment that is on the left or the right side of the primary chain
} }
if(is_done == 0) { if(is_done == 0) {
if((m.qe - m.qs) > (qlen*primary_fragment_cov_rate)) { if((m.qe - m.qs) > (qlen*primary_fragment_cov_rate)) {
dedup_second_chain(uref, idx->a, idx_n, m.ts, m.te, a, raw_idx, raw_chn, b, bw, diff_ec_ul, qlen);
om = m.qe - m.qs; om = m.qe - m.qs;
for (i = a_n = 0; i < idx_n; ++i) { for (i = a_n = 0; i < idx_n; ++i) {
s_idx = a[a_n].qs; a_n += ((uint32_t)idx->a[i]); e_idx = a[a_n-1].qe; s_idx = a[a_n].qs; a_n += ((uint32_t)idx->a[i]); e_idx = a[a_n-1].qe;
@@ -5723,7 +6298,7 @@ vec_sp_node_t *out, vec_mg_pathv_t *res)
if(ovlp == 0) continue; if(ovlp == 0) continue;
ok = e_idx - s_idx; ok = e_idx - s_idx;
if(ok > om) ok = om; if(ok > om) ok = om;
if((ovlp > ok*0.25) && ((int64_t)(idx->a[i]>>32)) > (m_sc*primary_fragment_second_score_rate)) break; if((ovlp > ok*0.1) && ((int64_t)(idx->a[i]>>32)) > (m_sc*primary_fragment_second_score_rate)) break;
} }
if(i >= idx_n) is_done = 3; if(i >= idx_n) is_done = 3;
} }
@@ -5734,13 +6309,14 @@ vec_sp_node_t *out, vec_mg_pathv_t *res)
// fprintf(stderr, "--[M::%s::id->%ld] [%u, %u), res->n:%lu\n", __func__, ulid, m.qs, m.qe, (uint64_t)res->n); // fprintf(stderr, "--[M::%s::id->%ld] [%u, %u), res->n:%lu\n", __func__, ulid, m.qs, m.qe, (uint64_t)res->n);
kv_resize(mg_lchain_t, *e, res->n); a = e->a; kv_resize(mg_lchain_t, *e, res->n); a = e->a;
for (i = ((int64_t)res->n)-1; i >= 0; i--) { for (i = ((int64_t)res->n)-1; i >= 0; i--) {
if(res->a[i].v == (uint32_t)-1) { if(res->a[i].v == (uint32_t)-1) {
a[i] = a[res->a[i].d]; a[i] = a[res->a[i].pre]; a[i].dist_pre = res->a[i].d;
// fprintf(stderr, "ulid:%ld\t%u\t%u\t%c\tutg%.6dl\t%u\t%u\n", ulid, a[i].qs, a[i].qe, "+-"[a[i].v&1], (int32_t)(a[i].v>>1)+1, a[i].rs, a[i].re); // fprintf(stderr, "ulid:%ld\t%u\t%u\t%c\tutg%.6dl\t%u\t%u\tdist_pre:%d\n", ulid, a[i].qs, a[i].qe, "+-"[a[i].v&1], (int32_t)(a[i].v>>1)+1, a[i].rs, a[i].re, a[i].dist_pre);
} }
else { else {
a[i].v = res->a[i].v; a[i].off = -1; a[i].v = res->a[i].v; a[i].off = -1; a[i].dist_pre = res->a[i].d;
// fprintf(stderr, "ulid:%ld\t*\t*\t%c\tutg%.6dl\t*\t*\n", ulid, "+-"[a[i].v&1], (int32_t)(a[i].v>>1)+1); // fprintf(stderr, "ulid:%ld\t*\t*\t%c\tutg%.6dl\t*\t*\tdist_pre:%d\n", ulid, "+-"[a[i].v&1], (int32_t)(a[i].v>>1)+1, a[i].dist_pre);
} }
} }
e->n = res->n; e->n = res->n;
@@ -5748,7 +6324,7 @@ vec_sp_node_t *out, vec_mg_pathv_t *res)
// debug_gchain(g, e->a, e->n); // debug_gchain(km, g, e->a, e->n, dst_done, out);
return 1; return 1;
@@ -5777,7 +6353,7 @@ void update_ul_vec_t(ul_vec_t *rch, vec_mg_lchain_t *u, const ul_idx_t *uref)
uint32_t direct_gchain(mg_tbuf_t *b, ul_vec_t *rch, glchain_t *ll, gdpchain_t *gdp, st_mt_t *sps, haplotype_evdience_alloc *hap, const ul_idx_t *uref, const ug_opt_t *uopt, uint32_t direct_gchain(mg_tbuf_t *b, ul_vec_t *rch, glchain_t *ll, gdpchain_t *gdp, st_mt_t *sps, haplotype_evdience_alloc *hap, const ul_idx_t *uref, const ug_opt_t *uopt,
int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid) int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid)
{ {
// if(ulid!=108) return 0; // if(ulid != 814) 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); gl_rg2ug_gen(rch, idx, uref, 1);
@@ -5797,22 +6373,24 @@ int64_t bw, double diff_ec_ul, int64_t max_skip, int64_t ulid)
///buffer ///buffer
kv_resize(uint64_t, ll->srt.a, gdp->l.n); kv_resize(uint64_t, hap->snp_srt, gdp->l.n); kv_resize(uint64_t, ll->srt.a, gdp->l.n); kv_resize(uint64_t, hap->snp_srt, gdp->l.n);
kv_resize(uint64_t, gdp->v, gdp->l.n); kv_resize(int64_t, gdp->f, gdp->l.n); kv_resize(uint64_t, gdp->v, gdp->l.n); kv_resize(int64_t, gdp->f, gdp->l.n);
max_idx = hc_gchain1_dp(b->km, uref->ug, &(gdp->l), &(gdp->swap), &(gdp->dst), &(gdp->out), &(gdp->path), rch->rlen, max_idx = hc_gchain1_dp(b->km, uref, uref->ug, &(gdp->l), &(gdp->swap), &(gdp->dst), &(gdp->out), &(gdp->path), rch->rlen,
uopt, bw, diff_ec_ul, ll->srt.a.a, sps, gdp->f.a, hap->snp_srt.a, gdp->v.a); uopt, bw, diff_ec_ul, N_GCHAIN_RATE, ll->srt.a.a, sps, gdp->f.a, hap->snp_srt.a, gdp->v.a);
// fprintf(stderr, "++++[M::%s::id->%ld, len->%u] gdp->l.n:%lu\n", __func__, ulid, rch->rlen, (uint64_t)gdp->l.n); // fprintf(stderr, "++++[M::%s::id->%ld, len->%u] gdp->l.n:%lu\n", __func__, ulid, rch->rlen, (uint64_t)gdp->l.n);
// fprintf(stderr, "+[M::%s::] gdp->l.n:%lu\n", __func__, (uint64_t)gdp->l.n); // fprintf(stderr, "+[M::%s::] gdp->l.n:%lu\n", __func__, (uint64_t)gdp->l.n);
//sps has the chain idx; gdp->l has the chain //sps has the chain idx; gdp->l has the chain
if(max_idx >= 0 && gen_max_gchain(b->km, ulid, sps, &(gdp->l), rch->rlen, P_CHAIN_COV, P_FRAGEMENT_PRIMARY_CHAIN_COV, if(max_idx >= 0 && gen_max_gchain(b->km, uref, ulid, sps, &(gdp->l), idx, init, rch->rlen, P_CHAIN_COV, P_FRAGEMENT_PRIMARY_CHAIN_COV,
P_FRAGEMENT_PRIMARY_SECOND_COV, uref->ug->g, &(gdp->dst_done), &(gdp->out), &(gdp->path))) { 0.1/**P_FRAGEMENT_PRIMARY_SECOND_COV**/, uref->ug->g, &(gdp->dst_done), &(gdp->out), &(gdp->path), ll->srt.a.a, bw, diff_ec_ul)) {
// update_ul_vec_t(rch, &(gdp->l), uref); // update_ul_vec_t(rch, &(gdp->l), uref);
// __ac_X31_hash_string("hehe"); // __ac_X31_hash_string("hehe");
return 1;
} else { } else {
uint64_t i; uint64_t i;
fprintf(stderr, "unsuccess->[M::%s::id->%ld, len->%u] gdp->l.n:%lu\n", __func__, ulid, rch->rlen, (uint64_t)gdp->l.n); fprintf(stderr, "unsuccess->[M::%s::id->%ld, len->%u] gdp->l.n:%lu\n", __func__, ulid, rch->rlen, (uint64_t)gdp->l.n);
for (i = 0; i < gdp->l.n; ++i) { for (i = 0; i < gdp->l.n; ++i) {
fprintf(stderr, "(%lu)\tutg%.6d%c(%u)\t%u\t%u\t%c\tsrc:%u\tscore:%d\n", fprintf(stderr, "(%lu)\t%u\t%u\t%c\tutg%.6d%c(%u)\t%u\t%u\tsrc:%u\tscore:%d\n",
i, (int32_t)(gdp->l.a[i].v>>1)+1, "lc"[uref->ug->u.a[gdp->l.a[i].v>>1].circ], uref->ug->u.a[gdp->l.a[i].v>>1].len, i, gdp->l.a[i].qs, gdp->l.a[i].qe, "+-"[gdp->l.a[i].v&1], (int32_t)(gdp->l.a[i].v>>1)+1, "lc"[uref->ug->u.a[gdp->l.a[i].v>>1].circ], uref->ug->u.a[gdp->l.a[i].v>>1].len,
gdp->l.a[i].qs, gdp->l.a[i].qe, "+-"[gdp->l.a[i].v&1], gdp->l.a[i].v^1, gdp->l.a[i].score); gdp->l.a[i].rs, gdp->l.a[i].re, gdp->l.a[i].v^1, gdp->l.a[i].score);
} }
} }
@@ -5837,7 +6415,8 @@ static void worker_for_ul_gchains_alignment(void *data, long i, int tid)
if(p->dd == 1) return; //fully aligned if(p->dd == 1) return; //fully aligned
if(p->bb.n == 1 && p->bb.a[0].base) return;///no alignment if(p->bb.n == 1 && p->bb.a[0].base) return;///no alignment
utepdat_t *s = (utepdat_t*)data; utepdat_t *s = (utepdat_t*)data;
direct_gchain(s->buf[tid], p, &(s->ll[tid]), &(s->gdp[tid]), &(s->sps[tid]), &(s->hab[tid]->hap), s->uu, s->uopt, G_CHAIN_BW, s->opt->diff_ec_ul, UG_SKIP, i); s->sum_len++;
s->n += direct_gchain(s->buf[tid], p, &(s->ll[tid]), &(s->gdp[tid]), &(s->sps[tid]), &(s->hab[tid]->hap), s->uu, s->uopt, G_CHAIN_BW, s->opt->diff_ec_ul, UG_SKIP, i);
// gl_chain_refine_advance(&b->olist, &b->correct, &b->hap, bl, s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, s->id+i, km); // gl_chain_refine_advance(&b->olist, &b->correct, &b->hap, bl, s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, s->id+i, km);
} }
@@ -5860,6 +6439,7 @@ void work_ul_gchains(uldat_t *sl)
} }
free(s.hab); free(s.buf); free(s.ll); free(s.gdp); free(s.mzs); free(s.sps); free(s.hab); free(s.buf); free(s.ll); free(s.gdp); free(s.mzs); free(s.sps);
fprintf(stderr, "[M::%s::] # try:%d, # done:%d\n", __func__, s.sum_len, s.n);
} }
void print_ul_ovlps(all_ul_t *x, int32_t prt_ovlp) void print_ul_ovlps(all_ul_t *x, int32_t prt_ovlp)