diff --git a/CommandLines.cpp b/CommandLines.cpp index c50258d..40cf1a9 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -93,6 +93,7 @@ static ko_longopt_t long_options[] = { { "ec-only", ko_no_argument, 375}, { "hyb-syn", ko_required_argument, 376}, { "simd-m", ko_required_argument, 377}, + { "del-hf", ko_no_argument, 378}, // { "path-round", ko_required_argument, 348}, { 0, 0, 0 } }; @@ -438,9 +439,11 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->hyb_syn = 1; asm_opt->step_rd = -1/**128**/; - asm_opt->dbg_bam = 0/**1**/; + asm_opt->dbg_bam = 0; asm_opt->simd_mm = -1; + + asm_opt->del_hf = 0; } void destory_enzyme(enzyme* f) @@ -1116,6 +1119,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) asm_opt->hyb_syn = atoi(opt.arg); } else if (c == 377) { asm_opt->simd_mm = atoi(opt.arg); + } else if (c == 378) { + asm_opt->del_hf = 1; } else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg); } diff --git a/CommandLines.h b/CommandLines.h index 7d2e16f..e71ba18 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.25.1-r891" +#define HA_VERSION "0.25.1-r897" #define VERBOSE 0 @@ -203,6 +203,7 @@ typedef struct { uint8_t dbg_bam; int8_t simd_mm; + int8_t del_hf; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Correct.cpp b/Correct.cpp index cde2e0d..aa772ac 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -2350,6 +2350,7 @@ inline char *return_str_seq(char *buf, int64_t s, int64_t pri_l, uint8_t rev, hp if(!hpc_g) { memset(buf, 'N', aux_beg); retrieve_u_seq(NULL, buf+aux_beg, &(uref->ug->u.a[id]), rev, s, pri_l, NULL); + // retrieve_u_seq_fast(NULL, buf, &(uref->ug->u.a[id]), rev, s, pri_l, NULL); memset(buf+aux_beg+pri_l, 'N', aux_end); return buf; } else { @@ -2369,6 +2370,7 @@ inline char *return_str_seq_exz(char *buf, int64_t s, int64_t pri_l, uint8_t rev { if(!hpc_g) { retrieve_u_seq(NULL, buf, &(uref->ug->u.a[id]), rev, s, pri_l, NULL); + // retrieve_u_seq_fast(NULL, buf, &(uref->ug->u.a[id]), rev, s, pri_l, NULL); return buf; } else { return hpc_str(*hpc_g, id, rev) + s; @@ -10732,8 +10734,9 @@ void generate_haplotypes_sv(haplotype_evdience_alloc* hap, overlap_region_alloc* -inline int64_t comput_sc_rphase(SnpStats *ai, uint64_t id, SnpStats *aj, uint64_t jd, haplotype_evdience *za, uint64_t occ0_cut) +inline int64_t comput_sc_rphase(SnpStats *ai, uint64_t id, SnpStats *aj, uint64_t jd, haplotype_evdience *za, uint64_t occ0_cut, uint32_t *oid, uint8_t *oph) { + (*oid) = ((uint32_t)-1); (*oph) = ((uint8_t)-1); if(ai->site == aj->site) return INT64_MIN; // if(ai->occ_0 < occ0_cut || aj->occ_0 < occ0_cut) return INT64_MIN; haplotype_evdience *iz = NULL, *jz = NULL; int64_t in, jn, ik, jk, nn[2]; uint8_t fi, fj; @@ -10750,6 +10753,7 @@ inline int64_t comput_sc_rphase(SnpStats *ai, uint64_t id, SnpStats *aj, uint64_ } else if(iz[ik].overlapSite == id){ fi = 1; } + (*oph) = fi; fj = 2; if(hh_tp(jz[jk]) == 0) { @@ -10763,7 +10767,12 @@ inline int64_t comput_sc_rphase(SnpStats *ai, uint64_t id, SnpStats *aj, uint64_ } if(fi == 2 || fj == 2) return INT64_MIN; - if(fi != fj) return INT64_MIN; + if(fi != fj) { + if(((*oph) == 0) || ((*oph) == 1)) { + (*oid) = iz[ik].overlapID; + } + return INT64_MIN; + } nn[fi]++; } } @@ -10873,7 +10882,7 @@ void gen_rphase_dp0_multiple_path(SnpStats *a, int64_t an, haplotype_evdience *z for (i = 0; i < an; ++i) { max_f = 1; ch_n = idx->n; for (j = i - 1; j >= 0; --j) { - sc = comput_sc_rphase(&a[i], i, &a[j], j, za, 0); + sc = comput_sc_rphase(&a[i], i, &a[j], j, za, 0, NULL, NULL); if (sc == INT64_MIN) continue; sc += f[j]; if (sc > max_f) { @@ -11000,7 +11009,7 @@ void gen_rphase_dp0_single_path(SnpStats *a, int64_t an, haplotype_evdience *za, max_f = 1; max_j = -1; st = 0; ///if(a[i].occ_0 < cc) st = i; for (j = i - 1; j >= st; --j) { - sc = comput_sc_rphase(&a[i], i, &a[j], j, za, 0/**cc**/); + sc = comput_sc_rphase(&a[i], i, &a[j], j, za, 0/**cc**/, NULL, NULL); if (sc == INT64_MIN) continue; sc += f[j]; if (sc > max_f) { @@ -11392,11 +11401,35 @@ void gen_rphase_dp0_single_path_hybrid_0_multi(SnpStats *a, int64_t an, haplotyp } } +void label_skip_mm(SnpStats *a, uint32_t oid, uint8_t hf, uint32_t **oidx, uint32_t **o_a, int32_t *ma, int32_t mm, uint32_t mek, int32_t max_n, int32_t rid) +{ + uint32_t sn = oidx[hf][oid+1] - oidx[hf][oid], *sa = o_a[hf] + oidx[hf][oid], k; int32_t i; + + if(hf) { + for (k = 0; k < sn; k++) { + if(sa[k] == ((uint32_t)-1)) continue; + if((sa[k] >= mek) && (a[sa[k]].site > a[mek].site)) break; + ma[sa[k]] = mm; + } + } else { + for (k = 0; k < sn; k++) { + if(sa[k] == ((uint32_t)-1)) continue; + if((sa[k] >= mek) && (a[sa[k]].site > a[mek].site)) break; + for (i = sa[k]; (i >= 0) && (a[sa[k]].site == a[i].site); i--) { + // if(i < 0 || i >= max_n) { + // /**if(rid == 307)**/ fprintf(stderr, "[M::%s::]\ti::%d\tmax_n::%d\trid::%d\toid::%u\n", __func__, i, max_n, rid, oid); + // } + + ma[i] = mm; + } + } + } +} -void gen_rphase_dp0_single_path_multi(SnpStats *a, int64_t an, haplotype_evdience *za, Chain_Data *dp, asg64_v *idx, int64_t het_cov, int64_t hom_cov, int64_t n_hap, double cut_rate, uint64_t cut_bd, asg64_v *res, uint8_t *qual_a, uint8_t site_sc, int64_t hap_cov_match, int64_t hap_cov_unmatch) +void gen_rphase_dp0_single_path_multi(overlap_region_alloc *ol, SnpStats *a, int64_t an, haplotype_evdience *za, Chain_Data *dp, asg64_v *idx, int64_t het_cov, int64_t hom_cov, int64_t n_hap, double cut_rate, uint64_t cut_bd, asg64_v *res, uint8_t *qual_a, uint8_t site_sc, int64_t hap_cov_match, int64_t hap_cov_unmatch, uint32_t **oidx, uint32_t **o_a) { if(an <= 0) return; - int64_t *p, ri, rj, st, max_f, max_j, sc, plus = 0; int32_t *f, *ii; uint64_t cc = 0, cci, cc_min; + int64_t *p, ri, rj, st, s0, max_f, max_j, /**sc,**/ plus = 0; int32_t *f, *ii; uint64_t cc = 0, cci, cc_min; uint32_t rid; uint8_t rh; resize_Chain_Data(dp, an, NULL); idx->n = res->n = 0; f = dp->score; p = dp->pre; ii = dp->occ; cc = ((het_cov > 0)?(het_cov):(hom_cov/n_hap)); cc *= cut_rate; if(cc < cut_bd) cc = cut_bd; @@ -11405,12 +11438,13 @@ void gen_rphase_dp0_single_path_multi(SnpStats *a, int64_t an, haplotype_evdienc // fprintf(stderr, "\n\n\n[M::%s]\tcc::%lu\n\n\n", __func__, cc); // fprintf(stderr, "[M::%s]\tan::%ld\n", __func__, an); + /** for (ri = 0; ri < an; ++ri) { // fprintf(stderr, "[M::%s::i->%ld]\tsite::%u\tocc0::%u\tocc1::%u\tocc2::%u\n", __func__, ri, a[ri].site, a[ri].occ_0, a[ri].occ_1, a[ri].occ_2); max_f = 1; max_j = -1; st = 0; ///if(a[i].occ_0 < cc) st = i; for (rj = ri - 1; rj >= st; --rj) { - sc = comput_sc_rphase(&a[ri], ri, &a[rj], rj, za, 0/**cc**/); + sc = comput_sc_rphase(&a[ri], ri, &a[rj], rj, za, 0); if (sc == INT64_MIN) continue; sc += f[rj]; if (sc > max_f) { @@ -11421,6 +11455,37 @@ void gen_rphase_dp0_single_path_multi(SnpStats *a, int64_t an, haplotype_evdienc f[ri] = max_f; p[ri] = max_j; ii[ri] = 0; if(f[ri] < plus) plus = f[ri]; } + **/ + for (ri = 0; ri < an; ++ri) { + max_f = 1; max_j = -1; st = 0; ii[ri] = 0; + + s0 = ol->list[za[a[ri].non_homopolymer_num].overlapID].x_pos_s; + for (rj = ri - 1; (rj >= st) && (a[rj].site >= s0); --rj) { + if (f[rj] + 1 <= max_f) continue; + if (ii[rj] == ri) continue; + + if (comput_sc_rphase(&a[ri], ri, &a[rj], rj, za, 0, &rid, &rh) == INT64_MIN) { + if(rid != ((uint32_t)-1)) { + label_skip_mm(a, rid, 1 - rh, oidx, o_a, ii, ri, rj, an, ol->list[0].x_id); + } + continue; + } + + ///f[rj] + 1 > max_f + max_f = f[rj] + 1; max_j = rj; + + ///f[rj] <= rj + 1 + ///f[max_f - 1] <= max_f -> f[st] <= max_f -> max(f[st]) = max_f + st = (max_f>=1)?(max_f-1):(0); + } + + // fprintf(stderr, "[M::%s::]\tst::%ld\trj::%ld\n", __func__, st, rj); + + f[ri] = max_f; p[ri] = max_j; ii[ri] = 0; + if(f[ri] < plus) plus = f[ri]; + } + + return; kv_resize(uint64_t, *res, ((uint64_t)an)); gen_rphase_dp0_single_path_hybrid_0_multi(a, an, za, f, p, ii, qual_a, res->a, idx, plus, cc, cci, cc_min, NULL, ((uint64_t)-1), ((uint64_t)-1), 0, n_hap, cut_rate, cut_bd, site_sc, hap_cov_match, hap_cov_unmatch); @@ -11535,12 +11600,16 @@ void get_wqual(uint64_t zid, uint64_t zpos, uint64_t zrev, asg8_v *v, uint8_t *v } -void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, double st_rate, uint64_t st_max, asg64_v *idx, asg64_v *res, uint64_t rid, uint8_t *qa, uint64_t tcut, int64_t hap_cov_match, int64_t hap_cov_unmatch) +void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, double st_rate, uint64_t st_max, asg64_v *idx, asg64_v *res, uint64_t rid, uint8_t *qa, uint64_t tcut, int64_t hap_cov_match, int64_t hap_cov_unmatch, asg32_v *b32, + uint32_t **ro_idx, uint32_t **ro_a) { + ro_idx[0] = ro_idx[1] = ro_a[0] = ro_a[1] = NULL; if(hl->length <= 0) return; uint64_t k, l, i, i0, o, ii, m_snp_stat, m_snp_stat0, m_list, m1, c0, c1, rev_n, tqual, wqual, hq_cut = 2; uint32_t hq[5], hp[4], is_st; SnpStats *s; haplotype_evdience ev; char mc; - uint16_t fl = ((res && idx)?(1):(0)); + uint16_t fl = ((res && idx)?(1):(0)); uint32_t *r_idx[2], *r_a[2], *rc = NULL, *rz = NULL, *rs = NULL, rl[2]; rl[0] = rl[1] = 0; + b32->n = (ol->length + 1)<<1; kv_resize(uint32_t, *b32, b32->n); memset(b32->a, 0, sizeof((*(b32->a)))*b32->n); + r_idx[0] = b32->a; r_idx[1] = b32->a + ol->length + 1; if(fl) { @@ -11553,13 +11622,24 @@ void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, doub assert(i < hl->length && hl->list[i].site == hl->snp_stat.a[l].site); for (o = l, m1 = 0; o < k; o++) { s = &(hl->snp_stat.a[o]); + // if(rid == 307 && hl->snp_stat.a[l].site == 3592) { + // fprintf(stderr, "s[M::%s]\tsite::%u\tn0::%u\tn1::%u\tbs::%u\n", __func__, + // s->site, s->occ_0, s->occ_1, is_st_bs((*s), st_rate, st_max)); + // } // fprintf(stderr, "+[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, s->site, s->occ_0, s->occ_1); if((s->occ_0 < 2 || s->occ_1 < 2) || (is_st_bs((*s), st_rate, st_max)) || (!(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch))) { continue; } + // if(rid == 307 && hl->snp_stat.a[l].site == 3592) { + // fprintf(stderr, "b[M::%s]\tsite::%u\tn0::%u\tn1::%u\tbs::%u\n", __func__, + // s->site, s->occ_0, s->occ_1, is_st_bs((*s), st_rate, st_max)); + // } // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, s->site, s->occ_0, s->occ_1); m1++; } + // if(rid == 307 && hl->snp_stat.a[l].site == 3592) { + // fprintf(stderr, "+[M::%s]\tsite::%u\tn0::%u\tn1::%u\tm1::%lu\n", __func__, hl->snp_stat.a[l].site, hl->snp_stat.a[l].occ_0, hl->snp_stat.a[l].occ_1, m1); + // } if(m1) { hq[0] = hq[1] = hq[2] = hq[3] = hq[4] = 0; hp[0] = hp[1] = hp[2] = hp[3] = ((uint32_t)-1); @@ -11582,22 +11662,22 @@ void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, doub } assert(i0 > i); if(qa) { - if(hp[0] != ((uint32_t)-1) && hq[0] < hq_cut) { + if(hp[0] != ((uint32_t)-1) && hq[0] < hq_cut) {///disable SNP // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, hl->snp_stat.a[hp[0]].site, hl->snp_stat.a[hp[0]].occ_0, hl->snp_stat.a[hp[0]].occ_1); hl->snp_stat.a[hp[0]].occ_1 = 0; } - if(hp[1] != ((uint32_t)-1) && hq[1] < hq_cut) { + if(hp[1] != ((uint32_t)-1) && hq[1] < hq_cut) {///disable SNP // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, hl->snp_stat.a[hp[1]].site, hl->snp_stat.a[hp[1]].occ_0, hl->snp_stat.a[hp[1]].occ_1); hl->snp_stat.a[hp[1]].occ_1 = 0; } - if(hp[2] != ((uint32_t)-1) && hq[2] < hq_cut) { + if(hp[2] != ((uint32_t)-1) && hq[2] < hq_cut) {///disable SNP // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, hl->snp_stat.a[hp[2]].site, hl->snp_stat.a[hp[2]].occ_0, hl->snp_stat.a[hp[2]].occ_1); hl->snp_stat.a[hp[2]].occ_1 = 0; } - if(hp[3] != ((uint32_t)-1) && hq[3] < hq_cut) { + if(hp[3] != ((uint32_t)-1) && hq[3] < hq_cut) {///disable SNP // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, hl->snp_stat.a[hp[3]].site, hl->snp_stat.a[hp[3]].occ_0, hl->snp_stat.a[hp[3]].occ_1); hl->snp_stat.a[hp[3]].occ_1 = 0; } @@ -11615,11 +11695,23 @@ void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, doub // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\thq[0]::%u\thq[1]::%u\thq[2]::%u\thq[3]::%u\thq[4]::%u\thq_cut::%lu\n", __func__, s->site, s->occ_0, s->occ_1, hq[0], hq[1], hq[2], hq[3], hq[4], hq_cut); if((hq[4] >= hq_cut) && (hq[0] >= hq_cut || hq[1] >= hq_cut || hq[2] >= hq_cut || hq[3] >= hq_cut)) { + // if(rid == 307 && hl->snp_stat.a[l].site == 3592) { + // fprintf(stderr, "-[M::%s]\tsite::%u\tn0::%u\tn1::%u\tbs::%u\thq[0]::%u\thq[1]::%u\thq[2]::%u\thq[3]::%u\thq[4]::%u\thq_cut::%lu\n", __func__, hl->snp_stat.a[l].site, hl->snp_stat.a[l].occ_0, hl->snp_stat.a[l].occ_1, is_st_bs(hl->snp_stat.a[l], st_rate, st_max), + // hq[0], hq[1], hq[2], hq[3], hq[4], hq_cut); + // } for (o = l, m_snp_stat0 = m_snp_stat; o < k; o++) { s = &(hl->snp_stat.a[o]); + // if(rid == 307 && hl->snp_stat.a[l].site == 3592) { + // fprintf(stderr, "+[M::%s]\tsite::%u\tn0::%u\tn1::%u\tbs::%u\n", __func__, + // s->site, s->occ_0, s->occ_1, is_st_bs((*s), st_rate, st_max)); + // } if((s->occ_0 < 2 || s->occ_1 < 2) || (is_st_bs((*s), st_rate, st_max)) || (!(s->occ_0 >= hap_cov_match && s->occ_1 >= hap_cov_unmatch))) { continue; } + // if(rid == 307 && hl->snp_stat.a[l].site == 3592) { + // fprintf(stderr, ">[M::%s]\tsite::%u\tn0::%u\tn1::%u\tbs::%u\tm_snp_stat::%lu\n", __func__, + // s->site, s->occ_0, s->occ_1, is_st_bs((*s), st_rate, st_max), m_snp_stat); + // } hl->snp_stat.a[m_snp_stat] = *s; for (ii = i, c1 = 0; ii < i0; ii++) { @@ -11637,22 +11729,31 @@ void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, doub hl->snp_stat.a[m_snp_stat++].score = -1; } - for (ii = i, m1 = m_list; ii < i0; ii++) { - if((hl->list[ii].site == ((uint32_t)-1))) { - hl->list[ii].overlapSite = ((uint32_t)-1); - hl->list[ii].site = hl->snp_stat.a[l].site; + if(m_snp_stat > m_snp_stat0) { + for (ii = i, m1 = m_list; ii < i0; ii++) { + if((hl->list[ii].site == ((uint32_t)-1))) { + hl->list[ii].overlapSite = ((uint32_t)-1); + hl->list[ii].site = hl->snp_stat.a[l].site; + } + if(fl) fill_incom(res, hl->list[ii].overlapID, hl->list[ii].site, idx->a, idx->n, ol->list[0].x_id); + if(hh_tp(hl->list[ii]) == 0 || hh_tp(hl->list[ii]) == 1) { + // if(rid == 307 && hl->list[ii].overlapID == 111 && hh_tp(hl->list[ii]) == 0) { + // fprintf(stderr, "-0-[M::%s]\ti::%lu\tsite::%u\thtp::%u\tsnp_stat->site::%u\n", __func__, m_list, hl->list[ii].site, hh_tp(hl->list[ii]), hl->snp_stat.a[l].site); + // // fprintf(stderr, "-0-[M::%s]\tsite::%u\n", __func__, hl->list[ii].site); + // } + r_idx[hh_tp(hl->list[ii])][hl->list[ii].overlapID]++; rl[hh_tp(hl->list[ii])]++; + } + hl->list[m_list++] = hl->list[ii]; } - if(fl) fill_incom(res, hl->list[ii].overlapID, hl->list[ii].site, idx->a, idx->n, ol->list[0].x_id); - hl->list[m_list++] = hl->list[ii]; - } - radix_sort_haplotype_evdience_id_srt(hl->list + m1, hl->list + m_list); - for (ii = m_snp_stat0; ii < m_snp_stat; ii++) { - hl->snp_stat.a[ii].homopolymer_num = m_list; - } - - i = i0; + radix_sort_haplotype_evdience_id_srt(hl->list + m1, hl->list + m_list); + for (ii = m_snp_stat0; ii < m_snp_stat; ii++) { + hl->snp_stat.a[ii].homopolymer_num = m_list; + } + + i = i0; - if(fl) kv_push(uint64_t, *idx, hl->snp_stat.a[l].site); + if(fl) kv_push(uint64_t, *idx, hl->snp_stat.a[l].site); + } } } @@ -11662,7 +11763,18 @@ void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, doub hl->snp_stat.n = m_snp_stat; hl->length = m_list; if(hl->snp_stat.n == 0 || hl->length == 0) return; - + // if(rid == 307) { + // for (ii = 0; ii < hl->length; ii++) { + // if(hl->list[ii].overlapID == 111) { + // fprintf(stderr, "-a-[M::%s]\ti::%lu\tsite::%u\thtp::%u\n", __func__, m_list, hl->list[ii].site, hh_tp(hl->list[ii])); + // // fprintf(stderr, "-0-[M::%s]\tsite::%u\n", __func__, hl->list[ii].site); + // } + // } + // for (ii = 0; ii < hl->snp_stat.n; ++ii) { + // fprintf(stderr, "-x-[M::%s]\tsite::%u\tocc_0::%u\tocc_1::%u\n", __func__, hl->snp_stat.a[ii].site, hl->snp_stat.a[ii].occ_0, hl->snp_stat.a[ii].occ_1); + // } + // } + // fprintf(stderr, "+[M::%s]\tres->n::%lu\tol->length::%lu\tsnp_stat.n::%lu\n", __func__, (uint64_t)res->n, ol->length, (uint64_t)hl->snp_stat.n); ///fill incomplete values if((fl) && (res->n > ol->length)) { @@ -11676,6 +11788,11 @@ void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, doub ev.cov = 1; addHaplotypeEvdience(hl, &ev, NULL); // fprintf(stderr, "[M::%s]\toid::%u\tsite::%u\n", __func__, ev.overlapID, ev.site); + r_idx[hh_tp(ev)][ev.overlapID]++; + rl[hh_tp(ev)]++; + // if(rid == 307 && ev.overlapID == 111 && hh_tp(ev) == 0) { + // fprintf(stderr, "-0-[M::%s]\tsite::%u\tsup\n", __func__, ev.site); + // } } radix_sort_haplotype_evdience_srt(hl->list, hl->list + hl->length); @@ -11705,6 +11822,7 @@ void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, doub for (ii = i, is_st = 1; ii < i0; ii++) { if(hh_tp(hl->list[ii]) == 0) { assert(hl->list[ii].overlapSite == ((uint32_t)-1) || hl->list[ii].overlapSite == m1); + // if(rid == 307) fprintf(stderr, "-1-[M::%s]\tm1::%lu\n", __func__, m1); hl->list[ii].overlapSite = m1; hl->list[ii].misBase = mc; } if((ii > i) && (hl->list[ii].overlapID <= hl->list[ii-1].overlapID)) { @@ -11735,6 +11853,76 @@ void call_rphase_sc(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, doub } } } + + // if(rid == 307) { + // for (ii = 0; ii < hl->length; ii++) { + // if(hl->list[ii].overlapID == 111) { + // fprintf(stderr, "-b-[M::%s]\ti::%lu\tsite::%u\thtp::%u\n", __func__, m_list, hl->list[ii].site, hh_tp(hl->list[ii])); + // // fprintf(stderr, "-0-[M::%s]\tsite::%u\n", __func__, hl->list[ii].site); + // } + // } + // for (ii = 0; ii < hl->snp_stat.n; ++ii) { + // fprintf(stderr, "-y-[M::%s]\tsite::%u\tocc_0::%u\tocc_1::%u\n", __func__, hl->snp_stat.a[ii].site, hl->snp_stat.a[ii].occ_0, hl->snp_stat.a[ii].occ_1); + // } + // } + + if(ol->length > 0) { + b32->n += rl[0] + rl[1]; kv_resize(uint32_t, *b32, b32->n); + r_idx[0] = b32->a; r_idx[1] = r_idx[0] + ol->length + 1; + r_a[0] = r_idx[1] + ol->length + 1; r_a[1] = r_a[0] + rl[0]; + + for (k = l = o = 0; k < ol->length; k++) { + o = r_idx[0][k]; + if(o > 0) r_a[0][l + o - 1] = 0; + r_idx[0][k] = l; l += o; + } + r_idx[0][k] = l; assert(l == rl[0]); + + for (k = l = o = 0; k < ol->length; k++) { + o = r_idx[1][k]; + if(o > 0) r_a[1][l + o - 1] = 0; + r_idx[1][k] = l; l += o; + } + r_idx[1][k] = l; assert(l == rl[1]); + + uint64_t sbn[2] = {0, 0}; + + for (k = 1, l = 0, i = m_snp_stat = m_list = 0; k <= hl->snp_stat.n; ++k) {///filter snps + if(k == hl->snp_stat.n || hl->snp_stat.a[k].site != hl->snp_stat.a[l].site) { + for (; i < hl->length && hl->list[i].site != hl->snp_stat.a[l].site; i++); + assert(i < hl->length && hl->list[i].site == hl->snp_stat.a[l].site); + for (; i < hl->length && hl->list[i].site == hl->snp_stat.a[l].site; i++) { + if((hh_tp(hl->list[i]) != 0) && (hh_tp(hl->list[i]) != 1)) continue; + + // if(rid == 307 && hl->list[i].overlapID == 111 && hh_tp(hl->list[i]) == 0) { + // fprintf(stderr, "-2-[M::%s]\ti::%lu\tzss::%u\thtp::%u\tsnp_stat->site::%u\n", __func__, i, hl->list[i].site, hh_tp(hl->list[i]), hl->snp_stat.a[l].site); + // // fprintf(stderr, "-2-[M::%s]\tsite::%u\n", __func__, hl->list[i].site); + // } + + rc = r_idx[hh_tp(hl->list[i])]; rz = r_a[hh_tp(hl->list[i])]; + o = rc[hl->list[i].overlapID + 1] - rc[hl->list[i].overlapID]; + assert(o > 0); + rs = &(rz[rc[hl->list[i].overlapID] + o - 1]); + assert((*rs) < o); + if(o - (*rs) > 1) { + rz[rc[hl->list[i].overlapID] + (*rs)] = hl->list[i].overlapSite; + // if(rz[rc[hl->list[i].overlapID] + (*rs)] == ((uint32_t)-1)) rz[rc[hl->list[i].overlapID] + (*rs)] = k - 1; + (*rs)++; + } else { + (*rs) = hl->list[i].overlapSite; + // if((*rs) == ((uint32_t)-1)) (*rs) = k - 1; + } + sbn[hh_tp(hl->list[i])]++; + } + l = k; + } + } + + ro_idx[0] = r_idx[0]; ro_idx[1] = r_idx[1]; + ro_a[0] = r_a[0]; ro_a[1] = r_a[1]; + + // if(rid == 307) fprintf(stderr, "-1-[M::%s] rl[0]::%u, rl[1]::%u, sbn[0]::%lu, sbn[1]::%lu\n", __func__, rl[0], rl[1], sbn[0], sbn[1]); + } } @@ -11811,19 +11999,19 @@ void recal_rphase(All_reads *rref, haplotype_evdience_alloc *hl, overlap_region_ void gen_rphase_dp_adv(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, UC_Read* g_read, double st_rate, uint64_t st_max, Chain_Data *dp, asg64_v *idx, asg64_v *res, uint64_t rid, uint8_t *qa, uint64_t tcut, uint64_t site_sc, - int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t n_hap, int64_t het_c, int64_t hom_c, double hf_rate) + int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t n_hap, int64_t het_c, int64_t hom_c, double hf_rate, asg32_v *b32) { if(hl->length <= 0) return; - uint64_t k, l, i, i0, o, ii, m_snp_stat, m_snp_stat0, m_list /**m_off**/, m1, c0, c1, rev_n; SnpStats *s; + uint64_t k, l, i, i0, o, ii, m_snp_stat, m_snp_stat0, m_list /**m_off**/, m1, c0, c1, rev_n; SnpStats *s; uint32_t *oidx[2], *o_a[2]; - call_rphase_sc(hl, ol, st_rate, st_max, idx, res, rid, qa, tcut, hap_cov_match, hap_cov_unmatch); + call_rphase_sc(hl, ol, st_rate, st_max, idx, res, rid, qa, tcut, hap_cov_match, hap_cov_unmatch, b32, oidx, o_a); // fprintf(stderr, "+[M::%s]\tyid::%u\thap_cov_match::%ld\thap_cov_unmatch::%ld\tn_hap::%ld\thet_c::%ld\thom_c::%ld\thf_rate::%f\n", // __func__, ol->list[0].y_id, hap_cov_match, hap_cov_unmatch, n_hap, het_c, hom_c, hf_rate); if(tcut == ((uint64_t)-1)) { // gen_rphase_dp0_single_path(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_a, hom_a, n_hap, 0.7, 6, res, qv->a, site_sc); - gen_rphase_dp0_single_path_multi(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_c, hom_c, n_hap, 0.6/**0.7**/, 6, res, qa, site_sc, hap_cov_match, hap_cov_unmatch); + gen_rphase_dp0_single_path_multi(ol, hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_c, hom_c, n_hap, 0.6/**0.7**/, 6, res, qa, site_sc, hap_cov_match, hap_cov_unmatch, oidx, o_a); } else { gen_rphase_dp0_single_path_hybrid(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_c, hom_c, n_hap, 0.6/**0.7**/, 0.4, 6, res, qa, ol, rid, tcut, site_sc, hap_cov_match, hap_cov_unmatch, hf_rate); } @@ -11899,7 +12087,7 @@ void gen_rphase_dp(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, UC_Re { if(hl->length <= 0) return; - uint64_t k, l, i, i0, o, ii, m_snp_stat, m_snp_stat0, m_list /**m_off**/, m1, c0, c1, rev_n, tqual, wqual, hq_cut = 2; uint32_t hq[5], hp[4], is_st; SnpStats *s; haplotype_evdience ev; char mc; + uint64_t k, l, i, i0, o, ii, m_snp_stat, m_snp_stat0, m_list /**m_off**/, m1, c0, c1, rev_n, tqual, wqual, hq_cut = 2; uint32_t hq[5], hp[4], is_st; SnpStats *s; haplotype_evdience ev; char mc; uint32_t *oidx[2], *o_a[2]; if(rid < tcut) { retrive_bqual(qv, NULL, rid, -1, -1, 0, sc_bn); } else { @@ -12103,7 +12291,7 @@ void gen_rphase_dp(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, UC_Re if(tcut == ((uint64_t)-1)) { // gen_rphase_dp0_single_path(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_a, hom_a, n_hap, 0.7, 6, res, qv->a, site_sc); - gen_rphase_dp0_single_path_multi(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_c, hom_c, n_hap, 0.7, 6, res, qv->a, site_sc, hap_cov_match, hap_cov_unmatch); + gen_rphase_dp0_single_path_multi(ol, hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_c, hom_c, n_hap, 0.7, 6, res, qv->a, site_sc, hap_cov_match, hap_cov_unmatch, oidx, o_a); } else { gen_rphase_dp0_single_path_hybrid(hl->snp_stat.a, hl->snp_stat.n, hl->list, dp, idx, het_c, hom_c, n_hap, 0.7, 0.4, 6, res, qv->a, ol, rid, tcut, site_sc, hap_cov_match, hap_cov_unmatch, hf_rate); } @@ -15354,6 +15542,50 @@ uint32_t align_hc_ed_post_extz(overlap_region *z, All_reads *rref, char* qstr, c return 1; } + +uint32_t align_gc_ed_post_extz(overlap_region *z, ul_idx_t *uref, char* qstr, char *tstr, bit_extz_t *exz, double e_rate, int64_t w_l, double ovlp_cut, int64_t force_aln, uint64_t *tot_b) +{ + int64_t q_s, q_e, nw, k, q_l, t_l, t_tot_l, aux_beg, aux_end, t_s, thre, aln_l, t_pri_l; ///uint64_t tot_b0 = 0; + char *q_string, *t_string; + z->w_list.n = 0; z->is_match = 0; z->align_length = 0; + nw = get_num_wins(z->x_pos_s, z->x_pos_e+1, w_l); + get_win_se_by_normalize_xs(z, (z->x_pos_s/w_l)*w_l, w_l, &q_s, &q_e); + + for (k = 0; k < nw; k++) { + aux_beg = aux_end = 0; q_l = 1 + q_e - q_s; + thre = q_l*e_rate; thre = Adjust_Threshold(thre, q_l); + if(thre > THRESHOLD_MAX_SIZE) thre = THRESHOLD_MAX_SIZE; + ///offset of y + t_s = (q_s - z->x_pos_s) + z->y_pos_s; + t_s += y_start_offset(q_s, &(z->f_cigar)); + + aln_l = q_l + (thre<<1); t_tot_l = uref->ug->u.a[z->y_id].len;; + if(init_waln(thre, t_s, t_tot_l, aln_l, &aux_beg, &aux_end, &t_s, &t_pri_l)) { + q_string = qstr+q_s; + + t_string = return_str_seq_exz(tstr, t_s, t_pri_l, z->y_pos_strand, NULL, uref, z->y_id); + t_l = t_pri_l; + + ed_band_cal_semi_64_w_absent_diag(t_string, t_l, q_string, q_l, thre, aux_beg, exz); + + if(tot_b) (*tot_b) += q_l; + + if (is_align(*exz)) { + ///t_s do not have aux_beg, while t_s + t_end (aka, te) has + if(!push_hc_wlst_exz(uref, NULL, NULL, z, qstr, tstr, exz, THRESHOLD_MAX_SIZE, q_s, q_e, t_s, t_s + exz->pe, + t_tot_l, aux_beg, aux_end, e_rate, w_l, ovlp_cut, force_aln, tot_b, 0)) { + return 0; + } + } + } + q_s = q_e + 1; q_e = q_s + w_l - 1; + if(q_e >= (int64_t)z->x_pos_e) q_e = z->x_pos_e; + } + + if((!force_aln) && (!pass_qovlp(z->x_pos_e+1-z->x_pos_s, z->align_length, ovlp_cut))) return 0; + return 1; +} + uint32_t align_hc_ed_post_extz_cut(overlap_region *z, overlap_region *zr, /**overlap_region *zuf,**/ All_reads *rref, char* qstr, char *tstr, bit_extz_t *exz, double e_rate, int64_t w_l, double ovlp_cut, int64_t force_aln, uint32_t *e_win, double nhit_cut, uint64_t *tot_b) { int64_t q_s, q_e, nw, k, kr = 0, nr = 0, /**nuf = 0, fqs = -1, fqe = -1,**/ q_l, t_l, t_tot_l, aux_beg, aux_end, t_s, thre, aln_l, t_pri_l, zol = z->x_pos_e+1-z->x_pos_s, zal, rual = 0; @@ -18064,6 +18296,7 @@ All_reads *rref, int64_t id) } if(uref) { retrieve_u_seq(NULL, str, &(uref->ug->u.a[id]), rev, ss, sl, NULL); + // retrieve_u_seq_fast(NULL, str, &(uref->ug->u.a[id]), rev, ss, sl, NULL); } else if(rref) { recover_UC_Read_sub_region(str, ss, sl, rev, rref, id); } @@ -20221,6 +20454,89 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u } } +void gc_ovlp_base_direct(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, int64_t wl, ul_idx_t *udb, 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, int64_t pre_mode) +{ + int64_t i, l, mode, q[2], t[2], qr, tr, is_done, zn, si, ei; + + 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) && + (z->w_list.a[i].y_end == (z->w_list.a[i-1].y_end + (z->w_list.a[i].x_end-z->w_list.a[i-1].x_end)))) { + continue; + } + break; + } + if(i >= zn) { + q[0] = z->w_list.a[0].x_start; q[1] = z->w_list.a[z->w_list.n-1].x_end; + t[1] = z->w_list.a[z->w_list.n-1].y_end; t[0] = z->w_list.a[0].y_end - (z->w_list.a[0].x_end-z->w_list.a[0].x_start); + + if(q[0] <= t[0]) { + t[0] -= q[0]; q[0] = 0; + } else { + q[0] -= t[0]; t[0] = 0; + } + + qr = ql-q[1]-1; tr = tl-t[1]-1; + if(qr <= tr) { + q[1] = ql-1; t[1] += qr; + } else { + t[1] = tl-1; q[1] += tr; + } + + if(q[0] == z->w_list.a[0].x_start && q[1] == z->w_list.a[z->w_list.n-1].x_end) { + // fprintf(stderr, "[M::%s::%u->%u::%c] ovlp::%u, w_list.n::%u\n", __func__, z->x_id, z->y_id+1, "+-"[z->y_pos_strand], z->x_pos_e+1-z->x_pos_s, (uint32_t)z->w_list.n); + set_exact_exz(exz, q[0], q[1] + 1, t[0], t[1] + 1); push_alnw(aux_o, exz); + return; + } + } + } + + 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; + } else { + q[0] = 0; + } + + if(i < ch_n) { + q[1] = ch_a[i].self_offset; t[1] = ch_a[i].offset; + } else { + q[1] = ql; + } + + if((t[0] != -1) && (t[1] != -1)) { + mode = 0;//global + } else if((t[0] != -1) && (t[1] == -1)) { + mode = 1;///forward extension + } else if((t[0] == -1) && (t[1] != -1)) { + mode = 2;///backward extension + } else { + mode = 3;///no primary hit within [ibeg, iend] + } + + if(mode == 1 || mode == 2) adjust_ext_offset(&(q[0]), &(q[1]), &(t[0]), &(t[1]), ql, tl, 0, mode); + is_done = hc_aln_exz_adv_hc(z, udb, NULL, NULL, qstr, tu, q[0], q[1], t[0], t[1], mode, wl, exz, ql, e_rate, + MAX_SIN_L, MAX_SIN_E, FORCE_SIN_L, -1, aux_o); + + if(!is_done) {///postprocess + push_unmap_alnw(aux_o, q[0], q[1]-1, t[0], t[1]-1, mode); + } + l = i; + } +} + 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, int64_t pre_mode) @@ -20457,6 +20773,78 @@ UC_Read *tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, // } } +void rechain_aln_gc(overlap_region *z, Candidates_list *cl, overlap_region *aux_o, int64_t aux_i, int64_t wl, + ul_idx_t *udb, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate, int64_t ql, int64_t tl, int64_t h_khit, int64_t rid) +{ + int64_t rcn = cl->length, ch_n, qs, qe, ts, te, mode, an0, an, todo; + k_mer_hit *ch_a; uint8_t q[2], t[2]; ///ul_ov_t idx; + ///[qs, qe) && [ts, te) + qs = aux_o->w_list.a[aux_i].x_start; qe = aux_o->w_list.a[aux_i].x_end+1; + ts = aux_o->w_list.a[aux_i].y_start; te = aux_o->w_list.a[aux_i].y_end+1; + if(qe - qs < FORCE_SIN_L || te - ts < FORCE_SIN_L) return; + mode = aux_o->w_list.a[aux_i].error_threshold; + ch_n = gen_win_chain(z, cl, qs, qe, ts, te, wl, udb, NULL, NULL, qstr, tu, exz, ql, tl, e_rate, h_khit, mode, rid, 1); + ch_a = cl->list + rcn; + if(ch_n) { + todo = 1; ///idx.ts = idx.te = (uint32_t)-1; idx.qs = 0; idx.qe = ql; + if(mode == 0) {//global + // idx.qn = 0; idx.tn = ch_n - 1; + // idx.qs = ch_a[idx.qn].self_offset; + // idx.ts = ch_a[idx.qn].offset; + // idx.qe = ch_a[idx.tn].self_offset; + // idx.te = ch_a[idx.tn].offset; + assert(ch_a[0].self_offset == qs && ch_a[0].offset == ts); + assert(ch_a[ch_n-1].self_offset == qe && ch_a[ch_n-1].offset == te); + if(ch_n <= 2) todo = 0; + } else if(mode == 1) {//forward ext + // idx.qn = 0; idx.tn = ch_n; + // idx.qs = ch_a[idx.qn].self_offset; + // idx.ts = ch_a[idx.qn].offset; + // idx.qe = ql; + assert(ch_a[0].self_offset == qs && ch_a[0].offset == ts); + if(ch_n <= 1) todo = 0; + } else if(mode == 2) {///backward ext + // idx.qn = (uint32_t)-1; idx.tn = ch_n-1; + // idx.qs = 0; + // idx.qe = ch_a[idx.tn].self_offset; + // idx.te = ch_a[idx.tn].offset; + assert(ch_a[ch_n-1].self_offset == qe && ch_a[ch_n-1].offset == te); + if(ch_n <= 1) todo = 0; + } + if(todo) { + an0 = aux_o->w_list.n; + // if(z->x_id == 29033 && z->y_id == 21307) { + // fprintf(stderr, "[M::%s]\tan0::%ld\tq::[%u,\t%u)\tt::[%u,\t%u)\tlw::%u\trw::%u\n", __func__, an0, + // 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); + gc_ovlp_base_direct(z, ch_a, ch_n, wl, udb, 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); + // } + // fprintf(stderr, "[M::%s::] awn0::%ld, awn::%lu\n", __func__, an0, an); + ///old unaligned window could be replaced by the new aligned window + if((an == (an0 + 1)) && (!(is_ualn_win(aux_o->w_list.a[an-1])))) { + if(aux_o->w_list.a[aux_i].x_start == aux_o->w_list.a[an-1].x_start) q[0] = 1; + if(aux_o->w_list.a[aux_i].x_end == aux_o->w_list.a[an-1].x_end) q[1] = 1; + if(aux_o->w_list.a[aux_i].y_start == aux_o->w_list.a[an-1].y_start) t[0] = 1; + if(aux_o->w_list.a[aux_i].y_end == aux_o->w_list.a[an-1].y_end) t[1] = 1; + if((mode == 0) && q[0] && q[1] && t[0] && t[1]) todo = 1; + if((mode == 1) && q[0] && t[0]) todo = 1; + if((mode == 2) && q[1] && t[1]) todo = 1; + if(todo) { + aux_o->w_list.a[aux_i] = aux_o->w_list.a[an-1]; aux_o->w_list.n--; + } + } + // if(an > an0) {///should always > 0 as there are unmapped windows + // } + // aux_o->w_list.n = an0; + } + } + cl->length = rcn;///must reset!!!! +} + void rechain_aln_hc(overlap_region *z, Candidates_list *cl, overlap_region *aux_o, int64_t aux_i, int64_t wl, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, double e_rate, int64_t ql, int64_t tl, int64_t h_khit, int64_t rid) { @@ -20618,6 +21006,78 @@ void push_trace_iter(asg16_v *exz, uint16_t c, uint32_t l) } **/ +/** +uint64_t gen_gc_fast_cigar0(overlap_region *z, Candidates_list *cl, uint64_t wl, ul_idx_t *udb, char* qstr, UC_Read *tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, uint64_t rid, int64_t h_khit, int64_t *re) +{ + return_t_chain(z, cl); + int64_t ch_idx = z->shared_seed, ch_n; + int64_t i, tl, id = z->y_id, m, tot_e, aln, xe, ye; + k_mer_hit *ch_a = cl->list + ch_idx; *re = INT64_MAX; + 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(ch_n <= 0) return 0; + + 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; + 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; + + gc_ovlp_base_direct(z, ch_a, ch_n, wl, udb, qstr, tu, exz, aux_o, e_rate, ql, tl, rid, -1); + + int64_t aux_n = aux_o->w_list.n; + if(z->w_list.n == 0) aux_n = 0; + for (i = 0; i < aux_n; i++) { + if(!(is_ualn_win(aux_o->w_list.a[i]))) continue; + //will overwrite ch_a; does not matter + rechain_aln_gc(z, cl, aux_o, i, wl, udb, qstr, tu, exz, e_rate, ql, tl, h_khit, rid); + } + + if(((int64_t)aux_o->w_list.n) > aux_n) { + for (i = m = 0; i < ((int64_t)aux_o->w_list.n); i++) { + if((i < aux_n) && (is_ualn_win(aux_o->w_list.a[i]))) continue; + aux_o->w_list.a[m++] = aux_o->w_list.a[i]; + } + aux_o->w_list.n = m; + radix_sort_window_list_xs_srt(aux_o->w_list.a, aux_o->w_list.a+aux_o->w_list.n); + } + + ///update z by aux_o + update_overlap_region(z, aux_o, ql, tl); + + aux_n = z->w_list.n; + for (i = tot_e = aln = 0; i < aux_n; i++) { + if(is_ualn_win(z->w_list.a[i])) { + xe = z->w_list.a[i].x_end + 1 - z->w_list.a[i].x_start; + ye = z->w_list.a[i].y_end + 1 - z->w_list.a[i].y_start; + tot_e += ((xe >= ye)?(xe):(ye)); + } else { + tot_e += z->w_list.a[i].error; aln += z->w_list.a[i].x_end + 1 - z->w_list.a[i].x_start; + } + } + *re = tot_e; + // fprintf(stderr, "[M::%s::%u->%u::%c] ovlp::%u, aln::%ld, tot_e::%ld, w_list.n::%u, ch_n::%ld\n", + // __func__, z->x_id, z->y_id+1, "+-"[z->y_pos_strand], z->x_pos_e+1-z->x_pos_s, aln, tot_e, (uint32_t)z->w_list.n, ch_n); + + // debug_overlap_region(aux_o, qstr, tu, NULL, NULL, rref); + + + // ch_a = cl->list + ch_idx; //update + // for (i = 0; i < wn; i++) z->w_list.a[i].clen = 0;///clean cigar + // if(on > 1) { + // fprintf(stderr, "[M::%s::] rid::%lu, on::%ld\n", __func__, rid, on); + // } + // if(z->y_id == 126) prt_k_mer_hit(ch_a, ch_n); + // for (i = ch_i = 0; i < on; i++) { + // assert((i<=0)||(ov[i].qs > ov[i-1].qe)); + // ov[i].sec = 16;///do not know the aln type + // ch_i = sub_base_aln(z, dp, ch_a, ch_n, pe, ov[i].qs, ov[i].qe, wl, uref, hpc_g, rref, qstr, tu, exz, e_rate, ql, tl, ch_i, rid); + // pe = ov[i].qe; + // } + + return 1; +} +**/ + uint64_t gen_hc_fast_cigar0(overlap_region *z, Candidates_list *cl, uint64_t wl, All_reads *rref, char* qstr, UC_Read *tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, uint64_t rid, int64_t h_khit, int64_t *re) { int64_t ch_idx = z->shared_seed, ch_n; @@ -20707,11 +21167,11 @@ uint64_t gen_hc_fast_cigar0(overlap_region *z, Candidates_list *cl, uint64_t wl, } -uint64_t gen_hc_fast_non_retrieve_cigar0(overlap_region *z, Candidates_list *cl, uint64_t wl, char* qstr, int64_t ql, char* tstr, int64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, uint64_t rid, int64_t h_khit) +uint64_t gen_hc_fast_non_retrieve_cigar0(overlap_region *z, Candidates_list *cl, uint64_t wl, char* qstr, int64_t ql, char* tstr, int64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, uint64_t rid, int64_t h_khit, int64_t *re) { int64_t ch_idx = z->shared_seed, ch_n; - int64_t i, m, tot_e, aln; - k_mer_hit *ch_a = cl->list + ch_idx; + int64_t i, m, tot_e, aln, xe, ye; + k_mer_hit *ch_a = cl->list + ch_idx; *re = INT64_MAX; for (i = ch_idx; i < cl->length && cl->list[i].readID == cl->list[ch_idx].readID; i++){;} ch_n = i-ch_idx; if(ch_n <= 0) return 0; @@ -20756,11 +21216,14 @@ uint64_t gen_hc_fast_non_retrieve_cigar0(overlap_region *z, Candidates_list *cl, aux_n = z->w_list.n; for (i = tot_e = aln = 0; i < aux_n; i++) { if(is_ualn_win(z->w_list.a[i])) { - tot_e += z->w_list.a[i].x_end + 1 - z->w_list.a[i].x_start; + xe = z->w_list.a[i].x_end + 1 - z->w_list.a[i].x_start; + ye = z->w_list.a[i].y_end + 1 - z->w_list.a[i].y_start; + tot_e += ((xe >= ye)?(xe):(ye)); } else { tot_e += z->w_list.a[i].error; aln += z->w_list.a[i].x_end + 1 - z->w_list.a[i].x_start; } } + *re = tot_e; // fprintf(stderr, "[M::%s::%u->%u::%c] ovlp::%u, aln::%ld, tot_e::%ld, w_list.n::%u, ch_n::%ld\n", // __func__, z->x_id, z->y_id+1, "+-"[z->y_pos_strand], z->x_pos_e+1-z->x_pos_s, aln, tot_e, (uint32_t)z->w_list.n, ch_n); @@ -24507,6 +24970,7 @@ void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_all } hp->length = t; + if(q8) gen_qvec_hvec(rref, q8, &qhf, ((dp)?(&qual):(NULL)), rid, tcut); //r829 @@ -24520,7 +24984,9 @@ void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_all // generate_haplotypes_naive(hap, overlap_list, R_INF, g_read->length, force_repeat); } else { // gen_rphase_dp(hp, ol, qu, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), dp, idx, buf, rid, q8, tcut, site_sc); - gen_rphase_dp_adv(hp, ol, qu, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), dp, idx, buf, rid, qual, tcut, site_sc, hap_cov_match, hap_cov_unmatch, n_hap, het_cov_a, hom_cov_a, hf_rate); + gen_rphase_dp_adv(hp, ol, qu, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), dp, idx, buf, rid, qual, tcut, site_sc, hap_cov_match, hap_cov_unmatch, n_hap, het_cov_a, hom_cov_a, hf_rate, b32); + + return; // if(!site_sc) generate_haplotypes_naive_HiFi_adv(hp, ol, 0.04, qu, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), b32, rid);//r835 if(!site_sc) generate_haplotypes_naive_HiFi_adv_hc(hp, ol, 0.04, qu, ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), b32, rid, hap_cov_match, hap_cov_unmatch);///r835 else generate_haplotypes_weight(hp, ol, 0.04, qu, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), 32, 2, hap_cov_match, hap_cov_unmatch); @@ -29857,7 +30323,6 @@ void ul_rid_lalign_adv(overlap_region_alloc* ol, Candidates_list *cl, const ul_i } } - uint64_t gen_hc_fast_cigar(overlap_region *z, Candidates_list *cl, All_reads *rref, int64_t wl, char *qstr, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t rid, int64_t khit, int64_t *re) { return_t_chain(z, cl); @@ -29865,10 +30330,10 @@ uint64_t gen_hc_fast_cigar(overlap_region *z, Candidates_list *cl, All_reads *rr return 1; } -uint64_t gen_hc_fast_non_retrieve_cigar(overlap_region *z, Candidates_list *cl, int64_t wl, char *qstr, int64_t ql, char *tstr, int64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t rid, int64_t khit) +uint64_t gen_hc_fast_non_retrieve_cigar(overlap_region *z, Candidates_list *cl, int64_t wl, char *qstr, int64_t ql, char *tstr, int64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t rid, int64_t khit, int64_t *re) { return_t_chain(z, cl); - gen_hc_fast_non_retrieve_cigar0(z, cl, wl, qstr, ql, tstr, tl, exz, aux_o, e_rate, rid, khit); + gen_hc_fast_non_retrieve_cigar0(z, cl, wl, qstr, ql, tstr, tl, exz, aux_o, e_rate, rid, khit, re); return 1; } @@ -32108,6 +32573,7 @@ void reassign_gaps(overlap_region *z, overlap_region *aux, char* qstr, int64_t q } + uint32_t is_ovlp_debug(int64_t s, int64_t e, int64_t ws, int64_t we, int64_t op) { int64_t os, oe, ovlp; @@ -33628,8 +34094,7 @@ uint8_t inline gen_hc_r_alin_flt_1_smp(overlap_region *z, Candidates_list *cl, A // z->is_match, z->non_homopolymer_errors, rr, err, re, f); // } } - - + if(f) { // if(z->y_id == 3621141) { @@ -33686,6 +34151,57 @@ uint8_t inline gen_hc_r_alin_flt_1_smp(overlap_region *z, Candidates_list *cl, A } +uint8_t inline gen_gc_r_alin_flt_1_smp(overlap_region *z, Candidates_list *cl, char *qstr, int64_t ql, char *tstr, int64_t tl, bit_extz_t *exz, overlap_region *aux_o, + double err, double e_max, double e_rate, int64_t wsl, int64_t rid, int64_t khit, uint8_t chem_drop, double align_gap_rate, int64_t align_gap_max, uint8_t *hpf, asg16_v* buf, uint8_t pre_win) +{ + uint8_t f = 1; double rr; int64_t re; + + z->non_homopolymer_errors = UINT32_MAX - 1;///primary chain that needs to be verfied + + if(pre_win == 0) { + f = align_hc_ed_post_non_retrieve_extz(z, qstr, ql, tstr, tl, exz, err, wsl, OVERLAP_THRESHOLD_HIFI_FILTER, 0, NULL); + // f = align_gc_ed_post_extz(z, uref, qu->seq, tu->seq, exz, err, wsl, OVERLAP_THRESHOLD_HIFI_FILTER, 0, tot_b); + } + + if(f) { + rr = gen_extend_err_non_retrieve_exz(z, qstr, tstr, tl, exz, NULL, wsl, -1, err, (e_max+0.000001), THRESHOLD_MAX_SIZE, 0, &re); + // rr = gen_extend_err_exz(z, uref, NULL, NULL, qu->seq, tu->seq, exz, NULL, wsl, -1, err, (e_max+0.000001), THRESHOLD_MAX_SIZE, 0, &re); + z->is_match = 0; + if (rr > err) f = 0; + } + + if(f) { + z->non_homopolymer_errors = re; + // if(!gen_hc_fast_cigar(z, cl, rref, wsl, qu->seq, tu, exz, aux_o, e_rate, ql, rid, khit, &re)) + if(!gen_hc_fast_non_retrieve_cigar(z, cl, wsl, qstr, ql, tstr, tl, exz, aux_o, e_rate, rid, khit, &re)) { + // if(!gen_gc_fast_cigar0(z, cl, wsl, uref, qu->seq, tu, exz, aux_o, e_rate, ql, rid, khit, &re)) { + f = 0; + } else { + z->non_homopolymer_errors = re; + } + } + + if(f) { + if((align_gap_max >= 0) && (!ff_lunalign(z, err, align_gap_rate, align_gap_max))) f = 0; + } + + + if(f) { + if(chem_drop && ff_tend(z, 384, 2000, 0.1, (((e_rate*10)<0.36)?(e_rate*10):(0.36)), 128)) f = 0; + } + + if(f) { + reassign_gaps(z, aux_o, qstr, ql, NULL, tstr, tl, NULL, NULL, buf, exz); + } + + if((!f)) { + z->w_list.n = z->w_list.c.n = 0; + } + + return f; +} + + ///need to consider coverage, this information is missing right now (currently only use numbers) void gen_hc_r_alin_flt(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, overlap_region *aux_b, 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 sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, @@ -35113,6 +35629,66 @@ uint64_t gen_hc_r_alin_adp_mmp_0(overlap_region_alloc* ol, Candidates_list *cl, } +uint64_t gen_gc_r_alin_adp_mmp_0(overlap_region_alloc* ol, Candidates_list *cl, ul_idx_t *uref, 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 sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, + asg64_v *sp, uint64_t ocw, uint8_t *hpf, asg32_v *v32, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min, uint8_t set_match) +{ + + uint64_t i, k, bs, ql = qu->length, *wsrt = NULL, wsrt_n = 0, spn0 = 0, tot_b = 0; Window_Pool w; double err, e_max; + overlap_region *z, t; ol->mapped_overlaps_length = 0; + if(ol->length <= 0 || ql <= 0) return tot_b; + + v32->n = ol->length<<1; kv_resize(uint32_t, *v32, v32->n); + for (i = 0; i < ol->length; i++) { + v32->a[i] = ol->list[i].align_length; + v32->a[i+ol->length] = ol->list[i].shared_seed; + ol->list[i].align_length = 0; + } + + uint32_t *ocn = v32->a, *osc = v32->a + ol->length; + + // prt_chain_cluster(ol, cl, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 0, NULL); + + ///base alignment + err = e_rate; e_max = err * 1.5; + init_Window_Pool(&w, ql, wl, (int)(1.0/err)); + bs = (w.window_length)+(THRESHOLD_MAX_SIZE<<1)+1; + resize_UC_Read(tu, bs<<1); spn0 = sp->n; + + wsrt = mmp_chn_select(ol, cl, sp, ocw, ocn, osc, ql, &wsrt_n, set_match, max_n_chain, max_n_chain_f, chain_cutoff, ave_cov_min, 0.333333, 16, 16); + + for (i = 0; i < wsrt_n; i++) { + z = &(ol->list[(uint32_t)wsrt[i]]); + if((z->is_match == 0) && (!gen_gc_r_alin_flt_1_smp(z, cl, qu->seq, ql, uref->ug->u.a[z->y_id].s, uref->ug->u.a[z->y_id].len, exz, aux_o, err, e_max, e_rate, w.window_length, rid, khit, chem_drop, align_gap_rate, align_gap_max, hpf, buf, 0))) { + continue; + } + z->is_match = 1; z->strong = z->without_large_indel = 0; + } + + + for (i = k = 0; i < ol->length; i++) {///primary chain + z = &(ol->list[i]); + if(z->is_match == 0) continue; + if(k != i) { + t = ol->list[k]; + ol->list[k] = ol->list[i]; + ol->list[i] = t; + } + k++; + } + // print_mm_wins_all(wcut, wcut_n, ocw, ql, 0); + // print_mm_wins_all(wcut, wcut_n, ocw, ql, 1); + + ol->length = k; + // prt_chain_cluster(ol, cl, a_cu, a_ci, ocn, osc, idx_cu, n_cu, 1, NULL); + // fprintf(stderr, "-[M::%s]\trid::%ld\ttot_b::%lu\tql::%lu\tmax_n_chain::%ld\ttot_b_cov:::%ld\n", __func__, rid, tot_b, ql, max_n_chain, tot_b/ql); + // exit(1); + sp->n = spn0; + // if(ol->length <= 0) return tot_b; + return tot_b; + +} + + ///[ws, we) uint8_t hc_aln_simd_avx512(overlap_region* ol, uint64_t *ffa, uint32_t *ia, uint64_t in, int64_t ws, int64_t we, int64_t wl, double e_rate, All_reads *rref, char *qu, UC_Read* tu, uint64_t *fi, int32_t *baux_beg, int32_t *baux_end, int32_t *bt_s, int32_t *bt_pri_l, bit_extz_t *exz, double ovlp_cut, int64_t force_aln, uint64_t *tot_b) @@ -36598,7 +37174,7 @@ uint64_t gen_hc_r_alin_re(overlap_region* z, Candidates_list *cl, char* qstr, ui // if(z->x_id == 3196 && z->y_id == 3199) fprintf(stderr, "-0-[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); - if(!gen_hc_fast_non_retrieve_cigar(z, cl, wl, qstr, ql, tstr, tl, exz, aux_o, e_rate, rid, khit)) return 0; + if(!gen_hc_fast_non_retrieve_cigar(z, cl, wl, qstr, ql, tstr, tl, exz, aux_o, e_rate, rid, khit, &re)) return 0; // 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); diff --git a/Correct.h b/Correct.h index 59cbda9..da081fe 100644 --- a/Correct.h +++ b/Correct.h @@ -1441,6 +1441,8 @@ uint64_t gen_hc_r_alin_adp_mmp_1(overlap_region_alloc* ol, Candidates_list *cl, asg64_v *sp, uint64_t ocw, uint8_t *hpf, asg32_v *v32, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min, uint8_t set_match); uint64_t gen_hc_r_alin_adp_mmp_0(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 sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *sp, uint64_t ocw, uint8_t *hpf, asg32_v *v32, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min, uint8_t set_match); +uint64_t gen_gc_r_alin_adp_mmp_0(overlap_region_alloc* ol, Candidates_list *cl, ul_idx_t *uref, 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 sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, + asg64_v *sp, uint64_t ocw, uint8_t *hpf, asg32_v *v32, asg64_v *bp, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, uint64_t ave_cov_min, uint8_t set_match); void gen_hc_r_alin_adv_adp_smp(gen_hc_aln_t *ez, uint32_t *a_cu, uint32_t *a_ci, uint32_t *ocn, uint32_t *osc, uint64_t *idx_cu, uint64_t n_cu, uint8_t set_match); void gen_hc_r_alin_adv_adp_smp_0(gen_hc_aln_t *ez, uint8_t set_match); void gen_hc_r_alin_adv_adp_smp_1(gen_hc_aln_t *ez, uint8_t set_match); @@ -1488,6 +1490,7 @@ void est_rep_err_rate(overlap_region_alloc* ol, asg64_v *ix, kv_ul_ov_t *c_idx, #define HPC_CC_Q 3 #define HC_MF_R 0.5 #define HC_AV_MIN 0.7 +#define GC_MF_N 12 // #define FORCE_CUT 1 diff --git a/Overlaps.cpp b/Overlaps.cpp index c0a3596..3f5f2e6 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -13918,9 +13918,9 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) ug = ma_ug_gen(sg); ma_ug_seq(ug, sg, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp, 0, 1); - // if(asm_opt.dbg_bam) { - // gen_gfa_bam(ug); - // } + if(asm_opt.dbg_bam) { + gen_gfa_bam(ug, R_INF.total_reads); + } fprintf(stderr, "Writing raw unitig GFA to disk... \n"); char* gfa_name = (char*)malloc(strlen(output_file_name)+25); @@ -40260,6 +40260,7 @@ long long bubble_dist, int read_graph, int write) for (j = 0; j < reverse_sources[i].length; j++) reverse_sources[i].buffer[j].del = 0; } if(asm_opt.is_ont) cmk = gen_chemical_arc_rf(asm_opt.thread_num, R_INF.total_reads); + if(asm_opt.del_hf) clean_arc_rf(asm_opt.thread_num, R_INF.total_reads); } try_rescue_overlaps(sources, reverse_sources, n_read, 4, asm_opt.is_ont); diff --git a/Process_Read.cpp b/Process_Read.cpp index b88d50d..717cb9e 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -1843,6 +1843,62 @@ void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_ } } +void retrieve_u_seq_fast(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_t s, int64_t l, void *km) +{ + if(u->m == 0 || u->n == 0) return; + if(l < 0) l = u->len; + char *r = NULL, *a = NULL; + int64_t e = s + l, ssp, sep, rs, re, des_i; + uint64_t k, rId, ori, r_l; + if(i_r) { + i_r->length = l; i_r->RID = 0; + if(i_r->length > i_r->size) { + i_r->size = i_r->length; + if(!km) REALLOC(i_r->seq, i_r->size); + else KREALLOC(km, i_r->seq, i_r->size); + // i_r->seq = (char*)realloc(i_r->seq,sizeof(char)*(i_r->size)); + } + r = i_r->seq; + } + if(i_s) r = i_s; + + if(u->s) { + memcpy(r, u->s + s, l * sizeof((*(u->s)))); + } else { + if(strand == 1) { + sep = u->len - s; + ssp = u->len - e; + s = ssp; e = sep; + } + for (k = l = des_i = 0; k < u->n; k++) { + rId = u->a[k]>>33; + ori = u->a[k]>>32&1; + r_l = (uint32_t)u->a[k]; + if(r_l == 0) continue; + ssp = l; sep = l + r_l; + l += r_l; + if(sep <= s) continue; + if(ssp >= e) break; + rs = MAX(ssp, s); re = MIN(sep, e); + a = r + des_i; des_i += re - rs; + recover_UC_Read_sub_region(a, rs-ssp, re-rs, ori, &R_INF, rId); + } + } + + if(strand == 1) { + char t; + re = (e - s); + l = re>>1; + for (k = 0; k < (uint64_t)l; k++) { + des_i = re - k - 1; + t = r[des_i]; + r[des_i] = RC_CHAR(r[k]); + r[k] = RC_CHAR(t); + } + if(re&1) r[l] = RC_CHAR(r[l]); + } +} + uint32_t retrieve_u_cov(const ul_idx_t *ul, uint64_t id, uint8_t strand, uint64_t pos, uint8_t dir, int64_t *pi) { uint64_t *a = ul->cc->interval.a + ul->cc->idx[id], cc = 0, ff = 0; diff --git a/Process_Read.h b/Process_Read.h index b8486f5..b2e10d8 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -268,6 +268,7 @@ int64_t load_compress_base_disk(FILE *fp, uint64_t *ul_rid, char *dest, uint32_t scaf_res_t *init_scaf_res_t(uint32_t n); void destroy_scaf_res_t(scaf_res_t *p); void read_ma(ma_hit_t* x, FILE* fp); +void retrieve_u_seq_fast(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_t s, int64_t l, void *km); const uint64_t sc_tb[8] = { 10, 20, 30, 40, 50, 60, 70, 80 diff --git a/anchor.cpp b/anchor.cpp index ab7f4f4..311fbff 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -2674,6 +2674,31 @@ void ul_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t 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, mcopy_rate, chain_cutoff, mcopy_khit_cut, sp); } +void p_ec_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, ul_idx_t *uref, overlap_region_alloc *ol, Candidates_list *cl, double bw_thres, + int apend_be, kvec_t_u8_warp* k_flag, 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, int64_t mcopy_num, double mcopy_rate, uint32_t mcopy_khit_cut) +{ + extern void *ha_flt_tab; + extern ha_pt_t *ha_idx; + 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); + + minimizers_qgen(ab, rs, rl, mz_w, mz_k, cl, k_flag, ha_flt_tab, ha_idx, NULL, uref, dbg_ct, sp, high_occ, low_occ); + + uint64_t k, l, m, cn = cl->length, yid; + clear_overlap_region_alloc(ol); + for (l = 0, k = 1, m = 0; k <= cn; k++) { + if((k == cn) || (cl->list[k].readID != cl->list[l].readID)) { + if(cl->list[l].readID != rid) { + yid = cl->list[l].readID; + m += lchain_qdp_mcopy_fast(cl, l, k-l, m, &(cl->chainDP), ol, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, + rid, rl, uref->ug->u.a[yid].len, quick_check, apend_be, gen_off, mcopy_num, mcopy_rate, mcopy_khit_cut, 1); + } + l = k; + } + } + cl->length = m; +} + void h_ec_lchain(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, All_reads *rref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, 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, int64_t mcopy_num, double mcopy_rate, uint32_t chain_cutoff, uint32_t mcopy_khit_cut, uint64_t ocv_w, uint8_t is_raw_chain) { diff --git a/ecovlp.cpp b/ecovlp.cpp index 579988e..686bec2 100644 --- a/ecovlp.cpp +++ b/ecovlp.cpp @@ -21,6 +21,11 @@ #define HC0_W 6 #define DBG_TIME 0 + + +void p_ec_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, ul_idx_t *uref, overlap_region_alloc *ol, Candidates_list *cl, double bw_thres, + int apend_be, kvec_t_u8_warp* k_flag, 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, int64_t mcopy_num, double mcopy_rate, uint32_t mcopy_khit_cut); + typedef struct { uint64_t fbs, faln; double chn_tm, aln_tm, phs_tm, cns_tm; @@ -59,35 +64,123 @@ typedef struct { kdq_t(uint32_t) *q; }cns_gfa; -typedef struct { - // chaining and overlapping related buffers - UC_Read self_read, ovlp_read; - Candidates_list clist; - overlap_region_alloc olist; - ha_abuf_t *ab; - // int64_t num_read_base, num_correct_base, num_recorrect_base; - uint64_t cnt[6], rr; - haplotype_evdience_alloc hap; - bit_extz_t exz; - kv_ul_ov_t pidx; - asg64_v v64; - asg32_v v32; - asg16_v v16; - asg8_v v8q, v8t; +void destroy_cns_gfa(cns_gfa *p) +{ + size_t k; + for (k = 0; k < p->m; k++) { + kv_destroy(p->a[k].arc); + } + free(p->a); kdq_destroy(uint32_t, p->q); +} - kvec_t_u8_warp k_flag; - st_mt_t sp; - cns_gfa cns; -} ec_ovec_buf_t0; +#define CAT2(a, b) a##b +#define XCAT2(a, b) CAT2(a, b) +#define CAT3(a, b, c) XCAT2(XCAT2(a, b), c) + +#define ECZ_INIT(sf, HType) \ +typedef struct {\ + UC_Read self_read, ovlp_read;\ + Candidates_list clist;\ + overlap_region_alloc olist;\ + HType##_t *ab;\ + uint64_t cnt[6], rr;\ + haplotype_evdience_alloc hap;\ + bit_extz_t exz;\ + kv_ul_ov_t pidx;\ + asg64_v v64;\ + asg32_v v32;\ + asg16_v v16;\ + asg8_v v8q, v8t;\ + kvec_t_u8_warp k_flag;\ + st_mt_t sp;\ + cns_gfa cns;\ +} sf##_ovec_buf_t0;\ +typedef struct {\ + sf##_ovec_buf_t0 *a;\ + uint32_t n, rev;\ + uint8_t *cr;\ +} sf##_ovec_buf_t;\ +sf##_ovec_buf_t* CAT3(gen_, sf, _ovec_buf_t)(uint32_t n) {\ + uint32_t k; sf##_ovec_buf_t0 *z = NULL;\ + sf##_ovec_buf_t *p = NULL; CALLOC(p, 1);\ + p->n = n; CALLOC(p->a, p->n);\ + for (k = 0; k < p->n; k++) {\ + z = &(p->a[k]);\ + init_UC_Read(&z->self_read);\ + init_UC_Read(&z->ovlp_read);\ + init_Candidates_list(&z->clist);\ + init_overlap_region_alloc(&z->olist);\ + kv_init(z->k_flag.a);\ + kv_init(z->sp);\ + kv_init(z->pidx);\ + kv_init(z->v64);\ + kv_init(z->v32);\ + kv_init(z->v16);\ + kv_init(z->v8q);\ + kv_init(z->v8t);\ + init_bit_extz_t(&(z->exz), 31);\ + z->ab = HType##_init();\ + InitHaplotypeEvdience(&z->hap);\ + z->cns.q = kdq_init(uint32_t);\ + }\ + return p;\ +}\ +void CAT3(destroy_, sf, _ovec_buf_t)(sf##_ovec_buf_t *p)\ +{\ + uint32_t k; sf##_ovec_buf_t0 *z = NULL;\ + for (k = 0; k < p->n; k++) {\ + z = &(p->a[k]); z->rr = 0;\ + destory_UC_Read(&z->self_read);\ + destory_UC_Read(&z->ovlp_read);\ + destory_Candidates_list(&z->clist);\ + destory_overlap_region_alloc(&z->olist);\ + kv_destroy(z->k_flag.a);\ + kv_destroy(z->sp);\ + kv_destroy(z->pidx);\ + kv_destroy(z->v64);\ + kv_destroy(z->v32);\ + kv_destroy(z->v16);\ + kv_destroy(z->v8q);\ + kv_destroy(z->v8t);\ + destroy_bit_extz_t(&(z->exz));\ + HType##_destroy(z->ab);\ + destoryHaplotypeEvdience(&z->hap);\ + destroy_cns_gfa(&(z->cns));\ + }\ + free(p->a); free(p->cr); free(p);\ +} + +ECZ_INIT(ec, ha_abuf) +ECZ_INIT(gc, ha_abufl) typedef struct { - ec_ovec_buf_t0 *a; - uint32_t n, rev; - uint8_t *cr; -} ec_ovec_buf_t; + uint64_t tns; uint32_t te; + uint32_t qn, qs, qe; +} raln_t; + +typedef struct { + size_t n, m; + raln_t *a; +} araln_t_v; + +typedef struct { + araln_t_v *res; gc_ovec_buf_t *p; + uint32_t rsi, rei, n_thre; + ul_idx_t *uu; uint64_t sum_len; +} rept_v; + +typedef struct { + uint32_t n_thre, rn, tot_r; + uint64_t sum_len, chunk_len; + ul_idx_t *uu; +} ec_rept_t; typedef struct { ec_ovec_buf_t *p; + uint32_t n_thre, rn; + uint64_t sum_len, chunk_len; + // araln_t_v *al; + // uint32_t al_n; asg64_v idx; ma_ug_t *ug; } ec_polish_buf_t; @@ -160,88 +253,6 @@ void h_ec_lchain_fast_new(ha_abuf_t *ab, uint32_t rid, UC_Read *qu, UC_Read *tu, void gen_chain_clus(ha_abuf_t *ab, overlap_region_alloc *ol, Candidates_list *cl, asg32_v *ik); void srt_olst(overlap_region_alloc* ol); -ec_ovec_buf_t* gen_ec_ovec_buf_t(uint32_t n) -{ - uint32_t k; ec_ovec_buf_t0 *z = NULL; - ec_ovec_buf_t *p = NULL; CALLOC(p, 1); - p->n = n; CALLOC(p->a, p->n); - for (k = 0; k < p->n; k++) { - z = &(p->a[k]); - init_UC_Read(&z->self_read); - init_UC_Read(&z->ovlp_read); - init_Candidates_list(&z->clist); - init_overlap_region_alloc(&z->olist); - - // init_fake_cigar(&(z->tmp.f_cigar)); - // memset(&(z->tmp.w_list), 0, sizeof(z->tmp.w_list)); - // CALLOC(z->tmp.w_list.a, 1); z->tmp.w_list.n = z->tmp.w_list.m = 1; - - // kv_init(z->b_buf.a); - // kv_init(z->r_buf.a); - kv_init(z->k_flag.a); - kv_init(z->sp); - kv_init(z->pidx); - kv_init(z->v64); - kv_init(z->v32); - kv_init(z->v16); - kv_init(z->v8q); - kv_init(z->v8t); - init_bit_extz_t(&(z->exz), 31); - - z->ab = ha_abuf_init(); - - InitHaplotypeEvdience(&z->hap); - z->cns.q = kdq_init(uint32_t); - } - - return p; -} - -void destroy_cns_gfa(cns_gfa *p) -{ - size_t k; - for (k = 0; k < p->m; k++) { - kv_destroy(p->a[k].arc); - } - free(p->a); kdq_destroy(uint32_t, p->q); -} - -void destroy_ec_ovec_buf_t(ec_ovec_buf_t *p) -{ - uint32_t k; ec_ovec_buf_t0 *z = NULL; - for (k = 0; k < p->n; k++) { - z = &(p->a[k]); z->rr = 0; - destory_UC_Read(&z->self_read); - destory_UC_Read(&z->ovlp_read); - destory_Candidates_list(&z->clist); - destory_overlap_region_alloc(&z->olist); - - // destory_fake_cigar(&(z->tmp.f_cigar)); - // free(z->tmp.w_list.a); free(z->tmp.w_list.c.a); - - // kv_destroy(z->r_buf.a); - kv_destroy(z->k_flag.a); - kv_destroy(z->sp); - kv_destroy(z->pidx); - kv_destroy(z->v64); - kv_destroy(z->v32); - kv_destroy(z->v16); - kv_destroy(z->v8q); - kv_destroy(z->v8t); - destroy_bit_extz_t(&(z->exz)); - - ha_abuf_destroy(z->ab); - - destoryHaplotypeEvdience(&z->hap); - destroy_cns_gfa(&(z->cns)); - - } - free(p->a); free(p->cr); free(p); - - // fprintf(stderr, "[M::%s-chains] #->%lld\n", __func__, asm_opt.num_bases); - // fprintf(stderr, "[M::%s-passed-chains-0] #->%lld\n", __func__, asm_opt.num_corrected_bases); - // fprintf(stderr, "[M::%s-cis-chains-1] #->%lld\n", __func__, asm_opt.num_recorrected_bases); -} inline void refresh_ec_ovec_buf_t0(ec_ovec_buf_t0 *z, uint64_t n) { @@ -274,6 +285,37 @@ inline void refresh_ec_ovec_buf_t0(ec_ovec_buf_t0 *z, uint64_t n) } } +inline void refresh_gc_ovec_buf_t0(gc_ovec_buf_t0 *z, uint64_t n) +{ + z->rr++; + if((z->rr%n) == 0) { + free(z->self_read.seq); memset(&(z->self_read), 0, sizeof(z->self_read)); + free(z->ovlp_read.seq); memset(&(z->ovlp_read), 0, sizeof(z->ovlp_read)); + + destory_Candidates_list(&z->clist); memset(&(z->clist), 0, sizeof(z->clist)); + destory_overlap_region_alloc(&z->olist); memset(&(z->olist), 0, sizeof(z->olist)); init_overlap_region_alloc(&z->olist); + + kv_destroy(z->k_flag.a); kv_init(z->k_flag.a); + kv_destroy(z->sp); kv_init(z->sp); + kv_destroy(z->pidx); kv_init(z->pidx); + kv_destroy(z->v64); kv_init(z->v64); + kv_destroy(z->v32); kv_init(z->v32); + kv_destroy(z->v16); kv_init(z->v16); + kv_destroy(z->v8q); kv_init(z->v8q); + kv_destroy(z->v8t); kv_init(z->v8t); + + destroy_bit_extz_t(&(z->exz)); init_bit_extz_t(&(z->exz), 31); + + ha_abufl_destroy(z->ab); z->ab = ha_abufl_init(); + + destoryHaplotypeEvdience(&z->hap); memset(&(z->hap), 0, sizeof(z->hap)); InitHaplotypeEvdience(&z->hap); + + destroy_cns_gfa(&(z->cns)); memset(&(z->cns), 0, sizeof(z->cns)); z->cns.q = kdq_init(uint32_t); + + // z->rr = 1; + } +} + void prt_chain(overlap_region_alloc *o) { @@ -3407,7 +3449,7 @@ uint64_t gen_hc_r_alin_ea_flt(ha_abuf_t *ab, overlap_region_alloc* ol, Candidate } -uint64_t gen_hc_r_alin_ea_flt_mmp(ha_abuf_t *ab, overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, 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, +uint64_t gen_hc_r_alin_ea_flt_mmp(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, uint64_t max_n_chain, uint64_t max_n_chain_f, uint64_t chain_cutoff, 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, uint64_t sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *bp, uint64_t ocw, asg8_v *hpz, asg32_v *v32) { if(ol->length <= 0) return 0; @@ -4390,6 +4432,8 @@ static void worker_hap_ec(void *data, long i, int tid) tt0 = tt1; } + + // stderr_phase_ovlp(&b->olist); ///debug for memory @@ -4407,7 +4451,7 @@ static void worker_hap_ec(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); // tot_b = gen_hc_r_alin_ea_flt(b->ab, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_n_chain, asm_opt.max_n_chain*HC_MF_R, asm_opt.chn_occ, 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]), asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1), (asm_opt.is_ont)?(512):(0), (asm_opt.is_ont)?(6):(0), (asm_opt.is_ont)?(1.5):(-1), (asm_opt.is_ont)?(0.1):(-1), &buf0, qw, &b->v8q, &b->v32, 1); - tot_b = gen_hc_r_alin_ea_flt_mmp(b->ab, &b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_n_chain, asm_opt.max_n_chain*HC_MF_R, asm_opt.chn_occ, asm_opt.max_ov_diff_ec, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), i, E_KHIT, + tot_b = gen_hc_r_alin_ea_flt_mmp(&b->olist, &b->clist, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, aux_o, asm_opt.max_n_chain, asm_opt.max_n_chain*HC_MF_R, asm_opt.chn_occ, 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]), asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1), (asm_opt.is_ont)?(512):(0), (asm_opt.is_ont)?(6):(0), (asm_opt.is_ont)?(1.5):(-1), (asm_opt.is_ont)?(0.1):(-1), &buf0, qw, &b->v8q, &b->v32); copy_asg_arr(b->sp, buf0); @@ -4466,6 +4510,8 @@ static void worker_hap_ec(void *data, long i, int tid) ///for debug indel // stderr_phase_ovlp(&b->olist); + return; + // est_rep_err_rate(&b->olist, &b->v64, &b->pidx, qlen, (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), NULL); if(DBG_TIME && dbg_a) { @@ -4583,10 +4629,75 @@ static void worker_hap_ec(void *data, long i, int tid) //fprintf(stderr, "-[M::%s]\trid::%ld\t%.*s\n", __func__, i, (int)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i)); } -// void dbg_gen_ovlp_lst() -// { -// ; -// } +void gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t rid); + +static void worker_gfa_ec(void *data, long i, int tid) +{ + rept_v *z = ((rept_v*)data); + gc_ovec_buf_t0 *b = &(z->p->a[tid]); i += z->rsi; + uint32_t high_occ = 2; uint64_t qid = i; int64_t het_a, hom_a, wl, max_chn, max_chn_f; + overlap_region *aux_o = NULL; double bw_h, e_h; + asg64_v buf0; uint32_t qlen = 0, qw = 0; + b->v8q.n = b->v8t.n = 0; set_ec_cov(asm_opt.het_cov, asm_opt.hom_cov, asm_opt.het_cov_set, asm_opt.polyploidy, het_a, hom_a); + + max_chn = asm_opt.max_n_chain/het_a; if(max_chn < GC_MF_N) max_chn = GC_MF_N; max_chn_f = max_chn*HC_MF_R; if(max_chn_f < 1) max_chn_f = 1; + e_h = asm_opt.max_ov_diff_ec; + if((asm_opt.is_ont) && (qid < R_INF.tqn)) { + bw_h = 0.05; wl = WINDOW_OHC; + } else { + bw_h = 0.02; wl = WINDOW_HC; if(asm_opt.is_ont) e_h = asm_opt.max_ov_diff_ec_sec;///HiFi reads in hybrid model + } + + recover_UC_Read(&b->ovlp_read, &R_INF, i); + gen_ori_seq0(b->ovlp_read.seq, b->ovlp_read.length, &b->self_read, &(scb.a[i]), i); + qlen = b->self_read.length; + + qw = ((qlen < (COV_W_AC<<1))?(qlen>>1):(COV_W_AC)); if(!qw) qw = 1; + + p_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, z->uu, &b->olist, &b->clist, bw_h, 1, NULL, NULL, &(b->sp), &high_occ, NULL, 1, 1, 3, 0.7, 32);///ONT high error + b->cnt[0] += b->self_read.length; + + aux_o = fetch_aux_ovlp(&b->olist, NULL);///must be here + + copy_asg_arr(buf0, b->sp); + gen_gc_r_alin_adp_mmp_0(&b->olist, &b->clist, z->uu, &b->self_read, &b->ovlp_read, &b->exz, aux_o, e_h, wl, qid, E_KHIT, + 1, &b->v16, asm_opt.is_ont, (asm_opt.is_ont)?(0.006):(-1), (asm_opt.is_ont)?(64):(-1), (asm_opt.is_ont)?(512):(0), (asm_opt.is_ont)?(6):(0), (asm_opt.is_ont)?(1.5):(-1), (asm_opt.is_ont)?(0.1):(-1), + &b->v64, qw, NULL, &b->v32, &buf0, max_chn, max_chn_f, asm_opt.chn_occ, ((hom_a*HC_AV_MIN)>0)?(hom_a*HC_AV_MIN):(1), 1); + copy_asg_arr(b->sp, buf0); + + /** + copy_asg_arr(buf0, b->sp); + rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, ((uint64_t)-1), 0, HC0_W, &b->v32, + asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, -1.0); + copy_asg_arr(b->sp, buf0); + + dedup_chains(&b->olist); + + + + copy_asg_arr(buf0, b->sp); + b->cnt[1] += wcns_gen(&b->olist, &R_INF, i, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, 512, b->self_read.length, 3, 0.500001, aux_o, &b->v32, &b->cns, 256, i, ((uint64_t)-1), + R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, NULL); + copy_asg_arr(b->sp, buf0); + + push_nec_re(aux_o, &(scc.a[i])); + // push_nec_re(aux_o, &(scb.a[i])); + if(asm_opt.dbg_bam) { + update_scb(&R_INF, &(scc.a[i]), &(scb.a[i]), &(b->v16), &b->self_read, &b->ovlp_read, &b->v64, &b->exz, i); + kv_resize(uint16_t, scb.a[i], b->v16.n); scb.a[i].n = b->v16.n; memcpy(scb.a[i].a, b->v16.a, b->v16.n*sizeof(*(scb.a[i].a))); + } + + push_ne_ovlp(&(R_INF.paf[i]), &b->olist, 1, &R_INF, &(scc.a[i])); + push_ne_ovlp(&(R_INF.reverse_paf[i]), &b->olist, 2, &R_INF, NULL); + + + check_well_cal(&(scc.a[i]), &b->v64, &(R_INF.paf[i].is_fully_corrected), &(R_INF.paf[i].is_abnormal), qlen, (MIN_COVERAGE_THRESHOLD*2), &(R_INF.paf[i])); + R_INF.trio_flag[i] = AMBIGU; + **/ + + refresh_gc_ovec_buf_t0(b, REFRESH_N); +} + void worker_hap_ec_back_dbg(void *data, long i, int tid) { @@ -6934,6 +7045,23 @@ static void worker_hap_dc_ec_chemical_arc_mark(void *data, long i, int tid) refresh_ec_ovec_buf_t0(b, REFRESH_N); } +static void worker_clean_hf(void *data, long i, int tid) +{ + uint64_t k; ma_hit_t_alloc *p = &(R_INF.paf[i]); + for (k = 0; k < p->length; k++) { + if((Get_qn(p->buffer[k]) >= R_INF.tqn) || (Get_tn(p->buffer[k]) >= R_INF.tqn)) { + p->buffer[k].del = 1; + } + } + + p = &(R_INF.reverse_paf[i]); + for (k = 0; k < p->length; k++) { + if((Get_qn(p->buffer[k]) >= R_INF.tqn) || (Get_tn(p->buffer[k]) >= R_INF.tqn)) { + p->buffer[k].del = 1; + } + } +} + uint64_t get_candidate_rrs(ma_utg_t *u, uint64_t rz) { // uint64_t rid, rs, re, rev, k, l[2], lr, ts, te; @@ -9118,6 +9246,7 @@ uint64_t cal_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, uint64 // fprintf(stderr, "[M::%s] # corrected bases->%lu\n", __func__, num_correct); // fprintf(stderr, "[M::%s::%.3f] running time\n", __func__, yak_realtime_0()-tt0); fprintf(stderr, "[M::pec::%.3f] # bases: %lu; # corrected bases: %lu\n", yak_realtime_0()-tt0, num_base, num_correct); + exit(1); (*r_base) = num_base; return num_correct; @@ -9460,8 +9589,74 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u // } } -void gen_gfa_bam(ma_ug_t *ug) +void ugfa_build(ma_ug_t *ug, int64_t ak, int64_t aw, int64_t max_cutoff) { + ha_flt_tab = ha_ft_ul_gen(&asm_opt, &(ug->u), ak, aw, -1, max_cutoff); + ha_idx = ha_pt_ul_gen(&asm_opt, ha_flt_tab, &(ug->u), ak, aw, max_cutoff); + fprintf(stderr, "[M::%s] Index has been built.\n", __func__); +} + +ec_polish_buf_t* init_ec_polish_buf_t(ec_polish_buf_t *z, int64_t n_thre, uint64_t chunk_size) +{ + if(z == NULL) { + //araln_t_v + ec_polish_buf_t *p; CALLOC(p, 1); + p->p = gen_ec_ovec_buf_t(n_thre); + p->n_thre = n_thre; p->chunk_len = chunk_size; + return p; + } else { + destroy_ec_ovec_buf_t(z->p); + free(z->idx.a); free(z); + return NULL; + } +} + +static void *worker_rec_pipeline(void *data, int step, void *in) // callback for kt_pipeline() +{ + ec_rept_t *p = (ec_rept_t *)data; + if (step == 0) { + rept_v *s = NULL; CALLOC(s, 1); uint32_t k; + s->n_thre = p->n_thre; s->uu = p->uu; + s->rsi = s->rei = p->rn; + CALLOC(s->res, s->n_thre); s->p = gen_gc_ovec_buf_t(s->n_thre); + + for (k = p->rn, s->sum_len = 0; k < p->tot_r && s->sum_len < p->chunk_len; k++) { + s->sum_len += Get_READ_LENGTH(R_INF, k); + } + s->rei = p->rn = k; + + kt_for(p->n_thre, worker_gfa_ec, s, s->rei - s->rsi); + + if (s->sum_len == 0) free(s); + else return s; + } else if (step == 1) { + rept_v *s = (rept_v*)in; uint32_t k; + destroy_gc_ovec_buf_t(s->p); + + + + + for (k = 0; k < s->n_thre; k++) free(s->res[k].a); + free(s->res); free(s); + } + return 0; +} + +void ugfa_aln(ma_ug_t *ug, uint64_t n_a) +{ + ec_rept_t *sl = NULL; CALLOC(sl, 1); sl->tot_r = n_a; + sl->chunk_len = asm_opt.thread_num*64000000; sl->n_thre = asm_opt.thread_num; sl->uu->ug = ug; + + kt_pipeline(3, worker_rec_pipeline, sl, 2); + + free(sl); +} + +///R_INF.total_reads +void gen_gfa_bam(ma_ug_t *ug, uint64_t n_a) +{ + // ec_ovec_buf_t *b = NULL; + if(ha_flt_tab) { ha_ft_destroy(ha_flt_tab); ha_flt_tab = NULL; } @@ -9469,8 +9664,10 @@ void gen_gfa_bam(ma_ug_t *ug) ha_pt_destroy(ha_idx); ha_idx = NULL; } - ha_flt_tab = ha_ft_ug_gen(&asm_opt, &(ug->u), 0, asm_opt.k_mer_length, asm_opt.mz_win, -1, -1); - // ha_idx = ha_pt_ug_gen(&asm_opt, ha_flt_tab, &(ug->u), hap_n); + ugfa_build(ug, asm_opt.k_mer_length, asm_opt.mz_win, asm_opt.max_n_chain); + ugfa_aln(ug, n_a); + fprintf(stderr, "[M::%s] gfa bam is done.\n", __func__); + exit(1); } void print_ov_dbg_paf(FILE *fp, char *ref_str, char *ref_id, int32_t ref_id_n, char *qry_str, char *qry_id, int32_t qry_id_n, uint64_t rs, uint64_t re, uint64_t rl, uint64_t qs, uint64_t qe, uint64_t ql, uint64_t rev, bit_extz_t *ez, char *ezh) @@ -9670,6 +9867,11 @@ uint8_t* gen_chemical_arc_rf(uint64_t n_thre, uint64_t n_a) return ra; } +void clean_arc_rf(uint64_t n_thre, uint64_t n_a) +{ + kt_for(n_thre, worker_clean_hf, NULL, n_a); +} + void gen_hc_polish(uint64_t n_thre, ma_ug_t *ug) { diff --git a/ecovlp.h b/ecovlp.h index 7ddf18a..60cf852 100644 --- a/ecovlp.h +++ b/ecovlp.h @@ -18,6 +18,7 @@ uint8_t* gen_chemical_arc_rf(uint64_t n_thre, uint64_t n_a); void cal_ec_r_dbg(uint64_t n_thre, uint64_t n_a); void write_ec_reads(const char *suffix_ou, cc_v *cvt, uint8_t is_rev); void destroy_cc_v(cc_v *z); -void gen_gfa_bam(ma_ug_t *ug); +void gen_gfa_bam(ma_ug_t *ug, uint64_t n_a); +void clean_arc_rf(uint64_t n_thre, uint64_t n_a); #endif \ No newline at end of file diff --git a/htab.cpp b/htab.cpp index e44fbf6..91da897 100644 --- a/htab.cpp +++ b/htab.cpp @@ -1120,13 +1120,26 @@ void debug_ct_index(void* q_ct_idx, void* r_ct_idx) /************************* * High-level interfaces * *************************/ -void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int cutoff) + + int64_t ha_ct_ug_cutoff(ha_ct_t *h, int64_t num_thre, double cut_rate) { - yak_ft_t *flt_tab; - ha_ct_t *h; + int64_t cnt[YAK_N_COUNTS], k, tot_n = 0, tot_cutn = 0; + ha_ct_hist(h, cnt, num_thre); + for (k = tot_n = 0; k < YAK_N_COUNTS; k++) tot_n += cnt[k]; + tot_cutn = tot_n - (tot_n*cut_rate); + + for (k = tot_n = 0; k < YAK_N_COUNTS && tot_n < tot_cutn; k++) tot_n += cnt[k]; + return k; +} + +void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int cutoff, int max_cutoff) +{ + yak_ft_t *flt_tab; ha_ct_t *h; ///HAF_COUNT_EXACT ---> no bf; HAF_COUNT_ALL ---> no minimizer h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_UG_READ, !(asm_opt->flag&HA_F_NO_HPC), k, w, NULL, NULL, NULL, us, 0, NULL, 0); - + + if(cutoff < 0) cutoff = ha_ct_ug_cutoff(h, asm_opt->thread_num, 0.0002); + if(cutoff > max_cutoff) cutoff = max_cutoff; // cutoff = (int)(asm_opt->hom_cov * asm_opt->high_factor); if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1; // fprintf(stderr, "[M::%s::] cutoff->%d\n\n", __func__, cutoff); @@ -1138,21 +1151,6 @@ void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int return (void*)flt_tab; } -// int64_t ha_ct_ug_cutoff(ha_ct_t *h, int64_t num_thre, double cut_rate, double cut_an_rate) -// { -// int64_t cnt[YAK_N_COUNTS], k, tot_n = 0, tot_cutn = 0, tot_an = 0, tot_cut_an = 0, tot_k = 0, tot_ak = 0; -// ha_ct_hist(h, cnt, num_thre); -// for (k = tot_n = tot_an = 0; k < YAK_N_COUNTS; k++) { -// tot_n += cnt[k]; tot_an += cnt[k]*k; -// } - -// tot_cutn = tot_n - (tot_n*cut_rate); -// tot_cut_an = tot_an - (tot_an*cut_an_rate); - - - -// } - void *ha_ft_ug_gen(hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq) { yak_ft_t *flt_tab; diff --git a/htab.h b/htab.h index d7ca3d6..0ae8ada 100644 --- a/htab.h +++ b/htab.h @@ -72,7 +72,7 @@ extern void *ha_flt_tab_hp; extern ha_pt_t *ha_idx_hp; extern void *ha_ct_table; -void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int cutoff); +void *ha_ft_ul_gen(const hifiasm_opt_t *asm_opt, ma_utg_v *us, int k, int w, int cutoff, int max_cutoff); void *ha_ft_ug_gen(hifiasm_opt_t *asm_opt, ma_utg_v *us, int is_HPC, int k, int w, int min_freq, int max_freq); void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs, int *hom_cov, int is_hp_mode, int read_from_store); int32_t ha_ft_cnt(const void *hh, uint64_t y); diff --git a/inter.cpp b/inter.cpp index e70d29b..b9e16a9 100644 --- a/inter.cpp +++ b/inter.cpp @@ -514,7 +514,7 @@ double diff_ec_ul, double diff_ec_ul_low, double diff_ec_ul_hpc, int ec_ul_round void uidx_l_build(ma_ug_t *ug, mg_idxopt_t *opt, int cutoff) { - ha_flt_tab = ha_ft_ul_gen(&asm_opt, &(ug->u), opt->k, opt->w, cutoff); + ha_flt_tab = ha_ft_ul_gen(&asm_opt, &(ug->u), opt->k, opt->w, cutoff, -1); ha_idx = ha_pt_ul_gen(&asm_opt, ha_flt_tab, &(ug->u), opt->k, opt->w, cutoff); fprintf(stderr, "[M::%s] Index has been built.\n", __func__); }