r513

CommandLines.cpp

@@ -1,3 +1,4 @@
+#define __STDC_LIMIT_MACROS
 #include <zlib.h>
 #include <stdlib.h>
 #include <stdio.h>
@@ -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);
         }

CommandLines.h

@@ -1,10 +1,11 @@
 #ifndef __COMMAND_LINE_PARSER__
 #define __COMMAND_LINE_PARSER__
 
+#define __STDC_LIMIT_MACROS
 #include <pthread.h>
 #include <stdint.h>
 
-#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;

Correct.cpp

@@ -1,3 +1,4 @@
+#define __STDC_LIMIT_MACROS
 #include <stdio.h>
 #include <stdlib.h>
 #include <stdint.h>
@@ -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;

Correct.h

@@ -1,5 +1,7 @@
 #ifndef __CORRECT__
 #define __CORRECT__
+
+#define __STDC_LIMIT_MACROS
 #include <stdint.h>
 #include "Hash_Table.h"
 #include "Levenshtein_distance.h"

Hash_Table.cpp

@@ -1,3 +1,4 @@
+#define __STDC_LIMIT_MACROS
 #include <stdio.h>
 #include <stdlib.h>
 #include <stdint.h>

Levenshtein_distance.h

@@ -1,5 +1,7 @@
 #ifndef __LEVENSHTEIN__
 #define __LEVENSHTEIN__
+
+#define __STDC_LIMIT_MACROS
 #include <stdint.h>
 #include "emmintrin.h"
 #include "nmmintrin.h"

Output.h

@@ -1,6 +1,7 @@
 #ifndef __OUTPUT__
 #define __OUTPUT__
 
+#define __STDC_LIMIT_MACROS
 #include<stdint.h>
 #include <string.h>
 #include <stdlib.h>

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);

Overlaps.h

@@ -1,5 +1,7 @@
 #ifndef __OVERLAPS__
 #define __OVERLAPS__
+
+#define __STDC_LIMIT_MACROS
 #include <stdio.h>
 #include <stdint.h>
 #include "kvec.h"

POA.h

@@ -1,5 +1,7 @@
 #ifndef __POA_PARSER__
 #define __POA_PARSER__
+
+#define __STDC_LIMIT_MACROS
 #include <stdint.h>
 #include "Hash_Table.h"
 #include "Process_Read.h"

Process_Read.h

@@ -1,6 +1,7 @@
 #ifndef __READ__
 #define __READ__
 
+#define __STDC_LIMIT_MACROS
 #include<stdint.h>
 #include <string.h>
 #include <stdlib.h>

Purge_Dups.h

@@ -1,5 +1,7 @@
 #ifndef __PURGEDUPS__
 #define __PURGEDUPS__
+
+#define __STDC_LIMIT_MACROS
 #include <stdio.h>
 #include <stdint.h>
 #include "kvec.h"

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", 

hic.h

@@ -1,5 +1,7 @@
 #ifndef __HIC__
 #define __HIC__
+
+#define __STDC_LIMIT_MACROS
 #include <stdint.h>
 #include "Overlaps.h"
 

horder.h

@@ -1,5 +1,7 @@
 #ifndef __HORDER__
 #define __HORDER__
+
+#define __STDC_LIMIT_MACROS
 #include <stdint.h>
 #include "hic.h"
 

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],

ksw2.h

@@ -1,6 +1,7 @@
 #ifndef KSW2_H_
 #define KSW2_H_
 
+#define __STDC_LIMIT_MACROS
 #include <stdint.h>
 
 #define KSW_NEG_INF -0x40000000

kthread.cpp

@@ -1,3 +1,4 @@
+#define __STDC_LIMIT_MACROS
 #include <pthread.h>
 #include <stdlib.h>
 #include <limits.h>

rcut.h

@@ -1,5 +1,7 @@
 #ifndef __RCUT__
 #define __RCUT__
+
+#define __STDC_LIMIT_MACROS
 #include <stdio.h>
 #include <stdint.h>
 #include "kvec.h"

tovlp.h

@@ -1,5 +1,7 @@
 #ifndef __TOVLP__
 #define __TOVLP__
+
+#define __STDC_LIMIT_MACROS
 #include <stdint.h>
 #include "Overlaps.h"