debug uc_block_t

Overlaps.h

@@ -70,7 +70,8 @@ typedef struct {
 typedef struct {
 	uint64_t qns;
 	uint32_t qe, tn, ts, te;
-	uint32_t ml:31, rev:1;
+	// uint32_t ml:31, rev:1;
+	uint32_t cc:30, ml:1, rev:1;
 	uint32_t bl:31, del:1;
 	uint8_t el;
 	uint8_t no_l_indel;
@@ -239,6 +240,7 @@ typedef struct {
 typedef struct {
 	ma_ug_t *ug;
 	ucov_t *cc;
+	ucov_t *cr;
 	ul_contain *ct;
 	// cvert_t *nug;
 	// kv_ul_ov_t *ov;

Process_Read.cpp

@@ -996,7 +996,7 @@ void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str,
 		}
 		
 
-		p->bb.n = p->N_site.n = p->r_base.n = 0;
+		p->bb.n = p->N_site.n = p->r_base.n = 0; p->h = 0;
 		p->rlen = str_l; 
 		
 		if(o == NULL || on == 0) on = 0;
@@ -1020,7 +1020,7 @@ void append_ul_t(all_ul_t *x, uint64_t *rid, char* id, int64_t id_l, char* str,
 
 			///push ovlp bases
 			kv_pushp(uc_block_t, p->bb, &b);
-			b->hid = z->tn; b->rev = z->rev;
+			b->hid = (z->tn<<1)>>1; b->rev = z->rev;
 			b->qs = z->qs; b->qe = z->qe;
 			b->ts = z->ts; b->te = z->te;
 			
@@ -1327,6 +1327,44 @@ uint64_t retrieve_u_cov_region(const ul_idx_t *ul, uint64_t id, uint8_t strand,
 }
 
 
+uint64_t retrieve_r_cov_region(const ul_idx_t *ul, uint64_t id, uint8_t strand, uint64_t s, uint64_t e, int64_t *pi)
+{
+    uint64_t *a = ul->cr->interval.a + ul->cr->idx[id], cc = 0, o = 0, tk, ts, te, tcc = 0;
+    int64_t a_n = ul->cr->idx[id+1]-ul->cr->idx[id], k = 0, cc_i = pi? *pi:0;
+    if(a_n == 0) return e>=s?e-s:0;
+    if(cc_i + 1 >= a_n || cc_i < 0) cc_i = 0;
+    if(strand) {
+		tk = s;
+		s = ul->ug->u.a[id].len - e;
+		e = ul->ug->u.a[id].len - tk;
+	}
+	// fprintf(stderr,"\nul->ug->u.a[id].len:%u, fs:%lu, fe:%lu\n", ul->ug->u.a[id].len, a[k]>>32, a[k+1]>>32);
+	k = cc_i; tk = s;
+	if(tk < (a[k]>>32)) {
+		for (; k >= 0; k--) {
+            if(tk>=(a[k]>>32) && tk<(a[k+1]>>32)) break;
+        }
+	} else if(tk >= (a[k+1]>>32)) {
+		for (; k + 1 < a_n; k++) {
+            if(tk>=(a[k]>>32) && tk<(a[k+1]>>32)) break;
+        }
+	}
+	if(pi) *pi = k;
+
+	
+
+	for (; k + 1 < a_n; k++) {
+		ts = a[k]>>32; te = a[k+1]>>32; cc = (uint32_t)a[k];
+		o = (MIN(e, te) > MAX(s, ts))?(MIN(e, te)-MAX(s, ts)):0;
+		tcc += o*cc;
+		// fprintf(stderr, ">>k:%ld, s:%lu, e:%lu, ts:%lu, te:%lu, o:%lu, cc:%lu\n", k, s, e, ts, te, o, cc);
+		if(e>=(a[k]>>32) && e<(a[k+1]>>32)) break;
+	}
+    
+    
+    return tcc + (e>=s?e-s:0);
+}
+
 uint32_t produce_u_cov(ul_idx_t *ul, uint64_t id, uint8_t strand, uint64_t pos, ma_hit_t_alloc* src, int64_t min_ovlp, int64_t max_hang, int64_t gap_fuzz,
 uint8_t *sset, kvec_t_u64_warp *buf)
 {

Process_Read.h

@@ -159,8 +159,9 @@ typedef struct
 
 typedef struct
 {
-    uint32_t hid:31, rev:1;
+    uint32_t hid;
     uint32_t qs, qe, ts, te;
+    uint8_t pchain:6, rev:1, base:1;
 } uc_block_t;
 
 typedef struct
@@ -177,10 +178,12 @@ typedef struct
 typedef struct
 {
     kvec_t(uint8_t) r_base;
-    uint32_t rlen;
+    uint32_t rlen; 
 
     kvec_t(uc_block_t) bb;
     N_t N_site;
+
+    uint8_t h;
 } ul_vec_t;
 
 typedef struct{
@@ -195,7 +198,7 @@ typedef struct
     ul_vec_t *a;
     size_t n, m;
     All_reads *hR;
-    uint32_t mm;
+    // uint32_t mm;
 } all_ul_t;
 
 extern all_ul_t UL_INF;
@@ -227,5 +230,6 @@ void retrieve_u_seq(UC_Read* i_r, char* i_s, ma_utg_t *u, uint8_t strand, int64_
 void debug_retrieve_rc_sub(const ug_opt_t *uopt, all_ul_t *ref, const All_reads *R_INF, ul_idx_t *ul, uint32_t n_step);
 uint32_t retrieve_u_cov(const ul_idx_t *ul, uint64_t id, uint8_t strand, uint64_t pos, uint8_t dir, int64_t *pi);
 uint64_t retrieve_u_cov_region(const ul_idx_t *ul, uint64_t id, uint8_t strand, uint64_t s, uint64_t e, int64_t *pi);
+uint64_t retrieve_r_cov_region(const ul_idx_t *ul, uint64_t id, uint8_t strand, uint64_t s, uint64_t e, int64_t *pi);
 
 #endif

gfa_ut.cpp

@@ -238,16 +238,16 @@ static void update_sg_uo_t(void *data, long i, int tid)
 	sset_aux *sl = (sset_aux *)data;
 	ma_hit_t_alloc *src = sl->src; asg_t *g = sl->g;
 	asg_arc_t *e = &(g->arc[i]); uint32_t k, qn, tn;
-
+    ma_hit_t_alloc *x = &(src[e->ul>>33]);
-	for (k = 0; k < src[i].length; k++) {
+	for (k = 0; k < x->length; k++) {
-		qn = Get_qn(src[i].buffer[k]);
+		qn = Get_qn(x->buffer[k]);
-		tn = Get_tn(src[i].buffer[k]);
+		tn = Get_tn(x->buffer[k]);
 		if(qn == (e->ul>>33) && tn == (e->v>>1)) {
-			e->ou = (src[i].buffer[k].bl&OU_MASK);
+			e->ou = (x->buffer[k].bl&OU_MASK);
 			break;
 		}
 	}
-	assert(k < src[i].length);
+	assert(k < x->length);
 }
 
 void update_sg_uo(asg_t *g, ma_hit_t_alloc *src)