diff --git a/Assembly.cpp b/Assembly.cpp index 39cc6a1..d012c5b 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -20,6 +20,7 @@ void ha_sort_list_by_anchor(overlap_region_alloc *overlap_list); All_reads R_INF; Debug_reads R_INF_FLAG; +all_ul_t UL_INF; void get_corrected_read_from_cigar(Cigar_record* cigar, char* pre_read, int pre_length, char* new_read, int* new_length) { @@ -417,31 +418,6 @@ long long push_final_overlaps_increment(ma_hit_t_alloc* paf, ma_hit_t_alloc* rev return available_overlaps; } -typedef struct { - int is_final, save_ov; - // chaining and overlapping related buffers - UC_Read self_read, ovlp_read; - Candidates_list clist; - overlap_region_alloc olist; - overlap_region_alloc olist_hp; - ha_abuf_t *ab; - ha_abufl_t *abl; - // error correction related buffers - int64_t num_read_base, num_correct_base, num_recorrect_base; - Cigar_record cigar1; - Graph POA_Graph; - Graph DAGCon; - Correct_dumy correct; - haplotype_evdience_alloc hap; - Round2_alignment round2; - kvec_t_u32_warp b_buf; - kvec_t_u64_warp r_buf; - kvec_t_u8_warp k_flag; - overlap_region tmp_region; - ma_utg_v *ua; - st_mt_t sp; -} ha_ovec_buf_t; - ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov, int is_ug) { ha_ovec_buf_t *b; diff --git a/Assembly.h b/Assembly.h index 7e61a28..7f37338 100644 --- a/Assembly.h +++ b/Assembly.h @@ -2,6 +2,9 @@ #define __ASSEMBLY__ #include "CommandLines.h" #include "Overlaps.h" +#include "Process_Read.h" +#include "Hash_Table.h" +#include "Correct.h" #define FORWARD 0 #define REVERSE_COMPLEMENT (0x8000000000000000) @@ -14,7 +17,34 @@ #define RESEED_LEN 2000 #define RESEED_HP_RATE 0.9 +typedef struct { + int is_final, save_ov; + // chaining and overlapping related buffers + UC_Read self_read, ovlp_read; + Candidates_list clist; + overlap_region_alloc olist; + overlap_region_alloc olist_hp; + ha_abuf_t *ab; + ha_abufl_t *abl; + // error correction related buffers + int64_t num_read_base, num_correct_base, num_recorrect_base; + Cigar_record cigar1; + Graph POA_Graph; + Graph DAGCon; + Correct_dumy correct; + haplotype_evdience_alloc hap; + Round2_alignment round2; + kvec_t_u32_warp b_buf; + kvec_t_u64_warp r_buf; + kvec_t_u8_warp k_flag; + overlap_region tmp_region; + ma_utg_v *ua; + st_mt_t sp; +} ha_ovec_buf_t; + int ha_assemble(void); void ug_idx_build(ma_ug_t *ug, int hap_n); +ha_ovec_buf_t *ha_ovec_init(int is_final, int save_ov, int is_ug); +void ha_ovec_destroy(ha_ovec_buf_t *b); #endif diff --git a/Makefile b/Makefile index 8da074f..8371724 100644 --- a/Makefile +++ b/Makefile @@ -76,5 +76,5 @@ hic.o: hic.h rcut.o: rcut.h horder.o: horder.h tovlp.o: tovlp.h -inter.o: inter.h +inter.o: inter.h Process_Read.h kalloc.o: kalloc.h \ No newline at end of file diff --git a/Overlaps.cpp b/Overlaps.cpp index e20a5d0..e5832bf 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -31102,6 +31102,13 @@ char *get_outfile_name(char* output_file_name) return buf; } +void create_ul_info() +{ + ug_opt_t opt; memset(&opt, 0, sizeof(opt)); + ul_load(&opt); + exit(1); +} + void clean_graph( int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, long long n_read, uint64_t* readLen, long long mini_overlap_length, @@ -31110,6 +31117,30 @@ float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name, long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr, ma_sub_t **coverage_cut_ptr, int debug_g) { + // if(asm_opt.ar) create_ul_info(); + + + + + + + + + + + + + + + + + + + + + + + char *o_file = get_outfile_name(output_file_name); ma_sub_t *coverage_cut = *coverage_cut_ptr; asg_t *sg = *sg_ptr; diff --git a/Overlaps.h b/Overlaps.h index d393ae1..2bf7342 100644 --- a/Overlaps.h +++ b/Overlaps.h @@ -53,7 +53,12 @@ #define CUT 11 #define CUT_DIF_HAP 12 - +///query is the read itself +typedef struct { + uint32_t qn, qs, qe; + uint32_t tn, ts, te; + uint8_t sec:6, el:1, rev:1; +} ul_ov_t; ///query is the read itself diff --git a/Process_Read.cpp b/Process_Read.cpp index c53e0dd..cf70f84 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -3,6 +3,8 @@ #include #include #include "Process_Read.h" +#include "htab.h" +#include "Correct.h" uint8_t seq_nt6_table[256] = { 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, @@ -735,4 +737,198 @@ void destory_Debug_reads(Debug_reads* x) free(x->read_name); fclose(x->fp); +} + + +void init_all_ul_t(all_ul_t *x, All_reads *hR) { + memset(x, 0, sizeof(*x)); + x->hR = hR; x->mm = 0x7fffffff; +} +void destory_all_ul_t(all_ul_t *x, All_reads *hR) { + 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); +} + +void ha_encode_base(uint8_t* dest, char* src, uint64_t src_l, N_t *nn) +{ + uint64_t i = 0; + uint64_t dest_i = 0; + uint8_t tmp = 0; + uint8_t c = 0; + + while (i + 4 <= src_l) { + tmp = 0; + + c = seq_nt6_table[(uint8_t)src[i]]; + if (c >= 4) { + c = 0; kv_push(uint32_t, *nn, i); + } + i++; + tmp = tmp | (c<<6); + + c = seq_nt6_table[(uint8_t)src[i]]; + if (c >= 4) { + c = 0; kv_push(uint32_t, *nn, i); + } + i++; + tmp = tmp | (c<<4); + + c = seq_nt6_table[(uint8_t)src[i]]; + if (c >= 4) { + c = 0; kv_push(uint32_t, *nn, i); + } + i++; + tmp = tmp | (c<<2); + + c = seq_nt6_table[(uint8_t)src[i]]; + if (c >= 4) { + c = 0; kv_push(uint32_t, *nn, i); + } + i++; + tmp = tmp | c; + + dest[dest_i] = tmp; + + dest_i++; + } + + //at most 3 bases here + uint64_t shift = 6; + if (i < src_l) { + tmp = 0; + while (i < src_l) { + c = seq_nt6_table[(uint8_t)src[i]]; + if (c >= 4) { + c = 0; kv_push(uint32_t, *nn, i); + } + i++; + tmp = tmp | (c << shift); + shift = shift -2; + } + + dest[dest_i] = tmp; + dest_i++; + } +} + +#define B4L(x) (((x)>>2)+(((x)&3)?1:0)) +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; + ul_ov_t *z = NULL, *zp = NULL; + uc_block_t *b = NULL; + if(rid == NULL) { + kv_pushp(ul_vec_t, *x, &p); + memset(p, 0, sizeof(*p)); + } + else { + p = &(x->a[*rid]); + } + + 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) { + p->bb.n = p->N_site.n = p->r_base.n = 0; + p->rlen = str_l; + + if(o == NULL || on == 0) on = 0; + for (i = end = 0, zp = NULL; i < on; i++) { + z = &(o[i]); + if(!z->el) continue; + if(zp) { + mine = MIN(zp->qe, z->qe); maxs = MAX(zp->qs, z->qs); + ovlp = mine - maxs; + if(zp->qe >= z->qe && zp->qs <= z->qs) continue; + } else { + ovlp = -z->qs; + } + if(ovlp < 0) {///push original bases + kv_pushp(uc_block_t, p->bb, &b); + b->hid = x->mm; b->rev = 0; + b->qs = end; b->qe = b->qs - ovlp; + b->ts = p->r_base.n; b->te = b->ts + B4L(-ovlp); + kv_resize(uint8_t, p->r_base, b->te); p->r_base.n = b->te; + ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site)); + } + + ///push ovlp bases + kv_pushp(uc_block_t, p->bb, &b); + b->hid = z->tn; b->rev = z->rev; + b->qs = z->qs + (ovlp>0?ovlp:0); b->qe = z->qe; + if(z->rev) { + b->ts = z->ts; b->te = z->ts + (b->qe - b->qs); + } else { + b->ts = z->te - (b->qe - b->qs); b->te = z->te; + } + + end = MAX(zp?zp->qe:0, z->qe); + } + + if(end < str_l) {///push original bases + kv_pushp(uc_block_t, p->bb, &b); + b->hid = x->mm; b->rev = 0; + b->qs = end; b->qe = str_l; + b->ts = p->r_base.n; b->te = b->ts + B4L(b->qe - b->qs); + kv_resize(uint8_t, p->r_base, b->te); p->r_base.n = b->te; + ha_encode_base(p->r_base.a+b->ts, str+b->qs, b->qe-b->qs, &(p->N_site)); + } + } +} + + +void retrieve_ul_t(UC_Read* r, all_ul_t *ref, uint64_t ID, uint8_t strand) { + ul_vec_t *p = &(ref->a[ID]); + uc_block_t *b = NULL; + char *a = NULL; + uint8_t *src = NULL; + int64_t k, i, a_n, l_chr, idx; + r->length = p->rlen; r->RID = ID; + + if(r->length + 4 > r->size) { + r->size = r->length + 4; + r->seq = (char*)realloc(r->seq,sizeof(char)*(r->size)); + } + + a = NULL; + if(strand == 0) { + for (k = 0; k < (int64_t)p->bb.n; k++) { + b = &(p->bb.a[k]); a = r->seq + b->qs; a_n = b->qe - b->qs; src = p->r_base.a + b->ts; + if(b->hid == ref->mm){///original bases + i = 0; + while (i < a_n) { + memcpy(a+i, bit_t_seq_table[src[i>>2]], 4); + i += 4; + } + } else {///ovlps + recover_UC_Read_sub_region(a, b->ts, b->te - b->ts, b->rev, ref->hR, b->hid); + } + } + for (k = 0; k < (int64_t)p->N_site.n; k++) r->seq[p->N_site.a[k]] = 'N'; + } else { + for (k = 0; k < (int64_t)p->bb.n; k++) { + b = &(p->bb.a[k]); a = r->seq + r->length - b->qe; a_n = b->qe - b->qs; src = p->r_base.a + b->ts; + if(b->hid == ref->mm){///original bases + l_chr = (b->qe - b->qs)&((uint32_t)3); + i = ((b->qe - b->qs)>>2)+(l_chr!= 0)-1; + idx = 0; + if(l_chr > 0) { + memcpy(a + idx, bit_t_seq_table_rc[src[i]]+4-l_chr, l_chr); + idx = l_chr; i--; + } + while (i >= 0) { + memcpy(a + idx, bit_t_seq_table_rc[src[i]], 4); + i--; idx += 4; + } + } else {///ovlps + recover_UC_Read_sub_region(a, b->ts, b->te - b->ts, b->rev==strand?0:1, ref->hR, b->hid); + } + } + for (k = 0; k < (int64_t)p->N_site.n; k++) r->seq[r->length - p->N_site.a[k] - 1] = 'N'; + } } \ No newline at end of file diff --git a/Process_Read.h b/Process_Read.h index 2ab7b42..49a7986 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -156,6 +156,41 @@ typedef struct pthread_mutex_t OutputMutex; } Debug_reads; + +typedef struct +{ + uint32_t hid:31, rev:1; + uint32_t qs, qe, ts, te; +} uc_block_t; + +typedef struct +{ + uint32_t *a; + ssize_t n ,m; +} N_t; + +typedef struct +{ + char *n_n; + uint32_t n_l; + + kvec_t(uint8_t) r_base; + uint32_t rlen; + + kvec_t(uc_block_t) bb; + N_t N_site; +} ul_vec_t; + +typedef struct +{ + ul_vec_t *a; + size_t n, m; + All_reads *hR; + uint32_t mm; +} all_ul_t; + +extern all_ul_t UL_INF; + 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); @@ -173,4 +208,10 @@ int destory_read_bin(All_reads* r); void init_Debug_reads(Debug_reads* x, const char* file); void destory_Debug_reads(Debug_reads* x); void recover_UC_sub_Read(UC_Read* i_r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID); + +void init_all_ul_t(all_ul_t *x, All_reads *hR); +void destory_all_ul_t(all_ul_t *x, All_reads *hR); +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); +void retrieve_ul_t(UC_Read* r, all_ul_t *ref, uint64_t ID, uint8_t strand); + #endif diff --git a/inter.cpp b/inter.cpp index 5399458..f723f8a 100644 --- a/inter.cpp +++ b/inter.cpp @@ -14,6 +14,8 @@ #include "htab.h" #include "Hash_Table.h" #include "Correct.h" +#include "Process_Read.h" +#include "Assembly.h" KSEQ_INIT(gzFile, gzread) #define MG_SEED_IGNORE (1ULL<<41) @@ -290,6 +292,7 @@ typedef struct { // data structure for each step in kt_pipeline() st_mt_t *sps; mg_gchains_t **gcs; mg_tbuf_t **buf; + ha_ovec_buf_t **hab; } utepdat_t; void init_mg_opt(mg_idxopt_t *opt, int is_HPC, int k, int w, int hap_n) @@ -2249,6 +2252,137 @@ int alignment_ul_pipeline(uldat_t* sl, const enzyme *fn) return 1; } +static void *worker_ul_scall_pipeline(void *data, int step, void *in) // callback for kt_pipeline() +{ + uldat_t *p = (uldat_t*)data; + ///uint64_t total_base = 0, total_pair = 0; + if (step == 0) { // step 1: read a block of sequences + int ret; + uint64_t l; + utepdat_t *s; + CALLOC(s, 1); + s->ha_flt_tab = p->ha_flt_tab; s->ha_idx = p->ha_idx; s->id = p->total_pair; + s->opt = p->opt; s->ug = p->ug; s->uopt = p->uopt; s->rg = p->rg; + while ((ret = kseq_read(p->ks)) >= 0) + { + if (p->ks->seq.l < (uint64_t)p->opt->k) continue; + if (s->n == s->m) { + s->m = s->m < 16? 16 : s->m + (s->n>>1); + REALLOC(s->len, s->m); + REALLOC(s->seq, s->m); + } + + append_ul_t(&UL_INF, NULL, p->ks->name.s, p->ks->name.l, NULL, 0, NULL, 0); + l = p->ks->seq.l; + MALLOC(s->seq[s->n], l); + s->sum_len += l; + memcpy(s->seq[s->n], p->ks->seq.s, l); + s->len[s->n++] = l; + if (s->sum_len >= p->chunk_size) break; + } + p->total_pair += s->n; + if (s->sum_len == 0) free(s); + else return s; + } + else if (step == 1) { // step 2: alignment + utepdat_t *s = (utepdat_t*)in; + + 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); + ///debug + /** + uint64_t i; + CALLOC(s->mzs, p->n_thread); + CALLOC(s->sps, p->n_thread); + CALLOC(s->gcs, s->n); + + s->buf = (mg_tbuf_t**)calloc(p->n_thread, sizeof(mg_tbuf_t*)); + for (i = 0; i < p->n_thread; ++i) s->buf[i] = mg_tbuf_init(); + + kt_for(p->n_thread, worker_for_ul_alignment, s, s->n); + for (i = 0; i < (uint64_t)s->n; ++i) { + free(s->seq[i]); + p->total_base += s->len[i]; + } + free(s->seq); free(s->len); + + for (i = 0; i < p->n_thread; ++i) { + mg_tbuf_destroy(s->buf[i]); + free(s->mzs[i].a); free(s->sps[i].a); + } + **/ + for (i = 0; i < p->n_thread; ++i) ha_ovec_destroy(s->hab[i]); + // free(s->buf); free(s->mzs); free(s->sps); + return s; + } + else if (step == 2) { // step 3: dump + utepdat_t *s = (utepdat_t*)in; + uint64_t i, rid; + for (i = 0; i < (uint64_t)s->n; ++i) { + ///debug + /** + // if(s->pos[i].s == (uint64_t)-1) continue; + // kv_push(pe_hit, p->hits.a, s->pos[i]); + if(!s->gcs[i]) continue; + dump_gaf(&(p->hits), s->gcs[i], 1); + free(s->gcs[i]->gc); free(s->gcs[i]->a); free(s->gcs[i]->lc); free(s->gcs[i]); + **/ + rid = s->id + i; + append_ul_t(&UL_INF, &rid, NULL, 0, s->seq[i], s->len[i], NULL, 0); + free(s->seq[i]); p->total_base += s->len[i]; + } + ///debug + /** + free(s->gcs); + **/ + free(s); + } + return 0; +} + +int scall_ul_pipeline(uldat_t* sl, const enzyme *fn) +{ + double index_time = yak_realtime(); + int i; + init_aux_table(); + + init_all_ul_t(&UL_INF, &R_INF); + for (i = 0; i < fn->n; i++){ + gzFile fp; + if ((fp = gzopen(fn->a[i], "r")) == 0) return 0; + sl->ks = kseq_init(fp); + kt_pipeline(3, worker_ul_scall_pipeline, sl, 3); + kseq_destroy(sl->ks); + gzclose(fp); + } + 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); + + return 1; +} + + +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, U_INF, i, 0); + fprintf(stderr, ">%s\n", p->n_n); + fprintf(stderr, "%.*s\n", (int)ur.length, ur.seq); + } + + destory_UC_Read(&ur); + return 1; +} + inline void get_ulname(mg_dbn_t *name, int32_t rid, char **rn, int32_t *rl) { (*rn) = name->cc.a + (rid>0?name->a[rid-1]:0); @@ -2728,4 +2862,33 @@ void ul_resolve(ma_ug_t *ug, const asg_t *rg, const ug_opt_t *uopt, int hap_n) if(exist == 0) uidx_write(ha_flt_tab, ha_idx, asm_opt.output_file_name); ul_align(&opt, uopt, rg, asm_opt.ar, ha_flt_tab, ha_idx, ug); uidx_destory(); +} + +int ul_v_call(mg_idxopt_t *opt, const ug_opt_t *uopt, const asg_t *rg, const enzyme *fn, void *ha_flt_tab, ha_pt_t *ha_idx, ma_ug_t *ug) +{ + uldat_t sl; memset(&sl, 0, sizeof(sl)); + sl.ha_flt_tab = ha_flt_tab; + sl.ha_idx = ha_idx; + sl.opt = opt; + sl.chunk_size = 500000000; + sl.n_thread = asm_opt.thread_num; + sl.ug = ug; + sl.rg = rg; + sl.uopt = uopt; + scall_ul_pipeline(&sl, fn); + print_ul_rs(&UL_INF); + // 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); + // } + + return 1; +} + +void ul_load(const ug_opt_t *uopt) +{ + fprintf(stderr, "[M::%s::] ==> UL\n", __func__); + mg_idxopt_t opt; + init_mg_opt(&opt, 0, 19, 10, 4095); + ul_v_call(&opt, uopt, /**rg**/NULL, asm_opt.ar, ha_flt_tab, ha_idx, /**ug**/NULL); } \ No newline at end of file diff --git a/inter.h b/inter.h index 6e65d07..d15cce8 100644 --- a/inter.h +++ b/inter.h @@ -3,5 +3,6 @@ #include "Overlaps.h" void ul_resolve(ma_ug_t *ug, const asg_t *rg, const ug_opt_t *uopt, int hap_n); +void ul_load(const ug_opt_t *uopt); #endif