From 7884b5ad88bf2053fd9819be02e47b43e2f638c9 Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Sun, 3 May 2026 03:34:23 -0400 Subject: [PATCH] regen_scb --- Assembly.cpp | 17 +- CommandLines.cpp | 5 + CommandLines.h | 4 +- Correct.cpp | 1044 ++++++++++++++++++++++++++++++++++++++-- Correct.h | 5 +- Hash_Table.cpp | 481 ++++++++++++++++++ Hash_Table.h | 10 + Levenshtein_distance.h | 33 +- Overlaps.cpp | 4 +- Process_Read.cpp | 31 +- anchor.cpp | 150 ++++++ ecovlp.cpp | 495 +++++++++++++++++-- gfa_ut.cpp | 3 +- htab.cpp | 26 +- htab.h | 2 +- 15 files changed, 2177 insertions(+), 133 deletions(-) diff --git a/Assembly.cpp b/Assembly.cpp index 02165dd..2821086 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1026,7 +1026,7 @@ void ha_ec(int64_t round, int num_pround, int des_idx, uint64_t *tot_b, uint64_t write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name); if((asm_opt.flag & HA_F_VERBOSE_GFA) && (asm_opt.bin_only == 1)) exit(1);///just for debug } - if (w_tmp) tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, round, asm_opt.number_of_round, 0); + if (w_tmp) tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, round, asm_opt.number_of_round, 0); // Output_corrected_fastq(); @@ -2077,7 +2077,7 @@ int ha_assemble(void) // quick_debug_phasing(MC_NAME); extern void ha_extract_print_list(const All_reads *rs, int n_rounds, const char *o); int r, r0 = -1, hom_cov = -1, ovlp_loaded = 0; uint64_t tot_b, tot_e; - if (asm_opt.load_index_from_disk && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) { + if ((asm_opt.load_index_from_disk) && (asm_opt.dbg_ec_rr < 0) && load_all_data_from_disk(&R_INF.paf, &R_INF.reverse_paf, asm_opt.output_file_name)) { ovlp_loaded = 1; fprintf(stderr, "[M::%s::%.3f*%.2f] ==> loaded corrected reads and overlaps from disk\n", __func__, yak_realtime(), yak_cpu_usage()); if (asm_opt.extract_list) { @@ -2096,12 +2096,18 @@ int ha_assemble(void) if((asm_opt.flag & HA_F_VERBOSE_GFA)) load_pt_index(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name), load_ct_index(&ha_ct_table, asm_opt.output_file_name); r = ha_idx?asm_opt.number_of_round-1:0; if((!ha_idx) && (asm_opt.restart)) { - for (r = asm_opt.number_of_round - 1; r >= 0; --r) { - if(tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name, r, asm_opt.number_of_round, 1)) { + r = asm_opt.number_of_round - 1; + if(asm_opt.dbg_ec_rr >= 0) r = asm_opt.dbg_ec_rr; + for (; r >= 0; --r) { + if(tmp_pt_pro(&ha_flt_tab, &ha_idx, &R_INF, &scb, &asm_opt, asm_opt.output_file_name, r, asm_opt.number_of_round, 1)) { load_ct_index(&ha_ct_table, asm_opt.output_file_name); r0 = r; break; } } + if((asm_opt.dbg_ec_rr >= 0) && (asm_opt.dbg_ec_rr != r)) { + fprintf(stderr, "[E::%s] no matching debug error-correction bins found\n", __func__); + exit(1); + } if(r < 0) r = 0; } @@ -2123,6 +2129,7 @@ int ha_assemble(void) // fprintf(stderr, "[M::%s] # bases: %lld; # corrected bases: %lld; # recorrected bases: %lld\n", __func__, // asm_opt.num_bases, asm_opt.num_corrected_bases, asm_opt.num_recorrected_bases); // fprintf(stderr, "[M::%s] size of buffer: %.3fGB\n", __func__, asm_opt.mem_buf / 1073741824.0); + if(asm_opt.dbg_ec_rr >= 0) exit(1); } if (asm_opt.flag & HA_F_WRITE_EC) { if(asm_opt.is_sc) Output_corrected_fastq(); @@ -2149,7 +2156,7 @@ int ha_assemble(void) build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf, R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round, asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, !ovlp_loaded); - destory_All_reads(&R_INF); if(asm_opt.dbg_bam) destroy_cc_v(&scb); + destory_All_reads(&R_INF); /**if(asm_opt.dbg_bam)**/ destroy_cc_v(&scb); return 0; } diff --git a/CommandLines.cpp b/CommandLines.cpp index 40cf1a9..debc942 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -94,6 +94,7 @@ static ko_longopt_t long_options[] = { { "hyb-syn", ko_required_argument, 376}, { "simd-m", ko_required_argument, 377}, { "del-hf", ko_no_argument, 378}, + { "dbg-rr", ko_required_argument, 379}, // { "path-round", ko_required_argument, 348}, { 0, 0, 0 } }; @@ -444,6 +445,8 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->simd_mm = -1; asm_opt->del_hf = 0; + + asm_opt->dbg_ec_rr = -1; } void destory_enzyme(enzyme* f) @@ -1121,6 +1124,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) asm_opt->simd_mm = atoi(opt.arg); } else if (c == 378) { asm_opt->del_hf = 1; + } else if (c == 379) { + asm_opt->dbg_ec_rr = atoi(opt.arg); } 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 4c4d24a..639a6ce 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -5,7 +5,7 @@ #include #include -#define HA_VERSION "0.25.1-r910" +#define HA_VERSION "0.25.1-r920" #define VERBOSE 0 @@ -204,6 +204,8 @@ typedef struct { int8_t simd_mm; int8_t del_hf; + + int64_t dbg_ec_rr; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Correct.cpp b/Correct.cpp index 7aff11b..2b39ad7 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -8847,7 +8847,7 @@ void prt_sub_cigar(overlap_region* z, uint64_t str_l, uint64_t site, uint64_t wi #define is_st_bs(s, rr, mm) (((mm) != ((uint64_t)-1)) && (((s).overlap_num + mm) >= ((s).occ_0)) && ((((s).occ_0*(rr) + (s).overlap_num)) >= ((s).occ_0))) -void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, UC_Read* g_read, uint64_t multi_check, double st_rate, uint64_t st_max, +void generate_haplotypes_naive_HiFi(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, uint64_t multi_check, double st_rate, uint64_t st_max, int64_t hap_cov_match, int64_t hap_cov_unmatch) { // fprintf(stderr, "[M::%s::] Done\n", __func__); @@ -9927,7 +9927,7 @@ void generate_haplotypes_naive_HiFi_adv(haplotype_evdience_alloc* hap, overlap_r -void generate_haplotypes_naive_HiFi_adv_hc(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, UC_Read* g_read, double st_rate, uint64_t st_max, asg32_v *b32, uint64_t rid, int64_t hap_cov_match, int64_t hap_cov_unmatch) +void generate_haplotypes_naive_HiFi_adv_hc(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, double st_rate, uint64_t st_max, asg32_v *b32, uint64_t rid, int64_t hap_cov_match, int64_t hap_cov_unmatch) { if(hap->length == 0) return; uint64_t k, l, i, o, obs, ii, m_snp_stat, m_list, m_off; @@ -9985,6 +9985,11 @@ void generate_haplotypes_naive_HiFi_adv_hc(haplotype_evdience_alloc* hap, overla s = &(hap->snp_stat.a[hap->list[i].overlapSite]); assert(s->site == hap->list[i].site); zf = is_plus_sc_obs(s, st_rate, st_max, hap_cov_match, hap_cov_unmatch); + // if(overlap_list->list[hap->list[l].overlapID].y_id == 5305) { + // fprintf(stderr, "[M::%s]\t%.*s\ts->site::%u\ts->occ_0::%u\ts->occ_1::%u\ts->overlap_num::%u\n", __func__, + // (int)Get_NAME_LENGTH(R_INF, overlap_list->list[hap->list[l].overlapID].y_id), Get_NAME(R_INF, overlap_list->list[hap->list[l].overlapID].y_id), + // s->site, s->occ_0, s->occ_1, s->overlap_num); + // } // if(s->site == 11114) { // fprintf(stderr, "[M::%s] s->site::%u, s->occ_0::%u, s->occ_1::%u, s->overlap_num::%u\n", __func__, s->site, s->occ_0, s->occ_1, s->overlap_num); // } @@ -10322,7 +10327,7 @@ void generate_haplotypes_naive_HiFi_adv_back(haplotype_evdience_alloc* hap, over -void generate_haplotypes_weight(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, UC_Read* g_read, uint64_t multi_check, double st_rate, uint64_t st_max, uint64_t snp_dis, int64_t snp_cut, +void generate_haplotypes_weight(haplotype_evdience_alloc* hap, overlap_region_alloc* overlap_list, double up, uint64_t multi_check, double st_rate, uint64_t st_max, uint64_t snp_dis, int64_t snp_cut, int64_t hap_cov_match, int64_t hap_cov_unmatch) { // fprintf(stderr, "-0-[M::%s::] Done\n", __func__); @@ -12610,7 +12615,7 @@ 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, +void gen_rphase_dp_adv(haplotype_evdience_alloc *hl, overlap_region_alloc *ol, 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, asg32_v *b32) { if(hl->length <= 0) return; @@ -13584,11 +13589,283 @@ uint8_t hpc_mask_ff_region(char *sa, int64_t sn, int64_t s0, int64_t e0, int64_t return 0; } + +char qck_cigar_chk(overlap_region *z, char *tstr, int64_t tl, int64_t qtarget, int64_t flk_len, int64_t *wk0, int64_t *ck0, int64_t *qk0, int64_t *tk0, uint8_t *r_op) +{ + (*r_op) = ((uint8_t)-1); + if((qtarget < z->x_pos_s) || (qtarget > z->x_pos_e)) return ((char)-1); + int64_t wk = *wk0, ck = *ck0, qk = *qk0, tk = *tk0, wn = z->w_list.n, cn = -1; uint16_t op; + if(wn <= 0) return ((char)-1); + + if((wk >= 0) && (wk < wn) && (qtarget >= z->w_list.a[wk].x_start) && (qtarget <= z->w_list.a[wk].x_end)) { + if(is_ualn_win(z->w_list.a[wk])) { + return ((char)-1); + } + cn = z->w_list.a[wk].clen; + if((ck < 0) || (ck > cn)) {//(*ck) == cn is allowed + ck = 0; qk = z->w_list.a[wk].x_start; tk = z->w_list.a[wk].y_start; + } + } else { + if((wk < 0) || (wk >= wn)) wk = 0; + + if(qtarget >= z->w_list.a[wk].x_start) { + for (; wk < wn; wk++) { + if((qtarget >= z->w_list.a[wk].x_start) && (qtarget <= z->w_list.a[wk].x_end)) break; + } + } else { + for (; wk >= 0; wk--) { + if((qtarget >= z->w_list.a[wk].x_start) && (qtarget <= z->w_list.a[wk].x_end)) break; + } + } + + if((wk >= 0) && (wk < wn) && (qtarget >= z->w_list.a[wk].x_start) && (qtarget <= z->w_list.a[wk].x_end)) { + if(is_ualn_win(z->w_list.a[wk])) { + (*wk0) = wk; (*ck0) = ck; (*qk0) = qk; (*tk0) = tk; + return ((char)-1); + } + cn = z->w_list.a[wk].clen; + if((ck < 0) || (ck > cn) || (wk != (*wk0))) {//(*ck) == cn is allowed + ck = 0; qk = z->w_list.a[wk].x_start; tk = z->w_list.a[wk].y_start; + } + } else { + (*wk0) = wk; (*ck0) = ck; (*qk0) = qk; (*tk0) = tk; + return ((char)-1); + } + } + + bit_extz_t ez; set_bit_extz_t(ez, (*z), wk); + if(!ez.cigar.n) return ((char)-1); + + int64_t ws, we, ttarget = -1; cn = ez.cigar.n; + + if((ck < 0) || (ck > cn)) {//(*ck) == cn is allowed + ck = 0; qk = ez.ts; tk = ez.ps; + } + + while ((ck > 0) && (qk > qtarget)) {///x -> t; y -> p + --ck; + op = ez.cigar.a[ck]>>14; + if(op!=2) qk -= (ez.cigar.a[ck]&(0x3fff)); + if(op!=3) tk -= (ez.cigar.a[ck]&(0x3fff)); + } + + // int64_t qs = qtarget - flk_len, qe = qtarget + flk_len, os, oe; + while (ck < cn && qk <= qtarget) {//[s, e) + ws = qk; + op = ez.cigar.a[ck]>>14; + if(op!=2) qk += (ez.cigar.a[ck]&(0x3fff)); + if(op!=3) tk += (ez.cigar.a[ck]&(0x3fff)); + ck++; we = qk; + + // os = MAX(qs, ws); oe = MIN(qe, we); + // ovlp = ((oe>os)? (oe-os):0); + + if(op != 0 && op != 1) continue;///only collect match/snp + + if((qtarget >= ws) && (qtarget < we)) { + (*r_op) = op; ttarget = qtarget-qk+tk; + (*wk0) = wk; (*ck0) = ck; (*qk0) = qk; (*tk0) = tk; + if((op == 1) && (tl > ttarget) && (ttarget >= 0)) { + + + + + + + return tstr[ttarget]; + } + return ((char)-1); + } + } + (*wk0) = wk; (*ck0) = ck; (*qk0) = qk; (*tk0) = tk; + return ((char)-1); +} + + +int push_info_ssl(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, haplotype_evdience* u_a, overlap_region *oa, asg8_v *v8, uint64_t tot_cov, uint64_t scw, uint64_t tcut, + char r1a, overlap_region *rchn, char *tstr, int64_t tlen, int64_t flk_len, int64_t *wk, int64_t *ck, int64_t *qk, int64_t *tk) +{ + if(a_n < 2) return 0; + uint64_t i, hi = 0, k, m, occ_0, occ_1[6], occ_2, diff, rev_n; char r0a = r1a; uint8_t ihpc = 0; const uint32_t HQ_MASK = 1u << 31; const uint32_t OD_MASK = ~HQ_MASK; + haplotype_evdience at; haplotype_evdience *ra = NULL; uint8_t op; uint32_t ft; + occ_0 = occ_2 = diff = rev_n = 0; memset(occ_1, 0, sizeof(uint64_t)*6); + + radix_sort_haplotype_evdience_id_srt(a, a + a_n); + if((rchn) && (tstr) && (tlen > 0)) { + r0a = qck_cigar_chk(rchn, tstr, tlen, a[0].site, flk_len, wk, ck, qk, tk, &op); + if((op != 1) || (r0a == ((char)-1))) r0a = r1a; + if(r0a != r1a) { + for (k = 0; k < a_n; k++) { + if((r0a == a[k].misBase) && ((k == 0) || (a[k].overlapID != a[k-1].overlapID))) break; + } + if(k >= a_n) r0a = r1a; + } + } + + if(tcut == ((uint64_t)-1)) { + for (k = 1, m = i = 0; k <= a_n; k++) { + if((k == a_n) || (a[k].overlapID != a[m].overlapID)) { + a[i] = a[m]; + + if(r0a == a[i].misBase) { + occ_0 += a[i].cov; + if((oa) && (oa[a[i].overlapID].y_pos_strand == 0)) { + rev_n += a[i].cov; + } + if((!ihpc) && (hh_hp(a[i]))) ihpc = 1; + ft = 0; + } else { + occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]] += a[i].cov; + diff += a[i].cov; + ft = 1; + } + if(hh_tp(a[i]) != ft) a[i].type ^= ((uint32_t)1); + + occ_2 += a[i].cov; + i++; m = k; + } + } + } else if(oa) { + for (k = 1, m = i = hi = 0; k <= a_n; k++) { + if((k == a_n) || (a[k].overlapID != a[m].overlapID)) { + a[i] = a[m]; + + if(r0a == a[i].misBase){ + occ_0 += a[i].cov; + if(oa[a[i].overlapID].y_pos_strand == 0) rev_n += a[i].cov; + if((!ihpc) && (hh_hp(a[i]))) ihpc = 1; + ft = 0; + } else { + occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]] += a[i].cov; + diff += a[i].cov; + ft = 1; + } + if(hh_tp(a[i]) != ft) a[i].type ^= ((uint32_t)1); + + occ_2 += a[i].cov; + if(oa[a[i].overlapID].y_id >= tcut) { + if(hi != i) { + at = a[i]; a[i] = a[hi]; a[hi] = at; + } + a[hi].overlapID |= HQ_MASK; hi++; + } + i++; m = k; + } + } + + if(i > hi + 1) { + ra = a + hi; hi = i - hi; + } + } + + a_n = i; + + // for (i = 0; i < a_n; i++) { + // if(hh_tp(a[i]) == 0){ + // occ_0 += a[i].cov; + // if((oa) && (oa[a[i].overlapID].y_pos_strand == 0)) { + // rev_n += a[i].cov; + // } + // if((!ihpc) && (hh_hp(a[i]))) ihpc = 1; + // }else if(hh_tp(a[i]) == 1){ + // occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]] += a[i].cov; + // diff += a[i].cov; + // } + // occ_2 += a[i].cov; + // } + + + /** + 1. if occ_0 = 0, that means all overlaps are different with this read at this site + 2. it is not possible that occ_1 = 0, + 3. if occ_1 = 1, there are only one difference. It must be a sequencing error. + (for repeat, it maybe a snp at repeat. but ...) + **/ + SnpStats *p = NULL; + if(occ_0 == 0 || diff <= 1) return 0; + if((oa) && (rev_n == occ_0)) return 0; + if(!oa) { + rev_n = occ_2; + } else { + rev_n++; + } + + for (i = m = 0; i < 4; i++) { + if(occ_1[i] >= 2){ + kv_pushp(SnpStats, h->snp_stat, &p); + p->id = h->snp_stat.n-1; + p->occ_0 = 1 + occ_0; + p->occ_1 = occ_1[i]; + if((tot_cov > 0) && (tot_cov >= (p->occ_0 + p->occ_1))) { + p->occ_2 = tot_cov - p->occ_0 - p->occ_1; + } else { + p->occ_2 = occ_2 + 1 - p->occ_0 - p->occ_1; + } + // fprintf(stderr, "[M::%s]\ttot_cov::%lu\tocc_2::%lu\n", __func__, tot_cov, occ_2); + p->site = a[0].site; + p->score = -1; + p->overlap_num = rev_n; + p->is_homopolymer = ihpc; + occ_1[i] = p->id; + m++; + // fprintf(stderr, "[M::%s]\tsite::%u\tn0::%u\tn1::%u\n", __func__, p->site, p->occ_0, p->occ_1); + } else { + occ_1[i] = (uint64_t)-1; + } + } + occ_1[4] = occ_1[5] = (uint64_t)-1; + if(m == 0) return 0; + + if(ra) { + radix_sort_haplotype_evdience_id_srt(ra, ra + hi); m = hi>>1; + for (k = 0; k < m; k++) { + at = ra[k]; ra[k] = ra[hi-k-1]; ra[hi-k-1] = at; + } + } + + if(!v8) { + for (i = m = 0; i < a_n; i++) { + // fprintf(stderr, "[M::%s] a[%lu].misBase->%c\n", __func__, i, a[i].misBase); + if(hh_tp(a[i]) == 0) { + a[i].overlapSite = h->snp_stat.n-1; + } else if(occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]]!=(uint64_t)-1){ + a[i].cov = a[i].overlapSite;///note: only renew cov here!!! + a[i].overlapSite = occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]]; + } else { + continue; + } + u_a[m++] = a[i]; + } + } else { + uint64_t tpos, tqual, wqual; + for (i = m = 0; i < a_n; i++) { + tpos = a[i].overlapSite; + if(hh_tp(a[i]) == 0) { + a[i].overlapSite = h->snp_stat.n-1; + } else if(occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]]!=(uint64_t)-1){ + a[i].cov = a[i].overlapSite;///note: only renew cov here!!! + a[i].overlapSite = occ_1[seq_nt6_table[(uint8_t)(a[i].misBase)]]; + } else { + continue; + } + tqual = wqual = 3; + if(oa[a[i].overlapID&OD_MASK].y_id < tcut) get_wqual(oa[a[i].overlapID&OD_MASK].y_id, tpos, oa[a[i].overlapID&OD_MASK].y_pos_strand, v8, NULL, scw, &tqual, &wqual); + a[i].type |= (tqual<<2); a[i].type |= (wqual<<(sc_bn+2)); + // assert(hh_bq(a[i]) == tqual); assert(hh_wq(a[i]) == wqual); + // fprintf(stderr, "[M::%s] site::%u, type::%u, occ0::%lu, occ2::%lu, tqual::%lu, wqual::%lu\n", __func__, a[i].site, hh_tp(a[i]), occ_0, occ_2, tqual, wqual); + u_a[m++] = a[i]; + } + } + + + return m; +} + int push_info(haplotype_evdience_alloc* h, haplotype_evdience* a, uint64_t a_n, haplotype_evdience* u_a, overlap_region *oa, asg8_v *v8, uint64_t tot_cov, uint64_t scw, uint64_t tcut) { uint64_t i, hi = 0, k, m, occ_0, occ_1[6], occ_2, diff, rev_n; uint8_t ihpc = 0; const uint32_t HQ_MASK = 1u << 31; const uint32_t OD_MASK = ~HQ_MASK; haplotype_evdience at; haplotype_evdience *ra = NULL; occ_0 = occ_2 = diff = rev_n = 0; memset(occ_1, 0, sizeof(uint64_t)*6); - + radix_sort_haplotype_evdience_id_srt(a, a + a_n); if(tcut == ((uint64_t)-1)) { for (k = 1, m = i = 0; k <= a_n; k++) { @@ -13789,7 +14066,7 @@ void partition_overlaps_advance(overlap_region_alloc* overlap_list, All_reads* R hap->length = m; // generate_haplotypes_naive_advance(hap, overlap_list, NULL); - generate_haplotypes_naive_HiFi(hap, overlap_list, 0.04, g_read, 1, 0, ((uint64_t)-1), asm_opt.s_hap_cov, asm_opt.infor_cov); + generate_haplotypes_naive_HiFi(hap, overlap_list, 0.04, 1, 0, ((uint64_t)-1), asm_opt.s_hap_cov, asm_opt.infor_cov); // generate_haplotypes_DP(hap, overlap_list, R_INF, g_read->length, force_repeat); // generate_haplotypes_naive(hap, overlap_list, R_INF, g_read->length, force_repeat); @@ -15953,6 +16230,9 @@ uint32_t push_hc_wlst_exz(const ul_idx_t *uref, hpc_t *hpc_g, All_reads *rref, o ol->align_length += qe + 1 - qs; ovl = ol->x_pos_e+1-ol->x_pos_s; ualn = (qe + 1 - ol->x_pos_s) - ol->align_length; aln = ovl-ualn; + // fprintf(stderr, "[M::%s::]\ttid::%u(%u)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\tovl::%ld\tcur_aln::%u\tun_aln::%ld\tinfer_aln::%ld\tovlp_cut::%f\n", __func__, + // ol->y_id, ol->y_pos_strand, (int)Get_NAME_LENGTH(R_INF, ol->y_id), Get_NAME(R_INF, ol->y_id), ol->x_pos_s, ol->x_pos_e + 1, ol->y_pos_s, ol->y_pos_e + 1, + // ovl, ol->align_length, ualn, aln, ovlp_cut); // if((!force_aln) && (!simi_pass(ovl, aln, 0, ovlp_cut, &e_rate)) && (!simi_pass(ovl, aln, sec_check, ovlp_cut, NULL))) { if((!force_aln) && (!pass_qovlp(ovl, aln, ovlp_cut))) { rf = 0; @@ -18926,6 +19206,260 @@ int64_t cal_estimate_err_hc(overlap_region *z, int64_t wl, int64_t qs, int64_t q return tot; } +int64_t cal_estimate_err_scc(asg16_v *scc, int64_t qs, int64_t qe, int64_t ts, int64_t te, int64_t *qk0, int64_t *tk0, int64_t *ck0) +{ + int64_t qk, tk, ck, cn, os, oe, ovlp, tot = 0, wq[2], wt[2]; uint16_t op, bq, bt; uint32_t cl; + if(qe <= qs || te <= ts) return (((qe>qs)?(qe-qs):(0)) + ((te>ts)?(te-ts):(0))); + + qk = *qk0; tk = *tk0; ck = *ck0; cn = scc->n; + + if((ck < 0) || (ck > cn)) {//(*ck) == cn is allowed + ck = qk = tk = 0; + } + + if((ck > 0) && ((qk > qs) || (tk > ts))) { + while ((ck > 0) && ((qk > qs) || (tk > ts))) { + wq[0] = qk; wt[0] = tk; + ck = pop_trace_bp_rev_f(scc, ck-1, &op, &bq, &bt, &cl)+1; + if(op != 2) qk -= cl; + if(op != 3) tk -= cl; + wq[1] = qk; wt[1] = tk; + } + } + + while ((ck < cn) && (qk < qe || tk < te)) {//[s, e) + wq[0] = qk; wt[0] = tk; + ck = pop_trace_bp_f(scc, ck, &op, &bq, &bt, &cl); + if(op != 2) qk += cl; + if(op != 3) tk += cl; + wq[1] = qk; wt[1] = tk; + + if(op == 0) { + continue; + } else if(op == 1 || op == 3) { + os = MAX(qs, wq[0]); oe = MIN(qe, wq[1]); + ovlp = ((oe>os)? (oe-os):0); + tot += ovlp; + } else if(op == 2) { + os = MAX(ts, wt[0]); oe = MIN(te, wt[1]); + ovlp = ((oe>os)? (oe-os):0); + tot += ovlp; + } + } + *qk0 = qk; *tk0 = tk; *ck0 = ck; + + cn = qe - qs; if(cn < (te - ts)) cn = te - ts; + if(tot > cn) tot = cn; + return tot; +} + +inline void push_khit(Candidates_list *res, int32_t xs, int32_t ys, uint32_t len, uint32_t h_khit, uint32_t *ic) +{ + uint32_t p, c; k_mer_hit *z; + c = ((len >= h_khit)?1:2); if(ic) c = *ic; c <<= 8; + if(len > 0) { + while (len >= (0xffu)) { + p = (c + (0xffu)); + kv_pushp_cl(k_mer_hit, (*res), &z); + memset(z, 0, sizeof((*z))); + z->self_offset = xs; z->offset = ys; z->cnt = p; + len -= (0xffu); + } + if(len) { + p = (c + len); + kv_pushp_cl(k_mer_hit, (*res), &z); + memset(z, 0, sizeof((*z))); + z->self_offset = xs; z->offset = ys; z->cnt = p; + } + } else { + p = (c + len); + kv_pushp_cl(k_mer_hit, (*res), &z); + memset(z, 0, sizeof((*z))); + z->self_offset = xs; z->offset = ys; z->cnt = p; + } +} + +inline uint8_t push_srt_khit(Candidates_list *res, int64_t *uk, int64_t uas, int64_t uae, int64_t pq, int64_t pt, int64_t l, int64_t khit) +{ + uint8_t ff = 1; int64_t kl, ml; + for (; ((*uk) <= uae) && (res->list[(*uk)].self_offset < pq); (*uk)++) { + kl = (res->list[(*uk)].cnt&((uint32_t)(0xffu))); if(kl > khit) kl = khit; + push_khit(res, res->list[(*uk)].self_offset, res->list[(*uk)].offset, kl, khit, NULL); + } + + for (; ((*uk) <= uae) && (res->list[(*uk)].self_offset == pq) && (res->list[(*uk)].offset < pt); (*uk)++) { + kl = (res->list[(*uk)].cnt&((uint32_t)(0xffu))); if(kl > khit) kl = khit; + push_khit(res, res->list[(*uk)].self_offset, res->list[(*uk)].offset, kl, khit, NULL); + } + + for (ml = -1; ((*uk) <= uae) && (res->list[(*uk)].self_offset == pq) && (res->list[(*uk)].offset == pt); (*uk)++) { + ff = 0; + kl = res->list[*uk].cnt & 0xffu; + if(ml < kl) ml = kl; + } + + if(!ff) { + kl = ml; if(kl > khit) kl = khit; if(kl < l) kl = l; + push_khit(res, pq, pt, kl, khit, NULL); + } else { + push_khit(res, pq, pt, l, khit, NULL); + } + + for (; ((*uk) <= uae) && (res->list[(*uk)].self_offset == pq); (*uk)++) { + kl = (res->list[(*uk)].cnt&((uint32_t)(0xffu))); if(kl > khit) kl = khit; + push_khit(res, res->list[(*uk)].self_offset, res->list[(*uk)].offset, kl, khit, NULL); + } + + return ff; +} + +int64_t extract_exact_cigar_scc(asg16_v *scc, int64_t qs, int64_t qe, int64_t ts, int64_t te, int64_t qk, int64_t tk, int64_t ck, int32_t minl, int64_t h_khit, Candidates_list *res, + int64_t uas, int64_t uae) ///[uas, uae] +{ + int64_t kl, cn = scc->n, wq[2], wt[2], oq[2], ot[2], occ = 0, l, pq, pt, uk = uas, os, oe; uint16_t op, bq, bt; uint32_t cl; + if(qe <= qs || te <= ts) return 0; + + if((ck < 0) || (ck > cn)) {//(*ck) == cn is allowed + ck = qk = tk = 0; + } + + if((ck > 0) && ((qk > qs) || (tk > ts))) { + while ((ck > 0) && ((qk > qs) || (tk > ts))) { + wq[0] = qk; wt[0] = tk; + ck = pop_trace_bp_rev_f(scc, ck-1, &op, &bq, &bt, &cl)+1; + if(op != 2) qk -= cl; + if(op != 3) tk -= cl; + wq[1] = qk; wt[1] = tk; + } + } + + while ((ck < cn) && (qk < qe || tk < te)) {//[s, e) + wq[0] = qk; wt[0] = tk; + ck = pop_trace_bp_f(scc, ck, &op, &bq, &bt, &cl); + if(op != 2) qk += cl; + if(op != 3) tk += cl; + wq[1] = qk; wt[1] = tk; + + if(op != 0) continue;///exact match + + oq[0] = MAX(qs, wq[0]); oq[1] = MIN(qe, wq[1]); + l = oq[1] - oq[0]; + if(l < minl) continue; + + ot[0] = wt[0] + oq[0] - wq[0]; ot[1] = wt[0] + oq[1] - wq[0]; + os = MAX(ts, ot[0]);oe = MIN(te, ot[1]); + l = oe - os; + if(l < minl) continue; + + /** clip query consistently with target clipping **/ + oq[0] += os - ot[0]; + oq[1] = oq[0] + l; + ot[0] = os; ot[1] = oe; + + if(l <= h_khit) { + pq = oq[1] - 1; pt = ot[1] - 1; + if(push_srt_khit(res, &uk, uas, uae, pq, pt, l, h_khit)) occ++; + } else { + pq = oq[1] - h_khit; pt = ot[1] - h_khit; + while ((pq - minl >= oq[0]) && (pt - minl >= ot[0])) { + pq -= minl; pt -= minl; + } + + pq = pq + h_khit - 1; pt = pt + h_khit - 1; + for (; pq < oq[1] && pt < ot[1]; pq += minl, pt += minl) { + if(push_srt_khit(res, &uk, uas, uae, pq, pt, h_khit, h_khit)) occ++; + } + } + } + + while(uk <= uae) { + kl = (res->list[uk].cnt&((uint32_t)(0xffu))); if(kl > h_khit) kl = h_khit; + push_khit(res, res->list[uk].self_offset, res->list[uk].offset, kl, h_khit, NULL); + uk++; + } + + return occ; +} + +int64_t extract_exact_cigar_scc_rev(asg16_v *scc, int64_t qs, int64_t qe, int64_t ts, int64_t te, int64_t qk, int64_t tk, int64_t ck, int32_t minl, int64_t h_khit, Candidates_list *res, + int64_t uas, int64_t uae) ///[uas, uae] +{ + int64_t kl, cn = scc->n, wq[2], wt[2], oq[2], ot[2], occ = 0, l, pq, pt, uk = uas, os, oe; uint16_t op, bq, bt; uint32_t cl; + if(qe <= qs || te <= ts) return 0; + + if((ck < 0) || (ck > cn)) {//(*ck) == cn is allowed + ck = qk = tk = 0; + } + + if((ck > 0) && ((qk > qs) || (tk > ts))) { + while ((ck > 0) && ((qk > qs) || (tk > ts))) { + wq[0] = qk; wt[0] = tk; + ck = pop_trace_bp_rev_f(scc, ck-1, &op, &bq, &bt, &cl)+1; + /** + if(op != 2) qk -= cl; + if(op != 3) tk -= cl; + **/ + if(op != 2) tk -= cl; + if(op != 3) qk -= cl; + wq[1] = qk; wt[1] = tk; + } + } + + while ((ck < cn) && (qk < qe || tk < te)) {//[s, e) + wq[0] = qk; wt[0] = tk; + ck = pop_trace_bp_f(scc, ck, &op, &bq, &bt, &cl); + /** + if(op != 2) qk += cl; + if(op != 3) tk += cl; + **/ + if(op != 2) tk += cl; + if(op != 3) qk += cl; + wq[1] = qk; wt[1] = tk; + + // if(op == 0) fprintf(stderr, "[M::%s::]\twq::[%ld,%ld)\twt::[%ld,%ld)\top::%u\n", __func__, wq[0], wq[1], wt[0], wt[1], op); + + if(op != 0) continue;///exact match + + oq[0] = MAX(qs, wq[0]); oq[1] = MIN(qe, wq[1]); + l = oq[1] - oq[0]; + if(l < minl) continue; + + ot[0] = wt[0] + oq[0] - wq[0]; ot[1] = wt[0] + oq[1] - wq[0]; + os = MAX(ts, ot[0]);oe = MIN(te, ot[1]); + l = oe - os; + if(l < minl) continue; + + /** clip query consistently with target clipping **/ + oq[0] += os - ot[0]; + oq[1] = oq[0] + l; + ot[0] = os; ot[1] = oe; + + + if(l <= h_khit) { + pq = oq[1] - 1; pt = ot[1] - 1; + if(push_srt_khit(res, &uk, uas, uae, pq, pt, l, h_khit)) occ++; + } else { + pq = oq[1] - h_khit; pt = ot[1] - h_khit; + while ((pq - minl >= oq[0]) && (pt - minl >= ot[0])) { + pq -= minl; pt -= minl; + } + + pq = pq + h_khit - 1; pt = pt + h_khit - 1; + for (; pq < oq[1] && pt < ot[1]; pq += minl, pt += minl) { + if(push_srt_khit(res, &uk, uas, uae, pq, pt, h_khit, h_khit)) occ++; + } + } + } + + while(uk <= uae) { + kl = (res->list[uk].cnt&((uint32_t)(0xffu))); if(kl > h_khit) kl = h_khit; + push_khit(res, res->list[uk].self_offset, res->list[uk].offset, kl, h_khit, NULL); + uk++; + } + + return occ; +} + ///[s, e); [ps, pe) inline char *retrieve_str_seq_exz(UC_Read *tu, int64_t s, int64_t l, @@ -19468,7 +20002,7 @@ void push_alnw(overlap_region *aux_o, bit_extz_t *exz) // exz->ts, exz->te, exz->ps, exz->pe, exz->err, p->clen); // assert((p->x_endts)&&(p->y_endps)); - if(p->clen > 0) { + if((p->clen > 0) && (p->cidx + p->clen == aux_o->w_list.c.n)) { t = ((int64_t)p->error) + ((int64_t)exz->err); ///note: t cannot be equal to INT16_MAX; otherwise it is unable to distiguish unaligned regions if(((p->x_end+1) == exz->ts) && ((p->y_end+1) == exz->ps) && (t < INT16_MAX)) { @@ -19502,6 +20036,25 @@ void push_unmap_alnw(overlap_region *aux_o, int32_t qs, int64_t qe, int64_t ts, p->cidx = p->clen = 0; } +///[qs, qe] && [ts, te] +void push_unmap_alnw_hold(overlap_region *aux_o, int32_t qs, int64_t qe, int64_t ts, int64_t te, int64_t mode, uint8_t hold, int64_t est_err) +{ + window_list *p = NULL; + kv_pushp(window_list, aux_o->w_list, &p); + p->x_start = qs; p->x_end = qe; + p->y_start = ts; p->y_end = te; + p->error_threshold = mode; p->error = INT16_MAX; + p->extra_begin = p->extra_end = -1; + if((hold) && (est_err < INT16_MAX)) { + p->extra_begin -= est_err; + p->extra_end -= est_err; + if((p->extra_begin >= -1) || (p->extra_end >= -1)) { + p->extra_begin = p->extra_end = -1; + } + } + p->cidx = p->clen = 0; +} + ///[qs, qe] && [ts, te] void push_replace_alnw(overlap_region *aux_o, int32_t qs, int64_t qe, int64_t ts, int64_t te, int64_t mode) { @@ -20341,32 +20894,6 @@ bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, u } -inline void push_khit(Candidates_list *res, int32_t xs, int32_t ys, uint32_t len, uint32_t h_khit, uint32_t *ic) -{ - uint32_t p, c; k_mer_hit *z; - c = ((len >= h_khit)?1:2); if(ic) c = *ic; c <<= 8; - if(len > 0) { - while (len >= (0xffu)) { - p = (c + (0xffu)); - kv_pushp_cl(k_mer_hit, (*res), &z); - memset(z, 0, sizeof((*z))); - z->self_offset = xs; z->offset = ys; z->cnt = p; - len -= (0xffu); - } - if(len) { - p = (c + len); - kv_pushp_cl(k_mer_hit, (*res), &z); - memset(z, 0, sizeof((*z))); - z->self_offset = xs; z->offset = ys; z->cnt = p; - } - } else { - p = (c + len); - kv_pushp_cl(k_mer_hit, (*res), &z); - memset(z, 0, sizeof((*z))); - z->self_offset = xs; z->offset = ys; z->cnt = p; - } -} - uint32_t extract_exact_cigar(asg16_v *ez, int32_t ps, int32_t ts, int32_t pmin, int32_t pmax, int32_t tmin, int32_t tmax, Candidates_list *res, int32_t minl, int64_t min_w_l, int64_t h_khit) { @@ -20817,17 +21344,21 @@ void merge_hc_fast_cigar0(overlap_region *des, overlap_region *src) int64_t k = 0, i, wn = src->w_list.n, ck_s, ck_e, zk, e0, e, e_tot = 0, mz; uint32_t op0 = (uint16_t)-1, ol0, op, ol, xs, ys, xk, yk, xe, ye; asg16_v b16; window_list *p; copy_asg_arr(b16, des->w_list.c); b16.n = 0; des->w_list.n = 0; + kv_resize(window_list, des->w_list, src->w_list.n); // fprintf(stderr, "\n[M::%s::aln::-xid->%u]\n", __func__, src->x_id); // /**if(z->y_id == 310315)**/ { - // for (i = 0; i < wn; i++) { - // fprintf(stderr, "[aln::-i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, cidx::%u, clen::%u, mode::%d\n", i, - // src->w_list.a[i].x_end+1-src->w_list.a[i].x_start, - // src->w_list.a[i].x_start, src->w_list.a[i].x_end+1, - // src->w_list.a[i].y_start, src->w_list.a[i].y_end+1, - // src->w_list.a[i].error, src->w_list.a[i].cidx, src->w_list.a[i].clen, src->w_list.a[i].error_threshold); - // } + // fprintf(stderr, "[M::%s::] des->w_list.n::%ld, src->w_list.n::%ld\n", + // __func__, (int64_t)des->w_list.n, (int64_t)src->w_list.n); + // for (i = 0; i < wn; i++) { + // fprintf(stderr, "[aln::-i->%ld::ql->%d] q::[%d, %d), t::[%d, %d), err::%d, cidx::%u, clen::%u, mode::%d, y_id::%u\n", i, + // src->w_list.a[i].x_end+1-src->w_list.a[i].x_start, + // src->w_list.a[i].x_start, src->w_list.a[i].x_end+1, + // src->w_list.a[i].y_start, src->w_list.a[i].y_end+1, + // src->w_list.a[i].error, src->w_list.a[i].cidx, src->w_list.a[i].clen, src->w_list.a[i].error_threshold, + // src->y_id); + // } // } while (k < wn) { @@ -21347,6 +21878,260 @@ void hc_ovlp_base_non_retrieve_direct(overlap_region *z, k_mer_hit *ch_a, int64_ } } + + +int64_t gen_scc_non_retrieve_chain( asg16_v* scc, uint8_t is_scc_rev, Candidates_list *cl, int64_t qs, int64_t qe, int64_t ts, int64_t te, + int64_t ql, int64_t tl, double e_rate, int64_t h_khit, int64_t mode, int64_t is_accurate, + int64_t qk, int64_t tk, int64_t ck, int64_t uas, int64_t uae) +{ + uint32_t occ = 0; int64_t rcn = cl->length, ncn; + ///global or forward + // uint32_t w = 1; + // if(mode == 0 || mode == 1) push_khit(cl, qs, ts, 0, 0, &w); + if(is_scc_rev) { + occ = extract_exact_cigar_scc_rev(scc, qs, qe, ts, te, qk, tk, ck, 10, h_khit, cl, uas, uae); + } else { + occ = extract_exact_cigar_scc(scc, qs, qe, ts, te, qk, tk, ck, 10, h_khit, cl, uas, uae); + } + + ///global or backward + // if(mode == 0 || mode == 2) push_khit(cl, qe-1, te-1, 0, 0, &w); + + if(!occ) { + cl->length = rcn; return 0; + } + ncn = cl->length; cl->length = rcn; + k_mer_hit *ch_a = cl->list + rcn; int64_t ch_n0 = ncn - rcn, ch_n; + int64_t max_skip, max_iter, max_dis, quick_check, k; double chn_pen_gap, chn_pen_skip; + set_lchain_dp_op(is_accurate, h_khit, &max_skip, &max_iter, &max_dis, &chn_pen_gap, &chn_pen_skip, &quick_check); + // max_dis = MAX_SIN_L>>1; + ch_n = lchain_qdp_fix_adv(ch_a, ch_n0, &(cl->chainDP), max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, + e_rate, ql, tl, quick_check, ((mode==0)||(mode==1))?1:0, ((mode==0)||(mode==2))?1:0, ch_a); + + if((ch_n > 0) && (ch_n == (uae + 1 - uas))) { + for (k = uas; k <= uae; k++) { + if(cl->list[k].offset != ch_a[k - uas].offset) return ch_n; + if(cl->list[k].self_offset != ch_a[k - uas].self_offset) return ch_n; + } + return 0; + } + + return ch_n; +} + +uint8_t hc_ovlp_base_non_retrieve_direct_scc(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, Candidates_list *cl, int64_t h_khit, char* qstr, char* tstr, + int64_t scc_qs, int64_t scc_ts, int64_t scc_ck, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, asg16_v* scc, uint8_t is_scc_rev, int32_t global_mode, uint64_t rid); + +void recal_non_retrieve_direct_scc(overlap_region *z, Candidates_list *cl, overlap_region *aux_o, int64_t aux_ei, + char* qstr, char* tstr, bit_extz_t *exz, double e_rate, int64_t ql, int64_t tl, int64_t h_khit, asg16_v* scc, uint8_t is_scc_rev, int64_t force_l, + int64_t qk, int64_t tk, int64_t ck, int64_t uas, int64_t uae) +{ + if((aux_ei < 0) || (aux_ei >= ((int64_t)aux_o->w_list.n))) return; + int64_t rcn = cl->length, ch_n, qs, qe, ts, te, mode, an0 = aux_o->w_list.n, cn0 = aux_o->w_list.c.n, an, mm, k, aux_si = INT64_MAX, err0, err1; + k_mer_hit *ch_a = NULL; uint8_t m[2]; + + for (k = aux_ei; (k >= 0) && (is_ualn_win(aux_o->w_list.a[k])); k--){;} aux_si = k + 1; + if(aux_si > aux_ei) return; + + ///[qs, qe) && [ts, te) + qs = aux_o->w_list.a[aux_si].x_start; qe = aux_o->w_list.a[aux_ei].x_end+1; + ts = aux_o->w_list.a[aux_si].y_start; te = aux_o->w_list.a[aux_ei].y_end+1; + + if(qe - qs < force_l) return; + + + mode = 3; + m[0] = aux_o->w_list.a[aux_si].error_threshold; m[1] = aux_o->w_list.a[aux_ei].error_threshold; + if(aux_si == aux_ei) { + mode = m[0]; + } else { + ///case 1: m[0] = 0 && m[1] = 0 + ///case 2: m[0] = 0 && m[1] = 1 + ///case 3: m[0] = 2 && m[1] = 0 + ///case 4: m[0] = 2 && m[1] = 1 + mode = m[0] + m[1]; + + assert(m[0] == 0 || m[0] == 2); + assert(m[1] == 0 || m[1] == 1); + } + + + ch_n = gen_scc_non_retrieve_chain(scc, is_scc_rev, cl, qs, qe, ts, te, ql, tl, e_rate, h_khit, mode, 1, qk, tk, ck, uas, uae); + if(ch_n) { + ch_a = cl->list + rcn; + if(mode == 0) { + assert((((int64_t)ch_a[0].offset) == aux_o->w_list.a[aux_si].y_start) && (((int64_t)ch_a[0].self_offset) == aux_o->w_list.a[aux_si].x_start)); + assert((((int64_t)ch_a[ch_n-1].offset) == (aux_o->w_list.a[aux_ei].y_end+1)) && (((int64_t)ch_a[ch_n-1].self_offset) == (aux_o->w_list.a[aux_ei].x_end+1))); + } else if(mode == 1) { + assert((((int64_t)ch_a[0].offset) == aux_o->w_list.a[aux_si].y_start) && (((int64_t)ch_a[0].self_offset) == aux_o->w_list.a[aux_si].x_start)); + } else if(mode == 2) { + assert((((int64_t)ch_a[ch_n-1].offset) == (aux_o->w_list.a[aux_ei].y_end+1)) && (((int64_t)ch_a[ch_n-1].self_offset) == (aux_o->w_list.a[aux_ei].x_end+1))); + } + hc_ovlp_base_non_retrieve_direct_scc(z, ch_a, ch_n, cl, h_khit, qstr, tstr, qk, tk, ck, exz, aux_o, e_rate, ql, tl, scc, 1, mode, (uint64_t)-1); + an = aux_o->w_list.n; + if(an > an0) { + err0 = qe - qs; if(err0 < te - ts) err0 = te - ts; + for (k = an0, err1 = 0; k < an; k++) { + if(is_ualn_win(aux_o->w_list.a[k])) { + mm = aux_o->w_list.a[k].x_end + 1 - aux_o->w_list.a[k].x_start; + if(mm < (aux_o->w_list.a[k].y_end + 1 - aux_o->w_list.a[k].y_start)) { + mm = aux_o->w_list.a[k].y_end + 1 - aux_o->w_list.a[k].y_start; + } + err1 += mm; + } else { + err1 += aux_o->w_list.a[k].error; + } + } + + if(err1 < err0) { + ///a::[aux_si, aux_ei + 1) -> removed + ///b::[aux_ei + 1, an0) + ///c::[an0, an) + int64_t wq[2], wt[2], bn = an0 - (aux_ei + 1), cn = an - an0; + wq[0] = aux_o->w_list.a[an0].x_start; wq[1] = aux_o->w_list.a[an-1].x_end + 1; + wt[0] = aux_o->w_list.a[an0].y_start; wt[1] = aux_o->w_list.a[an-1].y_end + 1; + + if(aux_si > 0) { + wq[0] = aux_o->w_list.a[aux_si-1].x_end+1; + wt[0] = aux_o->w_list.a[aux_si-1].y_end+1; + } + + if(aux_ei + 1 < an0) { + wq[1] = aux_o->w_list.a[aux_ei + 1].x_start; + wt[1] = aux_o->w_list.a[aux_ei + 1].y_start; + } + + assert((wq[0] == aux_o->w_list.a[an0].x_start) && (wq[1] == aux_o->w_list.a[an-1].x_end + 1) + && (wt[0] == aux_o->w_list.a[an0].y_start) && (wt[1] == aux_o->w_list.a[an-1].y_end + 1)); + /** function for a, b, c**/ + if(bn) { + kv_resize(window_list, aux_o->w_list, (uint64_t)(an + bn)); + memcpy(aux_o->w_list.a + an, aux_o->w_list.a + aux_ei + 1, sizeof(*(aux_o->w_list.a))*bn); + } + + if(an > an0) { + for (k = an0, mm = aux_si; k < an; k++, mm++) { + aux_o->w_list.a[mm] = aux_o->w_list.a[k]; + } + } + + if(bn) { + memcpy(aux_o->w_list.a + aux_si + cn, aux_o->w_list.a + an, sizeof(*(aux_o->w_list.a))*bn); + } + + aux_o->w_list.n = aux_si + bn + cn; + return; + } + } + } + aux_o->w_list.n = an0; aux_o->w_list.c.n = cn0; +} + +uint8_t hc_ovlp_base_non_retrieve_direct_scc(overlap_region *z, k_mer_hit *ch_a, int64_t ch_n, Candidates_list *cl, int64_t h_khit, char* qstr, char* tstr, + int64_t scc_qs, int64_t scc_ts, int64_t scc_ck, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t ql, int64_t tl, asg16_v* scc, uint8_t is_scc_rev, int32_t global_mode, uint64_t rid) +{ + if(ch_n < 1) return 0; + int64_t i, l, ls, is, ie, mode, q[2], t[2], is_done, qk = scc_qs, tk = scc_ts, ck = scc_ck, uq, ut, uc, wn0 = -1, cidx = ch_a - cl->list, uai[2], est_err, est_ll; uint8_t umap = 0, est_hd; + ///ch_n >= 1 + uq = qk; ut = tk; uc = ck; umap = 0; uai[0] = uai[1] = INT64_MIN; + ls = -1; is = 0; ie = ch_n + 1; + if(global_mode == 0) { + ls = 0; is = 1; ie = ch_n; + } else if(global_mode == 1) { + ls = 0; is = 1; + } else if(global_mode == 2) { + ie = ch_n; + } + if(is >= ie) return 0; + + for (l = ls, i = is; i < ie; 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; ///non-full-length + } + + if(i < ch_n) { + q[1] = ch_a[i].self_offset; t[1] = ch_a[i].offset; + } else { + q[1] = ql; ///non-full-length + } + + 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 + } + + if(mode == 1 || mode == 2) adjust_ext_offset(&(q[0]), &(q[1]), &(t[0]), &(t[1]), ql, tl, 0, mode); + + if(!umap) { + uq = qk; ut = tk; uc = ck; + } else { + wn0 = aux_o->w_list.n; + } + + if(is_scc_rev) { + est_err = cal_estimate_err_scc(scc, t[0], t[1], q[0], q[1], &tk, &qk, &ck); + } else { + est_err = cal_estimate_err_scc(scc, q[0], q[1], t[0], t[1], &qk, &tk, &ck); + } + is_done = hc_aln_exz_non_retrieve_adv_hc(z, qstr, tstr, q[0], q[1], t[0], t[1], mode, -1, exz, ql, tl, e_rate, MAX_SIN_L, MAX_SIN_E, FORCE_SIN_L, est_err, aux_o); + + if(!is_done) {///postprocess + est_hd = 0; + if((global_mode != -1)) {///no recal_non_retrieve_direct_scc anymore, directly get cigar from + est_ll = q[1] - q[0]; if(est_ll < (t[1] - t[0])) est_ll = t[1] - t[0]; + if((est_err > 16) && (est_err >= (est_ll*0.5))) {///bad window in anyway + est_hd = 1; + } + // fprintf(stderr, "[M::%s]\tq::[%ld,%ld)->qgap::%ld\tt::[%ld,%ld)->tgap::%ld\test_ll::%ld\test_err::%ld\test_hd::%u\n", __func__, q[0], q[1], q[1] - q[0], t[0], t[1], t[1] - t[0], est_ll, est_err, est_hd); + } + push_unmap_alnw_hold(aux_o, q[0], q[1]-1, t[0], t[1]-1, mode, est_hd, est_err); + if(global_mode == -1) { + umap = 1; + } + + if(uai[0] == INT64_MIN) { + uai[0] = l; + } + uai[1] = i; + } else if((umap) && (wn0 > 0)) { + if(is_ualn_win(aux_o->w_list.a[wn0-1])) { + assert(uai[0] != INT64_MIN && uai[1] != INT64_MIN); + assert(uai[1] > uai[0]); + if(uai[0] < 0) uai[0] = 0; + if(uai[1] >= ch_n) uai[1] = ch_n - 1; + recal_non_retrieve_direct_scc(z, cl, aux_o, wn0-1, qstr, tstr, exz, e_rate, ql, tl, h_khit, scc, is_scc_rev, FORCE_SIN_L, uq, ut, uc, + uai[0] + cidx, uai[1] + cidx); + ch_a = cl->list + cidx; + } + uq = qk; ut = tk; uc = ck; umap = 0; + uai[0] = uai[1] = INT64_MIN; wn0 = aux_o->w_list.n; + } + + l = i; + } + + wn0 = aux_o->w_list.n; //assert(wn0); + if(wn0 <= 0) return 0; + if((umap) && (wn0) && (is_ualn_win(aux_o->w_list.a[wn0-1]))) { + assert(uai[0] != INT64_MIN && uai[1] != INT64_MIN); + assert(uai[1] > uai[0]); + if(uai[0] < 0) uai[0] = 0; + if(uai[1] >= ch_n) uai[1] = ch_n - 1; + recal_non_retrieve_direct_scc(z, cl, aux_o, wn0-1, qstr, tstr, exz, e_rate, ql, tl, h_khit, scc, is_scc_rev, FORCE_SIN_L, uq, ut, uc, + uai[0] + cidx, uai[1] + cidx); + ch_a = cl->list + cidx; + } + return 1; +} + + void cigar_gen_by_chain_adv_local(overlap_region *z, Candidates_list *cl, ul_ov_t *ov, int64_t on, uint64_t wl, const ul_idx_t *uref, hpc_t *hpc_g, 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) { @@ -25448,8 +26233,8 @@ void est_rep_err_rate(overlap_region_alloc* ol, asg64_v *ix, kv_ul_ov_t *c_idx, ix->n = ixn0 - iin; } -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, uint64_t tcut, uint64_t site_sc, int64_t h0_w, asg32_v *b32, - int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t het_cov_a, int64_t hom_cov_a, int64_t n_hap, double hf_rate) +void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, asg16_v *qc0, 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, uint64_t tcut, uint64_t site_sc, int64_t h0_w, asg32_v *b32, + int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t het_cov_a, int64_t hom_cov_a, int64_t n_hap, double hf_rate, overlap_region *rchn, int64_t rref_len) { int64_t on = ol->length, k, i, zwn, q[2]; if(occ_thres + 1 < hap_cov_unmatch) occ_thres = hap_cov_unmatch-1; uint64_t m, l0, wi, wl0, si, ei, fi; overlap_region *z; ul_ov_t *cp; uint8_t *qhf = NULL, *qual = NULL; @@ -25595,6 +26380,7 @@ void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_all } **/ + int64_t wk = -1, ck = -1, qk = -1, tk = -1; SetSnpMatrix(hp, &(hp->nn_snp), &(ol->length), 0, NULL); srt_n = hp->length; z = (((std_bs)||(t8))?(ol->list):(NULL)); for (k = 1, m = /**ms =**/ i = t = 0; k <= srt_n; ++k) { @@ -25609,7 +26395,9 @@ void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_all **/ // if(m >= idx->n) then the region will be [ms, ql), odp = 0; // dbg_cp_cov(ol, c_idx, hp->list[i].site, odp); - t += push_info(hp, hp->list+i, k-i, hp->list+t, z, t8, 0/**odp**/, sc_wn, tcut); + // t += push_info(hp, hp->list+i, k-i, hp->list+t, z, t8, 0/**odp**/, sc_wn, tcut); + t += push_info_ssl(hp, hp->list+i, k-i, hp->list+t, z, t8, 0/**odp**/, sc_wn, tcut, qu->seq[hp->list[i].site], + rchn, qu->seq + qu->length, rref_len, &wk, &ck, &qk, &tk); i = k; } } @@ -25624,17 +26412,17 @@ void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_all if(!dp) { // generate_haplotypes_naive_advance(hap, overlap_list, NULL); - generate_haplotypes_naive_HiFi(hp, ol, 0.04, qu, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), hap_cov_match, hap_cov_unmatch); + generate_haplotypes_naive_HiFi(hp, ol, 0.04, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), hap_cov_match, hap_cov_unmatch); // generate_haplotypes_DP(hap, overlap_list, R_INF, g_read->length, force_repeat); // 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, b32); + gen_rphase_dp_adv(hp, ol, ((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); + if(!site_sc) generate_haplotypes_naive_HiFi_adv_hc(hp, ol, 0.04, ((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, ((std_bs)?(0):(1)), ((std_bs)?(0.05):(0)), ((std_bs)?(2):((uint64_t)-1)), 32, 2, hap_cov_match, hap_cov_unmatch); } if(lindel) { @@ -30982,6 +31770,129 @@ uint64_t gen_hc_fast_non_retrieve_cigar(overlap_region *z, Candidates_list *cl, return 1; } +uint64_t gen_hc_fast_non_retrieve_cigar_scc(overlap_region *z, Candidates_list *cl, int64_t wl, int64_t h_khit, 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, asg16_v* scc, int64_t *re) +{ + return_t_chain(z, cl); + + int64_t ch_idx = z->shared_seed, ch_n; + int64_t i, 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++){ + ;// fprintf(stderr, "[M::%s::]\tqpos::%u\ttpos::%u\n", __func__, cl->list[i].self_offset, cl->list[i].offset); + } + 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; + + hc_ovlp_base_non_retrieve_direct_scc(z, ch_a, ch_n, cl, h_khit, qstr, tstr, 0, 0, 0, exz, aux_o, e_rate, ql, tl, scc, 1, -1, rid); + + ///update z by aux_o + update_overlap_region(z, aux_o, ql, tl); + + int64_t wn = z->w_list.n; + for (i = tot_e = aln = 0; i < wn; 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; +} + +/** +uint8_t end_extend_l2g(uint16_t *sca, int64_t scn, int64_t iz_qs, int64_t iz_qe, int64_t iz_ts, int64_t iz_te, int64_t al_qs, int64_t al_qe, int64_t al_ts, int64_t al_te, + int64_t minl, int64_t ext_cc_bd, int64_t len_bd, uint8_t is_backward) +{ + int64_t ck, qk, tk, os, oe, wq[2], wt[2], oq[2], ot[2], qecut, tecut, qscut, tscut, l, err, qbd, tbd; asg16_v scc; scc.a = sca; scc.n = scc.m = scn; + uint16_t op, bq, bt; uint32_t cl; + + if(is_backward) { + ck = 0; qk = iz_qs; tk = iz_ts; qbd = iz_qs + ext_cc_bd; tbd = iz_ts + ext_cc_bd; err = 0; + qscut = MAX(iz_qs, al_qs); tscut = MAX(iz_ts, al_ts); + qecut = MIN(al_qe, iz_qe); if(qecut > qbd) qecut = qbd; if(qecut > len_bd) qecut = len_bd; + tecut = MIN(al_te, iz_qe); if(tecut > tbd) tecut = tbd; if(tecut > len_bd) tecut = len_bd; + + while ((ck < scn) && (qk < qecut) && (tk < tecut)) {//[s, e) + wq[0] = qk; wt[0] = tk; + ck = pop_trace_bp_f(&scc, ck, &op, &bq, &bt, &cl); + if(op != 2) qk += cl; + if(op != 3) tk += cl; + wq[1] = qk; wt[1] = tk; + + if(op != 0) { + err += cl; + continue;///exact match + } + + if((wq[1] <= qscut) || (wt[1] <= tscut)) continue; + + + oq[0] = MAX(qscut, wq[0]); oq[1] = MIN(qecut, wq[1]); + l = oq[1] - oq[0]; + if(l < minl) continue; + + ot[0] = wt[0] + oq[0] - wq[0]; ot[1] = wt[0] + oq[1] - wq[0]; + os = MAX(tscut, ot[0]);oe = MIN(tecut, ot[1]); + l = oe - os; + if(l < minl) continue; + + + qrch = wq[0] - iz_qs; trch = wt[0] - iz_ts; + + + + } + } +} + +uint64_t l2g_non_retrieve_cigar(overlap_region *z, Candidates_list *cl, int64_t h_khit, char *qstr, int64_t ql, char *tstr, int64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate) +{ + if(z->w_list.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; + + // ((b16).a = (des->w_list.c).a, (b16).n = (des->w_list.c).n, (b16).m = (des->w_list.c).m) + + + if((z->w_list.a[0].x_start > 0) || (z->w_list.a[0].y_start > 0)) { + end_extend_l2g(z->w_list.c.a + z->w_list.a[0].cidx, z->w_list.a[0].clen - z->w_list.a[0].cidx, + z->w_list.a[0].x_start, z->w_list.a[0].x_end + 1, z->w_list.a[0].y_start, z->w_list.a[0].y_end + 1, 0, ql, 0, tl, 16, 768, 2048, 1); + } + int64_t ch_idx = z->shared_seed, ch_n, i; + k_mer_hit *ch_a = cl->list + ch_idx; + 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; +} +**/ + void append_cigar(window_list *idx, window_list_alloc *res, uint16_t c, uint32_t l) { if(l <= 0) return; @@ -36033,6 +36944,8 @@ double fcov_rat, uint64_t ch_occ, uint64_t ch_sc) uint64_t spn0 = sp->n, *wsrt = NULL, *wcut = NULL, focv = 0, fcov_0 = 0, i, k, m, wcut_n = 0, wsrt_n = 0, wma_n = 0, bs, tot_b = 0, t_cov0, z_cov, tz; int32_t s[4]; uint32_t scn[4]; scn[0] = scn[1] = scn[2] = scn[3] = 0; s[0] = s[1] = s[2] = s[3] = 0; sp->n += ol->length; kv_resize(uint64_t, (*sp), sp->n); wsrt = sp->a + spn0; + + // fprintf(stderr, "[M::%s::]\tol->length::%lu\tmax_n_chain::%lu\n", __func__, ol->length, max_n_chain); if(fc) { permute_in_place_ss(ol->list, ol->length, wsrt, osc, ocn); wcut = infer_chn_bar_0(ol, max_n_chain, ocn, osc, ql, ocw, s, sp, &wcut_n, &focv); @@ -36068,6 +36981,9 @@ double fcov_rat, uint64_t ch_occ, uint64_t ch_sc) ff = 0; } } + + // fprintf(stderr, "-0-[M::%s::]\ttid::%u(%u)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\twf::%u\tosc::%u\tocn::%u\tff::%u\n", __func__, + // z->y_id, z->y_pos_strand, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1, wf, osc[i], ocn[i], ff); if((!ff) && (z->is_match == 0)) {///fitered out due to coverage wsrt[wsrt_n++] = (((uint64_t)osc[i])<<32)|(i); @@ -36089,6 +37005,9 @@ double fcov_rat, uint64_t ch_occ, uint64_t ch_sc) wsrt_n++; wma_n++; tot_b += z->x_pos_e + 1 - z->x_pos_s; z->non_homopolymer_errors = UINT32_MAX - 1;///primary chain that needs to be verfied + + // fprintf(stderr, "-1-[M::%s::]\ttid::%u(%u)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\twf::%u\tosc::%u\tocn::%u\tff::%u\tferr::%u\n", __func__, + // z->y_id, z->y_pos_strand, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1, wf, osc[i], ocn[i], ff, z->non_homopolymer_errors); } if(fc && wsrt_n > wma_n) { @@ -37538,6 +38457,8 @@ void gen_hc_r_alin_adv_adp_smp_1(gen_hc_aln_t *ez, uint8_t set_match) for (i = nol_1 = 0; i < wsrt_n; i++) { z = &(ez->ol->list[(uint32_t)wsrt[i]]); + // fprintf(stderr, "-a-[M::%s::]\ttid::%u(%u)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\n", __func__, + // z->y_id, z->y_pos_strand, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1); if(z->is_match == 0) { z->w_list.n = 0; z->align_length = 0; } @@ -37556,6 +38477,11 @@ void gen_hc_r_alin_adv_adp_smp_1(gen_hc_aln_t *ez, uint8_t set_match) ocn = ez->v32->a; osc = ez->v32->a + ez->ol->length; for (i = 0; i < wsrt_n; i++) { + // z = &(ez->ol->list[(uint32_t)wsrt[i]]); + // fprintf(stderr, "-0-[M::%s::]\ttid::%u(%u)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\tpass::%u\n", __func__, + // z->y_id, z->y_pos_strand, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1, + // pass_qovlp(z->x_pos_e+1-z->x_pos_s, z->align_length, OVERLAP_THRESHOLD_HIFI_FILTER)); + if((wsrt[i]>>32) == UINT32_MAX) { continue; } @@ -37568,12 +38494,20 @@ void gen_hc_r_alin_adv_adp_smp_1(gen_hc_aln_t *ez, uint8_t set_match) } e_max = err * 1.5; + // fprintf(stderr, "-1-[M::%s::]\ttid::%u(%u)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\tpass::%u\n", __func__, + // z->y_id, z->y_pos_strand, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1, + // pass_qovlp(z->x_pos_e+1-z->x_pos_s, z->align_length, OVERLAP_THRESHOLD_HIFI_FILTER)); + if(z->is_match == 0) { if((!pass_qovlp(z->x_pos_e+1-z->x_pos_s, z->align_length, OVERLAP_THRESHOLD_HIFI_FILTER)) || (!gen_hc_r_alin_flt_1_smp(z, ez->cl, ez->rref, ez->qu, ez->tu, ez->exz, ez->aux_o, err, e_max, err, w, ql, ez->rid, ez->khit, chem_drop, align_gap_rate, align_gap_max, NULL/**ez->hpz->a**/, ez->buf, 1, NULL/**&tot_b**/))) { continue; } } + + // fprintf(stderr, "-2-[M::%s::]\ttid::%u(%u)\t%.*s\tq::[%u,%u)\tt::[%u,%u)\terr::%u\n", __func__, + // z->y_id, z->y_pos_strand, (int)Get_NAME_LENGTH(R_INF, z->y_id), Get_NAME(R_INF, z->y_id), z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1, + // z->non_homopolymer_errors); // if((z->is_match == 0) && (!gen_hc_r_alin_flt_1_smp(z, ez->cl, ez->rref, ez->qu, ez->tu, ez->exz, ez->aux_o, err, e_max, err, w, ql, ez->rid, ez->khit, chem_drop, align_gap_rate, align_gap_max, NULL/**ez->hpz->a**/, ez->buf, 0, &tot_b))) { // continue; // } @@ -37833,6 +38767,22 @@ uint64_t gen_hc_r_alin_re(overlap_region* z, Candidates_list *cl, char* qstr, ui return 1; } +uint64_t gen_hc_r_alin_self(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, asg16_v* scc) +{ + int64_t re; + + z->shared_seed = z->non_homopolymer_errors;///for index + + if(!gen_hc_fast_non_retrieve_cigar_scc(z, cl, wl, khit, qstr, ql, tstr, tl, exz, aux_o, e_rate, rid, khit, scc, &re)) return 0; + + reassign_gaps(z, aux_o, qstr, ql, NULL, tstr, tl, NULL, NULL, buf, exz); + + z->is_match = 1; ///z->non_homopolymer_errors = re; + z->strong = z->without_large_indel = 0; + + return 1; +} + /** void ul_raw_lalign_adv(overlap_region_alloc* ol, Candidates_list *cl, const ul_idx_t *uref, All_reads *rdb, const ug_opt_t *uopt, char *qstr, uint64_t ql, UC_Read* qu, UC_Read* tu, Correct_dumy* dumy, bit_extz_t *exz, haplotype_evdience_alloc* hap, diff --git a/Correct.h b/Correct.h index da081fe..3889a72 100644 --- a/Correct.h +++ b/Correct.h @@ -1456,8 +1456,9 @@ void gen_hc_r_alin_adv(gen_hc_aln_t *ez); uint64_t 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 sec_aln_win, uint64_t sec_aln_cov, double sec_aln_err_rate, double sec_aln_max, asg64_v *kp, uint8_t *hpf); void gen_hc_r_alin_nec_adv(gen_hc_aln_t *ez); 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, uint64_t tcut, uint64_t site_sc, int64_t h0_w, asg32_v *b32, - int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t het_cov_a, int64_t hom_cov_a, int64_t n_hap, double hf_rate); +uint64_t gen_hc_r_alin_self(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, asg16_v* scc); +void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, asg16_v *qc0, 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, uint64_t tcut, uint64_t site_sc, int64_t h0_w, asg32_v *b32, + int64_t hap_cov_match, int64_t hap_cov_unmatch, int64_t het_cov_a, int64_t hom_cov_a, int64_t n_hap, double hf_rate, overlap_region *rchn, int64_t rref_len); void set_exact_exz(bit_extz_t *exz, int64_t qs, int64_t qe, int64_t ts, int64_t te); void push_alnw(overlap_region *aux_o, bit_extz_t *exz); void cal_exz_global(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez); diff --git a/Hash_Table.cpp b/Hash_Table.cpp index c44fcfe..9062998 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -1540,6 +1540,63 @@ inline int32_t comput_sc_ch_ec(const k_mer_hit *ai, const k_mer_hit *aj, double return sc; } +inline int32_t comput_sc_ch_ec_global(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, double chn_pen_gap, double chn_pen_skip, int64_t sl, int64_t ol, uint8_t no_adj) +{ + ///ai is the suffix of aj + int32_t dq, dr, dd, dg, q_span, sc; double dg_of; + dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset); + dr = (int64_t)(ai->offset) - (int64_t)(aj->offset); + if((dq < 0) || (dr < 0)) return INT32_MIN; + if((no_adj) && ((dq == 0) || (dr == 0))) return INT32_MIN; + + dd = dr > dq? dr - dq : dq - dr;//gap + if((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, sl, ol))) return INT32_MIN; + + dg = dr < dq? dr : dq;//len + q_span = ai->cnt&(0xffu); + sc = q_span < dg? q_span : dg; + sc = normal_w(sc, ((int32_t)(ai->cnt>>8))); + if (dd || (dg > q_span && dg > 0)) { + double lin_pen, a_pen; + lin_pen = (chn_pen_gap*(double)dd); + dg_of = (dg>0)?((double)dg):(0.333333); + a_pen = ((double)(sc))*((((double)dd)/dg_of)/bw_rate); + ///for long gap + // if(lin_pen > a_pen) lin_pen = a_pen; + if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen)); + else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen)); + lin_pen += (chn_pen_skip*dg_of); + sc -= (int32_t)lin_pen; + } + return sc; +} + +inline int32_t comput_sc_ff_adv(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, double chn_pen_gap, double chn_pen_skip, int64_t sl, int64_t ol) +{ + ///ai is the suffix of aj + int32_t dq, dr, dd, dg, q_span, sc; + dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset); + if(dq < 0) return INT32_MIN; + dr = (int64_t)(ai->offset) - (int64_t)(aj->offset); + if(dr < 0) return INT32_MIN; + dd = dr > dq? dr - dq : dq - dr;//gap + // if((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, sl, ol))) return INT32_MIN; + dg = dr < dq? dr : dq;//len + if(dg <= 0) return INT32_MIN; + q_span = ai->cnt&(0xffu); + sc = q_span < dg? q_span : dg; + sc = normal_w(sc, ((int32_t)(ai->cnt>>8))); + if (dd || (dg > q_span && dg > 0)) { + double lin_pen, a_pen; + lin_pen = (chn_pen_gap*(double)dd); + a_pen = ((double)(sc))*((((double)dd)/((double)dg))/bw_rate); + if(lin_pen > a_pen) lin_pen = a_pen; + lin_pen += (chn_pen_skip*(double)dg); + sc -= (int32_t)lin_pen; + } + return sc; +} + inline int32_t comput_sc_ff(const k_mer_hit *ai, const k_mer_hit *aj, double bw_rate, double chn_pen_gap, double chn_pen_skip, int64_t sl, int64_t ol) { ///ai is the suffix of aj @@ -2283,6 +2340,315 @@ uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n, } +/** +void quick_ck_lchain_global(k_mer_hit* a, int64_t a_n, int64_t xl, int64_t yl, double chn_pen_gap, double chn_pen_skip, double bw_rate, +int64_t *p, int64_t *t, int32_t *f, int32_t *ii, int32_t *gf, int64_t *gp, int64_t *si, int64_t *ei) +{ + *si = 0; *ei = a_n; gf[0] = gf[1] = gp[0] = gp[1] = INT32_MIN; + if((a_n <= 0) || (xl <= 0) || (yl <= 0)) return; + int64_t l, k, is_srt = 1, z; k_mer_hit *ai, *aj, ft, fz; + int64_t dq, dr, dd, dg, q_span, sc, csc, ddt; uint8_t ff; double lin_pen, a_pen, dg_of; + ft.cnt = fz.cnt = 0xFFFFFF00u; ft.readID = fz.readID = a[0].readID; + ft.offset = ft.self_offset = 0; fz.offset = yl-1; fz.self_offset = xl-1; + + for (k = 1, l = 0; k <= a_n; k++) { + if(k == a_n || a[k].strand != a[l].strand) { + t[k-1] = 0; ii[k-1] = 0; + + if(is_srt) { + ddt = 0; ff = 0; p[l] = f[l] = INT32_MIN; + + ft.strand = a[l].strand; + aj = &ft; z = l; ai = &a[z]; + dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset); + dr = (int64_t)(ai->offset) - (int64_t)(aj->offset); + dd = dr > dq? dr - dq : dq - dr;//gap + if((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, xl, yl))) ff = 1; + if(!ff) { + dg = dr < dq? dr : dq;//len + q_span = ai->cnt&(0xffu); + sc = q_span < dg? q_span : dg; + sc = normal_w(sc, ((int32_t)(ai->cnt>>8))); + if (dd || (dg > q_span && dg > 0)) { + lin_pen = (chn_pen_gap*(double)dd); + dg_of = ((dg>0)?(dg):(0.333333)); + a_pen = ((double)(sc))*((((double)dd)/dg_of)/bw_rate); + ///for long gap + // if(lin_pen > a_pen) lin_pen = a_pen; + if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen)); + else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen)); + lin_pen += (chn_pen_skip*dg_of); + sc -= (int32_t)lin_pen; + } + + csc = a[z].cnt&(0xffu); if(sc < csc) ff = 1; + if(!ff) { + p[z] = -1; f[z] = sc; ddt += dd; + } + } + + if(!ff) { + for (z = l + 1; z < k; z++) { + ///roughly same to comput_sc_ch(&a[z], &a[z-1]) + ai = &a[z]; aj = &a[z-1]; + dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset); + if(dq <= 0) break; + dr = (int64_t)(ai->offset) - (int64_t)(aj->offset); + if(dr <= 0) break; + dd = dr > dq? dr - dq : dq - dr;//gap + + if((dd > 16) && (dd > cal_bw(&(a[z]), &(a[z-1]), bw_rate, xl, yl))) break; + dg = dr < dq? dr : dq;//len + q_span = ai->cnt&(0xffu); + sc = q_span < dg? q_span : dg; + sc = normal_w(sc, ((int32_t)(ai->cnt>>8))); + if (dd || (dg > q_span && dg > 0)) { + lin_pen = (chn_pen_gap*(double)dd); + dg_of = ((dg>0)?(dg):(0.333333)); + a_pen = ((double)(sc))*((((double)dd)/dg_of)/bw_rate); + ///for long gap + // if(lin_pen > a_pen) lin_pen = a_pen; + if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen)); + else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen)); + lin_pen += (chn_pen_skip*dg_of); + sc -= (int32_t)lin_pen; + } + + sc += f[z-1]; csc = a[z].cnt&(0xffu); if(sc < csc) break; + p[z] = z - 1; f[z] = sc; ddt += dd; + } + if(z < k) ff = 1; + } + + if(!ff) { + fz.strand = a[l].strand; + ai = &fz; aj = &(a[k-1]); + dq = (int64_t)(ai->self_offset) - (int64_t)(aj->self_offset); + dr = (int64_t)(ai->offset) - (int64_t)(aj->offset); + dd = dr > dq? dr - dq : dq - dr;//gap + if(((dd > 16) && (dd > cal_bw(ai, aj, bw_rate, xl, yl)))|| + (((ddt + dd) > 16) && ((ddt + dd) > cal_bw(&fz, &ft, bw_rate, xl, yl)))) { + ff = 1; + } else { + dg = dr < dq? dr : dq;//len + sc = q_span = 0; + if (dd || (dg > q_span && dg > 0)) { + lin_pen = (chn_pen_gap*(double)dd); + dg_of = ((dg>0)?(dg):(0.333333)); + a_pen = ((double)(sc))*((((double)dd)/dg_of)/bw_rate); + ///for long gap + // if(lin_pen > a_pen) lin_pen = a_pen; + if(dd < 4) lin_pen = ((lin_pen > a_pen)?(a_pen):(lin_pen)); + else lin_pen = ((lin_pen < a_pen)?(a_pen):(lin_pen)); + lin_pen += (chn_pen_skip*dg_of); + sc -= (int32_t)lin_pen; + } + sc += f[k-1]; ///csc = a[z].cnt&(0xffu); if(sc < csc) break; + gp[fz.strand] = k-1; gf[fz.strand] = sc; ///ddt += dd; + if((*ei) > k) { + (*si) = k; + } else { + (*ei) = l; + } + } + } + } + l = k; is_srt = 1; + } else { + if((a[k].self_offset <= a[k-1].self_offset) || (a[k].offset <= a[k-1].offset)) is_srt = 0; + t[k-1] = 0; ii[k-1] = 0; + } + } +} + + +uint64_t lchain_qdp_global_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64_t des_idx, + Chain_Data* dp, overlap_region_alloc* res, int64_t max_skip, int64_t max_iter, + int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate, + uint32_t xid, int64_t xl, int64_t yl, int64_t quick_check, uint32_t apend_be, + int64_t gen_cigar, int64_t khit_n) +{ + if(a_n <= 0) return 0; + int64_t *p, *t, *gp, max_f, n_skip, st, max_j, end_j, sc, max_ii, ovl, min_sc, ch_n, si, ei; + int32_t *f, max, tmp, *ii, *gf; int64_t i, k, j, cL = 0; k_mer_hit* a; k_mer_hit* des; k_mer_hit *swap, ft; overlap_region *z; + resize_Chain_Data(dp, a_n + 2, NULL); ch_n = 1; + t = dp->tmp; f = dp->score; p = dp->pre; ii = dp->occ; gp = p + a_n; gf = f + a_n; gf[0] = gf[1] = gp[0] = gp[1] = INT32_MIN; + + a = cl->list + a_idx; des = cl->list + des_idx; + if(quick_check) { + quick_ck_lchain_global(a, a_n, xl, yl, chn_pen_gap, chn_pen_skip, bw_rate, p, t, f, ii, gf, gp, &si, &ei); + } else { + si = 0; ei = a_n; memset(t, 0, (a_n*sizeof((*t)))); + } + + ft.cnt = 0xFFFFFF00u; ft.readID = a[0].readID; ft.strand = 0; ft.offset = ft.self_offset = 0; + for (i = st = si, max_ii = -1; i < ei; ++i) { + ///max_f = a[i].cnt&(0xffu); + ft.strand = a[i].strand; + + + n_skip = 0; max_j = end_j = -1; + if ((i-st) > max_iter) st = i-max_iter; + while (a[i].strand != a[st].strand) ++st; + + for (j = i - 1; j >= st; --j) { + sc = comput_sc_ch_ec_global(&a[i], &a[j], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl); + if (sc == INT32_MIN) continue; + sc += f[j]; + if (sc > max_f) { + max_f = sc, max_j = j; + if (n_skip > 0) --n_skip; + } else if (t[j] == (int32_t)i) { + if (++n_skip > max_skip) + break; + } + if (p[j] >= 0) t[p[j]] = i; + } + end_j = j; + + if ((max_ii<0) || (a[i].self_offset>a[max_ii].self_offset+max_dis) || (a[i].strand!=a[max_ii].strand)) { + max = INT32_MIN; max_ii = -1; + for (j=i-1; (j>=st) && (a[i].self_offset<=max_dis+a[j].self_offset)&&(a[i].strand==a[j].strand); --j) { + if (max < f[j]) { + max = f[j], max_ii = j; + } + } + } + + if ((max_ii >= 0) && (max_ii < end_j) && (a[i].strand == a[max_ii].strand)) {///just have a try with a[i]<->a[max_ii] + tmp = comput_sc_ch_ec(&a[i], &a[max_ii], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl); + if (tmp != INT32_MIN && max_f < tmp + f[max_ii]) + max_f = tmp + f[max_ii], max_j = max_ii; + } + f[i] = max_f; p[i] = max_j; + if ((max_ii < 0) || ((a[i].self_offset<=max_dis+a[max_ii].self_offset)&&(a[i].strand==a[max_ii].strand)&&(f[max_ii]= msc) { + ovl = get_chainLen(a[i].self_offset, a[i].self_offset, xl, a[i].offset, a[i].offset, yl); + if(f[i] > msc || ovl < movl) { + msc = f[i]; msc_i = i; movl = ovl; + } + } + if(f[i] < plus) plus = f[i]; + ii[i] = 0;///for mcopy, not here + // if(a_n && (a[0].readID == 27105 || a[0].readID == 7603)) {///r833 + // fprintf(stderr, "i::%ld[M::%s::rid->%u::%c] q::%u, t::%u, st::%ld, max_ii::%ld, f[i]::%d, p[i]::%ld, msc_i::%ld, msc::%ld, movl::%ld\n", + // i, __func__, a[i].readID, "+-"[a[i].strand], + // a[i].self_offset, a[i].offset, st, max_ii, f[i], p[i], msc_i, msc, movl); + // } + } + + for (i = msc_i, cL = 0; i >= 0; i = p[i]) { ii[i] = 1; t[cL++] = i;}///label the best chain + + if(mcopy_num > 1) { + // if(a[0].readID == 4412344) { + // fprintf(stderr, "[M::%s::] msc::%ld, cL::%ld\n", __func__, msc, cL); + // } + if(cL >= mcopy_khit_cutoff) {///if there are too few k-mers, disable mcopy + msc -= plus; min_sc = msc*mcopy_rate; ii[msc_i] = 0; + for (i = ch_n = 0; i < a_n; ++i) {///make all f[] positive + f[i] -= plus; if(i >= ch_n) t[i] = 0; + if((!(ii[i])) && (f[i] >= min_sc)) {///!(ii[i]): skip the best chain + t[ch_n] = ((uint64_t)f[i])<<32; t[ch_n] += (i<<1); ch_n++; + } + } + // if(a[0].readID == 4412344) { + // fprintf(stderr, "[M::%s::] msc::%ld, min_sc::%ld, cL::%ld, ch_n::%ld, mcopy_num::%ld\n", __func__, msc, min_sc, cL, ch_n, mcopy_num); + // } + if(ch_n > 1) { + int64_t n_v, n_v0, ni, n_u, n_u0 = res->length; + radix_sort_hc64i(t, t + ch_n); + for (k = ch_n-1, n_v = n_u = 0; k >= 0 && n_u < mcopy_num; --k) { + n_v0 = n_v; + for (i = ((uint32_t)t[k])>>1; i >= 0 && (t[i]&1) == 0; ) { + ii[n_v++] = i; t[i] |= 1; i = p[i]; + } + if(n_v0 == n_v) continue; + sc = (i<0?(t[k]>>32):((t[k]>>32)-f[i])); + // if(a[0].readID == 4412344) { + // fprintf(stderr, "+[M::%s::] sc::%ld, n_a::%ld\n", __func__, sc, n_v-n_v0); + // } + if(sc >= min_sc) { + kv_pushp_ol(overlap_region, (*res), &z); + push_ovlp_chain_qgen(z, xid, xl, yl, sc+plus, &(a[ii[n_v-1]]), &(a[ii[n_v0]])); + // if(a[0].readID == 4412344) { + // fprintf(stderr, "-[M::%s::] sc::%ld, n_a::%ld, q::[%u,%u), t::[%u,%u), %c\n", __func__, sc, n_v-n_v0, z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1, "+-"[z->y_pos_strand]); + // } + ///mcopy_khit_cutoff <= 1: disable the mcopy_khit_cutoff filtering, for the realignment + // if((mcopy_khit_cutoff <= 1) || ((z->x_pos_e+1-z->x_pos_s) <= (movl<<2))) { + if((!n_u) || (n_v - n_v0 > 1)) { + z->align_length = n_v-n_v0; z->x_id = n_v0; + n_u++; + } else {///non-best is tiny + res->length--; n_v = n_v0; + } + } else { + n_v = n_v0; + } + } + + // if(n_u > 1) ks_introsort_or_sss(n_u, res->list + n_u0); + // res->length = n_u0 + filter_non_ovlp_xchains(res->list + n_u0, n_u, &n_v); + n_u = res->length; + if(n_u > n_u0 + 1) { + kv_resize_cl(k_mer_hit, (*cl), (n_v+cl->length)); + a = cl->list + a_idx; des = cl->list + des_idx; swap = cl->list + cl->length; + for (k = n_u0, i = n_v0 = n_v = 0; k < n_u; k++) { + z = &(res->list[k]); + z->non_homopolymer_errors = des_idx + i; + n_v0 = z->x_id; ni = z->align_length; + for (j = 0; j < ni; j++, i++) { + ///k0 + (ni - j - 1) + swap[i] = a[ii[n_v0 + (ni- j - 1)]]; + swap[i].readID = k; + } + z->x_id = xid; + if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, swap+i-ni, ni); + if(!khit_n) z->align_length = 0; + } + memcpy(des, swap, i*sizeof((*swap))); //assert(i == ch_n); + + // fprintf(stderr, "[M::%s::msc->%ld] msc_k_hits::%u, cL::%ld, min_sc::%ld, best_sc::%ld, n_u0_sc::%d, mcopy_rate::%f, # chains::%ld\n", + // __func__, msc, res->list[n_u0].align_length, cL, min_sc, msc+plus, res->list[n_u0].shared_seed, + // mcopy_rate, n_u-n_u0); + } else if(n_u == n_u0 + 1) { + z = &(res->list[n_u0]); k = n_u0; i = 0; + z->non_homopolymer_errors = des_idx + i; + n_v0 = z->x_id; ni = z->align_length; + for (j = 0; j < ni; j++, i++) { + ///k0 + (ni - j - 1) + des[i] = a[ii[n_v0 + (ni- j - 1)]]; + des[i].readID = k; + } + z->x_id = xid; + if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, des+i-ni, ni); + if(!khit_n) z->align_length = 0; + } + return i; + } else { + msc += plus; i = msc_i; cL = 0; + while (i >= 0) {t[cL++] = i; i = p[i];} + } + } + } + + + + ///a[] has been sorted by self_offset + // i = msc_i; cL = 0; + // while (i >= 0) {t[cL++] = i; i = p[i];} + kv_pushp_ol(overlap_region, (*res), &z); + push_ovlp_chain_qgen(z, xid, xl, yl, msc, &(a[t[cL-1]]), &(a[t[0]])); + for (i = 0; i < cL; i++) {des[i] = a[t[cL-i-1]]; des[i].readID = res->length-1;} + z->non_homopolymer_errors = des_idx; + if(gen_cigar) gen_fake_cigar(&(z->f_cigar), z, apend_be, des, cL); + if(khit_n) z->align_length = cL; + return cL; +} +**/ + + #define rev_khit(an, xl, yl) do { \ (an).self_offset = (xl)-1-((an).self_offset+1-((an).cnt&((uint32_t)(0xffu)))); \ (an).offset = (yl)-1-((an).offset+1-((an).cnt&((uint32_t)(0xffu))));\ @@ -2453,6 +2819,121 @@ uint64_t lchain_qdp_fix(k_mer_hit* a, int64_t a_n, Chain_Data* dp, int64_t max_s } +uint64_t lchain_qdp_fix_adv(k_mer_hit *a, int64_t a_n, Chain_Data* dp, int64_t max_skip, + int64_t max_iter, int64_t max_dis, double chn_pen_gap, double chn_pen_skip, + double bw_rate, int64_t xl, int64_t yl, int64_t quick_check, + int64_t left_fix, int64_t right_fix, k_mer_hit *res) +{ + if(a_n <= 0) return 0; + + int64_t *p, *t, max_f, n_skip, st, max_j, end_j, sc, msc, msc_i, bw, max_ii, ovl, movl; + int32_t *f, max, tmp; int64_t i, j, ret, cL = 0; + resize_Chain_Data(dp, a_n, NULL); + t = dp->tmp; f = dp->score; p = dp->pre; + bw = ((xl < yl)?xl:yl); bw *= bw_rate; + msc = msc_i = -1; movl = INT32_MAX; + + if(quick_check) { + ret = lchain_qcheck(a, a_n, dp, bw_rate); + if (ret > 0) { + a_n = ret; msc_i = a_n-1; msc = f[msc_i]; + goto skip_ldp; + } + } + + memset(t, 0, (a_n*sizeof((*t)))); + for (i = st = 0, max_ii = -1; i < a_n; ++i) { + max_f = a[i].cnt&(0xffu); if(left_fix && i > 0) max_f = INT32_MIN; + n_skip = 0; max_j = end_j = -1; + if ((i-st) > max_iter) st = i-max_iter; + + for (j = i - 1; j >= 0; --j) { + if(left_fix && f[j] == INT32_MIN)continue; + sc = comput_sc_ff_adv(&a[i], &a[j], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl); + if (sc == INT32_MIN) continue; + sc += f[j]; + if (sc > max_f) { + max_f = sc, max_j = j; + if (n_skip > 0) --n_skip; + } else if (t[j] == (int32_t)i) { + if ((++n_skip) > max_skip) { + if((max_j != -1) || (left_fix == 0)) break; + } + } + if (p[j] >= 0) t[p[j]] = i; + ///put it here will allow at least one prefix no matter max_dis + ///this is special for gap filling, not for chaining + if (a[i].self_offset > (max_dis + a[j].self_offset)) { + if((max_j != -1)) break; + } + if (j < st) { + if((max_j != -1) || (left_fix == 0)) break; + } + } + end_j = j; + + if (max_ii < 0 || ((int64_t)a[i].self_offset) - ((int64_t)a[max_ii].self_offset) > max_dis) { + max = INT32_MIN; max_ii = -1; + for (j = i - 1; (j >= st) && ((((int64_t)a[i].self_offset)-((int64_t)a[j].self_offset))<=max_dis); --j) { + if ((f[j] != INT32_MIN) && (max < f[j])) { + max = f[j], max_ii = j; + } + } + } + + if ((max_ii >= 0) && (max_ii < end_j) && (f[max_ii] != INT32_MIN)) {///just have a try with a[i]<->a[max_ii] + tmp = comput_sc_ff_adv(&a[i], &a[max_ii], bw_rate, chn_pen_gap, chn_pen_skip, xl, yl); + if (tmp != INT32_MIN && max_f < tmp + f[max_ii]) + max_f = tmp + f[max_ii], max_j = max_ii; + } + + p[i] = max_j; f[i] = max_f; + + if ((max_ii < 0) || (((((int64_t)a[i].self_offset)-((int64_t)a[max_ii].self_offset))<=max_dis) && (f[max_ii]= msc) { + ovl = get_chainLen(a[i].self_offset, a[i].self_offset, xl, a[i].offset, a[i].offset, yl); + if(f[i] > msc || ovl < movl) { + msc = f[i]; msc_i = i; movl = ovl; + } + } + } + + skip_ldp: + if(right_fix && f[a_n-1] == INT32_MIN) return 0; + if(right_fix) msc_i = a_n-1; + ///a[] has been sorted by self_offset + i = msc_i; cL = 0; + while (i >= 0) { + t[cL++] = i; msc_i = i; i = p[i]; + } + + n_skip = cL>>1; + for (i = 0; i < n_skip; i++) { + msc_i = t[i]; t[i] = t[cL-i-1]; t[cL-i-1] = msc_i; + } + + if((cL > 0) && (right_fix) && (t[cL-1] != (a_n-1))) { + cL = 0; + } + + if((cL > 0) && (left_fix) && (t[0] != 0)) { + cL = 0; + } + + if(cL > 0 && res) { + for (i = 0; i < cL; i++) { + res[i] = a[t[i]]; + } + } + + return cL; +} + + + uint64_t lchain_refine(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, int64_t max_skip, int64_t max_iter, int64_t max_dis, int64_t long_gap) { diff --git a/Hash_Table.h b/Hash_Table.h index 9004fda..7c692be 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -234,6 +234,10 @@ uint64_t lchain_qdp_fix(k_mer_hit* a, int64_t a_n, Chain_Data* dp, int64_t max_s int64_t max_iter, int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate, int64_t xl, int64_t yl, int64_t quick_check, int64_t left_fix, int64_t right_fix); +uint64_t lchain_qdp_fix_adv(k_mer_hit* a, int64_t a_n, Chain_Data* dp, int64_t max_skip, + int64_t max_iter, int64_t max_dis, double chn_pen_gap, double chn_pen_skip, + double bw_rate, int64_t xl, int64_t yl, int64_t quick_check, + int64_t left_fix, int64_t right_fix, k_mer_hit *res); uint64_t lchain_simple(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, int64_t max_skip, int64_t max_iter); uint64_t lchain_simple0(k_mer_hit* a, int64_t a_n, k_mer_hit* des, Chain_Data* dp, int64_t max_skip, int64_t max_iter); @@ -253,6 +257,12 @@ uint64_t lchain_qdp_mcopy_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64_t gen_cigar, int64_t enable_mcopy, double mcopy_rate, int64_t mcopy_khit_cutoff, int64_t khit_n); +uint64_t lchain_qdp_global_fast(Candidates_list *cl, int64_t a_idx, int64_t a_n, int64_t des_idx, + Chain_Data* dp, overlap_region_alloc* res, int64_t max_skip, int64_t max_iter, + int64_t max_dis, double chn_pen_gap, double chn_pen_skip, double bw_rate, + uint32_t xid, int64_t xl, int64_t yl, int64_t quick_check, uint32_t apend_be, + int64_t gen_cigar, int64_t khit_n); + #define kv_pushp_ol(type, v, p) do { \ if ((v).length == (v).size) { \ (v).list = (type*)realloc((v).list, sizeof(type)*((v).size?((v).size<<1):(2))); \ diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h index d53b2a9..015473d 100644 --- a/Levenshtein_distance.h +++ b/Levenshtein_distance.h @@ -712,19 +712,30 @@ inline uint32_t pop_trace_bp_f(asg16_v *res, uint32_t i, uint16_t *c, uint16_t * inline int64_t pop_trace_bp_rev_f(asg16_v *res, int64_t i, uint16_t *c, uint16_t *bq, uint16_t *bt, uint32_t *len) { - (*c) = (res->a[i]>>14); (*bq) = (*bt) = (uint16_t)-1; - if((*c) == 2 || (*c) == 3) { - (*bt) = ((res->a[i]>>12)&3); - (*len) = (res->a[i]&(0xfff)); - } else if((*c) == 1) { - (*bt) = ((res->a[i]>>12)&3); - (*bq) = ((res->a[i]>>10)&3); - (*len) = (res->a[i]&(0x3ff)); - } else { - (*len) = (res->a[i]&(0x3fff)); + uint32_t sl = 1; + (*c) = (*bq) = (*bt) = (uint16_t)-1; (*len) = 0; + if(i >= ((int64_t)res->n)) { + i = ((int64_t)res->n) - 1; sl = 0; + } + if(i >= 0) { + (*c) = (res->a[i]>>14); + (*bq) = (*bt) = (uint16_t)-1; + if((*c) == 2 || (*c) == 3) { + (*bt) = ((res->a[i]>>12)&3); + (*len) = (res->a[i]&(0xfff)); + } else if((*c) == 1) { + (*bt) = ((res->a[i]>>12)&3); + (*bq) = ((res->a[i]>>10)&3); + (*len) = (res->a[i]&(0x3ff)); + } else { + (*len) = (res->a[i]&(0x3fff)); + } + } + if(sl == 0) { + (*len) = 0; i = res->n; } - uint32_t sl; uint16_t sbq, sbt; + uint16_t sbq, sbt; for (i--; (i >= 0) && ((*c) == (res->a[i]>>14)); i--) { sbq = sbt = (uint16_t)-1; if((*c) == 2 || (*c) == 3) { diff --git a/Overlaps.cpp b/Overlaps.cpp index 3f5f2e6..2799503 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -23747,7 +23747,7 @@ void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sou sprintf(gfa_name, "%s.ovlp.reverse", output_file_name); write_ma_hit_ts(reverse_sources, RNF->total_reads, gfa_name); - if(asm_opt.dbg_bam) { + /**if(asm_opt.dbg_bam)**/ { sprintf(gfa_name, "%s.rec", output_file_name); write_cc_v(&scb, gfa_name); } @@ -23794,7 +23794,7 @@ int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_s return 0; } - if(asm_opt.dbg_bam) { + /**if(asm_opt.dbg_bam)**/ { sprintf(gfa_name, "%s.rec", output_file_name); load_cc_v(&scb, gfa_name); ///write_ec_reads("lec.raw.fa", &scb, 0); } diff --git a/Process_Read.cpp b/Process_Read.cpp index 717cb9e..5e3d18b 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -125,6 +125,11 @@ void write_All_reads(All_reads* r, char* read_file_name) } } + mm = 3; + fwrite(&mm, sizeof(mm), 1, fp); + fwrite(&(r->tr[0]), sizeof(r->tr[0]), 1, fp); + fwrite(&(r->tr[1]), sizeof(r->tr[1]), 1, fp); + free(index_name); fflush(fp); fclose(fp); @@ -235,11 +240,18 @@ int load_All_reads(All_reads* r, char* read_file_name) MALLOC(r->rsc[i], (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0)); f_flag += fread(r->rsc[i], sizeof(uint8_t), (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0), fp); } + + mm = 0; + if(!feof(fp)) fread(&mm, sizeof(mm), 1, fp); + } + + if(mm == 3) { + fread(&(r->tr[0]), sizeof(r->tr[0]), 1, fp); + fread(&(r->tr[1]), sizeof(r->tr[1]), 1, fp); } } - free(index_name); fclose(fp); fprintf(stderr, "Reads has been loaded.\n"); @@ -283,12 +295,17 @@ uint8_t load_cc_v(cc_v* r, char* read_file_name) ///typedef struct {size_t n, m; asg16_v *a; uint8_t *f; uint16_t *er; uint64_t bid;} cc_v; asg16_v *z; int f_flag = 0; uint64_t k, rn; uint32_t zn; f_flag += fread(&rn, sizeof(rn), 1, fp); - r->n = r->m = rn; MALLOC(r->a, r->n); - for (k = 0; k < r->n; k++) { - z = &(r->a[k]); - f_flag += fread(&zn, sizeof(zn), 1, fp); - z->n = z->m = zn; MALLOC(z->a, z->n); - f_flag += fread(z->a, sizeof((*(z->a))), zn, fp); + r->n = r->m = rn; + if(rn == 0) { + free(r->a); r->a = NULL; + } else { + MALLOC(r->a, r->n); + for (k = 0; k < r->n; k++) { + z = &(r->a[k]); + f_flag += fread(&zn, sizeof(zn), 1, fp); + z->n = z->m = zn; MALLOC(z->a, z->n); + f_flag += fread(z->a, sizeof((*(z->a))), zn, fp); + } } fflush(fp); fclose(fp); diff --git a/anchor.cpp b/anchor.cpp index 311fbff..2ec3cef 100644 --- a/anchor.cpp +++ b/anchor.cpp @@ -4014,6 +4014,156 @@ void get_pi_ec_chain(ha_abuf_t *ab, uint64_t rid, uint64_t rl, uint32_t tid, cha } +uint64_t gen_srt_chain(ha_abuf_t *ab, ha_mz1_t *rfa, uint64_t *rfi, uint64_t rfn, ha_mz1_t *qra, uint64_t *qri, uint64_t qrn, uint32_t *high_occ, uint32_t *low_occ, Candidates_list *cl, uint32_t qrid) +{ + uint64_t i, k, l, max_cnt = UINT32_MAX, min_cnt = 0, ri, qi, sn, x, rz, qz, n_a0, l0, k0, i0; anchor1_t *an; k_mer_hit *p; + if(high_occ) { + max_cnt = (*high_occ); + if(max_cnt < 2) max_cnt = 2; + } + if(low_occ) { + min_cnt = (*low_occ); + if(min_cnt < 2) min_cnt = 2; + } + + radix_sort_anc64(rfi, rfi + rfn); + radix_sort_anc64(qri, qri + qrn); + i = ab->n_a = 0; n_a0 = l0 = i0 = 0; k0 = 1; + for (k = 1, l = 0; k <= rfn; ++k) { + if (k == rfn || (rfi[k]>>32) != (rfi[l]>>32)) { + x = rfi[l]>>32; + for (; i < qrn && (qri[i]>>32) < x; i++); + if(ab->n_a < ab->m_a) { + n_a0 = ab->n_a; l0 = l; i0 = i; k0 = k; + } + if(i < qrn && (qri[i]>>32) == x) { + // sn = 0; + // if(ab->n_a < ab->m_a) sn = ab->seed[l].n;///ha_pt_cnt(ha_idx, ra[l].x); + for (qi = i; qi < qrn && (qri[qi]>>32) == x; qi++) { + qz = (uint32_t)qri[qi]; + for (ri = l; ri < k; ri++) { + rz = (uint32_t)rfi[ri]; + if(qra[qz].rev != rfa[rz].rev) continue; + if(ab->n_a < ab->m_a) { + an = &(ab->a[ab->n_a++]); + an->other_off = qra[qz].pos; + an->self_off = rfa[rz].pos; + ///an->cnt: cnt<<8|span + an->cnt = ab->seed[rz].n; if(an->cnt > ((uint32_t)(0xffffffu))) an->cnt = 0xffffffu; + an->cnt <<= 8; an->cnt |= ((rfa[rz].span <= ((uint32_t)(0xffu)))?rfa[rz].span:((uint32_t)(0xffu))); + an->srt = (((uint64_t)(an->self_off))<<32)|((uint64_t)(an->other_off)); + } else { + ab->n_a++; + } + } + } + i = qi; + } + l = k; + } + } + + if (ab->n_a > ab->m_a) { + ab->m_a = ab->n_a; REALLOC(ab->a, ab->m_a); + k = k0; i = i0; l = l0; ab->n_a = n_a0; + + for (; k <= rfn; ++k) { + if (k == rfn || (rfi[k]>>32) != (rfi[l]>>32)) { + x = rfi[l]>>32; + for (; i < qrn && (qri[i]>>32) < x; i++); + if(ab->n_a < ab->m_a) { + n_a0 = ab->n_a; l0 = l; i0 = i; k0 = k; + } + if(i < qrn && (qri[i]>>32) == x) { + // sn = 0; + // if(ab->n_a < ab->m_a) sn = ab->seed[l].n;///ha_pt_cnt(ha_idx, ra[l].x); + for (qi = i; qi < qrn && (qri[qi]>>32) == x; qi++) { + qz = (uint32_t)qri[qi]; + for (ri = l; ri < k; ri++) { + rz = (uint32_t)rfi[ri]; + if(qra[qz].rev != rfa[rz].rev) continue; + if(ab->n_a < ab->m_a) { + an = &(ab->a[ab->n_a++]); + an->other_off = qra[qz].pos; + an->self_off = rfa[rz].pos; + ///an->cnt: cnt<<8|span + an->cnt = ab->seed[rz].n; if(an->cnt > ((uint32_t)(0xffffffu))) an->cnt = 0xffffffu; + an->cnt <<= 8; an->cnt |= ((rfa[rz].span <= ((uint32_t)(0xffu)))?rfa[rz].span:((uint32_t)(0xffu))); + an->srt = (((uint64_t)(an->self_off))<<32)|((uint64_t)(an->other_off)); + } else { + ab->n_a++; + } + } + } + i = qi; + } + l = k; + } + } + } + + // copy over to _cl_ + sn = ab->n_a + cl->length; + if (sn > (uint64_t)cl->size) { + cl->size = sn; + REALLOC(cl->list, cl->size); + } + + radix_sort_ha_an1(ab->a, ab->a + ab->n_a); + for (i = 0; i < ab->n_a; i++) { + p = &cl->list[cl->length++]; + p->readID = qrid; + p->strand = 0; + p->offset = ab->a[i].other_off; + p->self_offset = ab->a[i].self_off; + + if(((ab->a[i].cnt>>8) < max_cnt) && ((ab->a[i].cnt>>8) > min_cnt)){ + p->cnt = 1; + } else if((ab->a[i].cnt>>8) <= min_cnt) { + p->cnt = 2; + } else{ + p->cnt = 1 + (((ab->a[i].cnt>>8) + (max_cnt<<1) - 1)/(max_cnt<<1)); + p->cnt = pow(p->cnt, 1.1); + } + if(p->cnt > ((uint32_t)(0xffffffu))) p->cnt = 0xffffffu; + p->cnt <<= 8; p->cnt |= (((uint32_t)(0xffu))&(ab->a[i].cnt)); + } + // cl->length = ab->n_a; + return ab->n_a; +} + +void gen_self_global_chain(ha_abuf_t *ab, Candidates_list *cl, uint32_t rid, uint64_t rl, uint32_t tid, char *ts, uint64_t tl, uint64_t mz_w, uint64_t mz_k, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, asg64_v *ix, + uint8_t is_accurate, double bw_thres, int apend_be, uint32_t gen_off, int64_t enable_mcopy, double mcopy_rate, uint32_t mcopy_khit_cut, overlap_region_alloc *ores) +{ + extern void *ha_flt_tab; + uint64_t k, rn = ab->mz.n, m = cl->length; + + mz1_ha_sketch(ts, tl, 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); + + kv_resize(uint64_t, *ix, ab->mz.n); + for (k = 0; k < rn; k++) { + ix->a[k] = ab->mz.a[k].x; + ix->a[k] <<= 32; ix->a[k] |= k; + } + + for (; k < ab->mz.n; k++) { + ix->a[k] = ab->mz.a[k].x; + ix->a[k] <<= 32; ix->a[k] |= (k-rn); + } + ix->n = ab->mz.n; + + if(gen_srt_chain(ab, ab->mz.a, ix->a, rn, ab->mz.a + rn, ix->a + rn, ab->mz.n - rn, high_occ, low_occ, cl, tid)) { + 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); + + m += lchain_qdp_mcopy_fast(cl, m, cl->length - m, m, &(cl->chainDP), ores, max_skip, max_iter, max_dis, chn_pen_gap, chn_pen_skip, bw_thres, + rid, rl, tl, quick_check, apend_be, gen_off, enable_mcopy, mcopy_rate, mcopy_khit_cut, 1); + cl->length = m; + } + ab->mz.n = rn; +} + + int64_t ug_map_lchain(ha_abufl_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, const ul_idx_t *uref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, double bw_thres_sec, int max_n_chain, int apend_be, kvec_t_u8_warp* k_flag, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, uint32_t is_accurate, uint32_t gen_off, double mcopy_rate, uint32_t mcopy_khit_cut, uint32_t is_hpc, ha_mzl_t *res, uint64_t res_n, ha_mzl_t *idx, uint64_t idx_n, uint64_t mzl_cutoff, uint64_t chain_cutoff, kv_u_trans_t *kov) diff --git a/ecovlp.cpp b/ecovlp.cpp index 9ad0666..03be9ec 100644 --- a/ecovlp.cpp +++ b/ecovlp.cpp @@ -13,6 +13,7 @@ #define del_cns_arc(z, arc_i) ((z).arc.a[(arc_i)].v == CNS_DEL_E) #define CNS_DEL_V (0x1fffffffu) #define del_cns_nn(z, nn_i) ((z).a[(nn_i)].sc == CNS_DEL_V) +#define is_cns_bb(z, nn_i) ((nn_i) >= (z).bb0 && (nn_i) < (z).bb1) #define REFRESH_N 128 #define COV_W 3072 #define COV_W_AC 512 @@ -242,6 +243,8 @@ uint64_t get_mz1(const char *str, int len, int w, int k, uint32_t rid, int is_hp void get_pi_ec_chain(ha_abuf_t *ab, uint64_t rid, uint64_t rl, uint32_t tid, char* ts, uint64_t tl, uint64_t mz_w, uint64_t mz_k, overlap_region_alloc *overlap_list, 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 enable_mcopy, double mcopy_rate, uint32_t mcopy_khit_cut, int64_t max_skip, int64_t max_iter, int64_t max_dis, int64_t quick_check, double chn_pen_gap, double chn_pen_skip); +void gen_self_global_chain(ha_abuf_t *ab, Candidates_list *cl, uint32_t rid, uint64_t rl, uint32_t tid, char *ts, uint64_t tl, uint64_t mz_w, uint64_t mz_k, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, asg64_v *ix, + uint8_t is_accurate, double bw_thres, int apend_be, uint32_t gen_off, int64_t enable_mcopy, double mcopy_rate, uint32_t mcopy_khit_cut, overlap_region_alloc *ores); void set_lchain_dp_op(uint32_t is_accurate, uint32_t mz_k, int64_t *max_skip, int64_t *max_iter, int64_t *max_dis, double *chn_pen_gap, double *chn_pen_skip, int64_t *quick_check); void h_ec_lchain_re_gen_srt(ha_abuf_t *ab, ha_pt_t *ha_idx, overlap_region_alloc *olst, Candidates_list *cl); uint64_t h_ec_lchain_re_gen_qry(ha_abuf_t *ab, uint64_t *k, uint64_t *l, uint64_t *i, uint64_t *idx_a, uint64_t idx_n, uint64_t *tid, uint64_t *trev); @@ -1430,7 +1433,174 @@ void del_cns_g_nn(cns_gfa *cns, uint32_t v) cns->a[v].arc.n = cns->a[v].arc.nou = 0; - cns->a[v].c = cns->a[v].f = 0; cns->a[v].sc = CNS_DEL_V; + /**cns->a[v].c =**/ cns->a[v].f = 0; cns->a[v].sc = CNS_DEL_V; +} + +void merge_cns_g_in_adv(cns_gfa *cns, uint32_t v0, asg32_v* b32) +{ + cns_t *av, *aw; uint32_t v, bp, vk, wk, wka, w, wn, nn, mn, wh, mn_k[2]; + + b32->n = 0; + kv_push(uint32_t, *b32, v0); + while (b32->n) { + v = b32->a[--b32->n]; + if(del_cns_nn((*cns), v)) continue; + + av = &((*cns).a[v]); + for (bp = 0; bp < 4; bp++) { + //nn: number of node; wh: weight + nn = wh = 0; mn = mn_k[0] = mn_k[1] = wka = (uint32_t)-1; + for (vk = av->arc.nou; vk < av->arc.n; vk++) {///in-edge of v + if(del_cns_arc((*av), vk)) continue; + w = av->arc.a[vk].v; aw = &((*cns).a[w]); + if(aw->c != bp) continue; + if(w == cns->si || w == cns->ei) continue; + + for (wk = wn = 0; wk < aw->arc.nou; wk++) {///out-edge of w + if(del_cns_arc((*aw), wk)) continue; + wn++; wka = wk; if(wn > 1) break; + } + + if(wn != 1) continue; + assert(aw->arc.a[wka].v == v); + + ///deal with out-edge of w + if((nn == 0) || is_cns_bb((*cns), w)) {///this is still risky as there might be multiple backbone + mn = w; mn_k[0] = vk; mn_k[1] = wka; + ///not sure if we should set these edges as visited + // av->arc.a[vk].f = 1; aw->arc.a[wka].f = 1; + } + wh += aw->arc.a[wka].sc; + nn++; + } + + + if(nn > 1) { + for (vk = av->arc.nou; vk < av->arc.n; vk++) {///in-edge of v + if(del_cns_arc((*av), vk)) continue; + w = av->arc.a[vk].v; aw = &((*cns).a[w]); + if(aw->c != bp) continue; + if(w == cns->si || w == cns->ei) continue; + + for (wk = wn = 0; wk < aw->arc.nou; wk++) {///out-edge of w + if(del_cns_arc((*aw), wk)) continue; + wn++; wka = wk; if(wn > 1) break; + } + + if(wn != 1) continue; + assert(aw->arc.a[wka].v == v); + + ///deal with in-edge of w + ///all edges to w, should be move to mn + if(mn != w) {///not sure if we should set these edges as visited; affect when nn == 0 + for (wk = aw->arc.nou; wk < aw->arc.n; wk++) { + if(del_cns_arc((*aw), wk)) continue; + ///previously, aw->arc.a[wk].v -> w + ///currently, aw->arc.a[wk].v -> mn + /// if(nn == 0), then mn = w + gen_mm_cns_arc(cns, aw->arc.a[wk].v, mn, aw->arc.a[wk].sc/**(nn?(aw->arc.a[wk].sc):(0))**/, aw->arc.a[wk].f);///not sure if we should set these edges as visited + } + del_cns_g_nn(cns, w); + } + } + } + + if(nn) { + aw = &((*cns).a[mn]); + av->arc.a[mn_k[0]].sc = aw->arc.a[mn_k[1]].sc = wh; + // merge_cns_g_in(cns_gfa *cns, uint32_t v, asg32_v* b32) + kv_push(uint32_t, *b32, mn); + } + } + } +} + +void merge_cns_g_ou_adv(cns_gfa *cns, uint32_t v0, asg32_v* b32) +{ + cns_t *av, *aw; uint32_t v, bp, vk, wk, wka, w, wn, nn, mn, wh, mn_k[2]; + + b32->n = 0; + kv_push(uint32_t, *b32, v0); + while (b32->n) { + v = b32->a[--b32->n]; + if(del_cns_nn((*cns), v)) continue; + + av = &((*cns).a[v]); + for (bp = 0; bp < 4; bp++) { + //nn: number of node; wh: weight + nn = wh = 0; mn = mn_k[0] = mn_k[1] = wka = (uint32_t)-1; + for (vk = 0; vk < av->arc.nou; vk++) {///ou-edge of v + if(del_cns_arc((*av), vk)) continue; + w = av->arc.a[vk].v; aw = &((*cns).a[w]); + if(aw->c != bp) continue; + if(w == cns->si || w == cns->ei) continue; + + for (wk = aw->arc.nou, wn = 0; wk < aw->arc.n; wk++) {///in-edge of w + if(del_cns_arc((*aw), wk)) continue; + wn++; wka = wk; if(wn > 1) break; + } + + if(wn != 1) continue; + + assert(aw->arc.a[wka].v == v); + + ///deal with out-edge of w + if((nn == 0) || is_cns_bb((*cns), w)) {///this is still risky as there might be multiple backbone + mn = w; mn_k[0] = vk; mn_k[1] = wka; + ///not sure if we should set these edges as visited + // av->arc.a[vk].f = 1; aw->arc.a[wka].f = 1; + } + wh += aw->arc.a[wka].sc; + nn++; + } + + + if(nn > 1) { + for (vk = 0; vk < av->arc.nou; vk++) {///ou-edge of v + if(del_cns_arc((*av), vk)) continue; + w = av->arc.a[vk].v; aw = &((*cns).a[w]); + if(aw->c != bp) continue; + if(w == cns->si || w == cns->ei) continue; + + for (wk = aw->arc.nou, wn = 0; wk < aw->arc.n; wk++) {///in-edge of w + if(del_cns_arc((*aw), wk)) continue; + wn++; wka = wk; if(wn > 1) break; + } + + if(wn != 1) continue; + + assert(aw->arc.a[wka].v == v); + + + + + ///deal with in-edge of w + ///all edges to w, should be move to mn + if(mn != w) {///not sure if we should set these edges as visited; affect when nn == 0 + for (wk = 0; wk < aw->arc.nou; wk++) { + if(del_cns_arc((*aw), wk)) continue; + ///previously, w -> aw->arc.a[wk].v + ///currently, mn -> aw->arc.a[wk].v + /// if(nn == 0), then mn = w + gen_mm_cns_arc(cns, mn, aw->arc.a[wk].v, aw->arc.a[wk].sc/**(nn?(aw->arc.a[wk].sc):(0))**/, aw->arc.a[wk].f);///not sure if we should set these edges as visited + } + del_cns_g_nn(cns, w); + } + } + } + + if(nn) { + aw = &((*cns).a[mn]); + // fprintf(stderr, "\n[M::%s] nn::%u, mn::%u\n", __func__, nn, mn); + // fprintf(stderr, "[M::%s] vi::%u, vn::%u\n", __func__, mn_k[0], (uint32_t)av->arc.n); + // fprintf(stderr, "[M::%s] wi::%u, wn::%u\n", __func__, mn_k[1], (uint32_t)aw->arc.n); + + av->arc.a[mn_k[0]].sc = aw->arc.a[mn_k[1]].sc = wh; + // merge_cns_g_in(cns_gfa *cns, uint32_t v, asg32_v* b32) + kv_push(uint32_t, *b32, mn); + } + } + } } void merge_cns_g_in(cns_gfa *cns, uint32_t v0, asg32_v* b32) @@ -1936,6 +2106,9 @@ uint64_t push_correct1_fhc_indel_exz(asg16_v *sc, int64_t sc0, window_list *idx, uint64_t push_correct1_fhc(window_list *idx, window_list_alloc *res, cns_gfa *cns, char* qstr, UC_Read* tu, bit_extz_t *exz, asg32_v *rc, uint32_t bl, uint32_t rid) { + // if(rid == 11206 && bl == 6) { + // fprintf(stderr, "[M::%s]\trc->n::%u\tbl::%u\trid::%u\n", __func__, (uint32_t)rc->n, bl, rid); + // } // fprintf(stderr, "[M::%s]\trc->n::%u\tbl::%u\n", __func__, (uint32_t)rc->n, bl); uint64_t nec = 0; uint32_t k, l, i, ff, sl, sk, bs = cns->off, be = bl + cns->off, bend = cns->off, is_i = 0, sc0 = res->c.n;///[bs, be) if(rc->n) {///it is possible that rc->n == 0, which means there is a deletion @@ -1983,13 +2156,17 @@ uint64_t push_correct1_fhc(window_list *idx, window_list_alloc *res, cns_gfa *cn l = k; } } - + + // if(rid == 11206 && bl == 6) { + // fprintf(stderr, "[M::%s]\trc->n::%u\tbl::%u\trid::%u\tbend::%u\tbe::%u\n", __func__, (uint32_t)rc->n, bl, rid, bend, be); + // } ///push remaining deletion if(be > bend) { if(is_i && exz) { nec += push_correct1_fhc_indel_exz(((asg16_v *)(&(res->c))), sc0, idx, cns, qstr, tu, exz, cns->off, 3, be - bend, bend-cns->off); } else { for (i = bend; i < be; i++) { + // fprintf(stderr, "%c(%u)\n", s_H[cns->a[i].c], i); push_trace_bp_f(((asg16_v *)(&(res->c))), 3, cns->a[i].c, (uint16_t)-1, 1, ((idx->clen>0)?1:0)); idx->clen = res->c.n - idx->cidx; nec++; } @@ -2017,6 +2194,23 @@ uint64_t cns_gen_full0(overlap_region* ol, All_reads *rref, uint64_t s, uint64_t } init_cns_g(cns, qstr + s, e - s, hw, rid); + + // uint64_t dk = 0; + // if(s <= 29788 && e > 29788) { + // fprintf(stderr, "[M::%s]\tqid::%u\tq::[%lu, %lu)\n", __func__, rid, s, e); + // fprintf(stderr, "-z-[M::%s]\toriginal::", __func__); + // for (dk = s; dk < e; dk++) { + // fprintf(stderr, "%c", qstr[dk]); + // } + // fprintf(stderr, "\n"); + + // fprintf(stderr, "-a-[M::%s]\tGraph::\t\t\t", __func__); + // for (dk = 2; dk < (*cns).n; dk++) { + // fprintf(stderr, "%c", s_H[(*cns).a[dk].c]); + // } + // fprintf(stderr, "\n"); + // } + id_n = iter_cc_idx_t(ol, idx, s, e, idx->rr, ((s==e)?1:0), &id_a); // debug_inter0(ol, idx->c_idx, idx->idx->a + idx->i0, idx->srt_n - idx->i0, id_a, id_n, s, e, ((s==e)?1:0), 0, "-1-"); uint64_t k, q[2], os, oe; ul_ov_t *p; overlap_region *z; idx->rr = 0; @@ -2049,16 +2243,42 @@ uint64_t cns_gen_full0(overlap_region* ol, All_reads *rref, uint64_t s, uint64_t } } + // if(s <= 29788 && e > 29788) { + // fprintf(stderr, "-b-[M::%s]\tGraph::\t\t\t", __func__); + // for (dk = 2; dk < (*cns).n; dk++) { + // fprintf(stderr, "%c(del::%u)", s_H[(*cns).a[dk].c], del_cns_nn((*cns), dk)); + // } + // fprintf(stderr, "\n"); + // } + // fprintf(stderr, "-2-[M::%s] cns->n::%u\n", __func__, (uint32_t)cns->n); // return; refine_cns_g(cns, b32); + + // if(s <= 29788 && e > 29788) { + // fprintf(stderr, "-c-[M::%s]\tGraph::\t\t\t", __func__); + // for (dk = 2; dk < (*cns).n; dk++) { + // fprintf(stderr, "%c(del::%u)", s_H[(*cns).a[dk].c], del_cns_nn((*cns), dk)); + // } + // fprintf(stderr, "\n"); + // } + // fprintf(stderr, "-3-[M::%s] cns->n::%u\n", __func__, (uint32_t)cns->n); gseq_cns_g(cns, b32, e - s); + + // if(s <= 29788 && e > 29788) { + // fprintf(stderr, "-d-[M::%s]\tGraph::\t\t\t", __func__); + // for (dk = 2; dk < (*cns).n; dk++) { + // fprintf(stderr, "%c(del::%u)", s_H[(*cns).a[dk].c], del_cns_nn((*cns), dk)); + // } + // fprintf(stderr, "\n"); + // } + // fprintf(stderr, "-4-[M::%s] cns->n::%u\n", __func__, (uint32_t)cns->n); // nec += push_correct1(ridx, res, cns, b32, e - s); @@ -2392,7 +2612,7 @@ uint64_t wcns_vote(overlap_region* ol, All_reads *rref, uint64_t q_hf, char* qst } ///a) pass coverage check; b) no enough coverage if((((oc[0] > (oc[1]*occ_exact)) && (oc[0] > (oc[1]-oc[0])) && (oc[1] >= occ_tot) && (oc[0] > 1) && ((!hf_idx) || ((ow[0] > (ow[1]*occ_exact)) && (ow[0] > (ow[1] - ow[0])))))) || (oc[1] < occ_tot)) { - fI = 0; + fI = 0;///keep q itself } if(fI) { @@ -2434,8 +2654,6 @@ uint64_t wcns_vote(overlap_region* ol, All_reads *rref, uint64_t q_hf, char* qst return rr; } - - void print_debug_ovlp_cigar(overlap_region_alloc* ol, asg64_v* idx, kv_ul_ov_t *c_idx) { uint64_t k, ci; uint32_t cl; ul_ov_t *cp; bit_extz_t ez; uint16_t c; char cm[4]; @@ -2533,7 +2751,7 @@ uint64_t wcns_gen(overlap_region_alloc* ol, All_reads *rref, uint64_t qid, UC_Re int64_t srt_n = idx->n, s, e, t, rr; i = 0; radix_sort_ec64(idx->a, idx->a+idx->n); - for (k = 1, i = 0; k < srt_n; k++) { + for (k = 1, i = 0; k <= srt_n; k++) { if (k == srt_n || (idx->a[k]>>32) != (idx->a[i]>>32)) { if(k - i > 1) { for (t = i; t < k; t++) { @@ -2555,10 +2773,11 @@ uint64_t wcns_gen(overlap_region_alloc* ol, All_reads *rref, uint64_t qid, UC_Re ///second index kv_resize(ul_ov_t, *c_idx, (c_idx->n<<1)); ul_ov_t *idx_a = NULL, *idx_b = NULL; - idx_a = c_idx->a; idx_b = c_idx->a; + idx_a = c_idx->a; idx_b = c_idx->a + c_idx->n;///update here? memcpy(idx_b, idx_a, c_idx->n * (sizeof((*(idx_a))))); kv_resize(uint64_t, *buf, ((wl<<1) + idx->n)); buf->n = ((wl<<1) + idx->n); + ///buf->a[wl<<1, idx->n): this is for sorted index memcpy(buf->a + (wl<<1), idx->a, idx->n * (sizeof((*(idx->a))))); memset(buf->a, 0, (wl<<1)*(sizeof((*(idx->a))))); if(hf_idx) { @@ -4340,15 +4559,157 @@ static void worker_init_ec_step(void *data, long i, int tid) void update_scb(All_reads *R_INF, asg16_v *scc, asg16_v *scb, asg16_v *scb_res, UC_Read *qu, UC_Read *tu, asg64_v *srt, bit_extz_t *exz, uint64_t rid); +uint64_t gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t rid); + +uint8_t refresh_check_scc(overlap_region *z, int64_t sc_e, int64_t ql, int64_t tl, int64_t gap_bd, double gap_rate, double err_rate, int64_t err_diff_bd, double err_diff_rate) +{ + int64_t k, wn = z->w_list.n, tot_g = 0, tot_e = 0, wq, wt; + if(wn <= 0) return 0; + + tot_g += z->w_list.a[0].x_start; tot_g += ql - z->w_list.a[wn-1].x_end - 1; + tot_g += z->w_list.a[0].y_start; tot_g += tl - z->w_list.a[wn-1].y_end - 1; + + for (k = 0; k < wn; k++) { + if(is_ualn_win(z->w_list.a[k])) { + if(z->w_list.a[k].extra_begin == -1 && z->w_list.a[k].extra_end == -1) { + return 0; + } + tot_e += (z->w_list.a[k].extra_begin*(-1)) - 1; + // fprintf(stderr, "+[M::%s]\tq::[%u,%u)\tt::[%u,%u)\ttot_e::%ld\n", __func__, 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, tot_e); + } else { + tot_e += z->w_list.a[k].error; + // fprintf(stderr, "-[M::%s]\tq::[%u,%u)\tt::[%u,%u)\ttot_e::%ld\n", __func__, 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, tot_e); + } + } + + // fprintf(stderr, "[M::%s]\tq::[%u,%u)\tt::[%u,%u)\ttot_e::%ld\trtot_g::%ld\n", __func__, z->x_pos_s, z->x_pos_e + 1, z->y_pos_s, z->y_pos_e + 1, tot_e, tot_g); + + if((tot_g > 0) && (tot_g > gap_bd)) { + if(tot_g > (ql*gap_rate)) return 0; + if(tot_g > (tl*gap_rate)) return 0; + } + + + wq = (z->w_list.a[wn-1].x_end+1) - z->w_list.a[0].x_start; + wt = (z->w_list.a[wn-1].y_end+1) - z->w_list.a[0].y_start; + if((tot_e > 0) && ((tot_e > (wq*err_rate)) || (tot_e > (wt*err_rate)))) return 0; + + // tot_e += tot_g; + if((tot_e > 0) && ((tot_e > (ql*err_rate)) || (tot_e > (tl*err_rate)))) return 0; + + if((tot_e > sc_e) && ((tot_e - sc_e) > err_diff_bd) && ((tot_e - sc_e) > (sc_e*err_diff_rate))) return 0; + + return 1; + +} + + +void dbg_prt_asg16_v_sc(uint64_t rid, asg16_v *sc) +{ + uint64_t ck, qk, tk, tot_e = 0; uint32_t len; uint16_t c, bq, bt; + + fprintf(stderr, "[M::%s]\trid::%lu\trlen::%lu\n", __func__, rid, Get_READ_LENGTH(R_INF, rid)); + ck = qk = tk = 0; + while (ck < sc->n) { + ck = pop_trace_bp_f(sc, ck, &c, &bq, &bt, &len); + if(c != 2) qk += len; + if(c != 3) tk += len; + if(c!=0) tot_e += len; + // fprintf(stderr, "%u(%c)\t", len, "MSID"[c]); + } + // fprintf(stderr, "\n"); + fprintf(stderr, "[M::%s]\trid::%lu\t#\talter::%lu\n", __func__, rid, tot_e); +} + +uint8_t regen_scb(ha_abuf_t *ab, Candidates_list *cl, uint32_t rid, asg16_v *sc, UC_Read *ia, UC_Read *ob, asg64_v *srt, + uint64_t mz_w, uint64_t mz_k, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* dbg_ct, st_mt_t *sp, uint32_t *high_occ, uint32_t *low_occ, + uint8_t is_accurate, double bw_thres, int apend_be, uint32_t gen_off, overlap_region_alloc *ol, bit_extz_t *exz, double erate, uint64_t wl, asg16_v *buf, + overlap_region **r_aux_o, overlap_region **rchn, int64_t *rref_len) +{ + // return 0; + // dbg_prt_asg16_v_sc(rid, sc); + uint64_t oln0 = ol->length, e0; (*rchn) = NULL; (*rref_len) = -1; + e0 = gen_ori_seq0(ia->seq, ia->length, ob, sc, rid); + // fprintf(stderr, "[M::%s]\tinitial\tlength::%ld\n", __func__, (int64_t)ob->length); + // fprintf(stderr, "\n[M::%s]\trid::%u\trlen::%ld\te0::%ld\n", __func__, rid, (int64_t)ia->length, e0); + if(e0 == 0) return 0; + cl->length = 0; + gen_self_global_chain(ab, cl, rid, ia->length, rid/**((uint32_t)-1)**/, ob->seq, ob->length, mz_w, mz_k, k_flag, dbg_ct, sp, high_occ, low_occ, srt, is_accurate, bw_thres, apend_be, gen_off, 1, -1, UINT32_MAX, ol); + (*r_aux_o) = fetch_aux_ovlp(ol, NULL);///refresh + if(ol->length > oln0) { + assert(ol->length == oln0 + 1); + if(gen_hc_r_alin_self(&(ol->list[oln0]), cl, ia->seq, ia->length, ob->seq, ob->length, exz, (*r_aux_o), erate, wl, rid, E_KHIT, 1, buf, sc)) { + if(refresh_check_scc(&(ol->list[oln0]), e0, ia->length, ob->length, 32, 0.012, 0.1, 32, 0.2)) { + (*rchn) = &(ol->list[oln0]); (*rref_len) = ob->length; + resize_UC_Read(ia, ia->length + ob->length); memcpy(ia->seq + ia->length, ob->seq, sizeof((*(ob->seq)))*ob->length); + ol->length = oln0; + return 1; + } + } + } + ol->length = oln0; + return 0; +} + +void cmp_smp_ac(UC_Read *ref, asg16_v *scz, UC_Read *res, uint64_t rid) +{ + uint64_t qn = ref->length; uint16_t c, bq, bt; uint32_t len, ck, qk, tk, tn, wq[2], wt[2], k/**, Nn = 0**/; char *qsr = NULL, *tsr = NULL; + + ck = qk = tk = 0; + while (ck < scz->n) { + ck = pop_trace_bp_f(scz, ck, &c, &bq, &bt, &len); + if(c != 3) tk += len; + } + tn = tk; resize_UC_Read(ref, qn + tn); + + qsr = ref->seq; tsr = ref->seq + qn; + + ck = qk = tk = 0; + while (ck < scz->n) { + wq[0] = qk; wt[0] = tk; + ck = pop_trace_bp_f(scz, ck, &c, &bq, &bt, &len); + if(c != 2) qk += len; + if(c != 3) tk += len; + wq[1] = qk; wt[1] = tk; + // fprintf(stderr, "%u(%c)\tq::[%u,%u)\tbq::%u\tt::[%u,%u)\tbt::%u\n", len, "MSID"[c], wq[0], wq[1], bq, wt[0], wt[1], bt); + if(c == 0) { + for (; wq[0] < wq[1]; wq[0]++, wt[0]++) { + tsr[wt[0]] = qsr[wq[0]]; + // if(p->a[wy[0]] == 'N') Nn++; + } + } else if(c == 1 || c == 2) { + for (k = wt[0]; k < wt[1]; k++) { + tsr[k] = s_H[bt]; + // if(p->a[k] == 'N') Nn++; + } + } + } + + gen_ori_seq0(tsr, tn, res, scz, rid); + + assert(res->length == ((int64_t)qn)); + if(memcmp(qsr, res->seq, qn) != 0) { + fprintf(stderr, "[M::%s]\trid::%lu(%.*s)\tres->length::%ld\tqn::%lu\n", __func__, rid, (int)Get_NAME_LENGTH(R_INF, rid), Get_NAME(R_INF, rid), (int64_t)res->length, qn); + for (k = 0; k < qn; k++) { + if(res->seq[k] != qsr[k]) { + fprintf(stderr, "[M::%s]\tk::%u\tqsr[%u]::%c\tres->seq[%u]::%c\n", __func__, k, k, qsr[k], k, res->seq[k]); + } + } + + } + // assert(memcmp(qsr, res->seq, qn) == 0); +} + static void worker_hap_ec(void *data, long i, int tid) { ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO); - uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; int64_t het_a, hom_a; ///gen_hc_aln_t ez; - overlap_region *aux_o = NULL/**, *rse_o = NULL**/; asg64_v buf0; uint32_t qlen = 0, qw = 0; uint64_t tot_b = 0; double tt0 = 0, tt1 = 0; + uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; int64_t het_a, hom_a, rl0 = -1; ///gen_hc_aln_t ez; + overlap_region *aux_o = NULL/**, *rse_o = NULL**/, *rcc = NULL; asg64_v buf0; uint32_t qlen = 0, qw = 0; uint64_t tot_b = 0; double tt0 = 0, tt1 = 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); + // if((i != 733166) && (i != 858708) && (i != 858732) && (i != 859819) && (i != 859899) && (i != 863486) && (i != 872165) && (i != 899887) && (i != 902298) && // (i != 906808) && (i != 946173) && (i != 952685) && (i != 983977) && (i != 1000227) && (i != 1011228) && (i != 1042858) && (i != 1045860) && (i != 1118558) && // (i != 1143886) && (i != 1155956) && (i != 1159490) && (i != 1179151) && (i != 1180199) && (i != 1230524) && (i != 1232338) && (i != 1244031) && (i != 1268467) && @@ -4502,11 +4863,15 @@ static void worker_hap_ec(void *data, long i, int tid) gen_reseed_re(&b->olist, &b->clist, aux_o, rse_o, &R_INF, &b->self_read, &b->ovlp_read, &b->exz, &b->pidx, &b->v64, &buf0, 0, asm_opt.mz_win, 19, i, asm_opt.max_ov_diff_ec, asm_opt.max_ov_diff_ec, &b->v16, R_INF.tqn, b->v8q.a); copy_asg_arr(b->sp, buf0); **/ + if(scb.a[i].n) { + regen_scb(b->ab, &b->clist, i, &(scb.a[i]), &b->self_read, &b->ovlp_read, &b->v64, asm_opt.mz_win, asm_opt.k_mer_length, NULL, NULL, &(b->sp), &high_occ, &low_occ, + 1, ((asm_opt.is_ont)?(0.05):(0.02)), 1, 1, &b->olist, &b->exz, ((asm_opt.max_ov_diff_ec>0.1)?(asm_opt.max_ov_diff_ec):(0.1)), (asm_opt.is_ont)?(WINDOW_OHC):(WINDOW_HC), &b->v16, &aux_o, &rcc, &rl0); + } copy_asg_arr(buf0, b->sp); //site_sc: r765 -> r766: 1 -> 0 - 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/**, 0**/, 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)/**&(b->v8t)**/, (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); + rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &(scb.a[i]), &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, 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)/**&(b->v8t)**/, (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, rcc, rl0); copy_asg_arr(b->sp, buf0); ///for debug indel // stderr_phase_ovlp(&b->olist); @@ -4540,7 +4905,7 @@ static void worker_hap_ec(void *data, long i, int tid) push_nec_re(aux_o, &(scc.a[i])); // push_nec_re(aux_o, &(scb.a[i])); - if(asm_opt.dbg_bam) { + /**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))); } @@ -4630,7 +4995,6 @@ 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 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) { @@ -4699,7 +5063,6 @@ static void worker_gfa_ec(void *data, long i, int tid) refresh_gc_ovec_buf_t0(b, REFRESH_N); } - void worker_hap_ec_back_dbg(void *data, long i, int tid) { ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); @@ -4817,8 +5180,8 @@ void worker_hap_ec_back_dbg(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); //site_sc: r765 -> r766: 1 -> 0 - 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/**, 0**/, 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)**/&(b->v8t), (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); + rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, 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)**/&(b->v8t), (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, NULL, -1); copy_asg_arr(b->sp, buf0); ///for debug indel // stderr_phase_ovlp(&b->olist); @@ -4940,8 +5303,6 @@ void worker_hap_ec_back_dbg(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)); } - - static void worker_hap_ec_step(void *data, long i, int tid) { ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); i += scc.bid; @@ -5085,8 +5446,8 @@ static void worker_hap_ec_step(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); //site_sc: r765 -> r766: 1 -> 0 - 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/**, 0**/, 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)**/&(b->v8t), (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); + rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, 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)**/&(b->v8t), (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, NULL, -1); copy_asg_arr(b->sp, buf0); ///for debug indel // stderr_phase_ovlp(&b->olist); @@ -5207,9 +5568,6 @@ static void worker_hap_ec_step(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)); } - - - static void worker_hap_ec_ss(void *data, long i, int tid) { ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); @@ -5303,8 +5661,8 @@ static void worker_hap_ec_ss(void *data, long i, int tid) copy_asg_arr(buf0, b->sp); //site_sc: r765 -> r766: 1 -> 0 - 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/**, 0**/, 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)**/&(b->v8t), (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); + rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, 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)**/&(b->v8t), (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, NULL, -1); copy_asg_arr(b->sp, buf0); ///for debug indel // stderr_phase_ovlp(&b->olist); @@ -5335,13 +5693,12 @@ static void worker_hap_ec_ss(void *data, long i, int tid) } - static void worker_hap_ec_hybrid(void *data, long i, int tid) { ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]); uint32_t high_occ = asm_opt.hom_cov * (2.0 - HA_KMER_GOOD_RATIO); int64_t het_a, hom_a; - uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; double bw_h, bw_l, e_h, e_l; - gen_hc_aln_t ez; overlap_region *aux_o = NULL, *rse_o = NULL; asg64_v buf0, buf1; uint64_t qlen = 0, qw = 0, qid = i; //uint64_t sk[2], ek[2], fn, qid = i, nec; + uint32_t low_occ = asm_opt.hom_cov * HA_KMER_GOOD_RATIO; double bw_h, bw_l, e_h, e_l; int64_t rl0 = -1; + gen_hc_aln_t ez; overlap_region *aux_o = NULL, *rse_o = NULL, *rcc = NULL; asg64_v buf0, buf1; uint64_t qlen = 0, qw = 0, qid = i; //uint64_t sk[2], ek[2], fn, qid = i, nec; if(qid < R_INF.tqn) {///ont bw_h = 0.05; bw_l = 0.035; e_h = asm_opt.max_ov_diff_ec; e_l = (asm_opt.max_ov_diff_ec + asm_opt.max_ov_diff_ec_sec)/2; } else { ///HiFi @@ -5353,16 +5710,20 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid) } 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); + // if(i != 10) return; + // if((i%16) != 0) return; - // if(i != 5966) return; + // if(i != 11206) return; // fprintf(stderr, "-a-[M::%s] rid::%ld\n", __func__, i); //id:i:21102 // if (memcmp("a59fab4a-892b-4ab7-bf4b-926bed57865b_1", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { //id:i:3504 - // if (memcmp("485f7963-eeb4-4745-ab74-1be4d61460c3", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { - // fprintf(stderr, "-a-[M::%s-beg] rid->%ld, rlen->%lu\n", __func__, i, Get_READ_LENGTH((R_INF),i)); + + + // if (memcmp("c7ecbd6b-e09d-4042-93ac-2400839feaf6", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) { + // fprintf(stderr, "-a-[M::%s-beg] rid->%ld, rlen->%lu, scb.a[i].n::%u\n", __func__, i, Get_READ_LENGTH((R_INF),i), (uint32_t)scb.a[i].n); // } else { // return; // } @@ -5434,12 +5795,24 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid) ///for debug indel // prt_ovlp_sam(&b->olist, &b->ovlp_read, b->self_read.seq, b->self_read.length); + if(scb.a[i].n) { + regen_scb(b->ab, &b->clist, i, &(scb.a[i]), &b->self_read, &b->ovlp_read, &b->v64, asm_opt.mz_win, asm_opt.k_mer_length, NULL, NULL, &(b->sp), &high_occ, &low_occ, + 1, bw_h, 1, 1, &b->olist, &b->exz, ((e_h>0.1)?(e_h):(0.1)), (qid < R_INF.tqn)?(WINDOW_OHC):(WINDOW_HC), &b->v16, &aux_o, &rcc, &rl0); + // if(rcc) { + // fprintf(stderr, "[M::%s]\t%.*s(qid::%u)\tql::%ld\tq::[%u,\t%u)\t%c\t%.*s(tid::%u)\ttl::%ld\tt::[%u,\t%u)\terr::%u\n", __func__, + // (int32_t)Get_NAME_LENGTH(R_INF, rcc->x_id), Get_NAME(R_INF, rcc->x_id), rcc->x_id, (int64_t)b->self_read.length, rcc->x_pos_s, rcc->x_pos_e + 1, "+-"[rcc->y_pos_strand], + // (int32_t)Get_NAME_LENGTH(R_INF, rcc->y_id), Get_NAME(R_INF, rcc->y_id), rcc->y_id, rl0, rcc->y_pos_s, rcc->y_pos_e + 1, rcc->non_homopolymer_errors); + // } else { + // fprintf(stderr, "[M::%s]\tunalined\n", __func__); + // } + } + 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/**, 0**/, 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)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, R_INF.tqn, 0/**1**/, HC0_W, &b->v32, - asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, (((double)R_INF.tr[1])/((double)(R_INF.tr[0] + R_INF.tr[1])))); + rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &(scb.a[i]), &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, 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)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, R_INF.tqn, 0/**1**/, HC0_W, &b->v32, + asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, (((double)R_INF.tr[1])/((double)(R_INF.tr[0] + R_INF.tr[1]))), rcc, rl0); copy_asg_arr(b->sp, buf0); ///for debug indel - // if(i == 23863) stderr_phase_ovlp(&b->olist); + // stderr_phase_ovlp(&b->olist); // exit(1); dedup_chains(&b->olist); @@ -5449,11 +5822,17 @@ static void worker_hap_ec_hybrid(void *data, long i, int tid) R_INF.tr[0], R_INF.tr[1], asm_opt.ont_rate, asm_opt.hf_rate, asm_opt.hf_rate_max, &buf1); copy_asg_arr(b->sp, buf0); copy_asg_arr(b->hap.snp_srt, buf1); - push_nec_re(aux_o, &(scc.a[i])); + // if(DBG_TIME && dbg_a) { + // dbg_a[i].faln = b->cnt[1]; + // } + + push_nec_re(aux_o, &(scc.a[i])); + // cmp_smp_ac(&b->self_read, &(scc.a[i]), &b->ovlp_read, i);///for debug // push_nec_re(aux_o, &(scb.a[i])); - if(asm_opt.dbg_bam) { + /**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); memcpy(scb.a[i].a, b->v16.a, b->v16.n); + kv_resize(uint16_t, scb.a[i], b->v16.n); memcpy(scb.a[i].a, b->v16.a, b->v16.n * sizeof((*(b->v16.a)))); scb.a[i].n = b->v16.n; + // fprintf(stderr, "-b-[M::%s-beg] rid->%ld, rlen->%lu, scb.a[i].n::%u\n", __func__, i, Get_READ_LENGTH((R_INF),i), (uint32_t)scb.a[i].n); } // if((asm_opt.is_ont) && is_chemical_r_qual(&b->olist, &b->v64, qlen, 1, 16, &(b->v8q), i)/**(is_uncorrected_read(&b->olist, &b->v64, qlen, 1600))**/) { @@ -5607,8 +5986,8 @@ static void worker_hap_ec_hybrid_sync(void *data, long i, int tid) 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/**, 0**/, 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)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, R_INF.tqn, 0/**1**/, HC0_W, &b->v32, - asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, hf_rate); + rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, 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)**/&(b->v8q), ((asm_opt.is_sc)?&(b->v8t):NULL), (asm_opt.is_ont)?1:0, R_INF.tqn, 0/**1**/, HC0_W, &b->v32, + asm_opt.s_hap_cov, asm_opt.infor_cov, het_a, hom_a, asm_opt.polyploidy, hf_rate, NULL, -1); copy_asg_arr(b->sp, buf0); ///for debug indel // if(i == 23863) stderr_phase_ovlp(&b->olist); @@ -8090,9 +8469,9 @@ overlap_region* h_ec_lchain_re3(ha_abuf_t *ab, uint32_t rid, UC_Read *qu, UC_Rea } -void gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t rid) +uint64_t gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t rid) { - uint64_t ck, qk, tk, k, wq[2], wt[2]; uint32_t len; uint16_t c, bq, bt; char *qstr = NULL; + uint64_t ck, qk, tk, k, wq[2], wt[2], tot_e = 0; uint32_t len; uint16_t c, bq, bt; char *qstr = NULL; ck = qk = tk = 0; while (ck < sc->n) { @@ -8101,7 +8480,14 @@ void gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t ri if(c != 2) qk += len; if(c != 3) tk += len; wq[1] = qk; wt[1] = tk; + if(c!=0) tot_e += len; + // if(rid == 24) { + // fprintf(stderr, "%u(%c)\t", len, "MSID"[c]); + // } } + // if(rid == 24) { + // fprintf(stderr, "\n"); + // } // if(!(tk == tl)) { // if(rid == 8) { // fprintf(stderr, "[M::%s] rid::%lu, tk::%lu, tl::%lu\n", __func__, rid, tk, tl); @@ -8116,6 +8502,9 @@ void gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t ri // } // } + // } + // if((!(tk == tl)) && (rid == 24)) { + // fprintf(stderr, "[M::%s]\trid::%lu\tqk::%lu\ttk::%lu\ttl::%lu\n", __func__, rid, qk, tk, tl); // } assert(tk == tl); @@ -8135,6 +8524,8 @@ void gen_ori_seq0(char *tstr, uint64_t tl, UC_Read *qu, asg16_v *sc, uint64_t ri } // fprintf(stderr, "%u%c(%c)(x::[%lu,%ld))(y::[%lu,%ld))\n", len, cm[c], ((c==1)||(c==2))?(cc[bt]):('*'), wx[0], wx[1], wy[0], wy[1]); // s_H } + + return tot_e; } void gen_cc_fly(asg16_v *sc, char *qstr, uint64_t ql, char *tstr, uint64_t tl, bit_extz_t *exz, double e_rate, uint64_t maxn, uint64_t maxe) @@ -8259,6 +8650,10 @@ void cal_updated_trace_len(asg16_v *sc, uint64_t *ql, uint64_t *tl) *ql = qk; *tl = tk; } +///qstr:: latest; tstr:: original; there is an intermidate string I between qstr and tstr +///tcc:: tstr -> I; +///qcc:: I -> qstr; +///tcc_res:: tstr -> qstr void gen_updated_trace(asg16_v *qcc, asg16_v *tcc, asg16_v *tcc_res, char *qstr, uint64_t ql, char *tstr, uint64_t tl, asg64_v *srt, bit_extz_t *exz, uint64_t rid) { uint64_t k, ck, qk, tk, wq[2], wt[2], old_dp, dp, s, e, srt_n, si, ei, so, os, oe, *qd, *td, qs, qe, ts, te, q0, t0; @@ -8466,8 +8861,8 @@ void update_scb(All_reads *R_INF, asg16_v *scc, asg16_v *scb, asg16_v *scb_res, // if(i == 700) fprintf(stderr, "|%u%c(%c)(x::%u)(y::%u)", len, cm[c], ((c==1)||(c==2))?(cc[b]):('*'), wx[1], wy[1]); // s_H } - qstr = tstr; ql = tl; - tstr = tu->seq; tl = tu->length; + qstr = tstr; ql = tl;///latest version + tstr = tu->seq; tl = tu->length;///orginal version // fprintf(stderr, "\n[M::%s] ql::%lu, tl::%lu, rid::%lu\n", __func__, ql, tl, rid); @@ -8575,8 +8970,8 @@ static void worker_hap_dc_ec0(void *data, long i, int tid) b->cnt[0] += b->self_read.length; 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/**, 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)**/&(b->v8q), ((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); + rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, NULL, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 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)**/&(b->v8q), ((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, NULL, -1); copy_asg_arr(b->sp, buf0); copy_asg_arr(buf0, b->sp); @@ -9299,7 +9694,7 @@ uint64_t cal_ec_multiple_step(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, u // 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); + // exit(1); (*r_base) = num_base; return num_correct; @@ -9567,7 +9962,7 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u // prt_nel_ovlp(R_INF.paf, n_a); // exit(1); - // dbg_write_ec_reads("ec12.fa", round, &scb, is_cr); + // dbg_write_ec_reads("ec12.fa", round, &scb, 0/**is_cr**/); if((!is_sv) || (is_sv && is_cr)) { kt_for(n_thre, worker_hap_post_rev, b, n_a); @@ -9582,7 +9977,7 @@ void cal_ec_r(uint64_t n_thre, uint64_t round, uint64_t n_round, uint64_t n_a, u fprintf(stderr, "-4-[M::%s]\t# tqn::%lu, Ont base::%lu, # HiFi bases::%lu\n", __func__, R_INF.tqn, R_INF.tr[0], R_INF.tr[1]); - // dbg_write_ec_reads("ec16.fa", round, &scb, !is_cr); + // dbg_write_ec_reads("ec16.fa", round, &scb, 0/**!is_cr**/); // exit(1); // uint64_t z; @@ -9753,12 +10148,12 @@ void cal_ov_r(uint64_t n_thre, uint64_t n_a, uint64_t new_idx) b = gen_ec_ovec_buf_t(n_thre); if(new_idx) { // kt_for(n_thre, worker_hap_dc_ec, b, n_a);///update overlaps - destroy_cc_v(&scc); if(!asm_opt.dbg_bam) destroy_cc_v(&scb); destroy_cc_v(&sca); + destroy_cc_v(&scc); /**if(!asm_opt.dbg_bam) destroy_cc_v(&scb);**/ destroy_cc_v(&sca); ha_print_ovlp_stat_0(b, n_thre, n_a); } else { ha_print_ovlp_stat_1(b, n_thre, n_a); - destroy_cc_v(&scc); if(!asm_opt.dbg_bam) destroy_cc_v(&scb); destroy_cc_v(&sca); + destroy_cc_v(&scc); /**if(!asm_opt.dbg_bam) destroy_cc_v(&scb);**/ destroy_cc_v(&sca); } destroy_ec_ovec_buf_t(b); diff --git a/gfa_ut.cpp b/gfa_ut.cpp index 54d5ab4..495c4fe 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -3053,13 +3053,14 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i // fprintf(stderr, "%.*s\tid::%u\tis_c::%u\n", // (int)Get_NAME_LENGTH(R_INF, 10819), Get_NAME(R_INF, 10819), 10819, is_contain_r((*rI), 10819)); // debug_info_of_specfic_node("m64011_190830_220126/47516220/ccs", sg, rI, "beg-0"); - // debug_info_of_specfic_node("bcb40bcc-d9cf-48e6-88ee-47ac3dde22ff", sg, rI, "beg-0"); + // debug_info_of_specfic_node("c7ecbd6b-e09d-4042-93ac-2400839feaf6", sg, rI, "beg-0"); if(asm_opt.is_ont) asg_arc_cut_weak(sg, &bu, max_tip, 0.975, 0, is_ou, 0, 1, 16, UL_COV_THRES-1, 0, rev, NULL, NULL); asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);///p_telo // fprintf(stderr, "[M::%s] count_edges_v_w(sg, 49778, 49847)->%ld\n", __func__, count_edges_v_w(sg, 49778, 49847)); // if(is_ou) dedup_contain_g(uopt, sg); + // debug_info_of_specfic_node("c7ecbd6b-e09d-4042-93ac-2400839feaf6", sg, rI, "beg-1"); for (i = 0; i < clean_round; i++, drop += step) { if(drop > max_ovlp_drop_ratio) drop = max_ovlp_drop_ratio; if(is_ou) { diff --git a/htab.cpp b/htab.cpp index 91da897..6061cfc 100644 --- a/htab.cpp +++ b/htab.cpp @@ -9,6 +9,7 @@ #include "ksort.h" #include "htab.h" #include "Process_Read.h" +#include "ecovlp.h" #define YAK_COUNTER_BITS 12 #define YAK_N_COUNTS (1<