Merge pull request #553 from chhylp123/hifiasm_dev_debug

fix bubble issue

CommandLines.cpp

@@ -65,6 +65,7 @@ static ko_longopt_t long_options[] = {
     { "dual-scaf",     ko_no_argument, 350},
     { "scaf-gap",     ko_required_argument, 351},
     { "sec-in",     ko_required_argument, 352},
+    { "somatic-cov",     ko_required_argument, 353},
     // { "path-round",     ko_required_argument, 348},
 	{ 0, 0, 0 }
 };
@@ -307,6 +308,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->self_scaf_reliable_min = 5000000;
     asm_opt->self_scaf_gap_max = 3000000;
     asm_opt->sec_in = NULL;
+    asm_opt->somatic_cov = -1;
 }
 
 void destory_enzyme(enzyme* f)
@@ -856,6 +858,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 350) asm_opt->self_scaf = 1;
         else if (c == 351) asm_opt->self_scaf_gap_max = atol(opt.arg);
         else if (c == 352) get_hic_enzymes(opt.arg, &(asm_opt->sec_in), 0);
+        else if (c == 353) asm_opt->somatic_cov = atol(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);
         }

CommandLines.h

@@ -5,7 +5,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.19.8-r602"
+#define HA_VERSION "0.19.8-r603"
 
 #define VERBOSE 0
 
@@ -152,6 +152,7 @@ typedef struct {
     uint64_t self_scaf_min;
     uint64_t self_scaf_reliable_min;
     int64_t self_scaf_gap_max;
+    int64_t somatic_cov;
 } hifiasm_opt_t;
 
 extern hifiasm_opt_t asm_opt;

Overlaps.cpp

@@ -58,6 +58,9 @@ KRADIX_SORT_INIT(u_trans_ts, u_trans_t, u_trans_ts_key, member_size(u_trans_t, t
 #define ha_mzl_t_key(p) ((p).x)
 KRADIX_SORT_INIT(ha_mzl_t_srt1, ha_mzl_t, ha_mzl_t_key, member_size(ha_mzl_t, x))
 
+#define Uc_beg(z) ((uint32_t)((z).a[0]>>32))
+#define Uc_end(z) ((uint32_t)((z).a[(z).n-1]>>32)^1)
+
 #define UL_COV_THRES 2
 #define PHASE_SEP 64
 #define PHASE_SEF 2
@@ -19409,10 +19412,11 @@ ma_ug_t* gen_fg(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc* src, ma_sub_t *cov, int3
     uint64_t i, k, l, m, rv, rw, uv, uw, zn, z, nist = 0; ma_utg_t *u;
     for (k = 0; k < fg->u.n; k++) {
         u = &(ug->u.a[k]); fg->g->seq[k].c = PRIMARY_LABLE;
-        if(u->circ) continue;
-        m = k<<1; m |= (((uint64_t)u->start)<<32); kv_push(uint64_t, srt, m);
-        m = (k<<1)+1; m |= (((uint64_t)u->end)<<32); kv_push(uint64_t, srt, m);
+        if((u->circ) || (u->n == 0)) continue;
+        m = k<<1; m |= (((uint64_t)Uc_beg((*u)))<<32); kv_push(uint64_t, srt, m);///((uint32_t)(((*u)).a[0]>>32));
+        m = (k<<1)+1; m |= (((uint64_t)Uc_end((*u)))<<32); kv_push(uint64_t, srt, m);
     }
+
 
     radix_sort_arch64(srt.a, srt.a+srt.n);
     for (k = 1, l = 0; k <= srt.n; k++) {
@@ -19427,9 +19431,9 @@ ma_ug_t* gen_fg(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc* src, ma_sub_t *cov, int3
     int32_t r; asg_arc_t t0, t1, *p;
     for (k = 0; k < fg->u.n; k++) {
         u = &(ug->u.a[k]);
-        if(u->circ) continue;
+        if((u->circ) || (u->n == 0)) continue;
 
-        uv = k<<1; rv = u->end^1;
+        uv = k<<1; rv = Uc_end((*u))^1;
         x = &(src[rv>>1]); 
         za = asg_arc_a(ug->g, uv); 
         zn = asg_arc_n(ug->g, uv); 
@@ -19446,6 +19450,8 @@ ma_ug_t* gen_fg(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc* src, ma_sub_t *cov, int3
             if((t0.ul>>32) != rv) continue;
             rw = t0.v;
             if(idx[rw>>1] == ((uint32_t)-1)) continue;
+            if(!(get_edge_from_source(src, cov, NULL, max_hang, min_ovlp, t0.ul>>32, t0.v, &t0))) continue;
+
             m = idx[rw>>1]; assert((srt.a[m]>>33) == (rw>>1));
             for (; m < srt.n && (srt.a[m]>>33) == (rw>>1); m++) {
                 if(rw == (srt.a[m]>>32)) {
@@ -19468,7 +19474,7 @@ ma_ug_t* gen_fg(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc* src, ma_sub_t *cov, int3
         }
 
 
-        uv = (k<<1)+1; rv = u->start^1;
+        uv = (k<<1)+1; rv = Uc_beg((*u))^1;
         x = &(src[rv>>1]); 
         za = asg_arc_a(ug->g, uv); 
         zn = asg_arc_n(ug->g, uv); 
@@ -19485,6 +19491,8 @@ ma_ug_t* gen_fg(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc* src, ma_sub_t *cov, int3
             if((t0.ul>>32) != rv) continue;
             rw = t0.v;
             if(idx[rw>>1] == ((uint32_t)-1)) continue;
+            if(!(get_edge_from_source(src, cov, NULL, max_hang, min_ovlp, t0.ul>>32, t0.v, &t0))) continue;
+
             m = idx[rw>>1]; assert((srt.a[m]>>33) == (rw>>1));
             for (; m < srt.n && (srt.a[m]>>33) == (rw>>1); m++) {
                 if(rw == (srt.a[m]>>32)) {
@@ -19522,6 +19530,7 @@ ma_ug_t* gen_fg(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc* src, ma_sub_t *cov, int3
     return fg;
 }
 
+
 rd_hamming_fly_simp_t* gen_rd_hamming_fly_simp_t(ma_ug_t *ug, asg_t *rg, ma_hit_t_alloc* src, ma_sub_t *cov, int32_t max_hang, int32_t min_ovlp, int32_t gap_fuzz, kvec_asg_arc_t_warp *ae)
 {
     rd_hamming_fly_simp_t *p; CALLOC(p, 1);
@@ -25997,8 +26006,9 @@ static inline void asg_arc_rest(asg_t *des, asg_t *src, uint32_t v0, uint32_t w0
     for (i = 0; i < nv && av[i].ul < arc->ul; ++i);
     if(i >= nv) i = nv - 1; av[i] = *arc;
 
-    rv = ((v0&1)?(ug->u.a[v0>>1].start^1):(ug->u.a[v0>>1].end^1));
-    rw = ((w0&1)?(ug->u.a[w0>>1].end):(ug->u.a[w0>>1].start));
+
+    rv = ((v0&1)?(Uc_beg(ug->u.a[v0>>1])^1):(Uc_end(ug->u.a[v0>>1])^1));
+    rw = ((w0&1)?(Uc_end(ug->u.a[w0>>1])):(Uc_beg(ug->u.a[w0>>1])));
     assert(get_edge_from_source(src_e, cov, NULL, max_hang, min_ovlp, rv, rw, &t));
     kv_push(asg_arc_t, ae->a, t);
 
@@ -26126,7 +26136,8 @@ uint32_t positive_flag, uint32_t negative_flag, hap_cov_t *cov, utg_trans_t *o,
         for (k = 0; k < an; k++) {
             ///s[av[k].v]&se:: in the existing graph
             if((!s[av[k].v]) || (s[av[k].v]&se)) continue; ///in the existing graph
-            if((av[k].v) == (v>>1)) continue;
+            // if((av[k].v) == (v>>1)) continue;///looks like a bug
+            if((av[k].v>>1) == (v>>1)) continue;
             av[k].del = 0;
             ra = asg_arc_a(fg->g, (av[k].v^1)); rn = asg_arc_n(fg->g, (av[k].v^1));
             for (ri = 0; ri < rn; ri++) {
@@ -37839,15 +37850,17 @@ void flat_bubbles_advance(asg_t *sg, ma_hit_t_alloc* sources, R_to_U* ruIndex, u
 void flat_soma_v(asg_t *sg, ma_hit_t_alloc* sources, R_to_U* ruIndex)
 {
     uint64_t dip_thre_max;
-    if(asm_opt.hom_global_coverage_set)
-    {
-        dip_thre_max = asm_opt.hom_global_coverage;
-    }
-    else
-    {
-        dip_thre_max = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE);
+    if(asm_opt.somatic_cov >= 0) {
+        dip_thre_max = asm_opt.somatic_cov;
+    } else {
+        if(asm_opt.hom_global_coverage_set) {
+            dip_thre_max = asm_opt.hom_global_coverage;
+        } else {
+            dip_thre_max = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE);
+        }
+        dip_thre_max = (((double)(dip_thre_max)*1.15)/asm_opt.polyploidy);
     }
-    flat_bubbles_advance(sg, sources, ruIndex, (((double)(dip_thre_max)*1.15)/asm_opt.polyploidy));
+    flat_bubbles_advance(sg, sources, ruIndex, dip_thre_max);
 }
 
 char *get_outfile_name(char* output_file_name)