From 0cc02035c54df59dc96e72ed55b5b552ea5ee06d Mon Sep 17 00:00:00 2001 From: chhylp123 Date: Mon, 26 Oct 2020 11:01:48 -0400 Subject: [PATCH] HiC full index --- CommandLines.cpp | 21 ++- CommandLines.h | 7 + Makefile | 3 +- Overlaps.cpp | 28 +++- hic.cpp | 328 +++++++++++++++++++++++++++++++++++++++++++++++ hic.h | 9 ++ 6 files changed, 390 insertions(+), 6 deletions(-) create mode 100644 hic.cpp create mode 100644 hic.h diff --git a/CommandLines.cpp b/CommandLines.cpp index 90f605a..3495719 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -25,6 +25,8 @@ static ko_longopt_t long_options[] = { { "pri-range", ko_required_argument, 310 }, { "high-het", ko_no_argument, 311 }, { "lowQ", ko_required_argument, 312 }, + { "h1", ko_required_argument, 313 }, + { "h2", ko_required_argument, 314 }, { 0, 0, 0 } }; @@ -97,8 +99,11 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->read_file_names = NULL; asm_opt->output_file_name = (char*)(DEFAULT_OUTPUT); asm_opt->required_read_name = NULL; + asm_opt->hic_files[0] = NULL; + asm_opt->hic_files[1] = NULL; asm_opt->thread_num = 1; asm_opt->k_mer_length = 51; + asm_opt->hic_mer_length = 31; asm_opt->mz_win = 51; asm_opt->bf_shift = 37; asm_opt->high_factor = 5.0; @@ -325,7 +330,7 @@ int check_option(hifiasm_opt_t* asm_opt) if(asm_opt->bed_inconsist_rate < 0 || asm_opt->bed_inconsist_rate > 100) { - fprintf(stderr, "[ERROR] inconsistency rate should be [0, 100] (--pb-range)\n"); + fprintf(stderr, "[ERROR] inconsistency rate should be [0, 100] (--lowQ)\n"); return 0; } @@ -335,6 +340,18 @@ int check_option(hifiasm_opt_t* asm_opt) if(asm_opt->fn_bin_list[0] != NULL && check_file(asm_opt->fn_bin_list[0], "LIST1") == 0) return 0; if(asm_opt->fn_bin_list[1] != NULL && check_file(asm_opt->fn_bin_list[1], "LIST2") == 0) return 0; if(asm_opt->required_read_name != NULL && check_file(asm_opt->required_read_name, "b") == 0) return 0; + if(asm_opt->hic_files[0] != NULL && check_file(asm_opt->hic_files[0], "HIC1") == 0) return 0; + if(asm_opt->hic_files[1] != NULL && check_file(asm_opt->hic_files[1], "HIC2") == 0) return 0; + if(asm_opt->hic_files[0] != NULL && asm_opt->hic_files[1] == NULL) + { + fprintf(stderr, "[ERROR] lack r2 of HiC reads (--h2)\n"); + return 0; + } + if(asm_opt->hic_files[1] != NULL && asm_opt->hic_files[0] == NULL) + { + fprintf(stderr, "[ERROR] lack r1 of HiC reads (--h1)\n"); + return 0; + } // fprintf(stderr, "input file num: %d\n", asm_opt->num_reads); // fprintf(stderr, "output file: %s\n", asm_opt->output_file_name); // fprintf(stderr, "number of threads: %d\n", asm_opt->thread_num); @@ -449,6 +466,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) } else if (c == 311) asm_opt->flag |= HA_F_HIGH_HET; else if (c == 312) asm_opt->bed_inconsist_rate = atoi(opt.arg); + else if (c == 313) asm_opt->hic_files[0] = opt.arg; + else if (c == 314) asm_opt->hic_files[1] = opt.arg; else if (c == 'l') { ///0: disable purge_dup; 1: purge containment; 2: purge overlap asm_opt->purge_level_primary = asm_opt->purge_level_trio = atoi(opt.arg); diff --git a/CommandLines.h b/CommandLines.h index b130949..1a600ed 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -30,9 +30,11 @@ typedef struct { char *fn_bin_yak[2]; char *fn_bin_list[2]; char *extract_list; + char *hic_files[2]; int extract_iter; int thread_num; int k_mer_length; + int hic_mer_length; int mz_win; int bf_shift; double high_factor; // coverage cutoff set to high_factor*hom_cov @@ -93,4 +95,9 @@ static inline int ha_opt_triobin(const hifiasm_opt_t *opt) return ((opt->fn_bin_yak[0] && opt->fn_bin_yak[1]) || (opt->fn_bin_list[0] && opt->fn_bin_list[1])); } +static inline int ha_opt_hic(const hifiasm_opt_t *opt) +{ + return ((opt->hic_files[0] && opt->hic_files[1])); +} + #endif diff --git a/Makefile b/Makefile index 7da5789..984e9a0 100644 --- a/Makefile +++ b/Makefile @@ -6,7 +6,7 @@ CPPFLAGS= INCLUDES= OBJS= CommandLines.o Process_Read.o Assembly.o Hash_Table.o \ POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o \ - htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o + htab.o hist.o sketch.o anchor.o extract.o sys.o ksw2_extz2_sse.o hic.o EXE= hifiasm LIBS= -lz -lpthread -lm @@ -71,3 +71,4 @@ main.o: CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h Assembly.h main.o: Levenshtein_distance.h htab.h sketch.o: kvec.h htab.h Process_Read.h Overlaps.h kdq.h CommandLines.h sys.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h +hic.o: hic.h diff --git a/Overlaps.cpp b/Overlaps.cpp index f86ac7c..2d563f3 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -9,6 +9,7 @@ #include "Hash_Table.h" #include "Correct.h" #include "Purge_Dups.h" +#include "hic.h" uint32_t debug_purge_dup = 0; @@ -11630,6 +11631,21 @@ ma_hit_t_alloc* sources, R_to_U* ruIndex, int max_hang, int min_ovlp) } +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) +{ + 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); + + ma_ug_destroy(ug); + kv_destroy(new_rtg_edges.a); +} + void merge_unitig_content(ma_utg_t* collection, ma_ug_t* ug, asg_t* read_g, kvec_asg_arc_t_warp* edge) { if(collection->m == 0) return; @@ -23738,7 +23754,7 @@ void pre_clean(ma_hit_t_alloc* sources, ma_sub_t* coverage_cut, asg_t *sg, long ///remove isoloated single read tri_flag += asg_arc_del_single_node_directly(sg, asm_opt.max_short_tip, sources); - if (!ha_opt_triobin(&asm_opt)) + if ((!ha_opt_triobin(&asm_opt))&&(!ha_opt_hic(&asm_opt))) { tri_flag += asg_arc_del_triangular_advance(sg, bubble_dist); ///remove the cross at the bubble carefully, just remove inexact cross @@ -25057,7 +25073,7 @@ kvec_t_u32_warp* new_rtg_nodes) int tri_flag = 0; ///remove very simple circle tri_flag += asg_arc_del_simple_circle_untig(sources, coverage_cut, r_g, 100, 0); - if (!ha_opt_triobin(&asm_opt)) + if ((!ha_opt_triobin(&asm_opt))&&(!ha_opt_hic(&asm_opt))) { ///remove isoloated single read tri_flag += asg_arc_del_triangular_advance(r_g, bubble_dist); @@ -25328,7 +25344,7 @@ kvec_asg_arc_t_warp* new_rtg_edges) int tri_flag = 0; ///remove very simple circle tri_flag += asg_arc_del_simple_circle_untig(sources, coverage_cut, r_g, 100, 0); - if (!ha_opt_triobin(&asm_opt)) + if ((!ha_opt_triobin(&asm_opt))&&(!ha_opt_hic(&asm_opt))) { ///remove isoloated single read tri_flag += asg_arc_del_triangular_advance(r_g, bubble_dist); @@ -25706,7 +25722,7 @@ uint32_t is_bubble_check, uint32_t is_primary_check) int tri_flag = 0; ///remove very simple circle tri_flag += asg_arc_del_simple_circle_untig(sources, coverage_cut, r_g, 100, 0); - if (!ha_opt_triobin(&asm_opt)) + if ((!ha_opt_triobin(&asm_opt))&&(!ha_opt_hic(&asm_opt))) { ///remove isoloated single read tri_flag += asg_arc_del_triangular_advance(r_g, bubble_dist); @@ -27450,6 +27466,10 @@ ma_sub_t **coverage_cut_ptr, int debug_g) reverse_sources, bubble_dist, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length, mini_overlap_length); } + else if(ha_opt_hic(&asm_opt)) + { + output_hic_graph(sg, coverage_cut, output_file_name, sources, ruIndex, max_hang_length, mini_overlap_length);; + } else { diff --git a/hic.cpp b/hic.cpp new file mode 100644 index 0000000..a80bfbc --- /dev/null +++ b/hic.cpp @@ -0,0 +1,328 @@ +#define __STDC_LIMIT_MACROS +#include "hic.h" +#include "htab.h" +#include "assert.h" +#include "Overlaps.h" +#include "khashl.h" + +#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) + +typedef struct { + hc_pt_t *h; + uint64_t n; + uint64_t *a; + khint_t end; +} hc_pt1_t; + +typedef struct { + ma_ug_t* ug; + uint64_t uID_bits; + uint64_t uID_mode; + uint64_t pos_bits; + uint64_t pos_mode; + uint64_t rev_mode; + hc_pt1_t idx; + uint64_t k; +} ha_ug_index; + +ha_ug_index* ug_index; + +inline uint64_t get_k_direction(uint64_t x[4]) +{ + if(x[1] != x[3]) + { + return x[1] < x[3]? 0 : 1; + } + else if(x[0] != x[2]) + { + return x[0] < x[2]? 0 : 1; + } + else + { + return (uint64_t)-1; + } +} + +inline uint64_t hc_hash_long(uint64_t x[4], uint64_t* skip) +{ + ///compare forward k-mer and reverse complementary strand + (*skip) = get_k_direction(x); + if((*skip) == (uint64_t)-1) return (*skip); + 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) +{ + uint64_t beg; + khint_t k; + k = hc_pt_get(idx->idx.h, key); + if (k == kh_end(idx->idx.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) + { + if (kh_exist(idx->idx.h, k)) + { + return kh_val(idx->idx.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 + } +} + +void test_hc_pt1(char* seq, uint64_t len, uint64_t uID, ha_ug_index* idx) +{ + uint64_t i, l, k, 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); + if(cnt == 0) fprintf(stderr, "ERROR cnt, uID: %lu\n", uID); + for (k = 0; k < cnt; k++) + { + if(pos_list[k]==pos) + { + pos_list[k] = (uint64_t)-1; + break; + } + } + if(k == cnt) fprintf(stderr, "ERROR k\n"); + + } + } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart + } +} + +void test_unitig_index(ha_ug_index* idx) +{ + uint32_t i; + ma_utg_t *u = NULL; + 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++) + { + if(idx->idx.a[i]!=(uint64_t)-1) + { + fprintf(stderr, "ERROR i\n"); + } + } +} + +void hc_pt_t_gen(hc_pt1_t* pt) +{ + khint_t k; + uint64_t c; + for (k = 0, pt->n = 0; k != kh_end(pt->h); ++k) { + if (kh_exist(pt->h, k)) { + c = kh_val(pt->h, k); + kh_val(pt->h, k) = pt->n; + pt->n += c; + pt->end = k; + } + } + 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); + sprintf(gfa_name, "%s.hc_tlb", file_name); + FILE* fp = fopen(gfa_name, "w"); + if (!fp) { + free(gfa_name); + return 0; + } + + fwrite(&idx->uID_bits, sizeof(idx->uID_bits), 1, fp); + fwrite(&idx->uID_mode, sizeof(idx->uID_mode), 1, fp); + fwrite(&idx->pos_bits, sizeof(idx->pos_bits), 1, fp); + 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); + + + fprintf(stderr, "[M::%s] Index has been written.\n", __func__); + free(gfa_name); + fclose(fp); + return 1; +} + + +int load_hc_pt_index(ha_ug_index** r_idx, char* file_name) +{ + char* gfa_name = (char*)malloc(strlen(file_name)+25); + sprintf(gfa_name, "%s.hc_tlb", file_name); + FILE* fp = fopen(gfa_name, "r"); + if (!fp) { + free(gfa_name); + return 0; + } + 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); + (*r_idx) = idx; + + fprintf(stderr, "[M::%s] Index has been loaded.\n", __func__); + free(gfa_name); + fclose(fp); + return 1; +} + + +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); + double index_time = yak_realtime(); + 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(); + + 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); + } + + hc_pt_t_gen(&(idx->idx)); + + 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] ==> HiC index has been built\n", __func__, yak_realtime()-index_time); + + return idx; +} + +void destory_hc_pt_index(ha_ug_index* r_idx) +{ + if(r_idx->idx.h) hc_pt_destroy(r_idx->idx.h); + if(r_idx->idx.a) free(r_idx->idx.a); +} + +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); + + destory_hc_pt_index(ug_index); +} \ No newline at end of file diff --git a/hic.h b/hic.h new file mode 100644 index 0000000..84c0114 --- /dev/null +++ b/hic.h @@ -0,0 +1,9 @@ +#ifndef __HIC__ +#define __HIC__ +#include +#include "Overlaps.h" + + +void hic_analysis(ma_ug_t *ug); + +#endif