diff --git a/Assembly.cpp b/Assembly.cpp index 407ee32..7c1b527 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -9,7 +9,7 @@ #include "POA.h" #include "Correct.h" #include "Output.h" -#include "Trio.h" +#include "yak.h" Total_Count_Table TCB; Total_Pos_Table PCB; diff --git a/Assembly.h b/Assembly.h index 557aa8a..81e1104 100644 --- a/Assembly.h +++ b/Assembly.h @@ -8,6 +8,7 @@ #define Get_Cigar_Type(RECORD) (RECORD&3) #define Get_Cigar_Length(RECORD) (RECORD>>2) +void ha_count_high(const hifiasm_opt_t *asm_opt); void Counting_multiple_thr(); void Build_hash_table_multiple_thr(); void Overlap_calculate_multipe_thr(); diff --git a/CommandLines.cpp b/CommandLines.cpp index ac424e6..da2e587 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -28,6 +28,7 @@ void Print_H(hifiasm_opt_t* asm_opt) fprintf(stderr, " -r INT round of correction [%d]\n", asm_opt->number_of_round); fprintf(stderr, " -a INT round of assembly cleaning [%d]\n", asm_opt->clean_round); fprintf(stderr, " -k INT k-mer length [%d] (must be < 64)\n", asm_opt->k_mer_length); + fprintf(stderr, " -f INT number of bits for bloom filter [%d]\n", asm_opt->bf_shift); ///fprintf(stderr, " -w write all overlaps to disk, can accelerate assembly next time [%d]\n", asm_opt->write_index_to_disk); ///fprintf(stderr, " -l load all overlaps from disk, can avoid overlap calculation [%d]\n", asm_opt->load_index_from_disk); ///fprintf(stderr, " -i ignore saved overlaps in *.ovlp*.bin files\n"); @@ -63,7 +64,8 @@ void init_opt(hifiasm_opt_t* asm_opt) asm_opt->pat_index = NULL; asm_opt->mat_index = NULL; asm_opt->thread_num = 1; - asm_opt->k_mer_length = 39; + asm_opt->k_mer_length = 41; + asm_opt->bf_shift = 37; asm_opt->k_mer_min_freq = 3; asm_opt->k_mer_max_freq = 66; asm_opt->load_index_from_disk = 1; @@ -295,7 +297,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) int c; - while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:lwm:n:r:a:b:z:x:y:p:c:d:M:P:i", 0)) >= 0) { + while ((c = ketopt(&opt, argc, argv, 1, "hvt:o:k:lwm:n:r:a:b:z:x:y:p:c:d:M:P:if:", 0)) >= 0) { if (c == 'h') { Print_H(asm_opt); @@ -305,7 +307,8 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt) { fprintf(stderr, "[Version] %s\n", VERSION); return 0; - } + } + else if (c == 'f') asm_opt->bf_shift = atoi(opt.arg); else if (c == 't') asm_opt->thread_num = atoi(opt.arg); else if (c == 'o') asm_opt->output_file_name = opt.arg; else if (c == 'r') asm_opt->number_of_round = atoi(opt.arg); diff --git a/CommandLines.h b/CommandLines.h index 7279cdb..b3f7f28 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -15,6 +15,7 @@ typedef struct { char* mat_index; int thread_num; int k_mer_length; + int bf_shift; int k_mer_min_freq; int k_mer_max_freq; int load_index_from_disk; @@ -52,4 +53,4 @@ void clear_opt(hifiasm_opt_t* asm_opt, int last_round); int CommandLine_process (int argc, char *argv[], hifiasm_opt_t* asm_opt); double Get_T(void); -#endif \ No newline at end of file +#endif diff --git a/Hash_Table.h b/Hash_Table.h index 458312d..c4ae1c8 100644 --- a/Hash_Table.h +++ b/Hash_Table.h @@ -2,6 +2,7 @@ #define __HASHTABLE__ #include "khashl.h" #include "kmer.h" +#include "yak.h" KHASHL_MAP_INIT(static inline, Count_Table, ha_ct, uint64_t, int, kh_hash_dummy, kh_eq_generic) KHASHL_MAP_INIT(static inline, Pos_Table, ha_pt, uint64_t, uint64_t, kh_hash_dummy, kh_eq_generic) diff --git a/Makefile b/Makefile index 4d00034..89b9d53 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,8 @@ CXXFLAGS= -g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall 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 + POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o \ + yak-bbf.o yak-count.o yak-sys.o EXE= hifiasm LIBS= -lz -lpthread -lm @@ -32,25 +33,28 @@ depend: # DO NOT DELETE Assembly.o: Assembly.h CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h -Assembly.o: kdq.h kmer.h Hash_Table.h khashl.h POA.h Correct.h -Assembly.o: Levenshtein_distance.h Output.h Trio.h +Assembly.o: kdq.h kmer.h Hash_Table.h khashl.h yak.h POA.h Correct.h +Assembly.o: Levenshtein_distance.h Output.h CommandLines.o: CommandLines.h ketopt.h Correct.o: Correct.h Hash_Table.h khashl.h kmer.h Process_Read.h kseq.h -Correct.o: Overlaps.h kvec.h kdq.h CommandLines.h Levenshtein_distance.h -Correct.o: POA.h Assembly.h +Correct.o: Overlaps.h kvec.h kdq.h CommandLines.h yak.h +Correct.o: Levenshtein_distance.h POA.h Assembly.h Hash_Table.o: Hash_Table.h khashl.h kmer.h Process_Read.h kseq.h Overlaps.h -Hash_Table.o: kvec.h kdq.h CommandLines.h Correct.h Levenshtein_distance.h -Hash_Table.o: POA.h ksort.h +Hash_Table.o: kvec.h kdq.h CommandLines.h yak.h Correct.h +Hash_Table.o: 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 kseq.h -Overlaps.o: CommandLines.h Hash_Table.h khashl.h kmer.h Correct.h +Overlaps.o: CommandLines.h Hash_Table.h khashl.h kmer.h yak.h Correct.h Overlaps.o: Levenshtein_distance.h POA.h POA.o: POA.h Hash_Table.h khashl.h kmer.h Process_Read.h kseq.h Overlaps.h -POA.o: kvec.h kdq.h CommandLines.h Correct.h Levenshtein_distance.h +POA.o: kvec.h kdq.h CommandLines.h yak.h Correct.h Levenshtein_distance.h Process_Read.o: Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h Trio.o: khashl.h kthread.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h -Trio.o: CommandLines.h Trio.h kmer.h +Trio.o: CommandLines.h yak.h kthread.o: kthread.h main.o: CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h main.o: Assembly.h Levenshtein_distance.h +yak-bbf.o: yak.h +yak-count.o: CommandLines.h yak.h khashl.h kthread.h kseq.h +yak-sys.o: yak.h diff --git a/Overlaps.cpp b/Overlaps.cpp index 72ab35e..fbd01d6 100644 --- a/Overlaps.cpp +++ b/Overlaps.cpp @@ -46,9 +46,6 @@ void ma_hit_sort_qns(ma_hit_t *a, long long n) radix_sort_hit_qns(a, a + n); } - - - asg_t *asg_init(void) { return (asg_t*)calloc(1, sizeof(asg_t)); diff --git a/Trio.cpp b/Trio.cpp index 4bf4b8b..ff2f7eb 100644 --- a/Trio.cpp +++ b/Trio.cpp @@ -6,45 +6,35 @@ #include "khashl.h" // hash table #include "kthread.h" #include "Process_Read.h" -#include "Trio.h" +#include "yak.h" #include "CommandLines.h" -#include "kmer.h" -#define CALLOC(ptr, len) ((ptr) = (__typeof__(ptr))calloc((len), sizeof(*(ptr)))) -#define MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr)))) -#define REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr)))) +#define YAK_MAX_KMER 31 +#define YAK_COUNTER_BITS 10 +#define YAK_N_COUNTS (1<>YAK_COUNTER_BITS == (b)>>YAK_COUNTER_BITS) // lower 8 bits for counts; higher bits for k-mer #define yak_ch_hash(a) ((a)>>YAK_COUNTER_BITS) -KHASHL_SET_INIT(, yak_ht_t, yak_ht, uint64_t, yak_ch_hash, yak_ch_eq) +KHASHL_SET_INIT(static klib_unused, yak_ht_t, yak_ht, uint64_t, yak_ch_hash, yak_ch_eq) -///#define CHUNK_SIZE 200000 +typedef struct { + struct yak_ht_t *h; +} yak_ch1_t; -unsigned char seq_nt4_table[256] = { // translate ACGT to 0123 - 0, 1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, - 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 -}; +typedef struct { + int k, pre, n_hash, n_shift; + uint64_t tot; + yak_ch1_t *h; +} yak_ch_t; -static inline uint64_t yak_hash_long(uint64_t x[4]) -{ - int j = x[1] < x[3]? 0 : 1; - return yak_hash64_64(x[j<<1|0]) + yak_hash64_64(x[j<<1|1]); -} - -int yak_ch_get(const yak_ch_t *h, uint64_t x) +static int yak_ch_get(const yak_ch_t *h, uint64_t x) { int mask = (1<pre) - 1; yak_ht_t *g = h->h[x&mask].h; @@ -53,21 +43,7 @@ int yak_ch_get(const yak_ch_t *h, uint64_t x) return k == kh_end(g)? -1 : kh_key(g, k)&YAK_MAX_COUNT; } -yak_bf_t *yak_bf_init(int n_shift, int n_hashes) -{ - yak_bf_t *b; - void *ptr = 0; - if (n_shift + YAK_BLK_SHIFT > 64 || n_shift < YAK_BLK_SHIFT) return 0; - b = (yak_bf_t*)calloc(1, sizeof(yak_bf_t)); - b->n_shift = n_shift; - b->n_hashes = n_hashes; - posix_memalign(&ptr, 1<<(YAK_BLK_SHIFT-3), 1ULL<<(n_shift-3)); - b->b = (uint8_t*)ptr; - bzero(b->b, 1ULL<<(n_shift-3)); - return b; -} - -yak_ch_t *yak_ch_init(int k, int pre, int n_hash, int n_shift) +static yak_ch_t *yak_ch_init(int k, int pre) { yak_ch_t *h; int i; @@ -77,32 +53,10 @@ yak_ch_t *yak_ch_init(int k, int pre, int n_hash, int n_shift) CALLOC(h->h, 1<pre); for (i = 0; i < 1<pre; ++i) h->h[i].h = yak_ht_init(); - if (n_hash > 0 && n_shift > h->pre) { - h->n_hash = n_hash, h->n_shift = n_shift; - for (i = 0; i < 1<pre; ++i) - h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash); - } return h; } -void yak_bf_destroy(yak_bf_t *b) -{ - if (b == 0) return; - free(b->b); free(b); -} - -void yak_ch_destroy_bf(yak_ch_t *h) -{ - int i; - for (i = 0; i < 1<pre; ++i) { - if (h->h[i].b) - yak_bf_destroy(h->h[i].b); - h->h[i].b = 0; - } -} - - -yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ...) +static yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ...) { va_list ap; FILE *fp; @@ -138,7 +92,7 @@ yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ...) return 0; } - ch = ch0 == 0? yak_ch_init(t[0], t[1], 0, 0) : ch0; + ch = ch0 == 0? yak_ch_init(t[0], t[1]) : ch0; assert((int)t[0] == ch->k && (int)t[1] == ch->pre); for (i = 0; i < 1<pre; ++i) { yak_ht_t *h = ch->h[i].h; @@ -172,17 +126,39 @@ yak_ch_t *yak_ch_restore_core(yak_ch_t *ch0, const char *fn, int mode, ...) return ch; } - -void yak_ch_destroy(yak_ch_t *h) +static void yak_ch_destroy(yak_ch_t *h) { int i; if (h == 0) return; - yak_ch_destroy_bf(h); for (i = 0; i < 1<pre; ++i) yak_ht_destroy(h->h[i].h); free(h->h); free(h); } +typedef struct { + int max; + uint32_t *s; +} tb_buf_t; + +typedef struct { + int k, n_threads, print_diff; + double ratio_thres; + const yak_ch_t *ch; + tb_buf_t *buf; + UC_Read *bseq; + All_reads* seq; +} tb_shared_t; + +typedef struct { + int c[16]; + int sc[2]; + int nk; +} tb_cnt_t; + +typedef struct { + int n_seq; + tb_shared_t *aux; +} tb_step_t; static char tb_classify(const int sc[2], const int *c, int k, double ratio_thres) { @@ -206,8 +182,6 @@ static char tb_classify(const int sc[2], const int *c, int k, double ratio_thres static void tb_worker(void *_data, long k, int tid) { - ///tb_step_t *t = (tb_step_t*)_data; - ///tb_shared_t *aux = t->aux; tb_shared_t *aux = (tb_shared_t*)_data; UC_Read *s = &aux->bseq[tid]; recover_UC_Read(s, aux->seq, k); @@ -243,12 +217,7 @@ static void tb_worker(void *_data, long k, int tid) if (++l >= aux->k) { int type = 0, c1, c2; uint64_t y; - - - //++t->cnt[k].nk; ++cnt.nk; - - if (aux->ch->k < 32) y = yak_hash64(x[0] < x[1]? x[0] : x[1], mask); else @@ -259,28 +228,18 @@ static void tb_worker(void *_data, long k, int tid) if (c1 == 2 && c2 == 0) type = 1; else if (c2 == 2 && c1 == 0) type = 2; b->s[i] = type; - - - ///++t->cnt[k].c[flag]; ++cnt.c[flag]; - - // if (aux->print_diff && (flag>>2&3) != (flag&3)) - // printf("D\t%s\t%d\t%d\t%d\n", s->name, i, flag&3, flag>>2&3); } } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; } for (l = 0, i = 1; i <= s->length; ++i) { if (i == s->length || b->s[i] != b->s[l]) { if (b->s[l] > 0 && i - l >= aux->k - 4) - { - ///t->cnt[k].sc[b->s[l] - 1] += i - l; cnt.sc[b->s[l] - 1] += i - l; - } l = i; } } - int *c = cnt.c; char type; type = tb_classify(cnt.sc, c, aux->k, aux->ratio_thres); @@ -289,12 +248,9 @@ static void tb_worker(void *_data, long k, int tid) if(type == 'm') aux->seq->trio_flag[k] = MOTHER; } - - -void trio_partition() +void trio_partition(void) { - if(asm_opt.pat_index == NULL || asm_opt.mat_index == NULL) - { + if (asm_opt.pat_index == NULL || asm_opt.mat_index == NULL) { memset(R_INF.trio_flag, AMBIGU, R_INF.total_reads*sizeof(uint8_t)); return; } @@ -303,32 +259,26 @@ void trio_partition() fprintf(stderr, "Start trio binning ...... \n"); yak_ch_t *ch; - int i/**, min_cnt = 2, mid_cnt = 5**/; + int i /**, min_cnt = 2, mid_cnt = 5**/; tb_shared_t aux; memset(&aux, 0, sizeof(tb_shared_t)); aux.n_threads = asm_opt.thread_num, aux.print_diff = 0; aux.ratio_thres = 0.33; aux.seq = &R_INF; - ch = yak_ch_restore_core(0, asm_opt.pat_index, YAK_LOAD_TRIOBIN1, asm_opt.min_cnt, asm_opt.mid_cnt); ch = yak_ch_restore_core(ch, asm_opt.mat_index, YAK_LOAD_TRIOBIN2, asm_opt.min_cnt, asm_opt.mid_cnt); - - aux.k = ch->k; aux.ch = ch; aux.buf = (tb_buf_t*)calloc(aux.n_threads, sizeof(tb_buf_t)); aux.bseq = (UC_Read*)calloc(aux.n_threads, sizeof(UC_Read)); for (i = 0; i < aux.n_threads; ++i) - { init_UC_Read(&aux.bseq[i]); - } kt_for(aux.n_threads, tb_worker, &aux, aux.seq->total_reads); - for (i = 0; i < aux.n_threads; ++i) - { + for (i = 0; i < aux.n_threads; ++i) { free(aux.buf[i].s); destory_UC_Read(&aux.bseq[i]); } diff --git a/Trio.h b/Trio.h deleted file mode 100644 index d5319f5..0000000 --- a/Trio.h +++ /dev/null @@ -1,71 +0,0 @@ -#ifndef __TRIO__ -#define __TRIO__ -#include - -#define YAK_MAX_KMER 31 -#define YAK_COUNTER_BITS 10 -#define YAK_N_COUNTS (1<x[3] = code->x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift; } -static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function -{ - key = (~key + (key << 21)) & mask; // key = (key << 21) - key - 1; - key = key ^ key >> 24; - key = ((key + (key << 3)) + (key << 8)) & mask; // key * 265 - key = key ^ key >> 14; - key = ((key + (key << 2)) + (key << 4)) & mask; // key * 21 - key = key ^ key >> 28; - key = (key + (key << 31)) & mask; - return key; -} - -static inline uint64_t yak_hash64_64(uint64_t key) -{ - key = ~key + (key << 21); - key = key ^ key >> 24; - key = (key + (key << 3)) + (key << 8); - key = key ^ key >> 14; - key = (key + (key << 2)) + (key << 4); - key = key ^ key >> 28; - key = key + (key << 31); - return key; -} - -inline void Hashcode_to_string(Hash_code* code, char* str, int k) // FIXME: not working -{ - uint8_t c; - int i; - for (i = 0; i < k; i++) - { - c = (code->x[1] >> (k - i - 1)) & ((uint64_t)1); - c = c << 1; - c = c | ((code->x[0] >> (k - i - 1)) & ((uint64_t)1)); - - str[i] = s_H[c]; - } - -} - -void init_HPC_seq(HPC_seq* seq, char* str, long long l); - #endif diff --git a/main.cpp b/main.cpp index 76fcfd8..c666d00 100644 --- a/main.cpp +++ b/main.cpp @@ -4,16 +4,26 @@ #include "Process_Read.h" #include "Assembly.h" #include "Levenshtein_distance.h" +#include "yak.h" int main(int argc, char *argv[]) { + int i; + init_opt(&asm_opt); if (!CommandLine_process(argc, argv, &asm_opt)) return 1; - Correct_Reads(asm_opt.number_of_round); - + yak_reset_realtime(); + ha_count_high(&asm_opt); + if (0) { + Correct_Reads(asm_opt.number_of_round); + } destory_opt(&asm_opt); - + fprintf(stderr, "[M::%s] Version: %s\n", __func__, "dummy"); + fprintf(stderr, "[M::%s] CMD:", __func__); + for (i = 0; i < argc; ++i) + fprintf(stderr, " %s", argv[i]); + fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, yak_realtime(), yak_cputime(), yak_peakrss() / 1024.0 / 1024.0 / 1024.0); return 0; } diff --git a/yak-bbf.cpp b/yak-bbf.cpp new file mode 100644 index 0000000..742b80a --- /dev/null +++ b/yak-bbf.cpp @@ -0,0 +1,42 @@ +#include +#include +#include "yak.h" + +yak_bf_t *yak_bf_init(int n_shift, int n_hashes) +{ + yak_bf_t *b; + void *ptr = 0; + if (n_shift + YAK_BLK_SHIFT > 64 || n_shift < YAK_BLK_SHIFT) return 0; + CALLOC(b, 1); + b->n_shift = n_shift; + b->n_hashes = n_hashes; + posix_memalign(&ptr, 1<<(YAK_BLK_SHIFT-3), 1ULL<<(n_shift-3)); + b->b = (uint8_t*)ptr; + bzero(b->b, 1ULL<<(n_shift-3)); + return b; +} + +void yak_bf_destroy(yak_bf_t *b) +{ + if (b == 0) return; + free(b->b); free(b); +} + +int yak_bf_insert(yak_bf_t *b, uint64_t hash) +{ + int x = b->n_shift - YAK_BLK_SHIFT; + uint64_t y = hash & ((1ULL<> x & YAK_BLK_MASK; + int h2 = hash >> b->n_shift & YAK_BLK_MASK; + uint8_t *p = &b->b[y<<(YAK_BLK_SHIFT-3)]; + int i, z = h1, cnt = 0; + if ((h2&31) == 0) h2 = (h2 + 1) & YAK_BLK_MASK; // otherwise we may repeatedly use a few bits + for (i = 0; i < b->n_hashes; z = (z + h2) & YAK_BLK_MASK) { + uint8_t *q = &p[z>>3], u; + u = 1<<(z&7); + cnt += !!(*q & u); + *q |= u; + ++i; + } + return cnt; +} diff --git a/yak-count.cpp b/yak-count.cpp new file mode 100644 index 0000000..3dead7b --- /dev/null +++ b/yak-count.cpp @@ -0,0 +1,383 @@ +#include +#include "CommandLines.h" +#include "yak.h" +#include "khashl.h" + +#define YAK_COUNTER_BITS 10 +#define YAK_N_COUNTS (1<>YAK_COUNTER_BITS == (b)>>YAK_COUNTER_BITS) // lower 8 bits for counts; higher bits for k-mer +#define yak_ch_hash(a) ((a)>>YAK_COUNTER_BITS) +KHASHL_SET_INIT(static klib_unused, yak_ht_t, yak_ht, uint64_t, yak_ch_hash, yak_ch_eq) + +typedef struct { + int32_t bf_shift, bf_n_hash; + int32_t k; + int32_t pre; + int32_t n_thread; + int64_t chunk_size; +} yak_copt_t; + +typedef struct { + yak_ht_t *h; + yak_bf_t *b; +} yak_ch1_t; + +typedef struct { + int k, pre, n_hash, n_shift; + uint64_t tot; + yak_ch1_t *h; +} yak_ch_t; + +const unsigned char seq_nt4_table[256] = { // translate ACGT to 0123 + 0, 1, 2, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 0, 4, 1, 4, 4, 4, 2, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, + 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4 +}; + +#include +#include +#include +#include "kthread.h" + +/*** hash table ***/ + +static yak_ch_t *yak_ch_init(int k, int pre, int n_hash, int n_shift) +{ + yak_ch_t *h; + int i; + if (pre < YAK_COUNTER_BITS) return 0; + CALLOC(h, 1); + h->k = k, h->pre = pre; + CALLOC(h->h, 1<pre); + for (i = 0; i < 1<pre; ++i) + h->h[i].h = yak_ht_init(); + if (n_hash > 0 && n_shift > h->pre) { + h->n_hash = n_hash, h->n_shift = n_shift; + for (i = 0; i < 1<pre; ++i) + h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash); + } + return h; +} + +static void yak_ch_destroy_bf(yak_ch_t *h) +{ + int i; + for (i = 0; i < 1<pre; ++i) { + if (h->h[i].b) + yak_bf_destroy(h->h[i].b); + h->h[i].b = 0; + } +} + +static void yak_ch_destroy(yak_ch_t *h) +{ + int i; + if (h == 0) return; + yak_ch_destroy_bf(h); + for (i = 0; i < 1<pre; ++i) + yak_ht_destroy(h->h[i].h); + free(h->h); free(h); +} + +static int yak_ch_insert_list(yak_ch_t *h, int create_new, int n, const uint64_t *a) +{ + int j, mask = (1<pre) - 1, n_ins = 0; + yak_ch1_t *g; + if (n == 0) return 0; + g = &h->h[a[0]&mask]; + for (j = 0; j < n; ++j) { + int ins = 1, absent; + uint64_t x = a[j] >> h->pre; + khint_t k; + if ((a[j]&mask) != (a[0]&mask)) continue; + if (create_new) { + if (g->b) + ins = (yak_bf_insert(g->b, x) == h->n_hash); + if (ins) { + k = yak_ht_put(g->h, x<h, k)&YAK_MAX_COUNT) < YAK_MAX_COUNT) + ++kh_key(g->h, k); + } + } else { + k = yak_ht_get(g->h, x<h) && (kh_key(g->h, k)&YAK_MAX_COUNT) < YAK_MAX_COUNT) + ++kh_key(g->h, k); + } + } + return n_ins; +} + +static int yak_ch_get(const yak_ch_t *h, uint64_t x) +{ + int mask = (1<pre) - 1; + yak_ht_t *g = h->h[x&mask].h; + khint_t k; + k = yak_ht_get(g, x >> h->pre << YAK_COUNTER_BITS); + return k == kh_end(g)? -1 : kh_key(g, k)&YAK_MAX_COUNT; +} + +/*** generate histogram ***/ + +typedef struct { + uint64_t c[YAK_N_COUNTS]; +} buf_cnt_t; + +typedef struct { + const yak_ch_t *h; + buf_cnt_t *cnt; +} hist_aux_t; + +static void worker_hist(void *data, long i, int tid) // callback for kt_for() +{ + hist_aux_t *a = (hist_aux_t*)data; + uint64_t *cnt = a->cnt[tid].c; + yak_ht_t *g = a->h->h[i].h; + khint_t k; + for (k = 0; k < kh_end(g); ++k) + if (kh_exist(g, k)) + ++cnt[kh_key(g, k)&YAK_MAX_COUNT]; +} + +static void yak_ch_hist(const yak_ch_t *h, int64_t cnt[YAK_N_COUNTS], int n_thread) +{ + hist_aux_t a; + int i, j; + a.h = h; + memset(cnt, 0, YAK_N_COUNTS * sizeof(uint64_t)); + CALLOC(a.cnt, n_thread); + kt_for(n_thread, worker_hist, &a, 1<pre); + for (i = 0; i < YAK_N_COUNTS; ++i) cnt[i] = 0; + for (j = 0; j < n_thread; ++j) + for (i = 0; i < YAK_N_COUNTS; ++i) + cnt[i] += a.cnt[j].c[i]; + free(a.cnt); +} + +/*** shrink a hash table ***/ + +typedef struct { + int min, max; + yak_ch_t *h; +} shrink_aux_t; + +static void worker_shrink(void *data, long i, int tid) // callback for kt_for() +{ + shrink_aux_t *a = (shrink_aux_t*)data; + yak_ch_t *h = a->h; + yak_ht_t *g = h->h[i].h, *f; + khint_t k; + f = yak_ht_init(); + yak_ht_resize(f, kh_size(g)); + for (k = 0; k < kh_end(g); ++k) { + if (kh_exist(g, k)) { + int absent, c = kh_key(g, k) & YAK_MAX_COUNT; + if (c >= a->min && c <= a->max) + yak_ht_put(f, kh_key(g, k), &absent); + } + } + yak_ht_destroy(g); + h->h[i].h = f; +} + +static void yak_ch_shrink(yak_ch_t *h, int min, int max, int n_thread) +{ + int i; + shrink_aux_t a; + a.h = h, a.min = min, a.max = max; + kt_for(n_thread, worker_shrink, &a, 1<pre); + for (i = 0, h->tot = 0; i < 1<pre; ++i) + h->tot += kh_size(h->h[i].h); +} + +#include +#include +#include "kseq.h" // FASTA/Q parser +KSEQ_INIT(gzFile, gzread) + +void yak_copt_init(yak_copt_t *o) +{ + memset(o, 0, sizeof(yak_copt_t)); + o->bf_shift = 0; + o->bf_n_hash = 4; + o->k = 31; + o->pre = YAK_COUNTER_BITS; + o->n_thread = 4; + o->chunk_size = 10000000; +} + +typedef struct { + int n, m; + uint64_t n_ins; + uint64_t *a; +} ch_buf_t; + +static inline void ch_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k-mer $y to a linear buffer +{ + int pre = y & ((1<n == b->m) { + b->m = b->m < 8? 8 : b->m + (b->m>>1); + REALLOC(b->a, b->m); + } + b->a[b->n++] = y; +} + +static void count_seq_buf(ch_buf_t *buf, int k, int p, int len, const char *seq) // insert k-mers in $seq to linear buffer $buf +{ + int i, l; + 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; + if (++l >= k) + ch_insert_buf(buf, p, yak_hash_long(x)); + } else l = 0, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart + } +} + +typedef struct { // global data structure for kt_pipeline() + const yak_copt_t *opt; + int create_new; + kseq_t *ks; + yak_ch_t *h; +} pldat_t; + +typedef struct { // data structure for each step in kt_pipeline() + pldat_t *p; + int n, m, sum_len, nk; + int *len; + char **seq; + ch_buf_t *buf; +} stepdat_t; + +static void worker_for(void *data, long i, int tid) // callback for kt_for() +{ + stepdat_t *s = (stepdat_t*)data; + ch_buf_t *b = &s->buf[i]; + yak_ch_t *h = s->p->h; + b->n_ins += yak_ch_insert_list(h, s->p->create_new, b->n, b->a); +} + +static void *worker_count_all(void *data, int step, void *in) // callback for kt_pipeline() +{ + pldat_t *p = (pldat_t*)data; + if (step == 0) { // step 1: read a block of sequences + int ret; + stepdat_t *s; + CALLOC(s, 1); + s->p = p; + while ((ret = kseq_read(p->ks)) >= 0) { + int l = p->ks->seq.l; + if (l < p->opt->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); + } + MALLOC(s->seq[s->n], l); + memcpy(s->seq[s->n], p->ks->seq.s, l); + s->len[s->n++] = l; + s->sum_len += l; + s->nk += l - p->opt->k + 1; + if (s->sum_len >= p->opt->chunk_size) + break; + } + if (s->sum_len == 0) free(s); + else return s; + } else if (step == 1) { // step 2: extract k-mers + stepdat_t *s = (stepdat_t*)in; + int i, n = 1<opt->pre, m; + CALLOC(s->buf, n); + m = (int)(s->nk * 1.2 / n) + 1; + for (i = 0; i < n; ++i) { + s->buf[i].m = m; + MALLOC(s->buf[i].a, m); + } + for (i = 0; i < s->n; ++i) { + count_seq_buf(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]); + free(s->seq[i]); + } + free(s->seq); free(s->len); + return s; + } else if (step == 2) { // step 3: insert k-mers to hash table + stepdat_t *s = (stepdat_t*)in; + int i, n = 1<opt->pre; + uint64_t n_ins = 0; + kt_for(p->opt->n_thread, worker_for, s, n); + for (i = 0; i < n; ++i) { + n_ins += s->buf[i].n_ins; + free(s->buf[i].a); + } + p->h->tot += n_ins; + free(s->buf); + fprintf(stderr, "[M::%s::%.3f*%.2f] processed %d sequences; %ld distinct k-mers in the hash table\n", __func__, + yak_realtime(), yak_cputime() / yak_realtime(), s->n, (long)p->h->tot); + free(s); + } + return 0; +} + +yak_ch_t *yak_count(const char *fn, const yak_copt_t *opt, yak_ch_t *h0) +{ + pldat_t pl; + gzFile fp; + if ((fp = gzopen(fn, "r")) == 0) return 0; + pl.ks = kseq_init(fp); + pl.opt = opt; + if (h0) { + pl.h = h0, pl.create_new = 0; + assert(h0->k == opt->k && h0->pre == opt->pre); + } else { + pl.create_new = 1; + pl.h = yak_ch_init(opt->k, opt->pre, opt->bf_n_hash, opt->bf_shift); + } + kt_pipeline(3, worker_count_all, &pl, 3); + kseq_destroy(pl.ks); + gzclose(fp); + return pl.h; +} + +yak_ch_t *yak_count_file(const yak_copt_t *opt, int n_fn, char **fn) +{ + int i; + yak_ch_t *h = 0; + for (i = 0; i < n_fn; ++i) + h = yak_count(fn[i], opt, h); + if (opt->bf_shift > 0) + yak_ch_destroy_bf(h); + return h; +} + +void ha_count_high(const hifiasm_opt_t *asm_opt) +{ + int64_t cnt[YAK_N_COUNTS]; + yak_copt_t opt; + yak_ch_t *h; + yak_copt_init(&opt); + opt.k = asm_opt->k_mer_length; + opt.n_thread = asm_opt->thread_num; + opt.bf_shift = asm_opt->bf_shift; + h = yak_count_file(&opt, asm_opt->num_reads, asm_opt->read_file_names); + yak_ch_hist(h, cnt, opt.n_thread); + yak_ch_destroy(h); +} diff --git a/yak-sys.cpp b/yak-sys.cpp new file mode 100644 index 0000000..d01c68c --- /dev/null +++ b/yak-sys.cpp @@ -0,0 +1,43 @@ +#include +#include +#include "yak.h" + +int yak_verbose = 3; + +static double yak_realtime0; + +double yak_cputime(void) +{ + struct rusage r; + getrusage(RUSAGE_SELF, &r); + return r.ru_utime.tv_sec + r.ru_stime.tv_sec + 1e-6 * (r.ru_utime.tv_usec + r.ru_stime.tv_usec); +} + +static inline double yak_realtime_core(void) +{ + struct timeval tp; + struct timezone tzp; + gettimeofday(&tp, &tzp); + return tp.tv_sec + tp.tv_usec * 1e-6; +} + +void yak_reset_realtime(void) +{ + yak_realtime0 = yak_realtime_core(); +} + +double yak_realtime(void) +{ + return yak_realtime_core() - yak_realtime0; +} + +long yak_peakrss(void) +{ + struct rusage r; + getrusage(RUSAGE_SELF, &r); +#ifdef __linux__ + return r.ru_maxrss * 1024; +#else + return r.ru_maxrss; +#endif +} diff --git a/yak.h b/yak.h new file mode 100644 index 0000000..31529b7 --- /dev/null +++ b/yak.h @@ -0,0 +1,68 @@ +#ifndef __YAK_H__ +#define __YAK_H__ +#include + +#define YAK_BLK_SHIFT 9 // 64 bytes, the size of a cache line +#define YAK_BLK_MASK ((1<<(YAK_BLK_SHIFT)) - 1) + +typedef struct { + int n_shift, n_hashes; + uint8_t *b; +} yak_bf_t; + +extern const unsigned char seq_nt4_table[256]; + +void trio_partition(void); + +double yak_cputime(void); +void yak_reset_realtime(void); +double yak_realtime(void); +long yak_peakrss(void); + +yak_bf_t *yak_bf_init(int n_shift, int n_hashes); +void yak_bf_destroy(yak_bf_t *b); +int yak_bf_insert(yak_bf_t *b, uint64_t hash); + +static inline uint64_t yak_hash64(uint64_t key, uint64_t mask) // invertible integer hash function +{ + key = (~key + (key << 21)) & mask; // key = (key << 21) - key - 1; + key = key ^ key >> 24; + key = ((key + (key << 3)) + (key << 8)) & mask; // key * 265 + key = key ^ key >> 14; + key = ((key + (key << 2)) + (key << 4)) & mask; // key * 21 + key = key ^ key >> 28; + key = (key + (key << 31)) & mask; + return key; +} + +static inline uint64_t yak_hash64_64(uint64_t key) +{ + key = ~key + (key << 21); + key = key ^ key >> 24; + key = (key + (key << 3)) + (key << 8); + key = key ^ key >> 14; + key = (key + (key << 2)) + (key << 4); + key = key ^ key >> 28; + key = key + (key << 31); + return key; +} + +static inline uint64_t yak_hash_long(uint64_t x[4]) +{ + int j = x[1] < x[3]? 0 : 1; + return yak_hash64_64(x[j<<1|0]) + yak_hash64_64(x[j<<1|1]); +} + +#define CALLOC(ptr, len) ((ptr) = (__typeof__(ptr))calloc((len), sizeof(*(ptr)))) +#define MALLOC(ptr, len) ((ptr) = (__typeof__(ptr))malloc((len) * sizeof(*(ptr)))) +#define REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr)))) + +#ifndef klib_unused +#if (defined __clang__ && __clang_major__ >= 3) || (defined __GNUC__ && __GNUC__ >= 3) +#define klib_unused __attribute__ ((__unused__)) +#else +#define klib_unused +#endif +#endif /* klib_unused */ + +#endif // __YAK_H__