diff --git a/Assembly.cpp b/Assembly.cpp index 6aa1cbd..39e4d63 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -1,10 +1,9 @@ -#include "Assembly.h" #include #include #include +#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(¤t_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(¤t_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, &l, 0.02, 1); + //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); + 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(¤t_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); diff --git a/Hash_Table.cpp b/Hash_Table.cpp index d387da7..a589a02 100644 --- a/Hash_Table.cpp +++ b/Hash_Table.cpp @@ -1,12 +1,11 @@ #include #include #include +#include #include "Hash_Table.h" #include "Process_Read.h" #include "Correct.h" #include "CommandLines.h" -#include "kmer.h" -#include #include "ksort.h" pthread_mutex_t output_mutex; diff --git a/Hash_Table.h b/Hash_Table.h index 1dd5ccd..279fecf 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -1,6 +1,5 @@ #ifndef __HASHTABLE__ #define __HASHTABLE__ -#include "kmer.h" #include "htab.h" #define PREFIX_BITS 16 diff --git a/Makefile b/Makefile index 0b221b4..2c22c33 100644 --- a/Makefile +++ b/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: kdq.h CommandLines.h htab.h Levenshtein_distance.h POA.h -Correct.o: Assembly.h -Hash_Table.o: Hash_Table.h kmer.h Process_Read.h Overlaps.h kvec.h kdq.h -Hash_Table.o: CommandLines.h htab.h Correct.h Levenshtein_distance.h POA.h -Hash_Table.o: ksort.h +Correct.o: Correct.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h +Correct.o: kdq.h CommandLines.h Levenshtein_distance.h POA.h Assembly.h +Hash_Table.o: Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.h +Hash_Table.o: CommandLines.h Correct.h Levenshtein_distance.h POA.h 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 -POA.o: POA.h Hash_Table.h kmer.h Process_Read.h Overlaps.h kvec.h kdq.h -POA.o: CommandLines.h htab.h Correct.h Levenshtein_distance.h +Overlaps.o: Hash_Table.h htab.h Correct.h Levenshtein_distance.h POA.h +POA.o: POA.h Hash_Table.h htab.h Process_Read.h Overlaps.h kvec.h kdq.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 diff --git a/anchor.cpp b/anchor.cpp new file mode 100644 index 0000000..445dd33 --- /dev/null +++ b/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); +} diff --git a/htab.h b/htab.h index 18567f8..180340a 100644 --- a/htab.h +++ b/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); diff --git a/kmer.h b/kmer.h deleted file mode 100644 index a06c822..0000000 --- a/kmer.h +++ /dev/null @@ -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