diff --git a/src/cli_parsing.cpp b/src/cli_parsing.cpp index b6a2707..087d4cc 100644 --- a/src/cli_parsing.cpp +++ b/src/cli_parsing.cpp @@ -72,7 +72,7 @@ static ko_longopt_t long_options[] = { // { (char*)"ava", ko_no_argument, 360 }, { (char*)"r10", ko_no_argument, 361 }, { (char*)"rev-query", ko_no_argument, 362 }, - { (char*)"map-model", ko_required_argument, 363 }, + // { (char*)"map-model", ko_required_argument, 363 }, { (char*)"fine-min", ko_required_argument, 364 }, { (char*)"fine-max", ko_required_argument, 365 }, { (char*)"fine-range", ko_required_argument, 366 }, @@ -390,7 +390,7 @@ CLIParsedArgs parse_args(int argc, char *argv[]) { opt.chain_gap_scale = 1.2f; } else if (c == 362) {ipt.flag |= RI_I_REV_QUERY;}// --rev-query, better name: --idx-dont-store-rev, or negate - else if (c == 363) {opt.model_map_path = o.arg;}// --map-model + // else if (c == 363) {opt.model_map_path = o.arg;}// --map-model else if (c == 364) {ipt.fine_min = atof(o.arg);}// --fine-min else if (c == 365) {ipt.fine_max = atof(o.arg);}// --fine-max else if (c == 366) {ipt.fine_range = atof(o.arg);}// --fine-range diff --git a/src/rawhash_wrapper.cpp b/src/rawhash_wrapper.cpp index da9d563..3ab707a 100644 --- a/src/rawhash_wrapper.cpp +++ b/src/rawhash_wrapper.cpp @@ -122,7 +122,7 @@ std::vector RawHashMapper::map(float* raw_signal, int signal_length) // .su_estimations = NULL, // .su_c_estimations = NULL, // .su_stop = 0, - .map_model = NULL + // .map_model = NULL }; ri_tbuf_t *buf = ri_tbuf_init(); char* read_name = "read123"; diff --git a/src/rmap.cpp b/src/rmap.cpp index ac3c97d..084857f 100644 --- a/src/rmap.cpp +++ b/src/rmap.cpp @@ -25,7 +25,7 @@ double ri_seedsorttime = 0.0; double ri_chainsorttime = 0.0; #endif -static ri_tbuf_t *ri_tbuf_init(void) +ri_tbuf_t *ri_tbuf_init(void) { ri_tbuf_t *b; b = (ri_tbuf_t*)calloc(1, sizeof(ri_tbuf_t)); @@ -33,7 +33,7 @@ static ri_tbuf_t *ri_tbuf_init(void) return b; } -static void ri_tbuf_destroy(ri_tbuf_t *b) +void ri_tbuf_destroy(ri_tbuf_t *b) { if (b == 0) return; ri_km_destroy(b->km); @@ -386,9 +386,9 @@ void ri_map_frag(const ri_idx_t *ri, reg->offset += n_events; } -static void map_worker_for(void *_data, - long i, - int tid) // kt_for() callback +void map_worker_for(void *_data, + long i, + int tid) // kt_for() callback { step_mt *s = (step_mt*)_data; //s->sig and s->n_sig (signals read in this step and num of them) const ri_mapopt_t *opt = s->p->opt; @@ -745,39 +745,15 @@ static void *map_worker_pipeline(void *shared, // fprintf(stderr, "%s %d %d\n", reg0->read_name, reg0->n_maps, reg0->maps[0].mapped); if(reg0->read_name){ - if(reg0->n_maps > 0 && (!p->su_stop || k < p->su_stop)){ - for(uint32_t m = 0; m < reg0->n_maps; ++m){ - if(reg0->maps[m].ref_id < ri->n_seq) - fprintf(stdout, "%s\t%u\t%u\t%u\t%c\t%s\t%u\t%u\t%u\t%u\t%u\t%u\t%s\n", - reg0->read_name, - reg0->maps[m].read_length, - reg0->maps[m].read_start_position, - reg0->maps[m].read_end_position, - reg0->maps[m].rev?'-':'+', - (ri->flag&RI_I_SIG_TARGET)?ri->sig[reg0->maps[m].ref_id].name:ri->seq[reg0->maps[m].ref_id].name, - (ri->flag&RI_I_SIG_TARGET)?ri->sig[reg0->maps[m].ref_id].l_sig:ri->seq[reg0->maps[m].ref_id].len, - reg0->maps[m].fragment_start_position, - reg0->maps[m].fragment_start_position + reg0->maps[m].fragment_length, - reg0->maps[m].read_end_position-reg0->maps[m].read_start_position-1, - reg0->maps[m].fragment_length, - reg0->maps[m].mapq, - reg0->maps[m].tags); - if(reg0->maps[m].tags) {free(reg0->maps[m].tags); reg0->maps[m].tags = NULL;} - } + if(!p->su_stop || k < p->su_stop){ + write_out_mappings_to_stdout(reg0, ri); }else{ - fprintf(stdout, "%s\t%u\t*\t*\t*\t*\t*\t*\t*\t*\t*\t%u\t%s\n", - reg0->read_name, - reg0->maps[0].read_length, - reg0->maps[0].mapq, - reg0->maps[0].tags); - - if(reg0->maps[0].tags) {free(reg0->maps[0].tags); reg0->maps[0].tags = NULL;} + write_out_mappings_to_stdout(reg0, ri, false); } + free_mappings_ri_reg1_t(reg0); + free(reg0); } - - // if(reg0->tags) {free(reg0->tags); reg0->tags = NULL;} - if(reg0->maps){free(reg0->maps); reg0->maps = NULL;} - free(reg0); s->reg[k] = NULL; + s->reg[k] = NULL; // was freed fflush(stdout); } } @@ -799,6 +775,56 @@ static void *map_worker_pipeline(void *shared, return 0; } +void free_mappings_ri_reg1_t(ri_reg1_t *reg0) { + for (int m = 0; m < reg0->n_maps; ++m) { + if (reg0->maps[m].tags) { + free(reg0->maps[m].tags); + reg0->maps[m].tags = NULL; + } + } + if(reg0->maps){free(reg0->maps); reg0->maps = NULL;} +} + +// todo4, todo1: refactor out reg0->maps,read_name into new structure so reg0 can be freed +void write_out_mappings_to_stdout(ri_reg1_t *reg0, const ri_idx_t *ri, bool was_mapped) { + was_mapped = was_mapped && (reg0->n_maps > 0); + if (was_mapped) { + for (int m = 0; m < reg0->n_maps; ++m) { + if (reg0->maps[m].ref_id < ri->n_seq) + fprintf(stdout, "%s\t%u\t%u\t%u\t%c\t%s\t%u\t%u\t%u\t%u\t%u\t%u\t%s\n", + reg0->read_name, + reg0->maps[m].read_length, + reg0->maps[m].read_start_position, + reg0->maps[m].read_end_position, + reg0->maps[m].rev ? '-' : '+', + (ri->flag & RI_I_SIG_TARGET) ? ri->sig[reg0->maps[m].ref_id].name : ri->seq[reg0->maps[m].ref_id].name, + (ri->flag & RI_I_SIG_TARGET) ? ri->sig[reg0->maps[m].ref_id].l_sig : ri->seq[reg0->maps[m].ref_id].len, + reg0->maps[m].fragment_start_position, + reg0->maps[m].fragment_start_position + reg0->maps[m].fragment_length, + reg0->maps[m].read_end_position - reg0->maps[m].read_start_position - 1, + reg0->maps[m].fragment_length, + reg0->maps[m].mapq, + reg0->maps[m].tags); + // if (reg0->maps[m].tags) { + // free(reg0->maps[m].tags); + // reg0->maps[m].tags = NULL; + // } + } + } else { + // maps[0] contains all necessary information + fprintf(stdout, "%s\t%u\t*\t*\t*\t*\t*\t*\t*\t*\t*\t%u\t%s\n", + reg0->read_name, + reg0->maps[0].read_length, + reg0->maps[0].mapq, + reg0->maps[0].tags); + + // if(reg0->maps[0].tags) {free(reg0->maps[0].tags); reg0->maps[0].tags = NULL;} + } + +// if(reg0->maps){free(reg0->maps); reg0->maps = NULL;} + // free(reg0); +} + int ri_map_file(const ri_idx_t *idx, const char *fn, const ri_mapopt_t *opt, diff --git a/src/rmap.h b/src/rmap.h index 5dcce6f..9964776 100644 --- a/src/rmap.h +++ b/src/rmap.h @@ -93,6 +93,27 @@ int ri_map_file(const ri_idx_t *idx, const char *fn, const ri_mapopt_t *opt, int */ int ri_map_file_frag(const ri_idx_t *idx, int n_segs, const char **fn, const ri_mapopt_t *opt, int n_threads, int io_n_threads = 1); + +// functions used in rmap.cpp and rawhash_wrapper.cpp + +// init thread-local buffer +ri_tbuf_t *ri_tbuf_init(void); +// destroy thread-local buffer +void ri_tbuf_destroy(ri_tbuf_t *b); + +// write out relevant mapping results to stdout, then free reg0 +// was_mapped: if false, write out "*" for mapping info even when read was mapped, e.g. when sequence until terminates the pipeline) +void write_out_mappings_to_stdout(ri_reg1_t *reg0, const ri_idx_t *ri, bool was_mapped=true); + +// frees mappings (works even if no mappings were set), which is not freed by "free_most_of_ri_reg1_t" and reg0 itself +void free_mappings_ri_reg1_t(ri_reg1_t* reg0); +// finally need to call free(reg0) + +// map a single read to a reference, called in parallel +void map_worker_for(void *_data, + long i, + int tid); + #ifdef __cplusplus } #endif diff --git a/src/roptions.h b/src/roptions.h index a6b6f36..16dc780 100644 --- a/src/roptions.h +++ b/src/roptions.h @@ -11,6 +11,7 @@ #define RI_I_SYNCMER 0x8 #define RI_I_STORE_SIG 0x10 #define RI_I_SIG_TARGET 0x20 +#define RI_I_REV_QUERY 0x40 // whether the index should NOT store the reverse strand, this assumes that we query the index with the revcomp of the read as well #define RI_I_NO_REV_TARGET 0x40 #define RI_I_OUT_QUANTIZE 0x80 #define RI_I_NO_EVENT_DETECTION 0x100 @@ -31,6 +32,7 @@ #define RI_M_LOG_NUM_ANCHORS 0x1000 //Overlapping related #define RI_M_ALL_CHAINS 0x2000 +#define RI_M_DONT_QUERY_REVCOMP 0x4000 // whether to transform the read to its revcomp and query index (defaults to RI_I_REV_QUERY) //Characterization related #define RI_M_OUT_ALL_CHAINS 0x4000 diff --git a/src/rseed.c b/src/rseed.c index 5bd2a6f..9af347b 100644 --- a/src/rseed.c +++ b/src/rseed.c @@ -7,8 +7,6 @@ void ri_seed_select(int32_t n, ri_seed_t *a, uint32_t len, int max_occ, int max_max_occ, int dist) { // for high-occ minimizers, choose up to max_high_occ in each high-occ streak - extern void ks_heapdown_uint64_t(size_t i, size_t n, uint64_t*); - extern void ks_heapmake_uint64_t(size_t n, uint64_t*); int32_t i, last0, m; uint64_t b[MAX_MAX_HIGH_OCC]; // this is to avoid a heap allocation diff --git a/src/rutils.h b/src/rutils.h index 6c0fbf9..d9e8e42 100644 --- a/src/rutils.h +++ b/src/rutils.h @@ -38,6 +38,10 @@ void liftrlimit(void); void load_pore(const char* fpore, const short k, const short lev_col, ri_pore_t* pore); +// Explicit definition for rseed.c +extern void ks_heapdown_uint64_t(size_t i, size_t n, uint64_t*); +extern void ks_heapmake_uint64_t(size_t n, uint64_t*); + #ifdef __cplusplus } #endif