backup

Correct.cpp

@@ -539,6 +539,236 @@ char* r_string)
     }
 }
 
+
+void verify_ul_window(long long window_start, long long window_end, overlap_region_alloc* overlap_list,Correct_dumy* dumy, const ma_ug_t *uref,
+char* r_string)
+{
+    long long i;
+    long long currentID;
+    long long x_start, y_start, o_len;
+    long long Window_Len = WINDOW + (THRESHOLD << 1);
+    char* x_string = NULL;
+    char* y_string = NULL;
+    long long x_end, x_len;
+    int end_site;
+    unsigned int error;
+    int groupLen = 0;
+    int return_sites[GROUP_SIZE];
+	unsigned int return_sites_error[GROUP_SIZE];
+    uint64_t overlapID[GROUP_SIZE];
+    uint64_t y_startGroup[GROUP_SIZE];
+    int y_extra_begin[GROUP_SIZE];
+    int y_extra_end[GROUP_SIZE];
+    int error_threshold[GROUP_SIZE];
+    int extra_begin;
+    int extra_end;
+
+    ///here are overlaps fully covered by WINDOW
+    for (i = 0; i < (long long)dumy->length; i++)
+    {
+        extra_begin = extra_end = 0;
+        ///if the window has been fully covered, the interval at x is [window_start, window_end]
+        x_len = WINDOW;
+        currentID = dumy->overlapID[i];
+        x_start = window_start;
+        ///offset of y
+        y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s;
+        /****************************may have bugs********************************/
+        y_start += y_start_offset(x_start, &(overlap_list->list[currentID].f_cigar));
+        /****************************may have bugs********************************/
+        
+  
+        if(!determine_overlap_region(THRESHOLD, y_start, overlap_list->list[currentID].y_id, Window_Len, R_INF,
+        &extra_begin, &extra_end, &y_start, &o_len))
+        {
+            append_window_list(&overlap_list->list[currentID], window_start, window_end, 
+            -1, -1, -1, -1, -1, -1);
+            continue;
+        }
+
+        fill_subregion(dumy->overlap_region_group[groupLen], y_start, o_len, overlap_list->list[currentID].y_pos_strand, 
+        R_INF, overlap_list->list[currentID].y_id, extra_begin, extra_end);
+
+        y_extra_begin[groupLen] = extra_begin;
+        y_extra_end[groupLen] = extra_end;
+        overlapID[groupLen] = currentID;
+        y_startGroup[groupLen] = y_start;
+        error_threshold[groupLen] = THRESHOLD;
+        x_string = r_string + x_start;
+        groupLen++;
+        
+
+        if (groupLen == GROUP_SIZE)
+        {
+            Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1], 
+            dumy->overlap_region_group[2], dumy->overlap_region_group[3], Window_Len, x_string, WINDOW,
+                return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE);
+            groupLen = 0;
+
+            
+            if (return_sites_error[0]!=(unsigned int)-1)
+            {
+                overlap_list->list[overlapID[0]].align_length += x_len;
+
+                append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, 
+                            y_startGroup[0], y_startGroup[0] + return_sites[0], (int)return_sites_error[0],
+                            y_extra_begin[0], y_extra_end[0], error_threshold[0]);
+            }
+            else
+            {
+                append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], -1, -1,
+                y_extra_begin[0], y_extra_end[0], error_threshold[0]);
+            }
+            
+
+            if (return_sites_error[1]!=(unsigned int)-1)
+            {
+                overlap_list->list[overlapID[1]].align_length += x_len;
+
+                append_window_list(&overlap_list->list[overlapID[1]], window_start, window_end, 
+                            y_startGroup[1], y_startGroup[1] + return_sites[1], (int)return_sites_error[1],
+                            y_extra_begin[1], y_extra_end[1], error_threshold[1]);
+            }
+            else
+            {
+                append_window_list(&overlap_list->list[overlapID[1]], window_start, window_end, y_startGroup[1], -1, -1,
+                y_extra_begin[1], y_extra_end[1], error_threshold[1]);
+            }
+            
+
+            if (return_sites_error[2]!=(unsigned int)-1)
+            {
+                overlap_list->list[overlapID[2]].align_length += x_len;
+
+                append_window_list(&overlap_list->list[overlapID[2]], window_start, window_end, 
+                            y_startGroup[2], y_startGroup[2] + return_sites[2], (int)return_sites_error[2],
+                            y_extra_begin[2], y_extra_end[2], error_threshold[2]);
+            }
+            else
+            {
+                append_window_list(&overlap_list->list[overlapID[2]], window_start, window_end, y_startGroup[2], -1, -1,
+                y_extra_begin[2], y_extra_end[2], error_threshold[2]);
+            }
+            
+
+            if (return_sites_error[3]!=(unsigned int)-1)
+            {
+                overlap_list->list[overlapID[3]].align_length += x_len;
+
+                append_window_list(&overlap_list->list[overlapID[3]], window_start, window_end, 
+                            y_startGroup[3], y_startGroup[3] + return_sites[3], (int)return_sites_error[3],
+                            y_extra_begin[3], y_extra_end[3], error_threshold[3]);
+            }
+            else
+            {
+                append_window_list(&overlap_list->list[overlapID[3]], window_start, window_end, y_startGroup[3], -1, -1,
+                y_extra_begin[3], y_extra_end[3], error_threshold[3]);
+            }    
+        }
+    }
+
+    if (groupLen == 1)
+    {
+        end_site = Reserve_Banded_BPM(dumy->overlap_region_group[0], Window_Len, x_string, WINDOW, THRESHOLD, &error);
+
+        if (error!=(unsigned int)-1)
+        {
+            overlap_list->list[overlapID[0]].align_length += x_len;
+
+            append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, 
+                                y_startGroup[0], y_startGroup[0] + end_site, (int)error,
+                                y_extra_begin[0], y_extra_end[0], error_threshold[0]);
+        }
+        else
+        {
+            append_window_list(&overlap_list->list[overlapID[0]], window_start, window_end, y_startGroup[0], -1, -1,
+            y_extra_begin[0], y_extra_end[0], error_threshold[0]);
+        }
+    }
+    else if (groupLen > 1)
+    {
+        Reserve_Banded_BPM_4_SSE_only(dumy->overlap_region_group[0], dumy->overlap_region_group[1], 
+                dumy->overlap_region_group[2], dumy->overlap_region_group[3], Window_Len, x_string, WINDOW,
+			        return_sites, return_sites_error, THRESHOLD, dumy->Peq_SSE);
+
+        for (i = 0; i < groupLen; i++)
+        {
+            if (return_sites_error[i]!=(unsigned int)-1)
+            {
+                overlap_list->list[overlapID[i]].align_length += x_len;
+                append_window_list(&overlap_list->list[overlapID[i]], window_start, window_end, 
+                                y_startGroup[i], y_startGroup[i] + return_sites[i], (int)return_sites_error[i],
+                                y_extra_begin[i], y_extra_end[i], error_threshold[i]);
+            }
+            else
+            {
+                append_window_list(&overlap_list->list[overlapID[i]], window_start, window_end, y_startGroup[i], -1, -1,
+                y_extra_begin[i], y_extra_end[i], error_threshold[i]);
+            }
+            
+        }
+
+        groupLen = 0;
+    }
+
+    long long reverse_i = dumy->size - 1;
+    int threshold;
+    
+    ///here are overlaps partially covered by WINDOW
+    for (i = 0; i < (long long)dumy->lengthNT; i++)
+    {
+        extra_begin = extra_end = 0;
+        currentID = dumy->overlapID[reverse_i--];
+        x_start = MAX(window_start, (long long)overlap_list->list[currentID].x_pos_s);
+        x_end = MIN(window_end, (long long)overlap_list->list[currentID].x_pos_e);
+
+        ///overlap length between [window_start, window_end]
+        x_len = x_end - x_start + 1;
+        threshold = x_len * asm_opt.max_ov_diff_ec;
+        /****************************may have bugs********************************/
+        threshold = Adjust_Threshold(threshold, x_len);
+        /****************************may have bugs********************************/
+        
+        ///offset of y
+        y_start = (x_start - overlap_list->list[currentID].x_pos_s) + overlap_list->list[currentID].y_pos_s;
+        /****************************may have bugs********************************/
+        y_start += y_start_offset(x_start, &(overlap_list->list[currentID].f_cigar));
+        /****************************may have bugs********************************/
+
+        Window_Len = x_len + (threshold << 1);
+
+        if(!determine_overlap_region(threshold, y_start, overlap_list->list[currentID].y_id, Window_Len, R_INF,
+        &extra_begin, &extra_end, &y_start, &o_len))
+        {
+            append_window_list(&overlap_list->list[currentID], x_start, x_end, 
+            -1, -1, -1, -1, -1, -1);
+            continue;
+        }
+        
+        fill_subregion(dumy->overlap_region, y_start, o_len, overlap_list->list[currentID].y_pos_strand, 
+        R_INF, overlap_list->list[currentID].y_id, extra_begin, extra_end);        
+
+        x_string = r_string + x_start;
+        y_string = dumy->overlap_region;
+
+        end_site = Reserve_Banded_BPM(y_string, Window_Len, x_string, x_len, threshold, &error);
+
+        if (error!=(unsigned int)-1)
+        {
+            overlap_list->list[currentID].align_length += x_len;
+            append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, y_start + end_site, (int)error,
+            extra_begin, extra_end, threshold);
+        }
+        else
+        {
+            append_window_list(&overlap_list->list[currentID], x_start, x_end, y_start, -1, -1,
+            extra_begin, extra_end, threshold);
+        }
+    }
+}
+
+
+
 inline double trim_error_rate(overlap_region_alloc* overlap_list, long long ID)
 {
     long long tLen, tError,i, subWinLen, subWinNum;
@@ -7745,6 +7975,63 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF,
 }
 
 
+void correct_ul_overlap(overlap_region_alloc* overlap_list, const ma_ug_t *uref, 
+                        UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, 
+                        Graph* g, Graph* DAGCon, Cigar_record* current_cigar, 
+                        haplotype_evdience_alloc* hap, Round2_alignment* second_round, 
+                        int force_repeat, int is_consensus, int* fully_cov, int* abnormal, double max_ov_diff_ec)
+{
+    
+    clear_Correct_dumy(dumy, overlap_list);
+
+    long long window_start, window_end;
+
+    Window_Pool w_inf;
+
+    init_Window_Pool(&w_inf, g_read->length, WINDOW, (int)(1.0/max_ov_diff_ec));
+    
+    int flag = 0;
+    
+    while(get_Window(&w_inf, &window_start, &window_end) && flag != -2)
+    {
+        dumy->length = 0;
+        dumy->lengthNT = 0;
+        flag = get_interval(window_start, window_end, overlap_list, dumy);
+
+        switch (flag)
+        {
+            case 1:    ///no match here
+                break;
+            case 0:    ///no match here
+                break;
+            case -2: ///if flag == -2, loop would be terminated
+                break;
+        }
+
+        
+        ///dumy->lengthNT represent how many overlaps that the length of them is not equal to WINDOW; may larger or less than WINDOW
+        ///dumy->length represent how many overlaps that the length of them is WINDOW
+        ///now the windows which are larger than WINDOW are verified one-by-one, to improve it, we can do it group-bygroup
+        verify_ul_window(window_start, window_end, overlap_list, dumy, uref, g_read->seq);
+    }
+
+    /**
+    // recalcate_window(overlap_list, R_INF, g_read, dumy, overlap_read);
+    // partition_overlaps(overlap_list, R_INF, g_read, dumy, hap, force_repeat);
+    recalcate_window_advance(overlap_list, R_INF, g_read, dumy, overlap_read);
+    partition_overlaps_advance(overlap_list, R_INF, g_read, overlap_read, dumy, hap, force_repeat);
+    
+    if(is_consensus)
+    {
+        generate_consensus(overlap_list, R_INF, g_read, dumy, g, DAGCon, current_cigar, second_round);
+    }
+
+
+    (*fully_cov) = check_if_fully_covered(overlap_list, R_INF, g_read, dumy, g, abnormal);
+    **/
+}
+
+
 void init_Cigar_record(Cigar_record* dummy)
 {
     dummy->length = 0;

Correct.h

@@ -1106,6 +1106,12 @@ void add_new_cell_to_cigar_record_with_different_base(Cigar_record* dummy, uint3
 void add_existing_cell_to_cigar_record_with_different_base(Cigar_record* dummy, uint32_t len, uint32_t type, char* seq);
 
 
+void correct_ul_overlap(overlap_region_alloc* overlap_list, const ma_ug_t *uref, 
+                        UC_Read* g_read, Correct_dumy* dumy, UC_Read* overlap_read, 
+                        Graph* g, Graph* DAGCon, Cigar_record* current_cigar, 
+                        haplotype_evdience_alloc* hap, Round2_alignment* second_round, 
+                        int force_repeat, int is_consensus, int* fully_cov, int* abnormal, double max_ov_diff_ec);
+
 /***
  type:
  0. match

Process_Read.cpp

@@ -754,10 +754,12 @@ void init_all_ul_t(all_ul_t *x, All_reads *hR) {
 void destory_all_ul_t(all_ul_t *x) {
 	uint64_t i;
 	for (i = 0; i < x->n; i++) {
-		free(x->a[i].n_n); free(x->a[i].N_site.a); 
-		free(x->a[i].r_base.a); free(x->a[i].bb.a);
+		free(x->a[i].N_site.a); free(x->a[i].r_base.a); free(x->a[i].bb.a);
 	}
 	free(x->a);
+
+	for (i = 0; i < x->nid.n; i++) free(x->nid.a[i].a);
+	free(x->nid.a);
 }
 
 void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn, uint64_t nn_offset)
@@ -843,23 +845,29 @@ void push_subblock_original_bases(char* str, all_ul_t *x, ul_vec_t *p, uint32_t
 void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str, int64_t str_l, ul_ov_t *o, int64_t on) {
 	int64_t i, mine, maxs, ovlp, end;
 	ul_vec_t *p = NULL;
+	nid_t *np = NULL;
 	ul_ov_t *z = NULL, *zp = NULL;
 	uc_block_t *b = NULL;
-	if(rid) fprintf(stderr, "rid:%lu\n", *rid);
-	if(rid == NULL) {
-		kv_pushp(ul_vec_t, *x, &p);
-		memset(p, 0, sizeof(*p));
-		fprintf(stderr, "x->n:%u\n", x->n);
-	}
-	else {
-		p = &(x->a[*rid]);
+
+	if(id) {
+		kv_pushp(nid_t, x->nid, &np);
+		np->n = id_l; MALLOC(np->a, np->n+1); memcpy(np->a, id, id_l); np->a[id_l] = '\0';
 	}
 
-	if(id && id_l > 0) {
-		free(p->n_n);
-		p->n_l = id_l; MALLOC(p->n_n, p->n_l+1); memcpy(p->n_n, id, id_l); p->n_n[id_l] = '\0';
-	}
-	if(str && str_l > 0) {
+	if(str) {
+		if(rid == NULL) {
+			kv_pushp(ul_vec_t, *x, &p);
+			memset(p, 0, sizeof(*p));
+		} else {
+			if((*rid) >= x->m) kv_resize(ul_vec_t, *x, (*rid) + 1);
+			if((*rid) >= x->n) {
+				memset(x->a+x->n, 0, sizeof(*p)*((*rid) + 1 - x->n));
+				x->n = (*rid) + 1;
+			}
+			p = &(x->a[(*rid)]);
+		}
+		
+
 		p->bb.n = p->N_site.n = p->r_base.n = 0;
 		p->rlen = str_l; 
 		

Process_Read.h

@@ -171,9 +171,11 @@ typedef struct
 
 typedef struct
 {
-    char *n_n;
-    uint32_t n_l;
+    char *a; uint32_t n;
+} nid_t;
 
+typedef struct
+{
     kvec_t(uint8_t) r_base;
     uint32_t rlen;
 
@@ -183,6 +185,7 @@ typedef struct
 
 typedef struct
 {
+    kvec_t(nid_t) nid;
     ul_vec_t *a;
     size_t n, m;
     All_reads *hR;

anchor.cpp

@@ -199,6 +199,133 @@ void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_reg
 	///ks_introsort_or_xs(overlap_list->length, overlap_list->list);
 }
 
+void ha_get_new_ul_candidates(ha_abufl_t *ab, int64_t rid, char* rs, int64_t rl, uint64_t mz_w, uint64_t mz_k, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain,
+						   kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, void *ha_flt_tab, ha_pt_t *ha_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp)
+{
+	uint32_t i;
+	uint64_t k, l;
+	uint32_t high_occ = asm_opt.hom_cov >= 1?asm_opt.hom_cov:1;
+	
+
+	// prepare
+    clear_Candidates_list(cl);
+    clear_overlap_region_alloc(overlap_list);
+	ab->mz.n = 0, ab->n_a = 0;
+
+	// get the list of anchors
+	mz2_ha_sketch(rs, rl, mz_w, mz_k, 0, !(asm_opt.flag & HA_F_NO_HPC), &ab->mz, ha_flt_tab, asm_opt.mz_sample_dist, k_flag, dbg_ct, NULL, -1, asm_opt.dp_min_len, -1, sp, asm_opt.mz_rewin, 0);
+	
+	// minimizer of queried read
+	if (ab->mz.m > ab->old_mz_m) {
+		ab->old_mz_m = ab->mz.m;
+		REALLOC(ab->seed, ab->old_mz_m);
+	}
+	for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) {
+		int n;
+		ab->seed[i].a = ha_ptl_get(ha_idx, ab->mz.a[i].x, &n);
+		ab->seed[i].n = n;
+		ab->n_a += n;
+	}
+	if (ab->n_a > ab->m_a) {
+		ab->m_a = ab->n_a;
+		kroundup64(ab->m_a);
+		REALLOC(ab->a, ab->m_a);
+	}
+	for (i = 0, k = 0; i < ab->mz.n; ++i) {
+		int j;
+		///z is one of the minimizer
+		ha_mzl_t *z = &ab->mz.a[i];
+		seedl_t *s = &ab->seed[i];
+		for (j = 0; j < s->n; ++j) {
+			const ha_idxposl_t *y = &s->a[j];
+			anchor1_t *an = &ab->a[k++];
+			uint8_t rev = z->rev == y->rev? 0 : 1;
+			an->other_off = y->pos;
+			an->self_off = rev? rl - 1 - (z->pos + 1 - z->span) : z->pos;
+			an->cnt = s->n;
+			an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off;
+		}
+	}
+	
+	// sort anchors
+	radix_sort_ha_an1(ab->a, ab->a + ab->n_a);
+	for (k = 1, l = 0; k <= ab->n_a; ++k) {
+		if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) {
+			if (k - l > 1)
+				radix_sort_ha_an2(ab->a + l, ab->a + k);
+			l = k;
+		}
+	}
+
+	
+	// copy over to _cl_
+	if (ab->m_a >= (uint64_t)cl->size) {
+		cl->size = ab->m_a;
+		REALLOC(cl->list, cl->size);
+	}
+	for (k = 0; k < ab->n_a; ++k) {
+		k_mer_hit *p = &cl->list[k];
+		p->readID = ab->a[k].srt >> 33;
+		p->strand = ab->a[k].srt >> 32 & 1;
+		p->offset = ab->a[k].other_off;
+		p->self_offset = ab->a[k].self_off;
+		if(ab->a[k].cnt <= high_occ){
+			p->cnt = 1;
+		}
+		else{
+			p->cnt = 1 + ((ab->a[k].cnt + (high_occ<<1) - 1)/(high_occ<<1));
+			p->cnt = pow(p->cnt, 1.1);
+		}
+	}
+	cl->length = ab->n_a;
+
+	calculate_overlap_region_by_chaining(cl, overlap_list, chain_idx, rid, rl, &R_INF, bw_thres, keep_whole_chain, f_cigar);
+	
+	#if 0
+	if (overlap_list->length > 0) {
+		fprintf(stderr, "B\t%ld\t%ld\t%d\n", (long)rid, (long)overlap_list->length, rlen);
+		for (int i = 0; i < (int)overlap_list->length; ++i) {
+			overlap_region *r = &overlap_list->list[i];
+			fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand],
+					(int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed, ha_ov_type(r, rlen));
+		}
+	}
+	#endif
+	
+	if ((int)overlap_list->length > max_n_chain) {
+		int32_t w, n[4], s[4];
+		n[0] = n[1] = n[2] = n[3] = 0, s[0] = s[1] = s[2] = s[3] = 0;
+		ks_introsort_or_ss(overlap_list->length, overlap_list->list);
+		for (i = 0; i < (uint32_t)overlap_list->length; ++i) {
+			const overlap_region *r = &overlap_list->list[i];
+			w = ha_ov_type(r, rl);
+			++n[w];
+			if ((int)n[w] == max_n_chain) s[w] = r->shared_seed;
+		}
+		if (s[0] > 0 || s[1] > 0 || s[2] > 0 || s[3] > 0) {
+			// n[0] = n[1] = n[2] = n[3] = 0;
+			for (i = 0, k = 0; i < (uint32_t)overlap_list->length; ++i) {
+				overlap_region *r = &overlap_list->list[i];
+				w = ha_ov_type(r, rl);
+				// ++n[w];
+				// if (((int)n[w] <= max_n_chain) || (r->shared_seed >= s[w] && s[w] >= (asm_opt.k_mer_length<<1))) {
+				if (r->shared_seed >= s[w]) {
+					if ((uint32_t)k != i) {
+						overlap_region t;
+						t = overlap_list->list[k];
+						overlap_list->list[k] = overlap_list->list[i];
+						overlap_list->list[i] = t;
+					}
+					++k;
+				}
+			}
+			overlap_list->length = k;
+		}
+	}
+
+	///ks_introsort_or_xs(overlap_list->length, overlap_list->list);
+}
+
 
 void calculate_ug_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list, kvec_t_u64_warp* chain_idx,
                                           uint64_t readID, ma_utg_v *ua, double band_width_threshold, int add_beg_end, overlap_region* f_cigar, long long mz_occ, double mz_rate)
@@ -626,6 +753,19 @@ void ha_get_candidates_interface(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overl
 	ks_introsort_or_xs(overlap_list->length, overlap_list->list);
 }
 
+
+void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, 
+								 int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp)
+{
+	extern void *ha_flt_tab;
+	extern ha_pt_t *ha_idx;
+
+	ha_get_new_ul_candidates(ab, rid, rs, rl, mz_w, mz_k, overlap_list, cl, bw_thres, max_n_chain, keep_whole_chain, k_flag, chain_idx, ha_flt_tab, ha_idx, f_cigar, dbg_ct, sp);
+	ks_introsort_or_xs(overlap_list->length, overlap_list->list);
+}
+
+
+
 void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list)
 {
 	ks_introsort_or_xs(overlap_list->length, overlap_list->list);

htab.cpp

@@ -1141,11 +1141,11 @@ ha_pt_t *ha_pt_ul_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, ma_utg_
 	///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen
 	if (flt_tab == 0) {
         if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
-        ha_ct_shrink(ct, 2, cutoff, asm_opt->thread_num);
+        ha_ct_shrink(ct, /**2**/1, cutoff, asm_opt->thread_num);
     } else {
         ///Note: here is just to remove minimizer appearing YAK_MAX_COUNT times
         ///minimizer with YAK_MAX_COUNT occ may apper > YAK_MAX_COUNT times, so it may lead to overflow at ha_pt_gen
-        ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num);
+        ha_ct_shrink(ct, /**2**/1, YAK_MAX_COUNT - 1, asm_opt->thread_num);
     }
 	
 	pt = ha_pt_gen(ct, asm_opt->thread_num, 1);

inter.cpp

@@ -18,6 +18,9 @@
 #include "Assembly.h"
 KSEQ_INIT(gzFile, gzread)
 
+void ha_get_ul_candidates_interface(ha_abufl_t *ab, int64_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, overlap_region_alloc *overlap_list, overlap_region_alloc *overlap_list_hp, Candidates_list *cl, double bw_thres, 
+								 int max_n_chain, int keep_whole_chain, kvec_t_u8_warp* k_flag, kvec_t_u64_warp* chain_idx, overlap_region* f_cigar, kvec_t_u64_warp* dbg_ct, st_mt_t *sp);
+
 #define MG_SEED_IGNORE     (1ULL<<41)
 #define MG_SEED_TANDEM     (1ULL<<42)
 #define MG_SEED_KEPT       (1ULL<<43)
@@ -68,6 +71,8 @@ typedef struct {
 	int max_lc_skip, max_lc_iter, min_lc_cnt, min_lc_score, max_gc_skip, ref_bonus;
 	int min_gc_cnt, min_gc_score, sub_diff, best_n;
     float chn_pen_gap, mask_level, pri_ratio;
+	///base-alignment
+	double bw_thres, diff_ec_ul; int max_n_chain;
 } mg_idxopt_t;
 
 typedef struct {
@@ -162,6 +167,7 @@ typedef struct { // global data structure for kt_pipeline()
     uint64_t total_pair;
 	mg_gres_a hits;
 	mg_dbn_t nn;
+	uint64_t num_bases, num_corrected_bases, num_recorrected_bases;
 } uldat_t;
 
 typedef struct {
@@ -293,9 +299,10 @@ typedef struct { // data structure for each step in kt_pipeline()
 	mg_gchains_t **gcs;///useless
     mg_tbuf_t **buf;///useless
 	ha_ovec_buf_t **hab;
+	uint64_t num_bases, num_corrected_bases, num_recorrected_bases;
 } utepdat_t;
 
-void init_mg_opt(mg_idxopt_t *opt, int is_HPC, int k, int w, int hap_n)
+void init_mg_opt(mg_idxopt_t *opt, int is_HPC, int k, int w, int hap_n, int max_n_chain, double bw_thres, double diff_ec_ul)
 {
     opt->k = k; 
     opt->w = w;
@@ -320,6 +327,9 @@ void init_mg_opt(mg_idxopt_t *opt, int is_HPC, int k, int w, int hap_n)
 	opt->sub_diff = 6;
 	opt->best_n = 5;
 	opt->pri_ratio = 0.8f;
+	opt->max_n_chain = max_n_chain;
+	opt->bw_thres = bw_thres;
+	opt->diff_ec_ul = diff_ec_ul;
 }
 
 void uidx_l_build(ma_ug_t *ug, mg_idxopt_t *opt, int cutoff)
@@ -2116,6 +2126,30 @@ static void worker_for_ul_alignment(void *data, long i, int tid) // callback for
     s->opt->is_HPC, asm_opt.mz_sample_dist, asm_opt.mz_rewin, s->opt, s->uopt, &(s->gcs[i]));
 }
 
+static void worker_for_ul_scall_alignment(void *data, long i, int tid) // callback for kt_for()
+{
+    utepdat_t *s = (utepdat_t*)data;
+	ha_ovec_buf_t *b = s->hab[tid];
+	int64_t rid = s->id+i;
+	int fully_cov, abnormal;
+	
+	ha_get_ul_candidates_interface(b->abl, rid, s->seq[i], s->len[i], s->opt->w, s->opt->k, &b->olist, &b->olist_hp, &b->clist, s->opt->bw_thres, 
+		s->opt->max_n_chain, 1, &(b->k_flag), &b->r_buf, &(b->tmp_region), NULL, &(b->sp));
+    // mg_map_frag(s->ha_flt_tab, s->ha_idx, s->ug, s->rg, s->id+i, s->len[i], s->seq[i], &(s->mzs[tid]), &(s->sps[tid]), s->buf[tid], s->opt->w, s->opt->k, 
+    // s->opt->is_HPC, asm_opt.mz_sample_dist, asm_opt.mz_rewin, s->opt, s->uopt, &(s->gcs[i]));
+	clear_Cigar_record(&b->cigar1);
+	clear_Round2_alignment(&b->round2);
+
+	b->self_read.seq = s->seq[i]; b->self_read.length = s->len[i]; b->self_read.size = 0;
+	correct_ul_overlap(&b->olist, s->ug, &b->self_read, &b->correct, &b->ovlp_read, &b->POA_Graph, &b->DAGCon,
+			&b->cigar1, &b->hap, &b->round2, 0, 1, &fully_cov, &abnormal, s->opt->diff_ec_ul);
+
+	b->num_read_base += b->self_read.length;
+	b->num_correct_base += b->correct.corrected_base;
+	b->num_recorrect_base += b->round2.dumy.corrected_base;
+	memset(&b->self_read, 0, sizeof(b->self_read));
+}
+
 void dump_gaf(mg_gres_a *hits, const mg_gchains_t *gs, uint32_t only_p)
 {
 	if (gs == NULL || gs->n_gc == 0 || gs->n_lc == 0) return;
@@ -2297,6 +2331,7 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac
 		uint64_t i;
 		CALLOC(s->hab, p->n_thread);
 		for (i = 0; i < p->n_thread; ++i) s->hab[i] = ha_ovec_init(0, 0, 1);
+		kt_for(p->n_thread, worker_for_ul_scall_alignment, s, s->n);
 		///debug
 		/**
 		uint64_t i;
@@ -2319,7 +2354,12 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac
 			free(s->mzs[i].a); free(s->sps[i].a);
 		}
 		**/
-		for (i = 0; i < p->n_thread; ++i) ha_ovec_destroy(s->hab[i]);
+		for (i = 0; i < p->n_thread; ++i) {
+			s->num_bases += s->hab[i]->num_read_base;
+			s->num_corrected_bases += s->hab[i]->num_correct_base;
+			s->num_recorrected_bases += s->hab[i]->num_recorrect_base;
+			ha_ovec_destroy(s->hab[i]);
+		}
 		free(s->hab);
 		// free(s->buf); free(s->mzs); free(s->sps);
 		return s;
@@ -2327,6 +2367,9 @@ static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callbac
     else if (step == 2) { // step 3: dump
         utepdat_t *s = (utepdat_t*)in;
 		uint64_t i, rid;
+		p->num_bases += s->num_bases;
+		p->num_corrected_bases += s->num_corrected_bases;
+		p->num_recorrected_bases += s->num_recorrected_bases;
         for (i = 0; i < (uint64_t)s->n; ++i) {
 			///debug
 			/**
@@ -2367,8 +2410,9 @@ int scall_ul_pipeline(uldat_t* sl, const enzyme *fn)
 	sl->hits.total_base = sl->total_base;
 	sl->hits.total_pair = sl->total_pair;
     fprintf(stderr, "[M::%s::%.3f] ==> Qualification\n", __func__, yak_realtime()-index_time);
-	fprintf(stderr, "[M::%s::%.3f] ==> # reads: %lu, # bases: %lu\n", __func__, yak_realtime()-index_time, 
-							UL_INF.n, sl->total_base);
+	fprintf(stderr, "[M::%s::] ==> # reads: %lu, # bases: %lu\n", __func__, UL_INF.n, sl->total_base);
+	fprintf(stderr, "[M::%s::] ==> # bases: %lu; # corrected bases: %lu; # recorrected bases: %lu\n", 
+	__func__, sl->num_bases, sl->num_corrected_bases, sl->num_recorrected_bases);
 	
     return 1;
 }
@@ -2377,13 +2421,11 @@ int scall_ul_pipeline(uldat_t* sl, const enzyme *fn)
 int print_ul_rs(all_ul_t *U_INF)
 {
 	uint32_t i;
-	ul_vec_t *p = NULL;
 	UC_Read ur;
 	init_UC_Read(&ur);
 	for (i = 0; i < U_INF->n; i++) {
-		p = &(U_INF->a[i]);
 		retrieve_ul_t(&ur, NULL, U_INF, i, 0, 0, -1);
-		fprintf(stderr, ">%s\n", p->n_n);
+		fprintf(stderr, ">%s\n", U_INF->nid.a[i].a);
 		fprintf(stderr, "%.*s\n", (int)ur.length, ur.seq);		
 	}
 
@@ -2864,7 +2906,7 @@ void ul_resolve(ma_ug_t *ug, const asg_t *rg, const ug_opt_t *uopt, int hap_n)
 {
 	fprintf(stderr, "[M::%s::] ==> UL\n", __func__);
     mg_idxopt_t opt;
-    init_mg_opt(&opt, 0, 19, 10, hap_n);
+    init_mg_opt(&opt, 0, 19, 10, hap_n, 0, 0, 0.05);
 	int exist = (asm_opt.load_index_from_disk? uidx_load(&ha_flt_tab, &ha_idx, asm_opt.output_file_name) : 0);
     if(exist == 0) uidx_build(ug, &opt);
 	if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name);
@@ -2884,7 +2926,7 @@ int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const enzyme *fn, void *ha
 	sl.uopt = uopt;
 	scall_ul_pipeline(&sl, fn);
 	// print_ul_rs(&UL_INF);
-	debug_retrieve_rc_sub(&UL_INF, &R_INF, &(ug->u), 100);
+	// debug_retrieve_rc_sub(&UL_INF, &R_INF, &(ug->u), 100);
 	// if(!load_ul_hits(&sl.hits, &sl.nn, asm_opt.output_file_name)) {
     // 	scall_ul_pipeline(&sl, fn);
 	// 	write_ul_hits(&sl.hits, &sl.nn, asm_opt.output_file_name);
@@ -2952,14 +2994,15 @@ void ul_load(const ug_opt_t *uopt)
 	fprintf(stderr, "[M::%s::] ==> UL\n", __func__);
 	mg_idxopt_t opt;
 	ma_ug_t *ug = dedup_HiFis(uopt->sources, uopt->min_ovlp, uopt->max_hang, uopt->gap_fuzz);
-	int cutoff = asm_opt.hom_cov * asm_opt.high_factor;
-	init_aux_table();
-	init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff);
-	/**
+	int cutoff;
+	init_aux_table(); ha_opt_update_cov(&asm_opt, asm_opt.hom_cov);
+	cutoff = asm_opt.max_n_chain;
+	init_mg_opt(&opt, !(asm_opt.flag&HA_F_NO_HPC), 19, 10, cutoff, asm_opt.max_n_chain, 0.05, 0.05);
+
 	int exist = (asm_opt.load_index_from_disk? uidx_load(&ha_flt_tab, &ha_idx, asm_opt.output_file_name) : 0);
     if(exist == 0) uidx_l_build(ug, &opt, cutoff);
-	if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name);
-	**/
+	if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name);	
+	
 	ul_v_call(&opt, uopt, asm_opt.ar, ha_flt_tab, ha_idx, ug);
 	ma_ug_destroy(ug); destory_all_ul_t(&UL_INF);
 }