From 286b897881fde283cff81e7088482a01763a318b Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Thu, 16 Feb 2023 00:04:06 -0500 Subject: [PATCH] r513 --- CommandLines.cpp | 4 ++++ CommandLines.h | 4 +++- Correct.cpp | 32 ++++++++++++++++++++++---------- Correct.h | 2 ++ Hash_Table.cpp | 1 + Levenshtein_distance.h | 2 ++ Output.h | 1 + Overlaps.cpp | 2 +- Overlaps.h | 2 ++ POA.h | 2 ++ Process_Read.h | 1 + Purge_Dups.h | 2 ++ gfa_ut.cpp | 10 ++++++++-- hic.h | 2 ++ horder.h | 2 ++ inter.cpp | 9 +++++---- ksw2.h | 1 + kthread.cpp | 1 + rcut.h | 2 ++ tovlp.h | 2 ++ 20 files changed, 66 insertions(+), 18 deletions(-) diff --git a/CommandLines.cpp b/CommandLines.cpp index 415f3c8..f002fe7 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -1,3 +1,4 @@ +#define __STDC_LIMIT_MACROS #include #include #include @@ -54,6 +55,7 @@ static ko_longopt_t long_options[] = { { "s-base", ko_required_argument, 340}, { "bin-only", ko_no_argument, 341}, { "ul-round", ko_required_argument, 342}, + { "prt-raw", ko_no_argument, 343}, { 0, 0, 0 } }; @@ -264,6 +266,7 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->is_bub_trans = 1; asm_opt->bin_only = 0; asm_opt->ul_clean_round = 1; + asm_opt->prt_dbg_gfa = 0; } void destory_enzyme(enzyme* f) @@ -797,6 +800,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) } else if (c == 341) asm_opt->bin_only = 1; else if (c == 342) asm_opt->ul_clean_round = atol(opt.arg); + else if (c == 343) asm_opt->prt_dbg_gfa = 1; else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg); } diff --git a/CommandLines.h b/CommandLines.h index eaef6ed..fb03367 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -1,10 +1,11 @@ #ifndef __COMMAND_LINE_PARSER__ #define __COMMAND_LINE_PARSER__ +#define __STDC_LIMIT_MACROS #include #include -#define HA_VERSION "0.18.6-r509" +#define HA_VERSION "0.18.6-r513" #define VERBOSE 0 @@ -138,6 +139,7 @@ typedef struct { uint8_t is_bub_trans; uint8_t bin_only; int32_t ul_clean_round; + int32_t prt_dbg_gfa; } hifiasm_opt_t; extern hifiasm_opt_t asm_opt; diff --git a/Correct.cpp b/Correct.cpp index 49e3814..9ee27d2 100644 --- a/Correct.cpp +++ b/Correct.cpp @@ -1,3 +1,4 @@ +#define __STDC_LIMIT_MACROS #include #include #include @@ -15796,7 +15797,8 @@ int64_t extract_sub_cigar_err_rr(overlap_region *z, int64_t s, int64_t e, ul_ov_ int64_t wk = ovlp_cur_wid(*p), xk = ovlp_cur_xoff(*p), ck = ovlp_cur_coff(*p); int64_t min_w = ovlp_min_wid(*p), max_w = ovlp_max_wid(*p);//[min_w, max_w] bit_extz_t ez; window_list *m; int64_t bd = ovlp_bd(*p), s0, e0; - s0 = z->w_list.a[min_w].x_start + bd; e0 = z->w_list.a[max_w].x_end+1-bd; + s0 = ((int64_t)(z->w_list.a[min_w].x_start)) + bd; + e0 = ((int64_t)(z->w_list.a[max_w].x_end))+1-bd; if(s < s0) s = s0; if(e > e0) e = e0;///exclude boundary if(s >= e) return -1; @@ -15811,7 +15813,7 @@ int64_t extract_sub_cigar_err_rr(overlap_region *z, int64_t s, int64_t e, ul_ov_ xk = z->w_list.a[wk].x_start; ck = 0; } - while(wk <= max_w && z->w_list.a[wk].x_start < e) {///[s, e) + while(wk <= max_w && z->w_list.a[wk].x_start < e) {///[s, e); [min_w, max_w] m = &(z->w_list.a[wk]); ws = m->x_start; we = m->x_end+1; os = MAX(s, ws); oe = MIN(e, we); @@ -15851,6 +15853,7 @@ int64_t extract_sub_cigar_err_rr(overlap_region *z, int64_t s, int64_t e, ul_ov_ return err; } +///[s, e) uint64_t gen_region_phase_robust_rr(overlap_region* ol, uint64_t *id_a, uint64_t id_n, uint64_t s, uint64_t e, uint64_t dp, ul_ov_t *c_idx, asg64_v *buf) { if(!id_n) return id_n; @@ -16412,7 +16415,7 @@ void rphase_rr(overlap_region_alloc* ol, const ul_idx_t *uref, const ug_opt_t *u m++; } z->w_list.n = m; - assert(zwn <= z->overlapLen); + assert(zwn <= z->overlapLen);///zwn is the length with secondary-best alignment z->align_length = z->overlapLen - zwn; } } @@ -16799,7 +16802,7 @@ int64_t get_pe_diff(overlap_region *q, uint64_t *qmask, uint64_t qmask_n, overla void retrieve_cigar_xcoord(bit_extz_t *ez, int64_t is, int64_t ie, int64_t *yk, int64_t *xk, int64_t *ck, int64_t *rxs, int64_t *rxe) { if(!ez->cigar.n) return; - int64_t cn = ez->cigar.n, op; int64_t ws, we, ws0, we0; + int64_t cn = ez->cigar.n, op; int64_t ws, we, ws0, we0, os, oe; if(((*ck) < 0) || ((*ck) > cn)) {//(*ck) == cn is allowed (*ck) = 0; (*yk) = ez->ps; (*xk) = ez->ts; } @@ -16818,9 +16821,17 @@ void retrieve_cigar_xcoord(bit_extz_t *ez, int64_t is, int64_t ie, int64_t *yk, if(op!=3) (*yk) += (ez->cigar.a[(*ck)]&(0x3fff)); if(op!=2) (*xk) += (ez->cigar.a[(*ck)]&(0x3fff)); we = (*yk); we0 = (*xk); - if(ws >= is && we <= ie) { + if(ws >= is && we <= ie) {///the cigar is fully contained if((*rxs) > ws0) (*rxs) = ws0; if((*rxe) < we0) (*rxe) = we0; + } else if(op == 0 || op == 1) {///overlap, it is hard to handle indels + os = MAX(is, ws); oe = MIN(ie, we); + if(oe > os) { + ws0 += (os - ws); + we0 = ws0 + (oe - os); + if((*rxs) > ws0) (*rxs) = ws0; + if((*rxe) < we0) (*rxe) = we0; + } } (*ck)++; } @@ -16904,7 +16915,7 @@ uint64_t extract_xcoordates(overlap_region *z, uint64_t *a, int64_t a_n, asg64_v ovlp_cur_xoff(m) = z->w_list.a[w[0]].y_start; ///cur ypos ovlp_cur_coff(m) = 0; ///cur cigar off in cur window - is = t[0]; ie = t[1]+1; + is = t[0]; ie = t[1]+1; ///[is, ie) for (a_z = a_i; a_z >= 0; a_z--) { as = a[a_z]>>32; ae = (uint32_t)a[a_z]; if(ae <= is) break; @@ -17067,7 +17078,7 @@ void update_masks(asg64_v *in, overlap_region *qi, overlap_region *ti, int64_t b } in_n1 = in->n; - radix_sort_bc64(in->a+in_n0, in->a+in_n1); + radix_sort_bc64(in->a+in_n0, in->a+in_n1);///(uint32_t)in->a[]: 0-> original; 1-> mask for (k = in_n0, dp = old_dp = dp_mask = old_dp_mask = 0, start = 0, end = -1; k < in_n1; ++k) {///[beg, end) but coordinates in idx is [, ] ///if idx->a.a[] is qe old_dp = dp; old_dp_mask = dp_mask; @@ -17150,7 +17161,8 @@ asg64_v* buf1, int64_t bd, ul_ov_t *res) { ul_ov_t *ref, r0, r1; uint32_t bn = buf->n, s, e; uint64_t is_exact = 0; overlap_region *q = &(a[qi]), *t = &(a[ti]); int64_t dis; - assert(q->y_id < t->y_id); assert(ovdb_idx < ov_db->n); + assert(q->y_id < t->y_id); + assert(ovdb_idx < ov_db->n);///quickly jump to the related qn-tn pair ov_db[] cal_x_ul_ovlp(ug, q, t, &r0); r1 = r0; r1.qn = r0.tn; r1.tn = r0.qn; @@ -17310,12 +17322,12 @@ int64_t rlen, mask_ul_ov_t *mk, idx_emask_t *mm, double len_diff) ///sort overlaps by yid for filtering kv_push(uint64_t, *idx, (((uint64_t)z->y_id)<<32)|((uint64_t)k)); } - radix_sort_bc64(idx->a, idx->a+idx->n); + radix_sort_bc64(idx->a, idx->a+idx->n); ///for filtering radix_sort_bc64(buf->a, buf->a+buf->n);///sort by error; smaller first for (l = 0, k = 1; k <= buf->n; k++) { if(k == buf->n || (buf->a[l]>>32) == (buf->a[k]>>32)) {///with equal number of normalized errors - if((k - l) >1) { + if((k - l) > 1) { for (i = l; i < k; i++) { z = &(ol->list[(uint32_t)buf->a[i]]); m0 = z->x_pos_e+1-z->x_pos_s; diff --git a/Correct.h b/Correct.h index c82755d..469ec2b 100644 --- a/Correct.h +++ b/Correct.h @@ -1,5 +1,7 @@ #ifndef __CORRECT__ #define __CORRECT__ + +#define __STDC_LIMIT_MACROS #include #include "Hash_Table.h" #include "Levenshtein_distance.h" diff --git a/Hash_Table.cpp b/Hash_Table.cpp index 796afad..0626305 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -1,3 +1,4 @@ +#define __STDC_LIMIT_MACROS #include #include #include diff --git a/Levenshtein_distance.h b/Levenshtein_distance.h index ad370fc..592f1b5 100644 --- a/Levenshtein_distance.h +++ b/Levenshtein_distance.h @@ -1,5 +1,7 @@ #ifndef __LEVENSHTEIN__ #define __LEVENSHTEIN__ + +#define __STDC_LIMIT_MACROS #include #include "emmintrin.h" #include "nmmintrin.h" diff --git a/Output.h b/Output.h index 9f95698..46271b3 100644 --- a/Output.h +++ b/Output.h @@ -1,6 +1,7 @@ #ifndef __OUTPUT__ #define __OUTPUT__ +#define __STDC_LIMIT_MACROS #include #include #include diff --git a/Overlaps.cpp b/Overlaps.cpp index 1b72381..6e2b4f1 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -33703,7 +33703,7 @@ bub_label_t *b_mask_t, uint32_t is_trio, int32_t ul_aln_round, char *o_file, con int32_t k, strl = strlen(bin_file)+1, kt, cl, sl; char *id = NULL; renew_g(src, rev_src, n_read, readLen, cov, ruIndex, sg, mini_overlap_length, max_hang_length, uopt, clean_round, min_ovlp_drop_ratio, max_ovlp_drop_ratio, asm_opt.max_short_tip, b_mask_t, - is_trio, o_file, bin_file, (ul_aln_round<=1)?1:0, 1, 1); + is_trio, o_file, bin_file, (ul_aln_round<=1)?1:0, 1, ((is_trio)?(0):(1))/**1**/); gen_ug_opt_t(uopt, *src, *rev_src, max_hang_length, mini_overlap_length, gap_fuzz, min_dp, *readLen, *cov, ruIndex, (asm_opt.max_short_tip*2), 0.15, 3, 0.05, 0.9, b_mask_t); ug_ext_gfa(uopt, *sg, ug_ext_len); diff --git a/Overlaps.h b/Overlaps.h index cef917e..21d7482 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -1,5 +1,7 @@ #ifndef __OVERLAPS__ #define __OVERLAPS__ + +#define __STDC_LIMIT_MACROS #include #include #include "kvec.h" diff --git a/POA.h b/POA.h index 00656f6..a19ed1c 100644 --- a/POA.h +++ b/POA.h @@ -1,5 +1,7 @@ #ifndef __POA_PARSER__ #define __POA_PARSER__ + +#define __STDC_LIMIT_MACROS #include #include "Hash_Table.h" #include "Process_Read.h" diff --git a/Process_Read.h b/Process_Read.h index fb56d5d..88362af 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -1,6 +1,7 @@ #ifndef __READ__ #define __READ__ +#define __STDC_LIMIT_MACROS #include #include #include diff --git a/Purge_Dups.h b/Purge_Dups.h index 018ab74..ce0d18a 100644 --- a/Purge_Dups.h +++ b/Purge_Dups.h @@ -1,5 +1,7 @@ #ifndef __PURGEDUPS__ #define __PURGEDUPS__ + +#define __STDC_LIMIT_MACROS #include #include #include "kvec.h" diff --git a/gfa_ut.cpp b/gfa_ut.cpp index f06235f..13e5ae4 100644 --- a/gfa_ut.cpp +++ b/gfa_ut.cpp @@ -2568,6 +2568,13 @@ void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc uopt->min_ovlp, 10, uopt->gap_fuzz, b_mask_t, no_trio_recover); } +void prt_dbg_gfa(ug_opt_t *uopt, asg_t *sg, char *o_file, const char *suffix) +{ + char* gfa_name; MALLOC(gfa_name, strlen(o_file)+strlen(suffix)+50); sprintf(gfa_name, "%s.%s", o_file, suffix); + print_debug_gfa(sg, NULL, uopt->coverage_cut, gfa_name, uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 0, 0, 1); + free(gfa_name); +} + void ul_clean_gfa(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, int64_t clean_round, double min_ovlp_drop_ratio, double max_ovlp_drop_ratio, double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, int32_t is_ou, int32_t is_trio, uint32_t ou_thres, char *o_file) { @@ -2585,8 +2592,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i // exit(1); // print_raw_u2rgfa_seq(&UL_INF, rI, 1); // exit(1); - - + if(asm_opt.prt_dbg_gfa) prt_dbg_gfa(uopt, sg, o_file, "raw"); // fprintf(stderr, "%.*s\tid::%u\tis_c::%u\n", // (int)Get_NAME_LENGTH(R_INF, 10785), Get_NAME(R_INF, 10785), 10785, is_contain_r((*rI), 10785)); // fprintf(stderr, "%.*s\tid::%u\tis_c::%u\n", diff --git a/hic.h b/hic.h index fc456e2..1c3226f 100644 --- a/hic.h +++ b/hic.h @@ -1,5 +1,7 @@ #ifndef __HIC__ #define __HIC__ + +#define __STDC_LIMIT_MACROS #include #include "Overlaps.h" diff --git a/horder.h b/horder.h index 8e136a9..8221b9f 100644 --- a/horder.h +++ b/horder.h @@ -1,5 +1,7 @@ #ifndef __HORDER__ #define __HORDER__ + +#define __STDC_LIMIT_MACROS #include #include "hic.h" diff --git a/inter.cpp b/inter.cpp index 9317e60..1fbda1e 100644 --- a/inter.cpp +++ b/inter.cpp @@ -17753,7 +17753,7 @@ void gen_trans_base_count_comp(ug_trans_t *p, kv_u_trans_t *res) double index_time = yak_realtime(); // ha_flt_tab = NULL; uint64_t i, k, l, occ, m, cc; kv_ul_ov_t *bl = NULL; - u_trans_t *z; ha_mzl_t *tz; + u_trans_t *z; ha_mzl_t *tz; double ww; clean_u_trans_t_idx_adv(res, p->ug, p->rg); p->filter = res; // fprintf(stderr, "[M::%s::] ==> 0\n", __func__); p->is_cnt = 1; p->is_ovlp = 0; @@ -17834,7 +17834,7 @@ void gen_trans_base_count_comp(ug_trans_t *p, kv_u_trans_t *res) } } // fprintf(stderr, "[M::%s::] ==> 5\n", __func__); - kt_for(p->n_thread, worker_for_sysm_trans_ovlp, p, p->ug->u.n); + // kt_for(p->n_thread, worker_for_sysm_trans_ovlp, p, p->ug->u.n);///not correct // assert(p->srt_a.n <= p->ug->u.n); // radix_sort_ha_mzl_t_srt(p->srt_a.a, p->srt_a.a + p->srt_a.n); kv_resize(u_trans_t, *res, occ); @@ -17846,13 +17846,14 @@ void gen_trans_base_count_comp(ug_trans_t *p, kv_u_trans_t *res) assert(bl->a[k].qn == (tz->x>>32)); for (; (k < bl->n) && (bl->a[k].qn == (tz->x>>32)); k++) { if(bl->a[k].qn == bl->a[k].tn) continue; - if(bl->a[k].qs == (uint32_t)-1 && bl->a[k].qe == (uint32_t)-1) continue; + ww = cal_trans_ov_w(&(bl->a[k])); + if(ww <= 0) continue; kv_pushp(u_trans_t, *res, &z); z->f = RC_3; z->rev = bl->a[k].rev; z->del = 0; z->qn = bl->a[k].qn; z->qs = bl->a[k].qs; z->qe = bl->a[k].qe; z->tn = bl->a[k].tn; z->ts = bl->a[k].ts; z->te = bl->a[k].te; - z->nw = cal_trans_ov_w(&(bl->a[k])); assert(z->nw > 0); + z->nw = ww; // if(z->qn == 56 || z->qn == 160 || z->tn == 56 || z->tn == 160) { // fprintf(stderr, ">>>utg%.6u%c\t%u\t%u\t%u\t%c\tutg%.6u%c\t%u\t%u\t%u\tnw::%f\n", // z->qn+1, "lc"[p->ug->u.a[z->qn].circ], p->ug->u.a[z->qn].len, z->qs, z->qe, "+-"[z->rev], diff --git a/ksw2.h b/ksw2.h index 213c27f..d8465aa 100644 --- a/ksw2.h +++ b/ksw2.h @@ -1,6 +1,7 @@ #ifndef KSW2_H_ #define KSW2_H_ +#define __STDC_LIMIT_MACROS #include #define KSW_NEG_INF -0x40000000 diff --git a/kthread.cpp b/kthread.cpp index ffdf940..d79d72c 100644 --- a/kthread.cpp +++ b/kthread.cpp @@ -1,3 +1,4 @@ +#define __STDC_LIMIT_MACROS #include #include #include diff --git a/rcut.h b/rcut.h index e2baf47..2679e3c 100644 --- a/rcut.h +++ b/rcut.h @@ -1,5 +1,7 @@ #ifndef __RCUT__ #define __RCUT__ + +#define __STDC_LIMIT_MACROS #include #include #include "kvec.h" diff --git a/tovlp.h b/tovlp.h index b2874ef..f009fe4 100644 --- a/tovlp.h +++ b/tovlp.h @@ -1,5 +1,7 @@ #ifndef __TOVLP__ #define __TOVLP__ + +#define __STDC_LIMIT_MACROS #include #include "Overlaps.h"