From ab073aa1fc4df0619d58d7e659306c8aaff7b17b Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Sat, 14 Nov 2020 23:35:32 -0500 Subject: [PATCH] hic first --- Overlaps.cpp | 89 +- Trio.cpp | 13 +- hic.cpp | 2369 +++++++++++++++++++++++++++++++++++++++++++++----- htab.cpp | 11 +- khashl.h | 11 +- 5 files changed, 2266 insertions(+), 227 deletions(-) diff --git a/Overlaps.cpp b/Overlaps.cpp index 2d563f3..bbea098 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -11630,18 +11630,74 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) kv_destroy(new_rtg_edges.a); } +void classify_untigs(ma_ug_t *ug, asg_t *sg, ma_sub_t* coverage_cut, +ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, +kvec_asg_arc_t_warp* new_rtg_edges, int max_hang, int min_ovlp) +{ + uint64_t i, dip_thres; + uint8_t* primary_flag = (uint8_t*)calloc(sg->n_seq, sizeof(uint8_t)); + int tmp_cov = asm_opt.hom_global_coverage; + asm_opt.hom_global_coverage = -1; + purge_dups(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, new_rtg_edges, + asm_opt.purge_simi_rate, asm_opt.purge_overlap_len, max_hang, min_ovlp, 0, 0, 0, 1); + dip_thres = ((double)asm_opt.hom_global_coverage)/((double)HOM_PEAK_RATE)*0.6; + asm_opt.hom_global_coverage = tmp_cov; + + for (i = 0; i < ug->g->n_seq; i++) + { + if(get_ug_coverage(&ug->u.a[i], sg, coverage_cut, sources, ruIndex, primary_flag)g->seq[i].c = 1; + } + else + { + ug->g->seq[i].c = 0; + } + } + free(primary_flag); + + ///fprintf(stderr, "[M::%s] diploid coverage threshold: %lu\n", __func__, dip_thres); +} void output_hic_graph(asg_t *sg, ma_sub_t* coverage_cut, char* output_file_name, -ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) -{ +ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, R_to_U* ruIndex, int max_hang, +int min_ovlp) +{ kvec_asg_arc_t_warp new_rtg_edges; kv_init(new_rtg_edges.a); ma_ug_t *ug = NULL; ug = ma_ug_gen(sg); ma_ug_seq(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, max_hang, min_ovlp); - hic_analysis(ug); + + + + + /** + fprintf(stderr, "Writing raw unitig GFA to disk... \n"); + char* gfa_name = (char*)malloc(strlen(output_file_name)+25); + sprintf(gfa_name, "%s.r_utg.gfa", output_file_name); + FILE* output_file = fopen(gfa_name, "w"); + ma_ug_print(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); + fclose(output_file); + sprintf(gfa_name, "%s.r_utg.noseq.gfa", output_file_name); + output_file = fopen(gfa_name, "w"); + ma_ug_print_simple(ug, &R_INF, sg, coverage_cut, sources, ruIndex, "utg", output_file); + fclose(output_file); + if(asm_opt.bed_inconsist_rate != 0) + { + sprintf(gfa_name, "%s.r_utg.lowQ.bed", output_file_name); + output_file = fopen(gfa_name, "w"); + ma_ug_print_bed(ug, sg, &R_INF, coverage_cut, sources, &new_rtg_edges, + max_hang, min_ovlp, asm_opt.bed_inconsist_rate, "utg", output_file); + fclose(output_file); + } + free(gfa_name); + **/ + classify_untigs(ug, sg, coverage_cut, sources, reverse_sources, ruIndex, &new_rtg_edges, + max_hang, min_ovlp); + hic_analysis(ug); ma_ug_destroy(ug); kv_destroy(new_rtg_edges.a); } @@ -27426,20 +27482,6 @@ ma_sub_t **coverage_cut_ptr, int debug_g) asg_arc_del_simple_circle_untig(sources, coverage_cut, sg, 100, 0); - if (asm_opt.flag & HA_F_VERBOSE_GFA) - { - /*******************************for debug***************************************/ - write_debug_graph(sg, sources, coverage_cut, output_file_name, n_read, reverse_sources, ruIndex); - debug_gfa:; - /*******************************for debug***************************************/ - } - /** - debug_ma_hit_t(sources, coverage_cut, n_read, max_hang_length, - mini_overlap_length); - debug_ma_hit_t(reverse_sources, coverage_cut, n_read, max_hang_length, - mini_overlap_length); - **/ - ///note: don't apply asg_arc_del_too_short_overlaps() after this function!!!! rescue_contained_reads_aggressive(NULL, sg, sources, coverage_cut, ruIndex, max_hang_length, mini_overlap_length, bubble_dist, 10, 1, 0, NULL, NULL); @@ -27452,6 +27494,14 @@ ma_sub_t **coverage_cut_ptr, int debug_g) // rescue_no_coverage_aggressive(sg, sources, reverse_sources, &coverage_cut, ruIndex, max_hang_length, // mini_overlap_length, bubble_dist, 10); + if (asm_opt.flag & HA_F_VERBOSE_GFA) + { + /*******************************for debug***************************************/ + write_debug_graph(sg, sources, coverage_cut, output_file_name, n_read, reverse_sources, ruIndex); + debug_gfa:; + /*******************************for debug***************************************/ + } + if (ha_opt_triobin(&asm_opt)) { char *buf = (char*)calloc(strlen(output_file_name) + 25, 1); @@ -27468,7 +27518,10 @@ ma_sub_t **coverage_cut_ptr, int debug_g) } else if(ha_opt_hic(&asm_opt)) { - output_hic_graph(sg, coverage_cut, output_file_name, sources, ruIndex, max_hang_length, mini_overlap_length);; + char *buf = (char*)calloc(strlen(output_file_name) + 25, 1); + sprintf(buf, "%s.hic", output_file_name); + output_hic_graph(sg, coverage_cut, buf, sources, reverse_sources, ruIndex, max_hang_length, mini_overlap_length);; + free(buf); } else { diff --git a/Trio.cpp b/Trio.cpp index b4934a9..75f20f4 100644 --- a/Trio.cpp +++ b/Trio.cpp @@ -67,7 +67,7 @@ static yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, .. { va_list ap; FILE *fp; - uint32_t t[3]; + uint32_t t[3], f_tmp = 0; char magic[4]; int i, j, absent, min_cnt = 0, mid_cnt = 0, mode_err = 0; uint64_t mask = (1ULL<k && (int)t[1] == ch->pre); for (i = 0; i < 1<pre; ++i) { yak_ht_t *h = ch->h[i].h; - fread(t, 4, 2, fp); + f_tmp += fread(t, 4, 2, fp); + ///t[0] = kh_capacity(h), t[1] = kh_size(h); if (ch0 == 0) yak_ht_resize(h, t[0]); for (j = 0; j < (int)t[1]; ++j) { uint64_t key; - fread(&key, 8, 1, fp); + f_tmp += fread(&key, 8, 1, fp); if (mode == YAK_LOAD_ALL) { ++n_ins; yak_ht_put(h, key, &absent); if (absent) ++n_new; } else if (mode == YAK_LOAD_TRIOBIN1 || mode == YAK_LOAD_TRIOBIN2) { int cnt = key & mask, x, shift = mode == YAK_LOAD_TRIOBIN1? 0 : 2; + //1. filter singleton k-mer; 2. label non-repeat and repeat if (cnt >= mid_cnt) x = 2<= min_cnt) x = 1<= 0) { khint_t k; + ///no need cnt at all key = (key & ~mask) | x; ++n_ins; k = yak_ht_put(h, key, &absent); diff --git a/hic.cpp b/hic.cpp index 4a44350..4f507de 100644 --- a/hic.cpp +++ b/hic.cpp @@ -3,10 +3,67 @@ #include "htab.h" #include "assert.h" #include "Overlaps.h" +#include "Hash_Table.h" +#include "Correct.h" #include "khashl.h" +#include "kthread.h" +#include "ksort.h" +#include "kseq.h" // FASTA/Q parser +#include "kdq.h" +KSEQ_INIT(gzFile, gzread) +KDQ_INIT(uint64_t) +#define kdq_clear(q) ((q)->count = (q)->front = 0) +#define kv_malloc(v, s) ((v).n = 0, (v).m = (s), MALLOC((v).a, (s))) -#define generic_hc_key(x) (x) -KHASHL_MAP_INIT(static klib_unused, hc_pt_t, hc_pt, uint64_t, uint64_t, generic_hc_key, kh_eq_generic) +#define HIC_COUNTER_BITS 12 +#define HIC_MAX_COUNT ((1<>HIC_COUNTER_BITS == (b)>>HIC_COUNTER_BITS) +#define hic_ct_hash(a) ((a)>>HIC_COUNTER_BITS) +KHASHL_MAP_INIT(static klib_unused, hc_pt_t, hc_pt, uint64_t, uint64_t, hic_ct_hash, hic_ct_eq) + +typedef struct{ + kvec_t(char) name; + kvec_t(uint64_t) name_Len; + kvec_t(char) r; + kvec_t(uint64_t) r_Len; +} reads_t; + +typedef struct{ + int weight; + uint32_t uID:31, del:1; + uint32_t enzyme; +} hc_edge; + +typedef struct{ + kvec_t(hc_edge) e; + kvec_t(hc_edge) f;//forbiden +} hc_linkeage; + +typedef struct{ + kvec_t(hc_linkeage) a; +} hc_links; + +typedef struct{ + ///kvec_t(hc_edge) a; + size_t n, m; + hc_edge *a; +}hc_edge_warp; + +typedef struct{ + kvec_t(hc_edge_warp) rGraph; + kvec_t(uint64_t) order; + kvec_t(uint8_t) rGraphSet; + kvec_t(uint8_t) rGraphVis; + kvec_t(uint8_t) utgVis; + kdq_t(uint64_t) *q; + kvec_t(uint64_t) parent; + uint64_t uID_mode, uID_shift, n, src, dest; +} min_cut_t; typedef struct { hc_pt_t *h; @@ -15,6 +72,13 @@ typedef struct { khint_t end; } hc_pt1_t; +typedef struct { + uint32_t* index; + ma_ug_t* ug; + kvec_t(uint32_t) list; + kvec_t(uint32_t) num; +} bubble_type; + typedef struct { ma_ug_t* ug; uint64_t uID_bits; @@ -22,12 +86,132 @@ typedef struct { uint64_t pos_bits; uint64_t pos_mode; uint64_t rev_mode; - hc_pt1_t idx; uint64_t k; + uint64_t max_cnt; + + + uint64_t pre; + uint64_t tot; + uint64_t tot_pos; + hc_pt1_t* idx_buf; } ha_ug_index; +typedef struct { // data structure for each step in kt_pipeline() + uint64_t key, pos; +} ch_buf_t; + +typedef struct { + kvec_t(uint64_t) a; +} kvec_cnt; + +typedef struct { + kvec_t(ch_buf_t) a; +} kvec_pos; + +typedef struct { // global data structure for kt_pipeline() + int is_cnt; + uint64_t buf_bytes; + ha_ug_index *h; + kvec_cnt* cnt; + kvec_pos* buf; + uint64_t n_thread; +} pldat_t; + +typedef struct { + uint64_t s, e, id; +} pe_hit; + +typedef struct { + kvec_t(pe_hit) a; +} kvec_pe_hit; + +#define pe_hit_an1_key(x) ((x).s<<1) +KRADIX_SORT_INIT(pe_hit_an1, pe_hit, pe_hit_an1_key, 8) +#define pe_hit_an2_key(x) ((x).e<<1) +KRADIX_SORT_INIT(pe_hit_an2, pe_hit, pe_hit_an2_key, 8) +#define generic_key(x) (x) +KRADIX_SORT_INIT(hc64, uint64_t, generic_key, 8) + +typedef struct { // global data structure for kt_pipeline() + const ha_ug_index* idx; + kseq_t *ks1, *ks2; + int64_t chunk_size; + uint64_t n_thread; + uint64_t total_base; + uint64_t total_pair; + kvec_pe_hit hits; +} sldat_t; + +typedef struct { + uint64_t ref; + uint64_t off_cnt; +} s_hit; + +typedef struct { + kvec_t(s_hit) a; +} kvec_vote; + +typedef struct { // data structure for each step in kt_pipeline() + const ha_ug_index* idx; + int n, m, sum_len; + uint64_t *len, id; + char **seq; + ch_buf_t *buf; + kvec_vote* pos_buf; + pe_hit* pos; +} stepdat_t; + +#define generic_key(x) (x) +KRADIX_SORT_INIT(b64, uint64_t, generic_key, 8) +#define ch_buf_t_key(a) ((a).key) +KRADIX_SORT_INIT(ch_buf, ch_buf_t, ch_buf_t_key, member_size(ch_buf_t, key)) +#define hc_pos_key(x) ((x)<<1) +KRADIX_SORT_INIT(hc_pos, uint64_t, hc_pos_key, 8) +#define hc_s_hit_an1_key(a) ((a).ref) +KRADIX_SORT_INIT(hc_s_hit_an1, s_hit, hc_s_hit_an1_key, 8) +#define hc_s_hit_an2_key(a) ((uint32_t)(a).off_cnt) +KRADIX_SORT_INIT(hc_s_hit_an2, s_hit, hc_s_hit_an2_key, 8) + +#define Get_bub_num(RECORD) ((RECORD).num.n-1) + +reads_t R1, R2; ha_ug_index* ug_index; +void init_ha_ug_index_opt(ha_ug_index* idx, ma_ug_t *ug, int k, pldat_t* p) +{ + uint64_t i, n; + for (idx->uID_bits=1; (uint64_t)(1<uID_bits)<(uint64_t)ug->u.n; idx->uID_bits++); + idx->pos_bits = 64 - idx->uID_bits - 1; + idx->uID_mode = (((uint64_t)-1) << (64-idx->uID_bits))>>1; + idx->pos_mode = ((uint64_t)-1) >> (64-idx->pos_bits); + idx->rev_mode = ((uint64_t)1) << 63; + idx->ug = ug; + idx->k = k; + idx->pre = HIC_COUNTER_BITS; + idx->tot = 1 << idx->pre; + idx->tot_pos = 0; + CALLOC(idx->idx_buf, idx->tot); + for (i = 0; i < idx->tot; i++) + { + idx->idx_buf[i].h = hc_pt_init(); + } + for (i = n = 0; i < ug->u.n; i++) + { + n += ug->u.a[i].len; + } + n = n << 3; + p->h = idx; + p->buf_bytes = n>>7; + CALLOC(p->cnt, idx->tot); + CALLOC(p->buf, idx->tot); + for (i = 0; i < idx->tot; i++) + { + kv_init(p->cnt[i].a); + kv_init(p->buf[i].a); + } + p->n_thread = asm_opt.thread_num; +} + inline uint64_t get_k_direction(uint64_t x[4]) { if(x[1] != x[3]) @@ -44,98 +228,38 @@ inline uint64_t get_k_direction(uint64_t x[4]) } } -inline uint64_t hc_hash_long(uint64_t x[4], uint64_t* skip) +inline uint64_t hc_hash_long(uint64_t x[4], uint64_t* skip, uint64_t k) { ///compare forward k-mer and reverse complementary strand (*skip) = get_k_direction(x); if((*skip) == (uint64_t)-1) return (*skip); + if (k <= 32) return ((x[(*skip)<<1|0]<<32)|(x[(*skip)<<1|1])); return yak_hash64_64(x[(*skip)<<1|0]) + yak_hash64_64(x[(*skip)<<1|1]); } -inline uint64_t get_hc_pt1_count(ha_ug_index* idx, uint64_t key, uint64_t** pos_list) +inline uint64_t get_hc_pt1_count(ha_ug_index* index, uint64_t key, uint64_t** pos_list) { + uint64_t bucket_mask = (1ULL<pre) - 1; + hc_pt1_t* h = &(index->idx_buf[key & bucket_mask]); uint64_t beg; khint_t k; - k = hc_pt_get(idx->idx.h, key); - if (k == kh_end(idx->idx.h)) + k = hc_pt_get(h->h, key); + if (k == kh_end(h->h)) { return 0; } - beg = kh_val(idx->idx.h, k); - if(pos_list) *pos_list = idx->idx.a + beg; - if(k == idx->idx.end) return idx->idx.n - beg; - for (k++; k != kh_end(idx->idx.h); ++k) + beg = kh_val(h->h, k); + if(pos_list) *pos_list = h->a + beg; + if((kh_key(h->h, k)&HIC_MAX_COUNT)h, k)&HIC_MAX_COUNT; + if(k == h->end) return h->n - beg; + for (k++; k != kh_end(h->h); ++k) { - if (kh_exist(idx->idx.h, k)) + if (kh_exist(h->h, k)) { - return kh_val(idx->idx.h, k) - beg; + return kh_val(h->h, k) - beg; } } - return idx->idx.n - beg; -} - -void count_hc_pt1(char* seq, uint64_t len, ha_ug_index* idx) -{ - uint64_t i, l; - khint_t key; - int absent; - uint64_t x[4], mask = (1ULL<k) - 1, shift = idx->k - 1, hash, skip; - for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) { - int c = seq_nt4_table[(uint8_t)seq[i]]; - ///c = 00, 01, 10, 11 - if (c < 4) { // not an "N" base - ///x[0] & x[1] are the forward k-mer - ///x[2] & x[3] are the reverse complementary k-mer - x[0] = (x[0] << 1 | (c&1)) & mask; - x[1] = (x[1] << 1 | (c>>1)) & mask; - x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift; - x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; - if (++l >= idx->k) - { - hash = hc_hash_long(x, &skip); - if(skip == (uint64_t)-1) continue; - key = hc_pt_put(idx->idx.h, hash, &absent); - if(absent) kh_val(idx->idx.h, key) = 0; - kh_val(idx->idx.h, key)++; - } - - } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart - } -} - -void fill_hc_pt1(char* seq, uint64_t len, uint64_t uID, ha_ug_index* idx) -{ - uint64_t i, l, pos, *pos_list = NULL, cnt; - uint64_t x[4], mask = (1ULL<k) - 1, shift = idx->k - 1, hash, skip; - for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) { - int c = seq_nt4_table[(uint8_t)seq[i]]; - ///c = 00, 01, 10, 11 - if (c < 4) { // not an "N" base - ///x[0] & x[1] are the forward k-mer - ///x[2] & x[3] are the reverse complementary k-mer - x[0] = (x[0] << 1 | (c&1)) & mask; - x[1] = (x[1] << 1 | (c>>1)) & mask; - x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift; - x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; - if (++l >= idx->k) - { - hash = hc_hash_long(x, &skip); - if(skip == (uint64_t)-1) continue; - pos = (skip << 63) | ((uID << (64-idx->uID_bits))>>1) | (i & idx->pos_mode); - cnt = get_hc_pt1_count(idx, hash, &pos_list); - ///assert(cnt != 0); - if(pos_list[cnt-1]!=cnt-1) - { - pos_list[pos_list[cnt-1]]=pos; - pos_list[cnt-1]++; - } - else - { - pos_list[pos_list[cnt-1]]=pos; - } - } - } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart - } + return h->n - beg; } void test_hc_pt1(char* seq, uint64_t len, uint64_t uID, ha_ug_index* idx) @@ -154,7 +278,7 @@ void test_hc_pt1(char* seq, uint64_t len, uint64_t uID, ha_ug_index* idx) x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; if (++l >= idx->k) { - hash = hc_hash_long(x, &skip); + hash = hc_hash_long(x, &skip, idx->k); if(skip == (uint64_t)-1) continue; pos = (skip << 63) | ((uID << (64-idx->uID_bits))>>1) | (i & idx->pos_mode); cnt = get_hc_pt1_count(idx, hash, &pos_list); @@ -174,68 +298,37 @@ void test_hc_pt1(char* seq, uint64_t len, uint64_t uID, ha_ug_index* idx) } } -void test_unitig_index(ha_ug_index* idx) +void test_unitig_index(ha_ug_index* idx, ma_ug_t *ug) { - uint32_t i; + double index_time = yak_realtime(); + uint32_t i, j; ma_utg_t *u = NULL; + hc_pt1_t *h = NULL; + idx->ug = ug; for (i = 0; i < idx->ug->u.n; i++) { - ///fprintf(stderr, "i: %u, n: %u\n", i, (uint32_t)idx->ug->u.n); u = &(idx->ug->u.a[i]); if(u->m == 0) continue; test_hc_pt1(u->s, u->len, i, idx); } - for (i = 0; i < idx->idx.n; i++) + + for (i = 0; i < idx->tot; i++) { - if(idx->idx.a[i]!=(uint64_t)-1) + h = &(idx->idx_buf[i]); + for (j = 0; j < h->n; j++) { - fprintf(stderr, "ERROR i\n"); - } - } -} - -void test_count1(char* seq, uint64_t len, uint64_t uID, ha_ug_index* idx, uint64_t direct_cnt) -{ - uint64_t i, l, cnt; - khint_t k; - uint64_t x[4], mask = (1ULL<k) - 1, shift = idx->k - 1, hash, skip; - for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) { - int c = seq_nt4_table[(uint8_t)seq[i]]; - ///c = 00, 01, 10, 11 - if (c < 4) { // not an "N" base - ///x[0] & x[1] are the forward k-mer - ///x[2] & x[3] are the reverse complementary k-mer - x[0] = (x[0] << 1 | (c&1)) & mask; - x[1] = (x[1] << 1 | (c>>1)) & mask; - x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift; - x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; - if (++l >= idx->k) + if(h->a[j] != (uint64_t)-1) { - hash = hc_hash_long(x, &skip); - if(skip == (uint64_t)-1) continue; - if(!direct_cnt) - { - cnt = get_hc_pt1_count(idx, hash, NULL); - } - else - { - k = hc_pt_get(idx->idx.h, hash); - if (k == kh_end(idx->idx.h)) - { - cnt = 0; - } - else - { - cnt = kh_val(idx->idx.h, k); - } - - } + fprintf(stderr, "ERROR j\n"); } - } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart + } + } + + fprintf(stderr, "[M::%s::%.3f] ==> Test has been passed\n", __func__, yak_realtime()-index_time); } -void hc_pt_t_gen(hc_pt1_t* pt) +void hc_pt_t_gen_single(hc_pt1_t* pt) { khint_t k; uint64_t c; @@ -250,8 +343,6 @@ void hc_pt_t_gen(hc_pt1_t* pt) CALLOC(pt->a, pt->n); } - - int write_hc_pt_index(ha_ug_index* idx, char* file_name) { char* gfa_name = (char*)malloc(strlen(file_name)+25); @@ -268,11 +359,18 @@ int write_hc_pt_index(ha_ug_index* idx, char* file_name) fwrite(&idx->pos_mode, sizeof(idx->pos_mode), 1, fp); fwrite(&idx->rev_mode, sizeof(idx->rev_mode), 1, fp); fwrite(&idx->k, sizeof(idx->k), 1, fp); - fwrite(&idx->idx.n, sizeof(idx->idx.n), 1, fp); - fwrite(&idx->idx.end, sizeof(idx->idx.end), 1, fp); - fwrite(idx->idx.a, sizeof(uint64_t), idx->idx.n, fp); - hc_pt_save(idx->idx.h, fp); + fwrite(&idx->pre, sizeof(idx->pre), 1, fp); + fwrite(&idx->tot, sizeof(idx->tot), 1, fp); + fwrite(&idx->tot_pos, sizeof(idx->tot_pos), 1, fp); + uint64_t i = 0; + for (i = 0; i < idx->tot; i++) + { + fwrite(&idx->idx_buf[i].n, sizeof(idx->idx_buf[i].n), 1, fp); + fwrite(&idx->idx_buf[i].end, sizeof(idx->idx_buf[i].end), 1, fp); + fwrite(idx->idx_buf[i].a, sizeof(uint64_t), idx->idx_buf[i].n, fp); + hc_pt_save(idx->idx_buf[i].h, fp); + } fprintf(stderr, "[M::%s] Index has been written.\n", __func__); free(gfa_name); @@ -280,9 +378,10 @@ int write_hc_pt_index(ha_ug_index* idx, char* file_name) return 1; } - int load_hc_pt_index(ha_ug_index** r_idx, char* file_name) { + uint64_t flag = 0; + double index_time = yak_realtime(); char* gfa_name = (char*)malloc(strlen(file_name)+25); sprintf(gfa_name, "%s.hc_tlb", file_name); FILE* fp = fopen(gfa_name, "r"); @@ -292,104 +391,1986 @@ int load_hc_pt_index(ha_ug_index** r_idx, char* file_name) } ha_ug_index* idx = NULL; CALLOC(idx, 1); - fread(&idx->uID_bits, sizeof(idx->uID_bits), 1, fp); - fread(&idx->uID_mode, sizeof(idx->uID_mode), 1, fp); - fread(&idx->pos_bits, sizeof(idx->pos_bits), 1, fp); - fread(&idx->pos_mode, sizeof(idx->pos_mode), 1, fp); - fread(&idx->rev_mode, sizeof(idx->rev_mode), 1, fp); - fread(&idx->k, sizeof(idx->k), 1, fp); - fread(&idx->idx.n, sizeof(idx->idx.n), 1, fp); - fread(&idx->idx.end, sizeof(idx->idx.end), 1, fp); - MALLOC(idx->idx.a, idx->idx.n); - fread(idx->idx.a, sizeof(uint64_t), idx->idx.n, fp); - hc_pt_load(&(idx->idx.h), fp); + flag += fread(&idx->uID_bits, sizeof(idx->uID_bits), 1, fp); + flag += fread(&idx->uID_mode, sizeof(idx->uID_mode), 1, fp); + flag += fread(&idx->pos_bits, sizeof(idx->pos_bits), 1, fp); + flag += fread(&idx->pos_mode, sizeof(idx->pos_mode), 1, fp); + flag += fread(&idx->rev_mode, sizeof(idx->rev_mode), 1, fp); + flag += fread(&idx->k, sizeof(idx->k), 1, fp); + flag += fread(&idx->pre, sizeof(idx->pre), 1, fp); + flag += fread(&idx->tot, sizeof(idx->tot), 1, fp); + flag += fread(&idx->tot_pos, sizeof(idx->tot_pos), 1, fp); + MALLOC(idx->idx_buf, idx->tot); + uint64_t i = 0; + for (i = 0; i < idx->tot; i++) + { + flag += fread(&idx->idx_buf[i].n, sizeof(idx->idx_buf[i].n), 1, fp); + flag += fread(&idx->idx_buf[i].end, sizeof(idx->idx_buf[i].end), 1, fp); + MALLOC(idx->idx_buf[i].a, idx->idx_buf[i].n); + flag += fread(idx->idx_buf[i].a, sizeof(uint64_t), idx->idx_buf[i].n, fp); + hc_pt_load(&(idx->idx_buf[i].h), fp); + } + + (*r_idx) = idx; - fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__); free(gfa_name); fclose(fp); + fprintf(stderr, "[M::%s::%.3f] ==> HiC index has been loaded\n", __func__, yak_realtime()-index_time); return 1; } +static void worker_for_sort(void *data, long i, int tid) // callback for kt_for() +{ + pldat_t *pl = (pldat_t*)data; + hc_pt1_t *h = &(pl->h->idx_buf[i]); + khint_t k; + uint64_t beg, cnt = 0; + uint64_t* pos_list; + for (k = 0; k != kh_end(h->h); ++k) { + if (kh_exist(h->h, k)) { + beg = kh_val(h->h, k); + pos_list = h->a + beg; + if((kh_key(h->h, k)&HIC_MAX_COUNT)h, k)&HIC_MAX_COUNT; + } + else if(k == h->end) + { + cnt = h->n - beg; + } + else + { + for (k++; k != kh_end(h->h); ++k) + { + if (kh_exist(h->h, k)) + { + cnt = kh_val(h->h, k) - beg; + break; + } + } + } + radix_sort_hc_pos(pos_list, pos_list+cnt); + } + } + +} + +void hc_pt_t_gen(ha_ug_index* idx, pldat_t* pl) +{ + if(pl == NULL) + { + uint64_t i; + for (i = 0; i < idx->tot; i++) + { + hc_pt_t_gen_single(&(idx->idx_buf[i])); + } + } + else + { + kt_for(pl->n_thread, worker_for_sort, pl, pl->h->tot); + } +} + +static void worker_for(void *data, long i, int tid) // callback for kt_for() +{ + pldat_t *pl = (pldat_t*)data; + hc_pt1_t *h = &(pl->h->idx_buf[i]); + uint64_t m = 0, beg, end, occ; + khint_t key; + int absent; + + + if(pl->is_cnt) + { + uint64_t* cnt = NULL; + if(pl->cnt[i].a.n > 2) radix_sort_b64(pl->cnt[i].a.a, pl->cnt[i].a.a + pl->cnt[i].a.n); + cnt = pl->cnt[i].a.a; + occ = pl->cnt[i].a.n; + for (m = beg = end = 0; m < occ; m++) + { + if(cnt[beg] == cnt[m]) + { + end = m; + } + else + { + key = hc_pt_put(h->h, cnt[beg], &absent); + if(absent) kh_val(h->h, key) = 0; + kh_val(h->h, key) += (end - beg + 1); + kh_key(h->h, key) = (kh_key(h->h, key)&HIC_KEY_MODE)| + (kh_val(h->h, key)h, key):HIC_MAX_COUNT); + beg = end = m; + } + } + if(occ > 0) + { + key = hc_pt_put(h->h, cnt[beg], &absent); + if(absent) kh_val(h->h, key) = 0; + kh_val(h->h, key) += (end - beg + 1); + kh_key(h->h, key) = (kh_key(h->h, key)&HIC_KEY_MODE)| + (kh_val(h->h, key)h, key):HIC_MAX_COUNT); + } + pl->cnt[i].a.n = 0; + } + + if(!pl->is_cnt) + { + ch_buf_t* pos = NULL; + uint64_t num, *pos_list = NULL, k, k_n, pos_k; + if(pl->buf[i].a.n > 2) radix_sort_ch_buf(pl->buf[i].a.a, pl->buf[i].a.a + pl->buf[i].a.n); + pos = pl->buf[i].a.a; + occ = pl->buf[i].a.n; + for (m = beg = end = 0; m < occ; m++) + { + if(pos[beg].key == pos[m].key) + { + end = m; + } + else + { + num = get_hc_pt1_count(pl->h, pos[beg].key, &pos_list); + + k_n=(end-beg+1);pos_k=pos_list[num-1];pos_list[num-1]+=k_n; + for (k = 0; k < k_n; k++) + { + pos_list[pos_k+k] = pos[beg+k].pos; + } + beg = end = m; + } + + } + if(occ > 0) + { + num = get_hc_pt1_count(pl->h, pos[beg].key, &pos_list); + k_n=(end-beg+1);pos_k=pos_list[num-1];pos_list[num-1]+=k_n; + for (k = 0; k < k_n; k++) + { + pos_list[pos_k+k] = pos[beg+k].pos; + } + } + pl->buf[i].a.n = 0; + } +} + +void parallel_count_hc_pt1(pldat_t* pl) +{ + uint64_t i, l = 0, uID, num_pos = 0, pos_thre; + uint64_t x[4], mask = (1ULL<h->k) - 1, shift = pl->h->k - 1, hash, pos, skip, bucket_mask = (1ULL<h->pre) - 1; + ma_utg_t *u = NULL; + ch_buf_t k_pos; + if(pl->is_cnt) l = ((pl->buf_bytes>>3)/pl->h->tot) + 1, pos_thre = pl->buf_bytes>>3; + if(!pl->is_cnt) l = ((pl->buf_bytes>>4)/pl->h->tot) + 1, pos_thre = pl->buf_bytes>>4; + for (i = 0; i < pl->h->tot; i++) + { + if(pl->is_cnt) + { + kv_resize(uint64_t, pl->cnt[i].a, l); + pl->cnt[i].a.n = 0; + } + + if(!pl->is_cnt) + { + kv_resize(ch_buf_t, pl->buf[i].a, l); + pl->buf[i].a.n = 0; + } + } + + + for (uID = 0; uID < pl->h->ug->u.n; uID++) + { + u = &(pl->h->ug->u.a[uID]); + if(u->m == 0) continue; + + for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < u->len; ++i) { + int c = seq_nt4_table[(uint8_t)u->s[i]]; + ///c = 00, 01, 10, 11 + if (c < 4) { // not an "N" base + ///x[0] & x[1] are the forward k-mer + ///x[2] & x[3] are the reverse complementary k-mer + x[0] = (x[0] << 1 | (c&1)) & mask; + x[1] = (x[1] << 1 | (c>>1)) & mask; + x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift; + x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; + if (++l >= pl->h->k) + { + hash = hc_hash_long(x, &skip, pl->h->k); + if(skip == (uint64_t)-1) continue; + if(pl->is_cnt) + { + kv_push(uint64_t, pl->cnt[hash & bucket_mask].a, hash); + } + else + { + pos = (skip << 63) | ((uID << (64-pl->h->uID_bits))>>1) | (i & pl->h->pos_mode); + k_pos.key = hash; k_pos.pos = pos; + kv_push(ch_buf_t, pl->buf[hash & bucket_mask].a, k_pos); + } + num_pos++; + + if(num_pos >= pos_thre) + { + num_pos = 0; + kt_for(pl->n_thread, worker_for, pl, pl->h->tot); + } + } + + } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart + } + } + + if(num_pos > 0) kt_for(pl->n_thread, worker_for, pl, pl->h->tot); + + for (i = 0; i < pl->h->tot; i++) + { + if(pl->cnt[i].a.m > 0) kv_destroy(pl->cnt[i].a), kv_init(pl->cnt[i].a); + if(pl->buf[i].a.m > 0) kv_destroy(pl->buf[i].a), kv_init(pl->buf[i].a); + } +} ha_ug_index* build_unitig_index(ma_ug_t *ug, int k) { - uint32_t i; - ma_utg_t *u = NULL; ha_ug_index* idx = NULL; CALLOC(idx, 1); + pldat_t pl; pl.h = idx; pl.is_cnt = 1; double index_time = yak_realtime(), beg_time; - for (idx->uID_bits=1; (uint64_t)(1<uID_bits)<(uint64_t)ug->u.n; idx->uID_bits++); - idx->pos_bits = 64 - idx->uID_bits - 1; - idx->uID_mode = (((uint64_t)-1) << (64-idx->uID_bits))>>1; - idx->pos_mode = ((uint64_t)-1) >> (64-idx->pos_bits); - idx->rev_mode = ((uint64_t)1) << 63; - idx->ug = ug; - idx->k = k; - idx->idx.h = hc_pt_init(); + init_ha_ug_index_opt(idx, ug, k, &pl); beg_time = yak_realtime(); - for (i = 0; i < idx->ug->u.n; i++) - { - u = &(idx->ug->u.a[i]); - if(u->m == 0) continue; - count_hc_pt1(u->s, u->len, idx); - } + pl.is_cnt = 1; + parallel_count_hc_pt1(&pl); fprintf(stderr, "[M::%s::%.3f] ==> Counting\n", __func__, yak_realtime()-beg_time); - // beg_time = yak_realtime(); - // for (i = 0; i < idx->ug->u.n; i++) - // { - // u = &(idx->ug->u.a[i]); - // if(u->m == 0) continue; - // test_count1(u->s, u->len, i, idx, 1); - // } - // fprintf(stderr, "[M::%s::%.3f] ==> Counting 1\n", __func__, yak_realtime()-beg_time); - beg_time = yak_realtime(); - hc_pt_t_gen(&(idx->idx)); + hc_pt_t_gen(pl.h, NULL); fprintf(stderr, "[M::%s::%.3f] ==> Memory allocating\n", __func__, yak_realtime()-beg_time); - - // beg_time = yak_realtime(); - // for (i = 0; i < idx->ug->u.n; i++) - // { - // u = &(idx->ug->u.a[i]); - // if(u->m == 0) continue; - // test_count1(u->s, u->len, i, idx, 0); - // } - // fprintf(stderr, "[M::%s::%.3f] ==> Counting 2\n", __func__, yak_realtime()-beg_time); - + beg_time = yak_realtime(); + pl.is_cnt = 0; + parallel_count_hc_pt1(&pl); + fprintf(stderr, "[M::%s::%.3f] ==> Filling pos\n", __func__, yak_realtime()-beg_time); beg_time = yak_realtime(); - for (i = 0; i < idx->ug->u.n; i++) - { - u = &(idx->ug->u.a[i]); - if(u->m == 0) continue; - fill_hc_pt1(u->s, u->len, i, idx); - } - fprintf(stderr, "[M::%s::%.3f] ==> Filling pos\n", __func__, yak_realtime()-beg_time); + hc_pt_t_gen(pl.h, &pl); + fprintf(stderr, "[M::%s::%.3f] ==> Sorting pos\n", __func__, yak_realtime()-beg_time); fprintf(stderr, "[M::%s::%.3f] ==> HiC index has been built\n", __func__, yak_realtime()-index_time); return idx; } -void destory_hc_pt_index(ha_ug_index* r_idx) +void destory_hc_pt_index(ha_ug_index* idx) { - if(r_idx->idx.h) hc_pt_destroy(r_idx->idx.h); - if(r_idx->idx.a) free(r_idx->idx.a); + if(idx->idx_buf) + { + uint64_t i = 0; + for (i = 0; i < idx->tot; i++) + { + if(idx->idx_buf[i].a) free(idx->idx_buf[i].a); + if(idx->idx_buf[i].h) hc_pt_destroy(idx->idx_buf[i].h); + } + free(idx->idx_buf); + } +} + +inline void interpret_pos(const ha_ug_index* idx, s_hit *p, uint64_t* rev, uint64_t* uID, +uint64_t* ref_p, uint64_t* self_p, uint64_t* exact_len, uint64_t* total_len) +{ + (*rev) = p->ref>>63; + (*uID) = (p->ref << 1) >> (64 - idx->uID_bits); + (*self_p) = (uint32_t)p->off_cnt; + (*exact_len) = p->off_cnt >> 32; + if(total_len != NULL) + { + (*exact_len) = (p->off_cnt>>32) & ((uint64_t)65535); + (*total_len) = (p->off_cnt>>48) + (*exact_len); + } + if((p->ref & idx->pos_mode)>>(idx->pos_bits - 1)) + { + (*ref_p) = (*self_p) - (p->ref&(idx->pos_mode>>1)); + } + else + { + (*ref_p) = (*self_p) + (p->ref&(idx->pos_mode)); + } +} + + +inline uint64_t check_exact_match(char* a, long long a_beg, long long a_total, char* b, long long b_beg, +long long b_total, long long Len, uint64_t rev, uint64_t dir) +{ + long long i = 0; + if(rev == 0) + { + if(dir == 0) + { + for (i = 0; i < Len && a_beg < a_total && b_beg < b_total; i++) + { + if(a[a_beg++] != b[b_beg++]) return i; + } + } + else + { + for (i = 0; i < Len && a_beg >= 0 && b_beg >= 0; i++) + { + if(a[a_beg--] != b[b_beg--]) return i; + } + } + } + else + { + + if(dir == 0) + { + for (i = 0; i < Len && a_beg < a_total && b_beg < b_total; i++) + { + if(a[a_beg] != b2rc[seq_nt4_table[(uint8_t)b[b_total - b_beg - 1]]]) return i; + a_beg++; b_beg++; + } + } + else + { + for (i = 0; i < Len && a_beg >= 0 && b_beg >= 0; i++) + { + if(a[a_beg] != b2rc[seq_nt4_table[(uint8_t)b[b_total - b_beg - 1]]]) return i; + a_beg--; b_beg--; + } + } + } + + return i; +} + +uint64_t debug_hash_value(char *r, uint64_t end, uint64_t k_mer) +{ + uint64_t i; + uint64_t x[4], mask = (1ULL<>1)) & mask; + x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift; + x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; + } + } + + return hc_hash_long(x, &skip, k_mer); + +} + +inline uint64_t collect_votes(s_hit* a, uint64_t n) +{ + if(n == 0) return 0; + if(n == 1) return (a[0].off_cnt>>32); + long long i = 0; + uint64_t cur_beg, cur_end, beg, end, ovlp = 0, tLen = 0; + cur_end = (uint32_t)a[n-1].off_cnt; + cur_beg = cur_end + 1 - (a[n-1].off_cnt>>32); + + for (i = n - 2; i >= 0; i--) + { + end = (uint32_t)a[i].off_cnt; + beg = end + 1 - (a[i].off_cnt>>32); + if(MAX(cur_beg, beg) <= MIN(cur_end, end)) + { + cur_beg = MIN(cur_beg, beg); + ///cur_end = MAX(cur_end, end); + } + else + { + ovlp += (cur_end + 1 - cur_beg); + cur_beg = beg; + cur_end = end; + } + } + ovlp += (cur_end + 1 - cur_beg); + tLen = (uint32_t)a[n-1].off_cnt + 1 - cur_beg; + tLen = tLen - ovlp; + tLen = tLen << 16; + return ovlp | tLen; +} + +inline void compress_mapped_pos(const ha_ug_index* idx, kvec_vote* buf, uint64_t buf_iter, uint64_t max_i, uint64_t thres) +{ + if(buf_iter >= buf->a.n) + { + buf->a.n = buf_iter; + return; + } + uint64_t rev, uID, ref_p, self_p, eLen, tLen, i, max_beg, max_end, cur_beg, cur_end, ovlp, max_eLen; + uint64_t secondLen = 0, second_i = (uint64_t)-1; + interpret_pos((ha_ug_index*)idx, &buf->a.a[max_i], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + max_end = self_p; + max_beg = self_p + 1 - tLen; + max_eLen = eLen; + for (i = buf_iter; i < buf->a.n; i++) + { + if(i == max_i) continue; + interpret_pos(idx, &buf->a.a[i], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + cur_end = self_p; + cur_beg = self_p + 1 - tLen; + if(MAX(cur_beg, max_beg) <= MIN(cur_end, max_end)) + { + ovlp = MIN(cur_end, max_end) - MAX(cur_beg, max_beg) + 1; + if(ovlp > thres) + { + if(eLen >= max_eLen * 0.8) + { + buf->a.n = buf_iter; + return; + } + continue; + } + } + if(secondLen < eLen) secondLen = eLen, second_i = i; + } + + if(second_i == (uint64_t)-1) + { + buf->a.a[buf_iter] = buf->a.a[max_i]; + buf->a.n = buf_iter + 1; + } + else + { + buf->a.a[buf_iter] = buf->a.a[MIN(max_i, second_i)]; + buf->a.a[buf_iter+1] = buf->a.a[MAX(max_i, second_i)]; + buf->a.n = buf_iter + 2; + } +} + +inline void print_pos_list(const ha_ug_index* idx, s_hit *l, uint64_t occ, uint64_t rid, uint64_t r1) +{ + if(rid == 33045391 || rid == 4239289 || rid == 5267597 || rid == 34474764 || rid == 35016489 + || rid == 36002255 || rid == 37811694 || rid == 46805824) + { + uint64_t i, rev, uID, ref_p, self_p, cnt; + for (i = 0; i < occ; i++) + { + interpret_pos(idx, &l[i], &rev, &uID, &ref_p, &self_p, &cnt, NULL); + fprintf(stderr, "(r%lu) rid: %lu, i: %lu, rev: %lu, uID: %lu, ref_p: %lu, self_p: %lu\n", + r1, rid, i, rev, uID, ref_p, self_p); + } + + } +} + +void get_alignment(char *r, uint64_t len, uint64_t k_mer, kvec_vote* buf, +const ha_ug_index* idx, uint64_t buf_iter, uint64_t rid) +{ + uint64_t i, j, l = 0, skip, *pos_list = NULL, cnt, rev, self_p, ref_p, u_len, uID; + uint64_t x[4], mask = (1ULL<a.n = 0; + for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) { + int c = seq_nt4_table[(uint8_t)r[i]]; + ///c = 00, 01, 10, 11 + if (c < 4) { // not an "N" base + ///x[0] & x[1] are the forward k-mer + ///x[2] & x[3] are the reverse complementary k-mer + x[0] = (x[0] << 1 | (c&1)) & mask; + x[1] = (x[1] << 1 | (c>>1)) & mask; + x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift; + x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; + if (++l >= k_mer) + { + hash = hc_hash_long(x, &skip, k_mer); + if(skip == (uint64_t)-1) continue; + /*******************************for debug************************************/ + // if(debug_hash_value(r, i, k_mer) != hash) + // { + // fprintf(stderr, "ERROR\n"); + // } + /*******************************for debug************************************/ + cnt = get_hc_pt1_count((ha_ug_index*)idx, hash, &pos_list); + if(cnt > idx->max_cnt) continue; + if(cnt != 1) continue; ///might be able to be disabled in future + + + for (j = 0; j < cnt; j++) + { + kv_pushp(s_hit, buf->a, &p); + rev = (pos_list[j]>>63) != skip; + self_p = i; + ref_p = pos_list[j] & idx->pos_mode; + uID = (pos_list[j] << 1) >> (64 - idx->uID_bits); + u_len = idx->ug->u.a[uID].len; + if(rev) ref_p = u_len - 1 - (ref_p + 1 - k_mer); + ///p->off_cnt = self_p | ((uint64_t)cnt << 32); + p->off_cnt = self_p | ((uint64_t)k_mer << 32); + + p->ref = ref_p >= self_p? (ref_p-self_p) + : (self_p-ref_p) + ((uint64_t)1 << (idx->pos_bits - 1)); + p->ref = (rev << 63)|(pos_list[j] & idx->uID_mode)|(p->ref&idx->pos_mode); + + + /*******************************for debug************************************/ + // if(check_exact_match(r, i + 1 - k_mer, len, + // idx->ug->u.a[uID].s, ref_p + 1 - k_mer, u_len, k_mer, rev, 0) != k_mer + // || + // check_exact_match(r, i, len, + // idx->ug->u.a[uID].s, ref_p, u_len, k_mer, rev, 1) != k_mer) + // { + // fprintf(stderr, "ERROR\n"); + // } + /*******************************for debug************************************/ + } + + if(cnt == 1) + { + ///uint64_t debug_right = 0, debug_left = 0, debug_len; + + j = check_exact_match(r, self_p + 1, len, idx->ug->u.a[uID].s, ref_p + 1, u_len, len, rev, 0); + + ///debug_right = j; + + ///if(j == 0) continue; + if((j + 1) >= k_mer) + { + l = 0, x[0] = x[1] = x[2] = x[3] = 0; + i = i + j - (k_mer - 1); + } + else + { + ///l = i - (i + j - (k_mer - 1)); + l = k_mer - j -1; + } + buf->a.a[buf->a.n-1].off_cnt += ((uint64_t)j << 32) + j; + + if(self_p >= k_mer && ref_p >= k_mer) + { + j = check_exact_match(r, self_p - k_mer, len, idx->ug->u.a[uID].s, + ref_p - k_mer, u_len, len, rev, 1); + buf->a.a[buf->a.n-1].off_cnt += ((uint64_t)j << 32); + ///debug_left = j; + } + + + // debug_len = check_exact_match(r, self_p + debug_right, len, idx->ug->u.a[uID].s, + // ref_p + debug_right, u_len, len, rev, 1); + // if(debug_len!= (debug_left + debug_right + k_mer)) + // { + // fprintf(stderr, "debug_len: %lu, debug_left: %lu, debug_right: %lu\n", + // debug_len, debug_left, debug_right); + // } + } + + } + + } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart + } + + ///if(buf->a.n - buf_iter <= 1) return; + if(buf->a.n - buf_iter == 0) return; + if(buf->a.n - buf_iter > 1) radix_sort_hc_s_hit_an1(buf->a.a + buf_iter, buf->a.a + buf->a.n); + + + + + /*******************************for debug************************************/ + // print_pos_list(idx, buf->a.a+buf_iter, buf->a.n - buf_iter, rid, (buf_iter != 0)); + // fprintf(stderr, "len0:%lu\n", buf->a.n - buf_iter); + // for (i = buf_iter; i < buf->a.n; i++) + // { + // interpret_pos(idx, &buf->a.a[i], &rev, &uID, &ref_p, &self_p, &cnt, NULL); + // fprintf(stderr, "(%lu) rev: %lu, uID: %lu, ref_p: %lu, self_p: %lu, len: %lu\n", + // i, rev, uID, ref_p, self_p, cnt); + // } + /*******************************for debug************************************/ + + + + + + uint64_t cur_ref_p, thres = (len * HIC_R_E_RATE) + 1, m, index_beg, ovlp, maxLen = 0, max_i = (uint64_t)-1; + i = m = buf_iter; + while (i < buf->a.n) + { + interpret_pos(idx, &buf->a.a[i], &rev, &uID, &ref_p, &self_p, &cnt, NULL); + /*******************************for debug************************************/ + // if(check_exact_match(r, self_p, len, idx->ug->u.a[uID].s, + // ref_p, idx->ug->u.a[uID].len, cnt, rev, 1) != cnt) + // { + // fprintf(stderr, "ERROR\n"); + // } + /*******************************for debug************************************/ + // if(self_p > ref_p) + // { + // i++; + // continue; ///fix this in future + // } + cur_ref_p = buf->a.a[i].ref; + index_beg = i; + while ((i < buf->a.n) && ((buf->a.a[i].ref>>idx->pos_bits) == (cur_ref_p>>idx->pos_bits)) && + (buf->a.a[i].ref - cur_ref_p <= thres)) + { + i++; + } + if(i - index_beg > 1) + { + radix_sort_hc_s_hit_an2(buf->a.a + index_beg, buf->a.a + i); + } + ovlp = collect_votes(buf->a.a + index_beg, i - index_beg); + buf->a.a[m] = buf->a.a[i - 1]; + buf->a.a[m].off_cnt = (buf->a.a[m].off_cnt << 32)>>32; + buf->a.a[m].off_cnt += ((uint64_t)ovlp<<32); + + if(maxLen < (ovlp&((uint64_t)65535))) maxLen = (ovlp&((uint64_t)65535)), max_i = m; + + m++; + } + buf->a.n = m; + + /*******************************for debug************************************/ + // for (i = buf_iter; i < buf->a.n; i++) + // { + // uint64_t eLen, tLen; + // interpret_pos(idx, &buf->a.a[i], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + // if(maxLen < eLen) fprintf(stderr, "ERROR1\n"); + // if(i == max_i && maxLen != eLen) fprintf(stderr, "ERROR2\n"); + // } + /*******************************for debug************************************/ + ///select the best alignment at [buf_iter, m) + + /*******************************for debug************************************/ + // fprintf(stderr, "len1:%lu, max_i: %lu\n", buf->a.n - buf_iter, max_i); + // for (i = buf_iter; i < buf->a.n; i++) + // { + // uint64_t eLen, tLen; + // interpret_pos(idx, &buf->a.a[i], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + // fprintf(stderr, "(%lu) rev: %lu, uID: %lu, ref_p: %lu, self_p: %lu, eLen: %lu, tLen: %lu\n", + // i, rev, uID, ref_p, self_p, eLen, tLen); + // } + /*******************************for debug************************************/ + + compress_mapped_pos(idx, buf, buf_iter, max_i, thres); + + /*******************************for debug************************************/ + // fprintf(stderr, "len2:%lu, max_i: %lu\n", buf->a.n - buf_iter, max_i); + // for (i = buf_iter; i < buf->a.n; i++) + // { + // uint64_t eLen, tLen; + // interpret_pos(idx, &buf->a.a[i], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + // fprintf(stderr, "(%lu) rev: %lu, uID: %lu, ref_p: %lu, self_p: %lu, eLen: %lu, tLen: %lu\n", + // i, rev, uID, ref_p, self_p, eLen, tLen); + // } + // if(buf->a.n != m) fprintf(stderr, "Changed\n"); + // fprintf(stderr, "\n"); + /*******************************for debug************************************/ +} + +inline void set_pe_pos(ha_ug_index* idx, s_hit *l1, uint64_t occ1, s_hit *l2, uint64_t occ2, +pe_hit* x, uint64_t rid) +{ + if(occ1 == 0 || occ2 == 0) return; + uint64_t rev1, rev2, uID1, uID2, ref_p1, ref_p2, self_p1, self_p2, eLen1, eLen2, tLen1, tLen2; + uint64_t rev_t, uID_t, ref_p_t, self_p_t, eLen_t, tLen_t; + ///5' end of r1 and r2 + interpret_pos(idx, &l1[0], &rev1, &uID1, &ref_p1, &self_p1, &eLen1, &tLen1); + interpret_pos(idx, &l2[0], &rev2, &uID2, &ref_p2, &self_p2, &eLen2, &tLen2); + /*******************************for debug************************************/ + // if(rid == 33045391 || rid == 4239289 || rid == 5267597 || rid == 34474764 || rid == 35016489 + // || rid == 36002255 || rid == 37811694 || rid == 46805824) + // { + // fprintf(stderr, "rid: %lu, rev1: %lu, uID1: %lu, ref_p1: %lu, self_p1: %lu, rev2: %lu, uID2: %lu, ref_p2: %lu, self_p2: %lu\n", + // rid, rev1, uID1, self_p1, ref_p1, rev2, uID2, ref_p2, self_p2); + // } + /*******************************for debug************************************/ + if(uID1 == uID2) return; + if(ref_p1 < self_p1 || ref_p2 < self_p2) return; + if(occ1 > 1) + { + interpret_pos(idx, &l1[1], &rev_t, &uID_t, &ref_p_t, &self_p_t, &eLen_t, &tLen_t); + if(uID_t != uID1 && uID_t != uID2) return; + } + + if(occ2 > 1) + { + interpret_pos(idx, &l2[1], &rev_t, &uID_t, &ref_p_t, &self_p_t, &eLen_t, &tLen_t); + if(uID_t != uID1 && uID_t != uID2) return; + } + + x->id = rid; + ref_p1 -= self_p1; + if(rev1) ref_p1 = idx->ug->u.a[uID1].len - 1 - ref_p1; + x->s = (rev1<<63) | ((uID1 << (64-idx->uID_bits))>>1) | (ref_p1 & idx->pos_mode); + + ref_p2 -= self_p2; + if(rev2) ref_p2 = idx->ug->u.a[uID2].len - 1 - ref_p2; + x->e = (rev2<<63) | ((uID2 << (64-idx->uID_bits))>>1) | (ref_p2 & idx->pos_mode); +} + +static void worker_for_alignment(void *data, long i, int tid) // callback for kt_for() +{ + stepdat_t *s = (stepdat_t*)data; + s->pos[i].id = s->pos[i].s = s->pos[i].e = (uint64_t)-1; + uint64_t len1 = s->len[i]>>32, len2 = (uint32_t)s->len[i], occ1, occ2; + char *r1 = s->seq[i], *r2 = s->seq[i] + len1; + /*******************************for debug************************************/ + // if(memcmp(r1, R1.r.a + R1.r_Len.a[s->id+i], len1) != 0) + // { + // fprintf(stderr, "haha1\n"); + // } + // if(memcmp(r2, R2.r.a + R2.r_Len.a[s->id+i], len2) != 0) + // { + // fprintf(stderr, "haha2\n"); + // } + /*******************************for debug************************************/ + s->pos_buf[tid].a.n = 0; + get_alignment(r1, len1, s->idx->k, &s->pos_buf[tid], s->idx, 0, s->id+i); + occ1 = s->pos_buf[tid].a.n; + if(occ1 == 0) return; + get_alignment(r2, len2, s->idx->k, &s->pos_buf[tid], s->idx, occ1, s->id+i); + occ2 = s->pos_buf[tid].a.n - occ1; + if(occ2 == 0) return; + + set_pe_pos((ha_ug_index*)s->idx, s->pos_buf[tid].a.a, occ1, s->pos_buf[tid].a.a + occ1, occ2, &(s->pos[i]), s->id+i); + + /*******************************for debug************************************/ + // if(memcmp(r1, R1.r.a + R1.r_Len.a[s->id+i], len1) != 0) + // { + // fprintf(stderr, "haha1\n"); + // } + // if(memcmp(r2, R2.r.a + R2.r_Len.a[s->id+i], len2) != 0) + // { + // fprintf(stderr, "haha2\n"); + // } + // uint64_t j, rev, uID, ref_p, self_p, eLen, tLen; + // char dir[2] = {'+', '-'}; + // fprintf(stderr, "(R1) %.*s\n", (int)(R1.name_Len.a[s->id + i + 1] - R1.name_Len.a[s->id+i]), + // R1.name.a + R1.name_Len.a[s->id+i]); + // for (j = 0; j < occ1; j++) + // { + // interpret_pos(s->idx, &s->pos_buf[tid].a.a[j], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + // fprintf(stderr, "utg%.6lu\t%c\t%lu\t%lu-%lu\n", uID+1, dir[rev], ref_p, self_p + 1 - tLen, self_p); + // } + // fprintf(stderr, "(R2) %.*s\n", (int)(R2.name_Len.a[s->id + i + 1] - R2.name_Len.a[s->id+i]), + // R2.name.a + R2.name_Len.a[s->id+i]); + // for (j = 0; j < occ2; j++) + // { + // interpret_pos(s->idx, &s->pos_buf[tid].a.a[j+occ1], &rev, &uID, &ref_p, &self_p, &eLen, &tLen); + // fprintf(stderr, "utg%.6lu\t%c\t%lu\t%lu-%lu\n", uID+1, dir[rev], ref_p, self_p + 1 - tLen, self_p); + // } + // fprintf(stderr, "\n"); + /*******************************for debug************************************/ +} + +static void *worker_pipeline(void *data, int step, void *in) // callback for kt_pipeline() +{ + sldat_t *p = (sldat_t*)data; + ///uint64_t total_base = 0, total_pair = 0; + if (step == 0) { // step 1: read a block of sequences + int ret1, ret2; + uint64_t l1, l2; + stepdat_t *s; + CALLOC(s, 1); + s->idx = p->idx; s->id = p->total_pair; + while (((ret1 = kseq_read(p->ks1)) >= 0)&&((ret2 = kseq_read(p->ks2)) >= 0)) + { + if (p->ks1->seq.l < p->idx->k || p->ks2->seq.l < p->idx->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); + } + + l1 = p->ks1->seq.l; l2 = p->ks2->seq.l; + MALLOC(s->seq[s->n], l1+l2); + s->sum_len += l1+l2; + memcpy(s->seq[s->n], p->ks1->seq.s, l1); + memcpy(s->seq[s->n]+l1, p->ks2->seq.s, l2); + s->len[s->n++] = (uint64_t)(l1<<32)|(uint64_t)l2; + + 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 + stepdat_t *s = (stepdat_t*)in; + CALLOC(s->pos_buf, p->n_thread); + MALLOC(s->pos, s->n); + int i; + kt_for(p->n_thread, worker_for_alignment, s, s->n); + for (i = 0; i < s->n; ++i) { + free(s->seq[i]); + p->total_base += (s->len[i]>>32) + (uint32_t)s->len[i]; + } + + free(s->seq); free(s->len); + for (i = 0; i < (int)p->n_thread; ++i) { + free(s->pos_buf[i].a.a); + } + free(s->pos_buf); + return s; + } + else if (step == 2) { // step 3: dump + stepdat_t *s = (stepdat_t*)in; + int i; + for (i = 0; i < s->n; ++i) { + if(s->pos[i].s == (uint64_t)-1) continue; + kv_push(pe_hit, p->hits.a, s->pos[i]); + } + free(s->pos); + free(s); + } + return 0; +} + + +void load_reads(reads_t* x, const char *fn) +{ + kv_init(x->name); + kv_init(x->name_Len); + kv_init(x->r); + kv_init(x->r_Len); + gzFile fp; + kseq_t *ks; + int ret; + uint64_t name_tot, base_total; + + name_tot = base_total = 0; + if ((fp = gzopen(fn, "r")) == 0) return; + ks = kseq_init(fp); + while (((ret = kseq_read(ks)) >= 0)) + { + kv_push(uint64_t, x->name_Len, name_tot); + kv_resize(char, x->name, name_tot + ks->name.l); + memcpy(x->name.a + name_tot, ks->name.s, ks->name.l); + name_tot += ks->name.l; + + kv_push(uint64_t, x->r_Len, base_total); + kv_resize(char, x->r, base_total + ks->seq.l); + memcpy(x->r.a + base_total, ks->seq.s, ks->seq.l); + base_total += ks->seq.l; + } + + kv_push(uint64_t, x->name_Len, name_tot); + kv_push(uint64_t, x->r_Len, base_total); + + kseq_destroy(ks); + gzclose(fp); +} + + +void test_reads(reads_t* x, const char *fn) +{ + gzFile fp; + kseq_t *ks; + int ret, i = 0; + + if ((fp = gzopen(fn, "r")) == 0) return; + ks = kseq_init(fp); + while (((ret = kseq_read(ks)) >= 0)) + { + if(memcmp(ks->name.s, x->name.a + x->name_Len.a[i], ks->name.l) != 0) + { + fprintf(stderr, "ERROR222: i: %d, len: %lu\n", i, x->name_Len.a[i]); + } + i++; + } + + kseq_destroy(ks); + gzclose(fp); +} + +void destory_reads(reads_t* x) +{ + kv_destroy(x->name); + kv_destroy(x->name_Len); + kv_destroy(x->r); + kv_destroy(x->r_Len); +} + +void print_hits(ha_ug_index* idx, kvec_pe_hit* hits, const char *fn) +{ + uint64_t k, shif = 64 - idx->uID_bits; + reads_t r1; + load_reads(&r1, fn); + char dir[2] = {'+', '-'}; + for (k = 0; k < hits->a.n; ++k) + { + fprintf(stderr, "%.*s\t%c\tutg%.6d\t%lu\t%c\tutg%.6d\t%lu\ti:%lu\n", + (int)(r1.name_Len.a[hits->a.a[k].id + 1] - r1.name_Len.a[hits->a.a[k].id]), + r1.name.a + r1.name_Len.a[hits->a.a[k].id], + dir[hits->a.a[k].s>>63], (int)((hits->a.a[k].s<<1)>>shif)+1, hits->a.a[k].s&idx->pos_mode, + dir[hits->a.a[k].e>>63], (int)((hits->a.a[k].e<<1)>>shif)+1, hits->a.a[k].e&idx->pos_mode, + hits->a.a[k].id); + } + destory_reads(&r1); +} + +void dedup_hits(kvec_pe_hit* hits) +{ + double index_time = yak_realtime(); + uint64_t k, l, m = 0, cur; + radix_sort_pe_hit_an1(hits->a.a, hits->a.a + hits->a.n); + for (k = 1, l = 0; k <= hits->a.n; ++k) + { + if (k == hits->a.n || (hits->a.a[k].s<<1) != (hits->a.a[l].s<<1)) + { + if (k - l > 1) radix_sort_pe_hit_an2(hits->a.a + l, hits->a.a + k); + cur = (uint64_t)-1; + while (l < k) + { + if(hits->a.a[l].e != cur) + { + cur = hits->a.a[l].e; + hits->a.a[m++] = hits->a.a[l]; + } + l++; + } + l = k; + } + } + hits->a.n = m; + fprintf(stderr, "[M::%s::%.3f] ==> Dedup\n", __func__, yak_realtime()-index_time); +} + +void sort_hits(kvec_pe_hit* hits) +{ + double index_time = yak_realtime(); + uint64_t k, l; + radix_sort_pe_hit_an1(hits->a.a, hits->a.a + hits->a.n); + for (k = 1, l = 0; k <= hits->a.n; ++k) + { + if (k == hits->a.n || (hits->a.a[k].s<<1) != (hits->a.a[l].s<<1)) + { + if (k - l > 1) radix_sort_pe_hit_an2(hits->a.a + l, hits->a.a + k); + l = k; + } + } + fprintf(stderr, "[M::%s::%.3f] ==> Sort\n", __func__, yak_realtime()-index_time); +} + +void destory_bubbles(bubble_type* bub) +{ + if(bub->index) free(bub->index); + kv_destroy(bub->list); + kv_destroy(bub->num); +} + +inline void get_bubbles(bubble_type* bub, uint64_t id, uint32_t* beg, uint32_t* sink, uint32_t** a, uint32_t* n) +{ + (*a) = bub->list.a + bub->num.a[id] + 2; + (*n) = bub->num.a[id+1] - bub->num.a[id] - 2; + (*beg) = bub->list.a[bub->num.a[id]]; + (*sink) = bub->list.a[bub->num.a[id] + 1]; +} + +void identify_bubbles(ma_ug_t* ug, bubble_type* bub) +{ + asg_cleanup(ug->g); + if (!ug->g->is_symm) asg_symm(ug->g); + memset(bub, 0, sizeof(bubble_type)); + uint32_t v, n_vtx = ug->g->n_seq * 2, tLen, i, mode = (((uint32_t)-1)<<2); + bub->ug = ug; + CALLOC(bub->index, n_vtx); + for (i = 0; i < ug->g->n_seq; i++) + { + if(ug->g->seq[i].c == 1) + { + ug->g->seq[i].c = 0; + bub->index[i] = 4; + } + } + kv_init(bub->list); kv_init(bub->num); + buf_t b; memset(&b, 0, sizeof(buf_t)); b.a = (binfo_t*)calloc(n_vtx, sizeof(binfo_t)); + for (i = 0, tLen = 1; i < ug->u.n; i++) tLen += ug->u.a[i].len; + + for (v = 0; v < n_vtx; ++v) + { + if(ug->g->seq[v>>1].del) continue; + if(asg_arc_n(ug->g, v) < 2) continue; + if((bub->index[v]&(uint32_t)3) != 0) continue; + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0)) + { + //beg is v, end is b.S.a[0] + //note b.b include end, does not include beg + for (i = 0; i < b.b.n; i++) + { + if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue; + bub->index[b.b.a[i]] &= mode; bub->index[b.b.a[i]] += 1; + bub->index[b.b.a[i]^1] &= mode; bub->index[b.b.a[i]^1] += 1; + } + bub->index[v] &= mode; bub->index[v] += 2; + bub->index[b.S.a[0]^1] &= mode;; bub->index[b.S.a[0]^1] += 3; + } + } + + + for (v = 0; v < n_vtx; ++v) + { + if((bub->index[v]&(uint32_t)3) !=2) continue; + kv_push(uint32_t, bub->num, bub->list.n); + if(asg_bub_pop1_primary_trio(ug->g, NULL, v, tLen, &b, (uint32_t)-1, (uint32_t)-1, 0)) + { + //beg is v, end is b.S.a[0] + kv_push(uint32_t, bub->list, v); + kv_push(uint32_t, bub->list, b.S.a[0]^1); + + //note b.b include end, does not include beg + for (i = 0; i < b.b.n; i++) + { + if(b.b.a[i]==v || b.b.a[i]==b.S.a[0]) continue; + kv_push(uint32_t, bub->list, b.b.a[i]); + } + } + } + + kv_push(uint32_t, bub->num, bub->list.n); + free(b.a); free(b.S.a); free(b.T.a); free(b.b.a); free(b.e.a); + + for (i = 0; i < ug->g->n_seq; i++) + { + if((bub->index[i]>>2) == 0) + { + bub->index[i] = (uint32_t)-1; + } + else + { + bub->index[i] = bub->num.n; + } + } + + uint32_t beg, sink, n, *a; + for (i = 0; i < bub->num.n-1; i++) + { + get_bubbles(bub, i, &beg, &sink, &a, &n); + for (v = 0; v < n; v++) + { + bub->index[(a[v]>>1)] = i; + } + } + + ///free(bub->index); bub->index = NULL; +} + + + +void print_bubbles(ma_ug_t* ug, bubble_type* bub, kvec_pe_hit* hits, hc_links* link, ha_ug_index* idx) +{ + uint64_t tLen, t_utg, i, k; + uint32_t beg, sink, n, *a; + for (i = 0, tLen = 0; i < bub->ug->u.n; i++) tLen += bub->ug->u.a[i].len; + fprintf(stderr, "[M::%s] # unitigs: %lu, # bases: %lu\n", __func__, bub->ug->u.n, tLen); + for (i = 0, tLen = 0, t_utg = 0; i < bub->num.n-1; i++) + { + get_bubbles(bub, i, &beg, &sink, &a, &n); + t_utg += n; + for (k = 0; k < n; k++) + { + tLen +=bub->ug->u.a[(a[k]>>1)].len; + } + } + fprintf(stderr, "[M::%s] # bubbles: %lu, # unitigs: %lu, # bases: %lu\n", __func__, + Get_bub_num(*bub), t_utg, tLen); + + + for (i = 0, tLen = 0, t_utg = 0; i < ug->g->n_seq; i++) + { + if(bub->index[i] < bub->num.n) + { + t_utg++; + tLen +=bub->ug->u.a[i].len; + } + } + fprintf(stderr, "[M::%s] # bubbles: %lu, # unitigs: %lu, # bases: %lu\n", __func__, + Get_bub_num(*bub), t_utg, tLen); + + for (i = 0, tLen = 0, t_utg = 0; i < ug->g->n_seq; i++) + { + if(bub->index[i] == bub->num.n) + { + t_utg++; + tLen +=bub->ug->u.a[i].len; + } + } + fprintf(stderr, "[M::%s] # het unitigs: %lu, # het bases: %lu\n", __func__, t_utg, tLen); + + uint8_t* flag; CALLOC(flag, ug->g->n_seq); + uint64_t s_uid, e_uid, shif = 64 - idx->uID_bits; + if(hits) + { + for (k = 0; k < hits->a.n; ++k) + { + s_uid = ((hits->a.a[k].s<<1)>>shif); + e_uid = ((hits->a.a[k].e<<1)>>shif); + if(bub->index[s_uid] == (uint32_t)-1 || bub->index[e_uid] == (uint32_t)-1) continue; + if(bub->index[s_uid] < bub->num.n && bub->index[e_uid] < bub->num.n) + { + flag[s_uid] |= 1; + flag[e_uid] |= 1; + continue; + } + if(bub->index[s_uid] == bub->num.n && bub->index[e_uid] == bub->num.n) + { + flag[s_uid] |= 4; + flag[e_uid] |= 4; + continue; + } + if(bub->index[s_uid] < bub->num.n) flag[s_uid] |= 2, flag[e_uid] |= 2; + if(bub->index[e_uid] < bub->num.n) flag[e_uid] |= 2, flag[s_uid] |= 2; + } + } + else if(link) + { + for (i = 0; i < link->a.n; i++) + { + for (k = 0; k < link->a.a[i].e.n; ++k) + { + if(link->a.a[i].e.a[k].del) continue; + s_uid = i; + e_uid = link->a.a[i].e.a[k].uID; + if(bub->index[s_uid] == (uint32_t)-1 || bub->index[e_uid] == (uint32_t)-1) continue; + if(bub->index[s_uid] < bub->num.n && bub->index[e_uid] < bub->num.n) + { + flag[s_uid] |= 1; + flag[e_uid] |= 1; + continue; + } + if(bub->index[s_uid] == bub->num.n && bub->index[e_uid] == bub->num.n) + { + flag[s_uid] |= 4; + flag[e_uid] |= 4; + continue; + } + if(bub->index[s_uid] < bub->num.n) flag[s_uid] |= 2, flag[e_uid] |= 2; + if(bub->index[e_uid] < bub->num.n) flag[e_uid] |= 2, flag[s_uid] |= 2; + } + } + } + + + + + for (i = 0, tLen = 0, t_utg = 0; i < ug->g->n_seq; i++) + { + if(flag[i] & (uint32_t)1) + { + t_utg++; + tLen +=bub->ug->u.a[i].len; + } + } + fprintf(stderr, "[M::%s] # bubble-chained unitigs: %lu, # bubble-chained bases: %lu\n", + __func__, t_utg, tLen); + + for (i = 0, tLen = 0, t_utg = 0; i < ug->g->n_seq; i++) + { + if((flag[i] & (uint32_t)1) || (flag[i] & (uint32_t)2)) + { + t_utg++; + tLen +=bub->ug->u.a[i].len; + } + } + fprintf(stderr, "[M::%s] # (bubble && het)-chained unitigs: %lu, # (bubble && het)-chained bases: %lu\n", + __func__, t_utg, tLen); + + + for (i = 0, tLen = 0, t_utg = 0; i < ug->g->n_seq; i++) + { + if((flag[i] & (uint32_t)1) || (flag[i] & (uint32_t)2) || (flag[i] & (uint32_t)4)) + { + t_utg++; + tLen +=bub->ug->u.a[i].len; + } + } + fprintf(stderr, "[M::%s] # (bubble || het)-chained unitigs: %lu, # (bubble || het)-chained bases: %lu\n", + __func__, t_utg, tLen); + free(flag); + + + // fprintf(stderr, "************bubble utgs************\n"); + // for (i = 0, tLen = 0, t_utg = 0; i < bub->num.n-1; i++) + // { + // get_bubbles(bub, i, &beg, &sink, &a, &n); + // t_utg += n; + // fprintf(stderr, "(%lu)\tbeg:utg%.6u\tsink:utg%.6u\n", i, (beg>>1)+1, (sink>>1)+1); + // for (k = 0; k < n; k++) + // { + // tLen +=bub->ug->u.a[(a[k]>>1)].len; + // fprintf(stderr, "utg%.6u,", (a[k]>>1)+1); + // } + // fprintf(stderr, "\n"); + // } + + // fprintf(stderr, "************het utgs************\n"); + // for (i = 0; i < ug->g->n_seq; i++) + // { + // if(bub->index[i] == bub->num.n) fprintf(stderr, "utg%.6lu\n", i+1); + // } + // fprintf(stderr, "************het utgs************\n"); +} + + + + +void init_hc_links(hc_links* link, uint64_t ug_num) +{ + kv_init(link->a); + kv_resize(hc_linkeage, link->a, ug_num); + link->a.n = ug_num; + uint64_t i; + for (i = 0; i < link->a.n; i++) + { + kv_init(link->a.a[i].e); + kv_init(link->a.a[i].f); + } +} + +void destory_hc_links(hc_links* link) +{ + uint64_t i; + for (i = 0; i < link->a.n; i++) + { + kv_destroy(link->a.a[i].e); + kv_destroy(link->a.a[i].f); + } + kv_destroy(link->a); +} + +void push_hc_edge(hc_linkeage* x, uint64_t uID, int weight, int dir) +{ + uint64_t k, n; + hc_edge* a = NULL; + hc_edge* p = NULL; + if(dir == 0) + { + a = x->e.a; + n = x->e.n; + } + else + { + a = x->f.a; + n = x->f.n; + } + + for (k = 0; k < n; k++) + { + if(a[k].del) continue; + if(a[k].uID == uID) + { + a[k].weight += weight; + return; + } + } + + if(dir == 0) + { + kv_pushp(hc_edge, x->e, &p); + } + else + { + kv_pushp(hc_edge, x->f, &p); + } + + p->del = p->enzyme = 0; + p->uID = uID; + p->weight = weight; +} + + +void collect_hc_links(const ha_ug_index* idx, kvec_pe_hit* hits, hc_links* link) +{ + uint64_t k, shif = 64 - idx->uID_bits, beg, end; + for (k = 0; k < hits->a.n; ++k) + { + beg = ((hits->a.a[k].s<<1)>>shif); + end = ((hits->a.a[k].e<<1)>>shif); + push_hc_edge(&(link->a.a[beg]), end, 1, 0); + push_hc_edge(&(link->a.a[end]), beg, 1, 0); + } +} + +void dfs_bubble(asg_t *g, kvec_t_u32_warp* stack, kvec_t_u32_warp* result, +uint32_t v, uint32_t beg, uint32_t sink) +{ + asg_arc_t *acur = NULL; + uint32_t cur, ncur, i; + stack->a.n = result->a.n = 0; + v = v << 1; + kv_push(uint32_t, stack->a, v); + while (stack->a.n > 0) + { + stack->a.n--; + cur = stack->a.a[stack->a.n]; + kv_push(uint32_t, result->a, cur>>1); + ncur = asg_arc_n(g, cur); + acur = asg_arc_a(g, cur); + for (i = 0; i < ncur; i++) + { + if(acur[i].del) continue; + if((acur[i].v>>1) == beg || (acur[i].v>>1) == sink) continue; + kv_push(uint32_t, stack->a, acur[i].v); + } + } + + + v = v + 1; + kv_push(uint32_t, stack->a, v); + while (stack->a.n > 0) + { + stack->a.n--; + cur = stack->a.a[stack->a.n]; + kv_push(uint32_t, result->a, cur>>1); + ncur = asg_arc_n(g, cur); + acur = asg_arc_a(g, cur); + for (i = 0; i < ncur; i++) + { + if(acur[i].del) continue; + if((acur[i].v>>1) == beg || (acur[i].v>>1) == sink) continue; + kv_push(uint32_t, stack->a, acur[i].v); + } + } +} + +void set_reverse_links(uint32_t* bub, uint32_t n, kvec_t_u32_warp* reach, uint32_t root, hc_links* link) +{ + uint64_t i, k; + uint32_t v; + for (i = 0; i < n; i++) + { + v = bub[i]>>1; + for (k = 0; k < reach->a.n; k++) + { + if(v == reach->a.a[k]) break; + } + + if(k == reach->a.n) + { + push_hc_edge(&(link->a.a[root]), v, 1, 1); + push_hc_edge(&(link->a.a[v]), root, 1, 1); + } + } + +} + +void collect_hc_reverse_links(hc_links* link, ma_ug_t* ug, bubble_type* bub) +{ + uint64_t i, k; + uint32_t beg, sink, n, v, *a = NULL; + kvec_t_u32_warp stack, result; + kv_init(stack.a); kv_init(result.a); + for (i = 0; i < bub->num.n-1; i++) + { + get_bubbles(bub, i, &beg, &sink, &a, &n); + if(n == 2) + { + push_hc_edge(&(link->a.a[a[0]>>1]), a[1]>>1, 1, 1); + push_hc_edge(&(link->a.a[a[1]>>1]), a[0]>>1, 1, 1); + continue; + } + beg = beg>>1; sink = sink>>1; + for (k = 0; k < n; k++) + { + v = a[k]>>1; + dfs_bubble(ug->g, &stack, &result, v, beg, sink); + set_reverse_links(a, n, &result, v, link); + } + } + kv_destroy(stack.a); kv_destroy(result.a); +} + +void write_hc_links(hc_links* link, const char *fn) +{ + uint64_t k; + char *buf = (char*)calloc(strlen(fn) + 25, 1); + sprintf(buf, "%s.hic.lk", fn); + FILE* fp = fopen(buf, "w"); + + fwrite(&link->a.n, sizeof(link->a.n), 1, fp); + for (k = 0; k < link->a.n; k++) + { + fwrite(&link->a.a[k].e.n, sizeof(link->a.a[k].e.n), 1, fp); + fwrite(link->a.a[k].e.a, sizeof(hc_edge), link->a.a[k].e.n, fp); + + fwrite(&link->a.a[k].f.n, sizeof(link->a.a[k].f.n), 1, fp); + fwrite(link->a.a[k].f.a, sizeof(hc_edge), link->a.a[k].f.n, fp); + } + + fclose(fp); + free(buf); +} + +int load_hc_links(hc_links* link, const char *fn) +{ + uint64_t k, flag = 0; + char *buf = (char*)calloc(strlen(fn) + 25, 1); + sprintf(buf, "%s.hic.lk", fn); + + FILE* fp = NULL; + fp = fopen(buf, "r"); + if(!fp) return 0; + + kv_init(link->a); + flag += fread(&link->a.n, sizeof(link->a.n), 1, fp); + link->a.m = link->a.n; CALLOC(link->a.a, link->a.n); + for (k = 0; k < link->a.n; k++) + { + flag += fread(&link->a.a[k].e.n, sizeof(link->a.a[k].e.n), 1, fp); + link->a.a[k].e.m = link->a.a[k].e.n; MALLOC(link->a.a[k].e.a, link->a.a[k].e.n); + flag += fread(link->a.a[k].e.a, sizeof(hc_edge), link->a.a[k].e.n, fp); + + flag += fread(&link->a.a[k].f.n, sizeof(link->a.a[k].f.n), 1, fp); + link->a.a[k].f.m = link->a.a[k].f.n; MALLOC(link->a.a[k].f.a, link->a.a[k].f.n); + flag += fread(link->a.a[k].f.a, sizeof(hc_edge), link->a.a[k].f.n, fp); + } + + + fclose(fp); + free(buf); + return 1; +} + +void print_hc_links(hc_links* link) +{ + uint64_t i, k; + for (i = 0; i < link->a.n; ++i) + { + for (k = 0; k < link->a.a[i].e.n; k++) + { + if(link->a.a[i].e.a[k].del) continue; + fprintf(stderr, "utg%.6d\tutg%.6d\t%d\t+\n", + (int)(i+1), (int)(link->a.a[i].e.a[k].uID+1), link->a.a[i].e.a[k].weight); + } + } + + + for (i = 0; i < link->a.n; ++i) + { + for (k = 0; k < link->a.a[i].f.n; k++) + { + if(link->a.a[i].f.a[k].del) continue; + fprintf(stderr, "utg%.6d\tutg%.6d\t%d\t-", + (int)(i+1), (int)(link->a.a[i].f.a[k].uID+1), link->a.a[i].f.a[k].weight); + if(link->a.a[i].f.a[k].weight > 2) fprintf(stderr,"\tcomplex"); + fprintf(stderr,"\n"); + } + } +} + +void init_min_cut_t(min_cut_t* x, hc_links* link, const bubble_type* bub, const ma_ug_t *ug) +{ + uint64_t utg_num = link->a.n, i, k; + x->n = utg_num; + + kv_malloc(x->rGraphSet, utg_num); x->rGraphSet.n = utg_num; + kv_malloc(x->rGraphVis, utg_num); x->rGraphVis.n = utg_num; + kv_malloc(x->utgVis, utg_num); x->utgVis.n = utg_num; + kv_malloc(x->order, utg_num); x->order.n = utg_num; + kv_malloc(x->parent, utg_num); x->parent.n = utg_num; + ///uresolved BUGs, if use kv_resize segfault; if use kv_malloc, work????? + kv_malloc(x->rGraph, utg_num); x->rGraph.n = utg_num; + // kv_init(x->rGraph); kv_resize(hc_edge_warp, x->rGraph, utg_num); x->rGraph.n = utg_num; + + //must be utg_num + 2 since we may need to add fake nodes + for (i = 1; (uint64_t)(1<uID_mode = ((uint64_t)-1) >> (64-i); + x->uID_shift = i; + for (i = 0; i < utg_num; i++) + { + ///x->order.a[i] = link->a.a[i].f.n; + x->order.a[i] = ug->u.a[i].len; + x->order.a[i] <<= x->uID_shift; + x->order.a[i] |= (uint64_t)(i & x->uID_mode); + + x->rGraphSet.a[i] = 0; + x->rGraphVis.a[i] = 0; + x->utgVis.a[i] = 0; + x->parent.a[i] = (uint64_t)-1; + + ///uresolved BUGs, if use kv_resize segfault; if use kv_malloc, work????? + // kv_init(x->rGraph.a[i]); kv_resize(hc_edge, x->rGraph.a[i], link->a.a[i].e.n); + kv_malloc(x->rGraph.a[i], link->a.a[i].e.n); + x->rGraph.a[i].n = link->a.a[i].e.n; + + if(x->rGraph.a[i].n) + { + for (k = 0; k < x->rGraph.a[i].n; k++) + { + ///kv_push(hc_edge, x->rGraph.a[i], link->a.a[i].e.a[k]); + x->rGraph.a[i].a[k] = link->a.a[i].e.a[k]; + if((bub->index[x->rGraph.a[i].a[k].uID] > bub->num.n) || + (bub->index[i] > bub->num.n)) + { + x->rGraph.a[i].a[k].del = 1; + } + } + } + } + + x->q = kdq_init(uint64_t); + radix_sort_hc64(x->order.a, x->order.a + x->order.n); + + ///fprintf(stderr, "[M::%s]\n", __func__); + ///exit(0); +} + +void destory_min_cut_t(min_cut_t* x) +{ + kv_destroy(x->order); + kv_destroy(x->parent); + kv_destroy(x->rGraphSet); + kv_destroy(x->rGraphVis); + kv_destroy(x->utgVis); + uint64_t i; + for (i = 0; i < x->rGraph.m; i++) + { + kv_destroy(x->rGraph.a[i]); + } + kv_destroy(x->rGraph); + kdq_destroy(uint64_t, x->q); +} + +void reset_min_cut_t(min_cut_t* x, hc_links* link) +{ + ///no need to reset parent[] and q + uint64_t i, j; + ///important to have this line + x->parent.n = x->order.n = x->rGraph.n = x->rGraphVis.n = x->rGraphSet.n = link->a.n; + kdq_clear(x->q); + + for (i = 0; i < x->rGraphSet.n; i++) + { + x->rGraphVis.a[i] = 0; + ///important to have this line + x->rGraph.a[i].n = link->a.a[i].e.n; + + if(x->rGraphSet.a[i] == 0) continue; + for (j = 0; j < x->rGraph.a[i].n; j++) + { + x->rGraph.a[i].a[j].weight = link->a.a[i].e.a[j].weight; + } + x->rGraphSet.a[i] = 0; + } +} + +uint64_t add_mul_convex(min_cut_t* x, uint64_t* a, uint64_t n) +{ + if(n == 0) return (uint64_t)-1; + if(n == 1) return a[0]; + kv_push(uint8_t, x->rGraphSet, 0); + kv_push(uint8_t, x->rGraphVis, 0); + kv_push(uint64_t, x->parent, 0); + kv_resize(hc_edge_warp, x->rGraph, x->rGraph.n+1); + kv_init(x->rGraph.a[x->rGraph.n]); + uint64_t i, k; + hc_edge t; + for (i = 0; i < n; i++) + { + t.uID = a[i]; t.del = t.enzyme = t.weight = 0; + for (k = 0; k < x->rGraph.a[a[i]].n; k++) + { + if(x->rGraph.a[a[i]].a[k].del) continue; + t.weight += x->rGraph.a[a[i]].a[k].weight; + } + kv_push(hc_edge, x->rGraph.a[x->rGraph.n], t); + t.uID = x->rGraph.n; + kv_push(hc_edge, x->rGraph.a[a[i]], t); + } + + x->rGraph.n++; + return x->rGraph.n - 1; +} + +void get_s_t(min_cut_t* x, hc_links* link, uint64_t uID, uint64_t* src, uint64_t* dest, kvec_t_u64_warp* buff) +{ + buff->a.n = 0; (*src) = (*dest) = (uint64_t)-1; + if(link->a.a[uID].f.n == 0) return; + (*src) = uID; + if(link->a.a[uID].f.n == 1) + { + (*dest) = link->a.a[uID].f.a[0].uID; + return; + } + + uint64_t i; + for (i = 0; i < link->a.a[uID].f.n; i++) + { + if(link->a.a[uID].f.a[i].del) continue; + kv_push(uint64_t, buff->a, link->a.a[uID].f.a[i].uID); + } + (*dest) = add_mul_convex(x, buff->a.a, buff->a.n); +} + + +uint64_t bfs_flow(uint64_t src, uint64_t dest, min_cut_t* x, kvec_t_u64_warp* buff) +{ + uint64_t *p = NULL, v, u, i; + if(dest != (uint64_t)-1) memset(x->rGraphVis.a, 0, x->rGraphVis.n); + + kdq_push(uint64_t, x->q, src); + if(buff) kv_push(uint64_t, buff->a, src); + + x->rGraphVis.a[src] = 1; + x->parent.a[src] = (uint64_t)-1; + + while (1) + { + p = kdq_shift(uint64_t, x->q); + if(!p) break; + v = *p; + if(v == dest) return 1; + for (i = 0; i < x->rGraph.a[v].n; i++) + { + if(x->rGraph.a[v].a[i].del) continue; + if(x->rGraph.a[v].a[i].weight == 0) continue; + u = x->rGraph.a[v].a[i].uID; + if(x->rGraphVis.a[u]) continue; + + ///x->parent.a[u] = v; + x->parent.a[u] = x->rGraph.a[v].a[i].weight; + x->parent.a[u] <<= x->uID_shift; + x->parent.a[u] |= v; + + kdq_push(uint64_t, x->q, u); + if(buff) kv_push(uint64_t, buff->a, u); + ///set u or v to be 1? doesn't matter + x->rGraphVis.a[u] = 1; + } + } + + return 0; +} + +hc_edge* get_rGraph_edge(min_cut_t* x, uint64_t src, uint64_t dest) +{ + if(src >= x->rGraph.n) return NULL; + uint64_t i; + for (i = 0; i < x->rGraph.a[src].n; i++) + { + if(x->rGraph.a[src].a[i].uID == dest) return &(x->rGraph.a[src].a[i]); + } + return NULL; +} + +uint64_t maxFlow(uint64_t src, uint64_t dest, min_cut_t* x) +{ + uint64_t flow = 0, max_flow = 0, v, u; + hc_edge *p; + + while (bfs_flow(src, dest, x, NULL)) + { + kdq_clear(x->q); + flow = (uint64_t)-1; + for (v = dest; v != src; v = x->parent.a[v]&x->uID_mode) + { + flow = MIN(flow, (x->parent.a[v]>>x->uID_shift)); + } + + fprintf(stderr, "***********flow: %lu*********\n", flow); + + for (v = dest; v != src; v = x->parent.a[v]&x->uID_mode) + { + u = x->parent.a[v]&x->uID_mode; + p = get_rGraph_edge(x, u, v); + + fprintf(stderr, "utg%.6lul (%d)\n", u+1, p->weight); + + + p->weight -= flow; + p = get_rGraph_edge(x, v, u); + p->weight += flow; + x->rGraphSet.a[u] = x->rGraphSet.a[v] = 1; + } + + max_flow += flow; + } + + return max_flow; +} + +void print_src_dest(uint64_t src, min_cut_t* x, const char* command) +{ + uint64_t i; + fprintf(stderr, "********************\n%s\n", command); + if(src >= x->n) + { + for (i = 0; i < x->rGraph.a[src].n; i++) + { + if(x->rGraph.a[src].a[i].del) continue; + fprintf(stderr, "utg%.6ul\n", x->rGraph.a[src].a[i].uID + 1); + } + + } + else + { + fprintf(stderr, "utg%.6lul\n", src+1); + } + fprintf(stderr, "!!!!!!!!!!!!!!!!!!!!\n"); + +} + +void graph_cut(uint64_t src, uint64_t dest, min_cut_t* x) +{ + if(maxFlow(src, dest, x)) + { + ///in the last time bfs of maxFlow, rGraphVis has already been set + uint64_t i, j, v, u; + hc_edge *p; + /*******************************for debug************************************/ + print_src_dest(src, x, "src utg:"); + print_src_dest(dest, x, "dest utg:"); + /*******************************for debug************************************/ + for (i = 0; i < x->rGraphVis.n; i++) + { + if(x->rGraphVis.a[i] == 0) continue; + v = i; + for (j = 0; j < x->rGraph.a[i].n; j++) + { + if(x->rGraph.a[i].a[j].del) continue; + u = x->rGraph.a[i].a[j].uID; + if(x->rGraphVis.a[u]) continue; + /*******************************for debug************************************/ + fprintf(stderr, "utg%.6lul\tutg%.6lul\t%d\n", v+1, u+1, x->rGraph.a[i].a[j].weight); + /*******************************for debug************************************/ + ///delete + x->rGraph.a[i].a[j].del = 1; + ///delete + p = get_rGraph_edge(x, u, v); + p->del = 1; + } + + } + } +} + +void check_connective(min_cut_t* x, hc_links* link) +{ + double index_time = yak_realtime(); + kvec_t_u64_warp buff; + kv_init(buff.a); + uint64_t i, k, uID; + for (i = 0; i < x->n; i++) + { + uID = x->order.a[i] & x->uID_mode; + if(link->a.a[uID].f.n == 0) continue; + for (k = 0; k < link->a.a[uID].f.n; k++) + { + if(link->a.a[uID].f.a[k].del) continue; + if(x->utgVis.a[link->a.a[uID].f.a[k].uID] == 0) break; + } + if(k == link->a.a[uID].f.n) continue; + reset_min_cut_t(x, link); + get_s_t(x, link, uID, &(x->src), &(x->dest), &buff); + + bfs_flow(x->src, x->dest, x, NULL); + + x->utgVis.a[uID] = 1; + } + + //reset x.utgVis + memset(x->utgVis.a, 0, x->utgVis.n); + kv_destroy(buff.a); + fprintf(stderr, "[M::%s::%.3f] \n", __func__, yak_realtime()-index_time); +} + +void get_Connected_Components(min_cut_t* x) +{ + double index_time = yak_realtime(); + uint64_t i, j, k = 0, uID, e; + kvec_t_u64_warp buff; + kv_init(buff.a); + while (1) + { + for (i = 0; i < x->n; i++) + { + uID = x->order.a[i] & x->uID_mode; + if(x->rGraphVis.a[uID] == 0) break; + } + if(i < x->n) + { + e = buff.a.n = 0; + bfs_flow(uID, (uint64_t)-1, x, &buff); + for (i = 0; i < buff.a.n; i++) + { + for (j = 0; j < x->rGraph.a[buff.a.a[i]].n; j++) + { + if(x->rGraph.a[buff.a.a[i]].a[j].del == 0) e++; + } + } + e >>= 1; + if(buff.a.n > 1) + { + fprintf(stderr, "(%lu) Component: # nodes: %lu, # edges: %lu\n", + k, (uint64_t)buff.a.n, e); + } + k++; + } + else + { + break; + } + } + + kv_destroy(buff.a); + fprintf(stderr, "[M::%s::%.3f] # Connected Components: %lu\n", + __func__, yak_realtime()-index_time, k); +} + +void print_rGraph(min_cut_t* x) +{ + uint64_t i, k; + for (i = 0; i < x->rGraph.n; ++i) + { + for (k = 0; k < x->rGraph.a[i].n; k++) + { + if(x->rGraph.a[i].a[k].del) continue; + fprintf(stderr, "src(utg%.6dl)\tdes(utg%.6dl)\t%d\n", + (int)(i+1), (int)(x->rGraph.a[i].a[k].uID+1), x->rGraph.a[i].a[k].weight); + } + } +} + + +int select_large_node(const ma_ug_t *ug, min_cut_t* x, +uint64_t src, uint64_t dest, uint64_t utg_thres, int weight_thres) +{ + if(src >= ug->u.n || dest >= ug->u.n) return 0; + if(ug->u.a[src].n < utg_thres || ug->u.a[dest].n < utg_thres) return 0; + + uint64_t k; + for (k = 0; k < x->rGraph.a[src].n; k++) + { + if(x->rGraph.a[src].a[k].del) continue; + if(x->rGraph.a[src].a[k].weight >= weight_thres) break; + } + if(k == x->rGraph.a[src].n) return 0; + + src = dest; + for (k = 0; k < x->rGraph.a[src].n; k++) + { + if(x->rGraph.a[src].a[k].del) continue; + if(x->rGraph.a[src].a[k].weight >= weight_thres) break; + } + if(k == x->rGraph.a[src].n) return 0; + + return 1; +} + +void clean_hap(hc_links* link, bubble_type* bub, const ma_ug_t *ug) +{ + min_cut_t x; + kvec_t_u64_warp buff; + kv_init(buff.a); + init_min_cut_t(&x, link, (const bubble_type*)bub, ug); + + // get_Connected_Components(&x); + // check_connective(&x, link); + // print_rGraph(&x); + long long i; + uint64_t k, uID; + ///for (i = x.n - 1; i >= 0; i--) + for (i = 0; (uint64_t)i < x.n; i++) + { + ///fprintf(stderr, "i: %lu\n", i); + uID = x.order.a[i] & x.uID_mode; + ///fprintf(stderr, "uID: %lu, f.n: %lu\n", uID, (uint64_t)link->a.a[uID].f.n); + if(link->a.a[uID].f.n == 0) continue; + for (k = 0; k < link->a.a[uID].f.n; k++) + { + if(link->a.a[uID].f.a[k].del) continue; + if(x.utgVis.a[link->a.a[uID].f.a[k].uID] == 0) break; + } + ///fprintf(stderr, "k: %lu\n", k); + if(k == link->a.a[uID].f.n) continue; + reset_min_cut_t(&x, link); + ///fprintf(stderr, "reset\n"); + get_s_t(&x, link, uID, &(x.src), &(x.dest), &buff); + ///fprintf(stderr, "x.src: %lu, x.dest: %lu\n", x.src, x.dest); + ///Note: should only consider edges betweem bubbles, ignore edges to homo untigs + + /*******************************for debug************************************/ + if(!select_large_node(ug, &x, x.src, x.dest, 10, 10)) continue; + /*******************************for debug************************************/ + + graph_cut(x.src, x.dest, &x); + ///fprintf(stderr, "graph_cut\n"); + x.utgVis.a[uID] = 1; + exit(0); + } + destory_min_cut_t(&x); + kv_destroy(buff.a); +} + +int hic_short_align(const char *fn1, const char *fn2, ha_ug_index* idx) +{ + double index_time = yak_realtime(); + hc_links link; + sldat_t sl; + gzFile fp1, fp2; + if ((fp1 = gzopen(fn1, "r")) == 0) return 0; + if ((fp2 = gzopen(fn2, "r")) == 0) return 0; + sl.ks1 = kseq_init(fp1); + sl.ks2 = kseq_init(fp2); + sl.idx = idx; + sl.chunk_size = 20000000; + sl.n_thread = asm_opt.thread_num; + sl.total_base = sl.total_pair = 0; + idx->max_cnt = 5; + kv_init(sl.hits.a); + fprintf(stderr, "u.n: %d, uID_bits: %lu, pos_bits: %lu\n", (uint32_t)idx->ug->u.n, idx->uID_bits, idx->pos_bits); + + if(!load_hc_links(&link, asm_opt.output_file_name)) + { + /*******************************for debug************************************/ + // load_reads(&R1, fn1); + // test_reads(&R1, fn1); + // load_reads(&R2, fn2); + // test_reads(&R1, fn1); + /*******************************for debug************************************/ + + kt_pipeline(3, worker_pipeline, &sl, 3); + /*******************************for debug************************************/ + // sort_hits(&sl.hits); + // print_hits(idx, &sl.hits, fn1); + /*******************************for debug************************************/ + dedup_hits(&sl.hits); + + init_hc_links(&link, sl.idx->ug->g->n_seq); + collect_hc_links(sl.idx, &sl.hits, &link); + write_hc_links(&link, asm_opt.output_file_name); + } + + /*******************************for debug************************************/ + // destory_reads(&R1); + // destory_reads(&R2); + /*******************************for debug************************************/ + bubble_type bub; + identify_bubbles(idx->ug, &bub); + print_bubbles(idx->ug, &bub, sl.hits.a.n?&sl.hits:NULL, &link, idx); + collect_hc_reverse_links(&link, idx->ug, &bub); + /*******************************for debug************************************/ + ///print_hc_links(&link); + /*******************************for debug************************************/ + clean_hap(&link, &bub, idx->ug); + + kv_destroy(sl.hits.a); + destory_hc_links(&link); + destory_bubbles(&bub); + kseq_destroy(sl.ks1); + kseq_destroy(sl.ks2); + gzclose(fp1); + gzclose(fp2); + + fprintf(stderr, "[M::%s::%.3f] processed %lu pairs; %lu bases\n", __func__, yak_realtime()-index_time, sl.total_pair, sl.total_base); + return 1; } void hic_analysis(ma_ug_t *ug) { + ug_index = NULL; - ///int exist = 0;//load_hc_pt_index(&ug_index, asm_opt.output_file_name); int exist = load_hc_pt_index(&ug_index, asm_opt.output_file_name); if(exist == 0) ug_index = build_unitig_index(ug, asm_opt.hic_mer_length); if(exist == 0) write_hc_pt_index(ug_index, asm_opt.output_file_name); ug_index->ug = ug; - - test_unitig_index(ug_index); + ///test_unitig_index(ug_index, ug); + hic_short_align(asm_opt.hic_files[0], asm_opt.hic_files[1], ug_index); destory_hc_pt_index(ug_index); } \ No newline at end of file diff --git a/htab.cpp b/htab.cpp index dcbf9c9..30f532d 100644 --- a/htab.cpp +++ b/htab.cpp @@ -1002,16 +1002,17 @@ int load_ct_index(void **i_ct_idx, char* file_name) } ha_ct_t** ct_idx = (ha_ct_t**)i_ct_idx; double index_time = 0; + uint64_t flag = 0; int i; ha_ct_t *h = 0; ha_ct1_t *g; CALLOC(h, 1); - fread(&h->k, sizeof(h->k), 1, fp); - fread(&h->pre, sizeof(h->pre), 1, fp); - fread(&h->n_hash, sizeof(h->n_hash), 1, fp); - fread(&h->n_shift, sizeof(h->n_shift), 1, fp); - fread(&h->tot, sizeof(h->tot), 1, fp); + flag += fread(&h->k, sizeof(h->k), 1, fp); + flag += fread(&h->pre, sizeof(h->pre), 1, fp); + flag += fread(&h->n_hash, sizeof(h->n_hash), 1, fp); + flag += fread(&h->n_shift, sizeof(h->n_shift), 1, fp); + flag += fread(&h->tot, sizeof(h->tot), 1, fp); CALLOC(h->h, 1<pre); diff --git a/khashl.h b/khashl.h index 4a8f2dc..4ecf294 100644 --- a/khashl.h +++ b/khashl.h @@ -147,13 +147,14 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26 SCOPE khint_t prefix##_load(HType **h, FILE* fp) { \ (*h) = prefix##_init(); \ khint_t n_buckets; \ - fread(&n_buckets, sizeof(n_buckets), 1, fp); \ - fread(&(*h)->bits, sizeof((*h)->bits), 1, fp); \ - fread(&(*h)->count, sizeof((*h)->count), 1, fp); \ + uint64_t flag = 0;\ + flag += fread(&n_buckets, sizeof(n_buckets), 1, fp); \ + flag += fread(&(*h)->bits, sizeof((*h)->bits), 1, fp); \ + flag += fread(&(*h)->count, sizeof((*h)->count), 1, fp); \ (*h)->used = (khint32_t*)kmalloc(__kh_fsize(n_buckets) * sizeof(khint32_t)); \ (*h)->keys = (khkey_t*)kmalloc(n_buckets * sizeof(khkey_t)); \ - fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \ - fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); \ + flag += fread((*h)->used, sizeof(khint32_t), __kh_fsize(n_buckets), fp); \ + flag += fread((*h)->keys, sizeof(khkey_t), n_buckets, fp); \ return 1; \ } \