ONT EC

Assembly.cpp

@@ -904,6 +904,29 @@ void Output_corrected_reads()
     fclose(output_file);
 }
 
+void Output_corrected_fastq()
+{
+    long long i; uint64_t k;
+    UC_Read g_read; asg8_v dv;
+    init_UC_Read(&g_read); kv_init(dv);
+    char* gfa_name = (char*)malloc(strlen(asm_opt.output_file_name)+35);
+    sprintf(gfa_name, "%s.ec.fq", asm_opt.output_file_name);
+    FILE* fp = fopen(gfa_name, "w");
+    free(gfa_name);
+
+    for (i = 0; i < (long long)R_INF.total_reads; i++) {
+        recover_UC_Read(&g_read, &R_INF, i);
+        fprintf(fp, "@%.*s\n", (int32_t)Get_NAME_LENGTH(R_INF, i), Get_NAME(R_INF, i));
+        fprintf(fp, "%.*s\n", (int32_t)g_read.length, g_read.seq);        
+        fprintf(fp, "+\n");
+        retrive_bqual(&dv, NULL, i, -1, -1, 0, sc_bn);
+        for (k = 0; k < dv.n; k++) fprintf(fp, "%c", (char)(sc_tb[dv.a[k]] + 33 - 1));
+        fprintf(fp, "\n");
+    }
+    destory_UC_Read(&g_read); kv_destroy(dv);
+    fclose(fp);
+}
+
 void debug_print_pob_regions()
 {
     uint64_t i, total = 0;
@@ -992,9 +1015,12 @@ void ha_ec(int64_t round, int num_pround, int des_idx, uint64_t *tot_b, uint64_t
 
     if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
 
+    // Output_corrected_fastq();
+
+
     cal_ec_r(asm_opt.thread_num, round, num_pround, R_INF.total_reads, (round == (asm_opt.number_of_round-1))?1:0, tot_b, tot_e);
 
-    // exit(1);
+    // exit(1);    
 
     // if (r_out) write_pt_index(ha_flt_tab, ha_idx, &R_INF, &asm_opt, asm_opt.output_file_name);
     if(des_idx) {
@@ -2014,7 +2040,10 @@ int ha_assemble(void)
 			ha_extract_print_list(&R_INF, asm_opt.extract_iter, asm_opt.extract_list);
 			exit(0);
 		}
-		if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads();
+		if (asm_opt.flag & HA_F_WRITE_EC) {
+            if(asm_opt.is_sc) Output_corrected_fastq();
+            else Output_corrected_reads();
+        }
 		if (asm_opt.flag & HA_F_WRITE_PAF) Output_PAF();
         if (asm_opt.het_cov == -1024) hap_recalculate_peaks(asm_opt.output_file_name), ovlp_loaded = 2;
 	}
@@ -2042,7 +2071,10 @@ int ha_assemble(void)
 			// 		asm_opt.num_bases, asm_opt.num_corrected_bases, asm_opt.num_recorrected_bases);
 			// fprintf(stderr, "[M::%s] size of buffer: %.3fGB\n", __func__, asm_opt.mem_buf / 1073741824.0);
 		}
-		if (asm_opt.flag & HA_F_WRITE_EC) Output_corrected_reads();
+		if (asm_opt.flag & HA_F_WRITE_EC) {
+            if(asm_opt.is_sc) Output_corrected_fastq();
+            else Output_corrected_reads();
+        }
 		// overlap between corrected reads
 		ha_opt_reset_to_round(&asm_opt, asm_opt.number_of_round);
 		// ha_overlap_final();

CommandLines.cpp

@@ -72,6 +72,7 @@ static ko_longopt_t long_options[] = {
     { "telo-s",     ko_required_argument, 357},
     { "ctg-n",     ko_required_argument, 358},
     { "ont",       ko_no_argument, 359},
+    { "sc-n",       ko_no_argument, 360},
     // { "path-round",     ko_required_argument, 348},
 	{ 0, 0, 0 }
 };
@@ -340,7 +341,8 @@ void init_opt(hifiasm_opt_t* asm_opt)
     asm_opt->telo_mic_sc = 500;
 
     asm_opt->is_ont = 0;
-}
+    asm_opt->is_sc = 0;
+}   
 
 void destory_enzyme(enzyme* f)
 {
@@ -912,7 +914,9 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
         else if (c == 357) asm_opt->telo_mic_sc = atol(opt.arg);
         else if (c == 358) asm_opt->max_contig_tip = atol(opt.arg);
         else if (c == 359) {
-            asm_opt->is_ont = 1; asm_opt->max_ov_diff_ec = 0.07;
+            asm_opt->is_ont = 1; asm_opt->max_ov_diff_ec = 0.07; ///asm_opt->mz_win = 37; asm_opt->k_mer_length = 37;
+        } else if (c == 360) {
+            asm_opt->is_sc = 1; 
         } else if (c == 'l') {   ///0: disable purge_dup; 1: purge containment; 2: purge overlap
             asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg);
         }

CommandLines.h

@@ -5,7 +5,7 @@
 #include <pthread.h>
 #include <stdint.h>
 
-#define HA_VERSION "0.20.0-r656"
+#define HA_VERSION "0.21.0-r666"
 
 #define VERBOSE 0
 
@@ -162,6 +162,7 @@ typedef struct {
     int64_t telo_mic_sc;
     
     uint64_t is_ont;
+    uint64_t is_sc;
 } hifiasm_opt_t;
 
 extern hifiasm_opt_t asm_opt;

Correct.h

@@ -144,7 +144,10 @@ typedef struct
     char misBase;
 }haplotype_evdience;
 
-
+#define hh_tp(z) (((z).type&1))
+#define hh_hp(z) ((((z).type>>1)&1))
+#define hh_bq(z) ((((z).type>>2))&sc_bm)
+#define hh_wq(z) ((((z).type>>(sc_bn+2)))&sc_bm)
 
 typedef struct
 {
@@ -1391,10 +1394,11 @@ bit_extz_t *exz, double e_rate, int64_t qs);
 void gen_hc_r_alin(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf);
 void gen_hc_r_alin_nec(overlap_region_alloc* ol, Candidates_list *cl, All_reads *rref, UC_Read* qu, UC_Read* tu, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf);
 uint64_t gen_hc_r_alin_re(overlap_region* z, Candidates_list *cl, char* qstr, uint64_t ql, char* tstr, uint64_t tl, bit_extz_t *exz, overlap_region *aux_o, double e_rate, int64_t wl, int64_t rid, int64_t khit, int64_t move_gap, asg16_v* buf);
-void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, int64_t wl, int64_t ql, uint8_t occ_thres/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs);
+void rphase_hc(overlap_region_alloc* ol, All_reads *rref, haplotype_evdience_alloc* hp, UC_Read* qu, UC_Read* tu, kv_ul_ov_t *c_idx, asg64_v* idx, asg64_v* buf, int64_t bd, int64_t wl, int64_t ql, uint8_t occ_thres/**, uint8_t is_dbg**/, uint64_t rid, uint64_t hpc_len, uint64_t std_bs, Chain_Data *dp, asg8_v *q8, asg8_v *t8);
 void set_exact_exz(bit_extz_t *exz, int64_t qs, int64_t qe, int64_t ts, int64_t te);
 void push_alnw(overlap_region *aux_o, bit_extz_t *exz);
 void cal_exz_global(char *pstr, int32_t pn, char *tstr, int32_t tn, int32_t thre, bit_extz_t *ez);
+void get_wqual(uint64_t zid, uint64_t zpos, uint64_t zrev, asg8_v *v, uint8_t *va, uint64_t scw, uint64_t *tqual, uint64_t *wqual);
 
 
 #define ovlp_id(x) ((x).tn)

Hash_Table.h

@@ -10,6 +10,7 @@
 #define WINDOW_BOUNDARY 375
 #define WINDOW_HC 775
 #define WINDOW_OHC 475
+// #define WINDOW_OHC 375 ///ONT high error
 #define WINDOW_HC_FAST 512
 ///for one side, the first or last WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY bases should not be corrected
 #define WINDOW_UNCORRECT_SINGLE_SIDE_BOUNDARY 25

Process_Read.cpp

@@ -51,6 +51,7 @@ void destory_All_reads(All_reads* r)
 		if (r->read_sperate[i]) free(r->read_sperate[i]);
 		if (r->paf && r->paf[i].buffer) free(r->paf[i].buffer);
 		if (r->reverse_paf && r->reverse_paf[i].buffer) free(r->reverse_paf[i].buffer);
+		if(r->rsc && r->rsc[i]) free(r->rsc[i]);
 		///if (r->pb_regions) kv_destroy(r->pb_regions[i].a);
 	}
 	free(r->paf);
@@ -61,6 +62,7 @@ void destory_All_reads(All_reads* r)
 	free(r->name_index);
 	free(r->read_length);
 	free(r->trio_flag);
+	free(r->rsc);
 	///if (r->pb_regions) free(r->pb_regions);
 }
 
@@ -108,6 +110,14 @@ void write_All_reads(All_reads* r, char* read_file_name)
 	fwrite(&(asm_opt.hom_cov), sizeof(asm_opt.hom_cov), 1, fp);
 	fwrite(&(asm_opt.het_cov), sizeof(asm_opt.het_cov), 1, fp);
 
+	uint64_t mm = 1;
+	if(asm_opt.is_sc) {
+		fwrite(&mm, sizeof(mm), 1, fp);
+		for (i = 0; i < r->total_reads; i++) {
+			fwrite(r->rsc[i], sizeof(uint8_t), ((r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0)), fp);
+		}
+	}
+
     free(index_name);    
 	fflush(fp);
     fclose(fp);
@@ -201,6 +211,19 @@ int load_All_reads(All_reads* r, char* read_file_name)
 	}
 	///r->pb_regions = NULL;
 
+	uint64_t mm = 0;
+	if (!feof(fp)) {
+		if((fread(&mm, sizeof(mm), 1, fp)) && (mm == 1)) {
+			MALLOC(r->rsc, r->total_reads);
+			for (i = 0; i < r->total_reads; i++) {
+				MALLOC(r->rsc[i], (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0));
+				f_flag += fread(r->rsc[i], sizeof(uint8_t), (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0), fp);
+			}
+		}
+	}
+
+
+
     free(index_name);    
     fclose(fp);
     fprintf(stderr, "Reads has been loaded.\n");
@@ -412,10 +435,13 @@ void malloc_All_reads(All_reads* r)
 	memcpy(r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads);
 
 	r->read_sperate = (uint8_t**)malloc(sizeof(uint8_t*)*r->total_reads);
+	if(asm_opt.is_sc) MALLOC(r->rsc, r->total_reads);
+
 	long long i = 0;
 	for (i = 0; i < (long long)r->total_reads; i++)
 	{
 		r->read_sperate[i] = (uint8_t*)malloc(sizeof(uint8_t)*(r->read_length[i]/4+1));
+		if(r->rsc) MALLOC(r->rsc[i], (r->read_length[i]/sc_bn) + ((r->read_length[i]%sc_bn)?1:0));
 	}
 
 	r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
@@ -823,6 +849,131 @@ void ha_compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_sit
 	}
 }
 
+void convert_qual(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitu, uint64_t rev, uint64_t sc_off)
+{
+	uint64_t i = 0; uint8_t c = 0, sc;
+	// fprintf(stderr, "\n[M::%s]\n", __func__);
+	for (i = 0; i < src_l; i++) {
+		for (c = 0, sc = ((uint8_t)src[i]) - sc_off; (c < bitu) && (sc_tb[c] < sc); c++);
+		if(c >= bitu) c = bitu - 1;
+		dest[(rev?(src_l-i-1):(i))] = c;
+		// fprintf(stderr, "%u->%u\n", sc, c);
+	}
+}
+
+void ha_compress_qual_bit(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitn)
+{
+
+	uint64_t i = 0, k, bit_r = 8/bitn, dest_i = 0;
+	uint8_t tmp = 0, c = 0;
+
+	for (i = 0; i + bit_r <= src_l;) {
+		for (k = tmp = 0; k < bit_r; k++) {
+			c = ((uint8_t)src[i]);
+			tmp <<= bitn; tmp |= c; i++;
+		}
+		dest[dest_i++] = tmp;
+	}
+
+	if(i < src_l) {
+		for (k = tmp = 0; i < src_l; k++) {
+			c = ((uint8_t)src[i]);
+			tmp <<= bitn; tmp |= c; i++;
+		}
+		 
+		dest[dest_i++] = (tmp<<(8-(bitn*k)));
+	}
+}
+
+void ha_compress_qual(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitn, uint64_t sc_off)
+{
+
+	uint64_t i = 0, k, bit_r = 8/bitn, dest_i = 0, bitu = (1<<bitn);
+	uint8_t tmp = 0, c = 0, sc;
+
+	for (i = 0; i + bit_r <= src_l;) {
+		for (k = tmp = 0; k < bit_r; k++) {
+			for (c = 0, sc = ((uint8_t)src[i]) - sc_off; (c < bitu) && (sc_tb[c] < sc); c++);
+			if(c >= bitu) c = bitu - 1;
+			tmp <<= bitn; tmp |= c; i++;
+		}
+		dest[dest_i++] = tmp;
+	}
+
+	if(i < src_l) {
+		for (k = tmp = 0; i < src_l; k++) {
+			for (c = 0, sc = ((uint8_t)src[i]) - sc_off; (c < bitu) && (sc_tb[c] < sc); c++);
+			if(c >= bitu) c = bitu - 1;
+			tmp <<= bitn; tmp |= c; i++;
+		}
+		 
+		dest[dest_i++] = (tmp<<(8-(bitn*k)));
+	}
+}
+
+///[s, e)
+int64_t retrive_bqual(asg8_v *dv, uint8_t *ds, uint64_t id, int64_t s, int64_t e, uint8_t rev, int64_t bitn)
+{
+	int64_t rl = Get_READ_LENGTH(R_INF, id), l; 
+	if(s < 0) s = 0; if(e < 0) e = rl;
+	if(s >= e || e > rl) return -1;
+
+	uint8_t *da = NULL, *src = Get_QUAL(R_INF, id), mm = (((uint8_t)1)<<bitn)-1, mlf = 8 - bitn, mrf;
+	int64_t bitr = 8/bitn, dk, sk, swk;
+	l = e - s;
+	if(dv) {
+		kv_resize(uint8_t, *dv, ((uint64_t)l)); da = dv->a;
+	} else {
+		da = ds;
+	}
+
+	if(!rev) {
+		dk = 0; sk = s;
+
+		mrf = ((s%bitr)*bitn);
+		// if(s == 21519 && e == 22332) {
+		// 	fprintf(stderr, "+[M::%s] id::%lu, in::[%ld, %ld), rev::%u, bitn::%ld, bitr::%ld, mrf::%u\n", __func__, id, s, e, rev, bitn, bitr, mrf);
+		// }
+		if(mrf) {
+			for (swk = sk/bitr; mrf < 8 && sk < e; mrf += bitn, sk++) da[dk++] = ((src[swk]<<mrf)>>mlf)&mm;
+		}
+		
+		for (swk = sk/bitr; (sk + bitr) <= e; sk += bitr, swk++) {
+			for (mrf = 0; mrf < 8; mrf += bitn) da[dk++] = ((src[swk]<<mrf)>>mlf)&mm;
+		}
+
+		if(sk < e) {
+			for (mrf = 0; sk < e; mrf += bitn, sk++) da[dk++] = ((src[swk]<<mrf)>>mlf)&mm;
+		}
+		// if(dk != l) {
+		// 	fprintf(stderr, "+[M::%s] id::%lu, in::[%ld, %ld), rev::%u, bitn::%ld, bitr::%ld\n", __func__, id, s, e, rev, bitn, bitr);
+		// }
+		assert(dk == l);
+	} else {
+		sk = s; s = e; e = sk;
+		s = rl - s; e = rl - e;
+		dk = l; sk = s; 
+
+		mrf = ((s%bitr)*bitn);
+		if(mrf) {
+			for (swk = sk/bitr; mrf < 8 && sk < e; mrf += bitn, sk++) da[--dk] = ((src[swk]<<mrf)>>mlf)&mm;
+		}
+		
+		for (swk = sk/bitr; (sk + bitr) <= e; sk += bitr, swk++) {
+			for (mrf = 0; mrf < 8; mrf += bitn) da[--dk] = ((src[swk]<<mrf)>>mlf)&mm;
+		}
+
+		if(sk < e) {
+			for (mrf = 0; sk < e; mrf += bitn, sk++) da[--dk] = ((src[swk]<<mrf)>>mlf)&mm;
+		}
+		
+		assert(dk == 0);
+	}
+	dv->n = l;
+
+	return l;
+}
+
 void reverse_complement(char* pattern, uint64_t length)
 {
 	uint64_t i = 0;
@@ -845,6 +996,24 @@ void reverse_complement(char* pattern, uint64_t length)
 	}
 }
 
+void print_fastq(FILE *fp, char *id, char *bs, char *qual, uint64_t bitu, uint64_t sc_off)
+{
+	uint64_t i = 0, ql = strlen(qual); uint8_t c = 0, sc;
+	
+	if(fp) fprintf(fp, "@%s\n%s\n+\n", id, bs);
+	else fprintf(stdout, "@%s\n%s\n+\n", id, bs);
+
+	for (i = 0; i < ql; i++) {
+		for (c = 0, sc = ((uint8_t)qual[i]) - sc_off; (c < bitu) && (sc_tb[c] < sc); c++);
+		if(c >= bitu) c = bitu - 1;
+		if(fp) fprintf(fp, "%u", c);
+		else fprintf(stdout, "%u", c);
+	}
+
+	if(fp) fprintf(fp, "\n");
+	else fprintf(stdout, "\n");
+}
+
 
 void init_Debug_reads(Debug_reads* x, const char* file)
 {

Process_Read.h

@@ -22,6 +22,7 @@
 #define Get_NAME_LENGTH(R_INF, ID) ((R_INF).name_index[(ID)+1] - (R_INF).name_index[(ID)])
 ///#define Get_READ(R_INF, ID) R_INF.read + (R_INF.index[ID]>>2) + ID
 #define Get_READ(R_INF, ID) (R_INF).read_sperate[(ID)]
+#define Get_QUAL(R_INF, ID) (R_INF).rsc[(ID)]
 #define Get_NAME(R_INF, ID) ((R_INF).name + (R_INF).name_index[(ID)])
 #define CHECK_BY_NAME(R_INF, NAME, ID) (Get_NAME_LENGTH((R_INF),(ID))==strlen((NAME)) && \
                                         memcmp((NAME), Get_NAME((R_INF), (ID)), Get_NAME_LENGTH((R_INF),(ID))) == 0)
@@ -38,6 +39,8 @@ extern char rc_Table[6];
 
 void init_aux_table();
 
+typedef struct { size_t n, m; uint8_t *a; } asg8_v;
+
 typedef struct
 {
     uint64_t x_id;
@@ -119,6 +122,7 @@ typedef struct
 	uint64_t* read_length;
 	uint64_t* read_size;
 	uint8_t* trio_flag;
+    uint8_t** rsc;
 
 	///seq start pos in uint8_t* read
 	///do not need it
@@ -223,6 +227,7 @@ void init_All_reads(All_reads* r);
 void malloc_All_reads(All_reads* r);
 void ha_insert_read_len(All_reads *r, int read_len, int name_len);
 void ha_compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ);
+void ha_compress_qual(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitn, uint64_t sc_off);
 void init_UC_Read(UC_Read* r);
 void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID);
 void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID);
@@ -254,4 +259,17 @@ scaf_res_t *init_scaf_res_t(uint32_t n);
 void destroy_scaf_res_t(scaf_res_t *p);
 void read_ma(ma_hit_t* x, FILE* fp);
 
+const uint64_t sc_tb[8] = { 
+	10, 20, 30, 40, 50, 60, 70, 80
+};
+
+#define sc_bn 2
+#define sc_bm ((((uint64_t)1)<<sc_bn)-1)
+#define sc_wn 5
+
+void convert_qual(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitu, uint64_t rev, uint64_t sc_off);
+int64_t retrive_bqual(asg8_v *dv, uint8_t *ds, uint64_t id, int64_t s, int64_t e, uint8_t rev, int64_t bitn);
+void print_fastq(FILE *fp, char *id, char *bs, char *qual, uint64_t bitu, uint64_t sc_off);
+void ha_compress_qual_bit(uint8_t* dest, char* src, uint64_t src_l, uint64_t bitn);
+
 #endif

ecovlp.cpp

@@ -27,7 +27,7 @@ cc_v scc = {0, 0, NULL, NULL};
 cc_v scb = {0, 0, NULL, NULL};
 cc_v sca = {0, 0, NULL, NULL};
 
-typedef struct {size_t n, m; char *a; UC_Read z;} sl_v;
+typedef struct {size_t n, m; char *a; UC_Read z; asg8_v q;} sl_v;
 
 
 void h_ec_lchain(ha_abuf_t *ab, uint32_t rid, char* rs, uint64_t rl, uint64_t mz_w, uint64_t mz_k, All_reads *rref, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, 
@@ -77,6 +77,8 @@ ec_ovec_buf_t* gen_ec_ovec_buf_t(uint32_t n)
 	    kv_init(z->v64);
         kv_init(z->v32);
         kv_init(z->v16);
+        kv_init(z->v8q);
+        kv_init(z->v8t);
         init_bit_extz_t(&(z->exz), 31);
 
         z->ab = ha_abuf_init();
@@ -117,6 +119,8 @@ void destroy_ec_ovec_buf_t(ec_ovec_buf_t *p)
 	    kv_destroy(z->v64);
         kv_destroy(z->v32);
         kv_destroy(z->v16);
+        kv_destroy(z->v8q);
+        kv_destroy(z->v8t);
         destroy_bit_extz_t(&(z->exz));
 
         ha_abuf_destroy(z->ab);
@@ -148,6 +152,8 @@ inline void refresh_ec_ovec_buf_t0(ec_ovec_buf_t0 *z, uint64_t n)
         kv_destroy(z->v64); kv_init(z->v64);
         kv_destroy(z->v32); kv_init(z->v32);
         kv_destroy(z->v16); kv_init(z->v16);
+        kv_destroy(z->v8q); kv_init(z->v8q);
+        kv_destroy(z->v8t); kv_init(z->v8t);
 
         destroy_bit_extz_t(&(z->exz)); init_bit_extz_t(&(z->exz), 31);
 
@@ -2930,6 +2936,56 @@ void dedup_chains(overlap_region_alloc* ol)
     }
 }
 
+void debug_retrive_bqual(asg8_v *vq, asg8_v *vt, uint64_t id, uint64_t rn)
+{
+    uint64_t k, n, z[2], s, e, rev;
+    retrive_bqual(vq, NULL, id, -1, -1, 0, sc_bn); n = vq->n;
+    retrive_bqual(vt, NULL, id, -1, -1, 1, sc_bn);
+    assert(vq->n == vt->n);
+    for (k = 0; k < vq->n && vq->a[k] == vt->a[vt->n - k - 1]; k++);
+    // if((k == vq->n)) {
+    //     fprintf(stderr, "[M::%s] id::%lu, k::%lu, n::%lu\n", __func__, id, k, ((uint64_t)vq->n));
+    // }
+    assert(k == vq->n);
+
+    // if(id == 0) {
+    //     s = 21519; e = 22332; rev = 0;
+    //     retrive_bqual(vt, NULL, id, s, e, rev, sc_bn);
+    //     if(memcmp(vq->a + s, vt->a, e - s)) {
+    //         fprintf(stderr, "+[M::%s] id::%lu, t::[%lu, %lu), rev::%lu\n", __func__, id, s, e, rev);
+    //         exit(1);
+    //     }
+
+    // }
+    // return;
+
+    for (k = 0, rev = 0; k < rn; k++) {
+        z[0] = rand()%(n + 1); 
+        z[1] = rand()%(n + 1); 
+        if(z[0] == z[1]) continue;
+        s = MIN(z[0], z[1]); e = MAX(z[0], z[1]);
+        retrive_bqual(vt, NULL, id, s, e, rev, sc_bn);
+        if(memcmp(vq->a + s, vt->a, e - s)) {
+            fprintf(stderr, "[M::%s] id::%lu, t::[%lu, %lu), rev::%lu\n", __func__, id, s, e, rev);
+            exit(1);
+        }
+    }
+    
+
+    retrive_bqual(vq, NULL, id, -1, -1, 1, sc_bn);
+    for (k = 0, rev = 1; k < rn; k++) {
+        z[0] = rand()%(n + 1); 
+        z[1] = rand()%(n + 1); 
+        if(z[0] == z[1]) continue;
+        s = MIN(z[0], z[1]); e = MAX(z[0], z[1]);
+        retrive_bqual(vt, NULL, id, s, e, rev, sc_bn);
+        if(memcmp(vq->a + s, vt->a, e - s)) {
+            fprintf(stderr, "[M::%s] id::%lu, t::[%lu, %lu), rev::%lu\n", __func__, id, s, e, rev);
+            exit(1);
+        }
+    }
+}
+
 static void worker_hap_ec(void *data, long i, int tid)
 {
 	ec_ovec_buf_t0 *b = &(((ec_ovec_buf_t*)data)->a[tid]);
@@ -2948,7 +3004,9 @@ static void worker_hap_ec(void *data, long i, int tid)
     **/
     // if(i < 1100000 || i > 1400000) return;
     // if(i % 100000 == 0) fprintf(stderr, "-a-[M::%s-beg] rid->%ld\n", __func__, i);
-    // if (memcmp("4da034b0-a94d-4576-8481-c0d9a9f96d40", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
+    // if (memcmp("c42804f3-0e13-43a0-8a71-b91b40accf9a", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
+    // if (memcmp("b2e68ecf-381a-439c-b676-c1e6831d6acf", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
+    // if (memcmp("0aec8c4f-c849-4c31-85ba-4ffb297eeb28", Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0) {
     //     fprintf(stderr, "-a-[M::%s-beg] rid->%ld\n", __func__, i);
     // } else {
     //     return;
@@ -2961,9 +3019,11 @@ static void worker_hap_ec(void *data, long i, int tid)
     // if(i != 2243244) return;
     // if(i != 19350) return;
 
+    // debug_retrive_bqual(D, &b->v8t, i, 256); return;
+
     recover_UC_Read(&b->self_read, &R_INF, i); qlen = b->self_read.length; 
 
-    h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, 0.02, asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32);
+    h_ec_lchain(b->ab, i, b->self_read.seq, b->self_read.length, asm_opt.mz_win, asm_opt.k_mer_length, &R_INF, &b->olist, &b->clist, /**((asm_opt.is_ont)?(0.05):(0.02))**/0.02, asm_opt.max_n_chain, 1, NULL, NULL, &(b->sp), &high_occ, &low_occ, 1, 1, 3, 0.7, 2, 32);///ONT high error
 
     // b->num_read_base += b->olist.length;
     b->cnt[0] += b->self_read.length;
@@ -2987,7 +3047,7 @@ static void worker_hap_ec(void *data, long i, int tid)
     // b->num_correct_base += b->olist.length;
 
     copy_asg_arr(buf0, b->sp); 
-    rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont);
+    rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 0**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), ((asm_opt.is_sc)?&(b->v8t):NULL));
     copy_asg_arr(b->sp, buf0);
 
     // stderr_phase_ovlp(&b->olist);
@@ -3495,6 +3555,12 @@ static void worker_hap_post_rev(void *data, long i, int tid)
     }
 
     ha_compress_base(Get_READ(R_INF, i), a, l, &R_INF.N_site[i], nn);
+
+    if(asm_opt.is_sc) {
+        retrive_bqual(&(b->v8q), NULL, i, -1, -1, 0, sc_bn);
+        for (k = 0; k < l; k++) a[l - k - 1] = b->v8q.a[k];
+        ha_compress_qual_bit(Get_QUAL(R_INF, i), a, l, sc_bn);
+    }
 }
 
 static void worker_hap_dc_ec_gen(void *data, long i, int tid)
@@ -5018,7 +5084,7 @@ static void worker_hap_dc_ec0(void *data, long i, int tid)
     b->cnt[0] += b->self_read.length;
 
     copy_asg_arr(buf0, b->sp); 
-    rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 1**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont);
+    rphase_hc(&b->olist, &R_INF, &b->hap, &b->self_read, &b->ovlp_read, &b->pidx, &b->v64, &buf0, 0, WINDOW_MAX_SIZE, b->self_read.length, 1/**, 1**/, i, (asm_opt.is_ont)?HPC_PL:0, asm_opt.is_ont, ((asm_opt.is_ont)?&(b->clist.chainDP):NULL), ((asm_opt.is_sc)?&(b->v8q):NULL), ((asm_opt.is_sc)?&(b->v8t):NULL));
     copy_asg_arr(b->sp, buf0); 
 
     copy_asg_arr(buf0, b->sp);
@@ -5165,7 +5231,7 @@ static void worker_sl_ec(void *data, long i, int tid)
 {
     // if(i != 0) return;
 
-	sl_v *p = &(((sl_v*)data)[tid]);
+	sl_v *p = &(((sl_v*)data)[tid]); uint8_t *oa = NULL; char *na = NULL; uint64_t tqual, wqual;
     uint32_t ci = 0, len, xk, yk, wx[2], wy[2], k, Nn, yn = 0, tot_e; uint16_t c, bq, bt; 
 
 
@@ -5213,11 +5279,13 @@ static void worker_sl_ec(void *data, long i, int tid)
     }
 
     // if(i == 700) fprintf(stderr, "|\n");
+    if(asm_opt.is_sc) retrive_bqual(&(p->q), NULL, i, -1, -1, 0, sc_bn);
 
 
     if (R_INF.read_size[i] < yn) {
 		R_INF.read_size[i] = yn;
 		REALLOC(R_INF.read_sperate[i], R_INF.read_size[i]/4+1);
+        if(asm_opt.is_sc) REALLOC(R_INF.rsc[i], ((R_INF.read_size[i]/sc_bn) + ((R_INF.read_size[i]%sc_bn)?1:0)));
 	}
 	R_INF.read_length[i] = yn;
     // if(Nn > 0) fprintf(stderr, "[M::%s] Nn->%u\n", __func__, Nn);
@@ -5236,7 +5304,25 @@ static void worker_sl_ec(void *data, long i, int tid)
 
 
 	ha_compress_base(Get_READ(R_INF, i), p->a, yn, &R_INF.N_site[i], Nn);
-
+    if(asm_opt.is_sc) {
+        oa = p->q.a; na = p->a; 
+        ci = 0; xk = yk = 0; Nn = 0;
+        while (ci < scc.a[i].n) {
+            wx[0] = xk; wy[0] = yk;
+            ci = pop_trace_bp_f(&scc.a[i], ci, &c, &bq, &bt, &len);
+            if(c != 2) xk += len;
+            if(c != 3) yk += len;
+            wx[1] = xk; wy[1] = yk;
+            if(c == 0 || c == 1) {
+                memcpy(na + wy[0], oa + wx[0], (wx[1]-wx[0])*sizeof((*oa)));
+            } else if(c == 2) {
+                get_wqual(i, wx[0], 0, NULL, oa, sc_wn, &tqual, &wqual);
+                for (k = wy[0]; k < wy[1]; k++) na[k] = wqual;
+            }
+        }
+        assert(yk == yn);
+        ha_compress_qual_bit(Get_QUAL(R_INF, i), na, yn, sc_bn);
+    }
 }
 
 uint64_t cal_ec_multiple(ec_ovec_buf_t *b, uint64_t n_thre, uint64_t n_a, uint64_t *r_base)
@@ -5496,7 +5582,7 @@ void sl_ec_r(uint64_t n_thre, uint64_t n_a)
     sl_v *b = NULL; uint64_t k; MALLOC(b, n_thre);
     for (k = 0; k < n_thre; k++) {
         b[k].a = NULL; b[k].n = b[k].m = 0;
-        init_UC_Read(&b[k].z);
+        init_UC_Read(&b[k].z); kv_init(b[k].q);
     }
 
     kt_for(n_thre, worker_sl_ec, b, n_a);///debug_for_fix

ecovlp.h

@@ -42,6 +42,7 @@ typedef struct {
 	asg64_v v64;
 	asg32_v v32;
 	asg16_v v16;
+	asg8_v v8q, v8t;
 
     kvec_t_u8_warp k_flag;
 	st_mt_t sp;

htab.cpp

@@ -692,6 +692,7 @@ static inline void sf##_pt_insert_buf(sf##_ch_buf_t *buf, int p, const HType *y)
 static void *sf##_worker_count(void *data, int step, void *in) /** callback for kt_pipeline()**/\
 {\
 	pl_data_t *p = (pl_data_t*)data;\
+	/**uint8_t src_a[1000000], des_a[1000000];**/\
 	if (step == 0) { /** step 1: read a block of sequences**/\
 		int ret;\
 		sf##_st_data_t *s;\
@@ -762,6 +763,16 @@ static void *sf##_worker_count(void *data, int step, void *in) /** callback for
 								++n_N;\
 						ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s+p->opt->adaLen, l, &p->rs_out->N_site[p->n_seq], n_N);\
 						memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);\
+						if(p->rs_out->rsc) {\
+							ha_compress_qual(Get_QUAL(*p->rs_out, p->n_seq), p->ks->qual.s+p->opt->adaLen, l, sc_bn, 33);\
+							/**print_fastq(NULL, p->ks->name.s, p->ks->seq.s, p->ks->qual.s, (1<<sc_bn), 33);**/\
+							/**if(l <= 1000000) {\
+								convert_qual(src_a, p->ks->qual.s+p->opt->adaLen, l, (1<<sc_bn), 0, 33);\
+								retrive_bqual(NULL, des_a, p->n_seq, -1, -1, 0, sc_bn);\
+								if(memcmp(src_a, des_a, l)!=0) fprintf(stderr, "ERROR: incorrect qual values\n");\
+								else fprintf(stderr, "Correct: correct qual values\n");\
+							}**/\
+						}\
 					}\
 				}\
 				if (s->n_seq == s->m_seq) {\