diff --git a/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch b/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch index 35b4a77..1c7d690 100644 Binary files a/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch and b/.vscode/ipch/4eee59091dc293a3/PROCESS_READ.ipch differ diff --git a/.vscode/ipch/756c62c8d5e4771d/COMMANDLINES.ipch b/.vscode/ipch/756c62c8d5e4771d/COMMANDLINES.ipch new file mode 100644 index 0000000..3caee36 Binary files /dev/null and b/.vscode/ipch/756c62c8d5e4771d/COMMANDLINES.ipch differ diff --git a/.vscode/ipch/856c69690c514632/KMER.ipch b/.vscode/ipch/856c69690c514632/KMER.ipch new file mode 100644 index 0000000..e65a7ca Binary files /dev/null and b/.vscode/ipch/856c69690c514632/KMER.ipch differ diff --git a/.vscode/ipch/856c69690c514632/mmap_address.bin b/.vscode/ipch/856c69690c514632/mmap_address.bin new file mode 100644 index 0000000..862b842 Binary files /dev/null and b/.vscode/ipch/856c69690c514632/mmap_address.bin differ diff --git a/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch b/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch new file mode 100644 index 0000000..71cb803 Binary files /dev/null and b/.vscode/ipch/b49c1e8b711f08aa/HASH_TABLE.ipch differ diff --git a/.vscode/ipch/b49c1e8b711f08aa/mmap_address.bin b/.vscode/ipch/b49c1e8b711f08aa/mmap_address.bin new file mode 100644 index 0000000..862b842 Binary files /dev/null and b/.vscode/ipch/b49c1e8b711f08aa/mmap_address.bin differ diff --git a/.vscode/ipch/b5c6b274f2611ee7/KMER.ipch b/.vscode/ipch/b5c6b274f2611ee7/KMER.ipch new file mode 100644 index 0000000..9221230 Binary files /dev/null and b/.vscode/ipch/b5c6b274f2611ee7/KMER.ipch differ diff --git a/.vscode/ipch/b5c6b274f2611ee7/mmap_address.bin b/.vscode/ipch/b5c6b274f2611ee7/mmap_address.bin new file mode 100644 index 0000000..862b842 Binary files /dev/null and b/.vscode/ipch/b5c6b274f2611ee7/mmap_address.bin differ diff --git a/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch b/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch index 2dd7dc6..62ad212 100644 Binary files a/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch and b/.vscode/ipch/c0e71cfe49f0fe81/ASSEMBLY.ipch differ diff --git a/.vscode/ipch/efd6dc00fdc3a0f/HASH_TABLE.ipch b/.vscode/ipch/efd6dc00fdc3a0f/HASH_TABLE.ipch new file mode 100644 index 0000000..5e161a1 Binary files /dev/null and b/.vscode/ipch/efd6dc00fdc3a0f/HASH_TABLE.ipch differ diff --git a/.vscode/ipch/efd6dc00fdc3a0f/mmap_address.bin b/.vscode/ipch/efd6dc00fdc3a0f/mmap_address.bin new file mode 100644 index 0000000..862b842 Binary files /dev/null and b/.vscode/ipch/efd6dc00fdc3a0f/mmap_address.bin differ diff --git a/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch b/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch index 456fe16..db8d684 100644 Binary files a/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch and b/.vscode/ipch/f2b98741f94c3e10/MAIN.ipch differ diff --git a/Assembly.cpp b/Assembly.cpp index 1a4489e..adbabb6 100644 --- a/Assembly.cpp +++ b/Assembly.cpp @@ -4,11 +4,15 @@ #include #include "Process_Read.h" #include "CommandLines.h" +#include "kmer.h" +#include "Hash_Table.h" + +Total_Count_Table TCB; - -void Counting() +/********************************for debug***************************************/ +void Verify_Counting() { @@ -17,46 +21,87 @@ void Counting() long long read_number = 0; + HPC_seq HPC_read; + uint64_t code; + Hash_code k_code; + int avalible_k = 0; + + fprintf(stdout, "Start Verifying ...\n"); + while (get_read(seq)) { - fprintf(stderr,"@%s\n",seq->name.s); - fprintf(stderr,"%s\n",seq->seq.s); - fprintf(stderr,"+\n"); - fprintf(stderr,"%s\n",seq->qual.s); + init_HPC_seq(&HPC_read, seq->seq.s, seq->seq.l); + init_Hash_code(&k_code); + + avalible_k = 0; + + while ((code = get_HPC_code(&HPC_read)) != 6) + { + if(code < 4) + { + k_mer_append(&k_code,code,k_mer_length); + avalible_k++; + if (avalible_k>=k_mer_length) + { + if(verify_Total_Count_Table(&TCB, &k_code, k_mer_length) == -1) + { + fprintf(stderr, "ERROR when subtracting!\n"); + } + } + + } + else + { + avalible_k = 0; + init_Hash_code(&k_code); + } + + } + + + read_number++; } fprintf(stdout, "read_number: %lld\n",read_number); + fprintf(stdout, "Start Traversing ...\n"); + + Traverse_Total_Count_Table(&TCB); + + fprintf(stdout, "Finish Traversing!\n"); } - - - - - - - - - - - +/********************************for debug***************************************/ void* Perform_Counting(void* arg) { + int thr_ID = *((int*)arg); + int i = 0; + HPC_seq HPC_read; + R_buffer_block curr_sub_block; init_R_buffer_block(&curr_sub_block); long long read_number = 0; - + long long k_mer_number = 0 ; + int file_flag = 1; + + uint64_t code; + + Hash_code k_code; + + int avalible_k = 0; + + while (file_flag != 0) { @@ -64,10 +109,46 @@ void* Perform_Counting(void* arg) read_number = read_number + curr_sub_block.num; + for (i = 0; i < curr_sub_block.num; i++) + { + + init_HPC_seq(&HPC_read, curr_sub_block.read[i].seq.s, curr_sub_block.read[i].seq.l); + init_Hash_code(&k_code); + + avalible_k = 0; + + while ((code = get_HPC_code(&HPC_read)) != 6) + { + if(code < 4) + { + k_mer_append(&k_code,code,k_mer_length); + avalible_k++; + if (avalible_k>=k_mer_length) + { + ///插入 + insert_Total_Count_Table(&TCB, &k_code, k_mer_length); + k_mer_number++; + } + + } + else + { + avalible_k = 0; + init_Hash_code(&k_code); + } + + } + + + } + + } - fprintf(stdout, "#########read_number: %lld\n",read_number); - fflush(stdout); + destory_R_buffer_block(&curr_sub_block); + free(arg); + + fprintf(stdout, "thr_ID: %d, read_number: %lld, k_mer_number: %lld\n",thr_ID, read_number, k_mer_number); } @@ -89,6 +170,15 @@ void* Perform_Counting(void* arg) void Counting_multiple_thr() { + + init_Total_Count_Table(k_mer_length, &TCB); + + fprintf(stdout, "TCB.prefix_bits: %d\n", TCB.prefix_bits); + fprintf(stdout, "TCB.suffix_bits: %d\n", TCB.suffix_bits); + fprintf(stdout, "TCB.size: %d\n", TCB.size); + + + pthread_t inputReadsHandle; init_R_buffer(thread_num); @@ -120,6 +210,11 @@ void Counting_multiple_thr() free(_r_threads); + destory_R_buffer(); + + + destory_Total_Count_Table(&TCB); + } @@ -134,6 +229,14 @@ void Counting_multiple_thr() + + + + + + + + diff --git a/Assembly.h b/Assembly.h index 66920a6..e7205de 100644 --- a/Assembly.h +++ b/Assembly.h @@ -3,7 +3,9 @@ -void Counting(); void Counting_multiple_thr(); +/********************************for debug***************************************/ +void Verify_Counting(); +/********************************for debug***************************************/ #endif diff --git a/CommandLines.cpp b/CommandLines.cpp index 54c19d1..5eb91fb 100644 --- a/CommandLines.cpp +++ b/CommandLines.cpp @@ -6,6 +6,7 @@ char* read_file_name = NULL; char* output_file_name = NULL; int thread_num = 1; +int k_mer_length = 40; void Print_H() { @@ -29,11 +30,12 @@ int CommandLine_process (int argc, char *argv[]) ketopt_t opt = KETOPT_INIT; int i, c; - while ((c = ketopt(&opt, argc, argv, 1, "ht:o:q:", longopts)) >= 0) { + while ((c = ketopt(&opt, argc, argv, 1, "ht:o:q:k:", longopts)) >= 0) { if (c == 100 || c == 'h') Print_H(); else if (c == 103 || c == 't') thread_num = atoi(opt.arg); else if (c == 102 || c == 'o') output_file_name = opt.arg; else if (c == 101 || c == 'q') read_file_name = opt.arg; + else if (c == 'k') k_mer_length = atoi(opt.arg); else if (c == '?') printf("unknown opt: -%c\n", opt.opt? opt.opt : ':'); else if (c == ':') printf("missing arg: -%c\n", opt.opt? opt.opt : ':'); } diff --git a/CommandLines.h b/CommandLines.h index cf430e1..f583693 100644 --- a/CommandLines.h +++ b/CommandLines.h @@ -7,6 +7,8 @@ extern char* read_file_name; extern char* output_file_name; extern int thread_num; +extern int k_mer_length; + int CommandLine_process (int argc, char *argv[]); diff --git a/Hash_Table.cpp b/Hash_Table.cpp new file mode 100644 index 0000000..ce5c2a5 --- /dev/null +++ b/Hash_Table.cpp @@ -0,0 +1,114 @@ +#include +#include +#include +#include "Hash_Table.h" + + +void init_Count_Table(Count_Table** table) +{ + *table = kh_init(COUNT64); +} + +void init_Total_Count_Table(int k, Total_Count_Table* TCB) +{ + if(k>64) + { + fprintf(stdout, "k-mer is too long. The length of k-mer must <= 64."); + fflush(stdout); + exit(0); + } + + int total_bits = k * 2; + TCB->prefix_bits = PREFIX_BITS; + TCB->suffix_bits = total_bits - TCB->prefix_bits; + if (TCB->suffix_bits > MAX_SUFFIX_BITS) + { + TCB->suffix_bits = MAX_SUFFIX_BITS; + TCB->prefix_bits = total_bits - TCB->suffix_bits; + } + ///TCB->suffix_mode = (1ULL<suffix_bits) - 1; + ///这个右移是安全的，因为TCB->suffix_bits不可能是0 + TCB->suffix_mode = ALL >> (64 - TCB->suffix_bits); + + ///number of small hash table + TCB->size = (1ULL<prefix_bits); + TCB->sub_h = (Count_Table**)malloc(sizeof(Count_Table*)*TCB->size); + TCB->sub_h_lock = (Hash_table_spin_lock*)malloc(sizeof(Hash_table_spin_lock)*TCB->size); + + int i = 0; + for (i = 0; i < TCB->size; i++) + { + init_Count_Table(&(TCB->sub_h[i])); + TCB->sub_h_lock[i].lock = 0; + } + + + +} + +void destory_Total_Count_Table(Total_Count_Table* TCB) +{ + int i; + for (i = 0; i < TCB->size; i++) + { + kh_destroy(COUNT64, TCB->sub_h[i]); + } + free(TCB->sub_h); + free(TCB->sub_h_lock); +} + + + + + +/********************************for debug***************************************/ +void test_COUNT64() +{ + + uint64_t sb[7] = {5000000000, 6000000000, 7000000000, 8000000000, 9000000000, 10000000000 ,5000000000}; + Count_Table* h; + ///构建哈希表 + init_Count_Table(&h); + + int i; + khint_t k; ///这就是个迭代器 + int absent; + + for (i = 0; i < 7; i++) + { + ///将5作为key插入到哈希表COUNT64中 + ///absent为0代表哈希表里面已经有这个key了 + k = kh_put(COUNT64, h, sb[i], &absent); + ///代表插入了一个哈希表中没有的元素 + ///则直接插入 + if (absent) + { + kh_value(h, k) = i; + } + else ///哈希表中已有的元素 + { + kh_value(h, k) = i; + } + } + + for (i = 0; i < 7; i++) + { + ///查询哈希表，key为k + k = kh_get(COUNT64, h, sb[i]); + + if (k != kh_end(h)) + { + fprintf(stderr, "value: %d\n", kh_value(h, k)); + } + else + { + fprintf(stderr, "not found!\n"); + } + } + + + kh_destroy(COUNT64, h); + + +} +/********************************for debug***************************************/ diff --git a/Hash_Table.h b/Hash_Table.h new file mode 100644 index 0000000..bd9e5da --- /dev/null +++ b/Hash_Table.h @@ -0,0 +1,189 @@ +#ifndef __HASHTABLE__ +#define __HASHTABLE__ +#include "khash.h" +#include "kmer.h" + +KHASH_MAP_INIT_INT64(COUNT64, int) + +typedef khash_t(COUNT64) Count_Table; + +#define PREFIX_BITS 16 +#define MAX_SUFFIX_BITS 64 + +typedef struct +{ + volatile int lock; + +}Hash_table_spin_lock; + +typedef struct +{ + Count_Table** sub_h; + Hash_table_spin_lock* sub_h_lock; + int prefix_bits; + int suffix_bits; + int size; + uint64_t suffix_mode; +} Total_Count_Table; + +/********************************for debug***************************************/ +inline void print_64bit(uint64_t x) +{ + int i; + for(i = 63; i >= 0; i--) + { + if(x & ((1ULL<x[0] | (code->x[1] << k); + low_key = code->x[0] | (code->x[1] << SAFE_SHIFT(k)); + //k不可能为0, 所以这个右移不会有问题 + h_key = code->x[1] >> (64 - k); + + ///注意suffix_bits最大就是64 + ///前一个右移不安全，因为TCB->suffix_bits有可能为64 + ///后一个左移安全，因为TCB->suffix_bits不可能为0 + //uint64_t sub_ID = (low_key >> TCB->suffix_bits) | (h_key << (64 - TCB->suffix_bits)); + uint64_t sub_ID = (low_key >> SAFE_SHIFT(TCB->suffix_bits)) | (h_key << (64 - TCB->suffix_bits)); + uint64_t sub_key = (low_key & TCB->suffix_mode); + + *get_sub_ID = sub_ID; + *get_sub_key = sub_key; + +} + + +inline void insert_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k) +{ + uint64_t sub_ID, sub_key; + get_sub_table(&sub_ID, &sub_key, TCB, code, k); + + khint_t t; ///这就是个迭代器 + int absent; + + + while (__sync_lock_test_and_set(&TCB->sub_h_lock[sub_ID].lock, 1)) + { + while (TCB->sub_h_lock[sub_ID].lock); + } + + t = kh_put(COUNT64, TCB->sub_h[sub_ID], sub_key, &absent); + if (absent) + { + kh_value(TCB->sub_h[sub_ID], t) = 1; + } + else ///哈希表中已有的元素 + { + //kh_value(TCB->sub_h[sub_ID], t) = kh_value(TCB->sub_h[sub_ID], t) + 1; + kh_value(TCB->sub_h[sub_ID], t)++; + } + + __sync_lock_release(&TCB->sub_h_lock[sub_ID].lock); +} + +inline int get_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k) +{ + + uint64_t sub_ID, sub_key; + get_sub_table(&sub_ID, &sub_key, TCB, code, k); + + khint_t t; ///这就是个迭代器 + int absent; + + ///查询哈希表，key为k + t = kh_get(COUNT64, TCB->sub_h[sub_ID], sub_key); + + if (t != kh_end(TCB->sub_h[sub_ID])) + { + return kh_value(TCB->sub_h[sub_ID], t); + } + else + { + return -1; + } + +} + +/********************************for debug***************************************/ +inline int verify_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k) +{ + uint64_t sub_ID, sub_key; + get_sub_table(&sub_ID, &sub_key, TCB, code, k); + + khint_t t; ///这就是个迭代器 + int absent; + + ///查询哈希表，key为k + t = kh_get(COUNT64, TCB->sub_h[sub_ID], sub_key); + + if (t != kh_end(TCB->sub_h[sub_ID])) + { + kh_value(TCB->sub_h[sub_ID], t)--; + if (kh_value(TCB->sub_h[sub_ID], t)<0) + { + return -1; + } + else + { + return 1; + } + } + else + { + return -1; + } +} +/********************************for debug***************************************/ + +/********************************for debug***************************************/ +inline int Traverse_Total_Count_Table(Total_Count_Table* TCB) +{ + int i; + Count_Table* h; + khint_t k; + + long long non_empty_k_mer = 0; + + for (i = 0; i < TCB->size; i++) + { + h = TCB->sub_h[i]; + for (k = kh_begin(h); k != kh_end(h); ++k) + { + if (kh_exist(h, k)) // test if a bucket contains data + { + non_empty_k_mer++; + + if (kh_value(h, k)!= 0) + { + fprintf(stderr, "ERROR when Traversing!\n"); + } + } + } + } + + fprintf(stdout, "non_empty_k_mer: %lld\n", non_empty_k_mer); +} +/********************************for debug***************************************/ + + +void init_Total_Count_Table(int k, Total_Count_Table* TCB); +void destory_Total_Count_Table(Total_Count_Table* TCB); + +void init_Count_Table(Count_Table** table); + +/********************************for debug***************************************/ +void test_COUNT64(); +/********************************for debug***************************************/ + +#endif \ No newline at end of file diff --git a/Makefile b/Makefile index dd44402..383a714 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ CC=g++ CFLAGS = -w -c -msse4.2 -mpopcnt -fomit-frame-pointer -Winline -O3 -lz LDFLAGS = -lm -lz -lpthread -O3 -mpopcnt -msse4.2 -lz -w -SOURCES = main.cpp CommandLines.cpp Process_Read.cpp Assembly.cpp +SOURCES = main.cpp CommandLines.cpp Process_Read.cpp Assembly.cpp kmer.cpp Hash_Table.cpp OBJECTS = $(SOURCES:.c=.o) EXECUTABLE = ccs_assembly diff --git a/Process_Read.cpp b/Process_Read.cpp index 2c266e4..32cc014 100644 --- a/Process_Read.cpp +++ b/Process_Read.cpp @@ -65,13 +65,8 @@ int get_read(kseq_t *s) void init_R_buffer_block(R_buffer_block* curr_sub_block) { - int i; - - for (i = 0; i < RDB.size; i++) - { - curr_sub_block->read = (kseq_t*)calloc(RDB.block_inner_size, sizeof(kseq_t)); - curr_sub_block->num = 0; - } + curr_sub_block->read = (kseq_t*)calloc(RDB.block_inner_size, sizeof(kseq_t)); + curr_sub_block->num = 0; } @@ -94,6 +89,27 @@ void init_R_buffer(int thread_num) } +void destory_R_buffer_block(R_buffer_block* curr_sub_block) +{ + + free(curr_sub_block->read); +} + + +void destory_R_buffer() +{ + int i = 0; + + for (i = 0; i < RDB.size; i++) + { + destory_R_buffer_block(&RDB.sub_block[i]); + } + + free(RDB.sub_block); + +} + + inline void load_read_block(R_buffer_block* read_batch, int batch_read_size, int* return_file_flag) { @@ -215,6 +231,8 @@ void* input_reads_muti_threads(void*) pthread_cond_signal(&i_readinputstallCond); //important pthread_mutex_unlock(&i_readinputMutex); + destory_R_buffer_block(&tmp_buf); + } diff --git a/Process_Read.h b/Process_Read.h index 76762d5..4734264 100644 --- a/Process_Read.h +++ b/Process_Read.h @@ -8,7 +8,7 @@ #include "kseq.h" #define READ_BLOCK_SIZE 64 -#define READ_BLOCK_NUM_PRE_THR 32 +#define READ_BLOCK_NUM_PRE_THR 100 #define IS_FULL(buffer) ((buffer.num >= buffer.size)?1:0) #define IS_EMPTY(buffer) ((buffer.num == 0)?1:0) @@ -46,5 +46,7 @@ void init_R_buffer_block(R_buffer_block* curr_sub_block); int get_reads_mul_thread(R_buffer_block* curr_sub_block); void Counting_block(); +void destory_R_buffer_block(R_buffer_block* curr_sub_block); +void destory_R_buffer(); #endif diff --git a/khash.h b/khash.h new file mode 100644 index 0000000..f75f347 --- /dev/null +++ b/khash.h @@ -0,0 +1,627 @@ +/* The MIT License + + Copyright (c) 2008, 2009, 2011 by Attractive Chaos + + Permission is hereby granted, free of charge, to any person obtaining + a copy of this software and associated documentation files (the + "Software"), to deal in the Software without restriction, including + without limitation the rights to use, copy, modify, merge, publish, + distribute, sublicense, and/or sell copies of the Software, and to + permit persons to whom the Software is furnished to do so, subject to + the following conditions: + + The above copyright notice and this permission notice shall be + included in all copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF + MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND + NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS + BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN + ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + SOFTWARE. +*/ + +/* + An example: + +#include "khash.h" +KHASH_MAP_INIT_INT(32, char) +int main() { + int ret, is_missing; + khiter_t k; + khash_t(32) *h = kh_init(32); + k = kh_put(32, h, 5, &ret); + kh_value(h, k) = 10; + k = kh_get(32, h, 10); + is_missing = (k == kh_end(h)); + k = kh_get(32, h, 5); + kh_del(32, h, k); + for (k = kh_begin(h); k != kh_end(h); ++k) + if (kh_exist(h, k)) kh_value(h, k) = 1; + kh_destroy(32, h); + return 0; +} +*/ + +/* + 2013-05-02 (0.2.8): + + * Use quadratic probing. When the capacity is power of 2, stepping function + i*(i+1)/2 guarantees to traverse each bucket. It is better than double + hashing on cache performance and is more robust than linear probing. + + In theory, double hashing should be more robust than quadratic probing. + However, my implementation is probably not for large hash tables, because + the second hash function is closely tied to the first hash function, + which reduce the effectiveness of double hashing. + + Reference: http://research.cs.vt.edu/AVresearch/hashing/quadratic.php + + 2011-12-29 (0.2.7): + + * Minor code clean up; no actual effect. + + 2011-09-16 (0.2.6): + + * The capacity is a power of 2. This seems to dramatically improve the + speed for simple keys. Thank Zilong Tan for the suggestion. Reference: + + - http://code.google.com/p/ulib/ + - http://nothings.org/computer/judy/ + + * Allow to optionally use linear probing which usually has better + performance for random input. Double hashing is still the default as it + is more robust to certain non-random input. + + * Added Wang's integer hash function (not used by default). This hash + function is more robust to certain non-random input. + + 2011-02-14 (0.2.5): + + * Allow to declare global functions. + + 2009-09-26 (0.2.4): + + * Improve portability + + 2008-09-19 (0.2.3): + + * Corrected the example + * Improved interfaces + + 2008-09-11 (0.2.2): + + * Improved speed a little in kh_put() + + 2008-09-10 (0.2.1): + + * Added kh_clear() + * Fixed a compiling error + + 2008-09-02 (0.2.0): + + * Changed to token concatenation which increases flexibility. + + 2008-08-31 (0.1.2): + + * Fixed a bug in kh_get(), which has not been tested previously. + + 2008-08-31 (0.1.1): + + * Added destructor +*/ + + +#ifndef __AC_KHASH_H +#define __AC_KHASH_H + +/*! + @header + + Generic hash table library. + */ + +#define AC_VERSION_KHASH_H "0.2.8" + +#include +#include +#include + +/* compiler specific configuration */ + +#if UINT_MAX == 0xffffffffu +typedef unsigned int khint32_t; +#elif ULONG_MAX == 0xffffffffu +typedef unsigned long khint32_t; +#endif + +#if ULONG_MAX == ULLONG_MAX +typedef unsigned long khint64_t; +#else +typedef unsigned long long khint64_t; +#endif + +#ifndef kh_inline +#ifdef _MSC_VER +#define kh_inline __inline +#else +#define kh_inline inline +#endif +#endif /* kh_inline */ + +#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 */ + +typedef khint32_t khint_t; +typedef khint_t khiter_t; + +#define __ac_isempty(flag, i) ((flag[i>>4]>>((i&0xfU)<<1))&2) +#define __ac_isdel(flag, i) ((flag[i>>4]>>((i&0xfU)<<1))&1) +#define __ac_iseither(flag, i) ((flag[i>>4]>>((i&0xfU)<<1))&3) +#define __ac_set_isdel_false(flag, i) (flag[i>>4]&=~(1ul<<((i&0xfU)<<1))) +#define __ac_set_isempty_false(flag, i) (flag[i>>4]&=~(2ul<<((i&0xfU)<<1))) +#define __ac_set_isboth_false(flag, i) (flag[i>>4]&=~(3ul<<((i&0xfU)<<1))) +#define __ac_set_isdel_true(flag, i) (flag[i>>4]|=1ul<<((i&0xfU)<<1)) + +#define __ac_fsize(m) ((m) < 16? 1 : (m)>>4) + +#ifndef kroundup32 +#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x)) +#endif + +#ifndef kcalloc +#define kcalloc(N,Z) calloc(N,Z) +#endif +#ifndef kmalloc +#define kmalloc(Z) malloc(Z) +#endif +#ifndef krealloc +#define krealloc(P,Z) realloc(P,Z) +#endif +#ifndef kfree +#define kfree(P) free(P) +#endif + +static const double __ac_HASH_UPPER = 0.77; + +#define __KHASH_TYPE(name, khkey_t, khval_t) \ + typedef struct kh_##name##_s { \ + khint_t n_buckets, size, n_occupied, upper_bound; \ + khint32_t *flags; \ + khkey_t *keys; \ + khval_t *vals; \ + } kh_##name##_t; + +#define __KHASH_PROTOTYPES(name, khkey_t, khval_t) \ + extern kh_##name##_t *kh_init_##name(void); \ + extern void kh_destroy_##name(kh_##name##_t *h); \ + extern void kh_clear_##name(kh_##name##_t *h); \ + extern khint_t kh_get_##name(const kh_##name##_t *h, khkey_t key); \ + extern int kh_resize_##name(kh_##name##_t *h, khint_t new_n_buckets); \ + extern khint_t kh_put_##name(kh_##name##_t *h, khkey_t key, int *ret); \ + extern void kh_del_##name(kh_##name##_t *h, khint_t x); + +#define __KHASH_IMPL(name, SCOPE, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) \ + SCOPE kh_##name##_t *kh_init_##name(void) { \ + return (kh_##name##_t*)kcalloc(1, sizeof(kh_##name##_t)); \ + } \ + SCOPE void kh_destroy_##name(kh_##name##_t *h) \ + { \ + if (h) { \ + kfree((void *)h->keys); kfree(h->flags); \ + kfree((void *)h->vals); \ + kfree(h); \ + } \ + } \ + SCOPE void kh_clear_##name(kh_##name##_t *h) \ + { \ + if (h && h->flags) { \ + memset(h->flags, 0xaa, __ac_fsize(h->n_buckets) * sizeof(khint32_t)); \ + h->size = h->n_occupied = 0; \ + } \ + } \ + SCOPE khint_t kh_get_##name(const kh_##name##_t *h, khkey_t key) \ + { \ + if (h->n_buckets) { \ + khint_t k, i, last, mask, step = 0; \ + mask = h->n_buckets - 1; \ + k = __hash_func(key); i = k & mask; \ + last = i; \ + while (!__ac_isempty(h->flags, i) && (__ac_isdel(h->flags, i) || !__hash_equal(h->keys[i], key))) { \ + i = (i + (++step)) & mask; \ + if (i == last) return h->n_buckets; \ + } \ + return __ac_iseither(h->flags, i)? h->n_buckets : i; \ + } else return 0; \ + } \ + SCOPE int kh_resize_##name(kh_##name##_t *h, khint_t new_n_buckets) \ + { /* This function uses 0.25*n_buckets bytes of working space instead of [sizeof(key_t+val_t)+.25]*n_buckets. */ \ + khint32_t *new_flags = 0; \ + khint_t j = 1; \ + { \ + kroundup32(new_n_buckets); \ + if (new_n_buckets < 4) new_n_buckets = 4; \ + if (h->size >= (khint_t)(new_n_buckets * __ac_HASH_UPPER + 0.5)) j = 0; /* requested size is too small */ \ + else { /* hash table size to be changed (shrink or expand); rehash */ \ + new_flags = (khint32_t*)kmalloc(__ac_fsize(new_n_buckets) * sizeof(khint32_t)); \ + if (!new_flags) return -1; \ + memset(new_flags, 0xaa, __ac_fsize(new_n_buckets) * sizeof(khint32_t)); \ + if (h->n_buckets < new_n_buckets) { /* expand */ \ + khkey_t *new_keys = (khkey_t*)krealloc((void *)h->keys, new_n_buckets * sizeof(khkey_t)); \ + if (!new_keys) { kfree(new_flags); return -1; } \ + h->keys = new_keys; \ + if (kh_is_map) { \ + khval_t *new_vals = (khval_t*)krealloc((void *)h->vals, new_n_buckets * sizeof(khval_t)); \ + if (!new_vals) { kfree(new_flags); return -1; } \ + h->vals = new_vals; \ + } \ + } /* otherwise shrink */ \ + } \ + } \ + if (j) { /* rehashing is needed */ \ + for (j = 0; j != h->n_buckets; ++j) { \ + if (__ac_iseither(h->flags, j) == 0) { \ + khkey_t key = h->keys[j]; \ + khval_t val; \ + khint_t new_mask; \ + new_mask = new_n_buckets - 1; \ + if (kh_is_map) val = h->vals[j]; \ + __ac_set_isdel_true(h->flags, j); \ + while (1) { /* kick-out process; sort of like in Cuckoo hashing */ \ + khint_t k, i, step = 0; \ + k = __hash_func(key); \ + i = k & new_mask; \ + while (!__ac_isempty(new_flags, i)) i = (i + (++step)) & new_mask; \ + __ac_set_isempty_false(new_flags, i); \ + if (i < h->n_buckets && __ac_iseither(h->flags, i) == 0) { /* kick out the existing element */ \ + { khkey_t tmp = h->keys[i]; h->keys[i] = key; key = tmp; } \ + if (kh_is_map) { khval_t tmp = h->vals[i]; h->vals[i] = val; val = tmp; } \ + __ac_set_isdel_true(h->flags, i); /* mark it as deleted in the old hash table */ \ + } else { /* write the element and jump out of the loop */ \ + h->keys[i] = key; \ + if (kh_is_map) h->vals[i] = val; \ + break; \ + } \ + } \ + } \ + } \ + if (h->n_buckets > new_n_buckets) { /* shrink the hash table */ \ + h->keys = (khkey_t*)krealloc((void *)h->keys, new_n_buckets * sizeof(khkey_t)); \ + if (kh_is_map) h->vals = (khval_t*)krealloc((void *)h->vals, new_n_buckets * sizeof(khval_t)); \ + } \ + kfree(h->flags); /* free the working space */ \ + h->flags = new_flags; \ + h->n_buckets = new_n_buckets; \ + h->n_occupied = h->size; \ + h->upper_bound = (khint_t)(h->n_buckets * __ac_HASH_UPPER + 0.5); \ + } \ + return 0; \ + } \ + SCOPE khint_t kh_put_##name(kh_##name##_t *h, khkey_t key, int *ret) \ + { \ + khint_t x; \ + if (h->n_occupied >= h->upper_bound) { /* update the hash table */ \ + if (h->n_buckets > (h->size<<1)) { \ + if (kh_resize_##name(h, h->n_buckets - 1) < 0) { /* clear "deleted" elements */ \ + *ret = -1; return h->n_buckets; \ + } \ + } else if (kh_resize_##name(h, h->n_buckets + 1) < 0) { /* expand the hash table */ \ + *ret = -1; return h->n_buckets; \ + } \ + } /* TODO: to implement automatically shrinking; resize() already support shrinking */ \ + { \ + khint_t k, i, site, last, mask = h->n_buckets - 1, step = 0; \ + x = site = h->n_buckets; k = __hash_func(key); i = k & mask; \ + if (__ac_isempty(h->flags, i)) x = i; /* for speed up */ \ + else { \ + last = i; \ + while (!__ac_isempty(h->flags, i) && (__ac_isdel(h->flags, i) || !__hash_equal(h->keys[i], key))) { \ + if (__ac_isdel(h->flags, i)) site = i; \ + i = (i + (++step)) & mask; \ + if (i == last) { x = site; break; } \ + } \ + if (x == h->n_buckets) { \ + if (__ac_isempty(h->flags, i) && site != h->n_buckets) x = site; \ + else x = i; \ + } \ + } \ + } \ + if (__ac_isempty(h->flags, x)) { /* not present at all */ \ + h->keys[x] = key; \ + __ac_set_isboth_false(h->flags, x); \ + ++h->size; ++h->n_occupied; \ + *ret = 1; \ + } else if (__ac_isdel(h->flags, x)) { /* deleted */ \ + h->keys[x] = key; \ + __ac_set_isboth_false(h->flags, x); \ + ++h->size; \ + *ret = 2; \ + } else *ret = 0; /* Don't touch h->keys[x] if present and not deleted */ \ + return x; \ + } \ + SCOPE void kh_del_##name(kh_##name##_t *h, khint_t x) \ + { \ + if (x != h->n_buckets && !__ac_iseither(h->flags, x)) { \ + __ac_set_isdel_true(h->flags, x); \ + --h->size; \ + } \ + } + +#define KHASH_DECLARE(name, khkey_t, khval_t) \ + __KHASH_TYPE(name, khkey_t, khval_t) \ + __KHASH_PROTOTYPES(name, khkey_t, khval_t) + +#define KHASH_INIT2(name, SCOPE, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) \ + __KHASH_TYPE(name, khkey_t, khval_t) \ + __KHASH_IMPL(name, SCOPE, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) + +#define KHASH_INIT(name, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) \ + KHASH_INIT2(name, static kh_inline klib_unused, khkey_t, khval_t, kh_is_map, __hash_func, __hash_equal) + +/* --- BEGIN OF HASH FUNCTIONS --- */ + +/*! @function + @abstract Integer hash function + @param key The integer [khint32_t] + @return The hash value [khint_t] + */ +#define kh_int_hash_func(key) (khint32_t)(key) +/*! @function + @abstract Integer comparison function + */ +#define kh_int_hash_equal(a, b) ((a) == (b)) +/*! @function + @abstract 64-bit integer hash function + @param key The integer [khint64_t] + @return The hash value [khint_t] + */ +#define kh_int64_hash_func(key) (khint32_t)((key)>>33^(key)^(key)<<11) +/*! @function + @abstract 64-bit integer comparison function + */ +#define kh_int64_hash_equal(a, b) ((a) == (b)) +/*! @function + @abstract const char* hash function + @param s Pointer to a null terminated string + @return The hash value + */ +static kh_inline khint_t __ac_X31_hash_string(const char *s) +{ + khint_t h = (khint_t)*s; + if (h) for (++s ; *s; ++s) h = (h << 5) - h + (khint_t)*s; + return h; +} +/*! @function + @abstract Another interface to const char* hash function + @param key Pointer to a null terminated string [const char*] + @return The hash value [khint_t] + */ +#define kh_str_hash_func(key) __ac_X31_hash_string(key) +/*! @function + @abstract Const char* comparison function + */ +#define kh_str_hash_equal(a, b) (strcmp(a, b) == 0) + +static kh_inline khint_t __ac_Wang_hash(khint_t key) +{ + key += ~(key << 15); + key ^= (key >> 10); + key += (key << 3); + key ^= (key >> 6); + key += ~(key << 11); + key ^= (key >> 16); + return key; +} +#define kh_int_hash_func2(key) __ac_Wang_hash((khint_t)key) + +/* --- END OF HASH FUNCTIONS --- */ + +/* Other convenient macros... */ + +/*! + @abstract Type of the hash table. + @param name Name of the hash table [symbol] + */ +#define khash_t(name) kh_##name##_t + +/*! @function + @abstract Initiate a hash table. + @param name Name of the hash table [symbol] + @return Pointer to the hash table [khash_t(name)*] + */ +#define kh_init(name) kh_init_##name() + +/*! @function + @abstract Destroy a hash table. + @param name Name of the hash table [symbol] + @param h Pointer to the hash table [khash_t(name)*] + */ +#define kh_destroy(name, h) kh_destroy_##name(h) + +/*! @function + @abstract Reset a hash table without deallocating memory. + @param name Name of the hash table [symbol] + @param h Pointer to the hash table [khash_t(name)*] + */ +#define kh_clear(name, h) kh_clear_##name(h) + +/*! @function + @abstract Resize a hash table. + @param name Name of the hash table [symbol] + @param h Pointer to the hash table [khash_t(name)*] + @param s New size [khint_t] + */ +#define kh_resize(name, h, s) kh_resize_##name(h, s) + +/*! @function + @abstract Insert a key to the hash table. + @param name Name of the hash table [symbol] + @param h Pointer to the hash table [khash_t(name)*] + @param k Key [type of keys] + @param r Extra return code: -1 if the operation failed; + 0 if the key is present in the hash table; + 1 if the bucket is empty (never used); 2 if the element in + the bucket has been deleted [int*] + @return Iterator to the inserted element [khint_t] + */ +#define kh_put(name, h, k, r) kh_put_##name(h, k, r) + +/*! @function + @abstract Retrieve a key from the hash table. + @param name Name of the hash table [symbol] + @param h Pointer to the hash table [khash_t(name)*] + @param k Key [type of keys] + @return Iterator to the found element, or kh_end(h) if the element is absent [khint_t] + */ +#define kh_get(name, h, k) kh_get_##name(h, k) + +/*! @function + @abstract Remove a key from the hash table. + @param name Name of the hash table [symbol] + @param h Pointer to the hash table [khash_t(name)*] + @param k Iterator to the element to be deleted [khint_t] + */ +#define kh_del(name, h, k) kh_del_##name(h, k) + +/*! @function + @abstract Test whether a bucket contains data. + @param h Pointer to the hash table [khash_t(name)*] + @param x Iterator to the bucket [khint_t] + @return 1 if containing data; 0 otherwise [int] + */ +#define kh_exist(h, x) (!__ac_iseither((h)->flags, (x))) + +/*! @function + @abstract Get key given an iterator + @param h Pointer to the hash table [khash_t(name)*] + @param x Iterator to the bucket [khint_t] + @return Key [type of keys] + */ +#define kh_key(h, x) ((h)->keys[x]) + +/*! @function + @abstract Get value given an iterator + @param h Pointer to the hash table [khash_t(name)*] + @param x Iterator to the bucket [khint_t] + @return Value [type of values] + @discussion For hash sets, calling this results in segfault. + */ +#define kh_val(h, x) ((h)->vals[x]) + +/*! @function + @abstract Alias of kh_val() + */ +#define kh_value(h, x) ((h)->vals[x]) + +/*! @function + @abstract Get the start iterator + @param h Pointer to the hash table [khash_t(name)*] + @return The start iterator [khint_t] + */ +#define kh_begin(h) (khint_t)(0) + +/*! @function + @abstract Get the end iterator + @param h Pointer to the hash table [khash_t(name)*] + @return The end iterator [khint_t] + */ +#define kh_end(h) ((h)->n_buckets) + +/*! @function + @abstract Get the number of elements in the hash table + @param h Pointer to the hash table [khash_t(name)*] + @return Number of elements in the hash table [khint_t] + */ +#define kh_size(h) ((h)->size) + +/*! @function + @abstract Get the number of buckets in the hash table + @param h Pointer to the hash table [khash_t(name)*] + @return Number of buckets in the hash table [khint_t] + */ +#define kh_n_buckets(h) ((h)->n_buckets) + +/*! @function + @abstract Iterate over the entries in the hash table + @param h Pointer to the hash table [khash_t(name)*] + @param kvar Variable to which key will be assigned + @param vvar Variable to which value will be assigned + @param code Block of code to execute + */ +#define kh_foreach(h, kvar, vvar, code) { khint_t __i; \ + for (__i = kh_begin(h); __i != kh_end(h); ++__i) { \ + if (!kh_exist(h,__i)) continue; \ + (kvar) = kh_key(h,__i); \ + (vvar) = kh_val(h,__i); \ + code; \ + } } + +/*! @function + @abstract Iterate over the values in the hash table + @param h Pointer to the hash table [khash_t(name)*] + @param vvar Variable to which value will be assigned + @param code Block of code to execute + */ +#define kh_foreach_value(h, vvar, code) { khint_t __i; \ + for (__i = kh_begin(h); __i != kh_end(h); ++__i) { \ + if (!kh_exist(h,__i)) continue; \ + (vvar) = kh_val(h,__i); \ + code; \ + } } + +/* More convenient interfaces */ + +/*! @function + @abstract Instantiate a hash set containing integer keys + @param name Name of the hash table [symbol] + */ +#define KHASH_SET_INIT_INT(name) \ + KHASH_INIT(name, khint32_t, char, 0, kh_int_hash_func, kh_int_hash_equal) + +/*! @function + @abstract Instantiate a hash map containing integer keys + @param name Name of the hash table [symbol] + @param khval_t Type of values [type] + */ +#define KHASH_MAP_INIT_INT(name, khval_t) \ + KHASH_INIT(name, khint32_t, khval_t, 1, kh_int_hash_func, kh_int_hash_equal) + +/*! @function + @abstract Instantiate a hash set containing 64-bit integer keys + @param name Name of the hash table [symbol] + */ +#define KHASH_SET_INIT_INT64(name) \ + KHASH_INIT(name, khint64_t, char, 0, kh_int64_hash_func, kh_int64_hash_equal) + +/*! @function + @abstract Instantiate a hash map containing 64-bit integer keys + @param name Name of the hash table [symbol] + @param khval_t Type of values [type] + */ +#define KHASH_MAP_INIT_INT64(name, khval_t) \ + KHASH_INIT(name, khint64_t, khval_t, 1, kh_int64_hash_func, kh_int64_hash_equal) + +typedef const char *kh_cstr_t; +/*! @function + @abstract Instantiate a hash map containing const char* keys + @param name Name of the hash table [symbol] + */ +#define KHASH_SET_INIT_STR(name) \ + KHASH_INIT(name, kh_cstr_t, char, 0, kh_str_hash_func, kh_str_hash_equal) + +/*! @function + @abstract Instantiate a hash map containing const char* keys + @param name Name of the hash table [symbol] + @param khval_t Type of values [type] + */ +#define KHASH_MAP_INIT_STR(name, khval_t) \ + KHASH_INIT(name, kh_cstr_t, khval_t, 1, kh_str_hash_func, kh_str_hash_equal) + +#endif /* __AC_KHASH_H */ diff --git a/kmer.cpp b/kmer.cpp new file mode 100644 index 0000000..ec3b380 --- /dev/null +++ b/kmer.cpp @@ -0,0 +1,16 @@ +#include +#include +#include "kmer.h" + +void init_HPC_seq(HPC_seq* seq, char* str, long long l) +{ + seq->i = 0; + seq->l = l; + seq->str = str; +} + +void init_Hash_code(Hash_code* code) +{ + code->x[0] = 0; + code->x[1] = 0; +} \ No newline at end of file diff --git a/kmer.h b/kmer.h new file mode 100644 index 0000000..55e35a3 --- /dev/null +++ b/kmer.h @@ -0,0 +1,89 @@ +#ifndef __KMER__ +#define __KMER__ +#include "Process_Read.h" + +#define ALL (0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff) + +#define SAFE_SHIFT(k) k & ((k < 64)?ALL:0) + + +static unsigned char seq_nt6_table[256] = { + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 0, 5, 1, 5, 5, 5, 2, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 0, 5, 1, 5, 5, 5, 2, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, + 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5 +}; + + + +typedef struct +{ + ///最大64-mer + ///x[0]低位 + ///x[1]高位 + uint64_t x[2]; + +} Hash_code; + + + +typedef struct +{ + char* str; + long long l; + long long i; + +} HPC_seq; + + +inline uint64_t get_HPC_code(HPC_seq* seq) +{ + + if(seq->i < seq ->l) + { + char code = seq->str[seq->i]; + + for (; seq->i < seq->l; seq->i++) + { + if (seq->str[seq->i] != code) + { + break; + } + } + + return (uint64_t)seq_nt6_table[(uint8_t)code]; + } + else + { + ///end + return 6; + } + +} + +inline void k_mer_append(Hash_code* code, uint64_t c, int k) +{ + + uint64_t mask = ALL >> (64 -k); + + code->x[0] = ((code->x[0]<<1) | (c&1)) & mask; + code->x[1] = ((code->x[1]<<1) | (c>>1)) & mask; +} + +void init_HPC_seq(HPC_seq* seq, char* str, long long l); +void init_Hash_code(Hash_code* code); + + +#endif \ No newline at end of file diff --git a/main.cpp b/main.cpp index 3484607..c704208 100644 --- a/main.cpp +++ b/main.cpp @@ -4,6 +4,17 @@ #include "Process_Read.h" #include "Assembly.h" +/********************************for debug***************************************/ +///使用这个函数的时候，必须把Counting_multiple_thr()里的destory_Total_Count_Table(&TCB)注释掉 +void debug_Counting() +{ + init_kseq(read_file_name); + Verify_Counting(); + fprintf(stderr, "debug over!\n"); + destory_kseq(); +} +/********************************for debug***************************************/ + int main(int argc, char *argv[]) { @@ -13,15 +24,15 @@ int main(int argc, char *argv[]) init_kseq(read_file_name); - + fprintf(stdout, "k-mer length: %d\n",k_mer_length); ///Counting(); Counting_multiple_thr(); - - destory_kseq(); + ///debug_Counting(); + return 1; }