Merge branch 'master' into dev-lh3

Assembly.cpp

@@ -1190,7 +1190,6 @@ int ha_assemble(void)
 	build_string_graph_without_clean(asm_opt.min_overlap_coverage, R_INF.paf, R_INF.reverse_paf, 
 			R_INF.total_reads, R_INF.read_length, asm_opt.min_overlap_Len, asm_opt.max_hang_Len, asm_opt.clean_round, 
 			asm_opt.gap_fuzz, asm_opt.min_drop_rate, asm_opt.max_drop_rate, asm_opt.output_file_name, asm_opt.large_pop_bubble_size, 0, !ovlp_loaded);
-	if (asm_opt.verbose_gfa) R_INF.paf = R_INF.reverse_paf = 0; // with --dbg-gfa, R_INF.paf and R_INF.reverse_paf are freed in load_debug_graph()
 	destory_All_reads(&R_INF);
 	return 0;
 }

Overlaps.cpp

@@ -22194,6 +22194,8 @@ char* output_file_name, ma_hit_t_alloc** reverse_sources, R_to_U* ruIndex)
         return 0;
     }
 
+    R_INF.paf = (*sources); R_INF.reverse_paf = (*reverse_sources);
+
     return 1;
 }
 
@@ -22276,7 +22278,8 @@ kvec_asg_arc_t_warp* new_rtg_edges)
 
 
 
-
+    
+    enable_debug_mode(1);
 
     purge_dups(*ug, read_g, coverage_cut, reverse_sources, ruIndex, new_rtg_edges, 0.75, 50, 50, 0.5, max_hang,
     min_ovlp, bubble_dist, drop_ratio, 0);
@@ -22286,28 +22289,25 @@ kvec_asg_arc_t_warp* new_rtg_edges)
 
 
 
-
     rescue_missing_overlaps_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang,
     min_ovlp, 0, 0, 1, NULL);
-
     renew_utg(ug, read_g, new_rtg_edges);
 
     rescue_contained_reads_aggressive(*ug, read_g, sources, coverage_cut, ruIndex, max_hang, 
     min_ovlp, 0, 10, 0, 1, NULL, NULL);
-
     renew_utg(ug, read_g, new_rtg_edges);
 
+
     ///debug_purge_dup = 1;
     ///deduplicate_advance(*ug, read_g, coverage_cut, sources, reverse_sources, 20, 100, 0.05, 0.2, ruIndex, 0);
     ///delete_useless_nodes(ug);
-    enable_debug_mode();
+    enable_debug_mode(0);
     purge_dups(*ug, read_g, coverage_cut, reverse_sources, ruIndex, new_rtg_edges, 0.75, 50, 50, 0.5, max_hang,
     min_ovlp, bubble_dist, drop_ratio, 0);
     delete_useless_nodes(ug);
     renew_utg(ug, read_g, new_rtg_edges);
     
 
-    ///set_drop_trio_flag(*ug);
 
     n_vtx = read_g->n_seq;
     for (v = 0; v < n_vtx; v++)

Purge_Dups.cpp

@@ -5,6 +5,7 @@
 #include "Purge_Dups.h"
 #include "Overlaps.h"
 #include "Correct.h"
+#include "kthread.h"
 
 #define Cal_Off(OFF) ((long long)((uint32_t)((OFF)>>32)) - (long long)((uint32_t)((OFF))))
 #define Get_match(x) ((x).weight)
@@ -98,7 +99,7 @@ void init_hap_alignment_struct(hap_alignment_struct* x, uint32_t size)
     kv_init(x->u_can.a);
 }
 
-void destory_hap_alignment_struct(hap_alignment_struct* x, uint32_t size)
+void destory_hap_alignment_struct(hap_alignment_struct* x)
 {
     free(x->vote_counting);
     free(x->visit);
@@ -107,6 +108,61 @@ void destory_hap_alignment_struct(hap_alignment_struct* x, uint32_t size)
     kv_destroy(x->u_can.a);
 }
 
+
+typedef struct {
+    hap_alignment_struct* buf;
+    uint32_t num_threads;
+
+    ma_ug_t *ug;
+    asg_t *read_g;
+    ma_hit_t_alloc* reverse_sources;
+    R_to_U* ruIndex;
+    ma_sub_t *coverage_cut; 
+    uint64_t* position_index;
+    float Hap_rate;
+    int max_hang;
+    int min_ovlp;
+    float chain_rate;
+    hap_overlaps_list* all_ovlp;
+}hap_alignment_struct_pip;
+
+void init_hap_alignment_struct_pip(hap_alignment_struct_pip* x, uint32_t num_threads, uint32_t n_seq,
+ma_ug_t *ug, asg_t *read_g,  ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, ma_sub_t *coverage_cut, 
+uint64_t* position_index, float Hap_rate, int max_hang, int min_ovlp, float chain_rate, hap_overlaps_list* all_ovlp)
+{
+    uint32_t i;
+    x->num_threads = num_threads;
+    x->buf = (hap_alignment_struct*)malloc(sizeof(hap_alignment_struct)*x->num_threads);
+    for (i = 0; i < x->num_threads; i++)
+    {
+        init_hap_alignment_struct(&(x->buf[i]), n_seq);
+    }
+
+    x->ug = ug;
+    x->read_g = read_g;
+    x->reverse_sources = reverse_sources;
+    x->ruIndex = ruIndex;
+    x->coverage_cut = coverage_cut;
+    x->position_index = position_index;
+    x->Hap_rate = Hap_rate;
+    x->max_hang = max_hang;
+    x->min_ovlp = min_ovlp;
+    x->chain_rate = chain_rate;
+    x->all_ovlp = all_ovlp;
+}
+
+
+void destory_hap_alignment_struct_pip(hap_alignment_struct_pip* x)
+{
+    uint32_t i;
+    for (i = 0; i < x->num_threads; i++)
+    {
+        destory_hap_alignment_struct(&(x->buf[i]));
+    }
+
+    free(x->buf);
+}
+
 void init_hap_overlaps_list(hap_overlaps_list* x, uint32_t num)
 {
     uint32_t i = 0;
@@ -118,9 +174,9 @@ void init_hap_overlaps_list(hap_overlaps_list* x, uint32_t num)
     }
 }
 
-void enable_debug_mode()
+void enable_debug_mode(uint32_t mode)
 {
-    debug_enable = 1;
+    debug_enable = mode;
 }
 
 void destory_hap_overlaps_list(hap_overlaps_list* x)
@@ -411,6 +467,7 @@ ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match,
         {
             max_count = 0;
             min_count = Len;
+            break;
         }
         qn = readIDs[i]>>33;
         if(reverse_sources[qn].length > 0) min_count++;
@@ -446,6 +503,60 @@ ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match,
     (*Match) = max_count;
     (*Total) = min_count;
 }
+/**
+void get_pair_hap_similarity_deduplicate(uint64_t* readIDs, uint32_t Len, uint32_t target_uId, 
+ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, double* Match, double* Total)
+{
+    get_pair_hap_similarity(readIDs, Len, target_uId, reverse_sources, read_g, ruIndex, Match, Total);
+    return;
+
+    #define CUTOFF_THRES 100
+    uint32_t i, j, qn, tn, is_Unitig, uId, min_count = 0, max_count = 0, cutoff = 0, is_found;
+    for (i = 0; i < Len; i++)
+    {
+        if(cutoff > CUTOFF_THRES)
+        {
+            max_count = 0;
+            min_count = Len;
+            break;
+        }
+        qn = readIDs[i]>>33;
+        is_found = 0;
+
+        for (j = 0; j < reverse_sources[qn].length; j++)
+        {
+            tn = Get_tn(reverse_sources[qn].buffer[j]);
+            if(read_g->seq[tn].del == 1)
+            {
+                get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
+                if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue;
+            }
+
+
+            get_R_to_U(ruIndex, tn, &uId, &is_Unitig);
+            if(uId!=(uint32_t)-1 && is_Unitig == 1 && uId == target_uId)
+            {
+                max_count++;
+            }
+            min_count++;
+            is_found = 1;
+        }
+
+        //means there is a match
+        if(is_found)
+        {
+            cutoff = 0;
+        }
+        else
+        {
+            cutoff++;
+        }
+    }
+
+    (*Match) = max_count;
+    (*Total) = min_count;
+}
+**/
 
 inline void check_hap_match(uint32_t qn, uint32_t targetBeg, uint32_t targetEnd, uint32_t targetID, 
 uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, uint32_t* is_found, uint32_t* is_match)
@@ -477,6 +588,42 @@ uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U
 
 }
 
+/**
+inline void check_hap_match_deduplicate(uint32_t qn, uint32_t targetBeg, uint32_t targetEnd, uint32_t targetID, 
+uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g, R_to_U* ruIndex, uint32_t* is_found, uint32_t* is_match)
+{
+    check_hap_match(qn, targetBeg, targetEnd, targetID, position_index, reverse_sources, read_g, 
+    ruIndex, is_found, is_match);
+    return;
+
+
+    uint32_t j, tn, uId, is_Unitig, offset;
+    (*is_found) = (*is_match) = 0;
+
+    for (j = 0; j < reverse_sources[qn].length; j++)
+    {
+        tn = Get_tn(reverse_sources[qn].buffer[j]);
+        if(read_g->seq[tn].del == 1)
+        {
+            get_R_to_U(ruIndex, tn, &tn, &is_Unitig);
+            if(tn == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[tn].del == 1) continue;
+        }
+
+
+        get_R_to_U(ruIndex, tn, &uId, &is_Unitig);
+        if(uId!=(uint32_t)-1 && is_Unitig == 1 && uId == targetID)
+        {
+            offset = (uint32_t)(position_index[tn]);
+            if(offset >= targetBeg && offset <= targetEnd)
+            {
+                (*is_match)++;
+            }                    
+        }
+        (*is_found)++;
+    }
+
+}
+**/
 
 void determin_hap_alignment_boundary_single_side(uint64_t* readIDs, long long queryLen, long long targetBeg, 
 long long targetEnd, long long targetID, long long eMatch, long long eTotal, long long dir, 
@@ -496,6 +643,7 @@ R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_co
         {
             check_hap_match(readIDs[i]>>33, targetBeg, targetEnd, targetID, position_index, reverse_sources, 
             read_g, ruIndex, &is_found, &is_match);
+
             min_count += is_found;
             max_count += is_match;
             if(max_count > min_count*Hap_rate) maxId = i;
@@ -506,6 +654,8 @@ R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_co
         {
             check_hap_match(readIDs[i]>>33, targetBeg, targetEnd, targetID, position_index, reverse_sources, 
             read_g, ruIndex, &is_found, &is_match);
+
+            ///if(is_found > 0 && is_match > 0 && is_match > is_found*Hap_rate)
             if(is_found == 1 && is_match == 1)
             {
                 break;
@@ -522,6 +672,7 @@ R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_co
         {
             check_hap_match(readIDs[i]>>33, targetBeg, targetEnd, targetID, position_index, reverse_sources, 
             read_g, ruIndex, &is_found, &is_match);
+
             min_count += is_found;
             max_count += is_match;
             if(max_count > min_count*Hap_rate) maxId = i;
@@ -531,6 +682,8 @@ R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_co
         {
             check_hap_match(readIDs[i]>>33, targetBeg, targetEnd, targetID, position_index, reverse_sources, 
             read_g, ruIndex, &is_found, &is_match);
+
+            ///if(is_found > 0 && is_match > 0 && is_match > is_found*Hap_rate)
             if(is_found == 1 && is_match == 1)
             {
                 break;
@@ -765,10 +918,10 @@ uint64_t* position_index, ma_utg_t* xReads, ma_utg_t* yReads)
 
 void get_base_boundary(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, 
 asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads,
-uint32_t xUid, uint32_t yUid, long long begIndex, long long endIndex, uint32_t dir, uint32_t rev,
+uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, long long yEndIndex, 
-uint32_t* x_off, uint32_t* y_off)
+uint32_t dir, uint32_t rev, uint32_t* x_off, uint32_t* y_off)
 {
-    long long k, j;
+    long long k, j, offset;
     ma_hit_t_alloc *xR = NULL;
     ma_hit_t *h = NULL;
     ma_sub_t *sq = NULL, *st = NULL;
@@ -779,7 +932,7 @@ uint32_t* x_off, uint32_t* y_off)
     (*x_off) = (*y_off) = (uint32_t)-1;
     if(dir == 1)
     {
-        for (k = endIndex; k >= begIndex; k--)
+        for (k = xEndIndex; k >= xBegIndex; k--)
         {
             xR = &(reverse_sources[xReads->a[k]>>33]);
             is_found = 0; oLen = 0;
@@ -817,6 +970,11 @@ uint32_t* x_off, uint32_t* y_off)
                 if(rev == 0 && v_dir != w_dir) continue;
                 if(rev == 1 && v_dir == w_dir) continue;
 
+                /****************************may have bugs********************************/
+                offset = (uint32_t)(position_index[rId]);
+                if(offset < yBegIndex || offset > yEndIndex) continue;
+                /****************************may have bugs********************************/
+
                 tmp = get_xy_pos(read_g, &t, v, w, xReads->len, yReads->len, position_index, &(t.el));
                 if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue;
 
@@ -835,7 +993,7 @@ uint32_t* x_off, uint32_t* y_off)
     }
     else
     {
-        for (k = begIndex; k <= endIndex; k++)
+        for (k = xBegIndex; k <= xEndIndex; k++)
         {
             xR = &(reverse_sources[xReads->a[k]>>33]);
             is_found = 0; oLen = 0;
@@ -872,8 +1030,11 @@ uint32_t* x_off, uint32_t* y_off)
                 w_dir = (t.v == w)?1:0;
                 if(rev == 0 && v_dir != w_dir) continue;
                 if(rev == 1 && v_dir == w_dir) continue;
-                
+            
-                
+                /****************************may have bugs********************************/
+                offset = (uint32_t)(position_index[rId]);
+                if(offset < yBegIndex || offset > yEndIndex) continue;
+                /****************************may have bugs********************************/
 
                 tmp = get_xy_pos(read_g, &t, v, w, xReads->len, yReads->len, position_index, &(t.el));
                 if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue;
@@ -905,11 +1066,12 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon
     uint32_t x_pos_beg, y_pos_beg, x_pos_end, y_pos_end;
     /*************************x***************************/
     get_base_boundary(ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, 
-    min_ovlp, xReads, yReads, xUid, yUid, Get_x_beg(*hap_can), Get_x_end(*hap_can), 0, 
+    min_ovlp, xReads, yReads, xUid, yUid, Get_x_beg(*hap_can), Get_x_end(*hap_can), 
-    Get_rev(*hap_can), &x_pos_beg, &y_pos_beg);
+    Get_y_beg(*hap_can), Get_y_end(*hap_can), 0, Get_rev(*hap_can), &x_pos_beg, &y_pos_beg);
+
     get_base_boundary(ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, 
-    min_ovlp, xReads, yReads, xUid, yUid, Get_x_beg(*hap_can), Get_x_end(*hap_can), 1, 
+    min_ovlp, xReads, yReads, xUid, yUid, Get_x_beg(*hap_can), Get_x_end(*hap_can), 
-    Get_rev(*hap_can), &x_pos_end, &y_pos_end);
+    Get_y_beg(*hap_can), Get_y_end(*hap_can), 1, Get_rev(*hap_can), &x_pos_end, &y_pos_end);
     /*************************x***************************/
 
     if(x_pos_beg == (uint32_t)-1 || y_pos_beg == (uint32_t)-1 
@@ -975,6 +1137,7 @@ long long* r_y_pos_beg, long long* r_y_pos_end)
     ///flag = classify_hap_overlap(xBasePos, xBasePos, xReads->len, yBasePos, yBasePos, yReads->len);
     flag = vote_overlap_type(u_buffer, hap_can, position_index, xReads, yReads);
 
+
     if(flag == XCY)
     {
         get_pair_hap_similarity(yReads->a, yLen, xUid, reverse_sources, read_g, ruIndex, 
@@ -1007,7 +1170,6 @@ long long* r_y_pos_beg, long long* r_y_pos_end)
         max_count = xLeftMatch + yRightMatch;
         min_count = xLeftTotal + yRightTotal;
     } else abort();
-    
 
     hap_can->weight = hap_can->index_beg = 0;
     if(min_count == 0) return NON_PLOID;
@@ -1042,6 +1204,8 @@ long long* r_y_pos_beg, long long* r_y_pos_end)
         hap_can->index_end = determine_hap_overlap_type(hap_can, xReads, yReads, ruIndex, 
         reverse_sources, coverage_cut, read_g, position_index, max_hang, min_ovlp, xUid, 
         yUid, r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end);
+        if(hap_can->index_end == XCY && yReads->len > (xReads->len*2)) return NON_PLOID;
+        if(hap_can->index_end == YCX && xReads->len > (yReads->len*2)) return NON_PLOID;
         if(hap_can->index_end == (uint32_t)-1) return NON_PLOID;
 
         return PLOID;
@@ -1053,9 +1217,11 @@ long long* r_y_pos_beg, long long* r_y_pos_end)
 
 void print_hap_paf(ma_ug_t *ug, hap_overlaps* ovlp)
 {
-    fprintf(stderr, "utg%.6d%c\t%u\t%u\t%u\t%c\tutg%.6d%c\t%u\t%u\t%u\t%u\t%u\n", 
+    fprintf(stderr, "utg%.6d%c\t%u(%u)\t%u(%u)\t%u(%u)\t%c\tutg%.6d%c\t%u(%u)\t%u(%u)\t%u(%u)\t%u\t%u\n", 
-    ovlp->xUid+1, "lc"[ug->u.a[ovlp->xUid].circ], ug->u.a[ovlp->xUid].len, ovlp->x_beg_pos, ovlp->x_end_pos, "+-"[ovlp->rev], 
+    ovlp->xUid+1, "lc"[ug->u.a[ovlp->xUid].circ], ug->u.a[ovlp->xUid].len, ug->u.a[ovlp->xUid].n,
-    ovlp->yUid+1, "lc"[ug->u.a[ovlp->yUid].circ], ug->u.a[ovlp->yUid].len, ovlp->y_beg_pos, ovlp->y_end_pos, ovlp->type, (uint32_t)ovlp->weight);
+    ovlp->x_beg_pos, ovlp->x_beg_id, ovlp->x_end_pos, ovlp->x_end_id, "+-"[ovlp->rev], 
+    ovlp->yUid+1, "lc"[ug->u.a[ovlp->yUid].circ], ug->u.a[ovlp->yUid].len, ug->u.a[ovlp->yUid].n,
+    ovlp->y_beg_pos, ovlp->y_beg_id, ovlp->y_end_pos, ovlp->y_end_id, ovlp->type, (uint32_t)ovlp->weight);
 }
 
 void hap_alignment(ma_ug_t *ug, asg_t *read_g,  ma_hit_t_alloc* reverse_sources, 
@@ -1139,18 +1305,291 @@ hap_overlaps_list* all_ovlp)
         xReads = &(ug->u.a[xUid]);
         yReads = &(ug->u.a[yUid]);
         u_buffer->a.n = 0;
-        
 
-        //if(xUid == 404 && yUid == 307)
+        for (k = 0; k < xReads->n; k++)
-        // if(xReads->n <= 100 && yReads->n <= 100)
+        {
-        // {
+            xR = &(reverse_sources[xReads->a[k]>>33]);
-        //     debug_enable = 1;
+            for (j = 0; j < xR->length; j++)
-        // }
+            {
-        // else
+                h = &(xR->buffer[j]);
-        // {
+                sq = &(coverage_cut[Get_qn(*h)]);
-        //     debug_enable = 0;
+                st = &(coverage_cut[Get_tn(*h)]);
-        // }
+                if(st->del || read_g->seq[Get_tn(*h)].del) continue;
-        // debug_enable = 1;
+
+                r = ma_hit2arc(h, sq->e - sq->s, st->e - st->s, max_hang, 
+                                asm_opt.max_hang_rate, min_ovlp, &t);
+                ///if it is a contained overlap, skip
+                if(r < 0) continue;
+
+                rId = t.v>>1;
+                if(read_g->seq[rId].del == 1) continue;
+                ///there are two cases: 
+                ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id  
+                ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1
+                get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig);
+                if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
+                if(Hap_uId != yUid) continue;
+
+                v = xReads->a[k]>>32;
+                get_R_to_U(ruIndex, v>>1, &Hap_uId, &is_Unitig);
+                if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
+                if(Hap_uId != xUid) continue;
+                if((uint32_t)(position_index[v>>1]) != k) continue;
+
+
+                if((prefilter((uint32_t)(position_index[v>>1]), (uint32_t)(position_index[rId]), 
+                    xReads->n, yReads->n, 0, Hap_rate, seedOcc)==NON_PLOID) && 
+                   (prefilter((uint32_t)(position_index[v>>1]), (uint32_t)(position_index[rId]), 
+                    xReads->n, yReads->n, 1, Hap_rate, seedOcc)==NON_PLOID))
+                {
+                    continue;
+                }
+            
+                t_offset.Off = get_xy_pos(read_g, &t, v, (yReads->a[(uint32_t)(position_index[rId])])>>32, 
+                xReads->len, yReads->len, position_index, &(t.el));
+                if(((t_offset.Off>>32) == (uint32_t)-1) || (((uint32_t)t_offset.Off) == (uint32_t)-1)) continue;
+
+                t_offset.x = t;
+                t_offset.weight = 1;
+
+                kv_push(asg_arc_t_offset, u_buffer->a, t_offset);
+            }
+
+            deduplicate_edge(u_buffer);
+        }
+
+        if(u_buffer->a.n == 0) continue;
+        
+        qsort(u_buffer->a.a, u_buffer->a.n, sizeof(asg_arc_t_offset), cmp_hap_alignment);
+        k = 0;
+        u_can->a.n = 0;
+        while (k < u_buffer->a.n)
+        {
+            hap_can.rev = u_buffer->a.a[k].x.el;
+            hap_can.index_beg = k;
+            hap_can.index_end = k;
+            hap_can.weight = u_buffer->a.a[k].weight;
+            hap_can.x_beg_pos = hap_can.x_end_pos = (uint32_t)(u_buffer->a.a[k].Off>>32);
+            hap_can.y_beg_pos = hap_can.y_end_pos = (uint32_t)(u_buffer->a.a[k].Off);
+            cur_offset = Cal_Off(u_buffer->a.a[k].Off);
+            interval_len = get_hap_overlapLen(hap_can.x_beg_pos, hap_can.x_end_pos, xReads->len, 
+            hap_can.y_beg_pos, hap_can.y_end_pos, yReads->len, NULL, NULL, NULL, NULL);
+
+
+            k++;
+            while (k < u_buffer->a.n)
+            {
+                new_offset = Cal_Off(u_buffer->a.a[k].Off);
+                if(u_buffer->a.a[k].x.el != hap_can.rev) break;
+                if((new_offset - cur_offset)>(interval_len*chain_rate)) break;
+                
+
+                hap_can.index_end = k;
+                hap_can.weight += u_buffer->a.a[k].weight;
+
+                is_update = 0;
+                xPos = (uint32_t)(u_buffer->a.a[k].Off>>32);
+                yPos = (uint32_t)(u_buffer->a.a[k].Off);
+                if(xPos < hap_can.x_beg_pos)
+                {
+                    hap_can.x_beg_pos = xPos;
+                    is_update = 1;
+                } 
+
+                if(xPos > hap_can.x_end_pos)
+                {
+                    hap_can.x_end_pos = xPos;
+                    is_update = 1;
+                }
+
+                if(yPos < hap_can.y_beg_pos)
+                {
+                    hap_can.y_beg_pos = yPos;
+                    is_update = 1;
+                }
+
+                if(yPos > hap_can.y_end_pos)
+                {
+                    hap_can.y_end_pos = yPos;
+                    is_update = 1;
+                }
+
+                if(new_offset == cur_offset) is_update = 0;
+
+                if(is_update)
+                {
+                    interval_len = get_hap_overlapLen(hap_can.x_beg_pos, hap_can.x_end_pos, xReads->len, 
+                    hap_can.y_beg_pos, hap_can.y_end_pos, yReads->len, NULL, NULL, NULL, NULL);
+                }
+
+                k++;
+            }
+
+            kv_push(hap_candidates, u_can->a, hap_can);
+        }
+
+        if(u_can->a.n == 0) continue;
+        
+        qsort(u_can->a.a, u_can->a.n, sizeof(hap_candidates), cmp_hap_candidates);
+
+        Get_match(hap_can) = Get_total(hap_can) = 0;
+        memset(&hap_align, 0, sizeof(hap_overlaps));
+        
+        for (k = 0; k < u_can->a.n; k++)
+        {
+            is_update = 0;
+            if(u_can->a.a[k].weight < Get_match(hap_can)*Hap_rate) continue;
+
+            if(calculate_pair_hap_similarity(u_buffer, &(u_can->a.a[k]), position_index, xUid, yUid, 
+            xReads, yReads, reverse_sources, read_g, ruIndex, coverage_cut, Hap_rate, max_hang, 
+            min_ovlp, &r_x_pos_beg, &r_x_pos_end, &r_y_pos_beg, &r_y_pos_end)!=PLOID)
+            {
+                continue;
+            }
+            
+            if(Get_match(hap_can) < Get_match(u_can->a.a[k]))
+            {
+                is_update = 1;
+            }
+            else if(Get_match(hap_can) == Get_match(u_can->a.a[k]) && 
+                                    Get_total(hap_can) > Get_total(u_can->a.a[k]))
+            {
+                is_update = 1;
+            }
+
+            if(is_update)
+            {
+                hap_can = u_can->a.a[k];
+                hap_align.rev = Get_rev(hap_can);
+                hap_align.type = Get_type(hap_can);
+                hap_align.x_beg_id = Get_x_beg(hap_can);
+                hap_align.x_end_id = Get_x_end(hap_can) + 1;
+                hap_align.y_beg_id = Get_y_beg(hap_can);
+                hap_align.y_end_id = Get_y_end(hap_can) + 1;
+                hap_align.weight = Get_match(hap_can);
+                hap_align.x_beg_pos = r_x_pos_beg;
+                hap_align.x_end_pos = r_x_pos_end + 1;
+                if(hap_align.rev == 0)
+                {
+                    hap_align.y_beg_pos = r_y_pos_beg;
+                    hap_align.y_end_pos = r_y_pos_end + 1;
+                }
+                else
+                {
+                    hap_align.y_beg_pos = yReads->len - r_y_pos_end - 1;
+                    hap_align.y_end_pos = yReads->len - r_y_pos_beg - 1 + 1;
+                }
+                hap_align.xUid = xUid;
+                hap_align.yUid = yUid;
+                hap_align.status = SELF_EXIST;
+            }
+        }
+
+        if(Get_match(hap_can) == 0 || Get_total(hap_can) == 0) continue;
+        
+        kv_push(hap_overlaps, all_ovlp->x[hap_align.xUid].a, hap_align);
+    }
+
+}
+
+
+static void hap_alignment_worker(void *_data, long eid, int tid)
+{
+    hap_alignment_struct_pip* hap_buf = (hap_alignment_struct_pip*)_data;
+    ma_ug_t *ug = hap_buf->ug;
+    asg_t *read_g = hap_buf->read_g;
+    ma_hit_t_alloc* reverse_sources = hap_buf->reverse_sources;
+    R_to_U* ruIndex = hap_buf->ruIndex;
+    ma_sub_t *coverage_cut = hap_buf->coverage_cut;
+    uint64_t* position_index = hap_buf->position_index;
+    float Hap_rate = hap_buf->Hap_rate;
+    int max_hang = hap_buf->max_hang;
+    int min_ovlp = hap_buf->min_ovlp;
+    float chain_rate = hap_buf->chain_rate;
+    hap_overlaps_list* all_ovlp = hap_buf->all_ovlp;
+    uint32_t Input_uId = eid;
+    uint64_t* vote_counting = hap_buf->buf[tid].vote_counting;
+    uint8_t* visit = hap_buf->buf[tid].visit;
+    kvec_t_u64_warp* u_vecs = &(hap_buf->buf[tid].u_vecs);
+    kvec_asg_arc_t_offset* u_buffer = &(hap_buf->buf[tid].u_buffer);
+    kvec_hap_candidates* u_can = &(hap_buf->buf[tid].u_can);
+
+
+    ma_utg_t *xReads = NULL, *yReads = NULL;
+    ma_hit_t_alloc *xR = NULL;
+    ma_hit_t *h = NULL;
+    ma_sub_t *sq = NULL, *st = NULL;
+    asg_t* nsg = ug->g;
+    uint32_t i, j, v, rId, k, is_Unitig, Hap_uId, xUid, yUid, seedOcc, xPos, yPos, is_update;
+    uint64_t tmp;
+    long long cur_offset, new_offset, interval_len;
+    long long r_x_pos_beg, r_x_pos_end, r_y_pos_beg, r_y_pos_end;
+    int32_t r;
+    asg_arc_t t;
+    asg_arc_t_offset t_offset;
+    hap_candidates hap_can;
+    hap_overlaps hap_align;
+    xUid = Input_uId;
+    if(nsg->seq[xUid].del || nsg->seq[xUid].c == ALTER_LABLE) return;
+    memset(vote_counting, 0, sizeof(uint64_t)*nsg->n_seq);
+    memset(visit, 0, nsg->n_seq);
+    u_vecs->a.n = 0;
+    u_can->a.n = 0;
+
+    xReads = &(ug->u.a[xUid]);
+    for (i = 0; i < xReads->n; i++)
+    {
+        xR = &(reverse_sources[xReads->a[i]>>33]);
+        
+        for (k = 0; k < xR->length; k++)
+        {
+            rId = Get_tn(xR->buffer[k]);
+        
+            if(read_g->seq[rId].del == 1)
+            {
+                ///get the id of read that contains it 
+                get_R_to_U(ruIndex, rId, &rId, &is_Unitig);
+                if(rId == (uint32_t)-1 || is_Unitig == 1 || read_g->seq[rId].del == 1) continue;
+            }
+
+            ///there are two cases: 
+            ///1. read at primary contigs, get_R_to_U() return its corresponding contig Id  
+            ///2. read at alternative contigs, get_R_to_U() return (uint32_t)-1
+            get_R_to_U(ruIndex, rId, &Hap_uId, &is_Unitig);
+            if(is_Unitig == 0 || Hap_uId == (uint32_t)-1) continue;
+            ///here rId is the id of the read coming from the different haplotype
+            ///Hap_cId is the id of the corresponding contig (note here is the contig, instead of untig)
+            if(visit[Hap_uId]!=0) continue;
+            visit[Hap_uId] = 1;
+            if(vote_counting[Hap_uId] < UINT64_MAX) vote_counting[Hap_uId]++;
+        }
+
+        clean_visit_flag(visit, read_g, ruIndex, nsg->n_seq, xR);
+    }
+
+
+
+    u_vecs->a.n = 0;
+    for (i = 0; i < nsg->n_seq; i++)
+    {
+        if(i == xUid) continue;
+        if(vote_counting[i] == 0) continue;
+        tmp = vote_counting[i]; tmp = tmp << 32; tmp = tmp | (uint64_t)i;
+        kv_push(uint64_t, u_vecs->a, tmp);
+    }
+
+    if(u_vecs->a.n == 0) return;
+    sort_kvec_t_u64_warp(u_vecs, 1);
+
+
+    ///scan each candidate unitig
+    for (i = 0; i < u_vecs->a.n; i++)
+    {
+        yUid = (uint32_t)u_vecs->a.a[i];
+        seedOcc = u_vecs->a.a[i]>>32;
+        xReads = &(ug->u.a[xUid]);
+        yReads = &(ug->u.a[yUid]);
+        u_buffer->a.n = 0;
 
         for (k = 0; k < xReads->n; k++)
         {
@@ -1344,7 +1783,6 @@ hap_overlaps_list* all_ovlp)
 
 }
 
-
 int inline get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn)
 {
     uint32_t i;
@@ -1587,10 +2025,10 @@ void clean_purge_graph(asg_t *purge_g, int max_dist, float drop_ratio)
 
 void get_node_boundary(R_to_U* ruIndex, ma_hit_t_alloc* reverse_sources, ma_sub_t *coverage_cut, 
 asg_t *read_g, uint64_t* position_index, int max_hang, int min_ovlp, ma_utg_t *xReads, ma_utg_t *yReads,
-uint32_t xUid, uint32_t yUid, long long begIndex, long long endIndex, uint32_t dir, uint32_t rev,
+uint32_t xUid, uint32_t yUid, long long xBegIndex, long long xEndIndex, long long yBegIndex, 
-asg_arc_t* reture_t_f, asg_arc_t* reture_t_r)
+long long yEndIndex, uint32_t dir, uint32_t rev, asg_arc_t* reture_t_f, asg_arc_t* reture_t_r)
 {
-    long long k, j;
+    long long k, j, offset;
     ma_hit_t_alloc *xR = NULL;
     ma_hit_t *h = NULL;
     ma_sub_t *sq = NULL, *st = NULL;
@@ -1600,7 +2038,7 @@ asg_arc_t* reture_t_f, asg_arc_t* reture_t_r)
     reture_t_f->del = reture_t_r->del = 1;
     if(dir == 1)
     {
-        for (k = endIndex; k >= begIndex; k--)
+        for (k = xEndIndex; k >= xBegIndex; k--)
         {
             xR = &(reverse_sources[xReads->a[k]>>33]);
             is_found = 0;
@@ -1637,6 +2075,12 @@ asg_arc_t* reture_t_f, asg_arc_t* reture_t_r)
                 if(rev == 0 && v_dir != w_dir) continue;
                 if(rev == 1 && v_dir == w_dir) continue;
                 if(v_dir == 1) continue;
+
+                /****************************may have bugs********************************/
+                offset = (uint32_t)(position_index[rId]);
+                if(offset < yBegIndex || offset > yEndIndex) continue;
+                /****************************may have bugs********************************/
+
                 /************************get reverse edge*************************/
                 index = get_specific_overlap(&(reverse_sources[Get_tn(*h)]), Get_tn(*h), Get_qn(*h));
                 if(index == -1) continue;
@@ -1663,7 +2107,7 @@ asg_arc_t* reture_t_f, asg_arc_t* reture_t_r)
     }
     else
     {
-        for (k = begIndex; k <= endIndex; k++)
+        for (k = xBegIndex; k <= xEndIndex; k++)
         {
             xR = &(reverse_sources[xReads->a[k]>>33]);
             is_found = 0; oLen = 0;
@@ -1701,6 +2145,11 @@ asg_arc_t* reture_t_f, asg_arc_t* reture_t_r)
                 if(rev == 0 && v_dir != w_dir) continue;
                 if(rev == 1 && v_dir == w_dir) continue;
                 if(v_dir == 0) continue;
+
+                /****************************may have bugs********************************/
+                offset = (uint32_t)(position_index[rId]);
+                if(offset < yBegIndex || offset > yEndIndex) continue;
+                /****************************may have bugs********************************/
             
                 /************************get reverse edge*************************/
                 index = get_specific_overlap(&(reverse_sources[Get_tn(*h)]), Get_tn(*h), Get_qn(*h));
@@ -1900,7 +2349,8 @@ uint64_t* position_index, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge,
             if(cut_end < endIndex) endIndex = cut_end;
 
             get_node_boundary(ruIndex, reverse_sources, coverage_cut, read_g, position_index, max_hang, 
-            min_ovlp, xReads, yReads, v>>1, w>>1, begIndex, endIndex, v&1, x->rev, &t_forward, &t_backward);
+            min_ovlp, xReads, yReads, v>>1, w>>1, begIndex, endIndex, x->y_beg_id, x->y_end_id-1, v&1, 
+            x->rev, &t_forward, &t_backward);
             if(t_forward.del || t_backward.del) break;
             
             kv_push(asg_arc_t, edge->a, t_forward);
@@ -2027,41 +2477,57 @@ uint64_t* position_index, int max_hang, int min_ovlp, kvec_asg_arc_t_warp* edge,
     }
     free(b_0.b.a);
 }
+
+void print_all_purge_ovlp(ma_ug_t *ug, hap_overlaps_list* all_ovlp)
+{
+    uint32_t v, uId, i;
+    for (v = 0; v < all_ovlp->num; v++)
+    {
+        uId = v;
+        if(uId != 96 && uId != 272) continue;
+        for (i = 0; i < all_ovlp->x[uId].a.n; i++)
+        {
+            print_hap_paf(ug, &(all_ovlp->x[uId].a.a[i]));
+        }
+    }
+
+}
 void purge_dups(ma_ug_t *ug, asg_t *read_g, ma_sub_t* coverage_cut, ma_hit_t_alloc* reverse_sources, 
 R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t bi_graph_Len, uint32_t long_hap_overlap, 
 float lable_match_rate, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio,
 uint32_t just_contain)
 {
-    fprintf(stderr, "*****************\n");
     asg_t *purge_g = NULL;
     purge_g = asg_init();
-    kvec_t_u64_warp u_vecs;
-    kv_init(u_vecs.a);
     asg_t* nsg = ug->g;
     uint32_t v, rId, uId, i, k, offset;
     ma_utg_t* reads = NULL;
-    uint8_t* visit = NULL;
+
-    visit = (uint8_t*)malloc(sizeof(uint8_t) * nsg->n_seq);
+    // kvec_t_u64_warp u_vecs;
-    ///int flag;
+    // kv_init(u_vecs.a);
-    memset(visit, 0, nsg->n_seq);
+    // uint8_t* visit = NULL;
+    // visit = (uint8_t*)malloc(sizeof(uint8_t) * nsg->n_seq);
+    // memset(visit, 0, nsg->n_seq);
+    // uint64_t* vote_counting = (uint64_t*)malloc(sizeof(uint64_t)*nsg->n_seq);
+    // memset(vote_counting, 0, sizeof(uint64_t)*nsg->n_seq);
+    // kvec_asg_arc_t_offset u_buffer;
+    // kv_init(u_buffer.a);
+    // kvec_hap_candidates u_can;
+    // kv_init(u_can.a);
     uint64_t* position_index = (uint64_t*)malloc(sizeof(uint64_t)*read_g->n_seq);
     memset(position_index, -1, sizeof(uint64_t)*read_g->n_seq);
-    uint64_t* vote_counting = (uint64_t*)malloc(sizeof(uint64_t)*nsg->n_seq);
+    
-    memset(vote_counting, 0, sizeof(uint64_t)*nsg->n_seq);
-    kvec_asg_arc_t_offset u_buffer;
-    kv_init(u_buffer.a);
-    kvec_hap_candidates u_can;
-    kv_init(u_can.a);
     hap_overlaps_list all_ovlp;
     init_hap_overlaps_list(&all_ovlp, nsg->n_seq);
     hap_overlaps_list back_all_ovlp;
     init_hap_overlaps_list(&back_all_ovlp, nsg->n_seq);
     ///uint32_t junk_cov, hap_cov, dip_cov, junk_occ, repeat_occ, single_cov;
-    ma_sub_t* purge_cut = NULL;
-    purge_cut = (ma_sub_t*)malloc(sizeof(ma_sub_t)*nsg->n_seq);
     asg_arc_t t;
     asg_arc_t* p = NULL;
     int r;
+    hap_alignment_struct_pip hap_buf;
+    
+
 
     /****************************may have bugs********************************/
     for (v = 0; v < nsg->n_seq; ++v) 
@@ -2099,24 +2565,29 @@ uint32_t just_contain)
             offset += (uint32_t)reads->a[i];
         }
 
-        purge_cut[uId].del = 0;
-        purge_cut[uId].s = 0;
-        purge_cut[uId].e = offset;
-        purge_cut[uId].c = PRIMARY_LABLE;
         asg_seq_set(purge_g, uId, offset, 0);
         purge_g->seq[uId].c = PRIMARY_LABLE;
     }
 
 
-    for (v = 0; v < nsg->n_seq; v++)
+    init_hap_alignment_struct_pip(&hap_buf, asm_opt.thread_num, nsg->n_seq, ug, read_g,  
-    {
+    reverse_sources, ruIndex, coverage_cut, position_index, density, max_hang, min_ovlp, 
-        uId = v;
+    0.05, &all_ovlp);
-        if(nsg->seq[uId].del || nsg->seq[uId].c == ALTER_LABLE) continue;
 
-        hap_alignment(ug, read_g,  reverse_sources, ruIndex, coverage_cut, position_index, 
+    kt_for(asm_opt.thread_num, hap_alignment_worker, &hap_buf, nsg->n_seq);
-        vote_counting, visit, &u_vecs, &u_buffer, &u_can, uId, density, max_hang, min_ovlp, 
+
-        0.05, &all_ovlp);
+    ///if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp);
-    }
+    
+
+    // for (v = 0; v < nsg->n_seq; v++)
+    // {
+    //     uId = v;
+    //     if(nsg->seq[uId].del || nsg->seq[uId].c == ALTER_LABLE) continue;
+
+    //     hap_alignment(ug, read_g,  reverse_sources, ruIndex, coverage_cut, position_index, 
+    //     vote_counting, visit, &u_vecs, &u_buffer, &u_can, uId, density, max_hang, min_ovlp, 
+    //     0.05, &all_ovlp);
+    // }
 
     normalize_hap_overlaps(&all_ovlp, &back_all_ovlp);
     ///debug_hap_overlaps(&all_ovlp, &back_all_ovlp);
@@ -2186,9 +2657,12 @@ uint32_t just_contain)
 		asg_symm(purge_g);
 
         clean_purge_graph(purge_g, bubble_dist, drop_ratio);
-        ///if(debug_enable) print_purge_gfa(ug, purge_g);
+
+        // if(debug_enable) print_purge_gfa(ug, purge_g);
+        // if(debug_enable) print_all_purge_ovlp(ug, &all_ovlp);
+
         link_unitigs(purge_g, ug, &all_ovlp, ruIndex, reverse_sources, coverage_cut, read_g, position_index, 
-        max_hang, min_ovlp, edge, visit);
+        max_hang, min_ovlp, edge, hap_buf.buf[0].visit);
     }
 
     for (v = 0; v < all_ovlp.num; v++)
@@ -2208,16 +2682,15 @@ uint32_t just_contain)
     }
 
     asg_cleanup(nsg);
-    kv_destroy(u_vecs.a);
-    kv_destroy(u_buffer.a);
-    kv_destroy(u_can.a);
     destory_hap_overlaps_list(&all_ovlp);
     destory_hap_overlaps_list(&back_all_ovlp);
     asg_destroy(purge_g);
-    free(purge_cut);
     free(position_index);
-    free(vote_counting);
+    // kv_destroy(u_vecs.a);
-    free(visit);
+    // kv_destroy(u_buffer.a);
-    fprintf(stderr, "#################\n");
+    // kv_destroy(u_can.a);
+    // free(vote_counting);
+    // free(visit);
+    destory_hap_alignment_struct_pip(&hap_buf);
 }
 

Purge_Dups.h

@@ -14,6 +14,6 @@ uint32_t just_contain);
 void fill_unitig(uint64_t* buffer, uint32_t bufferLen, asg_t* read_g, kvec_asg_arc_t_warp* edge,
 uint32_t is_circle, uint64_t* rLen);
 
-void enable_debug_mode();
+void enable_debug_mode(uint32_t mode);
 
 #endif