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

r721

Browse Source
This commit is contained in:
chhylp123
2025-03-14 01:40:21 -04:00
parent a96191560e
commit b3b18ab1d0
4 changed files with 246 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CommandLines.h
+1 -1
View File
@@ -5,7 +5,7 @@
#include <pthread.h>
#include <stdint.h>
#define HA_VERSION "0.25.0-r710"
#define HA_VERSION "0.25.0-r721"
#define VERBOSE 0
Correct.cpp
+103 -17
View File
@@ -17423,11 +17423,11 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u
void hc_ovlp_base_direct(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, int64_t wl, All_reads *rref, char* qstr, UC_Read *tu,
bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, uint64_t rid)
bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, uint64_t rid, int64_t pre_mode)
{
int64_t i, l, mode, q[2], t[2], qr, tr, is_done, zn;
int64_t i, l, mode, q[2], t[2], qr, tr, is_done, zn, si, ei;
if(z->non_homopolymer_errors == 0 && z->w_list.n) {
if((pre_mode < 0) && (z->non_homopolymer_errors == 0) && (z->w_list.n)) {
zn = z->w_list.n;
for (i = 1; i < zn; i++) {
if((z->w_list.a[i].error == 0 && z->w_list.a[i-1].error == 0) && (z->w_list.a[i].x_start == z->w_list.a[i-1].x_end + 1) &&
@@ -17461,7 +17461,16 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u
}
}
for (l = -1, i = 0; i <= ch_n; i++) {
si = 0; ei = ch_n;
if(pre_mode == 0) {
si = 1; ei = ch_n - 1;
} else if(pre_mode == 1) {
si = 1;
} else if(pre_mode == 2) {
ei = ch_n - 1;
}
for (l = si - 1, i = si; i <= ei; i++) {
q[0] = q[1] = t[0] = t[1] = mode = -1; is_done = 0;
if(l >= 0) {
q[0] = ch_a[l].self_offset; t[0] = ch_a[l].offset;
@@ -17702,7 +17711,7 @@ All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate, int64_
// idx.qs, idx.qe, idx.ts, idx.te, idx.qn, idx.tn);
// }
// ovlp_base_aln(z, ch_a, ch_n, &idx, wl, uref, hpc_g, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, (uint64_t)-1);
hc_ovlp_base_direct(z, ch_a, ch_n, wl, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, (uint64_t)-1);
hc_ovlp_base_direct(z, ch_a, ch_n, wl, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, (uint64_t)-1, mode);
an = aux_o->w_list.n; q[0] = q[1] = t[0] = t[1] = 0; todo = 0;
// if(z->x_id == 29033 && z->y_id == 21307) {
// fprintf(stderr, "[M::%s]\tan::%ld\n", __func__, an);
@@ -17820,11 +17829,17 @@ uint64_t gen_hc_fast_cigar0(overlap_region *z, Candidates_list *cl, uint64_t wl,
aux_o->x_pos_s = z->x_pos_s; aux_o->x_pos_e = z->x_pos_e;
aux_o->y_pos_s = z->y_pos_s; aux_o->y_pos_e = z->y_pos_e;
hc_ovlp_base_direct(z, ch_a, ch_n, wl, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, rid);
hc_ovlp_base_direct(z, ch_a, ch_n, wl, rref, qstr, tu, exz, aux_o, e_rate, ql, tl, rid, -1);
int64_t aux_n = aux_o->w_list.n;
for (i = 0; i < aux_n; i++) {
// if(z->y_id == 30129) {
// fprintf(stderr, "[aln::-i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, clen::%u, mode::%d\n", i,
// aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start,
// aux_o->w_list.a[i].x_start, aux_o->w_list.a[i].x_end+1,
// aux_o->w_list.a[i].y_start, aux_o->w_list.a[i].y_end+1,
// aux_o->w_list.a[i].error, aux_o->w_list.a[i].clen, aux_o->w_list.a[i].error_threshold);
// }
if(!(is_ualn_win(aux_o->w_list.a[i]))) continue;
// if((aux_o->w_list.a[i].x_end+1-aux_o->w_list.a[i].x_start) <= FORCE_CNS_L) {
// fprintf(stderr, "[aln::-i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, clen::%u, mode::%d\n", i,
@@ -19761,7 +19776,7 @@ inline uint64_t set_cgid(ul_ov_t *z, uint64_t *ia, uint64_t *iak, uint64_t ian,
}
///need to print some examples for double check
uint64_t is_get_group(ul_ov_t *a, uint64_t *ga, uint64_t *ia, uint64_t gi, uint64_t tn)
uint64_t is_get_group(ul_ov_t *a/**, uint64_t an, uint64_t rid**/, uint64_t *ga, uint64_t *ia, uint64_t gi, uint64_t tn)
{
// fprintf(stderr, "+[M::%s] gi::%lu\n", __func__, gi);
uint64_t k = ga[gi]>>32;
@@ -19771,7 +19786,10 @@ uint64_t is_get_group(ul_ov_t *a, uint64_t *ga, uint64_t *ia, uint64_t gi, uint6
// fprintf(stderr, "-[M::%s] k::%lu\n", __func__, k);
// }
for (k = ga[gi]>>32; (k != ((uint64_t)-1)) && (a[k].tn != tn); k = ia[k]);
if((k != ((uint32_t)-1)) && (a[k].tn == tn)) return 0;
// if((k != ((uint32_t)-1)) && (k >= an)) {
// fprintf(stderr, "-[M::%s] rid::%lu, an::%lu, k::%lu\n", __func__, rid, an, k);
// }
if((k != ((uint64_t)-1)) && (a[k].tn == tn)) return 0;
return 1;
}
@@ -19867,12 +19885,13 @@ int64_t rphase_lidel_cc(overlap_region_alloc* oa, ul_ov_t *a, int64_t an, double
ol = (a[k].qe - a[k].qs) * len_st; if(ol < len_w) ol = len_w;
if((a[k].qe - a[k].qs) < ol) continue;
mk = -1; msw = -1;
for (z = k + 1; (z < an) && (a[z].qs < a[k].qe); z++) {
for (z = 0/**k + 1**/; (z < an) && (a[z].qs < a[k].qe); z++) {
if(a[z].sec < c_sz) continue;
if(a[z].ts == ((uint32_t)-1)) continue;
if(z == k) continue;
sw = cal_lindel_dd(&(a[k]), &(a[z]), len_st, len_w, err_dif, c_sz);
if(sw < 0) continue;
if((sw > msw) && (is_get_group(a, ca, ia, a[z].ts, a[k].tn))) {
if((sw > msw) && (is_get_group(a, /**an, rid,**/ ca, ia, a[z].ts, a[k].tn))) {
mk = a[z].ts; msw = sw;
}
}
@@ -22979,7 +22998,7 @@ 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->y_id == 4378833 || z->y_id == 4378837) {
// if(z->y_id == 30129) {
// 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);
@@ -22998,14 +23017,14 @@ int64_t return_t_chain(overlap_region *z, Candidates_list *cl)
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->y_id == 30129) fprintf(stderr, "\n-a-[M::%s]\tcn::%ld\n", __func__, cn);
// if(z->y_id == 4378833 || z->y_id == 4378837) {
// if(z->y_id == 30129) {
// 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;
// i = z->shared_seed; pid = cl->list[i].readID; cn = cl->length;
// 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,
@@ -25532,7 +25551,70 @@ uint32_t inline ff_tend(overlap_region *z, int64_t wn, int64_t dn, double dr, do
return 0;
}
void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop)
uint32_t inline ff_lunalign(overlap_region *z, double erate, double gap_rate, int64_t max_gap)
{
if((z->w_list.n == 1) && (!(is_ualn_win(z->w_list.a[0])))
&& (z->w_list.a[0].x_start == ((int64_t)z->x_pos_s)) && (z->w_list.a[0].x_end == ((int64_t)z->x_pos_e))
&& (z->w_list.a[0].y_start == ((int64_t)z->y_pos_s)) && (z->w_list.a[0].y_end == ((int64_t)z->y_pos_e))) {
return 1;
}
int64_t k, zwn = z->w_list.n, zq, zt, wq, wt, tot_e, tot_g, ql, tl;
// fprintf(stderr, "[M::%s]\n", __func__);
zq = z->x_pos_s; zt = z->y_pos_s; tot_e = tot_g = 0;
for (k = 0; k < zwn; k++) {
// fprintf(stderr, "[M::%s]\twk::%ld\tq::[%d, %d)\tt::[%d, %d)\terr::%d\n", __func__,
// k, z->w_list.a[k].x_start, z->w_list.a[k].x_end + 1, z->w_list.a[k].y_start, z->w_list.a[k].y_end + 1, z->w_list.a[k].error);
if(is_ualn_win(z->w_list.a[k])) continue;
wq = z->w_list.a[k].x_start;
wt = z->w_list.a[k].y_start;
if(wq != zq) tot_g += ((wq>=zq)?(wq-zq):(zq-wq));
if(wt != zt) tot_g += ((wt>=zt)?(wt-zt):(zt-wt));
zq = z->w_list.a[k].x_end + 1;
zt = z->w_list.a[k].y_end + 1;
tot_e += z->w_list.a[k].error;
// fprintf(stderr, "[M::%s]\twk::%ld\tq::[%d, %d)\tt::[%d, %d)\terr::%d\n", __func__,
// k, z->w_list.a[k].x_start, z->w_list.a[k].x_end + 1, z->w_list.a[k].y_start, z->w_list.a[k].y_end + 1, z->w_list.a[k].error);
}
wq = z->x_pos_e + 1;
wt = z->y_pos_e + 1;
if(wq != zq) tot_g += ((wq>=zq)?(wq-zq):(zq-wq));
if(wt != zt) tot_g += ((wt>=zt)?(wt-zt):(zt-wt));
// fprintf(stderr, "[M::%s]\t%.*s(id::%u)\tq::[%u, %u)\t%.*s(id::%u)\tt::[%u, %u)\tre::%ld\trg::%ld\n", __func__, (int)Get_NAME_LENGTH(R_INF, z->x_id), Get_NAME(R_INF, z->x_id), z->x_id, z->x_pos_s, z->x_pos_e + 1,
// (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), z->y_id, z->y_pos_s, z->y_pos_e + 1, tot_e, tot_g);
if(!tot_g) return 1;
// fprintf(stderr, "-0-[M::%s]\n", __func__);
if(tot_g > max_gap) return 0;
// fprintf(stderr, "-1-[M::%s]\n", __func__);
ql = z->x_pos_e + 1 - z->x_pos_s;
tl = z->y_pos_e + 1 - z->y_pos_s;
if((tot_g > (ql*gap_rate)) || (tot_g > (tl*gap_rate))) return 0;
// fprintf(stderr, "-2-[M::%s]\n", __func__);
tot_e += tot_g;
if((tot_e > (ql*erate)) || (tot_e > (tl*erate))) return 0;
// fprintf(stderr, "-3-[M::%s]\n", __func__);
return 1;
}
void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max)
{
uint64_t i, bs, k, ql = qu->length; Window_Pool w; double err, e_max, rr; int64_t re;
overlap_region t; overlap_region *z; //asg64_v iidx, buf, buf1;
@@ -25566,6 +25648,8 @@ void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rre
if(!gen_hc_fast_cigar(z, cl, rref, w.window_length, qu->seq, tu, exz, aux_o, e_rate, ql, rid, khit, &re)) continue;
if((align_gap_max >= 0) && (!ff_lunalign(z, err, align_gap_rate, align_gap_max))) continue;
if(chem_drop && ff_tend(z, 384, 2000, 0.1, (((e_rate*10)<0.36)?(e_rate*10):(0.36)), 128)) continue;
@@ -25590,7 +25674,7 @@ void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rre
}
void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop)
void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max)
{
uint64_t i, bs, k, ql = qu->length; Window_Pool w; double err, e_max, rr; int64_t re;
overlap_region t; overlap_region *z; //asg64_v iidx, buf, buf1;
@@ -25632,6 +25716,8 @@ void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads
if(!gen_hc_fast_cigar(z, cl, rref, w.window_length, qu->seq, tu, exz, aux_o, e_rate, ql, rid, khit, &re)) continue;
if((align_gap_max >= 0) && (!ff_lunalign(z, err, align_gap_rate, align_gap_max))) continue;
if(chem_drop && ff_tend(z, 384, 2000, 0.1, (((e_rate*10)<0.36)?(e_rate*10):(0.36)), 128)) continue;
// if(z->x_id == 3196 && z->y_id == 3199) fprintf(stderr, "-1-[M::%s] tid::%u\t%.*s\trr::%f\tre::%ld\n", __func__, z->y_id, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), rr, re);
Correct.h
+2 -2
View File
@@ -1391,8 +1391,8 @@ int64_t get_rid_backward_cigar_err(rtrace_iter *it, ul_ov_t *aln, kv_rtrace_t *t
const ul_idx_t *uref, char* qstr, UC_Read *tu, overlap_region_alloc *ol, overlap_region *o,
bit_extz_t *exz, double e_rate, int64_t qs);
void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop);
void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop);
void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max);
void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max);
uint64_t gen_hc_r_alin_re(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf);
void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, int64_t wl, int64_t ql, uint8_t occ_thres/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs, Chain_Data *dp, asg8_v *q8, asg8_v *t8, uint8_t lindel);
void set_exact_exz(bit_extz_t *exz, int64_t qs, int64_t qe, int64_t ts, int64_t te);
ecovlp.cpp
+140 -9
View File
@@ -65,6 +65,12 @@ typedef struct {
uint8_t *cr;
} ec_ovec_buf_t;
typedef struct {
ec_ovec_buf_t *p;
asg64_v idx;
ma_ug_t *ug;
} ec_polish_buf_t;
typedef struct {
uint32_t n_thread, n_a, chunk_size, cn;
FILE *fp;
@@ -2801,7 +2807,7 @@ inline uint64_t exact_ec_check(char *qstr, uint64_t ql, char *tstr, uint64_t tl,
return 0;
}
void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in, uint8_t chem_drop)
void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v *buf, asg64_v *srt, ma_hit_t_alloc *in, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max)
{
if(ol->length <= 0) return;
@@ -2817,7 +2823,7 @@ void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads *
}
if(!(srt->n)) {
gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop);
gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max);
} else {
///debug for memory
// snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n);
@@ -2852,9 +2858,10 @@ void gen_hc_r_alin_ea(overlap_region_alloc* ol, Candidates_list *cl, All_reads *
// snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n);
if(on > nec) {
gen_hc_r_alin_nec(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop);
gen_hc_r_alin_nec(ol, cl, rref, qu, tu, exz, aux_o, e_rate, wl, rid, khit, move_gap, buf, chem_drop, align_gap_rate, align_gap_max);
}
// fprintf(stderr, "[M::%s] srt->n::%u, nec::%lu, on::%lu\n", __func__, (uint32_t)srt->n, nec, on);
///debug for memory
// snprintf(NULL, 0, "dwn::%u\tdcn::%u", (uint32_t)aux_o->w_list.n, (uint32_t)aux_o->w_list.c.n);
}
@@ -3257,6 +3264,8 @@ static void worker_hap_ec(void *data, long i, int tid)
///for debug indel
// if(i != 1238) return;
// if(i != 2410) return;
// if(i != 4198005) return;
// if(i != 682) return;
// debug_retrive_bqual(D, &b->v8t, i, 256); return;
@@ -3276,7 +3285,7 @@ static void worker_hap_ec(void *data, long i, int tid)
///mz1_ha_sketch(rs, rl, mz_w, mz_k, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0, NULL);
// if((asm_opt.is_ont) && (b->olist.length)) get_mz1(qu->seq, qu->length, RES_W, RES_K, 0, !(asm_opt.flag & HA_F_NO_HPC), b->ab, NULL, NULL, asm_opt.mz_sample_dist, NULL, NULL, NULL, -1, asm_opt.dp_min_len, -1, &(b->sp), asm_opt.mz_rewin, 0, NULL, 0);
gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT/**asm_opt.k_mer_length**/, 1, &b->v16, &b->v64, &(R_INF.paf[i]), asm_opt.is_ont);
gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT/**asm_opt.k_mer_length**/, 1, &b->v16, &b->v64, &(R_INF.paf[i]), asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1));
///for debug indel
// prt_ovlp_sam(&b->olist, &b->ovlp_read, b->self_read.seq, b->self_read.length);
@@ -3404,7 +3413,7 @@ static void worker_hap_ec_dbg_paf(void *data, long i, int tid)
aux_o = fetch_aux_ovlp(&b->olist);///must be here
// stderr_phase_ovlp(&b->olist);
gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, 1, &b->v16, &b->v64, &(R_INF.paf[i]), 0);
gen_hc_r_alin_ea(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, 1, &b->v16, &b->v64, &(R_INF.paf[i]), 0, -1, -1);
uint32_t k, m, tl; overlap_region *z; bit_extz_t ez; ma_hit_t *t;
for (k = 0; k < b->olist.length; k++) {
@@ -3796,6 +3805,42 @@ static void worker_hap_dc_ec(void *data, long i, int tid)
refresh_ec_ovec_buf_t0(b, REFRESH_N);
}
static void worker_update_dc_ec(void *data, long i, int tid)
{
ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
uint64_t k; ma_hit_t *z;
// fprintf(stderr, "-0-[M::%s-beg] rid->%ld\n", __func__, i);
// if (memcmp("m64012_190921_234837/139067658/ccs", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// fprintf(stderr, "-0-[M::%s-beg] rid->%ld\n", __func__, i);
// } else if (memcmp("m64012_190921_234837/28968323/ccs", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
// fprintf(stderr, "-1-[M::%s-beg] rid->%ld\n", __func__, i);
// } else {
// return;
// }
// if(i != 2851) return;
// if(scb.a[i].m < scc.a[i].n) {
// scb.a[i].m = scc.a[i].n;
// REALLOC(scb.a[i].a, scb.a[i].m);
// }
// scb.a[i].n = scc.a[i].n;
// memcpy(scb.a[i].a, scc.a[i].a, scc.a[i].n*sizeof((*(scb.a[i].a))));
if(!(R_INF.paf[i].length)) return;
recover_UC_Read(&b->self_read, &R_INF, i);
for (k = 0; k < R_INF.paf[i].length; k++) {
z = &(R_INF.paf[i].buffer[k]);
if((z->el) && (quick_exact_match(z, &R_INF, &b->self_read, &b->ovlp_read, &scc))) {
z->el = 1; b->cnt[0]++;
} else {
z->el = 0; b->cnt[1]++;
}
}
refresh_ec_ovec_buf_t0(b, REFRESH_N);
}
void flip_paf_rc(uint64_t rid, ma_hit_t_alloc *paf, All_reads *rref)
{
@@ -4280,6 +4325,52 @@ static void worker_hap_dc_ec_chemical_arc_mark(void *data, long i, int tid)
refresh_ec_ovec_buf_t0(b, REFRESH_N);
}
uint64_t get_candidate_rrs(ma_utg_t *u, uint64_t rz)
{
// uint64_t rid, rs, re, rev, k, l[2], lr, ts, te;
// ma_hit_t_alloc *z = NULL;
// ma_hit_t *h = NULL;
// rid = u->a[rz]>>33; rev = (u->a[rz]>>32)&1;
// rs = 0; re = (uint32_t)u->a[rz];
// if(!rev) {
// rs = Get_READ_LENGTH(R_INF, rid) - ((uint32_t)u->a[rz]);
// re = Get_READ_LENGTH(R_INF, rid);
// }
// for (k = lr = 0; k < rz; k++) lr += (uint32_t)u->a[k];
// ts = lr; te = lr + (uint32_t)u->a[k];
// z = &(sources[rid]);
// for (k = 0; k < z->length; k++) {
// h = &(z->buffer[k]);
// if ((Get_qs(*h) <= rs) && (Get_qe(*h) >= re)) {
// l[0] = l[1] = 0;
// if(!(h->rev&rev)) {
// l[0] = Get_ts(*h); l[1] = Get_READ_LENGTH(R_INF, Get_tn(*h)) - Get_te(*h);
// } else {
// l[1] = Get_ts(*h); l[0] = Get_READ_LENGTH(R_INF, Get_tn(*h)) - Get_te(*h);
// }
// }
// }
return 1;
}
static void worker_ec_polish(void *data, long i, int tid)
{
ec_ovec_buf_t0 *b = &(((ec_polish_buf_t*)data)->p->a[tid]);
uint64_t *idx = &(((ec_polish_buf_t*)data)->idx.a[i]);
if(get_candidate_rrs(&(((ec_polish_buf_t*)data)->ug->u.a[(*idx)>>32]), (uint32_t)(*idx))) {
*idx = (uint64_t)-1;
}
refresh_ec_ovec_buf_t0(b, REFRESH_N);
}
void gen_ovlst_paf(ma_hit_t_alloc *in_e, ma_hit_t_alloc *in_r, asg64_v *ou)
{
uint32_t n = 0, k;
@@ -4482,7 +4573,7 @@ overlap_region* h_ec_lchain_re1(ha_abuf_t *ab, uint32_t rid, UC_Read *qu, UC_Rea
// fprintf(stderr, "-0-[M::%s]\n", __func__);
gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf, 0); rs = qu->seq;
gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf, 0, -1, -1); rs = qu->seq;
// fprintf(stderr, "-1-[M::%s]\n", __func__);
@@ -5151,7 +5242,7 @@ overlap_region* h_ec_lchain_re3(ha_abuf_t *ab, uint32_t rid, UC_Read *qu, UC_Rea
// fprintf(stderr, "-0-[M::%s]\n", __func__);
gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf, 0); rs = qu->seq;
gen_hc_r_alin(ol, cl, rref, qu, tu, exz, aux_o, asm_opt.max_ov_diff_ec, w.window_length, rid, E_KHIT, 1, buf, 0, -1, -1); rs = qu->seq;
// fprintf(stderr, "-1-[M::%s]\n", __func__);
@@ -6001,6 +6092,22 @@ uint64_t cal_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, uint64
return num_correct;
}
void cal_update_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a)
{
double tt0 = yak_realtime_0();
uint64_t k, num_ec_o = 0, num_nec_o = 0;
for (k = 0; k < n_thre; ++k) b->a[k].cnt[0] = b->a[k].cnt[1] = 0;
kt_for(n_thre, worker_update_dc_ec, b, n_a);///debug_for_fix
for (k = 0; k < n_thre; ++k) {
num_ec_o += b->a[k].cnt[0]; num_nec_o += b->a[k].cnt[1];
}
fprintf(stderr, "[M::pec::%.3f] # exact o: %lu; # non-exact o: %lu\n", yak_realtime_0()-tt0, num_ec_o, num_nec_o);
}
void ha_print_ovlp_stat_1(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a)
{
@@ -6162,7 +6269,7 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u
{
// write_ec_reads("ec0.fa");
fprintf(stderr, "[M::%s]\tn_thre::%lu, round::%lu, n_round::%lu, n_a::%lu, is_sv::%lu\n", __func__, n_thre, round, n_round, n_a, is_sv);
// fprintf(stderr, "[M::%s]\tn_thre::%lu, round::%lu, n_round::%lu, n_a::%lu, is_sv::%lu\n", __func__, n_thre, round, n_round, n_a, is_sv);
ec_ovec_buf_t *b = NULL; uint64_t k, is_cr = (round&1);
(*tot_b) = (*tot_e) = 0;
@@ -6177,7 +6284,10 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u
sl_ec_r(n_thre, n_a);
}
if(is_sv) kt_for(n_thre, worker_hap_dc_ec, b, n_a);///update overlaps
cal_update_ec_multiple(b, n_thre, n_a);///update overlaps
// if(is_sv) kt_for(n_thre, worker_hap_dc_ec, b, n_a);///update overlaps
if((!is_sv) || (is_sv && is_cr)) {
kt_for(n_thre, worker_hap_post_rev, b, n_a);
@@ -6379,4 +6489,25 @@ uint8_t* gen_chemical_arc_rf(uint64_t n_thre, uint64_t n_a)
b->cr = NULL; destroy_ec_ovec_buf_t(b);
return ra;
}
void gen_hc_polish(uint64_t n_thre, ma_ug_t *ug)
{
ec_polish_buf_t b; memset(&b, 0, sizeof(b)); b.ug = ug;
uint64_t k, z, zn;
for (k = b.idx.n = 0; k < ug->u.n; k++) {
b.idx.n += ug->u.a[k].n;
}
MALLOC(b.idx.a, b.idx.n);
for (k = zn = 0; k < ug->u.n; k++) {
for (z = 0; z < ug->u.a[k].n; z++, zn++) {
b.idx.a[zn] = (k<<32)|z;
}
}
b.p = gen_ec_ovec_buf_t(n_thre);
kt_for(n_thre, worker_ec_polish, &b, b.idx.n);
destroy_ec_ovec_buf_t(b.p); free(b.idx.a);
}