ortholog-based cutting

Overlaps.cpp

@@ -2387,17 +2387,6 @@ uint32_t detect_single_path(asg_t *g, uint32_t begNode, uint32_t* endNode, long
     asg_arc_t *av;
     (*Len) = 0;
 
-    // if(begNode == 18001026)
-    // {
-    //     fprintf(stderr, "(*Len): %d, v: %d, begNode: %d\n", (*Len), v, begNode);
-    //     fflush(stderr);
-    // }
-
-    // if(begNode == 18011441)
-    // {
-    //     fprintf(stderr, "(*Len): %d, v: %d, begNode: %d\n", (*Len), v, begNode);
-    //     fflush(stderr);
-    // }
 
     while (1)
     {
@@ -2408,14 +2397,6 @@ uint32_t detect_single_path(asg_t *g, uint32_t begNode, uint32_t* endNode, long
 
         if(b) kv_push(uint32_t, b->b, v>>1);
 
-
-        // if(begNode == 18001026)
-        // {
-        //     fprintf(stderr, "!!!!!!!!!!!!(*Len): %d, v: %d, begNode: %d, nv: %d, nvr: %d\n", 
-        //     (*Len), v, begNode, nv, asg_arc_n(g, v^1));
-        //     fflush(stderr);
-        // }
-
         if(nv == 0)
         {
             return END_TIPS;
@@ -2742,25 +2723,20 @@ buf_t* b, uint32_t max_ext)
 
     if(b) kv_push(uint32_t, b->b, begNode>>1);
 
-    // if(begNode == 18021291)
-    // {
-    //     fprintf(stderr, "begNode: %d, vn: %d\n", begNode, asg_arc_n(g, begNode));
-    //     fprintf(stderr, "asg_arc_a(g, begNode)[0].v: %d\n", asg_arc_a(g, begNode)[0].v);
-    //     fprintf(stderr, "asg_arc_a(g, begNode)[1].v: %d\n", asg_arc_a(g, begNode)[1].v);
-    //     fflush(stderr);
-    // }
 
-
+    if(detect_single_path_with_single_bubbles(g, asg_arc_a(g, begNode)[0].v, &e1, &l1, b, max_ext) == TWO_INPUT)
-    if(detect_single_path_with_single_bubbles(g, asg_arc_a(g, begNode)[0].v, &e1, &l1, b, max_ext) == TWO_INPUT 
-       && 
-       detect_single_path_with_single_bubbles(g, asg_arc_a(g, begNode)[1].v, &e2, &l2, b, max_ext) == TWO_INPUT)
     {
-        if(e1 == e2)
+        b->b.n--;
+        if(detect_single_path_with_single_bubbles(g, asg_arc_a(g, begNode)[1].v, &e2, &l2, b, max_ext) == TWO_INPUT)
         {
-            (*endNode) = e1;
+            b->b.n--;
-            (*minLen) = (l1 <= l2)? l1: l2;
+            if(e1 == e2)
-            (*minLen)++;
+            {
-            return 1;
+                (*endNode) = e1;
+                (*minLen) = (l1 <= l2)? l1: l2;
+                (*minLen)++;
+                return 1;
+            }
         }
     }
 
@@ -2980,6 +2956,7 @@ uint32_t startNode, uint32_t endNode, int max_dist, buf_t* bub)
 
 
 
+
 int find_single_link(asg_t *g, uint32_t link_beg, int linkLen, uint32_t* link_end)
 {
     uint32_t v, w;
@@ -3609,45 +3586,6 @@ int asg_arc_del_triangular_advance(asg_t *g, long long max_dist)
 }
 
 
-int asg_arc_del_triangular_directly(asg_t *g, long long max_dist)
-{
-    double startTime = Get_T();
-    ///the reason is that each read has two direction (query->target, target->query)
-	uint32_t v, w, n_vtx = g->n_seq * 2, n_reduced = 0;
-
-
-    int flag0, flag1, node;
-    for (v = 0; v < n_vtx; ++v) 
-    {
-        uint32_t nv = asg_arc_n(g, v);
-        asg_arc_t *av = asg_arc_a(g, v);
-        if (g->seq[v>>1].del)
-        {
-            continue;
-        } 
-
-        if(nv < 2)
-        {
-            continue;
-        }
-
-
-        //test_triangular_directly();
-    }
-
-    
-
-    if (n_reduced) {
-		asg_cleanup(g);
-		asg_symm(g);
-	}
-
-    fprintf(stderr, "[M::%s] removed %d triangular overlaps\n", 
-    __func__, n_reduced);
-    fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime);
-
-    return n_reduced;
-}
 
 
 int asg_arc_del_triangular_advance_debug(asg_t *g, long long max_dist)
@@ -5324,11 +5262,98 @@ int asg_arc_del_short_diploid_unclean(asg_t *g, float drop_ratio, ma_hit_t_alloc
 	return n_short;
 }
 
-int asg_arc_del_too_short_overlaps(asg_t *g, long long dropLen)
+
+long long check_if_diploid(uint32_t v1, uint32_t v2, asg_t *g, 
+ma_hit_t_alloc* reverse_sources, long long min_edge_length)
+{
+    buf_t b_0, b_1;
+    memset(&b_0, 0, sizeof(buf_t));
+    memset(&b_1, 0, sizeof(buf_t));
+
+    uint32_t convex1, convex2, f1, f2;
+    long long l1, l2;
+
+    b_0.b.n = 0;
+    b_1.b.n = 0;
+    uint32_t flag1 = detect_single_path(g, v1, &convex1, &l1, &b_0);
+    uint32_t flag2 = detect_single_path(g, v2, &convex2, &l2, &b_1);
+    
+    if(flag1 == LOOP || flag2 == LOOP)
+    {
+        return -1;
+    }
+
+    if(flag1 != END_TIPS && flag1 != LONG_TIPS)
+    {
+        l1--;
+        b_0.b.n--;
+    }
+
+    if(flag2 != END_TIPS && flag2 != LONG_TIPS)
+    {
+        l2--;
+        b_1.b.n--;
+    }
+
+
+    if(l1 <= min_edge_length || l2 <= min_edge_length)
+    {
+        return 0;
+    }
+
+    buf_t* b_min;
+    buf_t* b_max;
+    if(l1<=l2)
+    {
+        b_min = &b_0;
+        b_max = &b_1;
+    }
+    else
+    {
+        b_min = &b_1;
+        b_max = &b_0;
+    }
+
+    long long i, j, k;
+    double max_count = 0;
+    double min_count = 0;
+    uint32_t qn, tn;
+    for (i = 0; i < b_min->b.n; i++)
+    {
+        qn = b_min->b.a[i];
+        for (j = 0; j < reverse_sources[qn].length; j++)
+        {
+            tn = Get_tn(reverse_sources[qn].buffer[j]);
+            if(g->seq[tn].del == 1) continue;
+            min_count++;
+            for (k = 0; k < b_max->b.n; k++)
+            {
+                if(b_max->b.a[k]==tn)
+                {
+                    max_count++;
+                    break;
+                }
+            }
+        }
+    }
+
+
+    free(b_0.b.a);
+    free(b_1.b.a);
+
+    if(min_count == 0) return -1;
+    if(max_count == 0) return 0;
+    if(max_count/min_count>0.3) return 1;
+    return 0;
+
+}
+
+int asg_arc_del_too_short_overlaps(asg_t *g, long long dropLen, float drop_ratio, ma_hit_t_alloc* reverse_sources)
 {
     double startTime = Get_T();
 
 	uint32_t v, n_vtx = g->n_seq * 2, n_short = 0;
+    long long drop_ratio_Len;
 	for (v = 0; v < n_vtx; ++v) 
     {
         if (g->seq[v>>1].del) continue;
@@ -5342,10 +5367,23 @@ int asg_arc_del_too_short_overlaps(asg_t *g, long long dropLen)
 		///remove short overlaps
         if(av[0].ol < dropLen) continue;
 
-		for (i = nv - 1; i >= 1 && av[i].ol < dropLen; --i);
+        drop_ratio_Len = av[0].ol * drop_ratio;
+        if(dropLen < drop_ratio_Len)
+        {
+            drop_ratio_Len = dropLen;
+        }
 
-		for (i = i + 1; i < nv; ++i)
+		for (i = nv - 1; i >= 1 && av[i].ol < drop_ratio_Len; --i);
-			av[i].del = 1, ++n_short;
+
+		// for (i = i + 1; i < nv; ++i)
+		// 	av[i].del = 1, ++n_short;
+        for (i = i + 1; i < nv; ++i)
+        {
+            if(check_if_diploid(av[0].v, av[i].v, g, reverse_sources, 2) != 1)
+            {
+                av[i].del = 1, ++n_short;
+            }
+        }
 	}
 
 
@@ -5497,6 +5535,13 @@ static int asg_topocut_aux(asg_t *g, uint32_t v, int max_ext)
 		}
 		v = asg_check_unambi1(g, v);
 	}
+    /**
+    if(v == (uint32_t)-1 && n_ext < max_ext) //return max_ext + 1;
+    {
+        fprintf(stderr, "v: %llu, n_ext: %llu, max_ext: %llu \n", v, n_ext, max_ext);
+    }
+    **/
+   if(v == (uint32_t)-1) return max_ext + 1;
 	return n_ext;
 }
 
@@ -5563,7 +5608,8 @@ int asg_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio, int max_ext)
 			++kw;
 		}
 		if (kw >= 2 && a->ol > ow_max * drop_ratio) continue;
-		if (kv == 1 && kw == 1) continue;
+		///if (kv == 1 && kw == 1) continue;
+        if (kv <= 1 && kw <= 1) continue;
 
 
         ///to see which one is the current edge (from v and w)
@@ -6536,7 +6582,8 @@ int asg_arc_del_short_diploid_by_exact(asg_t *g, int max_ext, ma_hit_t_alloc* so
 		}
 		if (kw >= 2 && a->ol == ow_max) continue;
 
-		if (kv == 1 && kw == 1) continue;
+		///if (kv == 1 && kw == 1) continue;
+        if (kv <= 1 && kw <= 1) continue;
 
 
         ///to see which one is the current edge (from v and w)
@@ -7735,15 +7782,17 @@ static uint64_t asg_bub_pop1(asg_t *g, uint32_t v0, int max_dist, buf_t *b)
 			if (av[i].del) continue;
 
 			///push the edge
+            ///high 32-bit of g->idx[v] is the start point of v's edges
+            //so here is the point of this specfic edge
 			kv_push(uint32_t, b->e, (g->idx[v]>>32) + i);
 
 			///find a too far path? directly terminate the whole bubble poping
 			if (d + l > max_dist) break; // too far
 
-
+            ///if this node
 			if (t->s == 0) { // this vertex has never been visited
 				kv_push(uint32_t, b->b, w); // save it for revert
-				///t->p means the in-node of w is v
+				///t->p is the parent node of 
 				///t->s = 1 means w has been visited
 				///d is len(v0->v), l is len(v->w), so t->d is len(v0->w)
 				t->p = v, t->s = 1, t->d = d + l;
@@ -7751,10 +7800,13 @@ static uint64_t asg_bub_pop1(asg_t *g, uint32_t v0, int max_dist, buf_t *b)
 				t->r = count_out(g, w^1);
 				++n_pending;
 			} else { // visited before
-				///c seems the max weight of node
+				///c is the weight (is very likely the number of node in this edge) of the parent node
-				if (c + 1 > t->c || (c + 1 == t->c && d + l > t->d)) t->p = v;
+				///select the longest edge (longest meams most reads/longest edge)
+                if (c + 1 > t->c || (c + 1 == t->c && d + l > t->d)) t->p = v;
 				if (c + 1 > t->c) t->c = c + 1;
 				///update len(v0->w)
+                ///node: t->d is not the length from this node's parent
+                ///it is the shortest edge
 				if (d + l < t->d) t->d = d + l; // update dist
 			}
 			///assert(t->r > 0);
@@ -7772,6 +7824,7 @@ static uint64_t asg_bub_pop1(asg_t *g, uint32_t v0, int max_dist, buf_t *b)
 	} while (b->S.n > 1 || n_pending);
 	asg_bub_backtrack(g, v0, b);
 	n_pop = 1 | (uint64_t)b->T.n<<32;
+    fprintf(stderr, "v>>1: %u\n", v0>>1);
 pop_reset:
 	for (i = 0; i < b->b.n; ++i) { // clear the states of visited vertices
 		binfo_t *t = &b->a[b->b.a[i]];
@@ -7809,6 +7862,212 @@ int asg_pop_bubble(asg_t *g, int max_dist)
 }
 
 
+
+
+
+
+
+
+
+
+int test_triangular_directly(asg_t *g, uint32_t v, 
+long long min_edge_length, ma_hit_t_alloc* reverse_sources)
+{
+    
+    uint32_t w;
+    int todel;
+    long long NodeLen_first[3];
+    long long NodeLen_second[3];
+
+    uint32_t Ns_first[3];
+    uint32_t Ns_second[3];
+    
+    uint32_t nv = asg_arc_n(g, v);
+    asg_arc_t *av = asg_arc_a(g, v);
+    if(av[0].v == av[1].v)
+    {
+        return 0;
+    }
+    /**********************test first node************************/
+    NodeLen_first[0] = NodeLen_first[1] = NodeLen_first[2] = -1;
+    if(asg_is_single_edge(g, av[0].v, v>>1) <= 2 && asg_is_single_edge(g, av[1].v, v>>1) <= 2)
+    {
+        NodeLen_first[asg_is_single_edge(g, av[0].v, v>>1)] = 0;
+        NodeLen_first[asg_is_single_edge(g, av[1].v, v>>1)] = 1;
+    }
+    ///one node has one out-edge, another node has two out-edges
+    if(NodeLen_first[1] == -1 || NodeLen_first[2] == -1)
+    {
+        return 0;
+    }
+    /**********************test first node************************/
+    ///if the potiential edge has already been removed 
+    if(av[NodeLen_first[2]].del == 1)
+    {
+        return 0;
+    }
+
+    /**********************test second node************************/
+    w = av[NodeLen_first[2]].v^1;
+    asg_arc_t *aw = asg_arc_a(g, w);
+    uint32_t nw = asg_arc_n(g, w);
+    if(nw != 2)
+    {
+        fprintf(stderr, "error\n");
+    }
+    NodeLen_second[0] = NodeLen_second[1] = NodeLen_second[2] = -1;
+    if(asg_is_single_edge(g, aw[0].v, w>>1) <= 2 && asg_is_single_edge(g, aw[1].v, w>>1) <= 2)
+    {
+        NodeLen_second[asg_is_single_edge(g, aw[0].v, w>>1)] = 0;
+        NodeLen_second[asg_is_single_edge(g, aw[1].v, w>>1)] = 1;
+    }
+    ///one node has one out-edge, another node has two out-edges
+    if(NodeLen_second[1] == -1 || NodeLen_second[2] == -1)
+    {
+        return 0;
+    }
+
+
+
+    todel = 0;
+    // if(check_if_diploid(av[0].v, av[1].v, g, reverse_sources, min_edge_length) && 
+    // check_if_diploid(aw[0].v, aw[1].v, g, reverse_sources, min_edge_length))
+    if(check_if_diploid(av[0].v, av[1].v, g, reverse_sources, min_edge_length) == 1||
+    check_if_diploid(aw[0].v, aw[1].v, g, reverse_sources, min_edge_length) == 1)
+    {
+        todel = 1;
+    }
+    
+    
+    
+    if(todel)
+    {
+        ///fprintf(stderr, "v: %u\n", v>>1);
+        av[NodeLen_first[2]].del = 1;
+        ///remove the reverse direction
+        asg_arc_del(g, av[NodeLen_first[2]].v^1, av[NodeLen_first[2]].ul>>32^1, 1);
+    }
+
+    return todel;
+}
+
+
+
+int asg_arc_del_triangular_directly(asg_t *g, long long min_edge_length, ma_hit_t_alloc* reverse_sources)
+{
+    double startTime = Get_T();
+    ///the reason is that each read has two direction (query->target, target->query)
+    uint32_t v, w, n_vtx = g->n_seq * 2, n_reduced = 0;
+
+
+    int flag0, flag1, node;
+    for (v = 0; v < n_vtx; ++v) 
+    {
+        uint32_t nv = asg_arc_n(g, v);
+        asg_arc_t *av = asg_arc_a(g, v);
+        if (g->seq[v>>1].del)
+        {
+            continue;
+        } 
+
+        if(nv < 2)
+        {
+            continue;
+        }
+
+
+        n_reduced += test_triangular_directly(g, v, min_edge_length, reverse_sources);
+    }
+
+    
+
+    if (n_reduced) {
+        asg_cleanup(g);
+        asg_symm(g);
+    }
+
+    fprintf(stderr, "[M::%s] removed %d triangular overlaps\n", 
+    __func__, n_reduced);
+    fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime);
+
+    return n_reduced;
+}
+
+
+
+
+int asg_arc_del_non_different_haps(asg_t *g, ma_hit_t_alloc* reverse_sources, float drop_ratio)
+{
+    double startTime = Get_T();
+    ///the reason is that each read has two direction (query->target, target->query)
+    uint32_t v, w, n_vtx = g->n_seq * 2, n_reduced = 0;
+    uint32_t idx[2];
+
+    for (v = 0; v < n_vtx; ++v) 
+    {
+        uint32_t i, n_arc = 0, nv = asg_arc_n(g, v);
+        asg_arc_t *av = asg_arc_a(g, v);
+        ///some node could be deleted
+        if (nv < 2 || g->seq[v>>1].del) continue;
+        ///some edges could be deleted
+        for (i = 0; i < nv; ++i) // asg_bub_pop1() may delete some edges/arcs
+            if (!av[i].del) ++n_arc;
+        if (n_arc != 2) continue;
+
+        for (i = 0, n_arc = 0; i < nv; i++)
+        {
+            if (!av[i].del)
+            {
+                idx[n_arc] = i;
+                n_arc++;
+            }
+        }
+
+        if(check_if_diploid(av[idx[0]].v, av[idx[1]].v, g, reverse_sources, 2) == 0)
+        {
+            float max = av[idx[0]].ol;
+            float min = av[idx[1]].ol;
+
+            if(min < drop_ratio * max)
+            {
+                av[idx[1]].del = 1;
+                asg_arc_del(g, av[idx[1]].v^1, av[idx[1]].ul>>32^1, 1);
+                n_reduced++;
+            }
+        }
+    }
+
+
+    
+    
+
+    if (n_reduced) {
+        asg_cleanup(g);
+        asg_symm(g);
+    }
+
+    fprintf(stderr, "[M::%s] removed %d different hap overlaps\n", 
+    __func__, n_reduced);
+    fprintf(stderr, "[M::%s] takes %0.2f s\n\n", __func__, Get_T()-startTime);
+
+    return n_reduced;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 void output_contig_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, long long n_read, long long bubble_dist)
 {
 
@@ -7921,11 +8180,13 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
 
     // debug_info_of_specfic_node("m64016_190918_162737/72220752/ccs", sg, "del_trans");
 
-    ///goto out;
+    
 
     asg_cut_tip(sg, MAX_SHORT_TIPS);
 
 
+    ///goto out;
+    
     // debug_info_of_specfic_node("m64016_190918_162737/72220752/ccs", sg, "cut_tip");
     
 
@@ -7987,6 +8248,8 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
                 ///asg_arc_del_single_node_bubble(sg, bubble_dist);
                 tri_flag += asg_arc_del_single_node_directly(sg, MAX_SHORT_TIPS, sources);
 
+                
+
                 if(tri_flag == 0)
                 {
                     break;
@@ -7994,9 +8257,12 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
             }
 
 
+
+            
             /****************************may have bugs********************************/
             //asg_arc_identify_simple_bubbles(sg);
             asg_arc_identify_simple_bubbles_multi(sg, 1);
+            //reomve edge between two chromesomes
             asg_arc_del_false_node(sg, MAX_SHORT_TIPS);
             asg_cut_tip(sg, MAX_SHORT_TIPS);
             /****************************may have bugs********************************/
@@ -8009,8 +8275,8 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
             asg_cut_tip(sg, MAX_SHORT_TIPS);
             /****************************may have bugs********************************/
 
-
+            asg_arc_del_non_different_haps(sg, reverse_sources, drop_ratio);
-        
+            asg_cut_tip(sg, MAX_SHORT_TIPS);
 
             //asg_arc_identify_simple_bubbles(sg);
             asg_arc_identify_simple_bubbles_multi(sg, 1);
@@ -8066,9 +8332,16 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
     asg_cut_tip(sg, MAX_SHORT_TIPS);
 
 
-    // asg_arc_identify_simple_bubbles_multi(sg, 0);
+    asg_arc_del_triangular_directly(sg, MAX_SHORT_TIPS, reverse_sources);
-    // asg_arc_del_too_short_overlaps(sg, 1000);
+
-    // asg_cut_tip(sg, MAX_SHORT_TIPS);
+    
+
+    
+    asg_arc_identify_simple_bubbles_multi(sg, 0);
+    asg_arc_del_too_short_overlaps(sg, 2000, min_ovlp_drop_ratio, reverse_sources);
+    asg_cut_tip(sg, MAX_SHORT_TIPS);
+
+    
     
 
 
@@ -8125,7 +8398,8 @@ char* output_file_name, long long bubble_dist, int read_graph, int write)
     output_unitig_graph(sg, coverage_cut, output_file_name, n_read);
     output_read_graph(sg, coverage_cut, output_file_name, n_read);
 
-    output_contig_graph(sg, coverage_cut, output_file_name, n_read, 50000);
+    output_contig_graph(sg, coverage_cut, output_file_name, n_read, 500000);
+    ///output_contig_graph(sg, coverage_cut, output_file_name, n_read, 10000000);
 
     
     asg_destroy(sg);