hic bubble

Overlaps.cpp

@@ -10981,71 +10981,7 @@ void set_pre_uid(buf_t* b, hc_links* link, ma_ug_t *ug)
 }
 
 
-void collect_reverse_unitigs_back(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug)
+inline uint64_t get_utg_len(buf_t* b, ma_ug_t *ug, asg_t *read_sg, uint64_t ignore_end, uint64_t* len_thre, uint64_t* occ)
-{
-    uint32_t k, m;
-    uint64_t d = (uint64_t)-1;
-    set_pre_uid(b_0, link, ug);
-    set_pre_uid(b_1, link, ug);
-    
-    for (k = 0; k < b_0->b.n; k++)
-    {
-        for (m = 0; m < b_1->b.n; m++)
-        {
-            if(b_0->b.a[k] == b_1->b.a[m]) continue;
-            push_hc_edge(&(link->a.a[b_0->b.a[k]]), b_1->b.a[m], 1, 1, &d);
-            push_hc_edge(&(link->a.a[b_1->b.a[m]]), b_0->b.a[k], 1, 1, &d);
-        }
-    }
-}
-
-
-void collect_reverse_unitigs_back(buf_t* b_0, uint32_t b_0_uLen, buf_t* b_1, uint32_t b_1_uLen, hc_links* link, ma_ug_t *ug)
-{
-    uint32_t b_0_i, b_0_i_l, b_0_k, b_1_i, b_1_i_l, b_1_k, pre_0, pre_1, rId_0, rId_1;
-    uint64_t d = (uint64_t)-1;
-    ma_utg_t* u_b_0 = NULL;
-    ma_utg_t* u_b_1 = NULL;
-    if(b_0->b.n == 0 || b_1->b.n == 0) return;
-    if(b_0_uLen == 0 || b_1_uLen == 0) return;
-
-    for (b_0_i = b_0_i_l = 0, pre_0 = (uint32_t)-1; b_0_i < b_0->b.n; b_0_i++)
-    {
-        u_b_0 = &(ug->u.a[b_0->b.a[b_0_i]>>1]);
-        if(u_b_0->n == 0) continue;
-        for (b_0_k = 0; b_0_k < u_b_0->n; b_0_k++)
-        {
-            rId_0 = u_b_0->a[b_0_k]>>33;
-            if(link->u_idx[rId_0] == (uint32_t)-1) continue;
-            if(pre_0 == link->u_idx[rId_0]) continue;
-            pre_0 = link->u_idx[rId_0];
-            b_0_i_l++;
-            if(b_0_i_l > b_0_uLen) return;
-
-            for (b_1_i = b_1_i_l = 0, pre_1 = (uint32_t)-1; b_1_i < b_1->b.n; b_1_i++)
-            {
-                u_b_1 = &(ug->u.a[b_1->b.a[b_1_i]>>1]);
-                if(u_b_1->n == 0) continue;
-                for (b_1_k = 0; b_1_k < u_b_1->n; b_1_k++)
-                {
-                    rId_1 = u_b_1->a[b_1_k]>>33;
-                    if(link->u_idx[rId_1] == (uint32_t)-1) continue;
-                    if(pre_1 == link->u_idx[rId_1]) continue;
-                    pre_1 = link->u_idx[rId_1];
-                    b_1_i_l++;
-                    if(b_1_i_l > b_1_uLen) goto b_1_i_end;
-
-                    push_hc_edge(&(link->a.a[pre_0]), pre_1, 1, 1, &d);
-                    push_hc_edge(&(link->a.a[pre_1]), pre_0, 1, 1, &d);
-                }
-            }
-
-            b_1_i_end:;
-        }
-    }
-}
-
-inline uint64_t get_utg_len(buf_t* b, ma_ug_t *ug, asg_t *read_sg, uint64_t* len_thre, uint64_t* occ)
 {
     if(len_thre && occ)(*occ) = (uint64_t)-1;
     uint32_t ori, uid, v, nv, l, k, idx;
@@ -11085,7 +11021,7 @@ inline uint64_t get_utg_len(buf_t* b, ma_ug_t *ug, asg_t *read_sg, uint64_t* len
                 }
                 if(k == nv) fprintf(stderr, "ERROR\n");
                 len += l;
-                if(len_thre && occ && len >= (*len_thre) && (*occ) != (uint64_t)-1)
+                if(len_thre && occ && len >= (*len_thre))
                 {
                     (*occ) = idx - 1;
                     return len;
@@ -11120,7 +11056,7 @@ inline uint64_t get_utg_len(buf_t* b, ma_ug_t *ug, asg_t *read_sg, uint64_t* len
                 }
                 if(k == nv) fprintf(stderr, "ERROR\n");
                 len += l;
-                if(len_thre && occ && len >= (*len_thre) && (*occ) != (uint64_t)-1)
+                if(len_thre && occ && len >= (*len_thre))
                 {
                     (*occ) = idx - 1;
                     return len;
@@ -11130,33 +11066,38 @@ inline uint64_t get_utg_len(buf_t* b, ma_ug_t *ug, asg_t *read_sg, uint64_t* len
         }
     }
 
-    if(p_v != (uint32_t)-1)
+    if(ignore_end == 0 && p_v != (uint32_t)-1)
     {
         len += read_sg->seq[p_v>>1].len;
-        if(len_thre && occ && len >= (*len_thre) && (*occ) != (uint64_t)-1)
+        if(len_thre && occ && len >= (*len_thre))
         {
             (*occ) = idx - 1;
             return len;
         }
     }
 
+    if(len_thre && occ)
+    {
+        (*occ) = idx;
+    }
+
     return len;
 }
 
 void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug, asg_t *read_sg)
 {
     uint32_t b_0_i, b_0_k, b_1_i, b_1_k, pre_0, pre_1, rId_0, rId_1, ori_0, ori_1;
-    uint64_t d = (uint64_t)-1, len_0, len_1, thre_0, thre_1;
+    uint64_t d = RC_1, len_0, len_1, thre_0, thre_1;
     ma_utg_t* u_b_0 = NULL;
     ma_utg_t* u_b_1 = NULL;
     if(b_0->b.n == 0 || b_1->b.n == 0) return;
 
-    len_0 = get_utg_len(b_0, ug, read_sg, NULL, NULL);
+    len_0 = get_utg_len(b_0, ug, read_sg, 1, NULL, NULL);
-    len_1 = get_utg_len(b_1, ug, read_sg, NULL, NULL);
+    len_1 = get_utg_len(b_1, ug, read_sg, 1, NULL, NULL);
 
     len_0 = MIN(len_0, len_1);
-    get_utg_len(b_0, ug, read_sg, &len_0, &thre_0);
+    get_utg_len(b_0, ug, read_sg, 0, &len_0, &thre_0);
-    get_utg_len(b_1, ug, read_sg, &len_0, &thre_1);
+    get_utg_len(b_1, ug, read_sg, 0, &len_0, &thre_1);
 
     
     for (b_0_i = len_0 = 0, pre_0 = (uint32_t)-1; b_0_i < b_0->b.n; b_0_i++)
@@ -11202,9 +11143,25 @@ void collect_reverse_unitigs(buf_t* b_0, buf_t* b_1, hc_links* link, ma_ug_t *ug
                     if(pre_1 == link->u_idx[rId_1]) continue;
                     pre_1 = link->u_idx[rId_1];
                 
-
+                    
                     push_hc_edge(&(link->a.a[pre_0]), pre_1, 1, 1, &d);
                     push_hc_edge(&(link->a.a[pre_1]), pre_0, 1, 1, &d);
+
+                    // if(pre_0 == 5 || pre_1 == 5)
+                    // {
+                    //     fprintf(stderr, "\npre_0: utg%.6ul, len_0: %lu, thre_0: %lu, pre_1: utg%.6ul, len_1: %lu, thre_1: %lu\n", 
+                    //     (int)(pre_0+1), len_0, thre_0, (int)(pre_1+1), len_1, thre_1);
+                    //     uint32_t xxx_i;
+                    //     for (xxx_i = 0; xxx_i < b_0->b.n; xxx_i++)
+                    //     {
+                    //         fprintf(stderr,"+: utg%.6ul\n", (int)((b_0->b.a[xxx_i]>>1)+1));
+                    //     }
+
+                    //     for (xxx_i = 0; xxx_i < b_1->b.n; xxx_i++)
+                    //     {
+                    //         fprintf(stderr,"-: utg%.6ul\n", (int)((b_1->b.a[xxx_i]>>1)+1));
+                    //     }
+                    // }
                 }
             }
 

Purge_Dups.cpp

@@ -3931,7 +3931,7 @@ kvec_t_i32_warp* prevIndex, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edg
 void collect_reverse_unitig_pair(hc_links* link, ma_ug_t *ug, hap_overlaps* t)
 {
     uint32_t i = 0, k = 0, rId_0, rId_1, pre_0, pre_1, b_0 = t->xUid, b_1 = t->yUid;
-    uint64_t d = (uint64_t)-1;
+    uint64_t d = RC_2;
     ma_utg_t* u_b_0 = &(ug->u.a[b_0]);
     ma_utg_t* u_b_1 = &(ug->u.a[b_1]);
     if(u_b_0->n == 0) return;

hic.cpp

@@ -114,6 +114,7 @@ typedef struct {
     kvec_t(uint32_t) list;
     kvec_t(uint32_t) num;
     kvec_t(uint64_t) pathLen;
+    uint64_t f_bub, b_bub;
 } bubble_type;
 
 typedef struct {
@@ -165,6 +166,11 @@ typedef struct {
     kvec_t(pe_hit) a;
 } kvec_pe_hit;
 
+typedef struct {
+    kvec_t(hc_edge) a;
+}kvec_hc_edge;
+
+
 #define pe_hit_an1_key(x) ((x).s)
 KRADIX_SORT_INIT(pe_hit_an1, pe_hit, pe_hit_an1_key, 8)
 #define pe_hit_an2_key(x) ((x).e)
@@ -229,6 +235,8 @@ typedef struct {
 reads_t R1, R2;
 ha_ug_index* ug_index;
 
+void build_bub_graph(ma_ug_t* ug, bubble_type* bub);
+
 void init_ha_ug_index_opt(ha_ug_index* idx, ma_ug_t *ug, int k, pldat_t* p)
 {
     uint64_t i, n;
@@ -1416,12 +1424,12 @@ inline void get_bubbles(bubble_type* bub, uint64_t id, uint32_t* beg, uint32_t*
     if(pathBase) (*pathBase) = bub->pathLen.a[id];
 }
 
-void identify_bubbles(ma_ug_t* ug, bubble_type* bub)
+void identify_bubbles(ma_ug_t* ug, bubble_type* bub, hc_links* link)
 {
     asg_cleanup(ug->g);
     if (!ug->g->is_symm) asg_symm(ug->g);
     memset(bub, 0, sizeof(bubble_type));
-    uint32_t v, n_vtx = ug->g->n_seq * 2, tLen, i, mode = (((uint32_t)-1)<<2);
+    uint32_t v, n_vtx = ug->g->n_seq * 2, tLen, i, k, mode = (((uint32_t)-1)<<2);
     uint64_t pathLen;
     bub->ug = ug; 
     CALLOC(bub->index, n_vtx);
@@ -1511,15 +1519,24 @@ void identify_bubbles(ma_ug_t* ug, bubble_type* bub)
 
         if(bub->index[(beg>>1)] != (bub->num.n + 1)) bub->index[(beg>>1)] = (uint32_t)-1;
         if(bub->index[(sink>>1)] != (bub->num.n + 1)) bub->index[(sink>>1)] = (uint32_t)-1;
-        ///bub->index[(beg>>1)] = bub->index[(sink>>1)] = (uint32_t)-1;
     }
 
     for (i = 0; i < ug->g->n_seq; i++)
     {
         if(bub->index[i] == bub->num.n + 1) bub->index[i] = bub->num.n;
+        if(bub->index[i] > bub->num.n)
+        {
+            for (k = 0; k < link->a.a[i].f.n; k++)
+            {
+                if(link->a.a[i].f.a[k].del || link->a.a[i].f.a[k].dis != RC_1) continue;
+                bub->index[i] = bub->num.n;
+                ///fprintf(stderr, "utg%.6ul\n", (int)(i+1));
+                break;
+            }
+        }
     }
 
-    ///free(bub->index); bub->index = NULL;
+    build_bub_graph(ug, bub);
 }
 
 
@@ -1922,12 +1939,13 @@ void pop_pdq(pdq* q, uint64_t* min_v, uint64_t* min_dis)
 }
 
 
-void get_shortest_path(uint32_t src, pdq* pq, asg_t *sg)
+void get_shortest_path(uint32_t src, pdq* pq, asg_t *sg, uint32_t* pre)
 {
     uint64_t v, u, i, nv, w;
     asg_arc_t *av = NULL;
     reset_pdq(pq);
     pq->dis.a[src] = 0;
+    if(pre) pre[src] = (uint32_t)-1;
     push_pdq(pq, src, 0);
     while (pdq_cnt(*pq) > 0)
     {
@@ -1947,6 +1965,7 @@ void get_shortest_path(uint32_t src, pdq* pq, asg_t *sg)
             {
                 pq->dis.a[u] = pq->dis.a[v] + w;
                 push_pdq(pq, u, pq->dis.a[u]);
+                if(pre) pre[u] = v;
             }
         }
     }
@@ -1969,7 +1988,7 @@ void all_pair_shortest_path(const ha_ug_index* idx, hc_links* link, MT* M)
         if (sg->seq[v>>1].del) continue;
         t = &(link->a.a[v>>1]);
         if (t->e.n == 0) continue;
-        get_shortest_path(v, &pq, sg);
+        get_shortest_path(v, &pq, sg, NULL);
         for (k = 0; k < pq.dis.n; k++)
         {
             if(pq.dis.a[k] == (uint64_t)-1) continue;
@@ -2477,7 +2496,7 @@ uint32_t v, uint32_t beg, uint32_t sink)
 
 void set_reverse_links(uint32_t* bub, uint32_t n, kvec_t_u32_warp* reach, uint32_t root, hc_links* link)
 {
-    uint64_t i, k, d = 0;
+    uint64_t i, k, d = RC_0;
     uint32_t v;
     for (i = 0; i < n; i++)
     {
@@ -2498,11 +2517,13 @@ void set_reverse_links(uint32_t* bub, uint32_t n, kvec_t_u32_warp* reach, uint32
 
 void collect_hc_reverse_links(hc_links* link, ma_ug_t* ug, bubble_type* bub)
 {
-    uint64_t i, j, k, d = 0, m, pre;
+    uint64_t i, j, k, d = RC_0, m, pre;
     uint32_t beg, sink, n, v, *a = NULL;
     kvec_t_u32_warp stack, result;
     hc_edge *e = NULL;
     kv_init(stack.a); kv_init(result.a);
+    ///clean all reverse overlaps within bubbles
+    ///might be wrong
     for (i = 0; i < bub->num.n-1; i++)
     {
         get_bubbles(bub, i, &beg, &sink, &a, &n, NULL);
@@ -2565,6 +2586,7 @@ void collect_hc_reverse_links(hc_links* link, ma_ug_t* ug, bubble_type* bub)
     }
     kv_destroy(stack.a); kv_destroy(result.a);
 
+
     for (i = 0; i < link->a.n; i++)
     {
         for (k = m = 0; k < link->a.a[i].f.n; k++)
@@ -2581,7 +2603,7 @@ void collect_hc_reverse_links(hc_links* link, ma_ug_t* ug, bubble_type* bub)
             if(link->a.a[i].f.a[k].del) continue;
             if(link->a.a[i].f.a[k].uID == pre)
             {
-                if(link->a.a[i].f.a[k].dis == 0) link->a.a[i].f.a[m-1].dis = 0;
+                if(link->a.a[i].f.a[k].dis == RC_0) link->a.a[i].f.a[m-1].dis = RC_0;
                 continue;
             }
              
@@ -2594,9 +2616,6 @@ void collect_hc_reverse_links(hc_links* link, ma_ug_t* ug, bubble_type* bub)
     }
 
 
-
-
-
     // hc_edge *e = NULL;
     // for (i = 0; i < link->a.n; i++)
     // {
@@ -4198,7 +4217,80 @@ void weight_edges(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_ty
     }
 }
 
-void init_hic_p(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub)
+void build_bub_graph(ma_ug_t* ug, bubble_type* bub)
+{
+    asg_t *sg = ug->g;
+    pdq pq;
+    init_pdq(&pq, sg->n_seq<<1);
+    uint32_t n_vtx = sg->n_seq<<1, v;
+    uint32_t *pre = NULL; MALLOC(pre, n_vtx);
+    uint64_t *flag = NULL; CALLOC(flag, sg->n_seq); 
+    uint64_t k, b_id_0, b_id_1;
+    uint32_t beg, sink, n, *a, pre_id, adjecent;
+    for (k = 0; k < bub->num.n-1; k++)
+    {
+        get_bubbles(bub, k, &beg, &sink, &a, &n, NULL);
+        pre_id = beg>>1;
+        if(bub->index[pre_id] > bub->num.n)
+        {
+            flag[pre_id] |= 1;
+            flag[pre_id] |= (uint64_t)((uint64_t)pre_id<<33);
+        }
+        
+        pre_id = sink>>1;
+        if(bub->index[pre_id] > bub->num.n)
+        {
+            flag[pre_id] |= 2;
+            flag[pre_id] |= (uint64_t)((uint64_t)pre_id<<2);
+        }
+    }
+
+    for (v = 0; v < n_vtx; ++v)
+    {
+        if(sg->seq[v>>1].del) continue;
+        if(flag[v>>1] == 0) continue;
+        if((flag[v>>1] & 1) && (flag[v>>1] & 2))
+        {
+            fprintf(stderr, "+utg%.6ul -> utg%.6ul\n", 
+                        (int)((flag[v>>1]>>33)+1), (int)((flag[v>>1]>>2)+1));
+            continue;
+        }
+
+        get_shortest_path(v, &pq, sg, pre);
+        for (k = 0; k < pq.dis.n; k++)
+        {
+            if(pq.dis.a[k] == (uint64_t)-1) continue;
+            if((flag[k>>1]&3) == 0) continue;
+            if(k == v) continue;
+            pre_id = pre[k];
+            adjecent = 0;
+            while (pre_id != v)
+            {
+                if(flag[pre_id>>1] > 0)
+                {
+                    adjecent = 1;
+                    break;
+                } 
+                pre_id = pre[pre_id];
+            }
+
+            if(adjecent == 0)
+            {
+                b_id_0 = flag[k>>1] & 1? flag[k>>1]>>33 : flag[k>>1]>>2;
+                b_id_1 = flag[v>>1] & 1? flag[v>>1]>>33 : flag[v>>1]>>2;
+                if((flag[k>>1] & 3) == 3) fprintf(stderr, "ERROR1\n");
+                if((flag[v>>1] & 3) == 3) fprintf(stderr, "ERROR2\n");
+                fprintf(stderr, "-utg%.6ul -> utg%.6ul\n", (int)(b_id_0+1), (int)(b_id_1+1));
+            }
+        }
+    }
+
+    free(pre);
+    free(flag);
+    destory_pdq(&pq);
+}
+
+void init_hic_p(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type* bub, kvec_hc_edge* back_hc_edge)
 {
     uint64_t k, i, m, beg, end, t_d, r_idx, f_idx, b_size, med = (uint64_t)-1, uID;
     uint32_t b_beg, b_end, n, *a, b_cnt;
@@ -4265,8 +4357,6 @@ void init_hic_p(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type
     {
         if((buf.a[k]&1) == 0) r_idx = k;
         if((buf.a[k]&1) == 1) f_idx = k;
-
-        ///fprintf(stderr, "%lu\t%lu\n", buf.a[k]>>1, buf.a[k]&1);
     }
     buf.n = MIN(r_idx, f_idx);
 
@@ -4372,19 +4462,29 @@ void init_hic_p(ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link, bubble_type
     fprintf(stderr, "idx->a: %f, idx->b: %f, idx->frac: %f, med: %lu\n", 
     (double)idx->a, (double)idx->b, (double)idx->frac, med);
     
-
+    back_hc_edge->a.n = 0;
     hc_edge *e = NULL;
     for (i = 0; i < link->a.n; i++)
     {
         for (k = 0; k < link->a.a[i].f.n; k++)
         {
             if(link->a.a[i].f.a[k].del) continue;
-            if(link->a.a[i].f.a[k].dis != 0) continue;
+            if(link->a.a[i].f.a[k].dis != RC_0) continue;
             uID = link->a.a[i].f.a[k].uID;
+
             e = get_hc_edge(link, i, uID, 0);
-            if(e) e->del = 1;
+            if(e) 
+            {
+                kv_push(hc_edge, back_hc_edge->a, *e);
+                e->del = 1;
+            }
+
             e = get_hc_edge(link, uID, i, 0);
-            if(e) e->del = 1;
+            if(e) 
+            {
+                kv_push(hc_edge, back_hc_edge->a, *e);
+                e->del = 1;
+            }
         }
     }
 
@@ -5355,6 +5455,8 @@ int hic_short_align(const char *fn1, const char *fn2, ha_ug_index* idx)
     double index_time = yak_realtime();
     sldat_t sl;
     gzFile fp1, fp2;
+    kvec_hc_edge back_hc_edge;
+    kv_init(back_hc_edge.a);
     if ((fp1 = gzopen(fn1, "r")) == 0) return 0;
     if ((fp2 = gzopen(fn2, "r")) == 0) return 0;
     sl.ks1 = kseq_init(fp1);
@@ -5368,7 +5470,7 @@ int hic_short_align(const char *fn1, const char *fn2, ha_ug_index* idx)
     fprintf(stderr, "u.n: %d, uID_bits: %lu, pos_bits: %lu\n", (uint32_t)idx->ug->u.n, idx->uID_bits, idx->pos_bits);
     
     bubble_type bub;
-    identify_bubbles(idx->ug, &bub);
+    identify_bubbles(idx->ug, &bub, idx->link);
     ///print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, idx->link, idx);
 
     if(!load_hc_hits(&sl.hits, asm_opt.output_file_name))
@@ -5394,7 +5496,7 @@ int hic_short_align(const char *fn1, const char *fn2, ha_ug_index* idx)
     
     ///debug_hc_links(idx, idx->link, &sl, &bub, fn1);
     
-    init_hic_p((ha_ug_index*)sl.idx, &sl.hits, idx->link, &bub);
+    init_hic_p((ha_ug_index*)sl.idx, &sl.hits, idx->link, &bub, &back_hc_edge);
     
     H_partition hap;
     init_contig_partition(&hap, idx, &bub);
@@ -5404,8 +5506,17 @@ int hic_short_align(const char *fn1, const char *fn2, ha_ug_index* idx)
     ///print_hc_links(idx->link, 0, &hap);
     ///print_contig_partition(&hap, "final");
 
-    destory_contig_partition(&hap);
+    // uint32_t i;
+    // for (i = 0; i < idx->ug->g->n_seq; i++)
+    // {
+    //     fprintf(stderr, "utg%.6ul, index: %u\n", (int)(i+1), bub.index[i]);
+    // }
+    
 
+
+
+    destory_contig_partition(&hap);
+    kv_destroy(back_hc_edge.a);
     return 1;
     
     /*******************************for debug************************************/

hic.h

@@ -5,6 +5,9 @@
 
 #define kdq_clear(q) ((q)->count = (q)->front = 0)
 #define kv_malloc(v, s) ((v).n = 0, (v).m = (s), MALLOC((v).a, (s)))
+#define RC_0 0
+#define RC_1 1
+#define RC_2 2
 
 hc_edge* get_hc_edge(hc_links* link, uint64_t src, uint64_t dest, uint64_t dir);
 void push_hc_edge(hc_linkeage* x, uint64_t uID, int weight, int dir, uint64_t* d);