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

Merge pull request #363 from chhylp123/hifiasm_dev_debug

Browse Source 
cleaner assembly graph
This commit is contained in:
chhylp123
2022-12-06 09:52:47 -05:00
committed by GitHub
parent 5b928dfe1d a3b72457f0
commit 6743dcf204
8 changed files with 426 additions and 38 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.cpp
+3
View File
@@ -172,8 +172,10 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->k_mer_length = 51;
asm_opt->hic_mer_length = 31;
asm_opt->ul_mer_length = 19;
asm_opt->trans_mer_length = 31;
asm_opt->mz_win = 51;
asm_opt->ul_mz_win = 19;
asm_opt->trans_win = 31;
asm_opt->mz_rewin = 1000;
asm_opt->ul_mz_rewin = 360;
asm_opt->mz_sample_dist = 500;
@@ -208,6 +210,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->purge_level_trio = 0;
asm_opt->purge_simi_rate_l2 = 0.75;
asm_opt->purge_simi_rate_l3 = 0.55;
asm_opt->trans_base_rate = 0.8;
asm_opt->purge_overlap_len = 1;
///asm_opt->purge_overlap_len_hic = 50;
asm_opt->recover_atg_cov_min = -1024;
CommandLines.h
+4 -1
View File
@@ -4,7 +4,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.18.0-r465"
#define HA_VERSION "0.18.1-r466"
#define VERBOSE 0
@@ -48,9 +48,11 @@ typedef struct {
int k_mer_length;
int hic_mer_length;
int ul_mer_length;
int trans_mer_length;
int bub_mer_length;
int mz_win;
int ul_mz_win;
int trans_win;
int mz_rewin;
int ul_mz_rewin;
int mz_sample_dist;
@@ -97,6 +99,7 @@ typedef struct {
float purge_simi_rate_l2;
float purge_simi_rate_l3;
float purge_simi_thres;
float trans_base_rate;
///float purge_simi_rate_hic;
Overlaps.cpp
+24 -1
View File
@@ -13823,6 +13823,8 @@ long long gap_fuzz, bub_label_t* b_mask_t)
if((asm_opt.flag & HA_F_VERBOSE_GFA)) write_trans_chain(cov->t_ch, output_file_name);
}
///for debug
// char* gfa_name = (char*)malloc(strlen(output_file_name)+50);
// sprintf(gfa_name, "%s.pre.clean_d_utg.noseq.gfa", output_file_name);
@@ -13830,7 +13832,11 @@ long long gap_fuzz, bub_label_t* b_mask_t)
// ma_ug_print_simple(ug, sg, coverage_cut, sources, ruIndex, "utg", output_file);
// fclose(output_file);
// free(gfa_name);
// exit(1);
///for debug
hic_analysis(ug, sg, cov?cov->t_ch:t_ch, &opt, 0, asm_opt.scffold?&rhits:NULL);
@@ -14876,12 +14882,29 @@ void filter_u_trans_t(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, uint32_t thr
kt_u_trans_t_idx(ta, ug->g->n_seq);
}
void dbg_prt_utg_trans(kv_u_trans_t *ta, ma_ug_t *ug, const char *o_n)
{
char* gfa_name = (char*)malloc(strlen(o_n)+100);
FILE *fn = NULL; sprintf(gfa_name, "%s.tran.dbg.ovlp.log", o_n);
u_trans_t *p = NULL; uint32_t i; fn = fopen(gfa_name, "w");
for (i = 0; i < ta->n; i++) {
p = &(ta->a[i]);
fprintf(fn, "utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\tw(%f)\tf(%u)\n",
p->qn+1, "lc"[ug->u.a[p->qn].circ], ug->u.a[p->qn].len, p->qs, p->qe, "+-"[p->rev],
p->tn+1, "lc"[ug->u.a[p->tn].circ], ug->u.a[p->tn].len, p->ts, p->te, p->nw, p->f);
}
fclose(fn); free(gfa_name); fprintf(stderr, "[M::%s::] done\n", __func__);
}
void clean_u_trans_t_idx(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g)
{
// dbg_prt_utg_trans(ta, ug, "pre");
kt_u_trans_t_idx(ta, ug->g->n_seq);
kt_u_trans_t_symm(ta, ug);
filter_u_trans_t(ta, ug, read_g, 3);
///debug_u_trans_t(ta);
// dbg_prt_utg_trans(ta, ug, "after");
}
@@ -32112,7 +32135,7 @@ ma_sub_t **coverage_cut_ptr, int debug_g)
// flat_bubbles(sg, ruIndex->is_het); free(ruIndex->is_het); ruIndex->is_het = NULL;
flat_soma_v(sg, sources, ruIndex);
**/
if(!(ha_opt_triobin(&asm_opt) && ha_opt_hic(&asm_opt))) {
if(!ha_opt_triobin(&asm_opt)) {
// output_unitig_graph(sg, coverage_cut, "pre_clean", sources, ruIndex, max_hang_length, mini_overlap_length);
hic_clean_adv(sg, &uopt);
}
Purge_Dups.cpp
+1 -1
View File
@@ -1788,7 +1788,7 @@ void chain_trans_ovlp(hap_cov_t *cover, utg_trans_t *o, ma_ug_t *ug, asg_t *read
}
p_v = v; len += l;
if(len >= targetBaseLen)
if(len >= targetBaseLen)///might be not right for centromeres
{
xOcc = aOcc;
break;
anchor.cpp
+13
View File
@@ -1278,3 +1278,16 @@ void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t
///no need to sort here, overlap_list has been sorted at lchain_gen
lchain_qgen_mcopy(cl, overlap_list, rid, rl, NULL, uref, apend_be, max_n_chain, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, quick_check, gen_off);
}
void ug_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres,
int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off, uint32_t is_hpc)
{
int64_t max_skip, max_iter, max_dis, quick_check; double chn_pen_gap, chn_pen_skip;
set_lchain_dp_op(is_accurate, mz_k, &max_skip, &max_iter, &max_dis, &chn_pen_gap, &chn_pen_skip, &quick_check);
if((!overlap_list) || (!cl)) {
ab->mz.n = 0;
mz2_ha_sketch(rs, rl, mz_w, mz_k, 0, is_hpc, &ab->mz, NULL, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0, NULL);
} else {
}
}
gfa_ut.cpp
+192 -31
View File
@@ -11510,18 +11510,20 @@ void debug_sysm_usg_t(usg_t *ng, const char *cmd)
}
}
void update_usg_t_threading_0(usg_t *ng, uint64_t *a, uint64_t a_n, uint32_t *occ, asg64_v *b)
void update_usg_t_threading_0(ul_resolve_t *uidx, usg_t *ng, uint64_t *a, uint64_t a_n, uint32_t *occ, asg64_v *b)
{
uint64_t k, i, *p, nid, nnid, v, w; b->n = 0; usg_seq_t *s;
assert(a_n > 1);
assert(occ[a[0]] == 1);
kv_pushp(uint64_t, *b, &p); (*p) = a[0]; (*p) <<= 32; (*p) |= a[0]; occ[a[0]]--;
// fprintf(stderr, "[M::%s::] a_n::%lu\n", __func__, a_n);
// for (k = 0; k < a_n; k++) {
// fprintf(stderr, "utg%.6dl(%c)\t", (int32_t)(a[k]>>1)+1, "+-"[a[k]&1]);
// fprintf(stderr, "utg%.6dl(%c)(hom::%u)\t", (int32_t)(a[k]>>1)+1, "+-"[a[k]&1], (IF_HOM((a[k]>>1), (*(uidx->bub)))));
// }
// fprintf(stderr, "\n");
assert(a_n > 1);
assert(occ[a[0]] == 1);
kv_pushp(uint64_t, *b, &p); (*p) = a[0]; (*p) <<= 32; (*p) |= a[0]; occ[a[0]]--;
for (k = 1; k + 1 < a_n; k++) {
assert(occ[a[k]] == occ[a[k]^1]); assert(occ[a[k]] > 0);
if(occ[a[k]] > 1) {///copy a node
@@ -11627,16 +11629,29 @@ void update_usg_t_threading(ul_resolve_t *uidx, usg_t *ng, uint64_t *arcs, uint6
{
uint64_t i, k, s, e, nvtx = ng->n<<1; memset(occ, 0, sizeof((*occ))*nvtx);
for (i = 0; i < arcs_gn; i++) {///set cluster
// prt = 0;
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]
occ[integ_seq[k]]++; occ[integ_seq[k]^1]++;
// if((integ_seq[k]>>1) == 2736) prt = 1;
}
occ[integ_seq[s]]++; occ[integ_seq[e]^1]++;
// if((integ_seq[s]>>1) == 2736) prt = 1;
// if(((integ_seq[e]^1)>>1) == 2736) prt = 1;
// if(occ[integ_seq[s]] > 1 || occ[integ_seq[e]^1] > 1) {
// fprintf(stderr, "s::utg%.6dl(%c)(occ::%u), e::utg%.6dl(%c)(occ::%u)\n",
// (int32_t)(integ_seq[s]>>1)+1, "+-"[integ_seq[s]&1], occ[integ_seq[s]],
// (int32_t)((integ_seq[e]^1)>>1)+1, "+-"[((integ_seq[e]^1)&1)], occ[integ_seq[e]^1]);
// }
// if(prt) {
// fprintf(stderr, "[M::%s::] a_n::%lu\n", __func__, e + 1 - s);
// for (k = s; k <= e; k++) {
// fprintf(stderr, "utg%.6dl(%c)(hom::%u)\t",
// (int32_t)(integ_seq[k]>>1)+1, "+-"[integ_seq[k]&1],
// (IF_HOM((integ_seq[k]>>1), (*(uidx->bub)))));
// }
// fprintf(stderr, "\n");
// }
}
for (i = 0; i < nvtx; i++) {
@@ -11644,7 +11659,7 @@ void update_usg_t_threading(ul_resolve_t *uidx, usg_t *ng, uint64_t *arcs, uint6
}
for (i = 0; i < arcs_gn; i++) {///arcs_g[]>>32 is the group id; (uint32_t)arcs_g[]
s = arcs[((uint32_t)arcs_g[i])]>>32; e = ((uint32_t)arcs[((uint32_t)arcs_g[i])]);
update_usg_t_threading_0(ng, integ_seq + s, e + 1 - s, occ, b);
update_usg_t_threading_0(uidx, ng, integ_seq + s, e + 1 - s, occ, b);
}
}
@@ -12245,39 +12260,68 @@ uint64_t get_arc_support(ul_resolve_t *uidx, uint64_t v, uint64_t w)
// }
// }
uint64_t get_arc_support_chain(ul_resolve_t *uidx, uint64_t *a, uint64_t a_n, usg_t *ng)
{
if(a_n <= 0) return 0;
uint64_t *hid_a, hid_n; ul_str_idx_t *str_idx = &(uidx->pstr);
uint64_t z, vz, occ = 0, v, w, ai; ul_str_t *str; int64_t s, s_n;
v = (ng->a[a[0]>>1].mm<<1)|(a[0]&1);
hid_a = str_idx->occ.a + str_idx->idx.a[v>>1];
hid_n = str_idx->idx.a[(v>>1)+1] - str_idx->idx.a[v>>1];
for (z = occ = 0; z < hid_n; z++) {
str = &(str_idx->str.a[hid_a[z]>>32]); s_n = str->cn;
if(s_n < 2) continue;
vz = (uint32_t)(str->a[(uint32_t)hid_a[z]]);
assert((v>>1) == (vz>>1));
if(v == vz) {
for(s = ((uint32_t)hid_a[z]) + 1, ai = 1; s < s_n && ai < a_n; s++, ai++) {
w = (ng->a[a[ai]>>1].mm<<1)|(a[ai]&1);
if(((uint32_t)(str->a[s])) != w) break;
}
if(ai >= a_n) occ++;
} else {
for(s = ((int32_t)((uint32_t)hid_a[z]))-1, ai = 1; s >= 0 && ai < a_n; s--, ai++) {
w = (ng->a[a[ai]>>1].mm<<1)|(a[ai]&1);
if(((uint32_t)(str->a[s])) != (w^1)) break;
}
if(ai >= a_n) occ++;
}
}
return occ;
}
static void worker_unique_bridge_check_s(void *data, long i, int tid) // callback for kt_for()
{
unique_bridge_check_t *uaux = (unique_bridge_check_t *)data;
uint64_t *integer_seq = uaux->integer_seq, s, e, v, w, nse, k;
uint64_t *integer_seq = uaux->integer_seq, s, e, vk, wk, nse, k;
uint64_t *gidx = uaux->gidx, *interval = uaux->interval_idx;
bubble_type *bub = uaux->uidx->bub; usg_t *ng = uaux->ng;
s = interval[((uint32_t)gidx[i])]>>32;
e = ((uint32_t)interval[((uint32_t)gidx[i])]);
assert(s < e);
for (k = s, v = w = (uint32_t)-1; k <= e; k++) {
for (k = s, vk = wk = (uint32_t)-1; k <= e; k++) {
if(k > s) {
nse = get_arc_support(uaux->uidx, (ng->a[integer_seq[k-1]>>1].mm<<1)|(integer_seq[k-1]&1),
(ng->a[integer_seq[k]>>1].mm<<1)|(integer_seq[k]&1));
// nse = get_arc_support(uaux->uidx, (ng->a[integer_seq[k-1]>>1].mm<<1)|(integer_seq[k-1]&1),
// (ng->a[integer_seq[k]>>1].mm<<1)|(integer_seq[k]&1));
nse = get_arc_support_chain(uaux->uidx, integer_seq+k-1, 2, ng);
if(nse < unique_bridge_occ) {
// fprintf(stderr, "+utg%.6dl(%c)\tutg%.6dl(%c)\tnse::%lu\n",
// ((int32_t)(integer_seq[k-1]>>1))+1, "+-"[integer_seq[k-1]&1],
// ((int32_t)(integer_seq[k]>>1))+1, "+-"[integer_seq[k]&1], nse);
gidx[i] |= ((uint64_t)0x8000000000000000); return;
}
}
if(IF_HOM((integer_seq[k]>>1), *bub)) continue;
w = (ng->a[integer_seq[k]>>1].mm<<1)|(integer_seq[k]&1);
if(v != (uint32_t)-1) {
nse = get_arc_support(uaux->uidx, v, w);
// w = (ng->a[integer_seq[k]>>1].mm<<1)|(integer_seq[k]&1);
wk = k;
if(vk != (uint32_t)-1) {
// nse = get_arc_support(uaux->uidx, v, w);
nse = get_arc_support_chain(uaux->uidx, integer_seq+vk, wk+1-vk, ng);
if(nse < unique_bridge_occ) {
// fprintf(stderr, "-utg%.6dl(%c)\tutg%.6dl(%c)\tnse::%lu\n",
// ((int32_t)(v>>1))+1, "+-"[v&1],
// ((int32_t)(w>>1))+1, "+-"[w&1], nse);
gidx[i] |= ((uint64_t)0x8000000000000000); return;
}
}
v = w;
vk = wk;
}
}
@@ -12382,7 +12426,7 @@ uint8_t *ff, uint32_t *ng_occ, uint32_t max_ext)
for (i = 0, ub64->n = 0; i < res.n; i++) {
s = res.a[i]>>32; e = ((uint32_t)res.a[i]);
update_usg_t_threading_0(ng, int_a + s, e + 1 - s, ng_occ, ub64);
update_usg_t_threading_0(uidx, ng, int_a + s, e + 1 - s, ng_occ, ub64);
}
}
kv_destroy(res);
@@ -13011,7 +13055,8 @@ void renew_usg_t_bub(ul_resolve_t *uidx, usg_t *ng, uint32_t *id_map, uint8_t *f
// fprintf(stderr, "[M::%s::k->%u] utg%.6dl(%c), ug->u.a[k].n::%u, rocc_cut::%u, is_hom::%u\n", __func__, k,
// i+1, "+-"[0], (uint32_t)ug->u.a[k].n, rocc_cut, IF_HOM(k, (*(uidx->bub))));
// }
if((ug->u.a[k].n >= rocc_cut) && (!(IF_HOM(k, (*(uidx->bub)))))) {
if((ug->u.a[k].n >= rocc_cut) &&
((!(IF_HOM(k, (*(uidx->bub))))) || (asm_opt.purge_level_primary == 0))) {
if(id_map[i]&p[0]) {
ff[i<<1] = 2;
// fprintf(stderr, "[M::%s::] utg%.6dl(%c), ug->u.a[k].n::%u, rocc_cut::%u\n", __func__,
@@ -13110,8 +13155,8 @@ asg64_v *b0, asg64_v *b1, double cutoff)
// uint64_t is_debug = 0;
// if(((v>>1) == 10600)/** && (b0->n > 1 && (b0->a[b0->n-1]>>1) == 34944)**/) {
// is_debug = 1;
// if((((v>>1) == 2736) && (v&1))/** && (b0->n > 1 && (b0->a[b0->n-1]>>1) == 34944)**/) {
// // is_debug = 1;
// fprintf(stderr, "[M::%s::] utg%.6dl(%c)\n", __func__,
// (int32_t)(v>>1)+1, "+-"[v&1]);
// prt_sub_integer_path(b0, it, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1, (uint64_t)-1);
@@ -13130,12 +13175,13 @@ asg64_v *b0, asg64_v *b1, double cutoff)
else return 0;///b0 is contained
if(z < b1->n) get_integer_seq_ovlps(uidx, b1->a, b1->n, z - 1, 0, NULL, &w1);
else continue;///b1 is contained
// if(((v>>1) == 8722) || ((v>>1) == 56768)) {
// if(((v>>1) == 2736) && (v&1)) {
// prt_sub_integer_path(b0, it, w0, w1, zn, z);
// prt_sub_integer_path(b1, it, w0, w1, zn, z);
// }
if(w0 == (uint64_t)-1) w0 = 0;
if(w1 == (uint64_t)-1) w1 = 0;
if((w0 <= w1) || (w1 > (w0*cutoff))) return 0;
// if(b0->n == zn) return 0;///b0 is shorter
// if(b1->n == zn) continue;///b1 is shorter
if((min_w0 == (uint64_t)-1) || (z == zn) || (min_w0 > w0) || (min_w0 == w0 && min_w1 < w1)) {
@@ -13143,6 +13189,10 @@ asg64_v *b0, asg64_v *b1, double cutoff)
}
is_contain = 0;
}
// if(((v>>1) == 2736) && (v&1)) {
// fprintf(stderr, "[M::%s::] min_w0::%lu, min_w1::%lu\n\n", __func__, min_w0, min_w1);
// }
// fprintf(stderr, "[M::%s::] utg%.6dl(%c)->utg%.6dl(%c)\n", __func__,
// (int32_t)(v>>1)+1, "+-"[v&1], (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1]);
if(alt_n == 0) return 1;
@@ -13165,13 +13215,13 @@ uint64_t *ridx, asg64_v *res)
v = int_a[k];
if((!f[v])&&(!f[v^1])) continue;
if((pi != (uint64_t)-1) && (f[v^1])) {
// if(((v>>1) == 8722) || ((v>>1) == 56768)) {
// if(((v>>1) == 2736)) {
// fprintf(stderr, "+[M::%s::ii[%lu, %lu)] utg%.6dl(%c), f[v^1]::%u, v^1::%lu, putg%.6dl(%c), f[pv]::%u, pv::%lu\n",
// __func__, s, e, (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1], f[v^1], v^1,
// (int32_t)(int_a[pi]>>1)+1, "+-"[int_a[pi]&1], f[int_a[pi]], int_a[pi]);
// }
if(is_best_path(uidx, ng, int_idx, int_a, s, e, k, v^1, ridx_a, ridx, b0, b1, 0.51)) {
// if(((v>>1) == 8722) || ((v>>1) == 56768)) {
// if(((v>>1) == 2736)) {
// fprintf(stderr, "-[M::%s::ii[%lu, %lu)] utg%.6dl(%c), f[v^1]::%u, v^1::%lu, putg%.6dl(%c), f[pv]::%u, pv::%lu\n",
// __func__, s, e, (int32_t)(int_a[k]>>1)+1, "+-"[int_a[k]&1], f[v^1], v^1,
// (int32_t)(int_a[pi]>>1)+1, "+-"[int_a[pi]&1], f[int_a[pi]], int_a[pi]);
@@ -13199,6 +13249,81 @@ uint64_t *ridx, asg64_v *res)
}
}
uint64_t gen_sub_integer_path_ff(uint64_t *int_a, int64_t s, int64_t e, int64_t it, uint64_t v, uint8_t *f, asg64_v *res)
{
int64_t k; res->n = 0; assert((int_a[it]>>1) == (v>>1));
if(int_a[it] == v) {///forward
kv_resize(uint64_t, *res, (uint64_t)(e - it));
for (k = it; k < e; k++) {
res->a[res->n++] = int_a[k];
if(res->n > 1 && f[res->a[res->n-1]^1]) return res->n;
}
} else {//reverse
kv_resize(uint64_t, *res, (uint64_t)(it - s));
for (k = it; k >= s; k--) {
res->a[res->n++] = int_a[k]^1;
if(res->n > 1 && f[res->a[res->n-1]^1]) return res->n;
}
}
res->n = 0; return 0;
}
uint64_t is_best_pair(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t *int_a,
uint64_t s, uint64_t e, uint64_t it, uint64_t v, uint64_t *ridx_a, uint64_t *ridx,
asg64_v *b0, asg64_v *b1, uint8_t *f, double cutoff)
{
uint64_t *arc_a, arc_n, k, is, ie, w0 = 0, w1 = 0, sup_cut = 3;
if(!gen_sub_integer_path_ff(int_a, s, e, it, v, f, b0)) return 0;
assert(b0->n >= 2); assert(f[b0->a[0]] && f[b0->a[b0->n-1]^1]);
w0 = get_arc_support_chain(uidx, b0->a, b0->n, ng);
if(w0 < sup_cut) return 0;
arc_a = ridx_a + (ridx[v>>1]>>32);
arc_n = (uint32_t)ridx[v>>1];
for (k = 0; k < arc_n; k++) {
if(it == ((uint32_t)arc_a[k])) continue;
is = int_idx[arc_a[k]>>32]>>32;
ie = is + ((uint32_t)(int_idx[arc_a[k]>>32]));
if(!gen_sub_integer_path_ff(int_a, is, ie, ((uint32_t)arc_a[k]), v, f, b1)) continue;
assert(b1->n >= 2); assert(f[b1->a[0]] && f[b1->a[b1->n-1]^1]);
if((b0->n == b1->n) && (memcmp(b0->a, b1->a, sizeof((*(b0->a)))*b0->n) == 0)) {
if(it < ((uint32_t)arc_a[k])) continue;
else return 0;///only keep 1 equal interval
}
w1 = get_arc_support_chain(uidx, b1->a, b1->n, ng);
// fprintf(stderr, "[M::%s::] w0::%lu, w1::%lu\n", __func__, w0, w1);
if((w0 <= w1) || (w1 > (w0*cutoff))) return 0;
}
return 1;
}
void get_best_pair(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t *int_a,
uint8_t *f, uint64_t s, uint64_t e, asg64_v *b0, asg64_v *b1, uint64_t *ridx_a,
uint64_t *ridx, asg64_v *res)
{
uint64_t pi, v, k, *pz;
b0->n = b1->n = 0;
for (k = s, pi = (uint64_t)-1; k < e; k++) {
v = int_a[k];
if((f[v^1])) {
if(pi != (uint64_t)-1) {
if(is_best_pair(uidx, ng, int_idx, int_a, s, e, k, v^1, ridx_a, ridx, b0, b1, f, 0.51)) {
kv_pushp(uint64_t, *res, &pz);
*pz = pi; (*pz) <<= 32; (*pz) |= k;
}
}
pi = (uint64_t)-1;
}
if(f[v]) {
pi = (uint64_t)-1;
if(is_best_pair(uidx, ng, int_idx, int_a, s, e, k, v, ridx_a, ridx, b0, b1, f, 0.51)) {
pi = k;//keep the shortest pi<->k
}
}
}
}
static void worker_ul_aln_path(void *data, long i, int tid) // callback for kt_for()
{
unique_bridge_check_t *u_aux = (unique_bridge_check_t*)data;
@@ -13226,6 +13351,33 @@ static void worker_ul_aln_path(void *data, long i, int tid) // callback for kt_f
buf->res_dump.a = res.a; buf->res_dump.n = res.n; buf->res_dump.m = res.m;
}
static void worker_ul_aln_pair(void *data, long i, int tid) // callback for kt_for()
{
unique_bridge_check_t *u_aux = (unique_bridge_check_t*)data;
ul_resolve_t *uidx = u_aux->uidx;
integer_t *buf = &(uidx->str_b.buf[tid]);
uint64_t s, e;
// uint64_t *x = &(uidx->uovl.iug_tra->a[i]);
// asg_arc_t *ve = &(uidx->uovl.i_ug->g->arc[*x]);
asg64_v b_v, b_r, res;
b_v.a = buf->u.a; b_v.n = buf->u.n; b_v.m = buf->u.m;
b_r.a = buf->o.a; b_r.n = buf->o.n; b_r.m = buf->o.m;
res.a = buf->res_dump.a; res.n = buf->res_dump.n; res.m = buf->res_dump.m;
b_v.n = b_r.n = 0;
s = u_aux->int_idx[i]>>32; ///the i-th integer contig/path
e = s + ((uint32_t)(u_aux->int_idx[i]));
assert(e > s + 1);//the length is at least 2
get_best_pair(uidx, u_aux->ng, u_aux->int_idx, u_aux->int_a, u_aux->f, s, e, &b_v, &b_r, u_aux->ridx_a, u_aux->ridx, &res);
// get_ul_arc_supports(uidx, ve, &b_v, &b_r, 1, &w_v, &w_r);
buf->u.a = b_v.a; buf->u.n = b_v.n; buf->u.m = b_v.m;
buf->o.a = b_r.a; buf->o.n = b_r.n; buf->o.m = b_r.m;
buf->res_dump.a = res.a; buf->res_dump.n = res.n; buf->res_dump.m = res.m;
}
uint32_t ug_ext_0(ul_resolve_t *uidx, usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint64_t *int_a, uint32_t max_ext,
uint8_t *ff, uint32_t *ng_occ, uint64_t *i_idx, asg64_v *b64, asg64_v *ub64, uint64_t a_n)
{
@@ -13294,7 +13446,8 @@ void debug_prt_renew_aln(usg_t *ng, uint64_t *int_idx, uint64_t int_idx_n, uint6
}
uint32_t ug_ext_free(asg64_v *ob, asg64_v *ub, ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext,
uint8_t **ff, uint32_t **ng_occ, uint64_t **i_idx, asg64_v *b64, asg64_v *ub64, uint32_t rocc_cut)
uint8_t **ff, uint32_t **ng_occ, uint64_t **i_idx, asg64_v *b64, asg64_v *ub64, uint32_t rocc_cut,
uint32_t thread_path)
{
// fprintf(stderr, "[M::%s::] rocc_cut::%u\n", __func__, rocc_cut);
uint32_t k, n_vtx = ng->n<<1, a_n; unique_bridge_check_t u_aux;
@@ -13313,7 +13466,13 @@ uint8_t **ff, uint32_t **ng_occ, uint64_t **i_idx, asg64_v *b64, asg64_v *ub64,
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;
}
kt_for(uidx->str_b.n_thread, worker_ul_aln_path, &u_aux, u_aux.int_idx_n);
if(thread_path) {
kt_for(uidx->str_b.n_thread, worker_ul_aln_path, &u_aux, u_aux.int_idx_n);
} else {
kt_for(uidx->str_b.n_thread, worker_ul_aln_pair, &u_aux, u_aux.int_idx_n);
}
for (k = b64->n = a_n = 0; k < uidx->str_b.n_thread; k++) {
a_n += uidx->str_b.buf[k].res_dump.n;
kv_resize(uint64_t, *b64, a_n);
@@ -14229,9 +14388,11 @@ void u2g_hybrid_detan_iter(ul_resolve_t *uidx, usg_t *ng, uint32_t max_ext, uint
// prt_usg_t(uidx, ng, "ng0");
for (k = 0; k < clean_round; k++) {
ncut += ug_ext_strict(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64);
ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 48);
ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 16);
ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 48, 1);
// prt_usg_t(uidx, ng, "ng_python");
ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 16, 1);
ncut += ug_ext_free(ob, ub, uidx, ng, max_ext, &ff, &ng_occ, &i_idx, &b64, &ub64, 0, 0);
///renew bubble for ug_ext_strict
n_vtx = ng->n<<1;
REALLOC(ff, n_vtx); REALLOC(ng_occ, n_vtx); REALLOC(i_idx, n_vtx);
hic.cpp
+1 -2
View File
@@ -285,7 +285,6 @@ KRADIX_SORT_INIT(k_trans_qs, u_trans_t, k_trans_qs_key, member_size(u_trans_t, q
#define get_hit_euid(x, k) (((x).a.a[(k)].e<<1)>>(64 - (x).uID_bits))
#define get_hit_epos(x, k) ((x).a.a[(k)].e & (x).pos_mode)
typedef struct {
kvec_t(kvec_t_u64_warp) matrix;
uint64_t uID_shift, dis_mode;
@@ -17551,4 +17550,4 @@ void hic_benchmark(ma_ug_t *ug, asg_t* read_g)
hic_short_align_bench(asm_opt.hic_reads[0], asm_opt.hic_reads[1], output_file_name, ug_index);
free(output_file_name);
}
}
inter.cpp
+188 -2
View File
@@ -25,6 +25,9 @@ void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint6
int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t high_occ, void *km);
void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres,
int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off);
void ug_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres,
int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off, uint32_t is_hpc);
#define MG_SEED_IGNORE (1ULL<<41)
#define MG_SEED_TANDEM (1ULL<<42)
#define MG_SEED_KEPT (1ULL<<43)
@@ -248,8 +251,6 @@ KRADIX_SORT_INIT(eg_srt_x, eg_srt_t, eg_srt_x_key, member_size(eg_srt_t, x))
#define eg_srt_d_key(p) ((p).d)
KRADIX_SORT_INIT(eg_srt_d, eg_srt_t, eg_srt_d_key, member_size(eg_srt_t, d))
// shortest path
typedef struct {
// input
@@ -374,6 +375,28 @@ typedef struct { // data structure for each step in kt_pipeline()
} utepdat_t;
typedef struct { // global data structure for kt_pipeline()
uint64_t x, i;
} ha_mzl_srt_t;
typedef struct { size_t n, m; ha_mzl_srt_t *a; } ha_mzl_srt_v;
typedef struct { // global data structure for kt_pipeline()
ma_ug_t *ug;
asg_t *rg;
ha_ovec_buf_t **hab;
glchain_t *ll;
int32_t w, k;
///base-alignment
double bw_thres, diff_ec_ul;
int32_t max_n_chain;
ha_mzl_v idx_a;
asg64_v idx_n;
ha_mzl_srt_v srt_a;
int32_t is_HPC, bw, max_gap, chn_pen_gap, n_thread, is_cnt;
ul_idx_t udb;
} ug_trans_t;
void hc_glchain_destroy(glchain_t *b)
{
if (!b) return;
@@ -9184,6 +9207,126 @@ uint32_t ck_ul_alignment(ul_vec_t *x)
}
static void worker_for_trans_ovlp(void *data, long i, int tid) // callback for kt_for()
{
/**
ug_trans_t *s = (ug_trans_t*)data;
ha_ovec_buf_t *b = s->hab[tid];
glchain_t *bl = &(s->ll[tid]);
int64_t winLen = MIN((((double)THRESHOLD_MAX_SIZE)/s->diff_ec_ul), WINDOW), ton = 0;
uint32_t high_occ = 2, phase = 1, k;
asg64_v b0, b1, b2; window_list p; memset(&p, 0, sizeof(p));
overlap_region *aux_o = NULL;
char *seq = s->ug->u.a[i].s;
int64_t len = s->ug->u.a[i].len;
**/
// uint64_t align = 0;
// if(UL_INF.a[s->id+i].rlen != s->len[i]) {
// fprintf(stderr, "[M::%s] rid:%ld, s->len:%lu, UL_INF->rlen:%u\n", __func__, s->id+i, s->len[i], UL_INF.a[s->id+i].rlen);
// }
// void *km = s->buf?(s->buf[tid]?s->buf[tid]->km:NULL):NULL;
// if(s->id+i!=3046) return;
// if((s->id+i!=871) && (s->id+i!=963) && (s->id+i!=980)) return;
// if(s->id+i!=944) return;
// if(s->id+i != 35437) return;
// if((s->id+i != 7086) && (s->id+i != 51705) && (s->id+i != 266022) && (s->id+i != 353608)
// && (s->id+i != 399416) && (s->id+i != 403014) && (s->id+i != 420915) && (s->id+i != 603855)
// && (s->id+i != 680134) && (s->id+i != 766261) && (s->id+i != 794527)) {
// return;
// }
// char *as = NULL;
// asprintf(&as, "\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);
// push_vlog(&(overall_zdbg->a[s->id+i]), as); free(as); as = NULL;
// 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);
// fprintf(stderr, ">%.*s\n%.*s\n", (int32_t)UL_INF.nid.a[s->id+i].n, UL_INF.nid.a[s->id+i].a,
// (int32_t)s->len[i], s->seq[i]);
// if (memcmp(UL_INF.nid.a[s->id+i].a, "d0aab024-b3a7-40fb-83cc-22c3d6d951f8", UL_INF.nid.a[s->id+i].n-1)) return;
// fprintf(stderr, "[M::%s::] ==> len: %lu\n", __func__, s->len[i]);
// ha_get_ul_candidates_interface(b->abl, i, s->seq[i], s->len[i], s->opt->w, s->opt->k, s->uu, &b->olist, &b->olist_hp, &b->clist, s->opt->bw_thres,
// s->opt->max_n_chain, 1, NULL, &b->r_buf, &(b->tmp_region), NULL, &(b->sp), 1, NULL);
/**
ug_map_lchain(b->abl, (uint32_t)-1, seq, len, s->w, s->k, &(s->udb), NULL, NULL, s->bw_thres,
s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, s->is_HPC);
**/
/**
clear_Cigar_record(&b->cigar1);
clear_Round2_alignment(&b->round2);
// return;
// b->num_correct_base += overlap_statistics(&b->olist, NULL, 0);
// int fully_cov, abnormal;
// b->self_read.seq = s->seq[i]; b->self_read.length = s->len[i]; b->self_read.size = 0;
// correct_ul_overlap(&b->olist, s->uu, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon,
// &b->cigar1, &b->hap, &b->round2, &b->r_buf, &(b->tmp_region.w_list), 0, 1, &fully_cov, &abnormal, s->opt->diff_ec_ul, winLen, NULL);
// memset(&b->self_read, 0, sizeof(b->self_read));
ul_lalign(&b->olist, &b->clist, s->uu, s->uopt, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read,
&b->correct, &b->exz, &b->hap, &b->r_buf, aux_o, s->opt->diff_ec_ul, winLen, NULL, s->id+i, s->opt->k, &(s->sps[tid]), NULL);
// ul_lalign_old_ed(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read,
// &b->correct, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, 1, NULL);
ton = b->olist.length;//all alignments pass similary check
aux_o = gen_aux_ovlp(&b->olist);///must be here
gl_chain_flter(&b->olist, &b->correct, &(s->sps[tid]), bl, s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, &phase);
// fprintf(stderr, "[M::%s] rid::%ld, len::%lu, name::%.*s, phase::%u\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, phase);
if(phase && gen_shared_intervals(&b->olist, s->uu, s->uopt, winLen, &b->r_buf, &(bl->lo))) {
filter_topN(&b->olist, &(bl->lo), s->len[i], winLen, UL_TOPN, bl);
// update_shared_intervals(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b->r_buf, &(s->sps[tid]), s->len[i], winLen, &(bl->lo), s->id+i);
copy_asg_arr(b0, b->hap.snp_srt); copy_asg_arr(b1, s->sps[tid]); copy_asg_arr(b2, b->r_buf.a);
// update_shared_intervals(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b0, &b1, &b2, s->len[i], winLen, &(bl->lo), s->id+i);
update_sketch_trace(&b->olist, s->uu, s->uopt, NULL, &b->ovlp_read, &b0, &b1, &b2, s->len[i], winLen, &(bl->lo), s->id+i, MAX_LGAP(s->len[i]), s->opt->diff_ec_ul);
copy_asg_arr(b->hap.snp_srt, b0); copy_asg_arr(s->sps[tid], b1); copy_asg_arr(b->r_buf.a, b2);
ul_lalign(&b->olist, &b->clist, s->uu, s->uopt, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read,
&b->correct, &b->exz, &b->hap, &b->r_buf, aux_o, s->opt->diff_ec_ul, winLen, &(bl->lo), s->id+i, s->opt->k, &(s->sps[tid]), NULL);
// ul_lalign_old_ed(&b->olist, &b->clist, s->uu, s->seq[i], s->len[i], &b->self_read, &b->ovlp_read,
// &b->correct, &b->hap, &b->r_buf, s->opt->diff_ec_ul, winLen, 0, NULL);
///recover alignments
for (k = b->olist.length; k < ton; k++) {
b->olist.list[k].w_list.n = 0;
p.x_start = b->olist.list[k].x_pos_s;
p.x_end = b->olist.list[k].x_pos_e+1;
p.clen = b->olist.list[k].non_homopolymer_errors;
kv_push(window_list, b->olist.list[k].w_list, p);
b->olist.list[k].align_length = 0;
b->olist.list[k].overlapLen = b->olist.list[k].x_pos_e+1-b->olist.list[k].x_pos_s;
}
b->olist.length = ton;
} else {
for (k = 0; k < b->olist.length; k++) {
b->olist.list[k].w_list.n = 0;
p.x_start = b->olist.list[k].x_pos_s;
p.x_end = b->olist.list[k].x_pos_e+1;
p.clen = 0;
kv_push(window_list, b->olist.list[k].w_list, p);
b->olist.list[k].align_length = b->olist.list[k].overlapLen =
b->olist.list[k].x_pos_e+1-b->olist.list[k].x_pos_s;
b->olist.list[k].non_homopolymer_errors = 0;
}
}
// prt_overlap_region_alloc_ol(&(b->olist), s->id+i);
gl_chain(s->buf[tid], &(UL_INF.a[s->id+i]), &b->olist, &(b->clist.chainDP), &b->hap, &(s->sps[tid]), bl, &(s->gdp[tid]), s->uu, s->opt->diff_ec_ul, winLen, s->len[i], s->uopt, s->id+i, tid, NULL);
if(UL_INF.a[s->id+i].dd == 1 || UL_INF.a[s->id+i].dd == 2) {
b->num_correct_base++;
}
if(UL_INF.a[s->id+i].dd != 3) {
free(s->seq[i]); s->seq[i] = NULL;
}
s->hab[tid]->num_read_base++;
**/
}
static void worker_for_ul_recorrect_alignment(void *data, long i, int tid) // callback for kt_for()
{
utepdat_t *s = (utepdat_t*)data;
@@ -14891,6 +15034,49 @@ void clear_all_ul_t(all_ul_t *x)
free(x->ridx.idx.a); free(x->ridx.occ.a); memset(&(x->ridx), 0, sizeof((x->ridx)));
}
void init_ug_trans_t(ug_trans_t *opt, int32_t is_HPC, int32_t k, int32_t w, int32_t max_n_chain,
double bw_thres, double diff_ec_ul, int32_t n_thread, ma_ug_t *ug, asg_t *sg)
{
int64_t i;
memset(opt, 0, sizeof((*opt)));
opt->k = k;
opt->w = w;
opt->is_HPC = is_HPC;
opt->max_n_chain = max_n_chain;
opt->bw_thres = bw_thres;
opt->diff_ec_ul = diff_ec_ul;
opt->bw = 10000;///2000 in minigraph
opt->max_gap = 500000;///5000 in minigraph
opt->ug = ug;
opt->rg = sg;
opt->n_thread = ((n_thread>=1)?n_thread:1);
CALLOC(opt->hab, opt->n_thread);
CALLOC(opt->ll, opt->n_thread);
for (i = 0; i < opt->n_thread; ++i) {
opt->hab[i] = ha_ovec_init(0, 0, 1);
}
opt->idx_n.n = opt->idx_n.m = ug->u.n;
MALLOC(opt->idx_n.a, opt->idx_n.n);
opt->udb.ug = ug;
}
void trans_base_count(ug_trans_t *p)
{
p->is_cnt = 1; ha_flt_tab = NULL;
kt_for(p->n_thread, worker_for_trans_ovlp, p, p->ug->u.n);
}
void trans_base_infer(ma_ug_t *ug, asg_t *sg)
{
ug_trans_t sl;
init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov);
init_ug_trans_t(&sl, 0, asm_opt.trans_mer_length, asm_opt.trans_win, asm_opt.max_n_chain, 1.0-asm_opt.trans_base_rate, 1.0-asm_opt.trans_base_rate, asm_opt.thread_num, ug, sg);
}
ma_ug_t *ul_realignment(const ug_opt_t *uopt, asg_t *sg, uint32_t double_check_cache)