From 2e5ba6057059b6b602556df168bb47dc195f295b Mon Sep 17 00:00:00 2001 From: vasudeva8 Date: Sun, 6 Sep 2026 15:39:29 +0100 Subject: [PATCH] filtering reads based on depth, review updates using last read for a position Signed-off-by: vasudeva8 --- Makefile | 7 +- hts.c | 736 +++++++++++++++--------- htslib/hts.h | 8 +- sam.c | 75 ++- sam_cache.c | 1069 +++++++++++++++++++++++++++++++++++ sam_cache.h | 132 +++++ test/cache.exp.1.T1T3T2.sam | 29 + test/cache.exp.1.sam | 31 + test/cache.exp.2.sam | 33 ++ test/cache.exp.T1.sam | 20 + test/cache.exp.w5.sam | 25 + test/cache.sam | 34 ++ test/test.pl | 53 ++ 13 files changed, 1974 insertions(+), 278 deletions(-) create mode 100644 sam_cache.c create mode 100644 sam_cache.h create mode 100644 test/cache.exp.1.T1T3T2.sam create mode 100644 test/cache.exp.1.sam create mode 100644 test/cache.exp.2.sam create mode 100644 test/cache.exp.T1.sam create mode 100644 test/cache.exp.w5.sam create mode 100644 test/cache.sam diff --git a/Makefile b/Makefile index 16ee1879f..d6c3502a6 100644 --- a/Makefile +++ b/Makefile @@ -228,6 +228,7 @@ LIBHTS_OBJS = \ regidx.o \ region.o \ sam.o \ + sam_cache.o \ sam_mods.o \ simd.o \ synced_bcf_reader.o \ @@ -292,6 +293,7 @@ hfile_internal_h = hfile_internal.h $(htslib_hts_defs_h) $(htslib_hfile_h) $(tex hts_internal_h = hts_internal.h $(htslib_hts_h) $(textutils_internal_h) hts_time_funcs_h = hts_time_funcs.h sam_internal_h = sam_internal.h $(htslib_sam_h) +sam_cache_h = sam_cache.h $(htslib_khash_h) $(htslib_sam_h) textutils_internal_h = textutils_internal.h $(htslib_kstring_h) thread_pool_internal_h = thread_pool_internal.h $(htslib_thread_pool_h) @@ -511,11 +513,12 @@ hfile.o hfile.pico: hfile.c config.h $(htslib_hfile_h) $(hfile_internal_h) $(hts hfile_gcs.o hfile_gcs.pico: hfile_gcs.c config.h $(htslib_hts_h) $(htslib_kstring_h) $(hfile_internal_h) hfile_libcurl.o hfile_libcurl.pico: hfile_libcurl.c config.h $(hfile_internal_h) $(htslib_hts_h) $(htslib_hts_alloc_h) $(htslib_kstring_h) $(htslib_khash_h) hfile_s3.o hfile_s3.pico: hfile_s3.c config.h $(hfile_internal_h) $(htslib_hts_h) $(htslib_hts_alloc_h) $(htslib_kstring_h) $(hts_time_funcs_h) -hts.o hts.pico: hts.c config.h os/lzma_stub.h $(htslib_hts_h) $(htslib_bgzf_h) $(cram_h) $(htslib_hfile_h) $(htslib_hts_endian_h) version.h config_vars.h $(hts_internal_h) $(hfile_internal_h) $(sam_internal_h) $(htslib_hts_alloc_h) $(htslib_hts_expr_h) $(htslib_hts_os_h) $(htslib_khash_h) $(htslib_kseq_h) $(htslib_ksort_h) $(htslib_tbx_h) $(htscodecs_htscodecs_h) +hts.o hts.pico: hts.c config.h os/lzma_stub.h $(htslib_hts_h) $(htslib_bgzf_h) $(cram_h) $(htslib_hfile_h) $(htslib_hts_endian_h) version.h config_vars.h $(hts_internal_h) $(hfile_internal_h) $(sam_internal_h) $(htslib_hts_alloc_h) $(htslib_hts_expr_h) $(htslib_hts_os_h) $(htslib_khash_h) $(htslib_kseq_h) $(htslib_ksort_h) $(htslib_tbx_h) $(htscodecs_htscodecs_h) $(sam_cache_h) hts_expr.o hts_expr.pico: hts_expr.c config.h $(htslib_hts_expr_h) $(htslib_hts_alloc_h) $(htslib_hts_log_h) $(textutils_internal_h) hts_os.o hts_os.pico: hts_os.c config.h $(htslib_hts_defs_h) os/rand.c vcf.o vcf.pico: vcf.c config.h $(fuzz_settings_h) $(htslib_vcf_h) $(htslib_bgzf_h) $(htslib_tbx_h) $(htslib_hfile_h) $(hts_internal_h) $(htslib_hts_alloc_h) $(htslib_hts_endian_h) $(htslib_khash_str2int_h) $(htslib_kstring_h) $(htslib_sam_h) $(htslib_khash_h) $(htslib_kseq_h) $(bgzf_internal_h) -sam.o sam.pico: sam.c config.h $(fuzz_settings_h) $(htslib_hts_defs_h) $(htslib_sam_h) $(htslib_bgzf_h) $(cram_h) $(hts_internal_h) $(sam_internal_h) $(htslib_hfile_h) $(htslib_hts_alloc_h) $(htslib_hts_endian_h) $(htslib_hts_expr_h) $(header_h) $(bgzf_internal_h) $(htslib_khash_h) $(htslib_kseq_h) $(htslib_kstring_h) +sam.o sam.pico: sam.c config.h $(fuzz_settings_h) $(htslib_hts_defs_h) $(htslib_sam_h) $(htslib_bgzf_h) $(cram_h) $(hts_internal_h) $(sam_internal_h) $(htslib_hfile_h) $(htslib_hts_alloc_h) $(htslib_hts_endian_h) $(htslib_hts_expr_h) $(header_h) $(bgzf_internal_h) $(htslib_khash_h) $(htslib_kseq_h) $(htslib_kstring_h) $(sam_cache_h) +sam_cache.o sam_cache.pico: sam_cache.c config.h $(sam_cache_h) $(htslib_hts_alloc_h) sam_mods.o sam_mods.pico: sam_mods.c config.h $(htslib_sam_h) $(textutils_internal_h) simd.o simd.pico: simd.c config.h $(htslib_sam_h) $(sam_internal_h) tbx.o tbx.pico: tbx.c config.h $(htslib_tbx_h) $(htslib_bgzf_h) $(htslib_hts_alloc_h) $(htslib_hts_endian_h) $(hts_internal_h) $(bgzf_internal_h) $(htslib_khash_h) diff --git a/hts.c b/hts.c index 8b3919591..7744a3669 100644 --- a/hts.c +++ b/hts.c @@ -62,6 +62,7 @@ DEALINGS IN THE SOFTWARE. */ #include "htslib/hts_alloc.h" #include "htslib/hts_expr.h" #include "htslib/hts_os.h" // drand48 +#include "sam_cache.h" #include "htslib/khash.h" #include "htslib/kseq.h" @@ -1222,6 +1223,13 @@ int hts_opt_add(hts_opt **opts, const char *c_arg) { strcmp(o->arg, "REMOVE_UR") == 0) o->opt = CRAM_OPT_RM_UR, o->val.i = atoi(val); + else if (strcmp(o->arg, "hts_flt_depth") == 0 || + strcmp(o->arg, "HTS_FLT_DEPTH") == 0) //todo make more SAM specific? + o->opt = HTS_OPT_CACHE_FLT_DEPTH, o->val.i = atoi(val); + else if (strcmp(o->arg, "hts_flt_wndsize") == 0 || + strcmp(o->arg, "HTS_FLT_WNDSIZE") == 0) //todo make more SAM specific? + o->opt = HTS_OPT_CACHE_FLT_SIZE, o->val.i = atoi(val); + else { hts_log_error("Unknown option '%s'", o->arg); free(o->arg); @@ -1705,6 +1713,7 @@ int hts_close(htsFile *fp) } save = errno; + destroy_readcache(fp); sam_hdr_destroy(fp->bam_header); hts_idx_destroy(fp->idx); hts_filter_free(fp->filter); @@ -1910,6 +1919,29 @@ int hts_set_opt(htsFile *fp, enum hts_fmt_option opt, ...) { break; } + case HTS_OPT_CACHE_FLT_DEPTH: { + va_start(args, opt); + int dpth = va_arg(args, int); + va_end(args); + if (dpth > 0) { + if(setup_readcache(fp, 0, dpth)) { + hts_log_warning("Failed to setup hts cache depth"); + } + } + return 0; + } + case HTS_OPT_CACHE_FLT_SIZE: { + va_start(args, opt); + int wndsz = va_arg(args, int); + va_end(args); + if (wndsz > 0) { + if(setup_readcache(fp, wndsz, 0)) { + hts_log_warning("Failed to setup hts cache window size"); + } + } + return 0; + } + default: break; } @@ -4276,11 +4308,20 @@ hts_itr_t *hts_itr_regions(const hts_idx_t *idx, hts_reglist_t *reglist, int cou return itr; } -int hts_itr_next(BGZF *fp, hts_itr_t *iter, void *r, void *data) +int hts_itr_next(BGZF *fp, hts_itr_t *iter, void *s, void *data) { int ret, tid; + cs sts = NOTREADY; hts_pos_t beg, end; + void *c = get_sam_readcache(iter, data); + void *e = NULL; + void *r = s; + if (iter == NULL || iter->finished) return -1; + //if sam caching in use, get from cache + if (c && getfrom_readcache_iter(c, s, &iter->curr_tid, &iter->curr_beg, &iter->curr_end)) { + return 0; + } if (iter->read_rest) { if (iter->curr_off) { // seek to the start if (bgzf_seek(fp, iter->curr_off, SEEK_SET) < 0) { @@ -4291,60 +4332,133 @@ int hts_itr_next(BGZF *fp, hts_itr_t *iter, void *r, void *data) } iter->curr_off = 0; // only seek once } - ret = iter->readrec(fp, data, r, &tid, &beg, &end); - if (ret < 0) iter->finished = 1; - iter->curr_tid = tid; - iter->curr_beg = beg; - iter->curr_end = end; - return ret; } // A NULL iter->off should always be accompanied by iter->finished. - assert(iter->off != NULL); + assert(iter->off != NULL || iter->read_rest); for (;;) { - if (iter->curr_off == 0 || iter->curr_off >= iter->off[iter->i].v) { // then jump to the next chunk - if (iter->i == iter->n_off - 1) { ret = -1; break; } // no more chunks - if (iter->i < 0 || iter->off[iter->i].v != iter->off[iter->i+1].u) { // not adjacent chunks; then seek - if (bgzf_seek(fp, iter->off[iter->i+1].u, SEEK_SET) < 0) { - hts_log_error("Failed to seek to offset %"PRIu64"%s%s", - iter->off[iter->i+1].u, - errno ? ": " : "", strerror(errno)); - return -2; - } - iter->curr_off = bgzf_tell(fp); - } - ++iter->i; + if (c) { //sam cache in use + if (!(e = get_cache_iter(data))) //get cached storage + return -3; + r = get_readbuffer_iter(e); //get bam pointer from retrived storage + sts = 0; } - if ((ret = iter->readrec(fp, data, r, &tid, &beg, &end)) >= 0) { - iter->curr_off = bgzf_tell(fp); - if (tid != iter->tid || beg >= iter->end) { // no need to proceed - ret = -1; break; - } else if (end > iter->beg && iter->end > beg) { + if (iter->read_rest) { + ret = iter->readrec(fp, data, r, &tid, &beg, &end); + if (ret < 0) { + if (c) { //end of file or iterator, mark in cache as well + notify_end_iter(c, e); + } else { + iter->curr_tid = tid; + iter->curr_beg = beg; + iter->curr_end = end; + } + break; + } else if (!c) { iter->curr_tid = tid; iter->curr_beg = beg; iter->curr_end = end; return ret; + } else { //cache in use, add read to cache + if (addto_readcache_iter(c, e, &sts)) { + return -3; + } + if (sts >= 2) //ready/wnd full/end + break; + continue; + } + } else { + if (iter->curr_off == 0 || iter->curr_off >= iter->off[iter->i].v) { // then jump to the next chunk + if (iter->i == iter->n_off - 1) { + ret = -1; + if (c) { //end of file or iterator, mark in cache as well + notify_end_iter(c, e); + } + break; + } // no more chunks + if (iter->i < 0 || iter->off[iter->i].v != iter->off[iter->i+1].u) { // not adjacent chunks; then seek + if (bgzf_seek(fp, iter->off[iter->i+1].u, SEEK_SET) < 0) { + hts_log_error("Failed to seek to offset %"PRIu64"%s%s", + iter->off[iter->i+1].u, + errno ? ": " : "", strerror(errno)); + return -2; + } + iter->curr_off = bgzf_tell(fp); + } + ++iter->i; + } + if ((ret = iter->readrec(fp, data, r, &tid, &beg, &end)) >= 0) { + iter->curr_off = bgzf_tell(fp); + if (tid != iter->tid || beg >= iter->end) { // no need to proceed + ret = -1; + if (c) { //end of file or iterator, mark in cache as well + notify_end_iter(c, e); + } + break; + } else if (end > iter->beg && iter->end > beg) { + if (c) { //cache in use, add read to cache + if (addto_readcache_iter(c, e, &sts)) { + return -3; + } + if (sts >= 2) //ready/wnd full/end + break; + continue; + } + iter->curr_tid = tid; + iter->curr_beg = beg; + iter->curr_end = end; + return ret; + } else { //non interested data? + if (c) + ret_cache(c,e); + } + } else { + if (c && ret == -1) { //eof + notify_end_iter(c, e); + } + break; // end of file or error } - } else break; // end of file or error + } + } + if (c && ret >= -1) { //sam cache in use + //process cached reads + if (process_readcache_iter(c) < 0) + return -3; + //get reads to be consumed + if (!getfrom_readcache_iter(c, s, &iter->curr_tid, &iter->curr_beg, &iter->curr_end)) { + iter->finished = 1; + reset_readcache_iter(c); + return -1; + } else if (ret < 0) + ret = 0; + return ret; + } else { + iter->finished = 1; } - iter->finished = 1; return ret; } -int hts_itr_multi_next(htsFile *fd, hts_itr_t *iter, void *r) +int hts_itr_multi_next(htsFile *fd, hts_itr_t *iter, void *s) { void *fp; - int ret, tid, i, cr, ci; + int ret, tid, i, cr, ci, used = 0; + cs sts = NOTREADY; hts_pos_t beg, end; hts_reglist_t *found_reg; + void *c = (void*)fd->c; + void *e = NULL; + void *r = s; if (iter == NULL || iter->finished) return -1; + //if sam caching in use, get from cache + if (c && getfrom_readcache_iter(c, s, &iter->curr_tid, &iter->curr_beg, &iter->curr_end)) { + return 0; + } if (iter->is_cram) { fp = fd->fp.cram; } else { fp = fd->fp.bgzf; } - if (iter->read_rest) { if (iter->curr_off) { // seek to the start if (iter->seek(fp, iter->curr_off, SEEK_SET) < 0) { @@ -4354,267 +4468,359 @@ int hts_itr_multi_next(htsFile *fd, hts_itr_t *iter, void *r) iter->curr_off = 0; // only seek once } - ret = iter->readrec(fp, fd, r, &tid, &beg, &end); - if (ret < 0) - iter->finished = 1; - - iter->curr_tid = tid; - iter->curr_beg = beg; - iter->curr_end = end; - - return ret; - } - // A NULL iter->off should always be accompanied by iter->finished. - assert(iter->off != NULL || iter->nocoor != 0); - - int next_range = 0; - for (;;) { - // Note that due to the way bam indexing works, iter->off may contain - // file chunks that are not actually needed as they contain data - // beyond the end of the requested region. These are filtered out - // by comparing the tid and index into hts_reglist_t::intervals - // (packed for reasons of convenience into iter->off[iter->i].max) - // associated with the file region with iter->curr_tid and - // iter->curr_intv. - - if (next_range - || iter->curr_off == 0 - || iter->i >= iter->n_off - || iter->curr_off >= iter->off[iter->i].v - || (iter->off[iter->i].max >> 32 == iter->curr_tid - && (iter->off[iter->i].max & 0xffffffff) < iter->curr_intv)) { - - // Jump to the next chunk. It may be necessary to skip more - // than one as the iter->off list can include overlapping entries. - do { - iter->i++; - } while (iter->i < iter->n_off - && (iter->curr_off >= iter->off[iter->i].v - || (iter->off[iter->i].max >> 32 == iter->curr_tid - && (iter->off[iter->i].max & 0xffffffff) < iter->curr_intv))); - - if (iter->is_cram && iter->i < iter->n_off) { - // Ensure iter->curr_reg is correct. - // - // We need this for CRAM as we shortcut some of the later - // logic by getting an end-of-range and continuing to the - // next offset. - // - // We cannot do this for BAM (and fortunately do not need to - // either) because in BAM world a query to genomic positions - // GX and GY leading to a seek offsets PX and PY may have - // GX > GY and PX < PY. (This is due to the R-tree and falling - // between intervals, bumping up to a higher bin.) - // CRAM strictly follows PX >= PY if GX >= GY, so this logic - // works. - int want_tid = iter->off[iter->i].max >> 32; - if (!(iter->curr_reg < iter->n_reg && - iter->reg_list[iter->curr_reg].tid == want_tid)) { - int j; - for (j = 0; j < iter->n_reg; j++) - if (iter->reg_list[j].tid == want_tid) - break; - if (j == iter->n_reg) - return -1; - iter->curr_reg = j; - iter->curr_tid = iter->reg_list[iter->curr_reg].tid; - }; - iter->curr_intv = iter->off[iter->i].max & 0xffffffff; + for(;;) { + if (c) { //sam cache in use + if (!(e = get_cache_iter(fd))) //get cached storage + return -3; + r = get_readbuffer_iter(e); //get bam pointer from retrived storage + sts = 0; } + ret = iter->readrec(fp, fd, r, &tid, &beg, &end); + if (ret < 0) { + if (c) { //end of file or iterator, mark in cache as well + notify_end_iter(c, e); + } + break; + } else if (!c) { + iter->curr_tid = tid; + iter->curr_beg = beg; + iter->curr_end = end; + return ret; + } else { //sam cache in use, add to read cache + if (addto_readcache_iter(c, e, &sts)) { + return -3; + } + if (sts >= 2) //ready/wnd full/end + break; + continue; + } + } + } else { + // A NULL iter->off should always be accompanied by iter->finished. + assert(iter->off != NULL || iter->nocoor != 0); + + int next_range = 0; + for (;;) { + // Note that due to the way bam indexing works, iter->off may contain + // file chunks that are not actually needed as they contain data + // beyond the end of the requested region. These are filtered out + // by comparing the tid and index into hts_reglist_t::intervals + // (packed for reasons of convenience into iter->off[iter->i].max) + // associated with the file region with iter->curr_tid and + // iter->curr_intv. + + if (c) { //sam cache in use + if (!(e = get_cache_iter(fd))) //get cached storage + return -3; + r = get_readbuffer_iter(e); //get bam pointer from retrived storage + used = 0; + sts = 0; + } + if (next_range + || iter->curr_off == 0 + || iter->i >= iter->n_off + || iter->curr_off >= iter->off[iter->i].v + || (iter->off[iter->i].max >> 32 == iter->curr_tid + && (iter->off[iter->i].max & 0xffffffff) < iter->curr_intv)) { + + // Jump to the next chunk. It may be necessary to skip more + // than one as the iter->off list can include overlapping entries. + do { + iter->i++; + } while (iter->i < iter->n_off + && (iter->curr_off >= iter->off[iter->i].v + || (iter->off[iter->i].max >> 32 == iter->curr_tid + && (iter->off[iter->i].max & 0xffffffff) < iter->curr_intv))); + + if (iter->is_cram && iter->i < iter->n_off) { + // Ensure iter->curr_reg is correct. + // + // We need this for CRAM as we shortcut some of the later + // logic by getting an end-of-range and continuing to the + // next offset. + // + // We cannot do this for BAM (and fortunately do not need to + // either) because in BAM world a query to genomic positions + // GX and GY leading to a seek offsets PX and PY may have + // GX > GY and PX < PY. (This is due to the R-tree and falling + // between intervals, bumping up to a higher bin.) + // CRAM strictly follows PX >= PY if GX >= GY, so this logic + // works. + int want_tid = iter->off[iter->i].max >> 32; + if (!(iter->curr_reg < iter->n_reg && + iter->reg_list[iter->curr_reg].tid == want_tid)) { + int j; + for (j = 0; j < iter->n_reg; j++) + if (iter->reg_list[j].tid == want_tid) + break; + if (j == iter->n_reg) + return -1; + iter->curr_reg = j; + iter->curr_tid = iter->reg_list[iter->curr_reg].tid; + }; + iter->curr_intv = iter->off[iter->i].max & 0xffffffff; + } - if (iter->i >= iter->n_off) { // no more chunks, except NOCOORs - if (iter->nocoor) { - next_range = 0; - if (iter->seek(fp, iter->nocoor_off, SEEK_SET) < 0) { - hts_log_error("Seek at offset %" PRIu64 " failed.", iter->nocoor_off); - return -2; - } - if (iter->is_cram) { - cram_range r = { HTS_IDX_NOCOOR }; - cram_set_option(fp, CRAM_OPT_RANGE_NOSEEK, &r); - } - - // The first slice covering the unmapped reads might - // contain a few mapped reads, so scroll - // forward until finding the first unmapped read. - do { - ret = iter->readrec(fp, fd, r, &tid, &beg, &end); - } while (tid >= 0 && ret >=0); - - if (ret < 0) - iter->finished = 1; - else - iter->read_rest = 1; + if (iter->i >= iter->n_off) { // no more chunks, except NOCOORs + if (iter->nocoor) { + next_range = 0; + if (iter->seek(fp, iter->nocoor_off, SEEK_SET) < 0) { + hts_log_error("Seek at offset %" PRIu64 " failed.", iter->nocoor_off); + return -2; + } + if (iter->is_cram) { + cram_range r = { HTS_IDX_NOCOOR }; + cram_set_option(fp, CRAM_OPT_RANGE_NOSEEK, &r); + } - iter->curr_off = 0; // don't seek any more - iter->curr_tid = tid; - iter->curr_beg = beg; - iter->curr_end = end; + // The first slice covering the unmapped reads might + // contain a few mapped reads, so scroll + // forward until finding the first unmapped read. + do { + ret = iter->readrec(fp, fd, r, &tid, &beg, &end); + } while (tid >= 0 && ret >=0); - return ret; - } else { - ret = -1; break; - } - } else if (iter->i < iter->n_off) { - // New chunk may overlap the last one, so ensure we - // only seek forwards. - if (iter->curr_off < iter->off[iter->i].u || next_range) { - iter->curr_off = iter->off[iter->i].u; - - // CRAM has the capability of setting an end location. - // This means multi-threaded decodes can stop once they - // reach that point, rather than pointlessly decoding - // more slices than we'll be using. - // - // We have to be careful here. Whenever we set the cram - // range we need a corresponding seek in order to ensure - // we can safely decode at that offset. We use next_range - // var to ensure this is always true; this is set on - // end-of-range condition. It's never modified for BAM. - if (iter->is_cram) { - // Next offset.[uv] tuple, but it's already been - // included in our cram range, so don't seek and don't - // reset range so we can efficiently multi-thread. - if (next_range || iter->curr_off >= iter->end) { - if (iter->seek(fp, iter->curr_off, SEEK_SET) < 0) { - hts_log_error("Seek at offset %" PRIu64 - " failed.", iter->curr_off); - return -2; + if (ret < 0) { + if (c) { //end of file or iterator, mark in cache as well + notify_end_iter(c, e); } - - // Find the genomic range matching this interval. - int j; - hts_reglist_t *rl = &iter->reg_list[iter->curr_reg]; - cram_range r = { - rl->tid, - rl->intervals[iter->curr_intv].beg, - rl->intervals[iter->curr_intv].end - }; - - // Expand it up to cover neighbouring intervals. - // Note we can only have a single chromosome in a - // range, so if we detect our blocks span chromosomes - // or we have a multi-ref mode slice, we just use - // HTS_IDX_START refid instead. This doesn't actually - // seek (due to CRAM_OPT_RANGE_NOSEEK) and is simply - // and indicator of decoding with no end limit. - // - // That isn't as efficient as it could be, but it's - // no poorer than before and it works. - int tid = r.refid; - int64_t end = r.end; - int64_t v = iter->off[iter->i].v; - j = iter->i+1; - while (j < iter->n_off) { - if (iter->off[j].u > v) + iter->curr_off = 0; // don't seek any more + iter->curr_tid = tid; + iter->curr_beg = beg; + iter->curr_end = end; + break; + } + else { + iter->read_rest = 1; + iter->curr_off = 0; // don't seek any more + if (!c) { + iter->curr_tid = tid; + iter->curr_beg = beg; + iter->curr_end = end; + return ret; + } else { //cache in use, add read to cache + if (addto_readcache_iter(c, e, &sts)) { + return -3; + } + if (sts >= 2) //ready/wnd full/end break; + continue; + } + } + } else { + ret = -1; + if (c) { //end of file or iterator, mark in cache as well + notify_end_iter(c, e); + } + break; + } + } else if (iter->i < iter->n_off) { + // New chunk may overlap the last one, so ensure we + // only seek forwards. + if (iter->curr_off < iter->off[iter->i].u || next_range) { + iter->curr_off = iter->off[iter->i].u; + + // CRAM has the capability of setting an end location. + // This means multi-threaded decodes can stop once they + // reach that point, rather than pointlessly decoding + // more slices than we'll be using. + // + // We have to be careful here. Whenever we set the cram + // range we need a corresponding seek in order to ensure + // we can safely decode at that offset. We use next_range + // var to ensure this is always true; this is set on + // end-of-range condition. It's never modified for BAM. + if (iter->is_cram) { + // Next offset.[uv] tuple, but it's already been + // included in our cram range, so don't seek and don't + // reset range so we can efficiently multi-thread. + if (next_range || iter->curr_off >= iter->end) { + if (iter->seek(fp, iter->curr_off, SEEK_SET) < 0) { + hts_log_error("Seek at offset %" PRIu64 + " failed.", iter->curr_off); + return -2; + } - uint64_t max = iter->off[j].max; - if ((max>>32) != tid) { - tid = HTS_IDX_START; // => no range limit - } else { - if (end < rl->intervals[max & 0xffffffff].end) - end = rl->intervals[max & 0xffffffff].end; + // Find the genomic range matching this interval. + int j; + hts_reglist_t *rl = &iter->reg_list[iter->curr_reg]; + cram_range r = { + rl->tid, + rl->intervals[iter->curr_intv].beg, + rl->intervals[iter->curr_intv].end + }; + + // Expand it up to cover neighbouring intervals. + // Note we can only have a single chromosome in a + // range, so if we detect our blocks span chromosomes + // or we have a multi-ref mode slice, we just use + // HTS_IDX_START refid instead. This doesn't actually + // seek (due to CRAM_OPT_RANGE_NOSEEK) and is simply + // and indicator of decoding with no end limit. + // + // That isn't as efficient as it could be, but it's + // no poorer than before and it works. + int tid = r.refid; + int64_t end = r.end; + int64_t v = iter->off[iter->i].v; + j = iter->i+1; + while (j < iter->n_off) { + if (iter->off[j].u > v) + break; + + uint64_t max = iter->off[j].max; + if ((max>>32) != tid) { + tid = HTS_IDX_START; // => no range limit + } else { + if (end < rl->intervals[max & 0xffffffff].end) + end = rl->intervals[max & 0xffffffff].end; + } + if (v < iter->off[j].v) + v = iter->off[j].v; + j++; } - if (v < iter->off[j].v) - v = iter->off[j].v; - j++; + r.refid = tid; + r.end = end; + + // Remember maximum 'v' here so we don't do + // unnecessary subsequent seeks for the next + // regions. We can't change curr_off, but + // beg/end are used only by single region iterator so + // we cache it there to avoid changing the struct. + iter->end = v; + + cram_set_option(fp, CRAM_OPT_RANGE_NOSEEK, &r); + next_range = 0; + } + } else { // Not CRAM + if (iter->seek(fp, iter->curr_off, SEEK_SET) < 0) { + hts_log_error("Seek at offset %" PRIu64 " failed.", + iter->curr_off); + return -2; } - r.refid = tid; - r.end = end; - - // Remember maximum 'v' here so we don't do - // unnecessary subsequent seeks for the next - // regions. We can't change curr_off, but - // beg/end are used only by single region iterator so - // we cache it there to avoid changing the struct. - iter->end = v; - - cram_set_option(fp, CRAM_OPT_RANGE_NOSEEK, &r); - next_range = 0; - } - } else { // Not CRAM - if (iter->seek(fp, iter->curr_off, SEEK_SET) < 0) { - hts_log_error("Seek at offset %" PRIu64 " failed.", - iter->curr_off); - return -2; } } } } - } - ret = iter->readrec(fp, fd, r, &tid, &beg, &end); - if (ret < 0) { - if (iter->is_cram && cram_eof(fp)) { - // Skip to end of range - // - // We should never be adjusting curr_off manually unless - // we also can guarantee we'll be doing a seek after to - // a new location. Otherwise we'll be reading wrong offset - // for the next container. - // - // We ensure this by adjusting our CRAM_OPT_RANGE - // accordingly above, but to double check we also - // set the skipped_block flag to enforce a seek also. - iter->curr_off = iter->off[iter->i].v; - next_range = 1; - - // Next region - if (++iter->curr_intv >= iter->reg_list[iter->curr_reg].count){ - if (++iter->curr_reg >= iter->n_reg) - break; - iter->curr_intv = 0; - iter->curr_tid = iter->reg_list[iter->curr_reg].tid; + ret = iter->readrec(fp, fd, r, &tid, &beg, &end); + if (ret < 0) { + if (iter->is_cram && cram_eof(fp)) { + // Skip to end of range + // + // We should never be adjusting curr_off manually unless + // we also can guarantee we'll be doing a seek after to + // a new location. Otherwise we'll be reading wrong offset + // for the next container. + // + // We ensure this by adjusting our CRAM_OPT_RANGE + // accordingly above, but to double check we also + // set the skipped_block flag to enforce a seek also. + iter->curr_off = iter->off[iter->i].v; + next_range = 1; + + // Next region + if (++iter->curr_intv >= iter->reg_list[iter->curr_reg].count){ + if (++iter->curr_reg >= iter->n_reg) { + if (c) { + notify_end_iter(c, e); + } + break; + } + iter->curr_intv = 0; + iter->curr_tid = iter->reg_list[iter->curr_reg].tid; + } + if (c) { //end of file or iterator, mark in cache as well + notify_end_iter(c, e); + } + continue; + } else { + if (c) { //end of file or iterator, mark in cache as well + notify_end_iter(c, e); + } + break; } - continue; - } else { - break; } - } - iter->curr_off = iter->tell(fp); + iter->curr_off = iter->tell(fp); - if (tid != iter->curr_tid) { - hts_reglist_t key; - key.tid = tid; + if (tid != iter->curr_tid) { + hts_reglist_t key; + key.tid = tid; - found_reg = (hts_reglist_t *)bsearch(&key, iter->reg_list, - iter->n_reg, - sizeof(hts_reglist_t), - compare_regions); - if (!found_reg) - continue; + found_reg = (hts_reglist_t *)bsearch(&key, iter->reg_list, + iter->n_reg, + sizeof(hts_reglist_t), + compare_regions); + if (!found_reg) { + if (c) { + ret_cache(c,e); + } + continue; + } - iter->curr_reg = (found_reg - iter->reg_list); - iter->curr_tid = tid; - iter->curr_intv = 0; - } + iter->curr_reg = (found_reg - iter->reg_list); + iter->curr_tid = tid; + iter->curr_intv = 0; + } - cr = iter->curr_reg; - ci = iter->curr_intv; + cr = iter->curr_reg; + ci = iter->curr_intv; + + for (i = ci; i < iter->reg_list[cr].count; i++) { + if (end > iter->reg_list[cr].intervals[i].beg && + iter->reg_list[cr].intervals[i].end > beg) { + if (c) { //add read to cache + used = 1; + sts = 0; + if (addto_readcache_iter(c, e, &sts)) { + return -3; + } + break; + } else { + iter->curr_beg = beg; + iter->curr_end = end; + iter->curr_intv = i; + return ret; + } + } - for (i = ci; i < iter->reg_list[cr].count; i++) { - if (end > iter->reg_list[cr].intervals[i].beg && - iter->reg_list[cr].intervals[i].end > beg) { - iter->curr_beg = beg; - iter->curr_end = end; - iter->curr_intv = i; + // Check if the read starts beyond intervals[i].end + // If so, the interval is finished so move on to the next. + if (beg > iter->reg_list[cr].intervals[i].end) + iter->curr_intv = i + 1; - return ret; + // No need to keep searching if the read ends before intervals[i].beg + if (end < iter->reg_list[cr].intervals[i].beg) { + break; + } + } + if (c) { //cache in use + if (used) { + if(sts >= WNDFULL) //ready/wnd full/end + break; + } else { //unused, return to cache + ret_cache(c,e); + } } - - // Check if the read starts beyond intervals[i].end - // If so, the interval is finished so move on to the next. - if (beg > iter->reg_list[cr].intervals[i].end) - iter->curr_intv = i + 1; - - // No need to keep searching if the read ends before intervals[i].beg - if (end < iter->reg_list[cr].intervals[i].beg) - break; } } - iter->finished = 1; + if (c && ret >= -1) { //sam cache in use + //prcess reads in cache + if (process_readcache_iter(c) < 0) + return -3; + //get reads to be consumed + if (!getfrom_readcache_iter(c, s, &iter->curr_tid, &iter->curr_beg, &iter->curr_end)) { + if (ret == -1) { + iter->finished = 1; + } + return -1; + } else + ret = 0; + return ret; + } else { + iter->finished = 1; + } return ret; } diff --git a/htslib/hts.h b/htslib/hts.h index 1cffe62f7..db2a24aa4 100644 --- a/htslib/hts.h +++ b/htslib/hts.h @@ -260,6 +260,7 @@ typedef struct htsFile { const char *fnidx; struct sam_hdr_t *bam_header; struct hts_filter_t *filter; + void *c; //for cache } htsFile; // A combined thread pool and queue allocation size. @@ -332,6 +333,8 @@ enum hts_fmt_option { HTS_OPT_BLOCK_SIZE, HTS_OPT_FILTER, HTS_OPT_PROFILE, + HTS_OPT_CACHE_FLT_DEPTH, //depth above which to discard reads + HTS_OPT_CACHE_FLT_SIZE, //size of caching window // Fastq @@ -876,7 +879,7 @@ typedef int64_t hts_tell_func(void *fp); */ typedef struct hts_itr_t { - uint32_t read_rest:1, finished:1, is_cram:1, nocoor:1, multi:1, dummy:27; + uint32_t read_rest:1, finished:1, is_cram:1, nocoor:1, multi:1, usecache:1, dummy:26; int tid, n_off, i, n_reg; hts_pos_t beg, end; hts_reglist_t *reg_list; @@ -1316,6 +1319,9 @@ typedef hts_itr_t *hts_itr_query_func(const hts_idx_t *idx, int tid, hts_pos_t b HTSLIB_EXPORT hts_itr_t *hts_itr_querys(const hts_idx_t *idx, const char *reg, hts_name2id_f getid, void *hdr, hts_itr_query_func *itr_query, hts_readrec_func *readrec); +/// use caching when available, invoked only for sequence data to differentiate between sequence and variant/tabix usage +static inline hts_itr_t* hts_itr_usecache(hts_itr_t *itr) { if (itr) itr->usecache = 1; return itr; } + /// Return the next record from an iterator /** @param fp Input file handle @param iter Iterator diff --git a/sam.c b/sam.c index 89e537e9d..e8af4ef96 100644 --- a/sam.c +++ b/sam.c @@ -58,6 +58,7 @@ DEALINGS IN THE SOFTWARE. */ #include "htslib/hts_expr.h" #include "header.h" #include "bgzf_internal.h" +#include "sam_cache.h" #include "htslib/khash.h" KHASH_DECLARE(s2i, kh_cstr_t, int64_t) @@ -1161,6 +1162,8 @@ static int sam_readrec(BGZF *ignored, void *fpv, void *bv, int *tid, hts_pos_t * htsFile *fp = (htsFile *)fpv; bam1_t *b = bv; fp->line.l = 0; + //mark iterator access to cache, if in use + set_iter_access(fp); int ret = sam_read1(fp, fp->bam_header, b); if (ret >= 0) { *tid = b->core.tid; @@ -1176,6 +1179,8 @@ static int sam_readrec_rest(BGZF *ignored, void *fpv, void *bv, int *tid, hts_po htsFile *fp = (htsFile *)fpv; bam1_t *b = bv; fp->line.l = 0; + //mark iterator access to cache, if in use + set_iter_access(fp); int ret = sam_read1(fp, fp->bam_header, b); return ret; } @@ -1737,11 +1742,11 @@ hts_itr_t *sam_itr_queryi(const hts_idx_t *idx, int tid, hts_pos_t beg, hts_pos_ { const hts_cram_idx_t *cidx = (const hts_cram_idx_t *) idx; if (idx == NULL) - return hts_itr_query(NULL, tid, beg, end, sam_readrec_rest); + return hts_itr_usecache(hts_itr_query(NULL, tid, beg, end, sam_readrec_rest)); else if (cidx->fmt == HTS_FMT_CRAI) - return cram_itr_query(idx, tid, beg, end, sam_readrec); + return hts_itr_usecache(cram_itr_query(idx, tid, beg, end, sam_readrec)); else - return hts_itr_query(idx, tid, beg, end, sam_readrec); + return hts_itr_usecache(hts_itr_query(idx, tid, beg, end, sam_readrec)); } static int cram_name2id(void *fdv, const char *ref) @@ -1753,9 +1758,9 @@ static int cram_name2id(void *fdv, const char *ref) hts_itr_t *sam_itr_querys(const hts_idx_t *idx, sam_hdr_t *hdr, const char *region) { const hts_cram_idx_t *cidx = (const hts_cram_idx_t *) idx; - return hts_itr_querys(idx, region, bam_name2id_wrapper, hdr, + return hts_itr_usecache(hts_itr_querys(idx, region, bam_name2id_wrapper, hdr, cidx->fmt == HTS_FMT_CRAI ? cram_itr_query : hts_itr_query, - sam_readrec); + sam_readrec)); } hts_itr_t *sam_itr_regarray(const hts_idx_t *idx, sam_hdr_t *hdr, char **regarray, unsigned int regcount) @@ -1781,6 +1786,7 @@ hts_itr_t *sam_itr_regarray(const hts_idx_t *idx, sam_hdr_t *hdr, char **regarra itr = hts_itr_regions(idx, r_list, r_count, bam_name2id_wrapper, hdr, hts_itr_multi_bam, sam_readrec, bgzf_pseek, bgzf_ptell); } + hts_itr_usecache(itr); if (!itr) hts_reglist_free(r_list, r_count); @@ -1796,11 +1802,11 @@ hts_itr_t *sam_itr_regions(const hts_idx_t *idx, sam_hdr_t *hdr, hts_reglist_t * return NULL; if (cidx->fmt == HTS_FMT_CRAI) - return hts_itr_regions(idx, reglist, regcount, cram_name2id, cidx->cram, - hts_itr_multi_cram, cram_readrec, cram_pseek, cram_ptell); + return hts_itr_usecache(hts_itr_regions(idx, reglist, regcount, cram_name2id, cidx->cram, + hts_itr_multi_cram, cram_readrec, cram_pseek, cram_ptell)); else - return hts_itr_regions(idx, reglist, regcount, bam_name2id_wrapper, hdr, - hts_itr_multi_bam, sam_readrec, bgzf_pseek, bgzf_ptell); + return hts_itr_usecache(hts_itr_regions(idx, reglist, regcount, bam_name2id_wrapper, hdr, + hts_itr_multi_bam, sam_readrec, bgzf_pseek, bgzf_ptell)); } /********************** @@ -4255,11 +4261,32 @@ static inline int sam_read1_sam(htsFile *fp, sam_hdr_t *h, bam1_t *b) { // Returns 0 on success, // -1 on EOF, // <-1 on error -int sam_read1(htsFile *fp, sam_hdr_t *h, bam1_t *b) +int sam_read1(htsFile *fp, sam_hdr_t *h, bam1_t *r) { int ret, pass_filter; + rc_t *c = NULL; + ce_t *e = NULL; + bam1_t *b = r; + + //if cache is in use and is not invoked thr' iterators, handle cache here itself + //otherwise handle in itr_nxt - no cache handling here! + if (!get_iter_access(fp)) { + c = (rc_t*)fp->c; + } + if (c) { //try to get cached reads + if ((ret = getfrom_readcache(c, r, NULL)) > 0) { + return 0; + } else if (ret < 0) + return -1; + //nothing cached or not ready yet + } do { + if (c) { //get cached storage + if (!(e = get_cache(fp))) + return -4; + b = e->r; //get bam record from storage + } switch (fp->format.format) { case bam: ret = sam_read1_bam(fp, h, b); @@ -4297,6 +4324,34 @@ int sam_read1(htsFile *fp, sam_hdr_t *h, bam1_t *b) pass_filter = (ret >= 0 && fp->filter) ? sam_passes_filter(h, b, fp->filter) : 1; + + if (c) { + //cache in use, add to cache if it is passed filtering + if (pass_filter) { + pass_filter = 0; + if (ret >= 0) { //successfull read + if (addto_readcache(c, e, NULL)) { + return -4; + } + } else { + if (ret == -1) { //end + notify_end(c, e); + } + } + if (get_readcache_status(c) >= WNDFULL) { //end/window full/ready + if (process_readcache(c) < 0) + return -4; + + pass_filter = getfrom_readcache(c, r, NULL); + if (-1 == ret && !pass_filter) { + pass_filter = 1; + } + else + ret = 0; + } + } else //return storage + ret_cache(c, e); + } } while (pass_filter == 0); return pass_filter < 0 ? -2 : ret; diff --git a/sam_cache.c b/sam_cache.c new file mode 100644 index 000000000..0ae6e6695 --- /dev/null +++ b/sam_cache.c @@ -0,0 +1,1069 @@ +/* sam_cache.c -- Functions to create a cache of reads for depth handling + + Copyright (C) 2026 Genome Research Ltd. + + Author: Vasudeva Sarma + +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. */ + +#include + +#include "sam_cache.h" +#include "htslib/hts_alloc.h" + + +#ifdef CACHE_DBG_LOG +FILE *cachelog = NULL; +#endif //CACHE_DBG_LOG + +#define CHUNK 1024 //no of items in a cache allocation + +/// @brief ensures depth buffer is big enough +/// @param c pointer to read cache +/// @param sz required size of depth buffer +/// @return 0 on success -ve on failure +static int ensure_depthbuffer(rc_t *c, hts_pos_t sz) +{ + if (c->dp_sz && sz <= c->dp_sz) + return 0; //big enough + int *dpth = hts_realloc_p(c->dpth, sizeof(int), sz); + if (!dpth) + return -1; + c->dpth = dpth; + c->dp_sz = sz; + return 0; +} +/// @brief setup read cache +/// @param fp htsFile pointer to which cache to be attached +/// @param wndsz size of cache window +/// @param maxdpth depth limit +/// @return 0 on success and non-zero on failure +int setup_readcache(htsFile *fp, int wndsz, int maxdpth) +{ + int i, j, first = 0; + rc_t *c = (rc_t*)fp->c; + ce_t *elem = NULL, *tail = NULL, **p = NULL; + const int def_wndsz = 3500, def_dpth = 1000; + /*create cache if it doesn't exists. set window and depth size during + initialisation. when free cache slots are 1, allocate next chunk.*/ + if (!c) { //create cache + first = 1; + wndsz = wndsz <= 0 ? def_wndsz : wndsz; + maxdpth = maxdpth <= 0 ? def_dpth : maxdpth; + if (!(c = hts_calloc(sizeof(rc_t), 1))) + goto fail; + fp->c = c; + } else if (c->cache.f > 1) { //re-init or retrieval + //update depth buffer / dpth settings if needed + if (wndsz || maxdpth) { + if (wndsz > 0) { //succeeding window size arg + //avoid retaining unnecessary buffer if wndsz is less than set + if (wndsz < c->wndsz) { + free(c->dpth); + c->dpth = NULL; + c->dp_sz = 0; + } + if (ensure_depthbuffer(c, wndsz+1)) + goto fail; + memset(c->dpth, 0, (wndsz+1) * sizeof(int)); + c->wndsz = wndsz; + } + if (maxdpth > 0) { //succeeding depth arg + c->maxdpth = maxdpth; + } +#ifdef CACHE_DBG_LOG + LG("setup_readcache wnd:%d dpth:%d\n", c->wndsz, c->maxdpth); +#endif //CACHE_DBG_LOG + } + return 0; + } + //make cache storage + if (!c->cache.m) { //initial + p = hts_calloc_ps(sizeof(ce_t*), c->cache.n, 1); + } else { //growing + p = hts_realloc_ps(c->cache.p, sizeof(ce_t*), c->cache.n, 1); + } + if (!p) + goto fail; + c->cache.p = p; //array holding chunks of storage + //allocate cache elements + if((elem = hts_calloc(sizeof(ce_t), CHUNK))) { + c->cache.p[c->cache.n++] = elem; + if (!c->cache.m) { //initial + if (!(elem->r = bam_init1())) + goto fail; + c->cache.head = tail = elem; + i = 1; +#ifdef CACHE_DBG_LOG + ks_initialize(&elem->log); +#endif //CACHE_DBG_LOG + } else { //growing + tail = c->cache.tail; + i = 0; + } + //initialize and add to tail + for (; i < CHUNK; ++i) { + if (!((elem + i )->r = bam_init1())) + goto fail; + tail->next = elem + i; + tail = tail->next; +#ifdef CACHE_DBG_LOG + ks_initialize(&elem->log); +#endif //CACHE_DBG_LOG + } + c->cache.m += CHUNK; + c->cache.f += CHUNK; + c->cache.tail = tail; + } else + goto fail; + + if (first) { //setup starting params + c->w_st = c->w_en = -1; + c->tid = -2; + if (!c->selpair && !(c->selpair = kh_init(pair))) + goto fail; + //todo check the +1 allocations + if (wndsz) { //window size arg in use + if (ensure_depthbuffer(c, wndsz+1)) + goto fail; + memset(c->dpth, 0, (wndsz+1) * sizeof(int)); + c->wndsz = wndsz; + } + if (maxdpth) { //depth arg in use + c->maxdpth = maxdpth; + } +#ifdef CACHE_DBG_LOG + LG("setup_readcache wnd:%d dpth:%d\n", c->wndsz, c->maxdpth); +#endif //CACHE_DBG_LOG + } + + return 0; + +fail: + if (c) { + for (i = 0; i < c->cache.n; ++i) { + elem = c->cache.p[i]; + for (j = 0; j < CHUNK; ++j) { + bam_destroy1(elem[j].r); +#ifdef CACHE_DBG_LOG + ks_free(&(elem[j].log)); +#endif //CACHE_DBG_LOG + } + free(elem); + c->cache.p[i] = NULL; + } + free(c->cache.p); + c->cache.p = NULL; + free(c->dpth); + c->dpth = NULL; + free(c); + fp->c = NULL; + } + return 1; +} +/// @brief destroys the read cache +/// @param fp htsFile pointer +void destroy_readcache(htsFile *fp) +{ + int i, j; + ce_t *elem = NULL; + rc_t *c = (rc_t*) fp->c; + khint_t iter; + + if(!c) //cache not in use + return; + + for (iter = kh_begin(c->selpair); iter != kh_end(c->selpair); ++iter) { + if (kh_exist(c->selpair, iter)) { + kh_del(pair, c->selpair, iter); + } + } + + while (c->head) { //clear remaining reads + elem = c->head->next; + LGlog(&c->head->log, "%s", "cleanup"); + ret_cache(c, c->head); + c->head = elem; + } + while (c->head_nsel) { //clear non-selected reads + elem = c->head_nsel->next; + LGlog(&c->head_nsel->log, "%s", "cleanup"); + ret_cache(c, c->head_nsel); + c->head_nsel = elem; + } + //cleanup cache + for (i = 0; i < c->cache.n; ++i) { + elem = c->cache.p[i]; + for (j = 0; j < CHUNK; ++j) { + bam_destroy1(elem[j].r); +#ifdef CACHE_DBG_LOG + ks_free(&elem[j].log); +#endif //CACHE_DBG_LOG + } + free(elem); + } + free(c->cache.p); + free(c->dpth); + free(c->inc); + kh_destroy(pair, c->selpair); + free(c); + fp->c = NULL; +} + +//implementation / internals +/// @brief marks the access as thr' iterator +/// @param fp htsFile pointer +void set_iter_access(htsFile *fp) { + /*when cache is used thr' iterator, cached read handling is done in itr_nxt + when it is used on whole file, it is done in sam_read1. this flag + helps to identify these scenarios and use cache appropriately*/ + if (fp->c) { //cache is in use + rc_t *c = (rc_t*)fp->c; + c->itr = 1; + } +} +/// @brief retrieves the how the cache is accessed, thr' iterator/for whole file +/// or not in use at all +/// @param fp htsFile pointer for cache access +/// @return 1 if thr' iterator and 0 if not in use / not thr' iterator +int get_iter_access(htsFile *fp) { + if (fp->c) { //cache in use + rc_t *c = (rc_t*) fp->c; + return c->itr; + } + return 0; //cache not in use +} +/// @brief cache's status +/// @param c read cache +/// @return cache status enum showing status of cache +cs get_readcache_status(rc_t *c) { + return c ? c->sts : NOTREADY; +} + +//todo htsopt3 to try +//-1 on failure and 1 when required and 0 on skip +/// @brief update the depth buffer based on reads length +/// @param c read cache +/// @param e cache element holding the read under processing +/// @param chk 1 checks whether the read is required or not; 0 to update depth +/// @return -ve - error, 0 - read not required, 1 - read required +static int update_depth(rc_t *c, ce_t *e, int chk) +{ + int *dpth = NULL; + uint32_t *cgr = bam_get_cigar(e->r), i, j; + int clen, off, req = 0; + hts_pos_t st, en, len = 0, h; + if (!c->dpth) { //setup depth buffer + if (ensure_depthbuffer(c, c->dp_sz + 1)) + goto fail; + c->dp_en = c->w_st + c->dp_sz; + off = 0; + } + st = c->w_st; en = c->dp_en; + if (st > e->r->core.pos) + goto fail; //not sorted? + //inc buffer is an attempt to force compiler to use intrinsics + if (e->len > c->inc_sz) { //grow increment buffer as required + int *inc = hts_realloc_p(c->inc, sizeof(int), e->len); + if (!inc) goto fail; + c->inc = inc; + c->inc_sz = e->len; + for(h = 0; h < c->inc_sz; ++h) + *(c->inc + h) = 1; + } + if (st < e->r->core.pos) { + off = e->r->core.pos - st; + len = off + e->len; + } else { + off = st - e->r->core.pos; + len = e->len - off; + off *= -1; + } + if (en < (c->w_st+len)) { //goes over the end of buffer, increase it + len = c->w_st + len - en; + hts_pos_t bkpsz = c->dp_sz; + if (ensure_depthbuffer(c, len + c->dp_sz)) + goto fail; + memset(c->dpth + bkpsz, 0, len * sizeof(int)); + c->dp_en = en = c->w_st + c->dp_sz; + } + len = 0; + dpth = c->dpth; + if (chk) { //check depth + for (i = 0; i < e->r->core.n_cigar; ++i) { + if (!(bam_cigar_type(bam_cigar_op(cgr[i])) & 2)) { + continue; //not consuming ref + } + //deletion is counted! + //check depth for each position is above required limit or not + //read required if depth is <= the limit + clen = bam_cigar_oplen(cgr[i]); + for (j = 0; j < clen; ++j) { + if (off + j >= 0) + req |= dpth[off + j] + 1 <= c->maxdpth; + } + off += clen; + if (req) { //required, no need to check further + break; + } + } + } else { //update depth + req = 1; + for (i = 0; i < e->r->core.n_cigar; ++i) { + if (!(bam_cigar_type(bam_cigar_op(cgr[i])) & 2)) { //not consuming ref + continue; + } + clen = bam_cigar_oplen(cgr[i]); + //inc buffer is an attempt to force compiler to use intrinsics + for (j = 0; j < clen; ++j) { + dpth[off + len + j] += c->inc[j]; + } + len += clen; + } + } + if (req) + return 1; + return 0; + +fail: + return -1; +} +/// @brief return the storage back to cache +/// @param c read cache +/// @param elem element/space in cache +void ret_cache(rc_t *c, ce_t* elem) +{ + if ((elem->r->core.flag & BAM_FPAIRED) && !(elem->r->core.flag & BAM_FUNMAP)) { + //paired and mate mapped, remove from expected pair + khiter_t it = kh_get(pair, c->selpair, bam_get_qname(elem->r)); + if (it != kh_end(c->selpair) && kh_exist(c->selpair, it)) { + kh_del(pair, c->selpair, it); + } + } + //add as head in cache + elem->prev = NULL; + elem->ord = 0; + elem->len = 0; + elem->next = c->cache.head; + c->cache.head->prev = elem; + c->cache.head = elem; + ++c->cache.f; +#ifdef CACHE_DBG_LOG + ks_clear(&elem->log); +#endif //CACHE_DBG_LOG +} + +/// @brief get cache element / storage from preallocated cache +/// @param fp htsFile pointer to setup / retrieve cache +/// @return ce_t* on success or NULL on failure +ce_t* get_cache(htsFile *fp) +{ + rc_t *c = (rc_t*)fp->c; + ce_t *ret = NULL; + + //create enough space + if (setup_readcache(fp, 0, 0)) //passing 0 to avoid re-initialization + goto fail; + + ret = c->cache.head; + c->cache.head = c->cache.head->next; + c->cache.head->prev = NULL; + ret->prev = NULL; + ret->next = NULL; + --c->cache.f; + return ret; +fail: + return NULL; +} +/// @brief mark end of input +/// @param p read cache pointer +/// @param e allocated and unused space +void notify_end(void *p, void *e) +{ + rc_t *c = (rc_t*)p; + c->sts = END; //end + c->tid = -3; //reset that it doesn't match to any/initial vals + //real end of input, unlike iterator where it could be just end of a region + ret_cache(c, (ce_t*)e); +} + +/// @brief add a read to cache +/// @param c pointer to read cache +/// @param e cache element containing the read to be cached +/// @param sts to return status of cache post caching +/// @return -1 on failure 0 on success +int addto_readcache(rc_t *c, ce_t *e, cs *sts) +{ + int unmap = 0; + if (!(e->r->core.flag & BAM_FUNMAP)) { + if (c->w_st == -1) { + //starting, use pos of 1st or one being added as start of window + c->w_st = c->head ? c->head->r->core.pos : e->r->core.pos; + c->w_en = c->w_st + c->wndsz; //end of wnd + c->dp_en = c->w_st + c->dp_sz; //end of depth buffer + } + } else { + unmap = 1; //unmapped, add w/o depth check + } + e->ord = ++(c->ord); + e->len = bam_cigar2rlen(e->r->core.n_cigar, bam_get_cigar(e->r)); + LG("+ %s %"PRIu64"\t\t%"PRIhts_pos" %"PRIhts_pos" %"PRIu64" %"PRIhts_pos"\n", bam_get_qname(e->r), e->ord, c->w_st, e->r->core.pos, e->len, c->w_en); + if (!c->head) { + c->head = c->tail = e; + } else { + ce_t *p = c->tail; + ce_t *tmpn = NULL; + //add to the tail + if (!unmap) { + if (p->r->core.tid == e->r->core.tid && p->r->core.pos > e->r->core.pos) { + hts_log_error("Unsorted data"); + return -1; //not sorted! + } + } + if (p) { //useful if read is sorted based on len and being inserted + tmpn = p->next; + p->next = e; + e->prev = p; + e->next = tmpn; + if (tmpn) + tmpn->prev = e; + if (p == c->tail) + c->tail = e; + } else { //either last or 1st + if (c->head == c->tail && !c->tail) { //none in list + c->tail = c->head = e; + e->prev = e->next = NULL; + } else { //insert 1st + tmpn = c->head; + c->head = e; + e->prev = NULL; + e->next = tmpn; + if(tmpn) + tmpn->prev = e; + } + } + } + //todo do we need a limit on max no of items that are cached? like the whole file is for same pos, probably cant be loaded! + if (c->tid == e->r->core.tid) { + if (c->w_en < e->r->core.pos) { //post window, process and advance + LG("wnd full\n"); + c->sts = WNDFULL; //wnd full, go for processing + } + else + c->sts = CACHING; //caching + } else if (c->tid != -2) { + LG("tid change\n"); + c->sts = READY; //ready for processing + } + else + c->sts = CACHING; //caching + + c->tid = e->r->core.tid; + if (sts) + *sts = c->sts; + LGlog(&e->log, "%s,%"PRIu64",added,%d,%d,%"PRIhts_pos",%"PRIhts_pos",%"PRIu64",%"PRIhts_pos",", bam_get_qname(e->r), e->ord,e->r->core.tid, e->r->core.flag,e->r->core.pos, e->r->core.mpos,e->len, e->len+e->r->core.pos); + return 0; +} + +/// @brief get read from processed cache +/// @param c pointer to read cache +/// @param b pointer to bam data, to which read data is copied +/// @param end end of read, for iterators +/// @return -1 on failure, 0 when nothing to retrieve and 1 with read retrieved +int getfrom_readcache(rc_t *c, bam1_t *b, hts_pos_t *end) +{ + if (!c || c->sts < READY) { //not ready! + return 0; + } + //todo at some point, removal from selpair need to be done based on pos as well + uint64_t sel = UINT64_MAX, ins = UINT64_MAX; + ce_t *e = c->head_sel, *f = c->head_ins, *p = NULL; + + //get from selected or inserted list, based on ordinal + if (e) + sel = e->ord; + if (f) + ins = f->ord; + if(sel < ins) + p = e; + else + p = f; + + if (p && (c->sts == READY || c->sts == END)) { + //send only upto start of wnd to maintain the order, except when it is end + if (!bam_copy1(b, p->r)) + return -1; + if (p == e) { //remove from sel list + c->head_sel = p->next; + if (!c->head_sel) { + c->tail_sel = NULL; + } + } else { //remove from ins list + c->head_ins = p->next; + if (!c->head_ins) { + c->tail_ins = NULL; + } + } + if (!c->head_sel && !c->head_ins && c->sts != END) { + c->sts = NOTREADY; //not ready + } else { + if (c->head_sel) hts_prefetch(c->head_sel); + if (c->head_ins) hts_prefetch(c->head_ins); + } + LG("- %s %"PRIu64"\n", bam_get_qname(b), p->ord); + LGlog(&p->log, "%s", ",retrieved"); + if (end) *end = p->len + p->r->core.pos; + ret_cache(c, p); //return storage to cache + return 1; + } + return 0; +} +/// @brief find a read matching to given one from non-selected list +/// @param c pointer to read cache +/// @param e pointer to read for which pair need to be found +/// @param ep pointer to previous one of the pair, to maintain list +/// @return NULL when not found and cache element pointer when found +static inline ce_t* find_nsel(rc_t *c, ce_t *e, ce_t **ep) +{ + ce_t *s = c->head_nsel; + *ep = NULL; + while (s) { + if (s->next) + hts_prefetch(s->next); + if (s->ord > e->ord) + break; //not found + if (s->r->core.pos == e->r->core.mpos && + s->r->core.mpos == e->r->core.pos && + s->r->core.tid == e->r->core.mtid && + s->r->core.mtid == e->r->core.tid && + !strcmp(bam_get_qname(s->r), bam_get_qname(e->r))) + return s; //found + *ep = s; + s = s->next; + } + return NULL; +} +/// @brief move the read/cache element from main list to selected/unselected/insert list +/// @param c pointer to cache +/// @param ep pointer to previous element to maintain the list +/// @param e element being moved +/// @param en next element +/// @param sel 1 to move read to selected list 0 to move to nselected list +/// @param ins 1 to move to insert list, relevant with sel = 1 +static inline void move_read(rc_t *c, ce_t *ep, ce_t *e, ce_t* en, int sel, int ins) +{ + int paired = (e->r->core.flag & BAM_FPAIRED) && + !(e->r->core.flag & BAM_FUNMAP) && !(e->r->core.flag & BAM_FMUNMAP) && + (e->r->core.mtid != -1) && (e->r->core.mpos != -1); + if (!ins) { //remove from cache + if (en) + en->prev = ep; + if (c->head == e) + c->head = en; + if (c->tail == e) + c->tail = ep; + if (!c->head) + c->tail = c->head; + if(ep) { + ep->next = en; + } + } else { //remove from nsel + if (ep) { + ep->next = en; + } else { + c->head_nsel = en; + } + if (en) + en->prev = ep; + else + c->tail_nsel = ep ? ep : NULL; + + } + e->next = NULL; + e->prev = NULL; + + if (sel) { //moving to sel/ins list + //insert in required pos, starting from tail, in order of ordinal + ce_t *s = ins? c->tail_ins : c->tail_sel, *p = NULL; + if (s && s->ord < e->ord) { //shortcut + s->next = e; + e->next = NULL; + e->prev = s; + if(ins) + c->tail_ins = e; + else + c->tail_sel = e; + return; + } + while (s) { + if (s->ord < e->ord) { //add in ascending order + break; + } + s = s->prev; + } + if (!s) { //as head + if (ins) { + p = c->head_ins; + c->head_ins = e; + } + else { + p = c->head_sel; + c->head_sel = e; + } + e->prev = NULL; + e->next = p; + if(p) + p->prev = e; + else { //update tail + if (ins) + c->tail_ins = e; + else + c->tail_sel = e; + } + return; + } else { + p = s->next; + s->next = e; + e->prev = s; + e->next = p; + if (p) + p->prev = e; + return; + } + return; + } else if (paired) { + //move to nsel if paired, otherwise discard and return cache + //add to non-selected list, for pair lookup + ce_t *s = c->tail_nsel, *p = NULL; + if (s && s->r->core.pos < e->r->core.pos && + s->r->core.tid == e->r->core.tid) { //shortcut + s->next = e; + e->next = NULL; + e->prev = s; + c->tail_nsel = e; + return; + } + //find pos and fit, in order of increasing pos, that it is easy to remove + while (s && (s->r->core.tid == e->r->core.tid)) { + if (s->r->core.pos < e->r->core.pos) { + break; + } + s = s->prev; + } + if (!s) { //add as head + p = c->head_nsel; + c->head_nsel = e; + e->prev = NULL; + e->next = p; + if(p) + p->prev = e; + if (!p) + c->tail_nsel = e; + return; + } else { + p = s->next; + s->next = e; + e->prev = s; + e->next = p; + if (p) + p->prev = e; + return; + } + return; + } else { //non selected, non paired reads, release them + LGlog(&e->log, "%s", "npair,disc"); + if(ep) + ep->next = en; + if (en) + en->prev = ep; + ret_cache(c, e); + return; + } +} +/// @brief reset cache status, for next tid/iterator... +/// @param c pointer to read cache +static inline void reset_depth(rc_t* c) +{ + c->w_st = -1; + if (c->dp_sz <= 0 || !c->dpth) + return; + memset(c->dpth, 0, c->dp_sz * sizeof(int)); + + ce_t *en = NULL; + //clear all from previous tid + while (c->head && c->head->r->core.tid != c->tid) { + en = c->head->next; + LGlog(&c->head->log, "%s", "h-reset"); + ret_cache(c, c->head); + c->head = en; + } + if (!c->head) c->tail = NULL; + else c->head->prev = NULL; + + //clear whole non selected ones + while (c->head_nsel) { + en = c->head_nsel->next; + LGlog(&c->head_nsel->log, "%s", "n-reset"); + ret_cache(c, c->head_nsel); + c->head_nsel = en; + } + + LG("reset: t %"PRIu64" s %"PRIu64" i %"PRIu64" n %"PRIu64"; nxt %d\n", c->rcnt, c->selcnt,c->inscnt, c->nselcnt, c->tid); + c->tail_nsel = NULL; +} +/// @brief process the cached reads and find required ones +/// @param c read cache +/// @return 0 on success and -ve on error +int process_readcache(rc_t *c) +{ + ce_t *e = NULL, *ep = NULL, *en = NULL; + hts_pos_t pos, off; + khiter_t pairitr; + if (!c->head) + return 0; + + hts_pos_t endpos = c->w_en < c->tail->r->core.pos ? c->tail->r->core.pos - 1 : c->w_en; + hts_pos_t lastpos = 0; + + if(c->sts == END) { + //fine tune endpos for last iteration, by looking for valid len which may not be the tail one! + e = c->head; + lastpos = c->head->r->core.pos + c->head->len; + while (e) { + pos = e->r->core.pos + e->len; + if (lastpos < pos) + lastpos = pos; + e = e->next; + } + if (endpos < lastpos) { + endpos = lastpos - 1; + } + } + pos = c->w_st; + while (pos <= endpos) { + if (!(e = c->head)) + break; + ep = NULL; + //discard any irrelevant ones + while (e && (e->r->core.pos <= pos) && ((c->tid != e->r->core.tid) || c->sts == WNDFULL)) { + en = e->next; + if (e->r->core.flag & BAM_FUNMAP) { //unmapped, nothing further to check + LGlog(&e->log, "%s%"PRIhts_pos, "sel-unmapped,", pos); + LG("* x %"PRIhts_pos" selunmapped\n", e->ord); + move_read(c, ep, e, en, 1, 0); + } else if (e->r->core.pos + e->len - 1 < pos) { + //not relevant for this pos or succeeding ones + //check whether pair is selected before discarding + pairitr = kh_get(pair, c->selpair, bam_get_qname(e->r)); + if (pairitr != kh_end(c->selpair)) { + pair_exp *p = &kh_val(c->selpair, pairitr); + if (p->mpos == e->r->core.pos && + p->mtid == e->r->core.tid && + p->pos == e->r->core.mpos && + p->tid == e->r->core.mtid) { //pair already selected + kh_del(pair, c->selpair, pairitr); //remove from expected pairs + LGlog(&e->log, "%s%"PRIhts_pos, "sel as paired,", pos); + move_read(c, ep, e, en, 1, 0); //select + LG("* x %"PRIhts_pos" selpaired\n", e->ord); + //no depth update! + } else { + LGlog(&e->log, "%s%"PRIhts_pos, "nsel,", pos); + move_read(c, ep, e, en, 0, 0); //move to unselected + LG("* x %"PRIhts_pos" nsel\n", e->ord); + } + } else { + LGlog(&e->log, "%s%"PRIhts_pos",", "nsel,", pos); + move_read(c, ep, e, en, 0, 0); //move to unselected + LG("* x %"PRIhts_pos" nsel2\n", e->ord); + } + } else { + ep = e; + } + e = en; + } + e = c->head; + ep = NULL; + off = pos - c->w_st; + if (c->head == c->tail && c->sts == WNDFULL) { + //lastone --> all from wnd are done and last one to be considered in nxt iteration + endpos = pos; + break; + } else { + if (c->dp_sz <= off) { + hts_pos_t bkp = c->dp_sz, ln = 100; + if (ensure_depthbuffer(c, c->dp_sz + ln)) { + goto fail; + } + memset(c->dpth + bkp, 0, ln * sizeof(int)); + c->dp_en += ln; + } + } + if (c->dpth[off] >= c->maxdpth) { //have enough depth + ++pos; + continue; + } + //find the last one covering the pos + while (e && (e->r->core.pos <= pos) && ((c->tid != e->r->core.tid) || c->sts == WNDFULL)) { + en = e->next; + if (e->r->core.flag & BAM_FUNMAP) { + LGlog(&e->log, "%s%"PRIhts_pos, "sel-unmapped,", pos); + LG("* x %"PRIhts_pos" selunmapped\n", e->ord); + move_read(c, ep, e, en, 1, 0); + e = en; + continue; + } + if (e->r->core.pos + e->len - 1 >= pos) { + ep = e; + } + e = en; + } + if (!ep) { //nothing! + ++pos; + continue; + } + LG("* x %"PRIhts_pos" sel @ %"PRIhts_pos"\n", ep->ord, pos); + LGlog(&ep->log, "%s%"PRIhts_pos, "sel,", pos); + move_read(c, ep->prev, ep, ep->next, 1, 0); + if (update_depth(c, ep, 0) < 0) { + goto fail; + } + //get/set for pair + if (ep->r->core.flag & BAM_FPAIRED && !(ep->r->core.flag & BAM_FMUNMAP)) { + if (ep->r->core.mpos >= pos) { //upcoming mate, add to hash + int r = -1; + pairitr = kh_put(pair, c->selpair, bam_get_qname(ep->r), &r); + if (r == -1) + goto fail; + pair_exp *p = &kh_val(c->selpair, pairitr); + p->pos = ep->r->core.pos; p->tid = ep->r->core.tid; + p->mpos = ep->r->core.mpos; p->mtid = ep->r->core.mtid; + } else { //mate already passed, selected or not? + int sel = 0; + pairitr = kh_get(pair, c->selpair, bam_get_qname(ep->r)); + if (pairitr != kh_end(c->selpair)) { + pair_exp *p = &kh_val(c->selpair, pairitr); + if (p->mpos == ep->r->core.pos && + p->mtid == ep->r->core.tid && + p->pos == ep->r->core.mpos && + p->tid == ep->r->core.mtid) { //pair already selected + sel = 1; + } + } + if (!sel) { //find and select from unselected ones + ce_t *o = NULL, *op = NULL; + if ((o = find_nsel(c, ep, &op))) { + LGlog(&o->log, "%s%"PRIhts_pos",%"PRIhts_pos",", "sel from nsel as paired ", ep->ord, pos); + move_read(c, op, o, o->next, 1, 1); //select + //no depth update! + } else { + LG("searching for pair failed, %"PRIu64"\n", ep->ord); + } + } + if (pairitr != kh_end(c->selpair)) { //remove from hash + kh_del(pair, c->selpair, pairitr); //remove from expected pairs + } + } + } + } + assert((c->head && c->tail) || (c->sts != WNDFULL)); + //check for any pairs discarded due to sufficient depth and select + LG("pos at lpexit %"PRIhts_pos", head %llu\n", pos, c->head?c->head->ord:0); + e = c->head; + ep = NULL; + if (c->head) + --pos; //out of loop --> ++pos + while (e && (e->r->core.pos <= pos) && ((c->tid != e->r->core.tid) || c->sts == WNDFULL)) { + en = e->next; + pairitr = kh_get(pair, c->selpair, bam_get_qname(e->r)); + if (pairitr != kh_end(c->selpair)) { + pair_exp *p = &kh_val(c->selpair, pairitr); + if (p->mpos == e->r->core.pos && + p->mtid == e->r->core.tid && + p->pos == e->r->core.mpos && + p->tid == e->r->core.mtid) { //pair already selected + kh_del(pair, c->selpair, pairitr); //remove from expected pairs + LGlog(&e->log, "%s%"PRIhts_pos, "sel as pair on wnd move,", pos); + move_read(c, ep, e, en, 1, 0); //select + //no depth update! + } else { + ep = e; + } + } else { + ep = e; + } + e = en; + } + if (c->sts != WNDFULL) { //when it is not wnd full, tid change or end + reset_depth(c); + } else { //update window + en = NULL; + hts_pos_t adj = c->tail ? c->tail->r->core.pos - c->w_en : 0; //last one, out of window - current end + hts_pos_t new_st = c->w_st + adj; + hts_pos_t bkp_st = c->w_st; + uint64_t last = c->tail_sel ? c->tail_sel->ord : 0; + uint64_t lasti = c->tail_ins ? c->tail_ins->ord : 0; + if (lasti > last) + last = lasti; + int rem = 0; + if (c->head_nsel) { + while ( c->head_nsel->r->core.pos + c->head_nsel->len - 1 < new_st || c->head_nsel->ord < last) { + rem = 1; + en = c->head_nsel->next; + LG("* nsel discarded %s %"PRIu64"\n", bam_get_qname(c->head_nsel->r), c->head_nsel->ord); + LGlog(&c->head_nsel->log,"%s","nsel-disc"); + ret_cache(c, c->head_nsel); + if(!(c->head_nsel = en)) { + c->tail_nsel = NULL; + break; + } else { + c->head_nsel->prev = NULL; + } + } + } + if (c->head) {////// + while ( c->head->r->core.pos + c->head->len - 1 < new_st || c->head->ord < last) { + rem = 1; + en = c->head->next; + LG("* discarded %s %"PRIu64" last %"PRIu64" new_st %"PRIhts_pos"\n", bam_get_qname(c->head->r), c->head->ord, last, new_st); + LGlog(&c->head->log,"%s,%"PRIhts_pos",%s", "wndchange", pos,"disc"); + ret_cache(c, c->head); + if(!(c->head = en)) { + c->tail = NULL; + break; + } else { + c->head->prev = NULL; + } + } + } + if (rem) { + LG("* wnd full, removed items from head_nsel\n"); + } + else { + LG("* wnd full, 0 removed items from head_nsel, [%"PRIhts_pos"-%"PRIhts_pos"] %"PRIhts_pos"\n", c->w_st, c->w_en, c->head_nsel?c->head_nsel->r->core.pos : 0); + } + c->w_st = c->head ? c->head->r->core.pos : new_st; //move wnd + c->w_en = c->w_st + c->wndsz; + adj = c->w_st - bkp_st; + if (adj >= c->dp_sz) { + memset(c->dpth, 0, c->dp_sz * sizeof(int)); + c->dp_en = c->w_st + c->dp_sz; + } else { + LG("adj %"PRIhts_pos", mv %"PRIhts_pos"-%"PRIhts_pos",", adj, c->w_st+adj, c->w_st+c->dp_sz); + LG("0 set %"PRIhts_pos" - %"PRIhts_pos"\n", c->w_st+c->dp_sz-adj,c->w_st+c->dp_sz); + memmove(c->dpth, c->dpth + adj, (c->dp_sz - adj) * sizeof(int)); + memset(c->dpth + c->dp_sz - adj, 0, adj * sizeof(int)); + c->dp_en += adj; + } + LG("* wnd moved, %"PRIhts_pos" - %"PRIhts_pos", dpth %"PRIhts_pos" - %"PRIhts_pos"; s %"PRIu64" i %"PRIu64" ns %"PRIu64"\n", c->w_st, c->w_en, c->w_st, c->dp_en, c->selcnt, c->inscnt, c->nselcnt); + c->sts = READY; //reset full status n get already processed + } + + return 0; +fail: + return -1; +} + +//wrappers for iterators +//these wrappers help to avoid complexity in hts / iterator code by moving the +//cache structure access to this file + +/// @brief get read cache +/// @param itr iterator which is in use +/// @param data custom data with iterator (htsFile/kstring based on invoker) +/// @return read cache pointer +void* get_sam_readcache(hts_itr_t *itr, void *data) +{ + htsFile *fp = NULL; + /*invoked from itr_nxt, which is used by utilities like tabix as well. + cache is in use only for sam data and to identify the invocation usecache + flag is used. this flag is set when iterator is used in sam context.*/ + if (itr && itr->usecache) + fp = (htsFile*)data; + return fp ? fp->c : NULL; +} +/// @brief wrapper to get read from cache +/// @param p read cache pointer +/// @param s bam storage for output +/// @param tid tid of output read +/// @param beg beg of output read +/// @param end end of output read +/// @return -1 on failure, 0 when nothing to retrieve and 1 with read retrieved +int getfrom_readcache_iter(void *p, void *s, int *tid, hts_pos_t *beg, hts_pos_t* end) +{ + rc_t *c = (rc_t*)p; + bam1_t *b = (bam1_t*)s; + int ret = getfrom_readcache(c, b, end); + if (ret > 0) { + *tid = b->core.tid; + *beg = b->core.pos; + } + return ret; +} +/// @brief wrapper to get cache storage +/// @param data custom data for iterator, htsfile pointer for sam data +/// @return storage released from cache or NULL +void *get_cache_iter(void *data) +{ + htsFile *fp = (htsFile*)data; + void *p = get_cache(fp); + return p; +} +/// @brief retrieves bam pointer from cache storage +/// @param p cache storage retrieved +/// @return bam pointer as void * +void *get_readbuffer_iter(void *p) +{ + ce_t* e = (ce_t*)p; + return e->r; +} +/// @brief mark end of input, end of region or file +/// @param p read cache pointer +/// @param e cache storage in use +void notify_end_iter(void *p, void *e) +{ + rc_t *c = (rc_t*)p; + c->sts = END; //end + c->tid = -2; //reset as in start + ret_cache(c, (ce_t*)e); +} +/// @brief wrapper to add read to cache +/// @param c cache +/// @param s read to be added +/// @param sts status of cache, output +/// @return -1 on failure 0 on success +int addto_readcache_iter(void *c, void *s, cs *sts) +{ + return addto_readcache(c, s, sts); +} +/// @brief wrapper to process cached data +/// @param c read cache +/// @return 0 on success and -ve on failure +int process_readcache_iter(void *c) +{ + return process_readcache((rc_t*)c); +} +/// @brief wrapper to reset cache status +/// @param c read cache +void reset_readcache_iter(rc_t *c) +{ + reset_depth(c); + c->sts = NOTREADY; + c->tid = -2; + c->w_st = c->w_en = -1; +} diff --git a/sam_cache.h b/sam_cache.h new file mode 100644 index 000000000..8632fb8ee --- /dev/null +++ b/sam_cache.h @@ -0,0 +1,132 @@ +/* sam_cache.h -- Functions to create a cache of reads for depth handling + + Copyright (C) 2026 Genome Research Ltd. + + Author: Vasudeva Sarma + +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. */ +#ifndef HTSLIB_SAMCACHE_H +#define HTSLIB_SAMCACHE_H + +#include "htslib/sam.h" +#include "htslib/khash.h" + +#ifdef __cplusplus +extern "C" { +#endif + +typedef struct ce_t {//cache element + uint64_t ord; //ordinal + bam1_t *r; + struct ce_t *next, *prev; + hts_pos_t len; +#ifdef CACHE_DBG_LOG + kstring_t log; +#endif //CACHE_DBG_LOG +} ce_t; + +typedef struct cache_t {//cache + int f, m, n; //free elements, max elements, count of elem chunks (of 1024) + ce_t **p; //array holding chunks of elements + struct ce_t *head, *tail; +} cache_t; + +typedef struct pair_exp { + int mtid, tid; + hts_pos_t mpos, pos; +} pair_exp; +KHASH_MAP_INIT_STR(pair, pair_exp) + +/// cache status +typedef enum cs {NOTREADY = 0, CACHING, WNDFULL, READY, END} cs; + +typedef struct rc_t {//read cache + cache_t cache; //cache of mem space + ce_t *head, *tail; //alignments + ce_t *head_sel, *tail_sel; //selected alignments + ce_t *head_nsel, *tail_nsel; //non-selected alignments + ce_t *head_ins, *tail_ins; //inserted alignments + uint64_t ord; //last ordinal + cs sts; + int wndsz, maxdpth, itr; //size of cache window, depth limit, thr' iterator or not + hts_pos_t w_st, w_en, dp_en, inc_sz; //wnd start, end, dpth buffer end, size of inc. buffer + khash_t(pair) *selpair; //hash holding name of selected reads for pair selection + int dp_sz, tid; + int *inc; //buffer holding inc val (1), attempt to force intrinsics + int *dpth; //depth buffer +} rc_t; + +/// @brief setup cache +/// @param fp file pointer to which cache is assigned and used +/// @param wndsz size of cache window +/// @param maxdpth depth limit +/// @return 0 on success others on failure +int setup_readcache(htsFile *fp, int wndsz, int maxdpth); +void destroy_readcache(htsFile *fp); +// return an element to cache +void ret_cache(rc_t *c, ce_t* elem); +// get a cached storage from cache +ce_t* get_cache(htsFile *fp); +//notify end of read +void notify_end(void *c, void *e); +//add a read to cache +int addto_readcache(rc_t *c, ce_t *e, cs *sts); +//retrieve a selected read from cached ones +int getfrom_readcache(rc_t *c, bam1_t *b, hts_pos_t *end); +//process cached reads and select required ones +int process_readcache(rc_t *c); +//get status of cache +cs get_readcache_status(rc_t *c); +//wrapper / for iterators +//get/set access status - thr' iterator or not; if thr' iterator, cache handling is done in itr_nxt +void set_iter_access(htsFile *fp); +int get_iter_access(htsFile *fp); +//get cache pointer, for use in iterator +void* get_sam_readcache(hts_itr_t *itr, void *data); +//retrieve a selected read from cached ones, wrapper for iterator +int getfrom_readcache_iter(void *c, void *s, int *tid, hts_pos_t *beg, hts_pos_t* end); +// get a cached storage from cache, wrapper for iterator +void *get_cache_iter(void *data); +//retrieve bam storage from cache +void *get_readbuffer_iter(void *e); +//notify end of read, wrapper for iterator +void notify_end_iter(void *c, void *e); +//add a read to cache, wrapper for iterator +int addto_readcache_iter(void *c, void *s, cs *sts); +//process cached reads and select required ones, wrapper for iterator +int process_readcache_iter(void *c); +//resets cache status and depth buffer, end of tid/region +void reset_readcache_iter(rc_t *c); + +#ifdef CACHE_DBG_LOG +extern FILE *cachelog; +//this is closed by system on exit! +#define LG(...) {if (!cachelog) cachelog = fopen("/tmp/op","w"); if (cachelog) { fprintf(cachelog, __VA_ARGS__);}} +#define LGlog(s,...) ksprintf(s,__VA_ARGS__) +#else +#define LG(...) ; +#define LGlog(s,...) ; +#endif //CACHE_DBG_LOG + +#ifdef __cplusplus +} +#endif + +#endif //HTSLIB_SAMCACHE_H + diff --git a/test/cache.exp.1.T1T3T2.sam b/test/cache.exp.1.T1T3T2.sam new file mode 100644 index 000000000..541091f9f --- /dev/null +++ b/test/cache.exp.1.T1T3T2.sam @@ -0,0 +1,29 @@ +@HD VN:1.17 SO:coordinate +@SQ SN:T1 LN:40 +@SQ SN:T2 LN:40 +@SQ SN:T3 LN:40 +@CO @SQ SN* LN* AH AN AS DS M5 SP TP UR +@CO @RG ID* BC CN DS DT FO KS LB PG PI PL PM PU SM +@CO @PG ID* PN CL PP DS VN +@CO this is a dummy alignment file to demonstrate different abilities of hts apis +@CO QNAME FLAG RNAME POS MAPQ CIGAR RNEXT PNEXT TLEN SEQ QUAL [TAG:TYPE:VALUE]… +@CO 1234567890123456789012345678901234567890 +@CO AAAAACTGAAAACCCCTTTTGGGGACTGTTAACAGTTTTT T1 +@CO TTTTCCCCACTGAAAACCCCTTTTGGGGACTGTTAACAGT T2 +ITR1 99 T1 5 40 4M = 33 10 ACTG ()() +UNMP1 73 T1 21 40 3M * 0 5 GGG &&1 +A1 99 T1 25 35 6M = 31 8 ACTGTT ****** +B1 99 T1 25 35 6M = 31 8 GCTATT ****** +B5 355 T1 25 35 4M = 33 5 AGTG PPPP +A1 147 T1 31 33 6M = 25 -8 ACTGTT ()()() +ITR1 147 T1 33 37 4M = 5 -10 ACTG $$$$ +A5 99 T3 12 50 3M = 23 5 GAA ()( +B5 99 T3 12 50 3M = 23 5 GAT ()( +ITR3 147 T3 23 49 2M = 35 -10 TT ** +B5 147 T3 23 47 2M1X = 12 -5 TTG ((( +A5 147 T3 23 47 2M1X = 12 -5 TAG ((( +ITR3 99 T3 35 51 2M = 23 10 AA && +B4 99 T2 12 50 3M = 23 5 GAT ()( +ITR2 147 T2 23 49 2M = 35 -10 TT ** +B4 147 T2 23 47 2M1X = 12 -5 TAG ((( +ITR2 99 T2 35 51 2M = 23 10 AA && diff --git a/test/cache.exp.1.sam b/test/cache.exp.1.sam new file mode 100644 index 000000000..2deafeffb --- /dev/null +++ b/test/cache.exp.1.sam @@ -0,0 +1,31 @@ +@HD VN:1.17 SO:coordinate +@SQ SN:T1 LN:40 +@SQ SN:T2 LN:40 +@SQ SN:T3 LN:40 +@CO @SQ SN* LN* AH AN AS DS M5 SP TP UR +@CO @RG ID* BC CN DS DT FO KS LB PG PI PL PM PU SM +@CO @PG ID* PN CL PP DS VN +@CO this is a dummy alignment file to demonstrate different abilities of hts apis +@CO QNAME FLAG RNAME POS MAPQ CIGAR RNEXT PNEXT TLEN SEQ QUAL [TAG:TYPE:VALUE]… +@CO 1234567890123456789012345678901234567890 +@CO AAAAACTGAAAACCCCTTTTGGGGACTGTTAACAGTTTTT T1 +@CO TTTTCCCCACTGAAAACCCCTTTTGGGGACTGTTAACAGT T2 +ITR1 99 T1 5 40 4M = 33 10 ACTG ()() +UNMP1 73 T1 21 40 3M * 0 5 GGG &&1 +A1 99 T1 25 35 6M = 31 8 ACTGTT ****** +B1 99 T1 25 35 6M = 31 8 GCTATT ****** +B5 355 T1 25 35 4M = 33 5 AGTG PPPP +A1 147 T1 31 33 6M = 25 -8 ACTGTT ()()() +ITR1 147 T1 33 37 4M = 5 -10 ACTG $$$$ +B4 99 T2 12 50 3M = 23 5 GAT ()( +ITR2 147 T2 23 49 2M = 35 -10 TT ** +B4 147 T2 23 47 2M1X = 12 -5 TAG ((( +ITR2 99 T2 35 51 2M = 23 10 AA && +A5 99 T3 12 50 3M = 23 5 GAA ()( +B5 99 T3 12 50 3M = 23 5 GAT ()( +ITR3 147 T3 23 49 2M = 35 -10 TT ** +B5 147 T3 23 47 2M1X = 12 -5 TTG ((( +A5 147 T3 23 47 2M1X = 12 -5 TAG ((( +ITR3 99 T3 35 51 2M = 23 10 AA && +UNMP2 141 * 0 0 * * 0 7 AA && +UNMP3 77 * 0 0 * * 0 5 GGG &&2 diff --git a/test/cache.exp.2.sam b/test/cache.exp.2.sam new file mode 100644 index 000000000..1764ab18e --- /dev/null +++ b/test/cache.exp.2.sam @@ -0,0 +1,33 @@ +@HD VN:1.17 SO:coordinate +@SQ SN:T1 LN:40 +@SQ SN:T2 LN:40 +@SQ SN:T3 LN:40 +@CO @SQ SN* LN* AH AN AS DS M5 SP TP UR +@CO @RG ID* BC CN DS DT FO KS LB PG PI PL PM PU SM +@CO @PG ID* PN CL PP DS VN +@CO this is a dummy alignment file to demonstrate different abilities of hts apis +@CO QNAME FLAG RNAME POS MAPQ CIGAR RNEXT PNEXT TLEN SEQ QUAL [TAG:TYPE:VALUE]… +@CO 1234567890123456789012345678901234567890 +@CO AAAAACTGAAAACCCCTTTTGGGGACTGTTAACAGTTTTT T1 +@CO TTTTCCCCACTGAAAACCCCTTTTGGGGACTGTTAACAGT T2 +ITR1 99 T1 5 40 4M = 33 10 ACTG ()() +UNMP1 73 T1 21 40 3M * 0 5 GGG &&1 +A1 99 T1 25 35 6M = 31 8 ACTGTT ****** +B1 99 T1 25 35 6M = 31 8 GCTATT ****** +B5 355 T1 25 35 4M = 33 5 AGTG PPPP +A1 147 T1 31 33 6M = 25 -8 ACTGTT ()()() +ITR1 147 T1 33 37 4M = 5 -10 ACTG $$$$ +A4 99 T2 12 50 3M = 23 5 GAA ()( +B4 99 T2 12 50 3M = 23 5 GAT ()( +ITR2 147 T2 23 49 2M = 35 -10 TT ** +A4 147 T2 23 47 2M1X = 12 -5 TTG ((( +B4 147 T2 23 47 2M1X = 12 -5 TAG ((( +ITR2 99 T2 35 51 2M = 23 10 AA && +A5 99 T3 12 50 3M = 23 5 GAA ()( +B5 99 T3 12 50 3M = 23 5 GAT ()( +ITR3 147 T3 23 49 2M = 35 -10 TT ** +B5 147 T3 23 47 2M1X = 12 -5 TTG ((( +A5 147 T3 23 47 2M1X = 12 -5 TAG ((( +ITR3 99 T3 35 51 2M = 23 10 AA && +UNMP2 141 * 0 0 * * 0 7 AA && +UNMP3 77 * 0 0 * * 0 5 GGG &&2 diff --git a/test/cache.exp.T1.sam b/test/cache.exp.T1.sam new file mode 100644 index 000000000..96af87a0b --- /dev/null +++ b/test/cache.exp.T1.sam @@ -0,0 +1,20 @@ +@HD VN:1.17 SO:coordinate +@SQ SN:T1 LN:40 +@SQ SN:T2 LN:40 +@SQ SN:T3 LN:40 +@CO @SQ SN* LN* AH AN AS DS M5 SP TP UR +@CO @RG ID* BC CN DS DT FO KS LB PG PI PL PM PU SM +@CO @PG ID* PN CL PP DS VN +@CO this is a dummy alignment file to demonstrate different abilities of hts apis +@CO QNAME FLAG RNAME POS MAPQ CIGAR RNEXT PNEXT TLEN SEQ QUAL [TAG:TYPE:VALUE]… +@CO 1234567890123456789012345678901234567890 +@CO AAAAACTGAAAACCCCTTTTGGGGACTGTTAACAGTTTTT T1 +@CO TTTTCCCCACTGAAAACCCCTTTTGGGGACTGTTAACAGT T2 +ITR1 99 T1 5 40 4M = 33 10 ACTG ()() +UNMP1 73 T1 21 40 3M * 0 5 GGG &&1 +A1 99 T1 25 35 6M = 31 8 ACTGTT ****** +A5 355 T1 25 35 4M = 33 5 ACTG PPPP +B1 99 T1 25 35 6M = 31 8 GCTATT ****** +B5 355 T1 25 35 4M = 33 5 AGTG PPPP +A1 147 T1 31 33 6M = 25 -8 ACTGTT ()()() +ITR1 147 T1 33 37 4M = 5 -10 ACTG $$$$ diff --git a/test/cache.exp.w5.sam b/test/cache.exp.w5.sam new file mode 100644 index 000000000..f35111458 --- /dev/null +++ b/test/cache.exp.w5.sam @@ -0,0 +1,25 @@ +@HD VN:1.17 SO:coordinate +@SQ SN:T1 LN:40 +@SQ SN:T2 LN:40 +@SQ SN:T3 LN:40 +@CO @SQ SN* LN* AH AN AS DS M5 SP TP UR +@CO @RG ID* BC CN DS DT FO KS LB PG PI PL PM PU SM +@CO @PG ID* PN CL PP DS VN +@CO this is a dummy alignment file to demonstrate different abilities of hts apis +@CO QNAME FLAG RNAME POS MAPQ CIGAR RNEXT PNEXT TLEN SEQ QUAL [TAG:TYPE:VALUE]… +@CO 1234567890123456789012345678901234567890 +@CO AAAAACTGAAAACCCCTTTTGGGGACTGTTAACAGTTTTT T1 +@CO TTTTCCCCACTGAAAACCCCTTTTGGGGACTGTTAACAGT T2 +ITR1 99 T1 5 40 4M = 33 10 ACTG ()() +UNMP1 73 T1 21 40 3M * 0 5 GGG &&1 +B1 99 T1 25 35 6M = 31 8 GCTATT ****** +B5 355 T1 25 35 4M = 33 5 AGTG PPPP +A1 147 T1 31 33 6M = 25 -8 ACTGTT ()()() +B4 99 T2 12 50 3M = 23 5 GAT ()( +B4 147 T2 23 47 2M1X = 12 -5 TAG ((( +ITR2 99 T2 35 51 2M = 23 10 AA && +B5 99 T3 12 50 3M = 23 5 GAT ()( +A5 147 T3 23 47 2M1X = 12 -5 TAG ((( +ITR3 99 T3 35 51 2M = 23 10 AA && +UNMP2 141 * 0 0 * * 0 7 AA && +UNMP3 77 * 0 0 * * 0 5 GGG &&2 diff --git a/test/cache.sam b/test/cache.sam new file mode 100644 index 000000000..183c81fa5 --- /dev/null +++ b/test/cache.sam @@ -0,0 +1,34 @@ +@HD VN:1.17 SO:coordinate +@SQ SN:T1 LN:40 +@SQ SN:T2 LN:40 +@SQ SN:T3 LN:40 +@CO @SQ SN* LN* AH AN AS DS M5 SP TP UR +@CO @RG ID* BC CN DS DT FO KS LB PG PI PL PM PU SM +@CO @PG ID* PN CL PP DS VN +@CO this is a dummy alignment file to demonstrate different abilities of hts apis +@CO QNAME FLAG RNAME POS MAPQ CIGAR RNEXT PNEXT TLEN SEQ QUAL [TAG:TYPE:VALUE]… +@CO 1234567890123456789012345678901234567890 +@CO AAAAACTGAAAACCCCTTTTGGGGACTGTTAACAGTTTTT T1 +@CO TTTTCCCCACTGAAAACCCCTTTTGGGGACTGTTAACAGT T2 +ITR1 99 T1 5 40 4M = 33 10 ACTG ()() +UNMP1 73 T1 21 40 3M * 0 5 GGG &&1 +A1 99 T1 25 35 6M = 31 8 ACTGTT ****** +A5 355 T1 25 35 4M = 33 5 ACTG PPPP +B1 99 T1 25 35 6M = 31 8 GCTATT ****** +B5 355 T1 25 35 4M = 33 5 AGTG PPPP +A1 147 T1 31 33 6M = 25 -8 ACTGTT ()()() +ITR1 147 T1 33 37 4M = 5 -10 ACTG $$$$ +A4 99 T2 12 50 3M = 23 5 GAA ()( +B4 99 T2 12 50 3M = 23 5 GAT ()( +ITR2 147 T2 23 49 2M = 35 -10 TT ** +A4 147 T2 23 47 2M1X = 12 -5 TTG ((( +B4 147 T2 23 47 2M1X = 12 -5 TAG ((( +ITR2 99 T2 35 51 2M = 23 10 AA && +A5 99 T3 12 50 3M = 23 5 GAA ()( +B5 99 T3 12 50 3M = 23 5 GAT ()( +ITR3 147 T3 23 49 2M = 35 -10 TT ** +B5 147 T3 23 47 2M1X = 12 -5 TTG ((( +A5 147 T3 23 47 2M1X = 12 -5 TAG ((( +ITR3 99 T3 35 51 2M = 23 10 AA && +UNMP2 141 * 0 0 * * 0 7 AA && +UNMP3 77 * 0 0 * * 0 5 GGG &&2 diff --git a/test/test.pl b/test/test.pl index 5fd6b5f7c..503a7227b 100755 --- a/test/test.pl +++ b/test/test.pl @@ -50,6 +50,8 @@ run_test('test_MD',$opts); +run_test('test_cache',$opts); + run_test('test_vcf_api',$opts,out=>'test-vcf-api.out',needed_by=>'test_vcf_sweep'); run_test('test_bcf2vcf',$opts); run_test('test_vcf_sweep',$opts,out=>'test-vcf-sweep.out'); @@ -1050,6 +1052,57 @@ sub test_MD } } +# Tests hts lib sam cache +sub test_cache +{ + my ($opts) = @_; + $test_view_failures = 0; + print "\ntest_cache:\n"; + #regular working + testv $opts, "./test_view -p $$opts{tmp}/cache.tmp.10.sam $$opts{path}/cache.sam"; + testv $opts, "./compare_sam.pl $$opts{tmp}/cache.tmp.10.sam $$opts{path}/cache.sam"; + #no impact with window size alone + testv $opts, "./test_view -i hts_flt_wndsize=5 -p $$opts{tmp}/cache.tmp.10.sam $$opts{path}/cache.sam"; + testv $opts, "./compare_sam.pl $$opts{tmp}/cache.tmp.10.sam $$opts{path}/cache.sam"; + # high depth or all in + testv $opts, "./test_view -i hts_flt_depth=10 -p $$opts{tmp}/cache.tmp.10.sam $$opts{path}/cache.sam"; + testv $opts, "./compare_sam.pl $$opts{tmp}/cache.tmp.10.sam $$opts{path}/cache.sam"; + #smaller window than the default + testv $opts, "./test_view -i hts_flt_wndsize=350 -i hts_flt_depth=10 -p $$opts{tmp}/cache.tmp.10.sam $$opts{path}/cache.sam"; + testv $opts, "./compare_sam.pl $$opts{tmp}/cache.tmp.10.sam $$opts{path}/cache.sam"; + #too small window - different result + testv $opts, "./test_view -i hts_flt_wndsize=5 -i hts_flt_depth=1 -p $$opts{tmp}/cache.tmp.10.sam $$opts{path}/cache.sam"; + testv $opts, "./compare_sam.pl $$opts{tmp}/cache.tmp.10.sam $$opts{path}/cache.exp.w5.sam"; + # filter as much as possible + testv $opts, "./test_view -i hts_flt_depth=1 -p $$opts{tmp}/cache.tmp.1.sam $$opts{path}/cache.sam"; + testv $opts, "./compare_sam.pl $$opts{tmp}/cache.tmp.1.sam $$opts{path}/cache.exp.1.sam"; + # different depth val + testv $opts, "./test_view -i hts_flt_depth=2 -p $$opts{tmp}/cache.tmp.2.sam $$opts{path}/cache.sam"; + testv $opts, "./compare_sam.pl $$opts{tmp}/cache.tmp.2.sam $$opts{path}/cache.exp.2.sam"; + + #using iterators + testv $opts, "./test_view -b -x $$opts{tmp}/cache.tmp.bam.bai -p $$opts{tmp}/cache.tmp.bam $$opts{path}/cache.sam"; + + testv $opts, "./test_view -i hts_flt_depth=10 -p $$opts{tmp}/cache.tmp.10.T1.sam $$opts{tmp}/cache.tmp.bam T1"; + testv $opts, "./compare_sam.pl $$opts{tmp}/cache.tmp.10.T1.sam $$opts{path}/cache.exp.T1.sam"; + # filter as much as possible + testv $opts, "./test_view -i hts_flt_depth=1 -p $$opts{tmp}/cache.tmp.1.T1T3T2.sam $$opts{tmp}/cache.tmp.bam T1 T3 T2"; + testv $opts, "./compare_sam.pl $$opts{tmp}/cache.tmp.1.T1T3T2.sam $$opts{path}/cache.exp.1.T1T3T2.sam"; + # different depth val + testv $opts, "./test_view -i hts_flt_depth=2 -p $$opts{tmp}/cache.tmp.2.sam $$opts{tmp}/cache.tmp.bam T1 T2 T3"; + testv $opts, "head -n 31 $$opts{path}/cache.exp.2.sam > $$opts{tmp}/cache.tmp.exp.sam"; + testv $opts, "./compare_sam.pl $$opts{tmp}/cache.tmp.2.sam $$opts{tmp}/cache.tmp.exp.sam"; + # multi region iterator + testv $opts, "./test_view -M -i hts_flt_depth=2 -p $$opts{tmp}/cache.tmp.2.T3T1T2.sam $$opts{tmp}/cache.tmp.bam T3 T1 T2"; + testv $opts, "./compare_sam.pl $$opts{tmp}/cache.tmp.2.T3T1T2.sam $$opts{tmp}/cache.tmp.exp.sam"; + + if ($test_view_failures == 0) { + passed($opts, "cache tests"); + } else { + failed($opts, "cache tests", "$test_view_failures subtests failed"); + } +} + sub test_index { my ($opts, $nthreads) = @_;