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

Merge pull request #429 from chhylp123/hifiasm_dev_debug

Browse Source 
avoid misassemblies; better polyploidy graph
This commit is contained in:
chhylp123
2023-03-22 12:24:31 -04:00
committed by GitHub
parent b763e1ff76 86b7dd424a
commit 00ad7458c3
12 changed files with 1931 additions and 79 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.cpp
+1 -1
View File
@@ -285,7 +285,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->min_path_drop_rate = 0.2;
asm_opt->max_path_drop_rate = 0.6;
asm_opt->hifi_pst_join = 1;
asm_opt->ul_pst_join = 0;
asm_opt->ul_pst_join = 1;
}
void destory_enzyme(enzyme* f)
CommandLines.h
+1 -1
View File
@@ -5,7 +5,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.19.2-r560"
#define HA_VERSION "0.19.2-r572"
#define VERBOSE 0
Correct.cpp
+70 -15
View File
@@ -15216,8 +15216,10 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u
else ibeg = ov->qn;
iend = ov->tn;
assert(iend>=ibeg+1);
// fprintf(stderr, "\n***[M::%s::rid->%lu] utg%.6dl(%c), s::%u, e::%u, z::[%u, %u), ibeg::%ld, iend::%ld, ch_n::%ld\n",
// __func__, rid, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], ov->qs, ov->qe, z->x_pos_s, z->x_pos_e+1, ibeg, iend, ch_n);
// if(z->x_id == 57 && z->y_id == 2175) {
// fprintf(stderr, "\n***[M::%s::rid->%lu] utg%.6dl(%c), s::%u, e::%u, z::[%u, %u), ibeg::%ld, iend::%ld, ch_n::%ld\n",
// __func__, rid, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], ov->qs, ov->qe, z->x_pos_s, z->x_pos_e+1, ibeg, iend, ch_n);
// }
for (l = ibeg, i = ibeg + 1; i <= iend; i++) {
q[0] = q[1] = t[0] = t[1] = mode = -1; is_done = 0;
if(l >= 0) {
@@ -15243,14 +15245,14 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u
}
if(mode == 1 || mode == 2) adjust_ext_offset(&(q[0]), &(q[1]), &(t[0]), &(t[1]), ql, tl, 0, mode);
// if(aux_o->x_id == 29033 && aux_o->y_id == 21307) {
// if(z->x_id == 57 && z->y_id == 2175) {
// 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_SIN_L, MAX_SIN_E, FORCE_SIN_L, -1, aux_o);
// if(aux_o->x_id == 29033 && aux_o->y_id == 21307) {
// if(z->x_id == 57 && z->y_id == 2175) {
// fprintf(stderr, "-is_done::%ld[M::%s::] utg%.6dl(%c), q::[%ld, %ld), t::[%ld, %ld), mode::%ld\n",
// is_done, __func__, (int32_t)z->y_id+1, "+-"[z->y_pos_strand], q[0], q[1], t[0], t[1], mode);
// }
@@ -19123,11 +19125,40 @@ uint64_t gen_sub_ov_adv(const ul_idx_t *udb, overlap_region* o, double o_rate, u
int64_t return_t_chain(overlap_region *z, Candidates_list *cl)
{
int64_t i, cn = cl->length, scn; uint64_t pid; k_mer_hit *ca;
// if(z->x_id == 57 && z->y_id == 2175) {
// fprintf(stderr, "\n-0-[M::%s]\tutg%.6ul\tx::[%u,\t%u)\t%c\tutg%.6ul\ty::[%u,\t%u)\n",
// __func__, z->x_id+1, z->x_pos_s, z->x_pos_e+1,
// "+-"[z->y_pos_strand], z->y_id+1, z->y_pos_s, z->y_pos_e+1);
// i = z->shared_seed; pid = cl->list[i].readID;
// for (; i < cn && cl->list[i].readID == pid && cl->list[i].readID != ((uint32_t)(0x7fffffff)); i++) {
// fprintf(stderr, "i::%ld[M::%s]\treadID::%u\tself_offset::%u\toffset::%u\t%c\n",
// i, __func__, cl->list[i].readID, cl->list[i].self_offset, cl->list[i].offset,
// "+-"[cl->list[i].strand]);
// }
// }
i = z->shared_seed; pid = cl->list[i].readID;
for (; i < cn && cl->list[i].readID == pid && cl->list[i].readID != ((uint32_t)(0x7fffffff)); i++);
scn = i - z->shared_seed; ca = cl->list+z->shared_seed;
i = lchain_refine(ca, scn, ca, &(cl->chainDP), 50, 5000, 512, 16); cn = i;
for (; i < scn; i++) ca[i].readID = ((uint32_t)(0x7fffffff));
// if(z->x_id == 57 && z->y_id == 2175) {
// fprintf(stderr, "\n-1-[M::%s]\tutg%.6ul\tx::[%u,\t%u)\t%c\tutg%.6ul\ty::[%u,\t%u)\n",
// __func__, z->x_id+1, z->x_pos_s, z->x_pos_e+1,
// "+-"[z->y_pos_strand], z->y_id+1, z->y_pos_s, z->y_pos_e+1);
// i = z->shared_seed; pid = cl->list[i].readID;
// for (; i < cn && cl->list[i].readID == pid && cl->list[i].readID != ((uint32_t)(0x7fffffff)); i++) {
// fprintf(stderr, "i::%ld[M::%s]\treadID::%u\tself_offset::%u\toffset::%u\t%c\n",
// i, __func__, cl->list[i].readID, cl->list[i].self_offset, cl->list[i].offset,
// "+-"[cl->list[i].strand]);
// }
// }
return cn;
}
@@ -21017,6 +21048,20 @@ bit_extz_t *exz, double e_rate, int64_t w_l, uint64_t ql, uint64_t rid, uint64_t
tl = udb->ug->u.a[id].len;
for (i = ch_idx; i < cl->length && cl->list[i].readID == cl->list[ch_idx].readID; i++); ch_n = i-ch_idx;
// if(z->x_id == 57 && z->y_id == 2175) {
// fprintf(stderr, "\n-*-[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\tch_idx::%ld\tch_n::%ld\n",
// __func__,
// z->x_id+1, "lc"[udb->ug->u.a[z->x_id].circ], udb->ug->u.a[z->x_id].len, z->x_pos_s, z->x_pos_e+1,
// "+-"[z->y_pos_strand],
// z->y_id+1, "lc"[udb->ug->u.a[z->y_id].circ], udb->ug->u.a[z->y_id].len, z->y_pos_s, z->y_pos_e+1,
// ch_idx, ch_n);
// for (i = ch_idx; i < cl->length && cl->list[i].readID == cl->list[ch_idx].readID; i++) {
// fprintf(stderr, "i::%ld[M::%s]\treadID::%u\tself_offset::%u\toffset::%u\t%c\n",
// i, __func__, cl->list[i].readID, cl->list[i].self_offset, cl->list[i].offset,
// "+-"[cl->list[i].strand]);
// }
// }
// on = fusion_chain_ovlp(z, ch_a, ch_n, ov, on, wl, ql, tl);
aux_o->w_list.n = aux_o->w_list.c.n = 0;
aux_o->y_id = z->y_id; aux_o->y_pos_strand = z->y_pos_strand;
@@ -21033,8 +21078,8 @@ bit_extz_t *exz, double e_rate, int64_t w_l, uint64_t ql, uint64_t rid, uint64_t
int64_t aux_n = aux_o->w_list.n;
// if(z->x_id == 77960 && z->y_id == 27340) {
// fprintf(stderr, "\n[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\n",
// if(z->x_id == 57 && z->y_id == 2175) {
// fprintf(stderr, "\n-0-[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\n",
// __func__,
// z->x_id+1, "lc"[udb->ug->u.a[z->x_id].circ], udb->ug->u.a[z->x_id].len, z->x_pos_s, z->x_pos_e+1,
// "+-"[z->y_pos_strand],
@@ -21057,8 +21102,8 @@ bit_extz_t *exz, double e_rate, int64_t w_l, uint64_t ql, uint64_t rid, uint64_t
rechain_aln(z, cl, aux_o, i, w_l, udb, NULL, NULL, qstr, tu, exz, e_rate, ql, tl, khit, rid);
}
// if(z->x_id == 77960 && z->y_id == 27340) {
// fprintf(stderr, "\n[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\n",
// if(z->x_id == 57 && z->y_id == 2175) {
// fprintf(stderr, "\n-1-[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\n",
// __func__,
// z->x_id+1, "lc"[udb->ug->u.a[z->x_id].circ], udb->ug->u.a[z->x_id].len, z->x_pos_s, z->x_pos_e+1,
// "+-"[z->y_pos_strand],
@@ -21085,6 +21130,7 @@ bit_extz_t *exz, double e_rate, int64_t w_l, uint64_t ql, uint64_t rid, uint64_t
///update z by aux_o
update_overlap_region(z, aux_o, ql, tl);
assert((z->x_pos_e>=z->x_pos_s) && (z->y_pos_e>=z->y_pos_s));
int64_t zwn = z->w_list.n, zerr = 0, zlen = z->x_pos_e+1-z->x_pos_s;
for (i = 0; i < zwn; i++) {
@@ -21097,8 +21143,8 @@ bit_extz_t *exz, double e_rate, int64_t w_l, uint64_t ql, uint64_t rid, uint64_t
}
z->non_homopolymer_errors = zerr;
// if(z->x_id == 77960 && z->y_id == 27340) {
// fprintf(stderr, "\n[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\n",
// if(z->x_id == 57 && z->y_id == 2175) {
// fprintf(stderr, "\n-2-[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\n",
// __func__,
// z->x_id+1, "lc"[udb->ug->u.a[z->x_id].circ], udb->ug->u.a[z->x_id].len, z->x_pos_s, z->x_pos_e+1,
// "+-"[z->y_pos_strand],
@@ -21114,6 +21160,15 @@ bit_extz_t *exz, double e_rate, int64_t w_l, uint64_t ql, uint64_t rid, uint64_t
// fprintf(stderr, "\n");
// }
// if(z->x_id == 57 && z->y_id == 2175) {
// fprintf(stderr, "\n-3-[M::%s]\tutg%.6u%c\txl::%u\tx::[%u,\t%u)\t%c\tutg%.6u%c\tyl::%u\ty::[%u,\t%u)\tzerr::%ld\tzlen::%ld\te_rate::%f\n",
// __func__,
// z->x_id+1, "lc"[udb->ug->u.a[z->x_id].circ], udb->ug->u.a[z->x_id].len, z->x_pos_s, z->x_pos_e+1,
// "+-"[z->y_pos_strand],
// z->y_id+1, "lc"[udb->ug->u.a[z->y_id].circ], udb->ug->u.a[z->y_id].len, z->y_pos_s, z->y_pos_e+1,
// zerr, zlen, e_rate);
// }
if(zerr >= zlen) return 0;
if(zerr <= 0 && zlen > 0) return 1;
if(zerr > (zlen*e_rate)) return 0;
@@ -21319,7 +21374,7 @@ void ug_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur
rr = gen_extend_err_exz(z, uref, NULL, NULL, in/**qu->seq**/, tu->seq, exz, NULL/**v_idx?v_idx->a.a:NULL**/, w.window_length, -1, err, (e_max+0.000001), &re);
z->is_match = 0;///must be here;
// if(z->x_id == 77960 && z->y_id == 27340) {
// if(z->x_id == 57 && z->y_id == 2175) {
// fprintf(stderr, "+utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\talign_length::%u\terr::%f\trr::%f\twn::%u\tk::%lu\ti::%lu\n",
// z->x_id+1, "lc"[uref->ug->u.a[z->x_id].circ], uref->ug->u.a[z->x_id].len,
// z->x_pos_s, z->x_pos_e+1, "+-"[z->y_pos_strand],
@@ -21339,7 +21394,7 @@ void ug_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur
}
ol->length = k;
// if(sid == 77960) {
// if(sid == 57) {
// fprintf(stderr, "+utg%.6ld%c\tol->length::%lu\n", sid+1, "lc"[uref->ug->u.a[sid].circ], ol->length);
// }
if(ol->length <= 0) return;
@@ -21351,7 +21406,7 @@ void ug_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur
}
for (i = k = 0; i < ol->length; i++) {
z = &(ol->list[i]);
// if(z->x_id == 77960 && z->y_id == 27340) {
// if(z->x_id == 57 && z->y_id == 2175) {
// fprintf(stderr, "-1-utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\talign_length::%u\terr::%f\twn::%u\tk::%lu\ti::%lu\n",
// z->x_id+1, "lc"[uref->ug->u.a[z->x_id].circ], uref->ug->u.a[z->x_id].len,
// z->x_pos_s, z->x_pos_e+1, "+-"[z->y_pos_strand],
@@ -21364,7 +21419,7 @@ void ug_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur
continue;
}
// if(z->x_id == 77960 && z->y_id == 27340) {
// if(z->x_id == 57 && z->y_id == 2175) {
// fprintf(stderr, "-2-utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\talign_length::%u\terr::%f\twn::%u\tk::%lu\ti::%lu\n",
// z->x_id+1, "lc"[uref->ug->u.a[z->x_id].circ], uref->ug->u.a[z->x_id].len,
// z->x_pos_s, z->x_pos_e+1, "+-"[z->y_pos_strand],
@@ -21382,7 +21437,7 @@ void ug_lalign(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *ur
}
ol->length = k;
// if(sid == 77960) {
// if(sid == 57) {
// fprintf(stderr, "-utg%.6ld%c\tol->length::%lu\n", sid+1, "lc"[uref->ug->u.a[sid].circ], ol->length);
// }
}
Overlaps.cpp
+996 -5
View File
File diff suppressed because it is too large Load Diff
Overlaps.h
+23
View File
@@ -1211,4 +1211,27 @@ int asg_arc_del_trans_ul(asg_t *g, int fuzz);
#define JUNK_COV 5
#define DISCARD_RATE 0.8
typedef struct {
uint32_t n, m, a;
} mmhap_status_t;
typedef struct {
kvec_t(mmhap_status_t) h;
kvec_t(uint32_t) a;
} mmhap_t;
typedef struct { // global data structure for kt_pipeline()
ma_ug_t *ug;
asg_t *rg;
uint64_t *idx;
asg64_v cov;
uint64_t hom_min, hom_max, hom_cov, het_cov;
} ug_rid_cov_t;
ug_rid_cov_t* gen_ug_rid_cov_t(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc *src);
void destory_ug_rid_cov_t(ug_rid_cov_t *p);
uint32_t append_cov_line_ug_rid_cov_t(uint64_t uid, uint64_t *qcc, u_trans_t *p, ug_rid_cov_t *idx, uint64_t hom_cut, double cut_rate);
uint64_t infer_mmhap_copy(ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc *src, uint8_t *ff, uint64_t uid, uint64_t het_cov, uint64_t n_hap);
uint64_t trans_sec_cut0(kv_u_trans_t *ta, asg64_v *srt, uint32_t id, double sec_rate, uint64_t bd, ma_ug_t *ug);
#endif
anchor.cpp
+16 -2
View File
@@ -1620,7 +1620,7 @@ inline uint64_t special_lchain(Candidates_list* cl, overlap_region_alloc* ol, ui
double chn_pen_skip, double bw_rate, int64_t quick_check, uint32_t gen_off, double mcopy_rate, uint32_t mcopy_khit_cut,
st_mt_t *sp, uint64_t *si, uint64_t m, uint64_t l, uint64_t k)
{
uint64_t z, s, e, yid, ol0 = ol->length, ol1;
uint64_t z, s, e, yid, ol0 = ol->length, ol1, m0 = m;
if(l >= k) return m;
yid = cl->list[l].readID;
for (z = l; z < k && cl->list[z].strand == cl->list[l].strand; z++);
@@ -1678,7 +1678,21 @@ inline uint64_t special_lchain(Candidates_list* cl, overlap_region_alloc* ol, ui
}
s++;
}
ol->length = s;
if(s != ol->length) {
ol->length = s;
for (z = ol0; z < ol->length; z++) {
s = ol->list[z].non_homopolymer_errors;
pi = cl->list[s].readID;
ol->list[z].non_homopolymer_errors = m0;
for (; s < m && cl->list[s].readID == pi; s++) {
cl->list[m0] = cl->list[s];
cl->list[m0].readID = z; m0++;
}
}
m = m0;
}
}
return m;
}
gfa_ut.cpp
+27 -14
View File
@@ -1262,7 +1262,7 @@ uint32_t minLen, asg64_v *t)
}
void asg_arc_cut_length(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_trio,
uint32_t is_topo, uint32_t min_diff, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len)
uint32_t is_topo, uint32_t min_diff, uint32_t min_ou, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len)
{
asg64_v tx = {0,0,0}, *b = NULL;
uint32_t i, k, v, w, n_vtx = g->n_seq<<1, nv, nw, kv, kw, trioF = (uint32_t)-1, ntrioF = (uint32_t)-1, ol_max, ou_max, to_del, cnt = 0, mm_ol, mm_ou;
@@ -1341,7 +1341,7 @@ uint32_t is_topo, uint32_t min_diff, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *
if (kv < 1) continue;
if (kv >= 2) {
if (mm_ol > ol_max*len_rat) continue;
if (is_ou && mm_ou > ou_max*ou_rat) continue;
if (is_ou && mm_ou > ou_max*ou_rat && mm_ou > min_ou) continue;
if ((mm_ol + min_diff) > ol_max) continue;
}
@@ -1356,7 +1356,7 @@ uint32_t is_topo, uint32_t min_diff, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *
if (kw < 1) continue;
if (kw >= 2) {
if (mm_ol > ol_max*len_rat) continue;
if (is_ou && mm_ou > ou_max*ou_rat) continue;
if (is_ou && mm_ou > ou_max*ou_rat && mm_ou > min_ou) continue;
if ((mm_ol + min_diff) > ol_max) continue;
}
@@ -2798,7 +2798,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
// prt_specfic_sge(sg, 10531, 10519, "--2--");
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 1);
asg_arc_cut_length(sg, &bu, max_tip, drop, ou_drop_rate, is_ou, is_trio, 1, min_diff, NULL, NULL, NULL);
asg_arc_cut_length(sg, &bu, max_tip, drop, ou_drop_rate, is_ou, is_trio, 1, min_diff, 1, NULL, NULL, NULL);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL);
// prt_specfic_sge(sg, 10531, 10519, "--3--");
@@ -2851,13 +2851,18 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
if(!is_ou) {
///asg_arc_del_triangular_directly might be unnecessary
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0);
asg_arc_cut_length(sg, &bu, max_tip, HARD_ORTHOLOGY_DROP/**min_ovlp_drop_ratio**/, ou_drop_rate, is_ou, 0/**is_trio**/, is_ou?1:0, min_diff, rev, rI, NULL);
asg_arc_cut_length(sg, &bu, max_tip, HARD_ORTHOLOGY_DROP/**min_ovlp_drop_ratio**/, ou_drop_rate, is_ou, 0/**is_trio**/, is_ou?1:0, min_diff, 1, rev, rI, NULL);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL);
// print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty7.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0);
asg_arc_cut_length(sg, &bu, max_tip, min_ovlp_drop_ratio, ou_drop_rate, is_ou, 0/**is_trio**/, is_ou?1:0, min_diff, rev, rI, &l_drop);
asg_arc_cut_length(sg, &bu, max_tip, min_ovlp_drop_ratio, ou_drop_rate, is_ou, 0/**is_trio**/, is_ou?1:0, min_diff, 1, rev, rI, &l_drop);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL);
} else {
min_diff = step_diff; l_drop = 6000;
asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0);
asg_arc_cut_length(sg, &bu, max_tip, 0.3, 0.9, is_ou, is_trio, 1, min_diff, 8, NULL, NULL, &l_drop);
asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL);
}
@@ -11923,16 +11928,18 @@ uint32_t ava_pass_unique_bridge(uint64_t *idx, uint64_t *integer_seq, uint64_t s
return 1;
}
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, double max_ext_rate, uint64_t ext_up)
{
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, nkeep = 0, ncut = 0; usg_arc_t *av;
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);
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;
nkeep += g->a[integer_seq[k]>>1].occ;
}
f[integer_seq[s]] = f[integer_seq[e]^1] = 1;
}
nkeep = nkeep*max_ext_rate;
///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++) {
@@ -11964,7 +11971,10 @@ uint32_t ava_pass_unique_bridge_tips(usg_t *g, asg64_v *b64, uint64_t g_s, uint6
}
}
while (b64->n > bn && n_ext < max_ext) {
ncut = MAX(nkeep, max_ext);
if(ncut > ext_up) ncut = ext_up;
if(ncut < max_ext) ncut = max_ext;
while (b64->n > bn && n_ext < ncut) {
v = b64->a[--b64->n]; if(f[v]) continue;
av = usg_arc_a(g, v^1); nv = usg_arc_n(g, v^1);
for (i = kv = 0; i < nv && kv < 1; i++) {
@@ -12014,7 +12024,7 @@ uint32_t ava_pass_unique_bridge_tips(usg_t *g, asg64_v *b64, uint64_t g_s, uint6
}
}
}
if((n_ext < ncut) && (n_ext > max_ext - 1)) n_ext = max_ext - 1;
return n_ext;
}
@@ -12954,7 +12964,7 @@ uint8_t *ff, uint32_t *ng_occ, uint32_t max_ext)
n_clus = 0;
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)) {
tip_l = ava_pass_unique_bridge_tips(ng, b64, i, k, int_a, ff, max_ext);
tip_l = ava_pass_unique_bridge_tips(ng, b64, i, k, int_a, ff, max_ext, 0.03, 16);
if(tip_l < max_ext) {///no long tip
if(/**(!tip_l) || (**/ava_pass_unique_bridge_cov(uidx, ng, b64, i, k, int_a, tip_l)) {
for (z = i; z < k; z++) {
@@ -15934,11 +15944,14 @@ void u2g_clean(ul_resolve_t *uidx, ulg_opt_t *ulopt, uint32_t keep_raw_utg, uint
// char sb[1000];
// output_integer_graph(uidx, iug, "ig_h0", 0);
// fprintf(stderr, "\n[M::%s::] max_path_drop_ratio::%f, min_path_drop_ratio::%f\n", __func__, ulopt->max_path_drop_ratio, ulopt->min_path_drop_ratio);
// fprintf(stderr, "\n[M::%s::] max_path_drop_ratio::%f, min_path_drop_ratio::%f, max_tip_hifi::%ld, max_tip::%ld\n",
// __func__, ulopt->max_path_drop_ratio, ulopt->min_path_drop_ratio,
// ulopt->max_tip_hifi, ulopt->max_tip);
for (ss = 0; ss < 2; ss++) {
for (i = 0, drop = ulopt->min_path_drop_ratio; i < ulopt->clean_round; i++, drop += step) {
if(drop > ulopt->max_path_drop_ratio) drop = ulopt->max_path_drop_ratio;
// fprintf(stderr, "[M::%s::] Starting round-%ld, drop::%f\n", __func__, i, drop);
// fprintf(stderr, "[M::%s::] Starting round-%ld, drop::%f, max_tip_hifi::%ld\n",
// __func__, i, drop, ulopt->max_tip_hifi);
cnt = 1; topo_level = 2; mm_tip = ulopt->max_tip;
while (cnt) {
cnt = 0;
gfa_ut.h
+1 -1
View File
@@ -22,7 +22,7 @@ void asg_iterative_semi_circ(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_
void asg_arc_cut_chimeric(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t ou_thres);
void asg_arc_cut_inexact(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, int32_t max_ext, uint32_t is_ou, uint32_t is_trio, uint32_t min_diff, float ou_rat/**, asg64_v *dbg**/);
void asg_arc_cut_length(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_trio,
uint32_t is_topo, uint32_t min_diff, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len);
uint32_t is_topo, uint32_t min_diff, uint32_t min_ou, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len);
void asg_arc_cut_bub_links(asg_t *g, asg64_v *in, float len_rat, float sec_len_rat, float ou_rat, uint32_t is_ou, uint64_t check_dist, ma_hit_t_alloc *rev, R_to_U* rI, int32_t max_ext);
void asg_arc_cut_complex_bub_links(asg_t *g, asg64_v *in, float len_rat, float ou_rat, uint32_t is_ou, bub_label_t *b_mask_t);
uint32_t asg_cut_large_indel(asg_t *g, asg64_v *in, int32_t max_ext, float ou_rat, uint32_t is_ou, uint32_t min_diff);
hic.cpp
+276 -3
View File
@@ -17430,8 +17430,280 @@ int hic_short_align_poy(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx,
return 1;
}
mmhap_t* gen_mmhap_t(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc *src)
{
uint64_t k, z, hom_cov, het_cov, s, *bs = NULL; uint8_t *ff; mmhap_t *p;
if(asm_opt.hom_global_coverage_set) {
hom_cov = asm_opt.hom_global_coverage;
} else {
hom_cov = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE);
}
het_cov = hom_cov/asm_opt.polyploidy;
CALLOC(ff, rg->n_seq); CALLOC(p, 1); CALLOC(bs, asm_opt.polyploidy+1);
p->h.n = p->h.m = ug->u.n; CALLOC(p->h.a, p->h.n);
for (k = p->a.n = s = 0; k < ug->u.n; k++) {
p->h.a[k].a = p->a.n;
p->h.a[k].n = 0;
p->h.a[k].m = infer_mmhap_copy(ug, rg, src, ff, k, het_cov, asm_opt.polyploidy);
if(p->h.a[k].m == (uint64_t)asm_opt.polyploidy) {
p->h.a[k].n = p->h.a[k].m; s = 1;
}
p->a.n += p->h.a[k].m; bs[p->h.a[k].m] += ug->g->seq[k].len;
}
free(ff);
p->a.m = p->a.n; MALLOC(p->a.a, p->a.n); memset(p->a.a, -1, sizeof((*(p->a.a)))*p->a.n);
if(s) {
for (k = p->a.n = 0; k < ug->u.n; k++) {
if(p->h.a[k].n == p->h.a[k].m) {
for (z = 0; z < p->h.a[k].m; z++) p->a.a[p->h.a[k].a+z] = z;
}
}
}
void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, uint32_t is_poy, kvec_pe_hit **rhits)
for (k = 1; k <= (uint64_t)asm_opt.polyploidy; k++) {
fprintf(stderr, "[M::stat] # %lu-copy bases: %lu\n", k, bs[k]);
}
free(bs);
return p;
}
bubble_type *gen_mmhap_bub(ma_ug_t* ug, uint8_t *r_het_flag, kv_u_trans_t *ref, mmhap_t *hh)
{
uint64_t k;
bubble_type *p; CALLOC(p, 1); p->n_round = asm_opt.n_weight; p->round_id = 0;
identify_bubbles(ug, p, r_het_flag, ref);
for (k = 0; k < ug->g->n_seq; k++) {
if(IF_BUB(k, (*p))) continue;
if(IF_HOM(k, (*p))) {
if(hh->h.a[k].n < hh->h.a[k].m) p->index[k] = p->num.n;
} else if(IF_HET(k, (*p))) {
if(hh->h.a[k].n >= hh->h.a[k].m) p->index[k] = (uint32_t)-1;
}
}
return p;
}
void purge_phase_0(ha_ug_index* idx, hc_links *link, bubble_type *bub, ps_t *s, kvec_pe_hit* hits, kv_u_trans_t *k_trans, uint8_t *del, mmhap_t *hh)
{
k_trans->idx.n = k_trans->n = 0; hits->idx.n = 0;
for (bub->round_id = 0; bub->round_id < bub->n_round; bub->round_id++) {
renew_kv_u_trans(k_trans, link, hits, &(idx->t_ch->k_trans), idx, bub, s->s, NULL, 1/**0**/);
mc_solve(NULL, NULL, k_trans, idx->ug, idx->read_g, 0.8, R_INF.trio_flag,
(bub->round_id==0?1:0), s->s, 1, bub, &(idx->t_ch->k_trans), 0,
/**(((bub.round_id+1) == bub.n_round)?1:0)**/0);
/*******************************for debug************************************/
// label_unitigs_sm(s->s, NULL, idx->ug);
}
uint64_t l[2], k, len; int64_t p; ma_ug_t *ug = idx->ug;
for (k = l[0] = l[1] = 0; k < ug->g->n_seq; k++) {
if((del[k] == 1) || (s->s[k] == 0)) continue;
if(s->s[k] > 0) l[0] += ug->g->seq[k].len;
else l[1] += ug->g->seq[k].len;
}
p = (l[0]>=l[1])?1:-1;
for (k = len = 0; k < ug->g->n_seq; k++) {
if((del[k] == 1) || (s->s[k] == 0)) {
if((del[k] == 2) ||
((hh->h.a[k].n == hh->h.a[k].m) && (hh->h.a[k].m == ((uint64_t)asm_opt.polyploidy)))) {
len += ug->g->seq[k].len;
}
continue;
}
if(s->s[k] == p) {
del[k] = 2; len += ug->g->seq[k].len;
} else {
del[k] = 1; bub->index[k] = (uint32_t)-1;
}
s->s[k] = 0;///reset
}
fprintf(stderr, "[M::%s::stat] # remaining bases: %lu\n", __func__, len);
}
void exchange_kv_u_trans_t(kv_u_trans_t *a, kv_u_trans_t *b)
{
uint64_t k, *ua; u_trans_t *u;
k = a->n; a->n = b->n; b->n = k;
k = a->m; a->m = b->m; b->m = k;
u = a->a; a->a = b->a; b->a = u;
k = a->idx.n; a->idx.n = b->idx.n; b->idx.n = k;
k = a->idx.m; a->idx.m = b->idx.m; b->idx.m = k;
ua = a->idx.a; a->idx.a = b->idx.a; b->idx.a = ua;
}
uint64_t cal_ave_ovlp(u_trans_t *a, uint64_t an, double top)
{
if(!an) return 0;
uint64_t k, len, cut, occ, tot;
for (k = len = 0; k < an; k++) {
len += a[k].qe - a[k].qs;
}
cut = len - (len*top);
for (k = occ = tot = 0; k < an; k++) {
if(a[k].qe - a[k].qs < cut) continue;
occ++; tot += a[k].qe - a[k].qs;
}
if(!occ) {
occ = an; tot = len;
}
return tot/occ;
}
void clean_trans_ovlp(bubble_type *bub, kv_u_trans_t *des, kv_u_trans_t *src, asg64_v *srt, ma_ug_t *ug)
{
uint64_t k, st, i, m = 0, ncut;
for (k = des->n = 0; k < src->n; k++) {
if(IF_HOM(src->a[k].qn, (*bub))) continue;
if(IF_HOM(src->a[k].tn, (*bub))) continue;
kv_push(u_trans_t, *des, src->a[k]);
}
kv_resize(uint64_t, des->idx, src->idx.n); des->idx.n = src->idx.n;
memset(des->idx.a, 0, des->idx.n*sizeof((*(des->idx.a))));
for (st = 0, i = 1; i <= des->n; ++i) {
if (i == des->n || des->a[i].qn != des->a[st].qn) {
des->idx.a[des->a[st].qn] = (((uint64_t)st)<<32)|(i-st); m++;
st = i;
}
}
if(srt && m) {
ncut = 0; kv_resize(uint64_t, *srt, m);
for (k = srt->n = 0; k < des->idx.n; k++) {
if(!((uint32_t)(des->idx.a[k]))) continue;
m = cal_ave_ovlp(des->a+(des->idx.a[k]>>32), ((uint32_t)(des->idx.a[k])), 0.9);
m = ((uint64_t)-1)-m; m <<= 32; m += k;
kv_push(uint64_t, *srt, m);
}
radix_sort_b64(srt->a, srt->a + srt->n);
for (k = 0; k < srt->n; k++) {
ncut += trans_sec_cut0(des, srt, (uint32_t)(srt->a[k]), 0.2, 256, ug);
}
if(ncut) {///renew idx
for (i = k = 0; i < des->n; i++) {
if(des->a[i].del) continue;
des->a[k++] = des->a[i];
}
des->n = k;
memset(des->idx.a, 0, des->idx.n*sizeof((*(des->idx.a))));
for (st = 0, i = 1; i <= des->n; ++i) {
if (i == des->n || des->a[i].qn != des->a[st].qn) {
des->idx.a[des->a[st].qn] = (((uint64_t)st)<<32)|(i-st); m++;
st = i;
}
}
}
}
}
void purge_phase(ha_ug_index* idx, mmhap_t *hh, uint64_t hapid, uint64_t max_round, hc_links *link, kvec_pe_hit* hits,
bubble_type *bub, ps_t *s, kv_u_trans_t *k_trans, uint8_t *del, uint32_t *bidx, kv_u_trans_t *ref, asg64_v *srt)
{
uint64_t k, len; uint32_t *bm;
ma_ug_t *ug = idx->ug;
if(max_round <= 0) {
for (k = 0; k < ug->g->n_seq; k++) {
if(hh->h.a[k].n >= hh->h.a[k].m) continue;
hh->a.a[hh->h.a[k].a+hh->h.a[k].n] = hapid;
hh->h.a[k].n++;
}
return;
}
for (k = 0; k < ug->g->n_seq; k++) {
del[k] = 0; s->s[k] = 0;
if(hh->h.a[k].n >= hh->h.a[k].m) {///all haplotypes have been set
bub->index[k] = (uint32_t)-1; del[k] = 1;
} else {
if(IF_HOM(k, (*bub))) bub->index[k] = bub->num.n;
}
bidx[k] = bub->index[k];
}
bm = bub->index; bub->index = bidx; bidx = bm;
kv_resize(u_trans_t, *ref, idx->t_ch->k_trans.n); ref->n = idx->t_ch->k_trans.n;
memcpy(ref->a, idx->t_ch->k_trans.a, ref->n*sizeof((*(ref->a))));
kv_resize(uint64_t, ref->idx, idx->t_ch->k_trans.idx.n); ref->idx.n = idx->t_ch->k_trans.idx.n;
memcpy(ref->idx.a, idx->t_ch->k_trans.idx.a, ref->idx.n*sizeof((*(ref->idx.a))));
exchange_kv_u_trans_t(ref, &(idx->t_ch->k_trans));
for (k = 0; k < max_round; k++) {
clean_trans_ovlp(bub, &(idx->t_ch->k_trans), ref, ((k+1)<max_round)?srt:NULL, ug);
purge_phase_0(idx, link, bub, s, hits, k_trans, del, hh);
}
for (k = len = 0; k < ug->g->n_seq; k++) {
if(del[k] != 2) {
if((hh->h.a[k].n == hh->h.a[k].m) && (hh->h.a[k].m == ((uint64_t)asm_opt.polyploidy))) {
len += ug->g->seq[k].len;
}
continue;
}
assert(hh->h.a[k].n < hh->h.a[k].m);
hh->a.a[hh->h.a[k].a+hh->h.a[k].n] = hapid;
hh->h.a[k].n++; len += ug->g->seq[k].len;
}
fprintf(stderr, "[M::%s::stat] # hap%lu bases: %lu\n", __func__, hapid+1, len);
bm = bub->index; bub->index = bidx; bidx = bm;
exchange_kv_u_trans_t(ref, &(idx->t_ch->k_trans));
}
int hic_short_align_mmhap(const enzyme *fn1, const enzyme *fn2, ha_ug_index* idx, ug_opt_t *opt, kvec_pe_hit **rhits, mmhap_t **rh)
{
if(rh) (*rh) = NULL;
sldat_t sl;
sl.idx = idx;
sl.t_ch = idx->t_ch;
sl.chunk_size = 20000000;
sl.n_thread = asm_opt.thread_num;
sl.total_base = sl.total_pair = 0;
idx->hap_cnt = asm_opt.hap_occ;
kv_init(sl.hits.a); kv_init(sl.hits.idx); kv_init(sl.hits.occ);
if(!load_hc_hits(&sl.hits, idx->ug, asm_opt.output_file_name)) {
alignment_worker_pipeline(&sl, fn1, fn2);
write_hc_hits(&sl.hits, idx->ug, asm_opt.output_file_name);
}
sl.hits.uID_bits = idx->uID_bits; sl.hits.pos_mode = idx->pos_mode;
if(sl.hits.idx.n == 0) idx_hc_links(&(sl.hits), idx, NULL);
mmhap_t *hh = gen_mmhap_t(idx->ug, idx->read_g, opt->sources);
bubble_type *bub = gen_mmhap_bub(idx->ug, idx->t_ch->ir_het, &(idx->t_ch->k_trans), hh);
hc_links link; init_hc_links(&link, idx->ug->g->n_seq, idx->t_ch);
measure_distance(idx, idx->ug, &sl.hits, &link, bub, &(idx->t_ch->k_trans));
kv_u_trans_t k_trans; kv_init(k_trans); kv_init(k_trans.idx);
ps_t *s = init_ps_t(11, idx->ug->g->n_seq); ///H_partition hap;
uint64_t k, n_hap = asm_opt.polyploidy; asg64_v srt; kv_init(srt);
uint8_t *ff; CALLOC(ff, idx->ug->g->n_seq);
uint32_t *bidx; MALLOC(bidx, idx->ug->g->n_seq);
kv_u_trans_t r_trans_buf; kv_init(r_trans_buf); kv_init(r_trans_buf.idx);
for (k = 0; k < n_hap; k++) {
purge_phase(idx, hh, k, ((n_hap>k)?(n_hap-k-1):(0)), &link, &sl.hits, bub, s, &k_trans, ff, bidx, &r_trans_buf, &srt);
}
if(rh) (*rh) = hh;
// if(rhits) (*rhits) = get_r_hits_order(&sl.hits, idx->uID_bits, idx->pos_mode, idx->read_g, idx->ug, &bub);
kv_destroy(sl.hits.a); kv_destroy(sl.hits.idx); kv_destroy(sl.hits.occ);
destory_hc_links(&link);
kv_destroy(k_trans); kv_destroy(k_trans.idx);
destory_ps_t(&s); destory_bubbles(bub); free(bub);
kv_destroy(srt); free(ff); free(bidx);
kv_destroy(r_trans_buf); kv_destroy(r_trans_buf.idx);
return 1;
}
void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, mmhap_t **rh, kvec_pe_hit **rhits)
{
ug_index = NULL;
int exist = (asm_opt.load_index_from_disk?
@@ -17442,8 +17714,9 @@ void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt,
ug_index->read_g = read_g;
ug_index->t_ch = t_ch;
///test_unitig_index(ug_index, ug);
if(!is_poy) hic_short_align(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt, rhits);
else hic_short_align_poy(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt);
if(!rh) hic_short_align(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt, rhits);
else hic_short_align_mmhap(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt, rhits, rh);
// else hic_short_align_poy(asm_opt.hic_reads[0], asm_opt.hic_reads[1], ug_index, opt);
destory_hc_pt_index(ug_index);free(ug_index);
hic.h
+1 -1
View File
@@ -113,7 +113,7 @@ pdq* pqw, uint32_t* path_w, buf_t *resw, asg_t *sg, uint8_t *dest, uint8_t df, u
long long *dis);
void set_utg_by_dis(uint32_t v, pdq* pq, asg_t *g, kvec_t_u32_warp *res, uint32_t dis);
void dedup_hits(kvec_pe_hit* hits, uint64_t is_dup);
void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, uint32_t is_poy, kvec_pe_hit **rhits);
void hic_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, mmhap_t **rh, kvec_pe_hit **rhits);
spg_t *hic_pre_analysis(ma_ug_t *ug, asg_t* read_g, trans_chain* t_ch, ug_opt_t *opt, kvec_pe_hit **rhits);
void prt_bubble_gfa_adv(FILE *fp, bubble_type *bub, const char* utg_pre, const char* bub_pre, const char* chain_pre);
void bp_solve(ug_opt_t *opt, kv_u_trans_t *ref, ma_ug_t *ug, asg_t *sg, bubble_type *bub, double cis_rate);
inter.cpp
+518 -36
View File
@@ -64,6 +64,8 @@ int64_t ug_map_lchain_simple(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl
// #define GBIN_L 15000
#define GBIN_L 256
#define FREE_BATCH 16
#define generic_key(x) (x)
KRADIX_SORT_INIT(gfa64, uint64_t, generic_key, 8)
@@ -417,6 +419,8 @@ typedef struct { // global data structure for kt_pipeline()
int32_t is_HPC, bw, max_gap, chn_pen_gap, n_thread, is_cnt, is_ovlp, mini_cut, chain_cut, keep_unsymm_arc;
ul_idx_t udb;
kv_u_trans_t *filter;
uint32_t *free_cnt;
ug_rid_cov_t *ccov;
} ug_trans_t;
typedef struct { // global data structure for kt_pipeline()
@@ -441,6 +445,7 @@ typedef struct { // global data structure for kt_pipeline()
// bit_mask_t *bm;
} ug_bin_t;
void hc_glchain_destroy(glchain_t *b)
{
if (!b) return;
@@ -9573,7 +9578,8 @@ static void worker_for_trans_ovlp(void *data, long i, int tid) // callback for k
}
void filter_by_reliable_ovlp_adv(uint32_t id, kv_u_trans_t *idx, st_mt_t *sp, overlap_region_alloc* ol, const ul_idx_t *udb, double sec_rate, uint64_t avoid_dup_aln, uint64_t *occ1)
void filter_by_reliable_ovlp_adv(uint32_t id, kv_u_trans_t *idx, st_mt_t *sp, overlap_region_alloc* ol, const ul_idx_t *udb, double sec_rate, uint64_t avoid_dup_aln,
uint64_t dedup_by_reliable_ovlp, uint64_t *occ1)
{
(*occ1) = 0;
u_trans_t *a; uint64_t n, k, l, s, e, s0, e0, z, ov, rr, r1, spn; overlap_region *m, t;
@@ -9629,25 +9635,28 @@ void filter_by_reliable_ovlp_adv(uint32_t id, kv_u_trans_t *idx, st_mt_t *sp, ov
}
}
for (k = sp->n = 0; k < n; k++) {
if(a[k].del) continue;
if(a[k].f == RC_0 || a[k].f == RC_1) {
kv_push(uint64_t, *sp, (((uint64_t)a[k].qs)<<32)|((uint64_t)a[k].qe));
}
}
if(sp->n > 1) {
radix_sort_gfa64(sp->a, sp->a + sp->n);
for (k = z = 0; k < sp->n; k++) {
s = sp->a[k]>>32; e = (uint32_t)sp->a[k];
if(z > 0 && s <= ((uint32_t)sp->a[z-1])) {
if(e > ((uint32_t)sp->a[z-1])) {
sp->a[z-1] >>= 32; sp->a[z-1] <<= 32; sp->a[z-1] |= e;
}
} else {
sp->a[z++] = sp->a[k];
sp->n = 0;
if(dedup_by_reliable_ovlp) {
for (k = sp->n = 0; k < n; k++) {
if(a[k].del) continue;
if(a[k].f == RC_0 || a[k].f == RC_1) {
kv_push(uint64_t, *sp, (((uint64_t)a[k].qs)<<32)|((uint64_t)a[k].qe));
}
}
sp->n = z;
if(sp->n > 1) {
radix_sort_gfa64(sp->a, sp->a + sp->n);
for (k = z = 0; k < sp->n; k++) {
s = sp->a[k]>>32; e = (uint32_t)sp->a[k];
if(z > 0 && s <= ((uint32_t)sp->a[z-1])) {
if(e > ((uint32_t)sp->a[z-1])) {
sp->a[z-1] >>= 32; sp->a[z-1] <<= 32; sp->a[z-1] |= e;
}
} else {
sp->a[z++] = sp->a[k];
}
}
sp->n = z;
}
}
for (k = rr = r1 = 0, spn = sp->n; k < ol->length; k++) {
@@ -9659,7 +9668,7 @@ void filter_by_reliable_ovlp_adv(uint32_t id, kv_u_trans_t *idx, st_mt_t *sp, ov
// ol->list[k].y_id+1, "lc"[udb->ug->u.a[ol->list[k].y_id].circ], udb->ug->u.a[ol->list[k].y_id].len,
// ol->list[k].y_pos_s, ol->list[k].y_pos_e+1, ol->list[k].x_pos_strand);
// }
if(m->x_pos_strand == 0) {
if((dedup_by_reliable_ovlp) && (m->x_pos_strand == 0)) {
s = m->x_pos_s; e = m->x_pos_e + 1; l = 0;
for (z = 0; z < spn; z++) {
s0 = sp->a[z]>>32; e0 = (uint32_t)sp->a[z];
@@ -9743,8 +9752,8 @@ uint64_t split_ug_lalign(uint64_t ol_h, overlap_region_alloc* ol, double errh, d
overlap_region *aux_o, int64_t sid, uint64_t khit, uint64_t chain_cut, void *km)
{
uint64_t ol_l = ol->length - ol_h, on0, k, m; int64_t wl; double erate; overlap_region t;
// if(sid == 160) fprintf(stderr, "[M::%s] errh::%f, errl::%f\n", __func__, errh, errl);
// if(sid == 160) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "st");
// if(sid == 57) fprintf(stderr, "[M::%s] errh::%f, errl::%f\n", __func__, errh, errl);
// if(sid == 57) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "st");
if(ol_h) {
erate = errh; on0 = ol->length;
wl = MIN((((double)THRESHOLD_MAX_SIZE)/erate), WINDOW);
@@ -9760,7 +9769,7 @@ uint64_t split_ug_lalign(uint64_t ol_h, overlap_region_alloc* ol, double errh, d
}
ol_h = ol->length; ol->length = m; ol_l = ol->length - ol_h;
}
// if(sid == 160) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "mi");
// if(sid == 57) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "mi");
if(ol_l) {
erate = errl; on0 = ol->length;
wl = MIN((((double)THRESHOLD_MAX_SIZE)/erate), WINDOW);
@@ -9774,10 +9783,10 @@ uint64_t split_ug_lalign(uint64_t ol_h, overlap_region_alloc* ol, double errh, d
m++;
}
}
// if(sid == 7) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "sw");
// if(sid == 57) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "sw");
ol->length = ol_l;
ug_lalign(ol, cl, uref, uopt, qstr, ql, qu, tu, dumy, exz, aux_o, erate, wl, sid, khit, chain_cut, km);
// if(sid == 7) {
// if(sid == 57) {
// fprintf(stderr, "[M::%s::]\ton0::%lu\tol_l::%lu\tol_h::%lu\tol->length::%lu\n", __func__,
// on0, ol_l, ol_h, ol->length);
// prt_split_ovs(ol, uref->ug, ol_h, ol_l, "u0");
@@ -9791,7 +9800,7 @@ uint64_t split_ug_lalign(uint64_t ol_h, overlap_region_alloc* ol, double errh, d
m++;
}
ol_l = ol->length; ol->length = m; ol_h = ol->length - ol_l;
// if(sid == 7) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "u1");
// if(sid == 57) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "u1");
if(ol_l) {///swap ol_h and ol_l
for (k = ol_l, m = 0; k < ol->length; k++) {
if(k != m) {
@@ -9803,7 +9812,7 @@ uint64_t split_ug_lalign(uint64_t ol_h, overlap_region_alloc* ol, double errh, d
}
}
}
// if(sid == 160) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "ed");
// if(sid == 57) prt_split_ovs(ol, uref->ug, ol_h, ol_l, "ed");
return ol_h;
}
@@ -9870,6 +9879,50 @@ uint32_t test_het_aln(ma_ug_t *ug, uint64_t rid, u_trans_t *a, uint64_t a_n, st_
return 0;
}
uint32_t is_mmhom_node(uint64_t *ca, ma_utg_t *u, asg_t *sg, uint64_t cov_bd, double cut_rate)
{
if(cut_rate < 0) cut_rate = 0; if(cut_rate > 1.0) cut_rate = 1.0;
uint64_t k, a, na, a_cut = u->n*cut_rate, na_cut = u->n*(1.0-cut_rate);
for (k = a = na = 0; k < u->n; k++) {
if(ca[k] > (cov_bd*((uint64_t)sg->seq[u->a[k]>>33].len))) {
a++; if((a) && (a>=a_cut)) return 1;
} else {
na++; if((na) && (na>=na_cut)) return 0;
}
}
if((a) && (a>=a_cut)) return 1;
return 0;
}
uint32_t test_het_aln_mmhap(uint64_t uid, ug_rid_cov_t *ccov, u_trans_t *a, uint64_t a_n, overlap_region_alloc* ol, st_mt_t *sp)
{
uint64_t k; u_trans_t p;
if(a_n == 0 && ol->length == 0) return 0;
kv_resize(uint64_t, *sp, ccov->ug->u.a[uid].n);
memcpy(sp->a, ccov->cov.a+ccov->idx[uid], sizeof((*(sp->a)))*ccov->ug->u.a[uid].n);
for (k = 0; k < a_n; k++) {
if(a[k].del) continue;
append_cov_line_ug_rid_cov_t(uid, sp->a, &(a[k]), ccov, ((uint64_t)-1), -1);
}
for (k = 0; k < ol->length; k++) {
p.qn = uid; p.tn = ol->list[k].y_id;
p.rev = ol->list[k].y_pos_strand; p.f = RC_3; p.nw = 0;
p.qs = ol->list[k].x_pos_s; p.qe = ol->list[k].x_pos_e+1;
if(p.rev) {
p.ts = ccov->ug->u.a[p.tn].len - (ol->list[k].y_pos_e+1);
p.te = ccov->ug->u.a[p.tn].len - ol->list[k].y_pos_s;
} else {
p.ts = ol->list[k].y_pos_s;
p.te = ol->list[k].y_pos_e+1;
}
append_cov_line_ug_rid_cov_t(uid, sp->a, &p, ccov, ((uint64_t)-1), -1);
}
return is_mmhom_node(sp->a, &(ccov->ug->u.a[uid]), ccov->rg, ccov->hom_min, 0.8);
}
void push_ul_ov_t(ul_idx_t *udb, u_trans_t *a, uint64_t a_n, uint64_t rid, st_mt_t *sp, overlap_region_alloc* ol, uint64_t len, uint64_t is_arc_filter, double max_err, kv_ul_ov_t *res)
{
uint64_t cnt, z, k, l, m, spn; ul_ov_t *p;
@@ -9921,11 +9974,10 @@ void push_ul_ov_t(ul_idx_t *udb, u_trans_t *a, uint64_t a_n, uint64_t rid, st_mt
// fprintf(stderr, ">0<[M::%s] utg%.6u%c -> utg%.6u%c\n", __func__,
// p->qn+1, "lc"[s->udb.ug->u.a[p->qn].circ],
// p->tn+1, "lc"[s->udb.ug->u.a[p->tn].circ]);
// if(i == 5)
// {
// fprintf(stderr, "***utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\ti::%ld\n",
// p->qn+1, "lc"[s->ug->u.a[p->qn].circ], s->ug->u.a[p->qn].len, p->qs, p->qe, "+-"[p->rev],
// p->tn+1, "lc"[s->ug->u.a[p->tn].circ], s->ug->u.a[p->tn].len, p->ts, p->te, i);
// if(p->ts >= p->te || p->qs >= p->qe) {
// fprintf(stderr, "+[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\n", __func__,
// p->qn+1, "lc"[udb->ug->u.a[p->qn].circ], udb->ug->u.a[p->qn].len, p->qs, p->qe, "+-"[p->rev],
// p->tn+1, "lc"[udb->ug->u.a[p->tn].circ], udb->ug->u.a[p->tn].len, p->ts, p->te);
// }
if((is_arc_filter) && (!trans_ovlp_connect(p, udb->ug))) res->n--;
// fprintf(stderr, ">1<[M::%s] utg%.6u%c -> utg%.6u%c\n", __func__,
@@ -9948,6 +10000,11 @@ void push_ul_ov_t(ul_idx_t *udb, u_trans_t *a, uint64_t a_n, uint64_t rid, st_mt
p->qs = a[sp->a[k]].qs; p->qe = a[sp->a[k]].qe;
p->ts = a[sp->a[k]].ts; p->te = a[sp->a[k]].te;
p->sec = (p->qe-p->qs)*max_err;
// if(p->ts >= p->te || p->qs >= p->qe) {
// fprintf(stderr, "-[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\n", __func__,
// p->qn+1, "lc"[udb->ug->u.a[p->qn].circ], udb->ug->u.a[p->qn].len, p->qs, p->qe, "+-"[p->rev],
// p->tn+1, "lc"[udb->ug->u.a[p->tn].circ], udb->ug->u.a[p->tn].len, p->ts, p->te);
// }
}
// if(is_sec_filter) {
@@ -10038,7 +10095,7 @@ uint64_t gen_trans_adaptive_aln(ug_trans_t *s, uint64_t rid, ha_ovec_buf_t *b, k
// fprintf(stderr, "-2-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length);
// }
filter_by_reliable_ovlp_adv(rid, s->filter, &(b->sp), &b->olist, &(s->udb), s->sec_cutoff, 1, &ol_h);
filter_by_reliable_ovlp_adv(rid, s->filter, &(b->sp), &b->olist, &(s->udb), s->sec_cutoff, 1, 1, &ol_h);
clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2);
if(!fi) ol_h = 0;
@@ -10064,6 +10121,47 @@ uint64_t gen_trans_adaptive_aln(ug_trans_t *s, uint64_t rid, ha_ovec_buf_t *b, k
return pass_aln;
}
void clear_count_buf(ug_trans_t *s, uint32_t tid, uint32_t free_count)
{
// fprintf(stderr, "[M::%s]\n", __func__);
ha_ovec_buf_t *b = s->hab[tid];
destory_fake_cigar(&(b->tmp_region.f_cigar));
free(b->tmp_region.w_list.a); free(b->tmp_region.w_list.c.a);
memset(&(b->tmp_region), 0, sizeof(b->tmp_region));
init_fake_cigar(&(b->tmp_region.f_cigar));
memset(&(b->tmp_region.w_list), 0, sizeof(b->tmp_region.w_list));
CALLOC(b->tmp_region.w_list.a, 1); b->tmp_region.w_list.n = b->tmp_region.w_list.m = 1;
ha_abufl_destroy(b->abl); b->abl = ha_abufl_init();
kv_destroy(b->sp); memset(&(b->sp), 0, sizeof((b->sp)));
if(free_count) return;
destory_Candidates_list(&b->clist);
memset((&(b->clist)), 0, sizeof(b->clist));
init_Candidates_list(&b->clist);
destory_overlap_region_alloc(&b->olist);
memset((&(b->olist)), 0, sizeof(b->olist));
init_overlap_region_alloc(&b->olist);
destory_UC_Read(&b->self_read);
memset((&(b->self_read)), 0, sizeof(b->self_read));
init_UC_Read(&b->self_read);
destory_UC_Read(&b->ovlp_read);
memset((&(b->ovlp_read)), 0, sizeof(b->ovlp_read));
init_UC_Read(&b->ovlp_read);
destory_Correct_dumy(&b->correct);
memset((&(b->correct)), 0, sizeof(b->correct));
init_Correct_dumy(&b->correct);
destroy_bit_extz_t(&(b->exz));
memset((&(b->exz)), 0, sizeof(b->exz));
init_bit_extz_t(&(b->exz), 31);
}
static void worker_for_trans_ovlp_adv(void *data, long i, int tid) // callback for kt_for()
{
ug_trans_t *s = (ug_trans_t*)data;
@@ -10085,6 +10183,10 @@ static void worker_for_trans_ovlp_adv(void *data, long i, int tid) // callback f
s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3, s->is_HPC, s->idx_a.a + s->idx_n.a[i], 0, NULL, 0, s->mini_cut, s->chain_cut, NULL);
assert(cnt == ((s->idx_n.a[i+1]-s->idx_n.a[i])));
}
if(s->free_cnt[tid] >= FREE_BATCH) {
clear_count_buf(s, tid, 1); s->free_cnt[tid] = 0;
}
s->free_cnt[tid]++;
return;
}
@@ -10096,6 +10198,240 @@ static void worker_for_trans_ovlp_adv(void *data, long i, int tid) // callback f
if(!gen_trans_adaptive_aln(s, i, b, bl, seq, len, s->filter, s->diff_ec_ul, s->diff_ec_ul_double, s->bw_thres, s->bw_thres_double)) {
gen_trans_adaptive_aln(s, i, b, bl, seq, len, NULL, s->diff_ec_ul_double, s->diff_ec_ul_double, s->bw_thres_double, s->bw_thres_double);
}
if(s->free_cnt[tid] >= FREE_BATCH) {
clear_count_buf(s, tid, 0); s->free_cnt[tid] = 0;
}
s->free_cnt[tid]++;
}
uint64_t *gen_reliable_cov_arr(uint32_t id, kv_u_trans_t *idx, ug_rid_cov_t *ccov, st_mt_t *sp)
{
u_trans_t *a; uint64_t n, k;
a = u_trans_a(*idx, id); n = u_trans_n(*idx, id);
kv_resize(uint64_t, *sp, ccov->ug->u.a[id].n);
memcpy(sp->a, ccov->cov.a+ccov->idx[id], sizeof((*(sp->a)))*ccov->ug->u.a[id].n);
for (k = 0; k < n; k++) {
if(a[k].del) continue;
if(a[k].f == RC_0 || a[k].f == RC_1) {
append_cov_line_ug_rid_cov_t(id, sp->a, &(a[k]), ccov, ((uint64_t)-1), -1);
}
}
return sp->a;
}
uint64_t is_above_cov(uint64_t uid, overlap_region *o, uint64_t *fc, ug_rid_cov_t *ccov, double sec_rate)
{
u_trans_t p;
p.qn = uid; p.tn = o->y_id;
p.rev = o->y_pos_strand; p.f = RC_3; p.nw = 0;
p.qs = o->x_pos_s; p.qe = o->x_pos_e+1;
if(p.rev) {
p.ts = ccov->ug->u.a[p.tn].len - (o->y_pos_e+1);
p.te = ccov->ug->u.a[p.tn].len - o->y_pos_s;
} else {
p.ts = o->y_pos_s;
p.te = o->y_pos_e+1;
}
if(append_cov_line_ug_rid_cov_t(uid, fc, &p, ccov, ccov->hom_max, sec_rate)) return 0;
return 1;
}
void filter_by_reliable_ovlp_mmhap_adv(uint32_t id, kv_u_trans_t *idx, st_mt_t *sp, overlap_region_alloc* ol, const ul_idx_t *udb, double sec_rate, uint64_t avoid_dup_aln,
uint64_t dedup_by_reliable_ovlp, ug_rid_cov_t *ccov, uint64_t *occ1)
{
(*occ1) = 0;
u_trans_t *a; uint64_t n, k, l, z, rr, r1, *fc; overlap_region *m, t;
a = u_trans_a(*idx, id); n = u_trans_n(*idx, id);
if(avoid_dup_aln) {
kv_resize(uint64_t, *sp, (ol->length)+n);
for (k = sp->n = 0; k < n; k++) {
if(a[k].del) continue;
if(a[k].f == RC_0 || a[k].f == RC_1) {
z = a[k].tn; z <<= 1; z |= a[k].rev; z <<= 32;
kv_push(uint64_t, *sp, z);
}
}
if(sp->n > 0) {
for (k = 0; k < ol->length; k++) {
z = ol->list[k].y_id; z <<= 1; z |= ol->list[k].y_pos_strand;
z <<= 32; z |= k; z |= ((uint64_t)0x80000000);
kv_push(uint64_t, *sp, z);
}
radix_sort_gfa64(sp->a, sp->a + sp->n);
for (k = 1, l = 0, rr = 0; k <= sp->n; k++) {
if(k == sp->n || (sp->a[l]>>32)!=(sp->a[k]>>32)) {
if((k - l > 1) && (!(sp->a[l]&((uint64_t)0x80000000)))) {///overlap within bck
for (z = l; z < k; z++) {
if(sp->a[z]&((uint64_t)0x80000000)) {
ol->list[(uint32_t)(sp->a[z]-((uint64_t)0x80000000))].y_id = ((uint32_t)-1);
rr++;
}
}
}
l = k;
}
}
if(rr > 0) {
for (k = rr = 0; k < ol->length; k++) {
if(ol->list[k].y_id == ((uint32_t)-1)) continue;
if(rr != k) {
t = ol->list[rr];
ol->list[rr] = ol->list[k];
ol->list[k] = t;
}
rr++;
}
ol->length = rr;
}
}
}
sp->n = 0; fc = NULL;
if(dedup_by_reliable_ovlp) {
for (k = 0; k < n; k++) {
if(a[k].del) continue;
if(a[k].f == RC_0 || a[k].f == RC_1) break;
}
if(k < n) fc = gen_reliable_cov_arr(id, idx, ccov, sp);
}
for (k = rr = r1 = 0; k < ol->length; k++) {
m = &(ol->list[k]);
if((dedup_by_reliable_ovlp) && (m->x_pos_strand == 0) && (fc)) {
if(is_above_cov(id, m, fc, ccov, sec_rate)) continue;
}
if(rr != k) {
t = ol->list[k];
ol->list[k] = ol->list[rr];
ol->list[rr] = t;
}
if(ol->list[rr].x_pos_strand) {
ol->list[rr].x_pos_strand = 0;
if(r1 != rr) {
t = ol->list[r1];
ol->list[r1] = ol->list[rr];
ol->list[rr] = t;
}
r1++;
}
rr++;
}
// if(id == 1576) {
// fprintf(stderr, "[M::%s] utg%.6ul, ol->length0::%lu, ol->length::%lu\n", __func__, id+1, ol->length, rr);
// }
ol->length = rr; (*occ1) = r1;
}
uint64_t gen_trans_adaptive_mmhap_aln(ug_trans_t *s, uint64_t rid, ha_ovec_buf_t *b, kv_ul_ov_t *bl, char *seq, uint64_t len, kv_u_trans_t *fi, double err_low, double err_high, double bw_low, double bw_high)
{
uint64_t cnt = ((s->idx_n.a[rid+1]-s->idx_n.a[rid])), ol_h = 0, pass_aln = 0;
uint32_t high_occ = asm_opt.polyploidy + 1; overlap_region *aux_o = NULL;
///note: high_occ is different
ug_map_lchain(b->abl, rid, seq, len, s->w, s->k, &(s->udb), &b->olist, &b->clist, bw_low, bw_high,
s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3,
s->is_HPC, s->idx_a.a + s->idx_n.a[rid], cnt, s->srt_a.a, s->srt_a.n, s->mini_cut, s->chain_cut, fi);
// if(rid == 57) {
// fprintf(stderr, "-1-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length);
// }
///remove candidate chains that have been calculated
if(!fi) backward_dedup_ol(rid, bl, &(b->sp), &b->olist);///it is ok
// if(rid == 57) {
// fprintf(stderr, "-2-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length);
// }
filter_by_reliable_ovlp_mmhap_adv(rid, s->filter, &(b->sp), &b->olist, &(s->udb), s->sec_cutoff, 1, 1, s->ccov, &ol_h);
clear_Cigar_record(&b->cigar1); clear_Round2_alignment(&b->round2);
if(!fi) ol_h = 0;
// if(rid == 57) {
// fprintf(stderr, "-3-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length);
// }
ol_h = split_ug_lalign(ol_h, &b->olist, err_high, err_low,
&b->clist, &(s->udb), s->uopt, seq, len, &b->self_read, &b->ovlp_read,
&b->correct, &b->exz, aux_o, rid, s->k, s->chain_cut, NULL);
// if(rid == 57) {
// fprintf(stderr, "-4-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length);
// }
aux_o = gen_aux_ovlp(&b->olist);///must be here
// if(rid == 57) {
// fprintf(stderr, "-5-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length);
// }
ol_h = split_ug_lalign(ol_h, &b->olist, err_high, err_low,
&b->clist, &(s->udb), s->uopt, seq, len, &b->self_read, &b->ovlp_read,
&b->correct, &b->exz, aux_o, rid, s->k, s->chain_cut, NULL);
// if(rid == 57) {
// fprintf(stderr, "-6-[M::%s] utg%.6lu%c, rid::%ld, b->olist->length::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, b->olist.length);
// }
if(fi) {///first round
pass_aln = test_het_aln_mmhap(rid, s->ccov, u_trans_a((*(s->filter)), rid), u_trans_n((*(s->filter)), rid), &b->olist, &(b->sp));
push_ul_ov_t(&(s->udb), u_trans_a((*(s->filter)), rid), u_trans_n((*(s->filter)), rid), rid, &(b->sp), &b->olist, len, pass_aln, err_high, bl);
// fprintf(stderr, "-1-[M::%s] utg%.6lu%c, rid::%lu, pass_aln::%lu\n",
// __func__, rid+1, "lc"[s->ug->u.a[rid].circ], rid, pass_aln);
} else {///second round
push_ul_ov_t(&(s->udb), NULL, 0, rid, &(b->sp), &b->olist, len, 0, err_high, bl);
remove_trans_ovlp_connect(s->udb.ug, rid, bl);
}
return pass_aln;
}
static void worker_for_trans_ovlp_mmhap_adv(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]; kv_ul_ov_t *bl = &(s->ll[tid].tk);
uint32_t high_occ = asm_opt.polyploidy + 1; uint64_t cnt;
char *seq = s->ug->u.a[i].s; int64_t len = s->ug->u.a[i].len;
if((!s->is_ovlp) && (s->is_cnt)) s->idx_n.a[i] = 0;
if(s->ug->g->seq[i].del) return;
if(is_mmhom_node(s->ccov->cov.a+s->ccov->idx[i], &(s->ug->u.a[i]), s->ccov->rg, s->ccov->hom_min, 0.9)) return;
// 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;
if(!s->is_ovlp) {
if(s->is_cnt) {
s->idx_n.a[i] = ug_map_lchain(b->abl, i, seq, len, s->w, s->k, &(s->udb), NULL, NULL, s->bw_thres, s->bw_thres_double,
s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3, s->is_HPC, NULL, 0, NULL, 0, s->mini_cut, s->chain_cut, NULL);
} else {
cnt = ug_map_lchain(b->abl, i, seq, len, s->w, s->k, &(s->udb), NULL, NULL, s->bw_thres, s->bw_thres_double,
s->max_n_chain, 1, NULL, &(b->tmp_region), NULL, &(b->sp), &high_occ, NULL, 0, 1, 0.2, 3, s->is_HPC, s->idx_a.a + s->idx_n.a[i], 0, NULL, 0, s->mini_cut, s->chain_cut, NULL);
assert(cnt == ((s->idx_n.a[i+1]-s->idx_n.a[i])));
}
if(s->free_cnt[tid] >= FREE_BATCH) {
clear_count_buf(s, tid, 1); s->free_cnt[tid] = 0;
}
s->free_cnt[tid]++;
return;
}
// if(i == 58) {
// fprintf(stderr, "\n-1-[M::%s] utg%.6u%c, rid::%ld, is_ovlp::%d, is_cnt::%d, len::%ld, str::%u\n",
// __func__, (uint32_t)i+1, "lc"[s->ug->u.a[i].circ], i, s->is_ovlp, s->is_cnt, len, (uint32_t)(!!seq));
// }
if(!gen_trans_adaptive_mmhap_aln(s, i, b, bl, seq, len, s->filter, s->diff_ec_ul, s->diff_ec_ul_double, s->bw_thres, s->bw_thres_double)) {
gen_trans_adaptive_mmhap_aln(s, i, b, bl, seq, len, NULL, s->diff_ec_ul_double, s->diff_ec_ul_double, s->bw_thres_double, s->bw_thres_double);
}
if(s->free_cnt[tid] >= FREE_BATCH) {
clear_count_buf(s, tid, 0); s->free_cnt[tid] = 0;
}
s->free_cnt[tid]++;
}
@@ -19562,7 +19898,7 @@ void clear_all_ul_t(all_ul_t *x)
void init_ug_trans_t(ug_trans_t *opt, ug_opt_t *uopt, int32_t is_HPC, int32_t k, int32_t w, int32_t max_n_chain,
double bw_thres, double diff_ec_ul, double bw_thres_double, double diff_ec_ul_double, double sec_cutoff, int32_t n_thread,
int32_t mini_cut, int32_t chain_cut, int32_t keep_unsymm_arc, ma_ug_t *ug, asg_t *sg, bubble_type *bub)
int32_t mini_cut, int32_t chain_cut, int32_t keep_unsymm_arc, ma_ug_t *ug, asg_t *sg, bubble_type *bub, uint8_t gen_bub)
{
int64_t i; uint8_t *bf = NULL;
memset(opt, 0, sizeof((*opt)));
@@ -19586,6 +19922,7 @@ int32_t mini_cut, int32_t chain_cut, int32_t keep_unsymm_arc, ma_ug_t *ug, asg_t
opt->rg = sg;
opt->n_thread = ((n_thread>=1)?n_thread:1);
CALLOC(opt->free_cnt, opt->n_thread);
CALLOC(opt->hab, opt->n_thread);
CALLOC(opt->ll, opt->n_thread);
for (i = 0; i < opt->n_thread; ++i) {
@@ -19598,7 +19935,7 @@ int32_t mini_cut, int32_t chain_cut, int32_t keep_unsymm_arc, ma_ug_t *ug, asg_t
opt->udb.ug = ug;
if(bub) {
opt->bub = bub;
} else {
} else if(gen_bub) {
opt->bub = gen_bubble_chain(sg, ug, uopt, &bf, 0); free(bf);
}
}
@@ -19727,6 +20064,7 @@ void gen_trans_base_count_comp(ug_trans_t *p, kv_u_trans_t *res)
clean_u_trans_t_idx_adv(res, p->ug, p->rg); p->filter = res;
// fprintf(stderr, "[M::%s::] ==> 0\n", __func__);
p->is_cnt = 1; p->is_ovlp = 0;
memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread);
kt_for(p->n_thread, worker_for_trans_ovlp_adv, p, p->ug->u.n);
for (i = l = 0; i < p->ug->u.n; i++) {
occ = p->idx_n.a[i]; p->idx_n.a[i] = l; l += occ;
@@ -19736,6 +20074,7 @@ void gen_trans_base_count_comp(ug_trans_t *p, kv_u_trans_t *res)
p->idx_a.n = p->idx_a.m = l; MALLOC(p->idx_a.a, p->idx_a.n);
// fprintf(stderr, "[M::%s::] ==> 1\n", __func__);
p->is_cnt = 0; p->is_ovlp = 0;
memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread);
kt_for(p->n_thread, worker_for_trans_ovlp_adv, p, p->ug->u.n);
p->srt_a.n = p->srt_a.m = p->idx_a.n; MALLOC(p->srt_a.a, p->srt_a.n);
// fprintf(stderr, "[M::%s::] p->idx_a.n::%lu \n", __func__, (uint64_t)p->idx_a.n);
@@ -19770,13 +20109,14 @@ void gen_trans_base_count_comp(ug_trans_t *p, kv_u_trans_t *res)
// fprintf(stderr, "[M::%s::] ==> 2\n", __func__);
p->is_cnt = 0; p->is_ovlp = 1;
memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread);
kt_for(p->n_thread, worker_for_trans_ovlp_adv, p, p->ug->u.n);
// fprintf(stderr, "[M::%s::] ==> 3\n", __func__);
for (i = 0; (int64_t)i < p->n_thread; i++) {
ha_ovec_destroy(p->hab[i]);
free(p->ll[i].lo.a); free(p->ll[i].srt.a.a); free(p->ll[i].tc.a);
}
free(p->idx_a.a); free(p->idx_n.a); free(p->hab);
free(p->idx_a.a); free(p->idx_n.a); free(p->hab); free(p->free_cnt);
// destory_bubbles(p->bub); free(p->bub);
// fprintf(stderr, "[M::%s::] ==> 4\n", __func__);
@@ -19838,13 +20178,144 @@ void gen_trans_base_count_comp(ug_trans_t *p, kv_u_trans_t *res)
fprintf(stderr, "[M::%s::%.3f] ==> Qualification\n", __func__, yak_realtime()-index_time);
}
void gen_trans_base_count_mmhap_comp(ug_trans_t *p, kv_u_trans_t *res)
{
double index_time = yak_realtime();
// ha_flt_tab = NULL;
uint64_t i, k, l, occ, m, cc; kv_ul_ov_t *bl = NULL;
u_trans_t *z; ha_mzl_t *tz; double ww;
p->ccov = gen_ug_rid_cov_t(p->ug, p->rg, p->uopt->sources);
clean_u_trans_t_idx_adv(res, p->ug, p->rg); p->filter = res;
// fprintf(stderr, "[M::%s::] ==> 0\n", __func__);
p->is_cnt = 1; p->is_ovlp = 0;
memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread);
kt_for(p->n_thread, worker_for_trans_ovlp_mmhap_adv, p, p->ug->u.n);
for (i = l = 0; i < p->ug->u.n; i++) {
occ = p->idx_n.a[i]; p->idx_n.a[i] = l; l += occ;
}
// fprintf(stderr, "[M::%s::] i::%lu, l::%lu\n", __func__, i, l);
p->idx_n.a[i] = l;
p->idx_a.n = p->idx_a.m = l; MALLOC(p->idx_a.a, p->idx_a.n);
// fprintf(stderr, "[M::%s::] ==> 1\n", __func__);
p->is_cnt = 0; p->is_ovlp = 0;
memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread);
kt_for(p->n_thread, worker_for_trans_ovlp_mmhap_adv, p, p->ug->u.n);
p->srt_a.n = p->srt_a.m = p->idx_a.n; MALLOC(p->srt_a.a, p->srt_a.n);
// fprintf(stderr, "[M::%s::] p->idx_a.n::%lu \n", __func__, (uint64_t)p->idx_a.n);
// memcpy(p->srt_a.a, p->idx_a.a, p->srt_a.n*sizeof((*(p->srt_a.a))));
for (i = 0; i < p->srt_a.n; i++) {
p->srt_a.a[i] = p->idx_a.a[i];
p->srt_a.a[i].pos = (uint32_t)i;
p->srt_a.a[i].rid = i>>32;
}
radix_sort_ha_mzl_t_srt(p->srt_a.a, p->srt_a.a + p->srt_a.n);
kvec_t(uint64_t) cut; kv_init(cut);
for (k = 1, l = 0; k <= p->srt_a.n; k++) {
if(k == p->srt_a.n || p->srt_a.a[l].x != p->srt_a.a[k].x) {
for (i = l; i < k; i++) {
m = p->srt_a.a[i].rid; m <<= 32; m |= p->srt_a.a[i].pos;
assert(p->srt_a.a[i].x == p->idx_a.a[m].x);
p->srt_a.a[i] = p->idx_a.a[m]; p->idx_a.a[m].x = i;
}
kv_push(uint64_t, cut, (k - l));
l = k;
}
}
if(cut.n > 0) {
radix_sort_gfa64(cut.a, cut.a + cut.n);
m = cut.n * 0.0002; cc = cut.a[cut.n-1] + 1;
if(m > 0 && m <= cut.n) cc = cut.a[cut.n-m] + 1;
if(cc < (uint64_t)p->mini_cut) p->mini_cut = cc;
}
kv_destroy(cut);
// fprintf(stderr, "[M::%s::] p->mini_cut::%d \n", __func__, p->mini_cut);
// fprintf(stderr, "[M::%s::] ==> 2\n", __func__);
p->is_cnt = 0; p->is_ovlp = 1;
memset(p->free_cnt, 0, sizeof((*(p->free_cnt)))*p->n_thread);
kt_for(p->n_thread, worker_for_trans_ovlp_mmhap_adv, p, p->ug->u.n);
// fprintf(stderr, "[M::%s::] ==> 3\n", __func__);
for (i = 0; (int64_t)i < p->n_thread; i++) {
ha_ovec_destroy(p->hab[i]);
free(p->ll[i].lo.a); free(p->ll[i].srt.a.a); free(p->ll[i].tc.a);
}
free(p->idx_a.a); free(p->idx_n.a); free(p->hab); free(p->free_cnt);
destory_ug_rid_cov_t(p->ccov); free(p->ccov);
// destory_bubbles(p->bub); free(p->bub);
// fprintf(stderr, "[M::%s::] ==> 4\n", __func__);
///make results consistent
kv_resize(ha_mzl_t, p->srt_a, p->ug->u.n); p->srt_a.n = p->ug->u.n;
for (i = 0; i < p->srt_a.n; i++) {
tz = &(p->srt_a.a[i]);
tz->x = (uint64_t)-1; tz->rev = 0;
tz->pos = tz->rid = tz->span = 0;
}
// memset(p->srt_a.a, 0, sizeof((*(p->srt_a.a)))*p->srt_a.n);
for (i = 0, occ = res->n; (int64_t)i < p->n_thread; i++) {
bl = &(p->ll[i].tk);
if(!(bl->n)) continue;
for (k = 1, l = 0; k <= bl->n; k++) {
if(k == bl->n || bl->a[k].qn != bl->a[l].qn) {
if(k > l) {
tz = &(p->srt_a.a[bl->a[l].qn]);
tz->x = bl->a[l].qn; tz->x <<= 32; tz->x |= i;
tz->rid = l>>32; tz->pos = (uint32_t)l; tz->rev = 1;
occ += (k - l);
}
l = k;
}
}
}
// fprintf(stderr, "[M::%s::] ==> 5\n", __func__);
// kt_for(p->n_thread, worker_for_sysm_trans_ovlp, p, p->ug->u.n);///not correct
// assert(p->srt_a.n <= p->ug->u.n);
// radix_sort_ha_mzl_t_srt(p->srt_a.a, p->srt_a.a + p->srt_a.n);
kv_resize(u_trans_t, *res, occ);
for (i = 0; i < p->srt_a.n; i++) {
tz = &(p->srt_a.a[i]);
if(!(tz->rev)) continue;
bl = &(p->ll[(uint32_t)(tz->x)].tk);
k = tz->rid; k <<= 32; k += tz->pos;
assert(bl->a[k].qn == (tz->x>>32));
for (; (k < bl->n) && (bl->a[k].qn == (tz->x>>32)); k++) {
if(bl->a[k].qn == bl->a[k].tn) continue;
ww = cal_trans_ov_w(&(bl->a[k]));
if(ww <= 0) continue;
kv_pushp(u_trans_t, *res, &z);
z->f = RC_3; z->rev = bl->a[k].rev; z->del = 0;
z->qn = bl->a[k].qn; z->qs = bl->a[k].qs; z->qe = bl->a[k].qe;
z->tn = bl->a[k].tn; z->ts = bl->a[k].ts; z->te = bl->a[k].te;
z->nw = ww;
// if(z->ts >= z->te || z->qs >= z->qe) {
// fprintf(stderr, "[M::%s]\tutg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\n", __func__,
// z->qn+1, "lc"[p->ug->u.a[z->qn].circ], p->ug->u.a[z->qn].len, z->qs, z->qe, "+-"[z->rev],
// z->tn+1, "lc"[p->ug->u.a[z->tn].circ], p->ug->u.a[z->tn].len, z->ts, z->te);
// }
// if(z->qn == 56 || z->qn == 160 || z->tn == 56 || z->tn == 160) {
// fprintf(stderr, ">>>utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\tnw::%f\n",
// z->qn+1, "lc"[p->ug->u.a[z->qn].circ], p->ug->u.a[z->qn].len, z->qs, z->qe, "+-"[z->rev],
// z->tn+1, "lc"[p->ug->u.a[z->tn].circ], p->ug->u.a[z->tn].len, z->ts, z->te, z->nw);
// }
// if(z->nw <= 0) res->n--;
}
}
// fprintf(stderr, "[M::%s::] ==> 6\n", __func__);
for (i = 0; (int64_t)i < p->n_thread; i++) free(p->ll[i].tk.a);
free(p->srt_a.a); free(p->ll);
fprintf(stderr, "[M::%s::%.3f] ==> Qualification\n", __func__, yak_realtime()-index_time);
}
void trans_base_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res, bubble_type *bub)
{
ug_trans_t sl;
init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov);
init_ug_trans_t(&sl, uopt, 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, 1.0-asm_opt.trans_base_rate_sec, 1.0-asm_opt.trans_base_rate_sec,
0.85, asm_opt.thread_num, 512, 3, 1, ug, sg, bub);
0.85, asm_opt.thread_num, 512, 3, 1, ug, sg, bub, 1);
// gen_trans_base_count(&sl, res);
gen_trans_base_count_comp(&sl, res);
if(!bub) {
@@ -19852,6 +20323,17 @@ void trans_base_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res,
}
}
void trans_base_mmhap_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res)
{
ug_trans_t sl;
init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov);
init_ug_trans_t(&sl, uopt, 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, 1.0-asm_opt.trans_base_rate_sec, 1.0-asm_opt.trans_base_rate_sec,
0.85, asm_opt.thread_num, 512, 3, 1, ug, sg, NULL, 0);
// gen_trans_base_count(&sl, res);
gen_trans_base_count_mmhap_comp(&sl, res);
}
void init_ug_bin_t(ug_bin_t *sl, const ug_opt_t *uopt, int32_t is_HPC, int32_t k, int32_t w, int32_t max_n_chain,
double bw_thres, double diff_ov, double diff_bin, uint64_t max_diff, uint64_t min_bin_len, int32_t n_thread,
int32_t mini_cut, int32_t chain_cut, int32_t keep_unsymm_arc, ma_ug_t *ug, asg_t *sg)
inter.h
+1
View File
@@ -126,5 +126,6 @@ uint32_t infer_se(uint32_t qs, uint32_t qe, uint32_t ts, uint32_t te, uint32_t r
uint32_t rqs, uint32_t rqe, uint32_t *rts, uint32_t *rte);
uint32_t clean_contain_g(const ug_opt_t *uopt, asg_t *sg, uint32_t push_trans);
void dedup_contain_g(const ug_opt_t *uopt, asg_t *sg);
void trans_base_mmhap_infer(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, kv_u_trans_t *res);
#endif