update for low coverge data

CommandLines.h

@@ -5,7 +5,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.23.0-r691"
+#define HA_VERSION "0.24.0-r702"
 
 #define VERBOSE 0
 

Correct.h

@@ -1415,4 +1415,6 @@ void get_wqual(uint64_t zid, uint64_t zpos, uint64_t zrev, asg8_v *v, uint8_t *v
 #define HPC_RR 4
 #define HPC_CC 2
 
+// #define FORCE_CUT 1
+
 #endif

Overlaps.h

@@ -632,7 +632,7 @@ static inline int count_out_without_del(const asg_t *g, uint32_t v)
 
 void build_string_graph_without_clean(
 int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
-long long n_read, uint64_t* readLen, long long mini_overlap_length, 
+uint64_t n_read, uint64_t* readLen, long long mini_overlap_length, 
 long long max_hang_length, long long clean_round, long long gap_fuzz,
 float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, 
 long long bubble_dist, int read_graph, int write);
@@ -927,7 +927,7 @@ ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, int m
 void rescue_bubble_by_chain(asg_t *sg, ma_sub_t *coverage_cut, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, 
 long long tipsLen, float tip_drop_ratio, long long stops_threshold, R_to_U* ruIndex, 
 float chimeric_rate, float drop_ratio, int max_hang, int min_ovlp, uint32_t chainLenThres, long long gap_fuzz, 
-bub_label_t* b_mask_t, long long no_trio_recover);
+bub_label_t* b_mask_t, long long no_trio_recover, uint8_t *cmk);
 
 typedef struct{
     double weight;
@@ -1242,6 +1242,9 @@ uint64_t infer_mmhap_copy(ma_ug_t *ug, asg_t *sg, ma_hit_t_alloc *src, uint8_t *
 uint64_t trans_sec_cut0(kv_u_trans_t *ta, asg64_v *srt, uint32_t id, double sec_rate, uint64_t bd, ma_ug_t *ug);
 void clean_u_trans_t_idx_filter_mmhap_adv(kv_u_trans_t *ta, ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* src, ug_rid_cov_t *in);
 void gen_ug_rid_cov_t_by_ovlp(kv_u_trans_t *ta, ug_rid_cov_t *cc);
+void rescue_chimeric_reads_aggressive(ma_ug_t *i_ug, asg_t *rg, ma_hit_t_alloc* sources, ma_sub_t *coverage_cut,
+R_to_U* ruIndex, int max_hang, int min_ovlp, uint32_t chainLenThres, uint32_t is_bubble_check, uint32_t is_primary_check, kvec_asg_arc_t_warp* new_rtg_edges, 
+kvec_t_u32_warp* new_rtg_nodes, bub_label_t* b_mask_t, uint8_t *cmk);
 
 #define UC_Read_resize(v, s) do {\
 		if ((v).size<(s)) {REALLOC((v).seq,(s));(v).size=(s);}\

ecovlp.cpp

@@ -60,6 +60,7 @@ typedef struct {
 typedef struct {
 	ec_ovec_buf_t0 *a;
 	uint32_t n, rev;
+    uint8_t *cr;
 } ec_ovec_buf_t;
 
 ec_ovec_buf_t* gen_ec_ovec_buf_t(uint32_t n);
@@ -180,7 +181,7 @@ void destroy_ec_ovec_buf_t(ec_ovec_buf_t *p)
         destroy_cns_gfa(&(z->cns));
 
     }
-    free(p->a); free(p);
+    free(p->a); free(p->cr); free(p);
 
     // fprintf(stderr, "[M::%s-chains] #->%lld\n", __func__, asm_opt.num_bases);
     // fprintf(stderr, "[M::%s-passed-chains-0] #->%lld\n", __func__, asm_opt.num_corrected_bases);
@@ -3931,10 +3932,11 @@ uint32_t is_chemical_r(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_t co
 }
 
 
-uint32_t is_chemical_r_adv(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_t cov, int64_t cut_len, double dup_rate)
+uint32_t is_chemical_r_adv(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_t cov, int64_t cut_len, double dup_rate, uint64_t is_del)
 {
     uint64_t k, s, e; int64_t dp, old_dp, st = 0, ed, s0, e0, rr, lt;
     for (k = idx->n = 0; k < ov->length; k++) {
+        if(is_del && ov->buffer[k].del) continue;
         s0 = (uint32_t)ov->buffer[k].qns; e0 = ov->buffer[k].qe;
         if(s0 > 0) s0 += cut_len;
         if(e0 < len) e0 -= cut_len;
@@ -3988,6 +3990,64 @@ uint32_t is_chemical_r_adv(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_
     return 0;
 }
 
+int64_t cal_chemical_r_adv(ma_hit_t_alloc *ov, asg64_v *idx, int64_t len, int64_t cut_len, double dup_rate, uint64_t is_del)
+{
+    uint64_t k, s, e; int64_t dp, old_dp, st = 0, ed, s0, e0, rr, lt, min_cov;
+    for (k = idx->n = 0; k < ov->length; k++) {
+        if(is_del && ov->buffer[k].del) continue;
+        s0 = (uint32_t)ov->buffer[k].qns; e0 = ov->buffer[k].qe;
+        if(s0 > 0) s0 += cut_len;
+        if(e0 < len) e0 -= cut_len;
+        if(e0 <= s0) continue;
+        s = s0; e = e0;
+
+        lt = Get_READ_LENGTH((R_INF), ov->buffer[k].tn);
+        rr = (lt >= len)?(lt - len):(len - lt);
+        if((rr <= (len*dup_rate)) && (rr <= (lt*dup_rate)) && (ov->buffer[k].rev)) {
+            dp = (ov->buffer[k].qe) - ((uint32_t)ov->buffer[k].qns); dp = len - dp;
+            old_dp = ov->buffer[k].te - ov->buffer[k].ts; old_dp = lt - old_dp;
+            if((dp <= (len*dup_rate)) && (old_dp <= (lt*dup_rate))) continue;
+        } 
+
+        kv_push(uint64_t, (*idx), (s<<1));
+        kv_push(uint64_t, (*idx), (e<<1)|1);
+    }
+
+    radix_sort_ec64(idx->a, idx->a + idx->n); s0 = e0 = rr = -1; min_cov = INT64_MAX;
+    for (k = 0, dp = 0, st = ed = 0; k < idx->n; ++k) {
+        old_dp = dp;
+        ///if a[j] is qe
+        if (idx->a[k]&1) --dp;
+        else ++dp;
+
+        ed = idx->a[k]>>1;
+        if(ed > st) {
+            // if(ov->length && ((ov->buffer[0].qns>>32) == 5045637)) {
+            //     fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\tid::%lu\n", __func__, st, ed, old_dp, len, ov->buffer[0].qns>>32);
+            // }
+            if(old_dp <= min_cov) {
+                // if(ov->length && (ov->buffer[0].qns>>32) == 22344) fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\n", __func__, st, ed, old_dp, len);
+                min_cov = old_dp;
+            } 
+        }
+        st = ed;
+    }
+
+
+    ed = len; old_dp = dp;
+    if(ed > st) {
+        // if(ov->length && ((ov->buffer[0].qns>>32) == 5045637)) {
+        //     fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\tid::%lu\n", __func__, st, ed, old_dp, len, ov->buffer[0].qns>>32);
+        // }
+        if(old_dp <= min_cov) {
+            // if(ov->length && (ov->buffer[0].qns>>32) == 22344) fprintf(stderr, "[M::%s]\tmd::[%ld,%ld)\tcov::%ld\tlen::%ld\n", __func__, st, ed, old_dp, len);
+            min_cov = old_dp;
+        }
+    }
+    
+    return min_cov;
+}
+
 void prt_dbg_rid_paf(ma_hit_t_alloc *ov, UC_Read *ra, asg8_v *qa)
 {
     if(!(ov->length)) return;
@@ -4051,7 +4111,7 @@ static void worker_hap_dc_ec_chemical_r(void *data, long i, int tid)
     if(b->cnt[1] == 0) {
         // if(i == 6204620) prt_dbg_rid_paf(&(R_INF.paf[i]), &(b->self_read), &(b->v8q));
         // if(is_chemical_r(&(R_INF.paf[i]), &b->v64, Get_READ_LENGTH((R_INF), i), 3, 16)) {
-        if(is_chemical_r_adv(&(R_INF.paf[i]), &b->v64, Get_READ_LENGTH((R_INF), i), asm_opt.chemical_cov, asm_opt.chemical_flank, 0.02)) {
+        if(is_chemical_r_adv(&(R_INF.paf[i]), &b->v64, Get_READ_LENGTH((R_INF), i), asm_opt.chemical_cov, asm_opt.chemical_flank, 0.02, 0)) {
             // fprintf(stderr, "-um-[M::%s]\tqn::%u::%.*s\n\n", __func__, (uint32_t)(i), (int)Get_NAME_LENGTH(R_INF, i), Get_NAME((R_INF), i));
             R_INF.paf[i].length = 0; b->cnt[0]++;
         } 
@@ -4075,11 +4135,56 @@ static void worker_hap_dc_ec_chemical_r(void *data, long i, int tid)
 static void worker_hap_dc_ec_chemical_arc(void *data, long i, int tid)
 {
     ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
-    ma_hit_t_alloc *paf = &(R_INF.paf[i]); uint64_t k;
+    ma_hit_t_alloc *paf = &(R_INF.paf[i]), *rev; uint64_t k, z;
 
-    if(is_chemical_r_adv(&(R_INF.paf[i]), &b->v64, Get_READ_LENGTH((R_INF), i), asm_opt.chemical_cov, asm_opt.chemical_flank, 0.02)) {
-        // fprintf(stderr, "-um-[M::%s]\tqn::%u::%.*s\n\n", __func__, (uint32_t)(i), (int)Get_NAME_LENGTH(R_INF, i), Get_NAME((R_INF), i));   
-        for (k = 0; k < paf->length; k++) paf->buffer[k].del = 1; b->cnt[0]++;
+    if(b->cnt[1] == 0) {
+        if(is_chemical_r_adv(&(R_INF.paf[i]), &b->v64, Get_READ_LENGTH((R_INF), i), asm_opt.chemical_cov, asm_opt.chemical_flank, 0.02, 1)) {
+            // fprintf(stderr, "-um-[M::%s]\tqn::%u::%.*s\n\n", __func__, (uint32_t)(i), (int)Get_NAME_LENGTH(R_INF, i), Get_NAME((R_INF), i));   
+            for (k = 0; k < paf->length; k++) paf->buffer[k].del = 1; b->cnt[0]++;
+        }
+    } else if(b->cnt[1] == 1) {
+        for (k = 0; k < paf->length; k++) {
+            if((Get_qn(paf->buffer[k])) > (Get_tn(paf->buffer[k]))) continue;
+            rev = &(R_INF.paf[paf->buffer[k].tn]);
+            for (z = 0; z < rev->length; z++) {
+                if((rev->buffer[z].tn == (Get_qn(paf->buffer[k])))) {
+                    if(paf->buffer[k].del != rev->buffer[z].del) {
+                        paf->buffer[k].del = rev->buffer[z].del = 1;
+                    }
+                }
+            }
+        }
+    }
+
+   refresh_ec_ovec_buf_t0(b, REFRESH_N);
+}
+
+static void worker_hap_dc_ec_chemical_arc_mark(void *data, long i, int tid)
+{
+    ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
+    ma_hit_t_alloc *paf = &(R_INF.paf[i]), *rev; uint64_t k, z; int64_t cov, msk_cut = asm_opt.chemical_cov;
+    uint8_t *msk = ((ec_ovec_buf_t*)data)->cr;
+
+    if(b->cnt[1] == 0) {
+        msk[i] = (uint8_t)-1;
+        cov = cal_chemical_r_adv(&(R_INF.paf[i]), &b->v64, Get_READ_LENGTH((R_INF), i), asm_opt.chemical_flank, 0.02, 1);
+        if(cov <= msk_cut) msk[i] = cov;
+        if(cov <= msk_cut/**FORCE_CUT**/) {
+            // fprintf(stderr, "-um-[M::%s]\tqn::%u::%.*s\n\n", __func__, (uint32_t)(i), (int)Get_NAME_LENGTH(R_INF, i), Get_NAME((R_INF), i));   
+            for (k = 0; k < paf->length; k++) paf->buffer[k].del = 1; b->cnt[0]++;
+        }
+    } else if(b->cnt[1] == 1) {
+        for (k = 0; k < paf->length; k++) {
+            if((Get_qn(paf->buffer[k])) > (Get_tn(paf->buffer[k]))) continue;
+            rev = &(R_INF.paf[paf->buffer[k].tn]);
+            for (z = 0; z < rev->length; z++) {
+                if((rev->buffer[z].tn == (Get_qn(paf->buffer[k])))) {
+                    if((paf->buffer[k].del != rev->buffer[z].del) || (msk[Get_qn(paf->buffer[k])] <= msk_cut/**FORCE_CUT**/) || (msk[Get_tn(paf->buffer[k])] <= msk_cut/**FORCE_CUT**/)) {
+                        paf->buffer[k].del = rev->buffer[z].del = 1;
+                    }
+                }
+            }
+        }
     }
 
    refresh_ec_ovec_buf_t0(b, REFRESH_N);
@@ -6066,7 +6171,7 @@ void handle_chemical_r(uint64_t n_thre, uint64_t n_a)
 
     kt_for(n_thre, worker_hap_dc_ec_chemical_r, b, n_a);
 
-    fprintf(stderr, "[M::%s] # chemical reads: %lu, # arcs:: %lu\n", __func__, chem_n, dedup);
+    fprintf(stderr, "[M::%s] # chimeric reads: %lu, # arcs:: %lu\n", __func__, chem_n, dedup);
 
     destroy_ec_ovec_buf_t(b);
 }
@@ -6076,16 +6181,44 @@ void handle_chemical_arc(uint64_t n_thre, uint64_t n_a)
     ec_ovec_buf_t *b = NULL; uint64_t k, chem_n = 0;
     b = gen_ec_ovec_buf_t(n_thre);
     for (k = 0; k < n_thre; ++k) {
-        b->a[k].cnt[0] = 0; 
+        b->a[k].cnt[0] = 0; b->a[k].cnt[1] = 0;
     }
 
     kt_for(n_thre, worker_hap_dc_ec_chemical_arc, b, n_a);
 
     for (k = 0; k < n_thre; ++k) {
         chem_n += b->a[k].cnt[0]; 
+        b->a[k].cnt[0] = 0; b->a[k].cnt[1] = 1;
     }
 
-    fprintf(stderr, "[M::%s] # chemical reads: %lu\n", __func__, chem_n);
+    kt_for(n_thre, worker_hap_dc_ec_chemical_arc, b, n_a);
+
+    fprintf(stderr, "[M::%s] # chimeric reads: %lu\n", __func__, chem_n);
 
     destroy_ec_ovec_buf_t(b);
+}
+
+uint8_t* gen_chemical_arc_rf(uint64_t n_thre, uint64_t n_a)
+{
+    ec_ovec_buf_t *b = NULL; uint64_t k, chem_n = 0; uint8_t *ra = NULL;
+    b = gen_ec_ovec_buf_t(n_thre);
+    for (k = 0; k < n_thre; ++k) {
+        b->a[k].cnt[0] = 0; b->a[k].cnt[1] = 0;
+    }
+    MALLOC(ra, n_a); ///memset(ra, -1, sizeof((*ra))*n_a); 
+    b->cr = ra;
+
+    kt_for(n_thre, worker_hap_dc_ec_chemical_arc_mark, b, n_a);
+
+    for (k = 0; k < n_thre; ++k) {
+        chem_n += b->a[k].cnt[0]; 
+        b->a[k].cnt[0] = 0; b->a[k].cnt[1] = 1;
+    }
+
+    kt_for(n_thre, worker_hap_dc_ec_chemical_arc_mark, b, n_a);
+
+    fprintf(stderr, "[M::%s] # chimeric reads: %lu\n", __func__, chem_n);
+
+    b->cr = NULL; destroy_ec_ovec_buf_t(b);
+    return ra;
 }

ecovlp.h

@@ -14,5 +14,6 @@ void sl_ec_r(uint64_t n_thre, uint64_t n_a);
 void cal_ov_r(uint64_t n_thre, uint64_t n_a, uint64_t new_idx);
 void handle_chemical_r(uint64_t n_thre, uint64_t n_a);
 void handle_chemical_arc(uint64_t n_thre, uint64_t n_a);
+uint8_t* gen_chemical_arc_rf(uint64_t n_thre, uint64_t n_a);
 
 #endif

gfa_ut.cpp

@@ -1574,6 +1574,107 @@ uint32_t is_topo, uint32_t min_diff, uint32_t min_ou, uint32_t test_bub, ma_hit_
 }
 
 
+void asg_arc_cut_chimeric_cmk(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_topo, uint32_t min_ou, uint32_t test_bub, uint8_t *cmk, uint32_t cmk_cut)
+{
+    asg64_v tx = {0,0,0}, *b = NULL; asg_arc_t *av, *aw, *ve, *we; 
+    uint32_t i, k, v, w, n_vtx = g->n_seq<<1, nv, nw, kv, kw, ol_max, ou_max, to_del, cnt = 0, mm_ol, mm_ou;
+    
+    if(in) b = in;
+    else b = &tx;
+	b->n = 0;
+
+    for (v = 0; v < n_vtx; ++v) {
+        // if((v>>1)==17078) fprintf(stderr, "[M::%s::] v:%u, del:%u, seq_vis:%u\n", __func__, v, g->seq[v>>1].del, g->seq_vis[v]);
+        if (g->seq[v>>1].del) continue;
+        if (cmk[v>>1] > cmk_cut) continue;
+        if((test_bub == 0) || (g->seq_vis[v] == 0)) {
+            av = asg_arc_a(g, v); nv = asg_arc_n(g, v);
+            if (nv < 2) continue;
+
+            for (i = kv = 0; i < nv; ++i) {
+                if(av[i].del) continue;
+                kv++;
+            }
+            if(kv < 2) continue;
+
+            for (i = 0; i < nv; ++i) {
+                if(av[i].del) continue;
+                kv_push(uint64_t, *b, (((uint64_t)av[i].ol)<<32) | ((uint64_t)(av-g->arc+i)));
+            }
+        }
+	}
+
+    radix_sort_srt64(b->a, b->a + b->n);
+    for (k = 0; k < b->n; k++) {
+        if(g->arc[(uint32_t)b->a[k]].del) continue;
+
+        v = g->arc[(uint32_t)b->a[k]].ul>>32; w = g->arc[(uint32_t)b->a[k]].v^1;
+        if(g->seq[v>>1].del || g->seq[w>>1].del) continue;
+        nv = asg_arc_n(g, v); nw = asg_arc_n(g, w);
+        av = asg_arc_a(g, v); aw = asg_arc_a(g, w);
+        if(nv<=1 && nw <= 1) continue;
+
+        ve = &(g->arc[(uint32_t)b->a[k]]);
+        for (i = 0; i < nw; ++i) {
+            if (aw[i].v == (v^1)) {
+                we = &(aw[i]);
+                break;
+            }
+        }
+        ///mm_ol and mm_ou are used to make edge with long indel more easy to be cutted
+        mm_ol = MIN(ve->ol, we->ol); mm_ou = MIN(ve->ou, we->ou);
+
+        for (i = kv = ol_max = ou_max = 0; i < nv; ++i) {
+            if(av[i].del) continue;
+            kv++;
+            
+            if(ol_max < av[i].ol) ol_max = av[i].ol;
+            if(ou_max < av[i].ou) ou_max = av[i].ou;
+        }
+        if (kv < 1) continue;
+        if (kv >= 2) {
+            if (mm_ol > ol_max*len_rat) continue;
+            if (is_ou && mm_ou > ou_max*ou_rat && mm_ou > min_ou) continue;
+        }
+        
+
+        for (i = kw = ol_max = ou_max = 0; i < nw; ++i) {
+            if(aw[i].del) continue;
+            kw++;
+
+            if(ol_max < aw[i].ol) ol_max = aw[i].ol;
+            if(ou_max < aw[i].ou) ou_max = aw[i].ou;
+        }
+        if (kw < 1) continue;
+        if (kw >= 2) {
+            if (mm_ol > ol_max*len_rat) continue;
+            if (is_ou && mm_ou > ou_max*ou_rat && mm_ou > min_ou) continue;
+        }
+
+        if (kv <= 1 && kw <= 1) continue;
+
+        to_del = 0;
+        if(is_topo) {
+            if (kv > 1 && kw > 1) {
+                to_del = 1;
+            } else if (kw == 1) {
+                if (asg_topocut_aux(g, w^1, max_ext) < max_ext) to_del = 1;
+            } else if (kv == 1) {
+                if (asg_topocut_aux(g, v^1, max_ext + 1) < max_ext + 1) to_del = 1;
+            }
+        }
+
+        if (to_del) {
+            asg_seq_del(g, v>>1); ++cnt;
+        }
+    }
+    // stats_sysm(g);
+    if(!in) free(tx.a);
+    if (cnt > 0) asg_cleanup(g);
+
+    // fprintf(stderr, "[M::%s::] cnt:%u\n", __func__, cnt);
+}
+
 
 void asg_arc_cut_length_adv(asg_t *g, asg64_v *in, int32_t max_ext, float len_rat, float ou_rat, uint32_t is_ou, uint32_t is_trio, 
 uint32_t is_topo, uint32_t min_diff, ma_hit_t_alloc *rev, R_to_U* rI, uint32_t *max_drop_len)
@@ -2909,22 +3010,22 @@ void print_raw_u2rgfa_seq(all_ul_t *aln, R_to_U* rI, uint32_t is_detail)
 }
 
 
-void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, bub_label_t *b_mask_t, long long no_trio_recover)
+void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, bub_label_t *b_mask_t, long long no_trio_recover, uint8_t *cmk)
 {
     rescue_contained_reads_aggressive(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 10, 1, 0, NULL, NULL, b_mask_t);
     rescue_missing_overlaps_aggressive(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 1, 0, NULL, b_mask_t);
     rescue_missing_overlaps_backward(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 10, 1, 0, b_mask_t);
+    if(cmk) rescue_chimeric_reads_aggressive(NULL, sg, src, uopt->coverage_cut, rI, uopt->max_hang, uopt->min_ovlp, 10, 1, 0, NULL, NULL, b_mask_t, cmk);
     // rescue_wrong_overlaps_to_unitigs(NULL, sg, sources, reverse_sources, coverage_cut, ruIndex, 
     // max_hang_length, mini_overlap_length, bubble_dist, NULL);
     // rescue_no_coverage_aggressive(sg, sources, reverse_sources, &coverage_cut, ruIndex, max_hang_length, 
     // mini_overlap_length, bubble_dist, 10);
     set_hom_global_coverage(&asm_opt, sg, uopt->coverage_cut, src, rev, rI, uopt->max_hang, uopt->min_ovlp);
-    rescue_bubble_by_chain(sg, uopt->coverage_cut, src, rev, (asm_opt.max_short_tip*2), 0.15, 3, rI, 0.05, 0.9, uopt->max_hang, 
-    uopt->min_ovlp, 10, uopt->gap_fuzz, b_mask_t, no_trio_recover);
+    rescue_bubble_by_chain(sg, uopt->coverage_cut, src, rev, (asm_opt.max_short_tip*2), 0.15, 3, rI, 0.05, 0.9, uopt->max_hang, uopt->min_ovlp, 10, uopt->gap_fuzz, b_mask_t, no_trio_recover, cmk);
 }
 
 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)
+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, uint8_t *cmk, char *o_file)
 {
     #define HARD_OU_DROP 0.75
     #define HARD_OL_DROP 0.6
@@ -3010,6 +3111,14 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
         //     print_debug_gfa(sg, NULL, uopt->coverage_cut, "UL.dirty3.debug", uopt->sources, uopt->ruIndex, uopt->max_hang, uopt->min_ovlp, 1, 0, 0);
         //     // exit(1);
         // }
+        /**
+        if(cmk && asm_opt.chemical_cov > FORCE_CUT) {
+            asg_arc_identify_simple_bubbles_multi(sg, b_mask_t, 0);
+            asg_arc_cut_chimeric_cmk(sg, &bu, max_tip, 1.1, 1.1, is_ou, 1, 1, 1, cmk, asm_opt.chemical_cov);
+            asg_arc_cut_tips(sg, max_tip, &bu, is_ou, is_ou?rI:NULL, uopt->te);
+        }
+        **/
+
 
         if(is_ou) {
             if(ul_refine_alignment(uopt, sg)) update_sg_uo(sg, src);
@@ -3074,7 +3183,7 @@ double ou_drop_rate, int64_t max_tip, int64_t gap_fuzz, bub_label_t *b_mask_t, i
     set_hom_global_coverage(&asm_opt, sg, uopt->coverage_cut, src, rev, rI, uopt->max_hang, uopt->min_ovlp);
     rescue_bubble_by_chain(sg, uopt->coverage_cut, src, rev, (asm_opt.max_short_tip*2), 0.15, 3, rI, 0.05, 0.9, uopt->max_hang, uopt->min_ovlp, 10, uopt->gap_fuzz, b_mask_t);
     **/
-    post_rescue(uopt, sg, src, rev, rI, b_mask_t, is_ou);
+    post_rescue(uopt, sg, src, rev, rI, b_mask_t, is_ou, cmk);
     
     ug_ext_gfa(uopt, sg, ug_ext_len);
 
@@ -17580,7 +17689,7 @@ ul_renew_t *ropt, const char *bin_file, uint64_t free_uld, uint64_t is_bridg, ui
     (*(ropt->src)) = R_INF.paf; (*(ropt->r_src)) = R_INF.reverse_paf;
     (*(ropt->n_read)) = R_INF.total_reads; (*(ropt->readLen)) = R_INF.read_length;
     renew_R_to_U(ng, (*(ropt->src)), (*(ropt->r_src)), (*(ropt->n_read)), (*(ropt->cov)), ropt->ruIndex, ropt->max_hang, ropt->mini_ovlp);
-    post_rescue(uopt, (*(ropt->sg)), (*(ropt->src)), (*(ropt->r_src)), ropt->ruIndex, ropt->b_mask_t, 0);
+    post_rescue(uopt, (*(ropt->sg)), (*(ropt->src)), (*(ropt->r_src)), ropt->ruIndex, ropt->b_mask_t, 0, NULL);
     // print_raw_uls_aln(uidx, asm_opt.output_file_name);
     // exit(0);
 }

gfa_ut.h

@@ -16,7 +16,7 @@ typedef struct {
 } sset_aux;
 
 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);
+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, uint8_t *cmk, char *o_file);
 uint32_t asg_arc_cut_tips(asg_t *g, uint32_t max_ext, asg64_v *in, uint32_t is_ou, R_to_U *ru, telo_end_t *te);
 void asg_iterative_semi_circ(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t normal_len, uint32_t pop_chimer, asg64_v *dbg, telo_end_t *te);
 void asg_arc_cut_chimeric(asg_t *g, ma_hit_t_alloc* src, asg64_v *in, uint32_t ou_thres, telo_end_t *te);
@@ -33,7 +33,7 @@ void recover_contain_g(asg_t *g, ma_hit_t_alloc *src, R_to_U* ruIndex, int64_t m
 void normalize_gou(asg_t *g);
 void prt_specfic_sge(asg_t *g, uint32_t src, uint32_t dst, const char* cmd);
 asg_t *gen_ng(ma_ug_t *ug, asg_t *sg, ug_opt_t *uopt, ma_sub_t **cov, R_to_U *ruI, uint64_t scaffold_len);
-void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, bub_label_t *b_mask_t, long long no_trio_recover);
+void post_rescue(ug_opt_t *uopt, asg_t *sg, ma_hit_t_alloc *src, ma_hit_t_alloc *rev, R_to_U* rI, bub_label_t *b_mask_t, long long no_trio_recover, uint8_t *cmk);
 // void print_raw_u2rgfa_seq(all_ul_t *aln, R_to_U* rI, uint32_t is_detail);
 bubble_type *gen_bubble_chain(asg_t *sg, ma_ug_t *ug, ug_opt_t *uopt, uint8_t **ir_het, uint8_t avoid_het);
 void filter_sg_by_ug(asg_t *rg, ma_ug_t *ug, ug_opt_t *uopt);