HiC full index

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);

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

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

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
     {
 

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<<idx->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<<idx->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<<idx->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<<idx->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);
+}

hic.h

@@ -0,0 +1,9 @@
+#ifndef __HIC__
+#define __HIC__
+#include <stdint.h>
+#include "Overlaps.h"
+
+
+void hic_analysis(ma_ug_t *ug);
+
+#endif