to replace get_new_candidates(); UNFINISHED!

Assembly.cpp

@@ -1,10 +1,9 @@
-#include "Assembly.h"
 #include <stdio.h>
 #include <stdlib.h>
 #include <zlib.h>
+#include "Assembly.h"
 #include "Process_Read.h"
 #include "CommandLines.h"
-#include "kmer.h"
 #include "Hash_Table.h"
 #include "POA.h"
 #include "Correct.h"
@@ -13,6 +12,7 @@
 
 void *ha_flt_tab;
 ha_pt_t *ha_idx;
+void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double band_width_threshold, int keep_whole_chain);
 
 All_reads R_INF;
 pthread_mutex_t statistics;
@@ -375,59 +375,6 @@ long long push_final_overlaps(ma_hit_t_alloc* paf, ma_hit_t_alloc* reverse_paf_l
     return available_overlaps;
 }
 
-void get_new_candidates(long long readID, UC_Read* g_read, overlap_region_alloc* overlap_list, k_mer_pos_list_alloc* array_list,
-						Candidates_list* l, double band_width_threshold, int keep_whole_chain)
-{
-    HPC_seq HPC_read;
-    Hash_code k_code;
-    long long avalible_k;
-    uint64_t code;
-    uint64_t end_pos;
-    k_mer_pos* list;
-    uint64_t list_length;
-    uint64_t sub_ID;
-
-    clear_Candidates_list(l);
-
-    clear_k_mer_pos_list_alloc(array_list);
-    clear_overlap_region_alloc(overlap_list);
-
-    recover_UC_Read(g_read, &R_INF, readID);
-
-    init_HPC_seq(&HPC_read, g_read->seq, g_read->length);
-    init_Hash_code(&k_code);
-    avalible_k = 0;
-
-    while ((code = get_HPC_code(&HPC_read, &end_pos)) != 6)
-    {
-        if(code < 4)
-        {
-            k_mer_append(&k_code, code, asm_opt.k_mer_length);
-            avalible_k++;
-            if (avalible_k >= asm_opt.k_mer_length)
-            {
-                //list_length = locate_Total_Pos_Table(&PCB, &k_code, &list, asm_opt.k_mer_length, &sub_ID);
-
-                if (list_length != 0)
-                {
-                    append_k_mer_pos_list_alloc(array_list, list, list_length, end_pos, 0);
-                }
-            }
-        }
-        else
-        {
-            avalible_k = 0;
-            init_Hash_code(&k_code);
-        }
-    }
-
-	// BIG CHANGES WILL GO HERE!!!
-
-    calculate_overlap_region_by_chaining(l, overlap_list, readID, g_read->length, &R_INF, 
-    band_width_threshold, keep_whole_chain);
-}
-
-
 void* Overlap_calculate_heap_merge(void* arg)
 {
     long long num_read_base = 0;
@@ -473,10 +420,14 @@ void* Overlap_calculate_heap_merge(void* arg)
     Round2_alignment second_round;
     init_Round2_alignment(&second_round);
 
+	ha_abuf_t *ab;
+	ab = ha_abuf_init();
+
     for (i = thr_ID; i < (long long)R_INF.total_reads; i = i + asm_opt.thread_num)
     {
         ///get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5);
-        get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.02, 1);
+        //get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.02, 1);
+		ha_get_new_candidates(ab, i, &g_read, &overlap_list, &l, 0.02, 1);
 
         clear_Cigar_record(&current_cigar);
         clear_Round2_alignment(&second_round);
@@ -507,7 +458,7 @@ void* Overlap_calculate_heap_merge(void* arg)
         push_overlaps(&(R_INF.reverse_paf[i]), &overlap_list, 2, &R_INF, asm_opt.roundID%2);
     }
 
-
+	ha_abuf_destroy(ab);
     finish_output_buffer();
     destory_buffer_sub_block(&current_sub_buffer);
     destory_Candidates_list(&l);
@@ -523,7 +474,6 @@ void* Overlap_calculate_heap_merge(void* arg)
     destory_Round2_alignment(&second_round);
 
 
-
     pthread_mutex_lock(&statistics);
     asm_opt.num_bases += num_read_base;
     asm_opt.num_corrected_bases += num_correct_base;
@@ -588,6 +538,8 @@ void* Output_related_reads(void* arg)
     Round2_alignment second_round;
     init_Round2_alignment(&second_round);
 
+	ha_abuf_t *ab;
+	ab = ha_abuf_init();
 
     long long required_read_name_length = strlen(asm_opt.required_read_name);
     for (i = thr_ID; i < (long long)R_INF.total_reads; i = i + asm_opt.thread_num)
@@ -597,8 +549,9 @@ void* Output_related_reads(void* arg)
             &&
            memcmp(asm_opt.required_read_name, Get_NAME((R_INF), i), Get_NAME_LENGTH((R_INF),i)) == 0)
         {
-            ////get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5);
+            //get_new_candidates(i, &g_read, &overlap_list, &array_list, &heap, &l, THRESHOLD_RATE*1.5);
-            get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.02, 1);
+            //get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.02, 1);
+			ha_get_new_candidates(ab, i, &g_read, &overlap_list, &l, 0.02, 1);
 
             fprintf(stderr, ">%.*s\n", (int)Get_NAME_LENGTH((R_INF), i),
             Get_NAME((R_INF), i));
@@ -618,6 +571,7 @@ void* Output_related_reads(void* arg)
         }
     }
 
+	ha_abuf_destroy(ab);
     finish_output_buffer();
 
     destory_buffer_sub_block(&current_sub_buffer);
@@ -1155,6 +1109,9 @@ void* Final_overlap_calculate_heap_merge(void* arg)
     Cigar_record_alloc cigarline;
     init_Cigar_record_alloc(&cigarline);
 
+	ha_abuf_t *ab;
+	ab = ha_abuf_init();
+
     uint8_t c2n[256];
     memset(c2n, 4, 256);
     c2n[(uint8_t)'A'] = c2n[(uint8_t)'a'] = 0; c2n[(uint8_t)'C'] = c2n[(uint8_t)'c'] = 1;
@@ -1163,7 +1120,9 @@ void* Final_overlap_calculate_heap_merge(void* arg)
     for (i = thr_ID; i < R_INF.total_reads; i = i + asm_opt.thread_num)
     {
 
-        get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0);
+        //get_new_candidates(i, &g_read, &overlap_list, &array_list, &l, 0.001, 0);
+		ha_get_new_candidates(ab, i, &g_read, &overlap_list, &l, 0.001, 0);
+
         /**
         correct_overlap(&overlap_list, &R_INF, &g_read, &correct, &overlap_read, &POA_Graph, &DAGCon,
         &matched_overlap_0, &matched_overlap_1, &potiental_matched_overlap_0, &potiental_matched_overlap_1,
@@ -1192,6 +1151,7 @@ void* Final_overlap_calculate_heap_merge(void* arg)
 
     }
 
+	ha_abuf_destroy(ab);
     finish_output_buffer();
 
     destory_Candidates_list(&l);

Hash_Table.cpp

@@ -1,12 +1,11 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <stdint.h>
+#include <pthread.h>
 #include "Hash_Table.h"
 #include "Process_Read.h"
 #include "Correct.h"
 #include "CommandLines.h"
-#include "kmer.h"
-#include <pthread.h>
 #include "ksort.h"
 pthread_mutex_t output_mutex;
 

Hash_Table.h

@@ -1,6 +1,5 @@
 #ifndef __HASHTABLE__
 #define __HASHTABLE__
-#include "kmer.h"
 #include "htab.h"
 
 #define PREFIX_BITS 16

Makefile

@@ -4,7 +4,7 @@ CPPFLAGS=
 INCLUDES=
 OBJS=		Output.o CommandLines.o Process_Read.o Assembly.o Hash_Table.o \
 			POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o \
-			htab.o hist.o sketch.o sys.o
+			htab.o hist.o sketch.o anchor.o sys.o
 EXE=		hifiasm
 LIBS=		-lz -lpthread -lm
 
@@ -33,24 +33,24 @@ depend:
 # DO NOT DELETE
 
 Assembly.o: Assembly.h CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h
-Assembly.o: kmer.h Hash_Table.h htab.h POA.h Correct.h Levenshtein_distance.h
+Assembly.o: Hash_Table.h htab.h POA.h Correct.h Levenshtein_distance.h
 Assembly.o: Output.h
 CommandLines.o: CommandLines.h ketopt.h
-Correct.o: Correct.h Hash_Table.h kmer.h Process_Read.h Overlaps.h kvec.h
+Correct.o: Correct.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h
-Correct.o: kdq.h CommandLines.h htab.h Levenshtein_distance.h POA.h
+Correct.o: kdq.h CommandLines.h Levenshtein_distance.h POA.h Assembly.h
-Correct.o: Assembly.h
+Hash_Table.o: Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h
-Hash_Table.o: Hash_Table.h kmer.h Process_Read.h Overlaps.h kvec.h kdq.h
+Hash_Table.o: CommandLines.h Correct.h Levenshtein_distance.h POA.h ksort.h
-Hash_Table.o: CommandLines.h htab.h Correct.h Levenshtein_distance.h POA.h
-Hash_Table.o: ksort.h
 Levenshtein_distance.o: Levenshtein_distance.h
 Output.o: Output.h CommandLines.h
 Overlaps.o: Overlaps.h kvec.h kdq.h ksort.h Process_Read.h CommandLines.h
-Overlaps.o: Hash_Table.h kmer.h htab.h Correct.h Levenshtein_distance.h POA.h
+Overlaps.o: Hash_Table.h htab.h Correct.h Levenshtein_distance.h POA.h
-POA.o: POA.h Hash_Table.h kmer.h Process_Read.h Overlaps.h kvec.h kdq.h
+POA.o: POA.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h
-POA.o: CommandLines.h htab.h Correct.h Levenshtein_distance.h
+POA.o: CommandLines.h Correct.h Levenshtein_distance.h
 Process_Read.o: Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
 Trio.o: khashl.h kthread.h Process_Read.h Overlaps.h kvec.h kdq.h
 Trio.o: CommandLines.h htab.h
+anchor.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
+anchor.o: ksort.h Hash_Table.h
 hist.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
 htab.o: kthread.h khashl.h kseq.h ksort.h htab.h Process_Read.h Overlaps.h
 htab.o: kvec.h kdq.h CommandLines.h

anchor.cpp

@@ -0,0 +1,96 @@
+#include "htab.h"
+#include "ksort.h"
+#include "Hash_Table.h"
+
+typedef struct {
+	uint64_t srt;
+	uint32_t self_off;
+	uint32_t other_off;
+} anchor1_t;
+
+#define an_key1(a) ((a).srt)
+#define an_key2(a) ((a).self_off)
+KRADIX_SORT_INIT(ha_an1, anchor1_t, an_key1, 8)
+KRADIX_SORT_INIT(ha_an2, anchor1_t, an_key2, 4)
+
+typedef struct {
+	int n;
+	const ha_idxpos_t *a;
+} seed1_t;
+
+struct ha_abuf_s {
+	uint64_t n_a, m_a;
+	uint32_t old_mz_m;
+	ha_mz1_v mz;
+	seed1_t *seed;
+	anchor1_t *a;
+};
+
+ha_abuf_t *ha_abuf_init(void)
+{
+	return (ha_abuf_t*)calloc(1, sizeof(ha_abuf_t));
+}
+
+void ha_abuf_destroy(ha_abuf_t *ab)
+{
+	free(ab->seed); free(ab->a); free(ab->mz.a); free(ab);
+}
+
+void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double band_width_threshold, int keep_whole_chain)
+{
+	extern void *ha_flt_tab;
+	extern ha_pt_t *ha_idx;
+	uint32_t i;
+	uint64_t k, l;
+
+    clear_Candidates_list(cl);
+    clear_overlap_region_alloc(overlap_list);
+	recover_UC_Read(ucr, &R_INF, rid);
+	ab->mz.n = 0, ab->n_a = 0;
+
+	ha_sketch(ucr->seq, ucr->length, asm_opt.mz_win, asm_opt.k_mer_length, 0, !asm_opt.no_HPC, &ab->mz, ha_flt_tab);
+	if (ab->mz.m > ab->old_mz_m) {
+		ab->old_mz_m = ab->mz.m;
+		REALLOC(ab->seed, ab->old_mz_m);
+	}
+	for (i = 0, ab->n_a = 0; i < ab->mz.n; ++i) {
+		ab->seed[i].a = ha_pt_get(ha_idx, ab->mz.a[i].x, &ab->seed[i].n);
+		ab->n_a += ab->seed[i].n;
+	}
+	if (ab->n_a > ab->m_a) {
+		ab->m_a = ab->n_a;
+		ab->m_a = ab->m_a > 16? ab->m_a + (ab->m_a>>1) : 16;
+		REALLOC(ab->a, ab->m_a);
+	}
+	for (i = 0, k = 0; i < ab->mz.n; ++i) {
+		uint32_t j;
+		ha_mz1_t *z = &ab->mz.a[i];
+		seed1_t *s = &ab->seed[i];
+		for (j = 0; j < s->n; ++j) {
+			const ha_idxpos_t *y = &s->a[j];
+			anchor1_t *an = &ab->a[k++];
+			uint8_t rev;
+			if (z->rev == y->rev) { // forward strand
+				rev = 0;
+				an->self_off = z->pos;
+				an->other_off = y->pos;
+			} else { // reverse strand
+				rev = 1;
+				an->self_off = ucr->length - 1 - (z->pos + 1 - z->span);
+				an->other_off = R_INF.read_length[y->rid] - 1 - (y->pos + 1 - y->span);
+			}
+			an->srt = (uint64_t)rev<<63 | (uint64_t)y->rid << 32 | (0x80000000ULL + (an->other_off - an->self_off));
+		}
+	}
+
+	radix_sort_ha_an1(ab->a, ab->a + ab->n_a);
+	for (k = 1, l = 0; k <= ab->n_a; ++k) {
+		if (k == ab->n_a || ab->a[k].srt != ab->a[l].srt) {
+			if (k - l > 1)
+				radix_sort_ha_an2(ab->a + l, ab->a + k);
+			l = k;
+		}
+	}
+
+	calculate_overlap_region_by_chaining(cl, overlap_list, rid, ucr->length, &R_INF, band_width_threshold, keep_whole_chain);
+}

htab.h

@@ -32,6 +32,9 @@ ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_f
 void ha_pt_destroy(ha_pt_t *h);
 const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n);
 
+ha_abuf_t *ha_abuf_init(void);
+void ha_abuf_destroy(ha_abuf_t *ab);
+
 double yak_cputime(void);
 void yak_reset_realtime(void);
 double yak_realtime(void);

kmer.h

@@ -1,81 +0,0 @@
-#ifndef __KMER__
-#define __KMER__
-#include "Process_Read.h"
-
-///#define ALL (0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff)
-#define ALL (0xffffffffffffffff)
-/****************************may have bugs********************************/
-#define SAFE_SHIFT(k) k & ((k < 64)?ALL:0)
-/****************************may have bugs********************************/
-
-typedef struct
-{
-    //can represent at most 64-mer
-    uint64_t x[4];
-} Hash_code;
-
-typedef struct
-{
-	char* str;
-	long long l;
-    long long i;
-    long long N_occ;
-
-} HPC_seq;
-
-inline void init_HPC_seq(HPC_seq* seq, char* str, long long l)
-{
-    seq->i = 0;
-    seq->l = l;
-    seq->N_occ = 0;
-    seq->str = str;
-}
-
-inline uint64_t get_HPC_code(HPC_seq* seq, uint64_t* end_pos)
-{
-
-    if(seq->i < seq ->l)
-    {
-        uint8_t code = seq_nt6_table[(uint8_t)seq->str[seq->i]];
-
-        (*end_pos) = seq->i;
-
-        for (; seq->i < seq->l; seq->i++)
-        {
-            ///number of Ns
-            if (seq_nt6_table[(uint8_t)seq->str[seq->i]] >= 4)
-            {
-                seq->N_occ++;
-            }
-            
-            if (seq_nt6_table[(uint8_t)seq->str[seq->i]] != code)
-            {
-                break;
-            }
-        }
-
-        return (uint64_t)code;
-    }
-    else
-    {
-        ///end
-        return 6;
-    }
-    
-}
-
-inline void init_Hash_code(Hash_code* code)
-{
-    code->x[0] = code->x[1] = code->x[2] = code->x[3] = 0;
-}
-
-inline void k_mer_append(Hash_code* code, uint64_t c, int k)
-{
-    uint64_t mask = ALL >> (64 - k), shift = k - 1;
-    code->x[0] = ((code->x[0]<<1) | (c&1))  & mask;
-    code->x[1] = ((code->x[1]<<1) | (c>>1)) & mask;
-	code->x[2] = code->x[2] >> 1 | (uint64_t)(1 - (c&1))  << shift;
-	code->x[3] = code->x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift;
-}
-
-#endif