xxay/hifiasm
1
0
Fork 0
mirror of https://github.com/chhylp123/hifiasm.git synced 2026-10-01 12:08:11 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' of github.com:chhylp123/hifiasm

Browse Source
This commit is contained in:
Heng Li
2020-04-06 08:25:21 -04:00
parent c57b63653d 87fc103e63
commit 7f580850e8
30 changed files with 2209 additions and 4602 deletions
Download Patch File Download Diff File Expand all files Collapse all files
Assembly.cpp
+301 -1005
View File
File diff suppressed because it is too large Load Diff
Assembly.h
+2 -4
View File
@@ -8,8 +8,6 @@
#define Get_Cigar_Type(RECORD) (RECORD&3)
#define Get_Cigar_Length(RECORD) (RECORD>>2)
void Counting_multiple_thr();
void Build_hash_table_multiple_thr();
void Overlap_calculate_multipe_thr();
void Correct_Reads(int last_round);
int ha_assemble(void);
#endif
CommandLines.cpp
+32 -14
View File
@@ -5,11 +5,16 @@
#include "ketopt.h"
#include <sys/time.h>
#define VERSION "0.3.0"
#define DEFAULT_OUTPUT "hifiasm.asm"
hifiasm_opt_t asm_opt;
static ko_longopt_t long_options[] = {
{ "version", ko_no_argument, 300 },
{ "dbg-gfa", ko_no_argument, 301 },
{ 0, 0, 0 }
};
double Get_T(void)
{
struct timeval t;
@@ -23,12 +28,14 @@ void Print_H(hifiasm_opt_t* asm_opt)
fprintf(stderr, "Options:\n");
fprintf(stderr, " Assembly:\n");
fprintf(stderr, " -o FILE prefix of output files [%s]\n", asm_opt->output_file_name);
///fprintf(stderr, " -c FILE file including trio information\n");
fprintf(stderr, " -t INT number of threads [%d]\n", asm_opt->thread_num);
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, " -w write all overlaps to disk, can accelerate assembly next time [%d]\n", asm_opt->write_index_to_disk);
fprintf(stderr, " -k INT k-mer length (must be <64) [%d]\n", asm_opt->k_mer_length);
fprintf(stderr, " -w INT minimizer window size [%d]\n", asm_opt->mz_win);
fprintf(stderr, " -f INT number of bits for bloom filter [%d]\n", asm_opt->bf_shift);
fprintf(stderr, " -D FLOAT drop k-mers occuring >FLOAT*coverage times [%.1f]\n", asm_opt->high_factor);
fprintf(stderr, " -N INT consider up to INT overlaps for each oriented read [%d]\n", asm_opt->max_n_chain);
///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");
fprintf(stderr, " -i ignore saved overlaps in *.ovlp* files\n");
@@ -63,7 +70,13 @@ 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 = 40;
asm_opt->k_mer_length = 51;
asm_opt->mz_win = 51;
asm_opt->bf_shift = 37;
asm_opt->high_factor = 5.0f;
asm_opt->no_HPC = 0;
asm_opt->no_kmer_flt = 0;
asm_opt->max_n_chain = 400;
asm_opt->k_mer_min_freq = 3;
asm_opt->k_mer_max_freq = 66;
asm_opt->load_index_from_disk = 1;
@@ -71,7 +84,6 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->number_of_round = 2;
asm_opt->adapterLen = 0;
asm_opt->clean_round = 4;
asm_opt->complete_threads = 0;
asm_opt->small_pop_bubble_size = 100000;
asm_opt->large_pop_bubble_size = 10000000;
asm_opt->min_drop_rate = 0.2;
@@ -84,6 +96,7 @@ void init_opt(hifiasm_opt_t* asm_opt)
asm_opt->max_short_tip = 3;
asm_opt->min_cnt = 2;
asm_opt->mid_cnt = 5;
asm_opt->verbose_gfa = 0;
}
void destory_opt(hifiasm_opt_t* asm_opt)
@@ -94,13 +107,13 @@ void destory_opt(hifiasm_opt_t* asm_opt)
}
}
void clear_opt(hifiasm_opt_t* asm_opt, int last_round)
void clear_opt(hifiasm_opt_t* asm_opt, int round)
{
asm_opt->complete_threads = 0;
asm_opt->num_bases = 0;
asm_opt->num_corrected_bases = 0;
asm_opt->num_recorrected_bases = 0;
asm_opt->roundID = asm_opt->number_of_round - last_round;
asm_opt->mem_buf = 0;
asm_opt->roundID = round;
}
int check_file(char* name, const char* opt)
@@ -295,24 +308,28 @@ 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:lw:m:n:r:a:b:z:x:y:p:c:d:M:P:if:D:FN:", long_options)) >= 0) {
if (c == 'h')
{
Print_H(asm_opt);
return 0;
}
else if (c == 'v')
else if (c == 'v' || c == 300)
{
fprintf(stderr, "[Version] %s\n", VERSION);
puts(HA_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);
else if (c == 'k') asm_opt->k_mer_length = atoi(opt.arg);
else if (c == 'i') asm_opt->load_index_from_disk = 0;
else if (c == 'l') asm_opt->load_index_from_disk = 1;
else if (c == 'w') asm_opt->write_index_to_disk = 1;
else if (c == 'w') asm_opt->mz_win = atoi(opt.arg);
else if (c == 'D') asm_opt->high_factor = atof(opt.arg);
else if (c == 'F') asm_opt->no_kmer_flt = 1;
else if (c == 'N') asm_opt->max_n_chain = atoi(opt.arg);
else if (c == 'a') asm_opt->clean_round = atoi(opt.arg);
else if (c == 'z') asm_opt->adapterLen = atoi(opt.arg);
else if (c == 'b') asm_opt->required_read_name = opt.arg;
@@ -325,6 +342,7 @@ int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt)
else if (c == 'p') asm_opt->small_pop_bubble_size = atoll(opt.arg);
else if (c == 'm') asm_opt->large_pop_bubble_size = atoll(opt.arg);
else if (c == 'n') asm_opt->max_short_tip = atoll(opt.arg);
else if (c == 301) asm_opt->verbose_gfa = 1;
else if (c == ':')
{
fprintf(stderr, "[ERROR] missing option argument in \"%s\"\n", argv[opt.i - 1]);
CommandLines.h
+13 -5
View File
@@ -3,8 +3,9 @@
#include <pthread.h>
#define HA_VERSION "0.3.0-dirty-r199"
#define VERBOSE 0
#define VERBOSE_GFA 1
typedef struct {
int num_reads;
@@ -15,6 +16,12 @@ typedef struct {
char* mat_index;
int thread_num;
int k_mer_length;
int mz_win;
int bf_shift;
float high_factor;
int no_HPC;
int no_kmer_flt;
int max_n_chain;
int k_mer_min_freq;
int k_mer_max_freq;
int load_index_from_disk;
@@ -22,7 +29,6 @@ typedef struct {
int number_of_round;
int adapterLen;
int clean_round;
int complete_threads;
int roundID;
int max_hang_Len;
int gap_fuzz;
@@ -31,6 +37,7 @@ typedef struct {
int max_short_tip;
int min_cnt;
int mid_cnt;
int verbose_gfa;
float max_hang_rate;
float min_drop_rate;
@@ -41,6 +48,7 @@ typedef struct {
long long num_bases;
long long num_corrected_bases;
long long num_recorrected_bases;
long long mem_buf;
long long coverage;
} hifiasm_opt_t;
@@ -48,8 +56,8 @@ extern hifiasm_opt_t asm_opt;
void init_opt(hifiasm_opt_t* asm_opt);
void destory_opt(hifiasm_opt_t* asm_opt);
void clear_opt(hifiasm_opt_t* asm_opt, int last_round);
int CommandLine_process (int argc, char *argv[], hifiasm_opt_t* asm_opt);
void clear_opt(hifiasm_opt_t* asm_opt, int round);
int CommandLine_process(int argc, char *argv[], hifiasm_opt_t* asm_opt);
double Get_T(void);
#endif
#endif
Correct.cpp
-3
View File
@@ -7213,8 +7213,6 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF,
haplotype_evdience_alloc* hap, Round2_alignment* second_round,
int force_repeat, int is_consensus, int* fully_cov, int* abnormal)
{
reverse_complement(g_read->seq, g_read->length);
clear_Correct_dumy(dumy, overlap_list);
long long window_start, window_end;
@@ -7261,7 +7259,6 @@ void correct_overlap(overlap_region_alloc* overlap_list, All_reads* R_INF,
(*fully_cov) = check_if_fully_covered(overlap_list, R_INF, g_read, dumy, g, abnormal);
}
Correct.h
+1 -1
View File
@@ -1173,4 +1173,4 @@ long long* max_t_pos, long long* max_q_pos, long long* score, long long* droped)
#define GAP_EXT_KSW 2
#define Z_DROP_KSW 400
#define BAND_KSW 50
#endif
#endif
Hash_Table.cpp
+116 -1107
View File
File diff suppressed because it is too large Load Diff
Hash_Table.h
+43 -400
View File
@@ -1,13 +1,6 @@
#ifndef __HASHTABLE__
#define __HASHTABLE__
#include "khash.h"
#include "kmer.h"
KHASH_MAP_INIT_INT64(COUNT64, int)
typedef khash_t(COUNT64) Count_Table;
KHASH_MAP_INIT_INT64(POS64, uint64_t)
typedef khash_t(POS64) Pos_Table;
#include "htab.h"
#define PREFIX_BITS 16
#define MAX_SUFFIX_BITS 64
@@ -37,26 +30,8 @@ typedef khash_t(POS64) Pos_Table;
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;
///number of subtable
int size;
uint64_t suffix_mode;
uint64_t non_unique_k_mer;
} Total_Count_Table;
typedef struct
{
uint64_t offset;
uint64_t readID;
uint32_t offset;
uint32_t readID:31, rev:1;
} k_mer_pos;
typedef struct
@@ -70,17 +45,9 @@ typedef struct
typedef struct
{
k_mer_pos_list* list;
uint64_t size;
uint64_t length;
} k_mer_pos_list_alloc;
typedef struct
{
int C_L[CIGAR_MAX_LENGTH];
char C_C[CIGAR_MAX_LENGTH];
int length;
int C_L[CIGAR_MAX_LENGTH];
char C_C[CIGAR_MAX_LENGTH];
int length;
} CIGAR;
typedef struct
@@ -97,407 +64,96 @@ typedef struct
CIGAR cigar;
} window_list;
typedef struct
{
window_list* buffer;
long long length;
long long size;
}window_list_alloc;
int32_t length;
int32_t size;
} window_list_alloc;
typedef struct
{
uint64_t* buffer;
uint64_t length;
uint64_t size;
}Fake_Cigar;
uint32_t length;
uint32_t size;
} Fake_Cigar;
typedef struct
{
uint64_t x_id;
uint32_t x_id;
///the begining and end of the whole overlap
uint64_t x_pos_s;
uint64_t x_pos_e;
uint64_t x_pos_strand;
uint32_t x_pos_s;
uint32_t x_pos_e;
uint32_t x_pos_strand;
uint64_t y_id;
uint64_t y_pos_s;
uint64_t y_pos_e;
uint64_t y_pos_strand;
uint32_t y_id;
uint32_t y_pos_s;
uint32_t y_pos_e;
uint32_t y_pos_strand;
uint64_t overlapLen;
uint64_t shared_seed;
uint64_t align_length;
///uint64_t total_errors;
uint32_t overlapLen;
uint32_t shared_seed;
uint32_t align_length;
uint8_t is_match;
uint8_t without_large_indel;
uint64_t non_homopolymer_errors;
int8_t strong;
uint32_t non_homopolymer_errors;
window_list* w_list;
uint64_t w_list_size;
uint64_t w_list_length;
int8_t strong;
uint32_t w_list_size;
uint32_t w_list_length;
Fake_Cigar f_cigar;
window_list_alloc boundary_cigars;
} overlap_region;
typedef struct
{
overlap_region* list;
uint64_t size;
uint64_t length;
///uint64_t mapped_overlaps_length;
long long mapped_overlaps_length;
int64_t mapped_overlaps_length;
} overlap_region_alloc;
typedef struct
{
///uint64_t offset;
long long offset;
///uint64_t self_offset;
long long self_offset;
uint64_t readID;
uint8_t strand;
uint32_t readID:31, strand:1;
uint32_t offset, self_offset;
} k_mer_hit;
typedef struct
{
k_mer_hit node;
uint64_t ID;
} ElemType;
typedef struct
{
ElemType* heap;
uint64_t* index_i;
int len;
int MaxSize;
} HeapSq;
typedef struct
{
long long* score;
long long* pre;
long long* indels;
long long* self_length;
long long length;
long long size;
typedef struct {
int32_t *score;
int64_t *pre;
int32_t *indels;
int32_t *self_length;
int64_t *tmp; // MUST BE 64-bit integer
int64_t length;
int64_t size;
} Chain_Data;
typedef struct
{
k_mer_hit* list;
k_mer_hit* tmp;
long long length;
long long size;
uint64_t foward_pos;
uint64_t rc_pos;
Chain_Data chainDP;
} Candidates_list;
typedef struct
{
Pos_Table** sub_h;
Hash_table_spin_lock* sub_h_lock;
int prefix_bits;
int suffix_bits;
///number of subtable
int size;
uint64_t suffix_mode;
k_mer_pos* pos;
uint64_t useful_k_mer;
uint64_t total_occ;
uint64_t* k_mer_index;
} Total_Pos_Table;
inline uint64_t mod_d(uint64_t h_key, uint64_t low_key, uint64_t d)
{
uint64_t result = (h_key >> 32) % d;
result = ((result << 32) + (h_key & (uint64_t)0xffffffff)) % d;
result = ((result << 32) + (low_key >> 32)) % d;
result = ((result << 32) + (low_key & (uint64_t)0xffffffff)) % d;
return result;
}
////suffix_bits = 64 in default
inline int recover_hash_code(uint64_t sub_ID, uint64_t sub_key, Hash_code* code,
uint64_t suffix_mode, int suffix_bits, int k)
{
uint64_t h_key, low_key;
h_key = low_key = 0;
low_key = sub_ID << SAFE_SHIFT(suffix_bits);
low_key = low_key | sub_key;
h_key = sub_ID >> (64 - suffix_bits);
code->x[0] = code->x[1] = 0;
uint64_t mask = ALL >> (64 - k);
code->x[0] = low_key & mask;
code->x[1] = h_key << (64 - k);
code->x[1] = code->x[1] | (low_key >> SAFE_SHIFT(k));
return 1;
}
///inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, Total_Count_Table* TCB, Hash_code* code, int k)
inline int get_sub_table(uint64_t* get_sub_ID, uint64_t* get_sub_key, uint64_t suffix_mode, int suffix_bits,
Hash_code* code, int k)
{
uint64_t h_key, low_key;
///k might be 64，so it is unsafe
///low_key = code->x[0] | (code->x[1] << k);
low_key = code->x[0] | (code->x[1] << SAFE_SHIFT(k));
//k cannot be 0, so this shift is safe
h_key = code->x[1] >> (64 - k);
if(mod_d(h_key, low_key, MODE_VALUE) > 3)
{
return 0;
}
uint64_t sub_ID = (low_key >> SAFE_SHIFT(suffix_bits)) | (h_key << (64 - suffix_bits));
uint64_t sub_key = (low_key & suffix_mode);
*get_sub_ID = sub_ID;
*get_sub_key = sub_key;
return 1;
}
inline int insert_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k)
{
uint64_t sub_ID, sub_key;
if(!get_sub_table(&sub_ID, &sub_key, TCB->suffix_mode, TCB->suffix_bits, code, k))
{
return 0;
}
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);
return 1;
}
inline int get_Total_Count_Table(Total_Count_Table* TCB, Hash_code* code, int k)
{
uint64_t sub_ID, sub_key;
if(!get_sub_table(&sub_ID, &sub_key, TCB->suffix_mode, TCB->suffix_bits, code, k))
{
return 0;
}
khint_t t;
///query hash table，key is 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 0;
}
}
inline uint64_t get_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k, uint64_t* r_sub_ID)
{
uint64_t sub_ID, sub_key;
if(!get_sub_table(&sub_ID, &sub_key, PCB->suffix_mode, PCB->suffix_bits, code, k))
{
return (uint64_t)-1;
}
khint_t t;
///query hash table，key is k
t = kh_get(POS64, PCB->sub_h[sub_ID], sub_key);
if (t != kh_end(PCB->sub_h[sub_ID]))
{
*r_sub_ID = sub_ID;
return kh_value(PCB->sub_h[sub_ID], t);
}
else
{
return (uint64_t)-1;
}
}
inline uint64_t count_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k)
{
uint64_t sub_ID;
uint64_t ret = get_Total_Pos_Table(PCB, code, k, &sub_ID);
if(ret != (uint64_t)-1)
{
return PCB->k_mer_index[ret + 1] - PCB->k_mer_index[ret];
}
else
{
return 0;
}
}
inline uint64_t locate_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, k_mer_pos** list, int k, uint64_t* r_sub_ID)
{
uint64_t ret = get_Total_Pos_Table(PCB, code, k, r_sub_ID);
if(ret != (uint64_t)-1)
{
*list = PCB->k_mer_index[ret] + PCB->pos;
return PCB->k_mer_index[ret + 1] - PCB->k_mer_index[ret];
}
else
{
*list = NULL;
return 0;
}
}
int cmp_k_mer_pos(const void * a, const void * b);
inline uint64_t insert_Total_Pos_Table(Total_Pos_Table* PCB, Hash_code* code, int k, uint64_t readID, uint64_t pos)
{
k_mer_pos* list;
int flag = 0;
uint64_t sub_ID;
uint64_t occ = locate_Total_Pos_Table(PCB, code, &list, k, &sub_ID);
if (occ)
{
while (__sync_lock_test_and_set(&PCB->sub_h_lock[sub_ID].lock, 1))
{
while (PCB->sub_h_lock[sub_ID].lock);
}
if (list[0].offset + 1 < occ)
{
list[0].offset++;
list[list[0].offset].readID = readID;
///list[list[0].offset].readID = readID|direction;
list[list[0].offset].offset = pos;
}
else
{
list[0].readID = readID;
///list[0].readID = readID|direction;
list[0].offset = pos;
flag = 1;
}
__sync_lock_release(&PCB->sub_h_lock[sub_ID].lock);
//if all pos has been saved, it is safe to sort
if (flag && occ>1)
{
qsort(list, occ, sizeof(k_mer_pos), cmp_k_mer_pos);
}
return 1;
}
else
{
return 0;
}
}
void init_Total_Count_Table(int k, Total_Count_Table* TCB);
void init_Total_Pos_Table(Total_Pos_Table* TCB, Total_Count_Table* pre_TCB);
void destory_Total_Count_Table(Total_Count_Table* TCB);
void init_Count_Table(Count_Table** table);
void init_Pos_Table(Count_Table** pre_table, Pos_Table** table);
void destory_Total_Pos_Table(Total_Pos_Table* TCB);
void write_Total_Pos_Table(Total_Pos_Table* TCB, char* read_file_name);
int load_Total_Pos_Table(Total_Pos_Table* TCB, char* read_file_name);
void Traverse_Counting_Table(Total_Count_Table* TCB, Total_Pos_Table* PCB, int k_mer_min_freq, int k_mer_max_freq);
void init_Candidates_list(Candidates_list* l);
void clear_Candidates_list(Candidates_list* l);
void destory_Candidates_list(Candidates_list* l);
void init_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list);
void destory_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list);
void clear_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list);
void append_k_mer_pos_list_alloc(k_mer_pos_list_alloc* list, k_mer_pos* n_list, uint64_t n_length,
uint64_t n_end_pos, uint8_t n_direction);
void merge_k_mer_pos_list_alloc_heap_sort(k_mer_pos_list_alloc* list, Candidates_list* candidates, HeapSq* HBT);
void Init_Heap(HeapSq* HBT);
void destory_Heap(HeapSq* HBT);
void clear_Heap(HeapSq* HBT);
void init_overlap_region_alloc(overlap_region_alloc* list);
void clear_overlap_region_alloc(overlap_region_alloc* list);
void destory_overlap_region_alloc(overlap_region_alloc* list);
void append_window_list(overlap_region* region, uint64_t x_start, uint64_t x_end, int y_start, int y_end, int error,
int extra_begin, int extra_end, int error_threshold);
void overlap_region_sort_y_id(overlap_region *a, long long n);
void calculate_overlap_region_by_chaining(Candidates_list* candidates, overlap_region_alloc* overlap_list,
uint64_t readID, uint64_t readLength, All_reads* R_INF, double band_width_threshold, int add_beg_end);
void init_fake_cigar(Fake_Cigar* x);
void destory_fake_cigar(Fake_Cigar* x);
void clear_fake_cigar(Fake_Cigar* x);
@@ -505,14 +161,14 @@ void add_fake_cigar(Fake_Cigar* x, uint32_t gap_site, int32_t gap_shift);
void resize_fake_cigar(Fake_Cigar* x, uint64_t size);
int get_fake_gap_pos(Fake_Cigar* x, int index);
int get_fake_gap_shift(Fake_Cigar* x, int index);
inline long long y_start_offset(long long x_start, Fake_Cigar* o)
static inline long long y_start_offset(long long x_start, Fake_Cigar* o)
{
if(x_start == get_fake_gap_pos(o, o->length - 1))
{
return get_fake_gap_shift(o, o->length - 1);
}
long long i;
for (i = 0; i < (long long)o->length; i++)
{
@@ -524,7 +180,7 @@ inline long long y_start_offset(long long x_start, Fake_Cigar* o)
if(i == 0 || i == (long long)o->length)
{
fprintf(stderr, "ERROR\n");
fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__);
exit(0);
}
@@ -532,24 +188,11 @@ inline long long y_start_offset(long long x_start, Fake_Cigar* o)
return get_fake_gap_shift(o, i - 1);
}
inline void print_fake_gap(Fake_Cigar* o)
{
long long i;
for (i = 0; i < (long long)o->length; i++)
{
fprintf(stderr, "**i: %lld, gap_pos_in_x: %d, gap_shift: %d\n",
i, get_fake_gap_pos(o, i),
get_fake_gap_shift(o, i));
}
}
void resize_Chain_Data(Chain_Data* x, long long size);
void init_window_list_alloc(window_list_alloc* x);
void clear_window_list_alloc(window_list_alloc* x);
void destory_window_list_alloc(window_list_alloc* x);
void resize_window_list_alloc(window_list_alloc* x, long long size);
void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result,
double band_width_threshold, int max_skip, int x_readLen, int y_readLen);
void chain_DP(k_mer_hit* a, long long a_n, Chain_Data* dp, overlap_region* result, double band_width_threshold, int max_skip, int x_readLen, int y_readLen);
#endif
#endif
Makefile
+32 -24
View File
@@ -1,11 +1,12 @@
CXX= g++
CXXFLAGS= -g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall #-fsanitize=address -fno-omit-frame-pointer#-Winline
CXXFLAGS= -g -O3 -msse4.2 -mpopcnt -fomit-frame-pointer -Wall
CPPFLAGS=
INCLUDES=
OBJS= Output.o CommandLines.o Process_Read.o Assembly.o kmer.o Hash_Table.o \
POA.o Correct.o Levenshtein_distance.o Overlaps.o Trio.o kthread.o Purge_Dups.o #ksw2_extz2_sse.o
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 Purge_Dups.o \
htab.o hist.o sketch.o anchor.o sys.o
EXE= hifiasm
LIBS= -lz -lpthread -lm #-fsanitize=address -fno-omit-frame-pointer
LIBS= -lz -lpthread -lm
ifneq ($(asan),)
CXXFLAGS+=-fsanitize=address
@@ -31,27 +32,34 @@ depend:
# DO NOT DELETE
Assembly.o: Assembly.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h
Assembly.o: CommandLines.h kmer.h Hash_Table.h khash.h POA.h Correct.h
Assembly.o: Levenshtein_distance.h Output.h Trio.h
Assembly.o: Assembly.h CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.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 khash.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 #ksw2.h
Hash_Table.o: Hash_Table.h khash.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
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 kseq.h
Overlaps.o: CommandLines.h Purge_Dups.h
POA.o: POA.h Hash_Table.h khash.h kmer.h Process_Read.h kseq.h Overlaps.h
POA.o: kvec.h kdq.h CommandLines.h Correct.h Levenshtein_distance.h
Process_Read.o: Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h
kmer.o: kmer.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h CommandLines.h
main.o: CommandLines.h Process_Read.h kseq.h Overlaps.h kvec.h kdq.h
main.o: Assembly.h Levenshtein_distance.h
Trio.o: Trio.h khashl.h kthread.h Process_Read.h CommandLines.h
Overlaps.o: Overlaps.h kvec.h kdq.h ksort.h Process_Read.h CommandLines.h
Overlaps.o: Hash_Table.h htab.h Correct.h Levenshtein_distance.h POA.h
Overlaps.o: Purge_Dups.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
Purge_Dups.o: ksort.h Purge_Dups.h kvec.h kdq.h Overlaps.h Hash_Table.h
Purge_Dups.o: htab.h Process_Read.h CommandLines.h Correct.h
Purge_Dups.o: Levenshtein_distance.h POA.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
kthread.o: kthread.h
Purge_Dups.o: Purge_Dups.h Overlaps.h Hash_Table.h Correct.h ksort.h
#ksw2_extz2_sse.o: ksw2.h
main.o: CommandLines.h Process_Read.h Overlaps.h kvec.h kdq.h Assembly.h
main.o: Levenshtein_distance.h htab.h
sketch.o: kvec.h htab.h Process_Read.h Overlaps.h kdq.h CommandLines.h
sys.o: htab.h Process_Read.h Overlaps.h kvec.h kdq.h CommandLines.h
Overlaps.cpp
+56 -71
View File
@@ -52,7 +52,6 @@ void ma_hit_sort_qns(ma_hit_t *a, long long n)
radix_sort_hit_qns(a, a + n);
}
void sort_kvec_t_u64_warp(kvec_t_u64_warp* u_vecs, uint32_t is_descend)
{
radix_sort_arch64(u_vecs->a.a, u_vecs->a.a + u_vecs->a.n);
@@ -68,7 +67,6 @@ void sort_kvec_t_u64_warp(kvec_t_u64_warp* u_vecs, uint32_t is_descend)
}
}
asg_t *asg_init(void)
{
return (asg_t*)calloc(1, sizeof(asg_t));
@@ -4356,7 +4354,7 @@ uint32_t startNode, uint32_t endNode)
}
else
{
fprintf(stderr, "ERROR\n");
fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__);
}
if(asg_arc_a(g, N_list[2])[0].v == (N_list[0]^1))
@@ -4369,7 +4367,7 @@ uint32_t startNode, uint32_t endNode)
}
else
{
fprintf(stderr, "ERROR\n");
fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__);
}
if(N_list[3] != N_list[4])
@@ -4920,12 +4918,10 @@ int asg_arc_del_short_diploid_unclean(asg_t *g, float drop_ratio, ma_hit_t_alloc
uint32_t detect_single_path_with_dels(asg_t *g, uint32_t begNode, uint32_t* endNode, long long* Len, buf_t* b)
{
uint32_t v = begNode, w;
uint32_t kv, kw;
(*Len) = 0;
while (1)
{
(*Len)++;
@@ -5594,6 +5590,7 @@ ma_hit_t_alloc* reverse_sources, long long min_edge_length, R_to_U* ruIndex)
if(flag1 == LOOP || flag2 == LOOP)
{
free(b_0.b.a); free(b_1.b.a);
return -1;
}
@@ -5612,6 +5609,7 @@ ma_hit_t_alloc* reverse_sources, long long min_edge_length, R_to_U* ruIndex)
if(l1 <= min_edge_length || l2 <= min_edge_length)
{
free(b_0.b.a); free(b_1.b.a);
return -1;
}
@@ -9051,7 +9049,7 @@ ma_sub_t *coverage_cut, int max_hang, int min_ovlp)
break;
}
}
if(k == nv) fprintf(stderr, "ERROR\n");
if(k == nv) fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__);
av = asg_arc_a(read_g, v);
nv = asg_arc_n(read_g, v);
@@ -9065,7 +9063,7 @@ ma_sub_t *coverage_cut, int max_hang, int min_ovlp)
}
}
if(k == nv) fprintf(stderr, "ERROR\n");
if(k == nv) fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__);
if(pE->el == 1 && aE->el == 1) continue;
@@ -9146,7 +9144,7 @@ ma_sub_t *coverage_cut, int max_hang, int min_ovlp)
}
l = asg_arc_len(t_f);
}
if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n");
if(l == (uint32_t)-1) fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__);
/*******************************for debug************************************/
@@ -9195,7 +9193,7 @@ ma_sub_t *coverage_cut, int max_hang, int min_ovlp)
}
l = asg_arc_len(t_f);
}
if(l == (uint32_t)-1) fprintf(stderr, "ERROR\n");
if(l == (uint32_t)-1) fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__);
/*******************************for debug************************************/
@@ -11671,7 +11669,7 @@ void print_untig(ma_ug_t *g, uint32_t uId, const char* info, uint32_t is_print_r
for (k = 0; k < u->n; k++)
{
fprintf(stderr, "%s: rId>>1: %lu, dir: %lu, name: %.*s\n",
info, u->a[k]>>33, (u->a[k]>>32)&1,
info, (unsigned long)(u->a[k]>>33), (unsigned long)((u->a[k]>>32)&1),
(int)Get_NAME_LENGTH((R_INF), (u->a[k]>>33)), Get_NAME((R_INF), (u->a[k]>>33)));
}
}
@@ -13624,6 +13622,9 @@ int load_ma_hit_ts(ma_hit_t_alloc** x, char* read_file_name)
f_flag += fread(&((*x)[i].length), sizeof((*x)[i].length), 1, fp);
(*x)[i].size = (*x)[i].length;
(*x)[i].buffer = NULL;
if((*x)[i].length == 0) continue;
(*x)[i].buffer = (ma_hit_t*)malloc(sizeof(ma_hit_t)*(*x)[i].length);
for (k = 0; k < (*x)[i].length; k++)
@@ -13640,8 +13641,6 @@ int load_ma_hit_ts(ma_hit_t_alloc** x, char* read_file_name)
}
void write_ma(ma_hit_t* x, FILE* fp)
{
fwrite(&(x->qns), sizeof(x->qns), 1, fp);
@@ -13692,51 +13691,43 @@ void write_ma_hit_ts(ma_hit_t_alloc* x, long long n_read, char* read_file_name)
fprintf(stderr, "ma_hit_ts has been written.\n");
}
void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
All_reads *RNF, char* output_file_name)
void write_all_data_to_disk(ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources, All_reads *RNF, char* output_file_name)
{
char* gfa_name = (char*)malloc(strlen(output_file_name)+25);
sprintf(gfa_name, "%s.ovlp", output_file_name);
write_All_reads(RNF, gfa_name);
sprintf(gfa_name, "%s.ovlp.source", output_file_name);
write_ma_hit_ts(sources, RNF->total_reads, gfa_name);
char* gfa_name = (char*)malloc(strlen(output_file_name)+25);
sprintf(gfa_name, "%s.ec", output_file_name);
write_All_reads(RNF, gfa_name);
sprintf(gfa_name, "%s.ovlp.reverse", output_file_name);
write_ma_hit_ts(reverse_sources, RNF->total_reads, gfa_name);
sprintf(gfa_name, "%s.ovlp.source", output_file_name);
write_ma_hit_ts(sources, RNF->total_reads, gfa_name);
free(gfa_name);
sprintf(gfa_name, "%s.ovlp.reverse", output_file_name);
write_ma_hit_ts(reverse_sources, RNF->total_reads, gfa_name);
free(gfa_name);
}
int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources,
char* output_file_name)
int load_all_data_from_disk(ma_hit_t_alloc **sources, ma_hit_t_alloc **reverse_sources, char* output_file_name)
{
char* gfa_name = (char*)malloc(strlen(output_file_name)+25);
sprintf(gfa_name, "%s.ovlp", output_file_name);
if(!load_All_reads(&R_INF, gfa_name))
{
return 0;
}
sprintf(gfa_name, "%s.ovlp.source", output_file_name);
if(!load_ma_hit_ts(sources, gfa_name))
{
return 0;
}
sprintf(gfa_name, "%s.ovlp.reverse", output_file_name);
if(!load_ma_hit_ts(reverse_sources, gfa_name))
{
return 0;
}
return 1;
char* gfa_name = (char*)malloc(strlen(output_file_name)+25);
sprintf(gfa_name, "%s.ec", output_file_name);
if (!load_All_reads(&R_INF, gfa_name)) {
free(gfa_name);
return 0;
}
sprintf(gfa_name, "%s.ovlp.source", output_file_name);
if (!load_ma_hit_ts(sources, gfa_name)) {
free(gfa_name);
return 0;
}
sprintf(gfa_name, "%s.ovlp.reverse", output_file_name);
if (!load_ma_hit_ts(reverse_sources, gfa_name)) {
free(gfa_name);
return 0;
}
free(gfa_name);
return 1;
}
// count the number of outgoing arcs, excluding reduced arcs
static inline int count_out(const asg_t *g, uint32_t v)
{
@@ -21006,7 +20997,7 @@ void merge_ug_nodes(ma_ug_t *ug, asg_t* read_g, kvec_t_u64_warp* array)
if(aw[i].del) continue;
if(aw[i].v == (v^1)) break;
}
if(i == nw) fprintf(stderr, "ERROR\n");
if(i == nw) fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__);
kmp = kmp | (uint64_t)(aw[i].ol);
@@ -21037,7 +21028,7 @@ void merge_ug_nodes(ma_ug_t *ug, asg_t* read_g, kvec_t_u64_warp* array)
if(aw[i].del) continue;
if(aw[i].v == (v^1)) break;
}
if(i == nw) fprintf(stderr, "ERROR\n");
if(i == nw) fprintf(stderr, "ERROR at %s:%d\n", __FILE__, __LINE__);
kmp = kmp | (uint64_t)(aw[i].ol);
@@ -26198,7 +26189,7 @@ void fix_binned_reads(ma_hit_t_alloc* paf, uint64_t n_read, ma_sub_t* coverage_c
}
fprintf(stderr, "n_read: %lu, binned_reads: %lu, binned_error_reads: %lu\n",
n_read, binned_reads, binned_error_reads);
(unsigned long)n_read, (unsigned long)binned_reads, (unsigned long)binned_error_reads);
}
@@ -26301,24 +26292,22 @@ int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
long long n_read, uint64_t* readLen, long long mini_overlap_length,
long long max_hang_length, long long clean_round, long long gap_fuzz,
float min_ovlp_drop_ratio, float max_ovlp_drop_ratio, char* output_file_name,
long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t* sg,
ma_sub_t* coverage_cut, int debug_g)
long long bubble_dist, int read_graph, R_to_U* ruIndex, asg_t **sg_ptr,
ma_sub_t **coverage_cut_ptr, int debug_g)
{
ma_sub_t *coverage_cut = *coverage_cut_ptr;
asg_t *sg = *sg_ptr;
if(debug_g) goto debug_gfa;
///just for debug
renew_graph_init(sources, reverse_sources, sg, coverage_cut, ruIndex, n_read);
///it's hard to say which function is better
///normalize_ma_hit_t_single_side(sources, n_read);
normalize_ma_hit_t_single_side_advance(sources, n_read);
normalize_ma_hit_t_single_side_advance(reverse_sources, n_read);
// debug_info_of_specfic_read(">m64011_190830_220126/117834372/ccs", sources, reverse_sources,
// -1, "clean");
@@ -26477,7 +26466,7 @@ ma_sub_t* coverage_cut, int debug_g)
asg_arc_del_simple_circle_untig(sources, coverage_cut, sg, 100, 0);
if(VERBOSE_GFA >= 1)
if (asm_opt.verbose_gfa)
{
/*******************************for debug***************************************/
write_debug_graph(sg, sources, coverage_cut, output_file_name, n_read, reverse_sources, ruIndex);
@@ -26532,13 +26521,13 @@ ma_sub_t* coverage_cut, int debug_g)
bubble_dist, (asm_opt.max_short_tip*2), 0.15, 3, ruIndex, 0.05, 0.9, max_hang_length,
mini_overlap_length);
output_contig_graph_alternative(sg, coverage_cut, output_file_name, sources, max_hang_length, mini_overlap_length);
}
*coverage_cut_ptr = coverage_cut;
*sg_ptr = sg;
}
void build_string_graph_without_clean(
int min_dp, ma_hit_t_alloc* sources, ma_hit_t_alloc* reverse_sources,
long long n_read, uint64_t* readLen, long long mini_overlap_length,
@@ -26551,13 +26540,12 @@ long long bubble_dist, int read_graph, int write)
asg_t *sg = NULL;
ma_sub_t* coverage_cut = NULL;
// debug_info_of_specfic_read("m64011_190329_072846/80545633/ccs",
// sources, reverse_sources, -1, "clean");
// debug_info_of_specfic_read("m64011_190329_072846/80545633/ccs", sources, reverse_sources, -1, "clean");
///actually min_thres = asm_opt.max_short_tip + 1 there are asm_opt.max_short_tip reads
min_thres = asm_opt.max_short_tip + 1;
if(VERBOSE_GFA >= 1)
if (asm_opt.verbose_gfa)
{
if(load_debug_graph(&sg, &sources, &coverage_cut, output_file_name, &reverse_sources, &ruIndex))
{
@@ -26565,7 +26553,7 @@ long long bubble_dist, int read_graph, int write)
clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length,
max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
output_file_name, bubble_dist, read_graph, &ruIndex, sg, coverage_cut, 1);
output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 1);
asg_destroy(sg);
free(coverage_cut);
destory_R_to_U(&ruIndex);
@@ -26573,20 +26561,17 @@ long long bubble_dist, int read_graph, int write)
}
}
if (asm_opt.write_index_to_disk && write)
{
write_all_data_to_disk(sources, reverse_sources,
&R_INF, output_file_name);
}
try_rescue_overlaps(sources, reverse_sources, n_read, 4);
clean_graph(min_dp, sources, reverse_sources, n_read, readLen, mini_overlap_length,
max_hang_length, clean_round, gap_fuzz, min_ovlp_drop_ratio, max_ovlp_drop_ratio,
output_file_name, bubble_dist, read_graph, &ruIndex, sg, coverage_cut, 0);
output_file_name, bubble_dist, read_graph, &ruIndex, &sg, &coverage_cut, 0);
asg_destroy(sg);
free(coverage_cut);
Overlaps.h
+1 -4
View File
@@ -47,9 +47,6 @@
#define CUT_DIF_HAP 12
///query is the read itself
typedef struct {
uint64_t qns;
@@ -1057,4 +1054,4 @@ int unitig_arc_del_short_diploid_by_length(asg_t *g, float drop_ratio);
#define JUNK_COV 5
#define DISCARD_RATE 0.8
#endif
#endif
POA.cpp
+73 -131
View File
@@ -1,29 +1,27 @@
#include "POA.h"
#include <stdlib.h>
#include <string.h>
#include "POA.h"
#include "Correct.h"
#include "Process_Read.h"
#define INIT_EDGE_SIZE 50
#define INCREASE_EDGE_SIZE 5
#define INIT_NODE_SIZE 16000
/********
* Edge *
********/
void init_Edge_alloc(Edge_alloc* list)
{
if (list->list == NULL)
{
list->size = INIT_EDGE_SIZE;
list->length = 0;
list->delete_length = 0;
list->list = (Edge*)malloc(sizeof(Edge)*list->size);
}
else
{
list->length = 0;
list->delete_length = 0;
}
if (list->list == NULL) {
list->size = INIT_EDGE_SIZE;
list->length = 0;
list->delete_length = 0;
list->list = (Edge*)malloc(sizeof(Edge)*list->size);
} else {
list->length = 0;
list->delete_length = 0;
}
}
void clear_Edge_alloc(Edge_alloc* list)
@@ -34,28 +32,29 @@ void clear_Edge_alloc(Edge_alloc* list)
void destory_Edge_alloc(Edge_alloc* list)
{
free(list->list);
if (list && list->list)
free(list->list);
}
void append_Edge_alloc(Edge_alloc* list, uint64_t in_node, uint64_t out_node, uint64_t weight, uint64_t length)
{
if (list->length + 1 > list->size)
{
list->size = list->size + INCREASE_EDGE_SIZE;
list->list = (Edge*)realloc(list->list, sizeof(Edge)*list->size);
}
if (list->length + 1 > list->size) {
uint64_t old_size = list->size;
list->size = list->size + INCREASE_EDGE_SIZE;
list->list = (Edge*)realloc(list->list, sizeof(Edge)*list->size);
memset(&list->list[old_size], 0, (list->size - old_size) * sizeof(Edge));
}
list->list[list->length].in_node = in_node;
list->list[list->length].out_node = out_node;
list->list[list->length].weight = weight;
list->list[list->length].length = length;
list->list[list->length].num_insertions = 0;
list->list[list->length].self_edge_ID = list->length;
list->list[list->length].in_node = in_node;
list->list[list->length].out_node = out_node;
list->list[list->length].weight = weight;
list->list[list->length].length = length;
list->list[list->length].num_insertions = 0;
list->list[list->length].self_edge_ID = list->length;
list->length++;
list->length++;
}
int add_and_check_bi_direction_edge(Graph* graph, Node* in_node, Node* out_node, uint64_t weight, uint64_t flag)
{
Edge* e_forward;
@@ -81,7 +80,6 @@ int add_and_check_bi_direction_edge(Graph* graph, Node* in_node, Node* out_node,
}
}
void add_bi_direction_edge(Graph* graph, Node* in_node, Node* out_node, uint64_t weight, uint64_t flag)
{
@@ -95,8 +93,6 @@ void add_bi_direction_edge(Graph* graph, Node* in_node, Node* out_node, uint64_t
= Output_Edges((*in_node)).length - 1;
}
int remove_and_check_bi_direction_edge_from_nodes(Graph* graph, Node* in_node, Node* out_node)
{
Edge* e_forward;
@@ -135,8 +131,6 @@ int remove_and_check_bi_direction_edge_from_nodes(Graph* graph, Node* in_node, N
}
}
int remove_and_check_bi_direction_edge_from_edge(Graph* graph, Edge* e)
{
Edge* e_forward;
@@ -171,109 +165,71 @@ int remove_and_check_bi_direction_edge_from_edge(Graph* graph, Edge* e)
{
return 0;
}
}
/********
* Node *
********/
void init_Node_alloc(Node_alloc* list)
{
list->size = INIT_NODE_SIZE;
list->length = 0;
list->delete_length = 0;
list->list = (Node*)malloc(sizeof(Node)*list->size);
list->sort.size = 0;
list->sort.list = NULL;
list->sort.visit = NULL;
list->sort.iterative_buffer = NULL;
list->sort.iterative_buffer_visit = NULL;
uint64_t i;
for (i = 0; i < list->size; i++)
{
list->list[i].insertion_edges.list=NULL;
list->list[i].mismatch_edges.list=NULL;
list->list[i].deletion_edges.list=NULL;
}
memset(list, 0, sizeof(Node_alloc));
list->size = INIT_NODE_SIZE;
list->list = (Node*)calloc(list->size, sizeof(Node));
}
void destory_Node_alloc(Node_alloc* list)
{
uint64_t i =0;
for (i = 0; i < list->length; i++)
{
destory_Edge_alloc(&list->list[i].deletion_edges);
destory_Edge_alloc(&list->list[i].insertion_edges);
destory_Edge_alloc(&list->list[i].mismatch_edges);
}
free(list->list);
free(list->sort.list);
free(list->sort.visit);
free(list->sort.iterative_buffer);
free(list->sort.iterative_buffer_visit);
///free(list->topo_order);
uint64_t i;
for (i = 0; i < list->size; i++) {
destory_Edge_alloc(&list->list[i].deletion_edges);
destory_Edge_alloc(&list->list[i].insertion_edges);
destory_Edge_alloc(&list->list[i].mismatch_edges);
}
free(list->list);
free(list->sort.list);
free(list->sort.visit);
free(list->sort.iterative_buffer);
free(list->sort.iterative_buffer_visit);
}
void clear_Node_alloc(Node_alloc* list)
{
uint64_t i =0;
for (i = 0; i < list->length; i++)
{
clear_Edge_alloc(&list->list[i].insertion_edges);
clear_Edge_alloc(&list->list[i].mismatch_edges);
clear_Edge_alloc(&list->list[i].deletion_edges);
}
list->length = 0;
list->delete_length = 0;
uint64_t i =0;
for (i = 0; i < list->length; i++) { // TODO: is this list->size or list->length? The original version is list->length.
clear_Edge_alloc(&list->list[i].insertion_edges);
clear_Edge_alloc(&list->list[i].mismatch_edges);
clear_Edge_alloc(&list->list[i].deletion_edges);
}
list->length = 0;
list->delete_length = 0;
}
uint64_t append_Node_alloc(Node_alloc* list, char base)
{
if (list->length + 1 > list->size)
{
uint64_t i = list->size;
if (list->length + 1 > list->size) {
uint64_t old_size = list->size;
list->size = list->size * 2;
list->list = (Node*)realloc(list->list, sizeof(Node) * list->size);
memset(&list->list[old_size], 0, (list->size - old_size) * sizeof(Node));
}
list->size = list->size * 2;
list->list = (Node*)realloc(list->list, sizeof(Node)*list->size);
///list->topo_order = (uint64_t*)realloc(list->topo_order, sizeof(uint64_t)*list->size);
for (; i < list->size; i++)
{
list->list[i].deletion_edges.list=NULL;
list->list[i].insertion_edges.list=NULL;
list->list[i].mismatch_edges.list=NULL;
}
}
list->list[list->length].ID = list->length;
list->list[list->length].base = base;
list->list[list->length].weight = 1;
list->list[list->length].num_insertions = 0;
init_Edge_alloc(&list->list[list->length].deletion_edges);
init_Edge_alloc(&list->list[list->length].insertion_edges);
init_Edge_alloc(&list->list[list->length].mismatch_edges);
list->length++;
list->list[list->length].ID = list->length;
list->list[list->length].base = base;
list->list[list->length].weight = 1;
list->list[list->length].num_insertions = 0;
init_Edge_alloc(&list->list[list->length].deletion_edges);
init_Edge_alloc(&list->list[list->length].insertion_edges);
init_Edge_alloc(&list->list[list->length].mismatch_edges);
return list->length - 1;
list->length++;
return list->length - 1;
}
/*********
* Graph *
*********/
void init_Graph(Graph* g)
{
@@ -310,12 +266,6 @@ void clear_Graph(Graph* g)
clear_Queue(&(g->node_q));
}
void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length, long long* startID, long long* endID)
{
long long firstID, lastID, nodeID, i;
@@ -356,8 +306,6 @@ void addUnmatchedSeqToGraph(Graph* g, char* g_read_seq, long long g_read_length,
}
void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_string, long long x_length,
char* y_string, long long y_length, CIGAR* cigar, long long backbone_start, long long backbone_end)
{
@@ -422,9 +370,3 @@ void addmatchedSeqToGraph(Graph* backbone, long long currentNodeID, char* x_stri
cigar_i++;
}
}
Process_Read.cpp
+26 -522
View File
@@ -1,25 +1,8 @@
#include "Process_Read.h"
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <fcntl.h>
#include <pthread.h>
gz_files fps;
R_buffer RDB;
static uint64_t total_reads;
pthread_mutex_t i_readinputMutex;
pthread_mutex_t i_queueMutex;
pthread_mutex_t i_terminateMutex;
pthread_cond_t i_flushCond;
pthread_cond_t i_readinputflushCond;
pthread_cond_t i_stallCond;
pthread_cond_t i_readinputstallCond;
pthread_mutex_t i_doneMutex;
#include "Process_Read.h"
uint8_t seq_nt6_table[256] = {
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
@@ -45,34 +28,27 @@ char bit_t_seq_table_rc[256][4] = {{0}};
char s_H[5] = {'A', 'C', 'G', 'T', 'N'};
char rc_Table[5] = {'T', 'G', 'C', 'A', 'N'};
void init_All_reads(All_reads* r)
{
memset(r, 0, sizeof(All_reads));
r->index_size = READ_INIT_NUMBER;
r->read_length = (uint64_t*)malloc(sizeof(uint64_t)*r->index_size);
r->read_sperate = NULL;
r->N_site = NULL;
r->total_reads_bases = 0;
r->name_index_size = READ_INIT_NUMBER;
r->name_index = (uint64_t*)malloc(sizeof(uint64_t)*r->name_index_size);
r->name_index[0] = 0;
r->name = NULL;
r->total_name_length = 0;
r->total_reads = 0;
r->trio_flag = NULL;
}
void destory_All_reads(All_reads* r)
{
uint64_t i = 0;
for (i = 0; i < r->total_reads; i++)
{
if (r->N_site[i] != NULL)
{
free(r->N_site[i]);
}
free(r->read_sperate[i]);
for (i = 0; i < r->total_reads; i++) {
if (r->N_site[i]) free(r->N_site[i]);
if (r->read_sperate[i]) free(r->read_sperate[i]);
if (r->paf&&r->paf[i].buffer) free(r->paf[i].buffer);
if (r->reverse_paf&&r->reverse_paf[i].buffer) free(r->reverse_paf[i].buffer);
}
free(r->paf);
free(r->reverse_paf);
free(r->N_site);
free(r->read_sperate);
free(r->name);
@@ -81,7 +57,6 @@ void destory_All_reads(All_reads* r)
free(r->trio_flag);
}
void write_All_reads(All_reads* r, char* read_file_name)
{
fprintf(stderr, "Writing reads to disk... \n");
@@ -112,9 +87,6 @@ void write_All_reads(All_reads* r, char* read_file_name)
{
fwrite(&zero, sizeof(zero), 1, fp);
}
}
fwrite(r->read_length, sizeof(uint64_t), r->total_reads, fp);
@@ -132,16 +104,13 @@ void write_All_reads(All_reads* r, char* read_file_name)
fprintf(stderr, "Reads has been written.\n");
}
int load_All_reads(All_reads* r, char* read_file_name)
{
fprintf(stderr, "Loading reads from disk... \n");
char* index_name = (char*)malloc(strlen(read_file_name)+15);
sprintf(index_name, "%s.bin", read_file_name);
FILE* fp = fopen(index_name, "r");
if (!fp)
{
if (!fp) {
free(index_name);
return 0;
}
int local_adapterLen;
@@ -164,12 +133,10 @@ int load_All_reads(All_reads* r, char* read_file_name)
r->N_site = (uint64_t**)malloc(sizeof(uint64_t*)*r->total_reads);
for (i = 0; i < r->total_reads; i++)
{
f_flag += fread(&zero, sizeof(zero), 1, fp);
if (zero)
{
r->N_site[i] = (uint64_t*)malloc(sizeof(uint64_t)*(zero + 1));
r->N_site[i][0] = zero;
if (r->N_site[i][0])
@@ -181,7 +148,6 @@ int load_All_reads(All_reads* r, char* read_file_name)
{
r->N_site[i] = NULL;
}
}
r->read_length = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads);
@@ -211,8 +177,6 @@ int load_All_reads(All_reads* r, char* read_file_name)
r->cigars = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
r->second_round_cigar = (Compressed_Cigar_record*)malloc(sizeof(Compressed_Cigar_record)*r->total_reads);
r->paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
r->reverse_paf = (ma_hit_t_alloc*)malloc(sizeof(ma_hit_t_alloc)*r->total_reads);
for (i = 0; i < r->total_reads; i++)
{
r->second_round_cigar[i].size = r->cigars[i].size = 0;
@@ -222,8 +186,6 @@ int load_All_reads(All_reads* r, char* read_file_name)
r->second_round_cigar[i].lost_base_size = r->cigars[i].lost_base_size = 0;
r->second_round_cigar[i].lost_base_length = r->cigars[i].lost_base_length = 0;
r->second_round_cigar[i].lost_base = r->cigars[i].lost_base = NULL;
init_ma_hit_t_alloc(&(r->paf[i]));
init_ma_hit_t_alloc(&(r->reverse_paf[i]));
}
free(index_name);
@@ -233,32 +195,28 @@ int load_All_reads(All_reads* r, char* read_file_name)
return 1;
}
inline void insert_read(All_reads* r, kstring_t* read, kstring_t* name)
void ha_insert_read_len(All_reads *r, int read_len, int name_len)
{
r->total_reads++;
r->total_reads_bases = r->total_reads_bases + read->l;
r->total_name_length = r->total_name_length + name->l;
r->total_reads_bases += (uint64_t)read_len;
r->total_name_length += (uint64_t)name_len;
///must +1
if (r->index_size < r->total_reads + 2)
{
// must +1
if (r->index_size < r->total_reads + 2) {
r->index_size = r->index_size * 2 + 2;
r->read_length = (uint64_t*)realloc(r->read_length,sizeof(uint64_t)*(r->index_size));
r->read_length = (uint64_t*)realloc(r->read_length, sizeof(uint64_t) * r->index_size);
r->name_index_size = r->name_index_size * 2 + 2;
r->name_index = (uint64_t*)realloc(r->name_index,sizeof(uint64_t)*(r->name_index_size));
r->name_index = (uint64_t*)realloc(r->name_index, sizeof(uint64_t) * r->name_index_size);
}
r->read_length[r->total_reads - 1] = read->l;
r->name_index[r->total_reads] = r->name_index[r->total_reads-1] + name->l;
r->read_length[r->total_reads - 1] = read_len;
r->name_index[r->total_reads] = r->name_index[r->total_reads - 1] + name_len;
}
void malloc_All_reads(All_reads* r)
{
r->read_size = (uint64_t*)malloc(sizeof(uint64_t)*r->total_reads);
memcpy (r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads);
memcpy(r->read_size, r->read_length, sizeof(uint64_t)*r->total_reads);
r->read_sperate = (uint8_t**)malloc(sizeof(uint8_t*)*r->total_reads);
long long i = 0;
@@ -313,7 +271,6 @@ void init_aux_table()
bit_t_seq_table_rc[i][2] = RC_CHAR(bit_t_seq_table[i][1]);
bit_t_seq_table_rc[i][3] = RC_CHAR(bit_t_seq_table[i][0]);
}
}
}
@@ -338,18 +295,12 @@ void init_UC_Read(UC_Read* r)
bit_t_seq_table_rc[i][2] = RC_CHAR(bit_t_seq_table[i][1]);
bit_t_seq_table_rc[i][3] = RC_CHAR(bit_t_seq_table[i][0]);
}
}
}
void recover_UC_Read_sub_region_begin_end
(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID, int extra_begin, int extra_end)
void recover_UC_Read_sub_region_begin_end(char* r, long long start_pos, long long length, uint8_t strand,
All_reads* R_INF, long long ID, int extra_begin, int extra_end)
{
long long readLen = Get_READ_LENGTH((*R_INF), ID);
uint8_t* src = Get_READ((*R_INF), ID);
@@ -357,17 +308,11 @@ void recover_UC_Read_sub_region_begin_end
long long copyLen;
long long end_pos = start_pos + length - 1;
if (strand == 0)
{
i = start_pos;
copyLen = 0;
long long initLen = start_pos % 4;
if (initLen != 0)
@@ -376,8 +321,6 @@ void recover_UC_Read_sub_region_begin_end
copyLen = copyLen + 4 - initLen;
i = i + copyLen;
}
while (copyLen < length)
{
memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4);
@@ -385,7 +328,6 @@ void recover_UC_Read_sub_region_begin_end
i = i + 4;
}
if (R_INF->N_site[ID])
{
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
@@ -400,17 +342,12 @@ void recover_UC_Read_sub_region_begin_end
}
}
}
}
else
{
start_pos = readLen - start_pos - 1;
end_pos = readLen - end_pos - 1;
///start_pos > end_pos
i = start_pos;
copyLen = 0;
@@ -436,7 +373,6 @@ void recover_UC_Read_sub_region_begin_end
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
{
if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos)
{
r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N';
@@ -447,21 +383,11 @@ void recover_UC_Read_sub_region_begin_end
}
}
}
}
}
void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID)
{
long long readLen = Get_READ_LENGTH((*R_INF), ID);
uint8_t* src = Get_READ((*R_INF), ID);
@@ -471,7 +397,6 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length,
if (strand == 0)
{
i = start_pos;
copyLen = 0;
@@ -483,8 +408,7 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length,
copyLen = copyLen + 4 - initLen;
i = i + copyLen;
}
while (copyLen < length)
{
memcpy(r+copyLen, bit_t_seq_table[src[i>>2]], 4);
@@ -492,7 +416,6 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length,
i = i + 4;
}
if (R_INF->N_site[ID])
{
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
@@ -507,17 +430,12 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length,
}
}
}
}
else
{
start_pos = readLen - start_pos - 1;
end_pos = readLen - end_pos - 1;
///start_pos > end_pos
i = start_pos;
copyLen = 0;
@@ -543,7 +461,6 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length,
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
{
if ((long long)R_INF->N_site[ID][i] >= end_pos && (long long)R_INF->N_site[ID][i] <= start_pos)
{
r[readLen - R_INF->N_site[ID][i] - 1 - offset] = 'N';
@@ -554,15 +471,11 @@ void recover_UC_Read_sub_region(char* r, long long start_pos, long long length,
}
}
}
}
}
void recover_UC_Read(UC_Read* r, All_reads* R_INF, uint64_t ID)
void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID)
{
r->length = Get_READ_LENGTH((*R_INF), ID);
uint8_t* src = Get_READ((*R_INF), ID);
@@ -623,7 +536,6 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID)
index = index + 4;
}
if (R_INF->N_site[ID])
{
for (i = 1; i <= (long long)R_INF->N_site[ID][0]; i++)
@@ -631,11 +543,8 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID)
r->seq[r->length - R_INF->N_site[ID][i] - 1] = 'N';
}
}
}
#define COMPRESS_BASE {c = seq_nt6_table[(uint8_t)src[i]];\
if (c >= 4)\
{\
@@ -645,9 +554,8 @@ void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID)
}\
i++;}\
void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ)
void ha_compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ)
{
///N_site_lis saves the pos of all Ns in this read
///N_site_lis[0] is the number of Ns
if (N_site_occ)
@@ -666,10 +574,8 @@ void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_l
uint8_t tmp = 0;
uint8_t c = 0;
while (i + 4 <= src_l)
{
tmp = 0;
COMPRESS_BASE;
@@ -705,351 +611,8 @@ void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_l
dest[dest_i] = tmp;
dest_i++;
}
}
void open_file(gz_files* nfps, char* name)
{
nfps->fp = gzopen(name, "r");
if(nfps->fp == 0)
{
fprintf(stderr, "[ERROR] Cannot find the input file: %s\n", name);
exit(0);
}
nfps->seq = kseq_init(nfps->fp);
}
void close_file(gz_files* nfps)
{
kseq_destroy(nfps->seq);
gzclose(nfps->fp);
}
void init_gz_files(hifiasm_opt_t* asm_opt)
{
fps.idx = 0;
fps.num_reads = asm_opt->num_reads;
fps.reads = asm_opt->read_file_names;
fps.seq = NULL;
fps.fp = NULL;
if(fps.num_reads > 0)
{
open_file(&fps, fps.reads[fps.idx]);
fps.idx++;
}
}
void destory_gz_files()
{
close_file(&fps);
}
int read_item()
{
int l = kseq_read(fps.seq);
if(l >= 0 || (l < 0 && fps.idx >= fps.num_reads))
{
return l;
}
///l < 0 && fps.idx < fps.num_reads
close_file(&fps);
open_file(&fps, fps.reads[fps.idx]);
fps.idx++;
return read_item();
}
inline void exchage_kstring_t(kstring_t* a, kstring_t* b)
{
kstring_t tmp;
tmp = *a;
*a = *b;
*b = tmp;
}
int get_read(kseq_t *s, int adapterLen)
{
int l;
///if ((l = kseq_read(seq)) >= 0)
if ((l = read_item()) >= 0)
{
exchage_kstring_t(&(fps.seq->comment), &s->comment);
exchage_kstring_t(&(fps.seq->name), &s->name);
exchage_kstring_t(&(fps.seq->qual), &s->qual);
exchage_kstring_t(&(fps.seq->seq), &s->seq);
if(adapterLen > 0)
{
if((int)s->seq.l <= adapterLen*2)
{
s->seq.l = 0;
}
else
{
long long i;
for (i = 0; i < ((int)s->seq.l - adapterLen*2); i++)
{
s->seq.s[i] = s->seq.s[i + adapterLen];
}
s->seq.l -= adapterLen*2;
}
}
return 1;
}
else
{
return 0;
}
}
void init_R_buffer_block(R_buffer_block* curr_sub_block)
{
curr_sub_block->read = (kseq_t*)calloc(RDB.block_inner_size, sizeof(kseq_t));
curr_sub_block->num = 0;
}
void clear_R_buffer()
{
RDB.all_read_end = 0;
RDB.num = 0;
}
void init_R_buffer(int thread_num)
{
RDB.all_read_end = 0;
RDB.num = 0;
RDB.block_inner_size = READ_BLOCK_SIZE;
RDB.size = thread_num*READ_BLOCK_NUM_PRE_THR;
RDB.sub_block = (R_buffer_block*)malloc(sizeof(R_buffer_block)*RDB.size);
int i = 0;
for (i = 0; i < RDB.size; i++)
{
init_R_buffer_block(&RDB.sub_block[i]);
}
}
void destory_R_buffer_block(R_buffer_block* curr_sub_block)
{
kseq_destroy(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, int is_insert, int adapterLen)
{
int inner_i = 0;
int file_flag = 1;
while (inner_i<batch_read_size)
{
file_flag = get_read(&read_batch->read[inner_i], adapterLen);
if (file_flag == 1)
{
read_batch->read[inner_i].ID = total_reads;
total_reads++;
if (is_insert)
{
insert_read(&R_INF, &read_batch->read[inner_i].seq,
&read_batch->read[inner_i].name);
}
inner_i++;
}
else if (file_flag == 0)
{
break;
}
}
if (inner_i || file_flag)
{
file_flag = 1;
}
*return_file_flag = file_flag;
read_batch->num = inner_i;
}
inline void push_R_block(R_buffer_block* tmp_sub_block)
{
///only exchange pointers
kseq_t *k1;
k1 = RDB.sub_block[RDB.num].read;
RDB.sub_block[RDB.num].read = tmp_sub_block->read;
tmp_sub_block->read = k1;
RDB.sub_block[RDB.num].num = tmp_sub_block->num;
tmp_sub_block->num = 0;
RDB.num++;
}
inline void pop_R_block(R_buffer_block* curr_sub_block)
{
RDB.num--;
///only exchange pointers
kseq_t *k1;
k1 = RDB.sub_block[RDB.num].read;
RDB.sub_block[RDB.num].read = curr_sub_block->read;
curr_sub_block->read = k1;
curr_sub_block->num = RDB.sub_block[RDB.num].num;
RDB.sub_block[RDB.num].num = 0;
}
void* input_reads_muti_threads(void* arg)
{
int is_insert = *((int*)arg);
total_reads = 0;
int file_flag = 1;
R_buffer_block tmp_buf;
init_R_buffer_block(&tmp_buf);
while (1)
{
load_read_block(&tmp_buf, RDB.block_inner_size, &file_flag, is_insert, asm_opt.adapterLen);
if (file_flag == 0)
{
break;
}
pthread_mutex_lock(&i_readinputMutex);
while (IS_FULL(RDB))
{
pthread_cond_signal(&i_readinputstallCond);
pthread_cond_wait(&i_readinputflushCond, &i_readinputMutex);
}
push_R_block(&tmp_buf);
pthread_cond_signal(&i_readinputstallCond);
pthread_mutex_unlock(&i_readinputMutex);
}
pthread_mutex_lock(&i_readinputMutex);
RDB.all_read_end = 1;
pthread_cond_signal(&i_readinputstallCond); //important
pthread_mutex_unlock(&i_readinputMutex);
destory_R_buffer_block(&tmp_buf);
fprintf(stderr, "Reads #: %lu\n", (unsigned long)total_reads);
fprintf(stderr, "Bases #: %lu\n", (unsigned long)R_INF.total_reads_bases);
return NULL;
}
int get_reads_mul_thread(R_buffer_block* curr_sub_block)
{
pthread_mutex_lock(&i_readinputMutex);
while (IS_EMPTY(RDB) && RDB.all_read_end == 0)
{
pthread_cond_signal(&i_readinputflushCond);
pthread_cond_wait(&i_readinputstallCond, &i_readinputMutex);
}
if (!IS_EMPTY(RDB))
{
pop_R_block(curr_sub_block);
pthread_cond_signal(&i_readinputflushCond);
pthread_mutex_unlock(&i_readinputMutex);
return 1;
}
else
{
curr_sub_block->num = 0;
pthread_cond_signal(&i_readinputstallCond); //important
pthread_mutex_unlock(&i_readinputMutex);
return 0;
}
}
void reverse_complement(char* pattern, uint64_t length)
{
uint64_t i = 0;
@@ -1069,63 +632,4 @@ void reverse_complement(char* pattern, uint64_t length)
{
pattern[end] = RC_CHAR(pattern[end]);
}
}
typedef struct {
char* tmp;
long long tmpSize;
char* dest;
long long destSize;
FILE* fp;
} LineReader;
int get_single_line(LineReader* line)
{
long long currentLen = 0, getLen = 0;
line->tmp[0] = '\0';
while (fgets(line->tmp, line->tmpSize, line->fp) != NULL)
{
getLen = strlen(line->tmp);
if(getLen + currentLen >= line->destSize)
{
line->destSize = getLen + currentLen + 1;
line->dest = (char*)realloc(line->dest, line->destSize);
}
memcpy(line->dest + currentLen, line->tmp, getLen+1);
currentLen = currentLen + getLen;
if(currentLen > 0 && line->dest[currentLen - 1] == '\n')
{
return 1;
}
}
if(currentLen > 0)
{
return 1;
}
else
{
return 0;
}
}
void get_trio_info(char* input, uint8_t* pm)
{
uint32_t i;
(*pm) = AMBIGU;
for (i = 0; input[i] != '\0'; i++)
{
if(input[i] == '\t')
{
break;
}
}
input[i] = '\0';
i++;
if(input[i] == 'p') (*pm) = FATHER;
if(input[i] == 'm') (*pm) = MOTHER;
}
Process_Read.h
+10 -65
View File
@@ -5,7 +5,6 @@
#include <string.h>
#include <stdlib.h>
#include <zlib.h>
#include "kseq.h"
#include "Overlaps.h"
#include "CommandLines.h"
///#include "Hash_Table.h"
@@ -18,16 +17,11 @@
#define IS_FULL(buffer) ((buffer.num >= buffer.size)?1:0)
#define IS_EMPTY(buffer) ((buffer.num == 0)?1:0)
///#define Get_READ_LENGTH(R_INF, ID) (R_INF.index[ID+1] - R_INF.index[ID])
#define Get_READ_LENGTH(R_INF, ID) R_INF.read_length[(ID)]
#define Get_NAME_LENGTH(R_INF, ID) (R_INF.name_index[(ID)+1] - R_INF.name_index[(ID)])
#define Get_READ_LENGTH(R_INF, ID) (R_INF).read_length[(ID)]
#define Get_NAME_LENGTH(R_INF, ID) ((R_INF).name_index[(ID)+1] - (R_INF).name_index[(ID)])
///#define Get_READ(R_INF, ID) R_INF.read + (R_INF.index[ID]>>2) + ID
#define Get_READ(R_INF, ID) R_INF.read_sperate[(ID)]
#define Get_NAME(R_INF, ID) R_INF.name + R_INF.name_index[(ID)]
KSEQ_INIT(gzFile, gzread)
#define Get_READ(R_INF, ID) (R_INF).read_sperate[(ID)]
#define Get_NAME(R_INF, ID) ((R_INF).name + (R_INF).name_index[(ID)])
extern uint8_t seq_nt6_table[256];
@@ -37,12 +31,9 @@ extern char s_H[5];
extern char rc_Table[5];
#define RC_CHAR(x) rc_Table[seq_nt6_table[(uint8_t)x]]
void init_aux_table();
int get_read(kseq_t *s, int adapterLen);
typedef struct
{
@@ -76,7 +67,6 @@ inline void init_PAF_alloc(PAF_alloc* list)
list->list = (PAF*)malloc(sizeof(PAF)*list->size);
}
inline void append_PAF_alloc(PAF_alloc* list, PAF* e)
{
if(list->length+1 > list->size)
@@ -89,9 +79,6 @@ inline void append_PAF_alloc(PAF_alloc* list, PAF* e)
list->length++;
}
typedef struct
{
/**[0-1] bits are type:**/
@@ -104,7 +91,7 @@ typedef struct
uint32_t lost_base_length;
uint32_t lost_base_size;
uint32_t new_length;
}Compressed_Cigar_record;
} Compressed_Cigar_record;
#define AMBIGU 0
#define FATHER 1
@@ -112,13 +99,13 @@ typedef struct
#define MIX_TRIO 3
#define NON_TRIO 4
#define DROP 5
typedef struct
{
uint64_t** N_site;
///uint8_t* read;
char* name;
uint8_t** read_sperate;
uint64_t* read_length;
uint64_t* read_size;
@@ -129,7 +116,6 @@ typedef struct
///uint64_t* index;
uint64_t index_size;
///name start pos in char* name
uint64_t* name_index;
uint64_t name_index_size;
@@ -142,32 +128,10 @@ typedef struct
ma_hit_t_alloc* paf;
ma_hit_t_alloc* reverse_paf;
ma_sub_t* coverage_cut;
} All_reads;
extern All_reads R_INF;
void malloc_All_reads(All_reads* r);
typedef struct
{
kseq_t* read;
long long num;
} R_buffer_block;
typedef struct
{
R_buffer_block* sub_block;
long long block_inner_size;
long long size;
long long num;
int all_read_end;
} R_buffer;
typedef struct
{
char* seq;
@@ -176,23 +140,12 @@ typedef struct
long long RID;
} UC_Read;
typedef struct
{
gzFile fp;
kseq_t *seq;
char** reads;
int num_reads;
int idx;
} gz_files;
void init_R_buffer(int thread_num);
void init_All_reads(All_reads* r);
void* input_reads_muti_threads(void*);
void init_R_buffer_block(R_buffer_block* curr_sub_block);
int get_reads_mul_thread(R_buffer_block* curr_sub_block);
void compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ);
void malloc_All_reads(All_reads* r);
void ha_insert_read_len(All_reads *r, int read_len, int name_len);
void ha_compress_base(uint8_t* dest, char* src, uint64_t src_l, uint64_t** N_site_lis, uint64_t N_site_occ);
void init_UC_Read(UC_Read* r);
void recover_UC_Read(UC_Read* r, All_reads* R_INF, uint64_t ID);
void recover_UC_Read(UC_Read* r, const All_reads *R_INF, uint64_t ID);
void recover_UC_Read_RC(UC_Read* r, All_reads* R_INF, uint64_t ID);
void recover_UC_Read_sub_region(char* r, long long start_pos, long long length, uint8_t strand, All_reads* R_INF, long long ID);
void destory_UC_Read(UC_Read* r);
@@ -201,12 +154,4 @@ void write_All_reads(All_reads* r, char* read_file_name);
int load_All_reads(All_reads* r, char* read_file_name);
void destory_All_reads(All_reads* r);
void destory_R_buffer_block(R_buffer_block* curr_sub_block);
void destory_R_buffer();
void clear_R_buffer();
void init_gz_files(hifiasm_opt_t* asm_opt);
void destory_gz_files();
#endif
Purge_Dups.cpp
+83 -61
View File
@@ -1,10 +1,10 @@
#include "Purge_Dups.h"
#include "Overlaps.h"
#include "Correct.h"
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#include <stdlib.h>
#include "ksort.h"
#include "Purge_Dups.h"
#include "Overlaps.h"
#include "Correct.h"
#define Cal_Off(OFF) ((long long)((uint32_t)((OFF)>>32)) - (long long)((uint32_t)((OFF))))
#define Get_match(x) ((x).weight)
@@ -76,6 +76,37 @@ typedef struct {
uint32_t num;
}hap_overlaps_list;
typedef struct {
uint64_t* vote_counting;
uint8_t* visit;
kvec_t_u64_warp u_vecs;
kvec_asg_arc_t_offset u_buffer;
kvec_hap_candidates u_can;
}hap_alignment_struct;
void init_hap_alignment_struct(hap_alignment_struct* x, uint32_t size)
{
x->vote_counting = (uint64_t*)malloc(sizeof(uint64_t)*size);
memset(x->vote_counting, 0, sizeof(uint64_t)*size);
x->visit = (uint8_t*)malloc(sizeof(uint8_t)*size);
memset(x->visit, 0, size);
kv_init(x->u_vecs.a);
kv_init(x->u_buffer.a);
kv_init(x->u_can.a);
}
void destory_hap_alignment_struct(hap_alignment_struct* x, uint32_t size)
{
free(x->vote_counting);
free(x->visit);
kv_destroy(x->u_vecs.a);
kv_destroy(x->u_buffer.a);
kv_destroy(x->u_can.a);
}
void init_hap_overlaps_list(hap_overlaps_list* x, uint32_t num)
{
uint32_t i = 0;
@@ -205,7 +236,20 @@ uint32_t xUnitigLen, uint32_t yUnitigLen, uint64_t* position_index, uint8_t* rev
}
(*rev) = x_dir^y_dir;
if((*rev)) y_pos = yUnitigLen - y_pos - 1;
if((*rev))
{
if(yUnitigLen <= y_pos)
{
y_pos = (uint32_t)-1;
}
else
{
y_pos = yUnitigLen - y_pos - 1;
}
}
if(x_pos>=xUnitigLen) x_pos = (uint32_t)-1;
if(y_pos>=yUnitigLen) y_pos = (uint32_t)-1;
tmp = x_pos; tmp = tmp << 32; tmp = tmp | y_pos;
return tmp;
@@ -439,6 +483,11 @@ long long targetEnd, long long targetID, long long eMatch, long long eTotal, lon
float Hap_rate, uint64_t* position_index, ma_hit_t_alloc* reverse_sources, asg_t *read_g,
R_to_U* ruIndex, uint32_t* n_matchLen, uint32_t* n_max_count, uint32_t* n_min_count)
{
if(queryLen == 0)
{
(*n_matchLen) = (*n_min_count) = (*n_max_count) = 0;
return;
}
long long i, maxId, min_count = eTotal, max_count = eMatch, matchLen = 0;
uint32_t is_found, is_match;
if(dir == 0)
@@ -577,7 +626,6 @@ long long* r_y_interval_beg, long long* r_y_interval_end)
long long x_interval_beg, x_interval_end, y_interval_beg, y_interval_end;
long long target_beg, target_end;
x_max_count = x_min_count = y_max_count = y_min_count = 0;
if(type == X2Y)
{
/********************x*********************/
@@ -675,7 +723,7 @@ long long* r_y_interval_beg, long long* r_y_interval_end)
xUid, 0, xReads->n - 1, Hap_rate, position_index, reverse_sources, read_g, ruIndex, rev,
&y_interval_beg, &y_interval_end);
/********************y*********************/
}
} else abort();
(*r_x_interval_beg) = x_interval_beg;
(*r_x_interval_end) = x_interval_end;
@@ -768,17 +816,17 @@ uint32_t* x_off, uint32_t* y_off)
w_dir = (t.v == w)?1:0;
if(rev == 0 && v_dir != w_dir) continue;
if(rev == 1 && v_dir == w_dir) continue;
tmp = get_xy_pos(read_g, &t, v, w, xReads->len, yReads->len, position_index, &(t.el));
if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue;
///if(is_found == 0 || ((uint32_t)(tmp>>32) > (*x_off) && ((uint32_t)tmp) > (*y_off)))
if(is_found == 0 || t.ol > oLen)
{
(*x_off) = tmp>>32;
(*y_off) = (uint32_t)tmp;
oLen = t.ol;
}
}
is_found = 1;
}
@@ -828,13 +876,15 @@ uint32_t* x_off, uint32_t* y_off)
tmp = get_xy_pos(read_g, &t, v, w, xReads->len, yReads->len, position_index, &(t.el));
if(((tmp>>32) == (uint32_t)-1) || (((uint32_t)tmp) == (uint32_t)-1)) continue;
///if(is_found == 0 || ((uint32_t)(tmp>>32) < (*x_off) && ((uint32_t)tmp) < (*y_off)))
if(is_found == 0 || t.ol > oLen)
{
(*x_off) = tmp>>32;
(*y_off) = (uint32_t)tmp;
oLen = t.ol;
}
}
is_found = 1;
}
@@ -868,15 +918,7 @@ long long* r_x_pos_beg, long long* r_x_pos_end, long long* r_y_pos_beg, long lon
return (uint32_t)-1;
}
if(x_pos_beg > x_pos_end || y_pos_beg > y_pos_end) return (uint32_t)-1;
/**
fprintf(stderr, "\nrev: %u, weight: %lu\n", Get_rev(*hap_can), Get_match(*hap_can));
fprintf(stderr, "xUid: %u, xLen: %u, xBase: %u, x_interval_beg: %u, x_interval_end: %u, x_pos_beg: %u, x_pos_end: %u\n",
xUid, xReads->n, xReads->len, Get_x_beg(*hap_can), Get_x_end(*hap_can), x_pos_beg, x_pos_end);
fprintf(stderr, "yUid: %u, yLen: %u, yBase: %u, y_interval_beg: %u, y_interval_end: %u, y_pos_beg: %u, y_pos_end: %u\n",
yUid, yReads->n, yReads->len, Get_y_beg(*hap_can), Get_y_end(*hap_can), y_pos_beg, y_pos_end);
**/
/**
#define X2Y 0
#define Y2X 1
@@ -900,8 +942,8 @@ long long* r_y_pos_beg, long long* r_y_pos_end)
uint32_t xLeftBeg, xLeftLen, yLeftBeg, yLeftLen;
uint32_t xRightBeg, xRightLen, yRightBeg, yRightLen;
uint64_t totalWeigth;
double xLeftMatch, xLeftTotal, yLeftMatch, yLeftTotal;
double xRightMatch, xRightTotal, yRightMatch, yRightTotal;
double xLeftMatch = 0, xLeftTotal = 0, yLeftMatch = 0, yLeftTotal = 0;
double xRightMatch = 0, xRightTotal = 0, yRightMatch = 0, yRightTotal = 0;
asg_arc_t_offset* arch = NULL;
for (i = hap_can->index_beg, totalWeigth = 0; i <= hap_can->index_end; i++)
@@ -964,7 +1006,7 @@ long long* r_y_pos_beg, long long* r_y_pos_end)
&yRightMatch, &yRightTotal);
max_count = xLeftMatch + yRightMatch;
min_count = xLeftTotal + yRightTotal;
}
} else abort();
hap_can->weight = hap_can->index_beg = 0;
@@ -972,6 +1014,7 @@ long long* r_y_pos_beg, long long* r_y_pos_end)
if(max_count > min_count*Hap_rate)
{
long long r_x_interval_beg, r_x_interval_end, r_y_interval_beg, r_y_interval_end;
///for containment, don't need to do anything
get_hap_alignment_boundary(xReads, yReads, flag, xLeftMatch, xLeftTotal,
yLeftMatch, yLeftTotal, xRightMatch, xRightTotal, yRightMatch, yRightTotal,
@@ -1007,6 +1050,14 @@ long long* r_y_pos_beg, long long* r_y_pos_end)
}
void print_hap_paf(ma_ug_t *ug, hap_overlaps* ovlp)
{
fprintf(stderr, "utg%.6d%c\t%u\t%u\t%u\t%c\tutg%.6d%c\t%u\t%u\t%u\t%u\t%u\n",
ovlp->xUid+1, "lc"[ug->u.a[ovlp->xUid].circ], ug->u.a[ovlp->xUid].len, ovlp->x_beg_pos, ovlp->x_end_pos, "+-"[ovlp->rev],
ovlp->yUid+1, "lc"[ug->u.a[ovlp->yUid].circ], ug->u.a[ovlp->yUid].len, ovlp->y_beg_pos, ovlp->y_end_pos, ovlp->type, (uint32_t)ovlp->weight);
}
void hap_alignment(ma_ug_t *ug, asg_t *read_g, ma_hit_t_alloc* reverse_sources,
R_to_U* ruIndex, ma_sub_t *coverage_cut, uint64_t* position_index, uint64_t* vote_counting,
uint8_t* visit, kvec_t_u64_warp* u_vecs, kvec_asg_arc_t_offset* u_buffer, kvec_hap_candidates* u_can,
@@ -1142,19 +1193,12 @@ hap_overlaps_list* all_ovlp)
t_offset.Off = get_xy_pos(read_g, &t, v, (yReads->a[(uint32_t)(position_index[rId])])>>32,
xReads->len, yReads->len, position_index, &(t.el));
if(((t_offset.Off>>32) == (uint32_t)-1) || (((uint32_t)t_offset.Off) == (uint32_t)-1)) continue;
t_offset.x = t;
t_offset.weight = 1;
kv_push(asg_arc_t_offset, u_buffer->a, t_offset);
// if(debug_enable)
// {
// fprintf(stderr, "xUid: %u, yUid: %u, (%u), x_index: %u, y_index: %u, dis: %lld\n\n",
// xUid, yUid, (uint32_t)(u_buffer->a.n-1),
// (uint32_t)(position_index[u_buffer->a.a[u_buffer->a.n-1].x.ul>>33]),
// (uint32_t)(position_index[u_buffer->a.a[u_buffer->a.n-1].x.v>>1]),
// (long long)((uint32_t)(u_buffer->a.a[u_buffer->a.n-1].Off>>32)) -
// (long long)((uint32_t)(u_buffer->a.a[u_buffer->a.n-1].Off)));
// }
}
deduplicate_edge(u_buffer);
@@ -1296,36 +1340,10 @@ hap_overlaps_list* all_ovlp)
if(Get_match(hap_can) == 0 || Get_total(hap_can) == 0) continue;
kv_push(hap_overlaps, all_ovlp->x[hap_align.xUid].a, hap_align);
/**
fprintf(stderr, "\nsplit\nxUid: %u, yUid: %u, u_buffer->a.n: %u\n", xUid, yUid, (uint32_t)u_buffer->a.n);
for (k = 0; k < u_buffer->a.n; k++)
{
v = u_buffer->a.a[k].x.ul>>33;
w = u_buffer->a.a[k].x.v>>1;
fprintf(stderr, "\n(%u), x_index: %u, x_pos: %u, x_real_pos: %u, dis: %lld, dir: %u, weight: %lu\n",
k, (uint32_t)(position_index[v]), (uint32_t)(position_index[v]>>32),
(uint32_t)(u_buffer->a.a[k].Off>>32), Cal_Off(u_buffer->a.a[k].Off),
u_buffer->a.a[k].x.el, u_buffer->a.a[k].weight);
fprintf(stderr, "(%u), y_index: %u, y_pos: %u, y_real_pos: %u\n", k,
(uint32_t)(position_index[w]), (uint32_t)(position_index[w]>>32),
(uint32_t)(u_buffer->a.a[k].Off));
}
**/
}
}
void print_hap_paf(ma_ug_t *ug, hap_overlaps* ovlp)
{
fprintf(stderr, "utg%.6d%c\t%u\t%u\t%u\t%c\tutg%.6d%c\t%u\t%u\t%u\t%u\t%u\n",
ovlp->xUid+1, "lc"[ug->u.a[ovlp->xUid].circ], ug->u.a[ovlp->xUid].len, ovlp->x_beg_pos, ovlp->x_end_pos, "+-"[ovlp->rev],
ovlp->yUid+1, "lc"[ug->u.a[ovlp->yUid].circ], ug->u.a[ovlp->yUid].len, ovlp->y_beg_pos, ovlp->y_end_pos, ovlp->type, (uint32_t)ovlp->weight);
}
int inline get_specific_hap_overlap(kvec_hap_overlaps* x, uint32_t qn, uint32_t tn)
{
@@ -2014,6 +2032,7 @@ R_to_U* ruIndex, kvec_asg_arc_t_warp* edge, float density, uint32_t bi_graph_Len
float lable_match_rate, int max_hang, int min_ovlp, long long bubble_dist, float drop_ratio,
uint32_t just_contain)
{
fprintf(stderr, "*****************\n");
asg_t *purge_g = NULL;
purge_g = asg_init();
kvec_t_u64_warp u_vecs;
@@ -2140,9 +2159,10 @@ uint32_t just_contain)
purge_g->seq[all_ovlp.x[uId].a.a[i].xUid].del||
purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].c == ALTER_LABLE||
purge_g->seq[all_ovlp.x[uId].a.a[i].yUid].del)
{
continue;
}
{
continue;
}
///print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
r = get_hap_arch(&(all_ovlp.x[uId].a.a[i]), ug->u.a[all_ovlp.x[uId].a.a[i].xUid].len,
ug->u.a[all_ovlp.x[uId].a.a[i].yUid].len, max_hang, asm_opt.max_hang_rate, min_ovlp, &t);
@@ -2156,7 +2176,8 @@ uint32_t just_contain)
else
{
print_hap_paf(ug, &(all_ovlp.x[uId].a.a[i]));
fprintf(stderr, "error\n");
fprintf(stderr, "error: uId: %u, i: %u, xUid: %u, yUid: %u\n",
uId, i, all_ovlp.x[uId].a.a[i].xUid, all_ovlp.x[uId].a.a[i].yUid);
}
}
}
@@ -2197,5 +2218,6 @@ uint32_t just_contain)
free(position_index);
free(vote_counting);
free(visit);
fprintf(stderr, "#################\n");
}
Trio.cpp
+56 -129
View File
@@ -6,68 +6,35 @@
#include "khashl.h" // hash table
#include "kthread.h"
#include "Process_Read.h"
#include "Trio.h"
#include "htab.h"
#include "CommandLines.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 // yak uses 10, but hifiasm uses 12; we have to copy over some yak code here due to this
#define YAK_N_COUNTS (1<<YAK_COUNTER_BITS)
#define YAK_MAX_COUNT ((1<<YAK_COUNTER_BITS)-1)
#define YAK_LOAD_ALL 1
#define YAK_LOAD_TRIOBIN1 2
#define YAK_LOAD_TRIOBIN2 3
#define YAK_MAGIC "YAK\2"
#define yak_ch_eq(a, b) ((a)>>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_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]);
}
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<<h->pre) - 1;
yak_ht_t *g = h->h[x&mask].h;
@@ -76,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;
@@ -100,32 +53,10 @@ yak_ch_t *yak_ch_init(int k, int pre, int n_hash, int n_shift)
CALLOC(h->h, 1<<h->pre);
for (i = 0; i < 1<<h->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<<h->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<<h->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;
@@ -161,13 +92,13 @@ 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<<ch->pre; ++i) {
yak_ht_t *h = ch->h[i].h;
fread(t, 4, 2, fp);
if (ch0 == 0) yak_ht_resize(h, t[0]);
for (j = 0; j < t[1]; ++j) {
for (j = 0; j < (int)t[1]; ++j) {
uint64_t key;
fread(&key, 8, 1, fp);
if (mode == YAK_LOAD_ALL) {
@@ -195,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<<h->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)
{
@@ -229,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);
@@ -266,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
@@ -282,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);
@@ -312,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;
}
@@ -326,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]);
}
@@ -361,4 +288,4 @@ void trio_partition()
fprintf(stderr, "Trio binning has been done.\n");
fprintf(stderr, "%-30s%18.2f\n\n", "Trio binning time:", Get_T() - start_time);
}
}
Trio.h
-71
View File
@@ -1,71 +0,0 @@
#ifndef __TRIO__
#define __TRIO__
#include <stdint.h>
#define YAK_MAX_KMER 31
#define YAK_COUNTER_BITS 10
#define YAK_N_COUNTS (1<<YAK_COUNTER_BITS)
#define YAK_MAX_COUNT ((1<<YAK_COUNTER_BITS)-1)
#define YAK_BLK_SHIFT 9 // 64 bytes, the size of a cache line
#define YAK_BLK_MASK ((1<<(YAK_BLK_SHIFT)) - 1)
#define YAK_LOAD_ALL 1
#define YAK_LOAD_TRIOBIN1 2
#define YAK_LOAD_TRIOBIN2 3
#define YAK_MAGIC "YAK\2"
typedef struct {
int n_shift, n_hashes;
uint8_t *b;
} yak_bf_t;
struct yak_ht_t;
typedef struct {
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;
typedef struct {
int max;
uint32_t *s;
} tb_buf_t;
typedef struct {
int k, n_threads, print_diff;
double ratio_thres;
///bseq_file_t *fp;
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;
///bseq1_t *seq;
///All_reads* seq;
///tb_cnt_t *cnt;
} tb_step_t;
void trio_partition();
#endif
anchor.cpp
+158
View File
@@ -0,0 +1,158 @@
#include <stdio.h>
#include "htab.h"
#include "ksort.h"
#include "Hash_Table.h"
typedef struct { // this struct is not strictly necessary; we can use k_mer_pos instead, with modifications
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)
#define oreg_xs_lt(a, b) (((uint64_t)(a).x_pos_s<<32|(a).x_pos_e) < ((uint64_t)(b).x_pos_s<<32|(b).x_pos_e))
KSORT_INIT(or_xs, overlap_region, oreg_xs_lt)
#define oreg_ss_lt(a, b) ((a).shared_seed > (b).shared_seed) // in the decending order
KSORT_INIT(or_ss, overlap_region, oreg_ss_lt)
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);
}
uint64_t ha_abuf_mem(const ha_abuf_t *ab)
{
return ab->m_a * sizeof(anchor1_t) + ab->mz.m * (sizeof(ha_mz1_t) + sizeof(seed1_t)) + sizeof(ha_abuf_t);
}
void ha_get_new_candidates(ha_abuf_t *ab, int64_t rid, UC_Read *ucr, overlap_region_alloc *overlap_list, Candidates_list *cl, double bw_thres, int max_n_chain, int keep_whole_chain)
{
extern void *ha_flt_tab;
extern ha_pt_t *ha_idx;
uint32_t i;
uint64_t k, l;
// prepare
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;
// get the list of anchors
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;
kroundup64(ab->m_a);
REALLOC(ab->a, ab->m_a);
}
for (i = 0, k = 0; i < ab->mz.n; ++i) {
int 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 = z->rev == y->rev? 0 : 1;
an->other_off = y->pos;
an->self_off = rev? ucr->length - 1 - (z->pos + 1 - z->span) : z->pos;
an->srt = (uint64_t)y->rid<<33 | (uint64_t)rev<<32 | an->other_off;
}
}
// sort anchors
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;
}
}
// copy over to _cl_
if (ab->m_a >= (uint64_t)cl->size) {
cl->size = ab->m_a;
REALLOC(cl->list, cl->size);
}
for (k = 0; k < ab->n_a; ++k) {
k_mer_hit *p = &cl->list[k];
p->readID = ab->a[k].srt >> 33;
p->strand = ab->a[k].srt >> 32 & 1;
p->offset = ab->a[k].other_off;
p->self_offset = ab->a[k].self_off;
}
cl->length = ab->n_a;
calculate_overlap_region_by_chaining(cl, overlap_list, rid, ucr->length, &R_INF, bw_thres, keep_whole_chain);
#if 0
if (overlap_list->length > 2000) {
fprintf(stderr, "B\t%ld\t%ld\t%ld\n", (long)rid, (long)overlap_list->length, (long)Get_READ_LENGTH(R_INF, rid));
for (int i = 0; i < (int)overlap_list->length; ++i) {
overlap_region *r = &overlap_list->list[i];
fprintf(stderr, "C\t%d\t%d\t%d\t%c\t%d\t%ld\t%d\t%d\t%c\t%d\n", (int)r->x_id, (int)r->x_pos_s, (int)r->x_pos_e, "+-"[r->x_pos_strand],
(int)r->y_id, (long)Get_READ_LENGTH(R_INF, r->y_id), (int)r->y_pos_s, (int)r->y_pos_e, "+-"[r->y_pos_strand], (int)r->shared_seed);
}
}
#endif
if ((int)overlap_list->length > max_n_chain) {
uint32_t n[2], s[2];
n[0] = n[1] = 0, s[0] = s[1] = 0;
for (i = 0; i < (uint32_t)overlap_list->length; ++i) {
const overlap_region *r = &overlap_list->list[i];
int dir = r->x_pos_s == 0? 0 : 1;
++n[dir];
if ((int)n[dir] == max_n_chain) s[dir] = r->shared_seed;
}
if (s[0] > 0 || s[1] > 0) {
for (i = 0, k = 0; i < (uint32_t)overlap_list->length; ++i) {
overlap_region *r = &overlap_list->list[i];
int dir = r->x_pos_s == 0? 0 : 1;
if (r->shared_seed > s[dir]) {
if ((uint32_t)k != i) {
overlap_region t;
t = overlap_list->list[k];
overlap_list->list[k] = overlap_list->list[i];
overlap_list->list[i] = t;
}
++k;
}
}
overlap_list->length = k;
}
}
ks_introsort_or_xs(overlap_list->length, overlap_list->list);
}
hifiasm.1
+4
View File
@@ -202,7 +202,11 @@ but occurs <
times in the other sample.
.SS Debugging options
.TP 10
.B --dbg-gfa
Write additional files to speed up the debugging of graph cleaning
.SH EXAMPLES
hist.cpp
+92
View File
@@ -0,0 +1,92 @@
#include <stdio.h>
#include "htab.h"
static void ha_hist_line(int c, int x, int exceed, int64_t cnt)
{
int j;
if (c >= 0) fprintf(stderr, "[M::%s] %5d: ", __func__, c);
else fprintf(stderr, "[M::%s] %5s: ", __func__, "rest");
for (j = 0; j < x; ++j) fputc('*', stderr);
if (exceed) fputc('>', stderr);
fprintf(stderr, " %lld\n", (long long)cnt);
}
int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het)
{
const int hist_max = 100;
int i, start, low_i, max_i, max2_i, max3_i;
int64_t max, max2, max3, min;
// find the low point from the left
*peak_het = -1;
start = cnt[1] > 0? 1 : 2;
low_i = start;
for (i = low_i + 1; i < n_cnt; ++i)
if (cnt[i] > cnt[i-1]) break;
low_i = i - 1;
fprintf(stderr, "[M::%s] lowest: count[%d] = %ld\n", __func__, low_i, (long)cnt[low_i]);
if (low_i == n_cnt - 1) return -1; // low coverage
// find the highest peak
max_i = low_i + 1, max = cnt[max_i];
for (i = low_i + 1; i < n_cnt; ++i)
if (cnt[i] > max)
max = cnt[i], max_i = i;
fprintf(stderr, "[M::%s] highest: count[%d] = %ld\n", __func__, max_i, (long)cnt[max_i]);
// print histogram
for (i = start; i < n_cnt; ++i) {
int x, exceed = 0;
x = (int)((double)hist_max * cnt[i] / cnt[max_i] + .499);
if (x > hist_max) exceed = 1, x = hist_max; // may happen if cnt[2] is higher
if (i > max_i && x == 0) break;
ha_hist_line(i, x, exceed, cnt[i]);
}
{
int x, exceed = 0;
int64_t rest = 0;
for (; i < n_cnt; ++i) rest += cnt[i];
x = (int)((double)hist_max * rest / cnt[max_i] + .499);
if (x > hist_max) exceed = 1, x = hist_max;
ha_hist_line(-1, x, exceed, rest);
}
// look for smaller peak on the low end
max2 = -1; max2_i = -1;
for (i = max_i - 1; i > low_i; --i) {
if (cnt[i] >= cnt[i-1] && cnt[i] >= cnt[i+1]) {
if (cnt[i] > max2) max2 = cnt[i], max2_i = i;
}
}
if (max2_i > low_i && max2_i < max_i) {
for (i = max2_i + 1, min = max; i < max_i; ++i)
if (cnt[i] < min) min = cnt[i];
if (max2 < max * 0.05 || min > max2 * 0.95)
max2 = -1, max2_i = -1;
}
if (max2 > 0) fprintf(stderr, "[M::%s] left: count[%d] = %ld\n", __func__, max2_i, (long)cnt[max2_i]);
else fprintf(stderr, "[M::%s] left: none\n", __func__);
// look for smaller peak on the high end
max3 = -1; max3_i = -1;
for (i = max_i + 1; i < n_cnt - 1; ++i) {
if (cnt[i] >= cnt[i-1] && cnt[i] >= cnt[i+1]) {
if (cnt[i] > max3) max3 = cnt[i], max3_i = i;
}
}
if (max3_i > max_i) {
for (i = max_i + 1, min = max; i < max3_i; ++i)
if (cnt[i] < min) min = cnt[i];
if (max3 < max * 0.05 || min > max3 * 0.95 || max3_i > max_i * 2.5)
max3 = -1, max3_i = -1;
}
if (max3 > 0) fprintf(stderr, "[M::%s] right: count[%d] = %ld\n", __func__, max3_i, (long)cnt[max3_i]);
else fprintf(stderr, "[M::%s] right: none\n", __func__);
if (max3_i > 0) {
*peak_het = max_i;
return max3_i;
} else {
if (max2_i > 0) *peak_het = max2_i;
return max_i;
}
}
htab.cpp
+832
View File
@@ -0,0 +1,832 @@
#include <stdint.h>
#include <zlib.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include "kthread.h"
#include "khashl.h"
#include "kseq.h"
#include "ksort.h"
#include "htab.h"
#define YAK_COUNTER_BITS 12
#define YAK_N_COUNTS (1<<YAK_COUNTER_BITS)
#define YAK_MAX_COUNT ((1<<YAK_COUNTER_BITS)-1)
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
};
void *ha_flt_tab;
ha_pt_t *ha_idx;
/***************************
* Yak specific parameters *
***************************/
typedef struct {
int32_t bf_shift, bf_n_hash;
int32_t k, w, is_HPC;
int32_t pre;
int32_t n_thread;
int64_t chunk_size;
} yak_copt_t;
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->w = 1;
o->pre = YAK_COUNTER_BITS;
o->n_thread = 4;
o->chunk_size = 20000000;
}
/************************
* Blocked bloom filter *
************************/
#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;
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) - 1);
int h1 = hash >> 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;
}
/********************
* Count hash table *
********************/
#define yak_ct_eq(a, b) ((a)>>YAK_COUNTER_BITS == (b)>>YAK_COUNTER_BITS) // lower 8 bits for counts; higher bits for k-mer
#define yak_ct_hash(a) ((a)>>YAK_COUNTER_BITS)
KHASHL_SET_INIT(static klib_unused, yak_ct_t, yak_ct, uint64_t, yak_ct_hash, yak_ct_eq)
typedef struct {
yak_ct_t *h;
yak_bf_t *b;
} ha_ct1_t;
typedef struct {
int k, pre, n_hash, n_shift;
uint64_t tot;
ha_ct1_t *h;
} ha_ct_t;
static ha_ct_t *ha_ct_init(int k, int pre, int n_hash, int n_shift)
{
ha_ct_t *h;
int i;
if (pre < YAK_COUNTER_BITS) return 0;
CALLOC(h, 1);
h->k = k, h->pre = pre;
CALLOC(h->h, 1<<h->pre);
for (i = 0; i < 1<<h->pre; ++i)
h->h[i].h = yak_ct_init();
if (n_hash > 0 && n_shift > h->pre) {
h->n_hash = n_hash, h->n_shift = n_shift;
for (i = 0; i < 1<<h->pre; ++i)
h->h[i].b = yak_bf_init(h->n_shift - h->pre, h->n_hash);
}
return h;
}
static void ha_ct_destroy_bf(ha_ct_t *h)
{
int i;
for (i = 0; i < 1<<h->pre; ++i) {
if (h->h[i].b)
yak_bf_destroy(h->h[i].b);
h->h[i].b = 0;
}
}
static void ha_ct_destroy(ha_ct_t *h)
{
int i;
if (h == 0) return;
ha_ct_destroy_bf(h);
for (i = 0; i < 1<<h->pre; ++i)
yak_ct_destroy(h->h[i].h);
free(h->h); free(h);
}
static int ha_ct_insert_list(ha_ct_t *h, int create_new, int n, const uint64_t *a)
{
int j, mask = (1<<h->pre) - 1, n_ins = 0;
ha_ct1_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_ct_put(g->h, x << YAK_COUNTER_BITS | (g->b? 1 : 0), &absent);
if (absent) ++n_ins;
if ((kh_key(g->h, k)&YAK_MAX_COUNT) < YAK_MAX_COUNT)
++kh_key(g->h, k);
}
} else {
k = yak_ct_get(g->h, x<<YAK_COUNTER_BITS);
if (k != kh_end(g->h) && (kh_key(g->h, k)&YAK_MAX_COUNT) < YAK_MAX_COUNT)
++kh_key(g->h, k);
}
}
return n_ins;
}
/*** generate histogram ***/
typedef struct {
uint64_t c[YAK_N_COUNTS];
} buf_cnt_t;
typedef struct {
const ha_ct_t *h;
buf_cnt_t *cnt;
} hist_aux_t;
static void worker_ct_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_ct_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 ha_ct_hist(const ha_ct_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_ct_hist, &a, 1<<h->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;
ha_ct_t *h;
} shrink_aux_t;
static void worker_ct_shrink(void *data, long i, int tid) // callback for kt_for()
{
shrink_aux_t *a = (shrink_aux_t*)data;
ha_ct_t *h = a->h;
yak_ct_t *g = h->h[i].h, *f;
khint_t k;
f = yak_ct_init();
yak_ct_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_ct_put(f, kh_key(g, k), &absent);
}
}
yak_ct_destroy(g);
h->h[i].h = f;
}
static void ha_ct_shrink(ha_ct_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_ct_shrink, &a, 1<<h->pre);
for (i = 0, h->tot = 0; i < 1<<h->pre; ++i)
h->tot += kh_size(h->h[i].h);
}
/***********************
* Position hash table *
***********************/
KHASHL_MAP_INIT(static klib_unused, yak_pt_t, yak_pt, uint64_t, uint64_t, yak_ct_hash, yak_ct_eq)
#define generic_key(x) (x)
KRADIX_SORT_INIT(ha64, uint64_t, generic_key, 8)
typedef struct {
yak_pt_t *h;
uint64_t n;
ha_idxpos_t *a;
} ha_pt1_t;
struct ha_pt_s {
int k, pre;
uint64_t tot, tot_pos;
ha_pt1_t *h;
};
typedef struct {
const ha_ct_t *ct;
ha_pt_t *pt;
} pt_gen_aux_t;
static void worker_pt_gen(void *data, long i, int tid) // callback for kt_for()
{
pt_gen_aux_t *a = (pt_gen_aux_t*)data;
ha_pt1_t *b = &a->pt->h[i];
yak_ct_t *g = a->ct->h[i].h;
khint_t k;
for (k = 0, b->n = 0; k != kh_end(g); ++k) {
if (kh_exist(g, k)) {
int absent;
khint_t l;
l = yak_pt_put(b->h, kh_key(g, k) >> a->ct->pre << YAK_COUNTER_BITS, &absent);
kh_val(b->h, l) = b->n;
b->n += kh_key(g, k) & YAK_MAX_COUNT;
}
}
yak_ct_destroy(g);
a->ct->h[i].h = 0;
CALLOC(b->a, b->n);
}
ha_pt_t *ha_pt_gen(ha_ct_t *ct, int n_thread)
{
pt_gen_aux_t a;
int i;
ha_pt_t *pt;
ha_ct_destroy_bf(ct);
CALLOC(pt, 1);
pt->k = ct->k, pt->pre = ct->pre, pt->tot = ct->tot;
CALLOC(pt->h, 1<<pt->pre);
for (i = 0; i < 1<<pt->pre; ++i) {
pt->h[i].h = yak_pt_init();
yak_pt_resize(pt->h[i].h, kh_size(ct->h[i].h));
}
a.ct = ct, a.pt = pt;
kt_for(n_thread, worker_pt_gen, &a, 1<<pt->pre);
free(ct->h); free(ct);
return pt;
}
int ha_pt_insert_list(ha_pt_t *h, int n, const ha_mz1_t *a)
{
int j, mask = (1<<h->pre) - 1, n_ins = 0;
ha_pt1_t *g;
if (n == 0) return 0;
g = &h->h[a[0].x&mask];
for (j = 0; j < n; ++j) {
uint64_t x = a[j].x >> h->pre;
khint_t k;
int n;
ha_idxpos_t *p;
if ((a[j].x&mask) != (a[0].x&mask)) continue;
k = yak_pt_get(g->h, x<<YAK_COUNTER_BITS);
if (k == kh_end(g->h)) continue;
n = kh_key(g->h, k) & YAK_MAX_COUNT;
assert(n < YAK_MAX_COUNT);
p = &g->a[kh_val(g->h, k) + n];
p->rid = a[j].rid, p->rev = a[j].rev, p->pos = a[j].pos, p->span = a[j].span;
//(uint64_t)a[j].rid<<36 | (uint64_t)a[j].rev<<35 | (uint64_t)a[j].pos<<8 | (uint64_t)a[j].span;
++kh_key(g->h, k);
++n_ins;
}
return n_ins;
}
/*
static void worker_pt_sort(void *data, long i, int tid)
{
ha_pt_t *h = (ha_pt_t*)data;
ha_pt1_t *g = &h->h[i];
khint_t k;
for (k = 0; k < kh_end(g->h); ++k) {
int n;
uint64_t *p;
if (!kh_exist(g->h, k)) continue;
n = kh_key(g->h, k) & YAK_MAX_COUNT;
p = &g->a[kh_val(g->h, k)];
radix_sort_ha64(p, p + n);
}
}
void ha_pt_sort(ha_pt_t *h, int n_thread)
{
kt_for(n_thread, worker_pt_sort, h, 1<<h->pre);
}
*/
void ha_pt_destroy(ha_pt_t *h)
{
int i;
if (h == 0) return;
for (i = 0; i < 1<<h->pre; ++i) {
yak_pt_destroy(h->h[i].h);
free(h->h[i].a);
}
free(h->h); free(h);
}
const ha_idxpos_t *ha_pt_get(const ha_pt_t *h, uint64_t hash, int *n)
{
khint_t k;
const ha_pt1_t *g = &h->h[hash & ((1ULL<<h->pre) - 1)];
*n = 0;
k = yak_pt_get(g->h, hash >> h->pre << YAK_COUNTER_BITS);
if (k == kh_end(g->h)) return 0;
*n = kh_key(g->h, k) & YAK_MAX_COUNT;
return &g->a[kh_val(g->h, k)];
}
/**********************************
* Buffer for counting all k-mers *
**********************************/
typedef struct {
int n, m;
uint64_t n_ins;
uint64_t *a;
ha_mz1_t *b;
} ch_buf_t;
static inline void ct_insert_buf(ch_buf_t *buf, int p, uint64_t y) // insert a k-mer $y to a linear buffer
{
int pre = y & ((1<<p) - 1);
ch_buf_t *b = &buf[pre];
if (b->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 inline void pt_insert_buf(ch_buf_t *buf, int p, const ha_mz1_t *y)
{
int pre = y->x & ((1<<p) - 1);
ch_buf_t *b = &buf[pre];
if (b->n == b->m) {
b->m = b->m < 8? 8 : b->m + (b->m>>1);
REALLOC(b->b, b->m);
}
b->b[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<<k) - 1, shift = k - 1;
for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) {
int c = seq_nt4_table[(uint8_t)seq[i]];
if (c < 4) { // not an "N" base
x[0] = (x[0] << 1 | (c&1)) & mask;
x[1] = (x[1] << 1 | (c>>1)) & mask;
x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift;
x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift;
if (++l >= k)
ct_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
}
}
static void count_seq_buf_HPC(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, last = -1;
uint64_t x[4], mask = (1ULL<<k) - 1, shift = k - 1;
for (i = l = 0, x[0] = x[1] = x[2] = x[3] = 0; i < len; ++i) {
int c = seq_nt4_table[(uint8_t)seq[i]];
if (c < 4) { // not an "N" base
if (c != last) {
x[0] = (x[0] << 1 | (c&1)) & mask;
x[1] = (x[1] << 1 | (c>>1)) & mask;
x[2] = x[2] >> 1 | (uint64_t)(1 - (c&1)) << shift;
x[3] = x[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift;
if (++l >= k)
ct_insert_buf(buf, p, yak_hash_long(x));
last = c;
}
} else l = 0, last = -1, x[0] = x[1] = x[2] = x[3] = 0; // if there is an "N", restart
}
}
/******************
* K-mer counting *
******************/
KSEQ_INIT(gzFile, gzread)
#define HAF_COUNT_EXACT 0x1
#define HAF_COUNT_ALL 0x2
#define HAF_RS_WRITE_LEN 0x4
#define HAF_RS_WRITE_SEQ 0x8
#define HAF_RS_READ 0x10
typedef struct { // global data structure for kt_pipeline()
const yak_copt_t *opt;
const void *flt_tab;
int flag, create_new, is_store;
uint64_t n_seq;
kseq_t *ks;
UC_Read ucr;
ha_ct_t *ct;
ha_pt_t *pt;
const All_reads *rs_in;
All_reads *rs_out;
} pl_data_t;
typedef struct { // data structure for each step in kt_pipeline()
pl_data_t *p;
uint64_t n_seq0;
int n_seq, m_seq, sum_len, nk;
int *len;
char **seq;
ha_mz1_v *mz_buf;
ha_mz1_v *mz;
ch_buf_t *buf;
} st_data_t;
static void worker_for_insert(void *data, long i, int tid) // callback for kt_for()
{
st_data_t *s = (st_data_t*)data;
ch_buf_t *b = &s->buf[i];
if (s->p->pt)
b->n_ins += ha_pt_insert_list(s->p->pt, b->n, b->b);
else
b->n_ins += ha_ct_insert_list(s->p->ct, s->p->create_new, b->n, b->a);
}
static void worker_for_mz(void *data, long i, int tid)
{
st_data_t *s = (st_data_t*)data;
ha_mz1_v *b = &s->mz_buf[tid];
s->mz_buf[tid].n = 0;
ha_sketch(s->seq[i], s->len[i], s->p->opt->w, s->p->opt->k, s->n_seq0 + i, s->p->opt->is_HPC, b, s->p->flt_tab);
s->mz[i].n = s->mz[i].m = b->n;
MALLOC(s->mz[i].a, b->n);
memcpy(s->mz[i].a, b->a, b->n * sizeof(ha_mz1_t));
}
static void *worker_count(void *data, int step, void *in) // callback for kt_pipeline()
{
pl_data_t *p = (pl_data_t*)data;
if (step == 0) { // step 1: read a block of sequences
int ret;
st_data_t *s;
CALLOC(s, 1);
s->p = p;
s->n_seq0 = p->n_seq;
if (p->rs_in && (p->flag & HAF_RS_READ)) {
while (p->n_seq < p->rs_in->total_reads) {
int l;
recover_UC_Read(&p->ucr, p->rs_in, p->n_seq);
l = p->ucr.length;
if (s->n_seq == s->m_seq) {
s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
REALLOC(s->len, s->m_seq);
REALLOC(s->seq, s->m_seq);
}
MALLOC(s->seq[s->n_seq], l);
memcpy(s->seq[s->n_seq], p->ucr.seq, l);
s->len[s->n_seq++] = l;
++p->n_seq;
s->sum_len += l;
s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
if (s->sum_len >= p->opt->chunk_size)
break;
}
} else {
while ((ret = kseq_read(p->ks)) >= 0) {
int l = p->ks->seq.l;
if (p->n_seq >= 1<<28) {
fprintf(stderr, "ERROR: this implementation supports no more than %d reads\n", 1<<28);
exit(1);
}
if (p->rs_out) {
if (p->flag & HAF_RS_WRITE_LEN) {
assert(p->n_seq == p->rs_out->total_reads);
ha_insert_read_len(p->rs_out, l, p->ks->name.l);
} else if (p->flag & HAF_RS_WRITE_SEQ) {
int i, n_N;
assert(l == (int)p->rs_out->read_length[p->n_seq]);
for (i = n_N = 0; i < l; ++i) // count number of ambiguous bases
if (seq_nt4_table[(uint8_t)p->ks->seq.s[i]] >= 4)
++n_N;
ha_compress_base(Get_READ(*p->rs_out, p->n_seq), p->ks->seq.s, l, &p->rs_out->N_site[p->n_seq], n_N);
memcpy(&p->rs_out->name[p->rs_out->name_index[p->n_seq]], p->ks->name.s, p->ks->name.l);
}
}
if (s->n_seq == s->m_seq) {
s->m_seq = s->m_seq < 16? 16 : s->m_seq + (s->m_seq>>1);
REALLOC(s->len, s->m_seq);
REALLOC(s->seq, s->m_seq);
}
MALLOC(s->seq[s->n_seq], l);
memcpy(s->seq[s->n_seq], p->ks->seq.s, l);
s->len[s->n_seq++] = l;
++p->n_seq;
s->sum_len += l;
s->nk += l >= p->opt->k? l - p->opt->k + 1 : 0;
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
st_data_t *s = (st_data_t*)in;
int i, n_pre = 1<<p->opt->pre, m;
// allocate the k-mer buffer
CALLOC(s->buf, n_pre);
m = (int)(s->nk * 1.2 / n_pre) + 1;
for (i = 0; i < n_pre; ++i) {
s->buf[i].m = m;
if (p->pt) MALLOC(s->buf[i].b, m);
else MALLOC(s->buf[i].a, m);
}
// fill the buffer
if (p->opt->w == 1) { // enumerate all k-mers
for (i = 0; i < s->n_seq; ++i) {
if (p->opt->is_HPC)
count_seq_buf_HPC(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
else
count_seq_buf(s->buf, p->opt->k, p->opt->pre, s->len[i], s->seq[i]);
if (!p->is_store) free(s->seq[i]);
}
} else { // minimizers only
uint32_t j;
// compute minimizers
CALLOC(s->mz, s->n_seq);
CALLOC(s->mz_buf, p->opt->n_thread);
kt_for(p->opt->n_thread, worker_for_mz, s, s->n_seq);
for (i = 0; i < p->opt->n_thread; ++i)
free(s->mz_buf[i].a);
free(s->mz_buf);
// insert minimizers
if (p->pt) {
for (i = 0; i < s->n_seq; ++i)
for (j = 0; j < s->mz[i].n; ++j)
pt_insert_buf(s->buf, p->opt->pre, &s->mz[i].a[j]);
} else {
for (i = 0; i < s->n_seq; ++i)
for (j = 0; j < s->mz[i].n; ++j)
ct_insert_buf(s->buf, p->opt->pre, s->mz[i].a[j].x);
}
for (i = 0; i < s->n_seq; ++i) {
free(s->mz[i].a);
if (!p->is_store) free(s->seq[i]);
}
free(s->mz);
}
free(s->seq); free(s->len);
s->seq = 0, s->len = 0;
return s;
} else if (step == 2) { // step 3: insert k-mers to hash table
st_data_t *s = (st_data_t*)in;
int i, n = 1<<p->opt->pre;
uint64_t n_ins = 0;
kt_for(p->opt->n_thread, worker_for_insert, s, n);
for (i = 0; i < n; ++i) {
n_ins += s->buf[i].n_ins;
if (p->pt) free(s->buf[i].b);
else free(s->buf[i].a);
}
if (p->ct) p->ct->tot += n_ins;
if (p->pt) p->pt->tot_pos += n_ins;
free(s->buf);
#if 0
fprintf(stderr, "[M::%s::%.3f*%.2f] processed %ld sequences; %ld %s in the hash table\n", __func__,
yak_realtime(), yak_cpu_usage(), (long)s->n_seq0 + s->n_seq,
(long)(p->pt? p->pt->tot_pos : p->ct->tot), p->pt? "positions" : "distinct k-mers");
#endif
free(s);
}
return 0;
}
static ha_ct_t *yak_count(const yak_copt_t *opt, const char *fn, int flag, ha_pt_t *p0, ha_ct_t *c0, const void *flt_tab, All_reads *rs)
{
int read_rs = (rs && (flag & HAF_RS_READ));
pl_data_t pl;
gzFile fp = 0;
memset(&pl, 0, sizeof(pl_data_t));
if (read_rs) {
pl.rs_in = rs;
init_UC_Read(&pl.ucr);
} else {
if ((fp = gzopen(fn, "r")) == 0) return 0;
pl.ks = kseq_init(fp);
}
if (rs && (flag & (HAF_RS_WRITE_LEN|HAF_RS_WRITE_SEQ)))
pl.rs_out = rs;
pl.flt_tab = flt_tab;
pl.opt = opt;
pl.flag = flag;
if (p0) {
pl.pt = p0, pl.create_new = 0;
assert(p0->k == opt->k && p0->pre == opt->pre);
} else if (c0) {
pl.ct = c0, pl.create_new = 0;
assert(c0->k == opt->k && c0->pre == opt->pre);
} else {
pl.create_new = 1;
pl.ct = ha_ct_init(opt->k, opt->pre, opt->bf_n_hash, opt->bf_shift);
}
kt_pipeline(3, worker_count, &pl, 3);
if (read_rs) {
destory_UC_Read(&pl.ucr);
} else {
kseq_destroy(pl.ks);
gzclose(fp);
}
return pl.ct;
}
ha_ct_t *ha_count(const hifiasm_opt_t *asm_opt, int flag, ha_pt_t *p0, const void *flt_tab, All_reads *rs)
{
int i;
yak_copt_t opt;
ha_ct_t *h = 0;
assert(!(flag & HAF_RS_WRITE_LEN) || !(flag & HAF_RS_WRITE_SEQ)); // not both
if (rs) {
if (flag & HAF_RS_WRITE_LEN)
init_All_reads(rs);
else if (flag & HAF_RS_WRITE_SEQ)
malloc_All_reads(rs);
}
yak_copt_init(&opt);
opt.k = asm_opt->k_mer_length;
opt.is_HPC = !asm_opt->no_HPC;
opt.w = flag & HAF_COUNT_ALL? 1 : asm_opt->mz_win;
opt.bf_shift = flag & HAF_COUNT_EXACT? 0 : asm_opt->bf_shift;
opt.n_thread = asm_opt->thread_num;
for (i = 0; i < asm_opt->num_reads; ++i)
h = yak_count(&opt, asm_opt->read_file_names[i], flag, p0, h, flt_tab, rs);
if (h && opt.bf_shift > 0)
ha_ct_destroy_bf(h);
return h;
}
/***************************
* High count filter table *
***************************/
KHASHL_SET_INIT(static klib_unused, yak_ft_t, yak_ft, uint64_t, kh_hash_dummy, kh_eq_generic)
static yak_ft_t *gen_hh(const ha_ct_t *h)
{
int i;
yak_ft_t *hh;
hh = yak_ft_init();
yak_ft_resize(hh, h->tot * 2);
for (i = 0; i < 1<<h->pre; ++i) {
yak_ct_t *ht = h->h[i].h;
khint_t k;
for (k = 0; k < kh_end(ht); ++k) {
if (kh_exist(ht, k)) {
uint64_t y = kh_key(ht, k) >> h->pre << YAK_COUNTER_BITS | i;
int absent;
yak_ft_put(hh, y, &absent);
}
}
}
return hh;
}
int ha_ft_isflt(const void *hh, uint64_t y)
{
yak_ft_t *h = (yak_ft_t*)hh;
khint_t k;
k = yak_ft_get(h, y);
return k == kh_end(h)? 0 : 1;
}
void ha_ft_destroy(void *h)
{
if (h) yak_ft_destroy((yak_ft_t*)h);
}
/*************************
* High-level interfaces *
*************************/
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs)
{
yak_ft_t *flt_tab;
int64_t cnt[YAK_N_COUNTS];
int peak_hom, peak_het, cutoff;
ha_ct_t *h;
h = ha_count(asm_opt, HAF_COUNT_ALL|HAF_RS_WRITE_LEN, NULL, NULL, rs);
ha_ct_hist(h, cnt, asm_opt->thread_num);
peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
cutoff = (int)(peak_hom * asm_opt->high_factor);
if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
ha_ct_shrink(h, cutoff, YAK_MAX_COUNT, asm_opt->thread_num);
flt_tab = gen_hh(h);
ha_ct_destroy(h);
fprintf(stderr, "[M::%s::%.3f*%.2f@%.3fGB] ==> filtered out %ld k-mers occurring %d or more times\n", __func__,
yak_realtime(), yak_cpu_usage(), yak_peakrss_in_gb(), (long)kh_size(flt_tab), cutoff);
return (void*)flt_tab;
}
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs)
{
int64_t cnt[YAK_N_COUNTS], tot_cnt;
int peak_hom, peak_het, i, extra_flag1, extra_flag2;
ha_ct_t *ct;
ha_pt_t *pt;
if (read_from_store) {
extra_flag1 = extra_flag2 = HAF_RS_READ;
} else if (rs->total_reads == 0) {
extra_flag1 = HAF_RS_WRITE_LEN;
extra_flag2 = HAF_RS_WRITE_SEQ;
} else {
extra_flag1 = HAF_RS_WRITE_SEQ;
extra_flag2 = HAF_RS_READ;
}
ct = ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag1, NULL, flt_tab, rs);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> counted %ld distinct minimizer k-mers\n", __func__,
yak_realtime(), yak_cpu_usage(), (long)ct->tot);
ha_ct_hist(ct, cnt, asm_opt->thread_num);
fprintf(stderr, "[M::%s] count[%d] = %ld (for sanity check)\n", __func__, YAK_MAX_COUNT, (long)cnt[YAK_MAX_COUNT]);
peak_hom = ha_analyze_count(YAK_N_COUNTS, cnt, &peak_het);
if (peak_hom > 0) fprintf(stderr, "[M::%s] peak_hom: %d; peak_het: %d\n", __func__, peak_hom, peak_het);
if (flt_tab == 0) {
int cutoff = (int)(peak_hom * asm_opt->high_factor);
if (cutoff > YAK_MAX_COUNT - 1) cutoff = YAK_MAX_COUNT - 1;
ha_ct_shrink(ct, 2, cutoff, asm_opt->thread_num);
for (i = 2, tot_cnt = 0; i <= cutoff; ++i) tot_cnt += cnt[i] * i;
} else {
ha_ct_shrink(ct, 2, YAK_MAX_COUNT - 1, asm_opt->thread_num);
for (i = 2, tot_cnt = 0; i <= YAK_MAX_COUNT - 1; ++i) tot_cnt += cnt[i] * i;
}
pt = ha_pt_gen(ct, asm_opt->thread_num);
ha_count(asm_opt, HAF_COUNT_EXACT|extra_flag2, pt, flt_tab, rs);
assert((uint64_t)tot_cnt == pt->tot_pos);
//ha_pt_sort(pt, asm_opt->thread_num);
fprintf(stderr, "[M::%s::%.3f*%.2f] ==> indexed %ld positions\n", __func__,
yak_realtime(), yak_cpu_usage(), (long)pt->tot_pos);
return pt;
}
htab.h
+103
View File
@@ -0,0 +1,103 @@
#ifndef __HA_HTAB_H__
#define __HA_HTAB_H__
#define __STDC_LIMIT_MACROS
#include <stdint.h>
#include "Process_Read.h"
#include "CommandLines.h"
typedef struct {
uint64_t x;
uint64_t rid:28, pos:27, rev:1, span:8;
} ha_mz1_t;
typedef struct {
uint64_t rid:28, pos:27, rev:1, span:8; // actually it is not necessary to keep span in the index
} ha_idxpos_t;
typedef struct { uint32_t n, m; ha_mz1_t *a; } ha_mz1_v;
struct ha_pt_s;
typedef struct ha_pt_s ha_pt_t;
struct ha_abuf_s;
typedef struct ha_abuf_s ha_abuf_t;
extern const unsigned char seq_nt4_table[256];
extern void *ha_flt_tab;
extern ha_pt_t *ha_idx;
void *ha_ft_gen(const hifiasm_opt_t *asm_opt, All_reads *rs);
int ha_ft_isflt(const void *hh, uint64_t y);
void ha_ft_destroy(void *h);
ha_pt_t *ha_pt_gen(const hifiasm_opt_t *asm_opt, const void *flt_tab, int read_from_store, All_reads *rs);
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);
uint64_t ha_abuf_mem(const ha_abuf_t *ab);
double yak_cputime(void);
void yak_reset_realtime(void);
double yak_realtime(void);
long yak_peakrss(void);
double yak_peakrss_in_gb(void);
double yak_cpu_usage(void);
void trio_partition(void);
void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf);
int ha_analyze_count(int n_cnt, const int64_t *cnt, int *peak_het);
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 kroundup32
#define kroundup32(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, ++(x))
#endif
#ifndef kroundup64
#define kroundup64(x) (--(x), (x)|=(x)>>1, (x)|=(x)>>2, (x)|=(x)>>4, (x)|=(x)>>8, (x)|=(x)>>16, x|=(x)>>32, ++(x))
#endif
#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__
ketopt.h
+1 -1
View File
@@ -17,7 +17,7 @@ typedef struct {
} ketopt_t;
typedef struct {
char *name;
const char *name;
int has_arg;
int val;
} ko_longopt_t;
khash.h
-669
View File
@@ -1,669 +0,0 @@
/* The MIT License
Copyright (c) 2008, 2009, 2011 by Attractive Chaos <attractor@live.co.uk>
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 <stdlib.h>
#include <string.h>
#include <limits.h>
#include <stdio.h>
/* 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);\
extern void kh_write_##name(kh_##name##_t *h, FILE* fp);\
extern void kh_load_##name(kh_##name##_t *h, FILE* fp);
#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; \
} \
} \
SCOPE void kh_write_##name(kh_##name##_t *h, FILE* fp)\
{\
fwrite(&(h->n_buckets), sizeof(khint_t), 1, fp);\
fwrite(&(h->size), sizeof(khint_t), 1, fp);\
fwrite(&(h->n_occupied), sizeof(khint_t), 1, fp);\
fwrite(&(h->upper_bound), sizeof(khint_t), 1, fp);\
if (h->n_buckets)\
{\
fwrite(h->flags, sizeof(khint32_t), __ac_fsize(h->n_buckets), fp);\
fwrite(h->keys, sizeof(khkey_t), h->n_buckets, fp);\
fwrite(h->vals, sizeof(khval_t), h->n_buckets, fp);\
}\
} \
SCOPE void kh_load_##name(kh_##name##_t *h, FILE* fp)\
{\
int f_flag;\
f_flag = fread(&(h->n_buckets), sizeof(khint_t), 1, fp);\
f_flag += fread(&(h->size), sizeof(khint_t), 1, fp);\
f_flag += fread(&(h->n_occupied), sizeof(khint_t), 1, fp);\
f_flag += fread(&(h->upper_bound), sizeof(khint_t), 1, fp);\
if (h->n_buckets)\
{\
h->flags = (khint32_t*)kmalloc(__ac_fsize(h->n_buckets) * sizeof(khint32_t));\
f_flag += fread(h->flags, sizeof(khint32_t), __ac_fsize(h->n_buckets), fp);\
h->keys = (khkey_t*)kmalloc(sizeof(khkey_t)*h->n_buckets);\
f_flag += fread(h->keys, sizeof(khkey_t), h->n_buckets, fp);\
h->vals = (khval_t*)kmalloc(sizeof(khval_t)*h->n_buckets);\
f_flag += fread(h->vals, sizeof(khval_t), h->n_buckets, fp);\
}\
}
#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)
#define kh_write(name, h, fp) kh_write_##name(h, fp)
#define kh_load(name, h, fp) kh_load_##name(h, fp)
#endif /* __AC_KHASH_H */
khashl.h
+1
View File
@@ -276,6 +276,7 @@ static kh_inline khint_t __kh_h2b(khint_t hash, khint_t bits) { return hash * 26
KHASHL_INIT(KH_LOCAL, HType, prefix##_m, HType##_m_bucket_t, prefix##_m_hash, prefix##_m_eq) \
SCOPE HType *prefix##_init(void) { return prefix##_m_init(); } \
SCOPE void prefix##_destroy(HType *h) { prefix##_m_destroy(h); } \
SCOPE void prefix##_resize(HType *h, khint_t new_n_buckets) { prefix##_m_resize(h, new_n_buckets); } \
SCOPE khint_t prefix##_get(const HType *h, khkey_t key) { HType##_m_bucket_t t; t.key = key; return prefix##_m_getp(h, &t); } \
SCOPE int prefix##_del(HType *h, khint_t k) { return prefix##_m_del(h, k); } \
SCOPE khint_t prefix##_put(HType *h, khkey_t key, int *absent) { HType##_m_bucket_t t; t.key = key; return prefix##_m_putp(h, &t, absent); }
kmer.cpp
-195
View File
@@ -1,195 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include "kmer.h"
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;
}
void init_Hash_code(Hash_code* code)
{
code->x[0] = 0;
code->x[1] = 0;
}
void init_small_hash_table(small_hash_table* x)
{
x->size = 0;
x->buffer = NULL;
x->length = 0;
}
void clear_small_hash_table(small_hash_table* x)
{
x->length = 0;
}
void resize_small_hash_table(small_hash_table* x, uint64_t size)
{
if(size > x->size)
{
x->size = size;
x->buffer = (k_v*)realloc(x->buffer, x->size*sizeof(k_v));
}
}
void destory_small_hash_table(small_hash_table* x)
{
free(x->buffer);
}
void add_small_hash_table(small_hash_table* x, k_v* element)
{
if(x->length + 1 > x->size)
{
x->size = (x->length + 1) * 2;
x->buffer = (k_v*)realloc(x->buffer, x->size*sizeof(k_v));
}
x->buffer[x->length] = (*element);
x->length++;
}
//x > y, return 1; x < y, return -1, x == y, return 0
int compare_k_mer(k_v* x, k_v* y)
{
if(x->key.x[1] != y->key.x[1])
{
return x->key.x[1] > y->key.x[1] ? 1: -1;
}
else
{
if(x->key.x[0] != y->key.x[0])
{
return x->key.x[0] > y->key.x[0] ? 1: -1;
}
else
{
return 0;
}
}
}
int cmp_k_mer_kv(const void * a, const void * b)
{
int flag = compare_k_mer((k_v*)a, (k_v*)b);
if(flag == 0)
{
if ((*(k_v*)a).value != (*(k_v*)b).value)
{
return (*(k_v*)a).value > (*(k_v*)b).value ? 1: -1;
}
else
{
return 0;
}
}
else
{
return flag;
}
}
void sort_small_hash_table(small_hash_table* x)
{
qsort(x->buffer, x->length, sizeof(k_v), cmp_k_mer_kv);
}
inline long long firstEqual(k_v* arr, long long arrLen, k_v* key)
{
long long L = 0, R = arrLen - 1; //[L, R]
long long mid;
int flag;
while( L <= R)
{
mid = L + (R - L)/2;
flag = compare_k_mer(&(arr[mid]), key);
///arr[mid] >= key
if(flag >= 0)
{
R = mid - 1;
}
else
{
L = mid + 1;
}
}
if(L < arrLen && (flag = compare_k_mer(&(arr[L]), key) == 0))
{
return L;
}
return -1;
}
inline long long lastEqual(k_v* arr, long long arrLen, k_v* key)
{
long long L = 0, R = arrLen - 1; //[L, R]
long long mid;
int flag;
while( L <= R)
{
mid = L + (R - L)/2;
flag = compare_k_mer(&(arr[mid]), key);
///arr[mid] <= key
if(flag <= 0)
{
L = mid + 1;
}
else
{
R = mid - 1;
}
}
if(R >= 0 && ((flag = compare_k_mer(&(arr[R]), key)) == 0))
{
return R;
}
return -1;
}
int query_small_hash_table(small_hash_table* target, k_v* query, long long* l_end, long long* r_end)
{
(*l_end) = -1;
(*r_end) = -1;
long long left_end;
long long right_end;
left_end = firstEqual(target->buffer, target->length, query);
if(left_end != -1)
{
right_end = lastEqual(target->buffer + left_end, target->length - left_end, query) + left_end;
(*l_end) = left_end;
(*r_end) = right_end;
if(right_end == -1)
{
fprintf(stderr, "error\n");
}
return right_end - left_end + 1;
}
return 0;
}
kmer.h
-114
View File
@@ -1,114 +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[2];
} Hash_code;
typedef struct {
Hash_code key; ///k-mer itself
uint64_t value; ///offset
} k_v;
typedef struct {
k_v* buffer;
uint32_t size;
uint32_t length;
} small_hash_table;
void init_small_hash_table(small_hash_table* x);
void clear_small_hash_table(small_hash_table* x);
void resize_small_hash_table(small_hash_table* x, uint64_t size);
void destory_small_hash_table(small_hash_table* x);
void add_small_hash_table(small_hash_table* x, k_v* element);
void sort_small_hash_table(small_hash_table* x);
int compare_k_mer(k_v* x, k_v* y);
int query_small_hash_table(small_hash_table* target, k_v* query, long long* l_end, long long* r_end);
typedef struct
{
char* str;
long long l;
long long i;
long long N_occ;
} HPC_seq;
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 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;
}
inline void Hashcode_to_string(Hash_code* code, char* str, int k)
{
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);
void init_Hash_code(Hash_code* code);
#endif
main.cpp
+10 -6
View File
@@ -4,16 +4,20 @@
#include "Process_Read.h"
#include "Assembly.h"
#include "Levenshtein_distance.h"
#include "htab.h"
int main(int argc, char *argv[])
{
int i, ret;
yak_reset_realtime();
init_opt(&asm_opt);
if (!CommandLine_process(argc, argv, &asm_opt)) return 1;
Correct_Reads(asm_opt.number_of_round);
ret = ha_assemble();
destory_opt(&asm_opt);
return 0;
fprintf(stderr, "[M::%s] Version: %s\n", __func__, HA_VERSION);
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_in_gb());
return ret;
}
sketch.cpp
+110
View File
@@ -0,0 +1,110 @@
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <string.h>
#include "kvec.h"
#include "htab.h"
typedef struct { // a simplified version of kdq
int front, count;
int a[64];
} tiny_queue_t;
static inline void tq_push(tiny_queue_t *q, int x)
{
q->a[((q->count++) + q->front) & 0x3f] = x;
}
static inline int tq_shift(tiny_queue_t *q)
{
int x;
if (q->count == 0) return -1;
x = q->a[q->front++];
q->front &= 0x3f;
--q->count;
return x;
}
/**
* Find symmetric (w,k)-minimizers on a DNA sequence
*
* @param str DNA sequence
* @param len length of $str
* @param w find a minimizer for every $w consecutive k-mers
* @param k k-mer size
* @param rid reference ID; will be copied to the output $p array
* @param is_hpc homopolymer-compressed or not
* @param p minimizers
*/
void ha_sketch(const char *str, int len, int w, int k, uint32_t rid, int is_hpc, ha_mz1_v *p, const void *hf)
{
static const ha_mz1_t dummy = { UINT64_MAX, 0, 0, 0 };
uint64_t shift1 = k - 1, mask = (1ULL<<k) - 1, kmer[4] = {0,0,0,0};
int i, j, l, buf_pos, min_pos, kmer_span = 0;
ha_mz1_t buf[256], min = dummy;
tiny_queue_t tq;
assert(len > 0 && len < 1<<27 && rid < 1<<28 && (w > 0 && w < 256) && (k > 0 && k <= 63));
memset(buf, 0xff, w * 16);
memset(&tq, 0, sizeof(tiny_queue_t));
kv_resize(ha_mz1_t, *p, p->n + len/w);
for (i = l = buf_pos = min_pos = 0; i < len; ++i) {
int c = seq_nt4_table[(uint8_t)str[i]];
ha_mz1_t info = dummy;
if (c < 4) { // not an ambiguous base
int z;
if (is_hpc) {
int skip_len = 1;
if (i + 1 < len && seq_nt4_table[(uint8_t)str[i + 1]] == c) {
for (skip_len = 2; i + skip_len < len; ++skip_len)
if (seq_nt4_table[(uint8_t)str[i + skip_len]] != c)
break;
i += skip_len - 1; // put $i at the end of the current homopolymer run
}
tq_push(&tq, skip_len);
kmer_span += skip_len;
if (tq.count > k) kmer_span -= tq_shift(&tq);
} else kmer_span = l + 1 < k? l + 1 : k;
kmer[0] = (kmer[0] << 1 | (c&1)) & mask; // forward k-mer
kmer[1] = (kmer[1] << 1 | (c>>1)) & mask;
kmer[2] = kmer[2] >> 1 | (uint64_t)(1 - (c&1)) << shift1; // reverse k-mer
kmer[3] = kmer[3] >> 1 | (uint64_t)(1 - (c>>1)) << shift1;
if (kmer[1] == kmer[3]) continue; // skip "symmetric k-mers" as we don't know it strand
z = kmer[1] < kmer[3]? 0 : 1; // strand
++l;
if (l >= k && kmer_span < 256) {
uint64_t y;
y = yak_hash64_64(kmer[z<<1|0]) + yak_hash64_64(kmer[z<<1|1]);
if (hf == 0 || ha_ft_isflt(hf, y) == 0)
info.x = y, info.rid = rid, info.pos = i, info.rev = z, info.span = kmer_span;
}
} else l = 0, tq.count = tq.front = 0, kmer_span = 0;
buf[buf_pos] = info; // need to do this here as appropriate buf_pos and buf[buf_pos] are needed below
if (l == w + k - 1 && min.x != UINT64_MAX) { // special case for the first window - because identical k-mers are not stored yet
for (j = buf_pos + 1; j < w; ++j)
if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
for (j = 0; j < buf_pos; ++j)
if (min.x == buf[j].x && buf[j].pos != min.pos) kv_push(ha_mz1_t, *p, buf[j]);
}
if (info.x <= min.x) { // a new minimum; then write the old min
if (l >= w + k && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
min = info, min_pos = buf_pos;
} else if (buf_pos == min_pos) { // old min has moved outside the window
if (l >= w + k - 1 && min.x != UINT64_MAX) kv_push(ha_mz1_t, *p, min);
for (j = buf_pos + 1, min.x = UINT64_MAX; j < w; ++j) // the two loops are necessary when there are identical k-mers
if (min.x >= buf[j].x) min = buf[j], min_pos = j; // >= is important s.t. min is always the closest k-mer
for (j = 0; j <= buf_pos; ++j)
if (min.x >= buf[j].x) min = buf[j], min_pos = j;
if (l >= w + k - 1 && min.x != UINT64_MAX) { // write identical k-mers
for (j = buf_pos + 1; j < w; ++j) // these two loops make sure the output is sorted
if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
for (j = 0; j <= buf_pos; ++j)
if (min.x == buf[j].x && min.pos != buf[j].pos) kv_push(ha_mz1_t, *p, buf[j]);
}
}
if (++buf_pos == w) buf_pos = 0;
}
if (min.x != UINT64_MAX)
kv_push(ha_mz1_t, *p, min);
}
sys.cpp
+53
View File
@@ -0,0 +1,53 @@
#include <sys/resource.h>
#include <sys/time.h>
#include "htab.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
}
double yak_peakrss_in_gb(void)
{
return yak_peakrss() / 1073741824.0;
}
double yak_cpu_usage(void)
{
return (yak_cputime() + 1e-9) / (yak_realtime() + 1e-9);
}