From f911f7f0f2d2d8c4b5c8d17ea72608b5fa9a4d94 Mon Sep 17 00:00:00 2001 From: mbstadler Date: Tue, 22 Sep 2026 14:29:04 +0200 Subject: [PATCH 01/13] initial draft for getBaseCoverageForBam --- R/RcppExports.R | 32 ++++++++++ src/RcppExports.cpp | 15 +++++ src/utils.cpp | 142 ++++++++++++++++++++++++++++++++++++++++++++ src/utils.h | 13 ++++ 4 files changed, 202 insertions(+) diff --git a/R/RcppExports.R b/R/RcppExports.R index 82a56fd..56d7fef 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -711,6 +711,11 @@ NULL #' @keywords internal NULL +#' Helper function for getBaseCoverageForBam +#' @noRd +#' @keywords internal +NULL + #' Concatenate files #' #' @param input_files Character vector with input file names to concatenate. @@ -781,3 +786,30 @@ complement <- function(n) { .Call(`_SingleMoleculeGenomicsIO_complement`, n) } +#' Get base coverage histogram for a BAM file +#' +#' @param bamfile A character scalar with the bam file name (and path). +#' @param regions Character vector specifying the region(s) for which +#' to calculate coverage, in the form \code{"."}, \code{"chr"} or +#' \code{"chr:start-end"}. +#' @param maxDepth An integer scalar defining the maximal depth to consider. +#' @param nThreads A numeric scalar with the number of threads used for +#' decompressing BAM records. +#' +#' @details CIGAR operations are considered, thus not counting the genomic +#' bases in a read-insertion as covered. Secondary and supplementary +#' alignments are not included. +#' +#' @reference The algorithm was described in Pedersen BS and Quinlan AR. +#' "Mosdepth: quick coverage calculation for genomes and exomes". +#' Bioinformatics. 2018; 34(5):867-868. doi: 10.1093/bioinformatics/btx699 +#' +#' @return A numeric vector of length \code{maxDepth + 1}, with values at +#' index \code{i} giving the number of positions that were overlapped by +#' \code{maxDepth} are also added to the value for \code{maxDepth} +#' at index \code{maxDepth + 1}. +#' +getBaseCoverageForBam <- function(bamfile, regions = NULL, maxDepth = 200L, nThreads = 3L) { + .Call(`_SingleMoleculeGenomicsIO_getBaseCoverageForBam`, bamfile, regions, maxDepth, nThreads) +} + diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 8ad93c6..b904184 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -247,6 +247,20 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// getBaseCoverageForBam +Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions, const uint maxDepth, int nThreads); +RcppExport SEXP _SingleMoleculeGenomicsIO_getBaseCoverageForBam(SEXP bamfileSEXP, SEXP regionsSEXP, SEXP maxDepthSEXP, SEXP nThreadsSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const std::string >::type bamfile(bamfileSEXP); + Rcpp::traits::input_parameter< Rcpp::Nullable> >::type regions(regionsSEXP); + Rcpp::traits::input_parameter< const uint >::type maxDepth(maxDepthSEXP); + Rcpp::traits::input_parameter< int >::type nThreads(nThreadsSEXP); + rcpp_result_gen = Rcpp::wrap(getBaseCoverageForBam(bamfile, regions, maxDepth, nThreads)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_SingleMoleculeGenomicsIO_estimateNoise", (DL_FUNC) &_SingleMoleculeGenomicsIO_estimateNoise, 4}, @@ -263,6 +277,7 @@ static const R_CallMethodDef CallEntries[] = { {"_SingleMoleculeGenomicsIO_getChromosomeNamesFromBam", (DL_FUNC) &_SingleMoleculeGenomicsIO_getChromosomeNamesFromBam, 1}, {"_SingleMoleculeGenomicsIO_get_unmodified_base", (DL_FUNC) &_SingleMoleculeGenomicsIO_get_unmodified_base, 1}, {"_SingleMoleculeGenomicsIO_complement", (DL_FUNC) &_SingleMoleculeGenomicsIO_complement, 1}, + {"_SingleMoleculeGenomicsIO_getBaseCoverageForBam", (DL_FUNC) &_SingleMoleculeGenomicsIO_getBaseCoverageForBam, 4}, {NULL, NULL, 0} }; diff --git a/src/utils.cpp b/src/utils.cpp index 937ba99..9eb336b 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -264,6 +264,7 @@ Rcpp::CharacterVector getChromosomeNamesFromBam(const std::string bamfile) { } } + //' Get target and text lines from loaded BAM header //' //' @param inbamhdr A loaded (populated) sam_hdr_t*. @@ -1091,3 +1092,144 @@ int check_bam_format(samFile *infile, } return result; } + + +//' Helper function for getBaseCoverageForBam +//' @noRd +//' @keywords internal +static void region_cov(samFile *fp, hts_itr_t *it, bam1_t *b, + hts_pos_t beg, hts_pos_t end, int32_t *diff, + uint32_t *hist, uint maxd) { + hts_pos_t len = end - beg; + memset(diff, 0, (len + 1) * sizeof(int32_t)); + + while (sam_itr_next(fp, it, b) >= 0) { + if (b->core.flag & (BAM_FUNMAP | BAM_FSECONDARY | BAM_FSUPPLEMENTARY)) continue; + hts_pos_t pos = b->core.pos; + const uint32_t *cig = bam_get_cigar(b); + for (uint32_t i = 0; i < b->core.n_cigar; i++) { + int op = bam_cigar_op(cig[i]); + hts_pos_t l = bam_cigar_oplen(cig[i]); + if (op == BAM_CMATCH || op == BAM_CEQUAL || op == BAM_CDIFF) { + hts_pos_t s = pos - beg, e = s + l; + if (e > 0 && s < len) { // clip to region + diff[s < 0 ? 0 : s]++; + diff[e > len ? len : e]--; + } + pos += l; + } else if (op == BAM_CDEL || op == BAM_CREF_SKIP) { + pos += l; // gap: no coverage + } + } + } + int32_t d = 0; + for (hts_pos_t i = 0; i < len; i++) { + d += diff[i]; + hist[d < maxd ? d : maxd]++; + } +} + +//' Get base coverage histogram for a BAM file +//' +//' @param bamfile A character scalar with the bam file name (and path). +//' @param regions Character vector specifying the region(s) for which +//' to calculate coverage, in the form \code{"."}, \code{"chr"} or +//' \code{"chr:start-end"}. +//' @param maxDepth An integer scalar defining the maximal depth to consider. +//' @param nThreads A numeric scalar with the number of threads used for +//' decompressing BAM records. +//' +//' @details CIGAR operations are considered, thus not counting the genomic +//' bases in a read-insertion as covered. Secondary and supplementary +//' alignments are not included. +//' +//' @reference The algorithm was described in Pedersen BS and Quinlan AR. +//' "Mosdepth: quick coverage calculation for genomes and exomes". +//' Bioinformatics. 2018; 34(5):867-868. doi: 10.1093/bioinformatics/btx699 +//' +//' @return A numeric vector of length \code{maxDepth + 1}, with values at +//' index \code{i} giving the number of positions that were overlapped by +//' \code{maxDepth} are also added to the value for \code{maxDepth} +//' at index \code{maxDepth + 1}. +//' +// [[Rcpp::export]] +Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, + Rcpp::Nullable> regions = R_NilValue, + const uint maxDepth = 200, + int nThreads = 3) { + + int buffer_len = 2000, success = 0; + char buffer[2000]; + bool had_error = false; + std::vector regionsvect; + samFile *inbamfile = NULL; + const char* inname = bamfile.c_str(); + bam1_t *bamdata = NULL; + hts_idx_t *idx = NULL; + hts_itr_t *iter = NULL; + sam_hdr_t *inbamhdr = NULL; + int32_t *diff = NULL; + uint32_t *hist = (uint32_t*)calloc(maxDepth + 1, sizeof(uint32_t)); // hist[maxDepth] = overflow bin + Rcpp::NumericVector histvect(maxDepth + 1); + hts_pos_t maxlen = 0; + + // set default regions if NULL + if (regions.isNotNull()) { + regionsvect = Rcpp::as>(regions); + } else { + regionsvect = {"."}; // default: whole genome + } + + // turn htslib logging off -> handle via Rcpp::warning or Rcpp::stop + hts_set_log_level(HTS_LOG_OFF); + + // open bam file and read index and header + success = open_bam_and_read_index_and_header(bamdata, inname, + inbamfile, idx, + inbamhdr, nThreads, + had_error, buffer_len, buffer); + if (success != 0) { + goto end; + } + + // create multi-region iterator + success = create_multi_region_iterator(regionsvect, iter, idx, inbamhdr, + had_error, buffer_len, buffer); + if (success != 0) { + goto end; + } + + // get length of longest chromosome and allocate per-base start/end vector `diff` + for (int t = 0; t < sam_hdr_nref(inbamhdr); t++) + if (sam_hdr_tid2len(inbamhdr, t) > maxlen) maxlen = sam_hdr_tid2len(inbamhdr, t); + diff = (int32_t*)malloc((maxlen + 1) * sizeof(int32_t)); + if (diff == NULL) { + had_error = true; + snprintf(buffer, buffer_len, "Failed to allocate memory for coverage vector\n"); + goto end; + } + + // populate `diff` and `hist` + region_cov(inbamfile, iter, bamdata, 0, maxlen, diff, hist, maxDepth); + for (int i = 0; i <= maxDepth; i++) { + histvect[i] = hist[i]; + } + + end: + //cleanup + if (diff) free(diff); + if (hist) free(hist); + if (bamdata) bam_destroy1(bamdata); + if (iter) sam_itr_destroy(iter); + if (idx) hts_idx_destroy(idx); + if (inbamhdr) sam_hdr_destroy(inbamhdr); + if (inbamfile) sam_close(inbamfile); + if (hist) free(hist); + + if (had_error) { + // we encountered an error (message in `buffer`) --> stop + Rcpp::stop(buffer); + } + + return histvect; +} diff --git a/src/utils.h b/src/utils.h index 121a169..b649cf9 100644 --- a/src/utils.h +++ b/src/utils.h @@ -87,3 +87,16 @@ int check_bam_format(samFile *infile, bool &had_error, char *buffer, int &buffer_len); + +static void region_cov(samFile*, + hts_itr_t*, + bam1_t*, + hts_pos_t, + hts_pos_t, + int32_t*, + uint32_t*, + uint); +Rcpp::NumericVector getBaseCoverageForBam(const std::string, + std::vector, + const uint, + int); From a4bce43462fe3c46677dca03a5c9e78016217e4c Mon Sep 17 00:00:00 2001 From: mbstadler Date: Tue, 22 Sep 2026 14:39:01 +0200 Subject: [PATCH 02/13] update docs --- R/RcppExports.R | 8 ++++--- man/getBaseCoverageForBam.Rd | 41 ++++++++++++++++++++++++++++++++++++ src/utils.cpp | 8 ++++--- 3 files changed, 51 insertions(+), 6 deletions(-) create mode 100644 man/getBaseCoverageForBam.Rd diff --git a/R/RcppExports.R b/R/RcppExports.R index 56d7fef..f1b826b 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -800,13 +800,15 @@ complement <- function(n) { #' bases in a read-insertion as covered. Secondary and supplementary #' alignments are not included. #' -#' @reference The algorithm was described in Pedersen BS and Quinlan AR. +#' @references The algorithm was described in Pedersen BS and Quinlan AR. #' "Mosdepth: quick coverage calculation for genomes and exomes". -#' Bioinformatics. 2018; 34(5):867-868. doi: 10.1093/bioinformatics/btx699 +#' Bioinformatics. 2018; 34(5):867-868. +#' \url{https://doi.org/10.1093/bioinformatics/btx699} #' #' @return A numeric vector of length \code{maxDepth + 1}, with values at #' index \code{i} giving the number of positions that were overlapped by -#' \code{maxDepth} are also added to the value for \code{maxDepth} +#' \code{i - 1} alignments. Positions overlapped by more than +#' \code{maxDepth} alignments are also added to the value for \code{maxDepth} #' at index \code{maxDepth + 1}. #' getBaseCoverageForBam <- function(bamfile, regions = NULL, maxDepth = 200L, nThreads = 3L) { diff --git a/man/getBaseCoverageForBam.Rd b/man/getBaseCoverageForBam.Rd new file mode 100644 index 0000000..aefa521 --- /dev/null +++ b/man/getBaseCoverageForBam.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{getBaseCoverageForBam} +\alias{getBaseCoverageForBam} +\title{Get base coverage histogram for a BAM file} +\usage{ +getBaseCoverageForBam(bamfile, regions = NULL, maxDepth = 200L, nThreads = 3L) +} +\arguments{ +\item{bamfile}{A character scalar with the bam file name (and path).} + +\item{regions}{Character vector specifying the region(s) for which +to calculate coverage, in the form \code{"."}, \code{"chr"} or +\code{"chr:start-end"}.} + +\item{maxDepth}{An integer scalar defining the maximal depth to consider.} + +\item{nThreads}{A numeric scalar with the number of threads used for +decompressing BAM records.} +} +\value{ +A numeric vector of length \code{maxDepth + 1}, with values at +index \code{i} giving the number of positions that were overlapped by +\code{i - 1} alignments. Positions overlapped by more than +\code{maxDepth} alignments are also added to the value for \code{maxDepth} +at index \code{maxDepth + 1}. +} +\description{ +Get base coverage histogram for a BAM file +} +\details{ +CIGAR operations are considered, thus not counting the genomic +bases in a read-insertion as covered. Secondary and supplementary +alignments are not included. +} +\references{ +The algorithm was described in Pedersen BS and Quinlan AR. +"Mosdepth: quick coverage calculation for genomes and exomes". +Bioinformatics. 2018; 34(5):867-868. +\url{https://doi.org/10.1093/bioinformatics/btx699} +} diff --git a/src/utils.cpp b/src/utils.cpp index 9eb336b..5021362 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1143,13 +1143,15 @@ static void region_cov(samFile *fp, hts_itr_t *it, bam1_t *b, //' bases in a read-insertion as covered. Secondary and supplementary //' alignments are not included. //' -//' @reference The algorithm was described in Pedersen BS and Quinlan AR. +//' @references The algorithm was described in Pedersen BS and Quinlan AR. //' "Mosdepth: quick coverage calculation for genomes and exomes". -//' Bioinformatics. 2018; 34(5):867-868. doi: 10.1093/bioinformatics/btx699 +//' Bioinformatics. 2018; 34(5):867-868. +//' \url{https://doi.org/10.1093/bioinformatics/btx699} //' //' @return A numeric vector of length \code{maxDepth + 1}, with values at //' index \code{i} giving the number of positions that were overlapped by -//' \code{maxDepth} are also added to the value for \code{maxDepth} +//' \code{i - 1} alignments. Positions overlapped by more than +//' \code{maxDepth} alignments are also added to the value for \code{maxDepth} //' at index \code{maxDepth + 1}. //' // [[Rcpp::export]] From a3f780c8607a00c7a9d556eff72ce7652b67a231 Mon Sep 17 00:00:00 2001 From: mbstadler Date: Tue, 22 Sep 2026 14:44:50 +0200 Subject: [PATCH 03/13] return IntegerVector from getBaseCoverageForBam --- src/RcppExports.cpp | 2 +- src/utils.cpp | 4 ++-- src/utils.h | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index b904184..747d2aa 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -248,7 +248,7 @@ BEGIN_RCPP END_RCPP } // getBaseCoverageForBam -Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions, const uint maxDepth, int nThreads); +Rcpp::IntegerVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions, const uint maxDepth, int nThreads); RcppExport SEXP _SingleMoleculeGenomicsIO_getBaseCoverageForBam(SEXP bamfileSEXP, SEXP regionsSEXP, SEXP maxDepthSEXP, SEXP nThreadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; diff --git a/src/utils.cpp b/src/utils.cpp index 5021362..9042a68 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1155,7 +1155,7 @@ static void region_cov(samFile *fp, hts_itr_t *it, bam1_t *b, //' at index \code{maxDepth + 1}. //' // [[Rcpp::export]] -Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, +Rcpp::IntegerVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions = R_NilValue, const uint maxDepth = 200, int nThreads = 3) { @@ -1172,7 +1172,7 @@ Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, sam_hdr_t *inbamhdr = NULL; int32_t *diff = NULL; uint32_t *hist = (uint32_t*)calloc(maxDepth + 1, sizeof(uint32_t)); // hist[maxDepth] = overflow bin - Rcpp::NumericVector histvect(maxDepth + 1); + Rcpp::IntegerVector histvect(maxDepth + 1); hts_pos_t maxlen = 0; // set default regions if NULL diff --git a/src/utils.h b/src/utils.h index b649cf9..a954834 100644 --- a/src/utils.h +++ b/src/utils.h @@ -96,7 +96,7 @@ static void region_cov(samFile*, int32_t*, uint32_t*, uint); -Rcpp::NumericVector getBaseCoverageForBam(const std::string, +Rcpp::IntegerVector getBaseCoverageForBam(const std::string, std::vector, const uint, int); From c3bc9801d35e0f9afd8cd74bcf8bc4ec17a8d526 Mon Sep 17 00:00:00 2001 From: mbstadler Date: Tue, 22 Sep 2026 17:52:58 +0200 Subject: [PATCH 04/13] only count zeros inside of regions Co-authored-by: Charlotte Soneson Co-authored-by: Michael Stadler --- R/RcppExports.R | 5 +- man/getBaseCoverageForBam.Rd | 5 +- src/RcppExports.cpp | 2 +- src/utils.cpp | 164 +++++++++++++++++++++++++++-------- src/utils.h | 6 +- 5 files changed, 140 insertions(+), 42 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index f1b826b..f0f1e21 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -791,7 +791,8 @@ complement <- function(n) { #' @param bamfile A character scalar with the bam file name (and path). #' @param regions Character vector specifying the region(s) for which #' to calculate coverage, in the form \code{"."}, \code{"chr"} or -#' \code{"chr:start-end"}. +#' \code{"chr:start-end"}. If \code{NULL}, the whole genome +#' (\code{"."}) is used by default. #' @param maxDepth An integer scalar defining the maximal depth to consider. #' @param nThreads A numeric scalar with the number of threads used for #' decompressing BAM records. @@ -807,7 +808,7 @@ complement <- function(n) { #' #' @return A numeric vector of length \code{maxDepth + 1}, with values at #' index \code{i} giving the number of positions that were overlapped by -#' \code{i - 1} alignments. Positions overlapped by more than +#' exactly \code{i} alignments. Positions overlapped by more than #' \code{maxDepth} alignments are also added to the value for \code{maxDepth} #' at index \code{maxDepth + 1}. #' diff --git a/man/getBaseCoverageForBam.Rd b/man/getBaseCoverageForBam.Rd index aefa521..a8a3f66 100644 --- a/man/getBaseCoverageForBam.Rd +++ b/man/getBaseCoverageForBam.Rd @@ -11,7 +11,8 @@ getBaseCoverageForBam(bamfile, regions = NULL, maxDepth = 200L, nThreads = 3L) \item{regions}{Character vector specifying the region(s) for which to calculate coverage, in the form \code{"."}, \code{"chr"} or -\code{"chr:start-end"}.} +\code{"chr:start-end"}. If \code{NULL}, the whole genome +(\code{"."}) is used by default.} \item{maxDepth}{An integer scalar defining the maximal depth to consider.} @@ -21,7 +22,7 @@ decompressing BAM records.} \value{ A numeric vector of length \code{maxDepth + 1}, with values at index \code{i} giving the number of positions that were overlapped by -\code{i - 1} alignments. Positions overlapped by more than +exactly \code{i} alignments. Positions overlapped by more than \code{maxDepth} alignments are also added to the value for \code{maxDepth} at index \code{maxDepth + 1}. } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 747d2aa..b904184 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -248,7 +248,7 @@ BEGIN_RCPP END_RCPP } // getBaseCoverageForBam -Rcpp::IntegerVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions, const uint maxDepth, int nThreads); +Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions, const uint maxDepth, int nThreads); RcppExport SEXP _SingleMoleculeGenomicsIO_getBaseCoverageForBam(SEXP bamfileSEXP, SEXP regionsSEXP, SEXP maxDepthSEXP, SEXP nThreadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; diff --git a/src/utils.cpp b/src/utils.cpp index 9042a68..d31f7c5 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include "utils.h" @@ -1099,27 +1100,34 @@ int check_bam_format(samFile *infile, //' @keywords internal static void region_cov(samFile *fp, hts_itr_t *it, bam1_t *b, hts_pos_t beg, hts_pos_t end, int32_t *diff, - uint32_t *hist, uint maxd) { + uint64_t *hist, uint maxd) { hts_pos_t len = end - beg; memset(diff, 0, (len + 1) * sizeof(int32_t)); while (sam_itr_next(fp, it, b) >= 0) { - if (b->core.flag & (BAM_FUNMAP | BAM_FSECONDARY | BAM_FSUPPLEMENTARY)) continue; + if (b->core.flag & (BAM_FUNMAP | BAM_FSECONDARY | BAM_FSUPPLEMENTARY)) { + continue; + } hts_pos_t pos = b->core.pos; const uint32_t *cig = bam_get_cigar(b); + for (uint32_t i = 0; i < b->core.n_cigar; i++) { + uint32_t opl = bam_cigar_oplen(cig[i]); int op = bam_cigar_op(cig[i]); - hts_pos_t l = bam_cigar_oplen(cig[i]); if (op == BAM_CMATCH || op == BAM_CEQUAL || op == BAM_CDIFF) { - hts_pos_t s = pos - beg, e = s + l; - if (e > 0 && s < len) { // clip to region - diff[s < 0 ? 0 : s]++; - diff[e > len ? len : e]--; + // soft-clips before this op consumed no reference bases, + // so the reference position of this op's start is still `pos` + hts_pos_t s_rel = pos - beg; + hts_pos_t e_rel = s_rel + opl; + if (e_rel > 0 && s_rel < len) { + diff[s_rel < 0 ? 0 : s_rel]++; + diff[e_rel > len ? len : e_rel]--; } - pos += l; + pos += opl; } else if (op == BAM_CDEL || op == BAM_CREF_SKIP) { - pos += l; // gap: no coverage + pos += opl; } + // S / I / P do not advance the reference position } } int32_t d = 0; @@ -1129,12 +1137,78 @@ static void region_cov(samFile *fp, hts_itr_t *it, bam1_t *b, } } +// Resolve a vector of region strings into (tid, beg, end) spans, where +// `beg` is 0-based and `end` is exclusive. "." or "*" expands to all +// non-empty contigs. Returns 0 on success, -1 on failure (message in buffer). +static int resolve_regions_to_spans(const std::vector ®ionsvect, + sam_hdr_t *hdr, + std::vector> &spans, + bool &had_error, int buffer_len, char *buffer) { + int n = sam_hdr_nref(hdr); + for (const std::string &rs : regionsvect) { + if (rs == "." || rs == "*") { + spans.clear(); + for (int t = 0; t < n; t++) { + hts_pos_t L = sam_hdr_tid2len(hdr, t); + if (L > 0) spans.push_back(std::make_tuple(t, (hts_pos_t)0, L)); + } + break; + } + size_t colon = rs.find(':'); + if (colon == std::string::npos) { + int tid = sam_hdr_name2tid(hdr, rs.c_str()); + if (tid < 0) { + had_error = true; + snprintf(buffer, buffer_len, "Unknown contig: %s\n", rs.c_str()); + return -1; + } + hts_pos_t L = sam_hdr_tid2len(hdr, tid); + if (L > 0) spans.push_back(std::make_tuple(tid, (hts_pos_t)0, L)); + continue; + } + std::string name = rs.substr(0, colon); + std::string posstr = rs.substr(colon + 1); + size_t dash = posstr.find('-'); + if (dash == std::string::npos) { + had_error = true; + snprintf(buffer, buffer_len, + "Malformed region (expected 'name:beg-end'): %s\n", rs.c_str()); + return -1; + } + errno = 0; + long long beg_in = strtoll(posstr.substr(0, dash).c_str(), NULL, 10); + errno = 0; + long long end_in = strtoll(posstr.substr(dash + 1).c_str(), NULL, 10); + if (errno != 0 || beg_in < 1 || end_in < beg_in) { + had_error = true; + snprintf(buffer, buffer_len, + "Invalid region position in: %s\n", rs.c_str()); + return -1; + } + int tid = sam_hdr_name2tid(hdr, name.c_str()); + if (tid < 0) { + had_error = true; + snprintf(buffer, buffer_len, "Unknown contig: %s\n", name.c_str()); + return -1; + } + hts_pos_t hdr_len = sam_hdr_tid2len(hdr, tid); + hts_pos_t beg = (hts_pos_t)beg_in - 1; // 0-based + hts_pos_t end = (hts_pos_t)end_in > hdr_len ? hdr_len : (hts_pos_t)end_in; // exclusive, clamped + if (end > beg) { + spans.push_back(std::make_tuple(tid, beg, end)); + } + } + return 0; +} + + //' Get base coverage histogram for a BAM file //' //' @param bamfile A character scalar with the bam file name (and path). //' @param regions Character vector specifying the region(s) for which //' to calculate coverage, in the form \code{"."}, \code{"chr"} or -//' \code{"chr:start-end"}. +//' \code{"chr:start-end"}. If \code{NULL}, the whole genome +//' (\code{"."}) is used by default. //' @param maxDepth An integer scalar defining the maximal depth to consider. //' @param nThreads A numeric scalar with the number of threads used for //' decompressing BAM records. @@ -1150,12 +1224,12 @@ static void region_cov(samFile *fp, hts_itr_t *it, bam1_t *b, //' //' @return A numeric vector of length \code{maxDepth + 1}, with values at //' index \code{i} giving the number of positions that were overlapped by -//' \code{i - 1} alignments. Positions overlapped by more than +//' exactly \code{i} alignments. Positions overlapped by more than //' \code{maxDepth} alignments are also added to the value for \code{maxDepth} //' at index \code{maxDepth + 1}. //' // [[Rcpp::export]] -Rcpp::IntegerVector getBaseCoverageForBam(const std::string bamfile, +Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions = R_NilValue, const uint maxDepth = 200, int nThreads = 3) { @@ -1164,16 +1238,15 @@ Rcpp::IntegerVector getBaseCoverageForBam(const std::string bamfile, char buffer[2000]; bool had_error = false; std::vector regionsvect; + std::vector> spans; samFile *inbamfile = NULL; const char* inname = bamfile.c_str(); bam1_t *bamdata = NULL; hts_idx_t *idx = NULL; hts_itr_t *iter = NULL; sam_hdr_t *inbamhdr = NULL; - int32_t *diff = NULL; - uint32_t *hist = (uint32_t*)calloc(maxDepth + 1, sizeof(uint32_t)); // hist[maxDepth] = overflow bin - Rcpp::IntegerVector histvect(maxDepth + 1); - hts_pos_t maxlen = 0; + uint64_t *hist = (uint64_t*)calloc(maxDepth + 1, sizeof(uint64_t)); // hist[maxDepth] = overflow bin + Rcpp::NumericVector histvect(maxDepth + 1); // set default regions if NULL if (regions.isNotNull()) { @@ -1194,39 +1267,62 @@ Rcpp::IntegerVector getBaseCoverageForBam(const std::string bamfile, goto end; } - // create multi-region iterator - success = create_multi_region_iterator(regionsvect, iter, idx, inbamhdr, - had_error, buffer_len, buffer); + // resolve regions to (tid, beg, end) spans in header coordinates + success = resolve_regions_to_spans(regionsvect, inbamhdr, spans, + had_error, buffer_len, buffer); if (success != 0) { goto end; } - // get length of longest chromosome and allocate per-base start/end vector `diff` - for (int t = 0; t < sam_hdr_nref(inbamhdr); t++) - if (sam_hdr_tid2len(inbamhdr, t) > maxlen) maxlen = sam_hdr_tid2len(inbamhdr, t); - diff = (int32_t*)malloc((maxlen + 1) * sizeof(int32_t)); - if (diff == NULL) { - had_error = true; - snprintf(buffer, buffer_len, "Failed to allocate memory for coverage vector\n"); - goto end; + // process each span independently, accumulating into `hist` + for (const std::tuple &sp : spans) { + int tid = std::get<0>(sp); + hts_pos_t beg = std::get<1>(sp); + hts_pos_t end = std::get<2>(sp); + hts_pos_t len = end - beg; + if (len == 0) { + continue; + } + + int32_t *diff = (int32_t*)malloc((len + 1) * sizeof(int32_t)); + if (diff == NULL) { + had_error = true; + snprintf(buffer, buffer_len, + "Failed to allocate memory for coverage vector\n"); + goto end; + } + + iter = sam_itr_queryi(idx, tid, beg, end); + if (iter == NULL) { + free(diff); + diff = NULL; + had_error = true; + snprintf(buffer, buffer_len, "Failed to create region iterator\n"); + goto end; + } + + region_cov(inbamfile, iter, bamdata, beg, end, diff, hist, maxDepth); + + if (iter) { + sam_itr_destroy(iter); + } + iter = NULL; + free(diff); + diff = NULL; } - // populate `diff` and `hist` - region_cov(inbamfile, iter, bamdata, 0, maxlen, diff, hist, maxDepth); + // copy result into return vector before releasing `hist` for (int i = 0; i <= maxDepth; i++) { - histvect[i] = hist[i]; + histvect[i] = (double) hist[i]; } end: - //cleanup - if (diff) free(diff); + // cleanup if (hist) free(hist); if (bamdata) bam_destroy1(bamdata); - if (iter) sam_itr_destroy(iter); if (idx) hts_idx_destroy(idx); if (inbamhdr) sam_hdr_destroy(inbamhdr); if (inbamfile) sam_close(inbamfile); - if (hist) free(hist); if (had_error) { // we encountered an error (message in `buffer`) --> stop diff --git a/src/utils.h b/src/utils.h index a954834..6f0c532 100644 --- a/src/utils.h +++ b/src/utils.h @@ -94,9 +94,9 @@ static void region_cov(samFile*, hts_pos_t, hts_pos_t, int32_t*, - uint32_t*, + uint64_t*, uint); -Rcpp::IntegerVector getBaseCoverageForBam(const std::string, - std::vector, +Rcpp::NumericVector getBaseCoverageForBam(const std::string, + Rcpp::Nullable>, const uint, int); From 4b23565bb64355dfe4923f9192a369c63481dff9 Mon Sep 17 00:00:00 2001 From: mbstadler Date: Wed, 23 Sep 2026 15:16:41 +0200 Subject: [PATCH 05/13] add unit test for getBaseCoverageForBam Co-authored-by: Charlotte Soneson Co-authored-by: Michael Stadler --- tests/testthat/test-utils.R | 82 +++++++++++++++++++++++++++++++++++++ 1 file changed, 82 insertions(+) diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R index 1fc8f0a..2b086cf 100644 --- a/tests/testthat/test-utils.R +++ b/tests/testthat/test-utils.R @@ -374,3 +374,85 @@ test_that("refargToDNAStringSet works", { unlink(rlibdir) }) + +## -------------------------------------------------------------------------- ## +## Checks, getBaseCoverageForBam +## -------------------------------------------------------------------------- ## +test_that("getBaseCoverageForBam works", { + bamfile <- system.file("extdata", "6mA_1_10reads.bam", + package = "SingleMoleculeGenomicsIO") + expect_true(file.exists(bamfile)) + + ## -- error cases -------------------------------------------------------- ## + expect_error(getBaseCoverageForBam(bamfile = "NO_SUCH_FILE.bam"), + "Could not open input file") + expect_error(getBaseCoverageForBam(bamfile, regions = "notarealcontig"), + "Unknown contig") + expect_error(getBaseCoverageForBam(bamfile, regions = "chr1:6940000"), + "Malformed region") + expect_error(getBaseCoverageForBam(bamfile, regions = "chr1:5-3"), + "Invalid region position") + expect_error(getBaseCoverageForBam(bamfile, regions = "chr1:-10"), + "Invalid region position") + + ## -- whole genome (NULL == ".") ------------------------------------------ ## + vNULL <- getBaseCoverageForBam(bamfile, regions = NULL) + vDot <- getBaseCoverageForBam(bamfile, regions = ".") + expect_identical(vNULL, vDot) + + # structure + expect_length(vNULL, 201L) + expect_true(all(vNULL >= 0)) + expect_false(any(is.na(vNULL))) + + # ground truth (verified against samtools depth -a) + expect_identical(head(vNULL, 11), + c(2728206677, 3841, 811, 563, 973, 596, 675, 463, 420, 1066, 6366)) + expect_true(all(vNULL[12:201] == 0)) # max depth is 10 for 10 reads + + ## -- single contig (no coordinates) -------------------------------------- ## + # this file only has reads on chr1, so the per-depth histogram matches + # the whole genome, but with zero-count restricted to chr1 length. + vChr1 <- getBaseCoverageForBam(bamfile, regions = "chr1") + expect_identical(head(vChr1, 11), + c(195138505, 3841, 811, 563, 973, 596, 675, 463, 420, 1066, 6366)) + expect_true(all(vChr1[12:201] == 0)) + expect_equal(sum(vChr1), 195154279) # contig length + + ## -- bounded region ------------------------------------------------------ ## + vRegion <- getBaseCoverageForBam(bamfile, regions = "chr1:6940000-6960000") + expect_identical(head(vRegion, 5), + c(18370, 841, 651, 139, 0)) + expect_true(all(vRegion[6:201] == 0)) + expect_equal(sum(vRegion), 20001) # total positions in region + + ## -- region end beyond contig end is clamped ------------------------------ ## + vClamp <- getBaseCoverageForBam( + bamfile, regions = "chr1:6940000-999999999") + expect_identical(head(vClamp, 4), + c(188212649, 841, 651, 139)) + # clamped to chr1 length, so sum = 195154279 - 6939999 + expect_equal(sum(vClamp), 195154279 - 6939999) + + ## -- multi-region (disjoint union) ---------------------------------------- ## + vA <- getBaseCoverageForBam(bamfile, regions = "chr1:6940000-6950000") + vB <- getBaseCoverageForBam(bamfile, regions = "chr1:6950001-6960000") + vMulti <- getBaseCoverageForBam( + bamfile, + regions = c("chr1:6940000-6950000", "chr1:6950001-6960000")) + # sum of disjoint regions == single region + expect_identical(vA + vB, vMulti) + expect_identical(vMulti, vRegion) + + ## -- overflow bin (maxDepth < actual max coverage) ------------------------ ## + vSmall <- getBaseCoverageForBam( + bamfile, regions = "chr1:6940000-6960000", maxDepth = 2L) + expect_length(vSmall, 3L) + # index 2 = overflow: positions covered by depth 3 AND 4 + expect_identical(vSmall, c(18370, 841, 651 + 139)) + + ## -- nThreads parameter does not affect result ---------------------------- ## + v2t <- getBaseCoverageForBam( + bamfile, regions = "chr1:6940000-6960000", nThreads = 2L) + expect_identical(v2t, vRegion) +}) From 1005aae1c2cbad62e987edc8a50cbaf38ffbb039 Mon Sep 17 00:00:00 2001 From: mbstadler Date: Thu, 24 Sep 2026 08:15:05 +0200 Subject: [PATCH 06/13] use sam_parse_region (htslib/sam.h) instead of own parser Co-authored-by: Charlotte Soneson Co-authored-by: Michael Stadler --- R/RcppExports.R | 7 +++-- man/getBaseCoverageForBam.Rd | 7 +++-- src/utils.cpp | 61 +++++++++++++----------------------- tests/testthat/test-utils.R | 18 +++++++---- 4 files changed, 44 insertions(+), 49 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index f0f1e21..556baeb 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -790,8 +790,11 @@ complement <- function(n) { #' #' @param bamfile A character scalar with the bam file name (and path). #' @param regions Character vector specifying the region(s) for which -#' to calculate coverage, in the form \code{"."}, \code{"chr"} or -#' \code{"chr:start-end"}. If \code{NULL}, the whole genome +#' to calculate coverage. Each region uses the grammar supported by the +#' \code{htslib} function \code{sam_parse_region}, for example +#' \code{"."} (all contigs), \code{"chr"} (whole contig), +#' \code{"chr:START"}, \code{"chr:-END"} or +#' \code{"chr:START-END"}. If \code{NULL}, the whole genome #' (\code{"."}) is used by default. #' @param maxDepth An integer scalar defining the maximal depth to consider. #' @param nThreads A numeric scalar with the number of threads used for diff --git a/man/getBaseCoverageForBam.Rd b/man/getBaseCoverageForBam.Rd index a8a3f66..612717f 100644 --- a/man/getBaseCoverageForBam.Rd +++ b/man/getBaseCoverageForBam.Rd @@ -10,8 +10,11 @@ getBaseCoverageForBam(bamfile, regions = NULL, maxDepth = 200L, nThreads = 3L) \item{bamfile}{A character scalar with the bam file name (and path).} \item{regions}{Character vector specifying the region(s) for which -to calculate coverage, in the form \code{"."}, \code{"chr"} or -\code{"chr:start-end"}. If \code{NULL}, the whole genome +to calculate coverage. Each region uses the grammar supported by the +\code{htslib} function \code{sam_parse_region}, for example +\code{"."} (all contigs), \code{"chr"} (whole contig), +\code{"chr:START"}, \code{"chr:-END"} or +\code{"chr:START-END"}. If \code{NULL}, the whole genome (\code{"."}) is used by default.} \item{maxDepth}{An integer scalar defining the maximal depth to consider.} diff --git a/src/utils.cpp b/src/utils.cpp index d31f7c5..d1b6260 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1139,13 +1139,21 @@ static void region_cov(samFile *fp, hts_itr_t *it, bam1_t *b, // Resolve a vector of region strings into (tid, beg, end) spans, where // `beg` is 0-based and `end` is exclusive. "." or "*" expands to all -// non-empty contigs. Returns 0 on success, -1 on failure (message in buffer). +// non-empty contigs. Individual region strings are parsed with htslib's +// sam_parse_region(), which accepts the forms: +// REF -> entire contig +// REF:START -> REF:START to end of contig +// REF:-END -> begin of contig to END +// REF:START-END +// Returns 0 on success, -1 on failure (message in buffer). static int resolve_regions_to_spans(const std::vector ®ionsvect, sam_hdr_t *hdr, std::vector> &spans, bool &had_error, int buffer_len, char *buffer) { int n = sam_hdr_nref(hdr); for (const std::string &rs : regionsvect) { + // htslib treats "." and "*" as "all contigs" via sam_parse_region, + // but we only want non-empty contigs, so handle them explicitly if (rs == "." || rs == "*") { spans.clear(); for (int t = 0; t < n; t++) { @@ -1154,46 +1162,18 @@ static int resolve_regions_to_spans(const std::vector ®ionsvect, } break; } - size_t colon = rs.find(':'); - if (colon == std::string::npos) { - int tid = sam_hdr_name2tid(hdr, rs.c_str()); - if (tid < 0) { - had_error = true; - snprintf(buffer, buffer_len, "Unknown contig: %s\n", rs.c_str()); - return -1; - } - hts_pos_t L = sam_hdr_tid2len(hdr, tid); - if (L > 0) spans.push_back(std::make_tuple(tid, (hts_pos_t)0, L)); - continue; - } - std::string name = rs.substr(0, colon); - std::string posstr = rs.substr(colon + 1); - size_t dash = posstr.find('-'); - if (dash == std::string::npos) { - had_error = true; - snprintf(buffer, buffer_len, - "Malformed region (expected 'name:beg-end'): %s\n", rs.c_str()); - return -1; - } - errno = 0; - long long beg_in = strtoll(posstr.substr(0, dash).c_str(), NULL, 10); - errno = 0; - long long end_in = strtoll(posstr.substr(dash + 1).c_str(), NULL, 10); - if (errno != 0 || beg_in < 1 || end_in < beg_in) { - had_error = true; - snprintf(buffer, buffer_len, - "Invalid region position in: %s\n", rs.c_str()); - return -1; - } - int tid = sam_hdr_name2tid(hdr, name.c_str()); - if (tid < 0) { + int tid = -1; + hts_pos_t beg = 0, end = 0; + const char *ok = sam_parse_region(hdr, rs.c_str(), &tid, &beg, &end, 0); + if (ok == NULL || tid < 0) { had_error = true; - snprintf(buffer, buffer_len, "Unknown contig: %s\n", name.c_str()); + snprintf(buffer, buffer_len, "Invalid region: %s\n", rs.c_str()); return -1; } + // Clamp end to contig length (htslib may return INT64_MAX for + // the "open end" forms REF and REF:START) hts_pos_t hdr_len = sam_hdr_tid2len(hdr, tid); - hts_pos_t beg = (hts_pos_t)beg_in - 1; // 0-based - hts_pos_t end = (hts_pos_t)end_in > hdr_len ? hdr_len : (hts_pos_t)end_in; // exclusive, clamped + if (end > hdr_len) end = hdr_len; if (end > beg) { spans.push_back(std::make_tuple(tid, beg, end)); } @@ -1206,8 +1186,11 @@ static int resolve_regions_to_spans(const std::vector ®ionsvect, //' //' @param bamfile A character scalar with the bam file name (and path). //' @param regions Character vector specifying the region(s) for which -//' to calculate coverage, in the form \code{"."}, \code{"chr"} or -//' \code{"chr:start-end"}. If \code{NULL}, the whole genome +//' to calculate coverage. Each region uses the grammar supported by the +//' \code{htslib} function \code{sam_parse_region}, for example +//' \code{"."} (all contigs), \code{"chr"} (whole contig), +//' \code{"chr:START"}, \code{"chr:-END"} or +//' \code{"chr:START-END"}. If \code{NULL}, the whole genome //' (\code{"."}) is used by default. //' @param maxDepth An integer scalar defining the maximal depth to consider. //' @param nThreads A numeric scalar with the number of threads used for diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R index 2b086cf..1f6b6f2 100644 --- a/tests/testthat/test-utils.R +++ b/tests/testthat/test-utils.R @@ -387,13 +387,9 @@ test_that("getBaseCoverageForBam works", { expect_error(getBaseCoverageForBam(bamfile = "NO_SUCH_FILE.bam"), "Could not open input file") expect_error(getBaseCoverageForBam(bamfile, regions = "notarealcontig"), - "Unknown contig") - expect_error(getBaseCoverageForBam(bamfile, regions = "chr1:6940000"), - "Malformed region") + "Invalid region") expect_error(getBaseCoverageForBam(bamfile, regions = "chr1:5-3"), - "Invalid region position") - expect_error(getBaseCoverageForBam(bamfile, regions = "chr1:-10"), - "Invalid region position") + "Invalid region") ## -- whole genome (NULL == ".") ------------------------------------------ ## vNULL <- getBaseCoverageForBam(bamfile, regions = NULL) @@ -455,4 +451,14 @@ test_that("getBaseCoverageForBam works", { v2t <- getBaseCoverageForBam( bamfile, regions = "chr1:6940000-6960000", nThreads = 2L) expect_identical(v2t, vRegion) + + ## -- htslib grammar: REF:START (open end) & REF:-END (open start) -------- ## + # "chr1:6940000" means chr1:6940000 to end of contig (clamped to length) + vOpenEnd <- getBaseCoverageForBam(bamfile, regions = "chr1:6940000") + # identical to the (6940000, contig-length) span + expect_identical(vOpenEnd, vClamp) + # "chr1:-10" means positions 1..10 of chr1 + vOpenStart <- getBaseCoverageForBam(bamfile, regions = "chr1:-10") + expect_equal(sum(vOpenStart), 10) + }) From 1111c836abad6b16468a6984ebbc1016609395ef Mon Sep 17 00:00:00 2001 From: mbstadler Date: Fri, 25 Sep 2026 08:05:28 +0200 Subject: [PATCH 07/13] add region_cov_simple method to getBaseCoverageForBam New method 'simple' ignores CIGAR strings and covers the whole reference span between the first and last aligned position of each alignment (via htslib's bam_endpos), so soft-clip and inserted bases also count as covered. 'full' remains the default and is unchanged. Benchmark on a 102 GB Dorado 6mA BAM (GRCm39, ~100x): simple is ~1.2x faster on a 1 Mb region and roughly tied at 10-100 Mb, so the win is modest for long ONT reads. Co-authored-by: Charlotte Soneson Co-authored-by: Michael Stadler --- R/RcppExports.R | 32 +++++++++++++--- man/getBaseCoverageForBam.Rd | 29 ++++++++++++-- src/RcppExports.cpp | 9 +++-- src/utils.cpp | 64 +++++++++++++++++++++++++++++-- src/utils.h | 9 +++++ tests/testthat/test-utils.R | 74 ++++++++++++++++++++++++++++++++++++ 6 files changed, 200 insertions(+), 17 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 556baeb..b17a595 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -716,6 +716,15 @@ NULL #' @keywords internal NULL +#' Simple coverage: a read covers the whole reference span between its first +#' and last aligned position, regardless of CIGAR operations (soft-clip bases +#' and read-inserted bases are also counted as covered). htslib's +#' bam_endpos() gives exactly this reference span (pos-1 .. endpos-2, 0-based +#' inclusive; it equals pos-1 if the CIGAR ends in a soft-clip). +#' @noRd +#' @keywords internal +NULL + #' Concatenate files #' #' @param input_files Character vector with input file names to concatenate. @@ -797,25 +806,38 @@ complement <- function(n) { #' \code{"chr:START-END"}. If \code{NULL}, the whole genome #' (\code{"."}) is used by default. #' @param maxDepth An integer scalar defining the maximal depth to consider. +#' @param method Character scalar with the method used for coverage +#' calculation. \code{"full"} (the default) considers CIGAR +#' operations, thus not counting soft-clip bases and read-inserted +#' positions as covered. \code{"simple"} ignores the CIGAR strings and +#' covers the whole reference span between the first and the last +#' aligned position of each alignment (using \code{htslib}'s +#' \code{bam_endpos}); this is slightly faster but overestimates +#' coverage in the soft-clipped ends and around indels. #' @param nThreads A numeric scalar with the number of threads used for #' decompressing BAM records. #' -#' @details CIGAR operations are considered, thus not counting the genomic -#' bases in a read-insertion as covered. Secondary and supplementary -#' alignments are not included. +#' @details Secondary and supplementary alignments and unmapped reads are +#' not included. #' #' @references The algorithm was described in Pedersen BS and Quinlan AR. #' "Mosdepth: quick coverage calculation for genomes and exomes". #' Bioinformatics. 2018; 34(5):867-868. #' \url{https://doi.org/10.1093/bioinformatics/btx699} #' +#' @examples +#' modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") +#' getBaseCoverageForBam(modbamfile, "chr1", 12L, "full") +#' getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "full") +#' getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "simple") +#' #' @return A numeric vector of length \code{maxDepth + 1}, with values at #' index \code{i} giving the number of positions that were overlapped by #' exactly \code{i} alignments. Positions overlapped by more than #' \code{maxDepth} alignments are also added to the value for \code{maxDepth} #' at index \code{maxDepth + 1}. #' -getBaseCoverageForBam <- function(bamfile, regions = NULL, maxDepth = 200L, nThreads = 3L) { - .Call(`_SingleMoleculeGenomicsIO_getBaseCoverageForBam`, bamfile, regions, maxDepth, nThreads) +getBaseCoverageForBam <- function(bamfile, regions = NULL, maxDepth = 200L, method = "full", nThreads = 3L) { + .Call(`_SingleMoleculeGenomicsIO_getBaseCoverageForBam`, bamfile, regions, maxDepth, method, nThreads) } diff --git a/man/getBaseCoverageForBam.Rd b/man/getBaseCoverageForBam.Rd index 612717f..b8c6dca 100644 --- a/man/getBaseCoverageForBam.Rd +++ b/man/getBaseCoverageForBam.Rd @@ -4,7 +4,13 @@ \alias{getBaseCoverageForBam} \title{Get base coverage histogram for a BAM file} \usage{ -getBaseCoverageForBam(bamfile, regions = NULL, maxDepth = 200L, nThreads = 3L) +getBaseCoverageForBam( + bamfile, + regions = NULL, + maxDepth = 200L, + method = "full", + nThreads = 3L +) } \arguments{ \item{bamfile}{A character scalar with the bam file name (and path).} @@ -19,6 +25,15 @@ to calculate coverage. Each region uses the grammar supported by the \item{maxDepth}{An integer scalar defining the maximal depth to consider.} +\item{method}{Character scalar with the method used for coverage +calculation. \code{"full"} (the default) considers CIGAR +operations, thus not counting soft-clip bases and read-inserted +positions as covered. \code{"simple"} ignores the CIGAR strings and +covers the whole reference span between the first and the last +aligned position of each alignment (using \code{htslib}'s +\code{bam_endpos}); this is slightly faster but overestimates +coverage in the soft-clipped ends and around indels.} + \item{nThreads}{A numeric scalar with the number of threads used for decompressing BAM records.} } @@ -33,9 +48,15 @@ at index \code{maxDepth + 1}. Get base coverage histogram for a BAM file } \details{ -CIGAR operations are considered, thus not counting the genomic -bases in a read-insertion as covered. Secondary and supplementary -alignments are not included. +Secondary and supplementary alignments and unmapped reads are +not included. +} +\examples{ +modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") +getBaseCoverageForBam(modbamfile, "chr1", 12L, "full") +getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "full") +getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "simple") + } \references{ The algorithm was described in Pedersen BS and Quinlan AR. diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index b904184..e00a3aa 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -248,16 +248,17 @@ BEGIN_RCPP END_RCPP } // getBaseCoverageForBam -Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions, const uint maxDepth, int nThreads); -RcppExport SEXP _SingleMoleculeGenomicsIO_getBaseCoverageForBam(SEXP bamfileSEXP, SEXP regionsSEXP, SEXP maxDepthSEXP, SEXP nThreadsSEXP) { +Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions, const uint maxDepth, const std::string method, int nThreads); +RcppExport SEXP _SingleMoleculeGenomicsIO_getBaseCoverageForBam(SEXP bamfileSEXP, SEXP regionsSEXP, SEXP maxDepthSEXP, SEXP methodSEXP, SEXP nThreadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const std::string >::type bamfile(bamfileSEXP); Rcpp::traits::input_parameter< Rcpp::Nullable> >::type regions(regionsSEXP); Rcpp::traits::input_parameter< const uint >::type maxDepth(maxDepthSEXP); + Rcpp::traits::input_parameter< const std::string >::type method(methodSEXP); Rcpp::traits::input_parameter< int >::type nThreads(nThreadsSEXP); - rcpp_result_gen = Rcpp::wrap(getBaseCoverageForBam(bamfile, regions, maxDepth, nThreads)); + rcpp_result_gen = Rcpp::wrap(getBaseCoverageForBam(bamfile, regions, maxDepth, method, nThreads)); return rcpp_result_gen; END_RCPP } @@ -277,7 +278,7 @@ static const R_CallMethodDef CallEntries[] = { {"_SingleMoleculeGenomicsIO_getChromosomeNamesFromBam", (DL_FUNC) &_SingleMoleculeGenomicsIO_getChromosomeNamesFromBam, 1}, {"_SingleMoleculeGenomicsIO_get_unmodified_base", (DL_FUNC) &_SingleMoleculeGenomicsIO_get_unmodified_base, 1}, {"_SingleMoleculeGenomicsIO_complement", (DL_FUNC) &_SingleMoleculeGenomicsIO_complement, 1}, - {"_SingleMoleculeGenomicsIO_getBaseCoverageForBam", (DL_FUNC) &_SingleMoleculeGenomicsIO_getBaseCoverageForBam, 4}, + {"_SingleMoleculeGenomicsIO_getBaseCoverageForBam", (DL_FUNC) &_SingleMoleculeGenomicsIO_getBaseCoverageForBam, 5}, {NULL, NULL, 0} }; diff --git a/src/utils.cpp b/src/utils.cpp index d1b6260..3c28aa3 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include "utils.h" @@ -1137,6 +1138,37 @@ static void region_cov(samFile *fp, hts_itr_t *it, bam1_t *b, } } +//' Simple coverage: a read covers the whole reference span between its first +//' and last aligned position, regardless of CIGAR operations (soft-clip bases +//' and read-inserted bases are also counted as covered). htslib's +//' bam_endpos() gives exactly this reference span (pos-1 .. endpos-2, 0-based +//' inclusive; it equals pos-1 if the CIGAR ends in a soft-clip). +//' @noRd +//' @keywords internal +static void region_cov_simple(samFile *fp, hts_itr_t *it, bam1_t *b, + hts_pos_t beg, hts_pos_t end, int32_t *diff, + uint64_t *hist, uint maxd) { + hts_pos_t len = end - beg; + memset(diff, 0, (len + 1) * sizeof(int32_t)); + + while (sam_itr_next(fp, it, b) >= 0) { + if (b->core.flag & (BAM_FUNMAP | BAM_FSECONDARY | BAM_FSUPPLEMENTARY)) { + continue; + } + hts_pos_t s_rel = b->core.pos - beg; + hts_pos_t e_rel = bam_endpos(b) - beg; + if (e_rel > 0 && s_rel < len) { + diff[s_rel < 0 ? 0 : s_rel]++; + diff[e_rel > len ? len : e_rel]--; + } + } + int32_t d = 0; + for (hts_pos_t i = 0; i < len; i++) { + d += diff[i]; + hist[d < maxd ? d : maxd]++; + } +} + // Resolve a vector of region strings into (tid, beg, end) spans, where // `beg` is 0-based and `end` is exclusive. "." or "*" expands to all // non-empty contigs. Individual region strings are parsed with htslib's @@ -1193,18 +1225,31 @@ static int resolve_regions_to_spans(const std::vector ®ionsvect, //' \code{"chr:START-END"}. If \code{NULL}, the whole genome //' (\code{"."}) is used by default. //' @param maxDepth An integer scalar defining the maximal depth to consider. +//' @param method Character scalar with the method used for coverage +//' calculation. \code{"full"} (the default) considers CIGAR +//' operations, thus not counting soft-clip bases and read-inserted +//' positions as covered. \code{"simple"} ignores the CIGAR strings and +//' covers the whole reference span between the first and the last +//' aligned position of each alignment (using \code{htslib}'s +//' \code{bam_endpos}); this is slightly faster but overestimates +//' coverage in the soft-clipped ends and around indels. //' @param nThreads A numeric scalar with the number of threads used for //' decompressing BAM records. //' -//' @details CIGAR operations are considered, thus not counting the genomic -//' bases in a read-insertion as covered. Secondary and supplementary -//' alignments are not included. +//' @details Secondary and supplementary alignments and unmapped reads are +//' not included. //' //' @references The algorithm was described in Pedersen BS and Quinlan AR. //' "Mosdepth: quick coverage calculation for genomes and exomes". //' Bioinformatics. 2018; 34(5):867-868. //' \url{https://doi.org/10.1093/bioinformatics/btx699} //' +//' @examples +//' modbamfile <- system.file("extdata", "6mA_1_10reads.bam", package = "SingleMoleculeGenomicsIO") +//' getBaseCoverageForBam(modbamfile, "chr1", 12L, "full") +//' getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "full") +//' getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "simple") +//' //' @return A numeric vector of length \code{maxDepth + 1}, with values at //' index \code{i} giving the number of positions that were overlapped by //' exactly \code{i} alignments. Positions overlapped by more than @@ -1215,8 +1260,15 @@ static int resolve_regions_to_spans(const std::vector ®ionsvect, Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions = R_NilValue, const uint maxDepth = 200, + const std::string method = "full", int nThreads = 3) { + std::vector allowed_methods = {"full", "simple"}; + if (std::find(allowed_methods.begin(), allowed_methods.end(), method) == allowed_methods.end()) { + Rcpp::stop(std::string("Invalid method: '") + method + + "' (only 'full' and 'simple' are supported)"); + } + int buffer_len = 2000, success = 0; char buffer[2000]; bool had_error = false; @@ -1284,7 +1336,11 @@ Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, goto end; } - region_cov(inbamfile, iter, bamdata, beg, end, diff, hist, maxDepth); + if (method == "full") { + region_cov(inbamfile, iter, bamdata, beg, end, diff, hist, maxDepth); + } else { // method == "simple" + region_cov_simple(inbamfile, iter, bamdata, beg, end, diff, hist, maxDepth); + } if (iter) { sam_itr_destroy(iter); diff --git a/src/utils.h b/src/utils.h index 6f0c532..de17f13 100644 --- a/src/utils.h +++ b/src/utils.h @@ -96,7 +96,16 @@ static void region_cov(samFile*, int32_t*, uint64_t*, uint); +static void region_cov_simple(samFile*, + hts_itr_t*, + bam1_t*, + hts_pos_t, + hts_pos_t, + int32_t*, + uint64_t*, + uint); Rcpp::NumericVector getBaseCoverageForBam(const std::string, Rcpp::Nullable>, const uint, + const std::string, int); diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R index 1f6b6f2..0fcb376 100644 --- a/tests/testthat/test-utils.R +++ b/tests/testthat/test-utils.R @@ -462,3 +462,77 @@ test_that("getBaseCoverageForBam works", { expect_equal(sum(vOpenStart), 10) }) + +## Checks, getBaseCoverageForBam method = "simple" +## -------------------------------------------------------------------------- ## +test_that("getBaseCoverageForBam method='simple' works", { + bamfile <- system.file("extdata", "6mA_1_10reads.bam", + package = "SingleMoleculeGenomicsIO") + expect_true(file.exists(bamfile)) + + ## -- default method is "full" -------------------------------------------- ## + vDefault <- getBaseCoverageForBam(bamfile, regions = "chr1:6940000-6960000") + vFull <- getBaseCoverageForBam( + bamfile, regions = "chr1:6940000-6960000", method = "full") + expect_identical(vDefault, vFull) + + ## -- invalid method ------------------------------------------------------- ## + expect_error(getBaseCoverageForBam(bamfile, method = "bogus"), + "Invalid method") + + ## -- "simple" on a region spanned by reads with soft-clips ---------------- ## + # reads aligned in this window have non-zero soft-clips, so the simple + # method (which covers the full first..last aligned reference span, + # including soft-clip bases and inserted bases) gives a slightly + # different histogram than the full method. + # Verified ground truth: + # full -> c(18370, 841, 651, 139, 0) + # simple -> c(18369, 829, 660, 143, 0) + vRegionFull <- getBaseCoverageForBam( + bamfile, regions = "chr1:6940000-6960000", method = "full") + vRegionSimple <- getBaseCoverageForBam( + bamfile, regions = "chr1:6940000-6960000", method = "simple") + expect_identical(vRegionFull, c(18370, 841, 651, 139, + rep(0, 197))) + expect_identical(vRegionSimple, c(18369, 829, 660, 143, + rep(0, 197))) + # both cover all positions in the region + expect_equal(sum(vRegionSimple), 20001) + # simple does not create depth > full in this case (only soft-clip + # positions change, and max depth in full is already 4 here) + expect_true(max(which(vRegionSimple > 0)) >= max(which(vRegionFull > 0))) + + ## -- "simple" on whole contig -------------------------------------------- ## + vChr1Full <- getBaseCoverageForBam( + bamfile, regions = "chr1", method = "full") + vChr1Simple <- getBaseCoverageForBam( + bamfile, regions = "chr1", method = "simple") + # both cover the whole contig + expect_equal(sum(vChr1Full), 195154279) + expect_equal(sum(vChr1Simple), 195154279) + # simple differs from full (reads have soft-clips) + expect_false(isTRUE(all.equal(vChr1Full, vChr1Simple))) + # specific first 5 bins + expect_identical(head(vChr1Simple, 5), + c(195138477, 3839, 823, 538, 991)) + + ## -- "simple" on maxDepth overflow bin ------------------------------------ ## + vSmall <- getBaseCoverageForBam( + bamfile, regions = "chr1:6940000-6960000", + method = "simple", maxDepth = 2L) + # simple's max depth is also 4; index 2 catches depth 3+4 bins + expect_identical(vSmall, c(18369, 829, 660 + 143)) + + ## -- "simple" on a region fully outside any read span -------------------- ## + # positions never covered by any read -> all positions go into the + # zero-depth bin (index 1), same as full + vNoneFull <- getBaseCoverageForBam( + bamfile, regions = "chr1:10000000-10001000") + vNoneSimple <- getBaseCoverageForBam( + bamfile, regions = "chr1:10000000-10001000", method = "simple") + expect_identical(vNoneFull, vNoneSimple) + expect_equal(vNoneSimple[1L], 1001) # all positions at zero depth + expect_true(all(vNoneSimple[-1L] == 0)) + expect_equal(sum(vNoneSimple), 1001) # region length +}) + From 9f73f3fae17428c922a5a4aeb9a782ac1c6442e1 Mon Sep 17 00:00:00 2001 From: mbstadler Date: Fri, 25 Sep 2026 08:37:13 +0200 Subject: [PATCH 08/13] add names to return value of getBaseCoverageForBam and update unit tests Co-authored-by: Charlotte Soneson Co-authored-by: Michael Stadler --- R/RcppExports.R | 4 ++-- man/getBaseCoverageForBam.Rd | 4 ++-- src/utils.cpp | 7 +++++-- src/utils.h | 16 ---------------- tests/testthat/test-utils.R | 24 ++++++++++++------------ 5 files changed, 21 insertions(+), 34 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index b17a595..c521642 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -831,9 +831,9 @@ complement <- function(n) { #' getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "full") #' getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "simple") #' -#' @return A numeric vector of length \code{maxDepth + 1}, with values at +#' @return A named numeric vector of length \code{maxDepth + 1}, with values at #' index \code{i} giving the number of positions that were overlapped by -#' exactly \code{i} alignments. Positions overlapped by more than +#' exactly \code{i-1} alignments. Positions overlapped by more than #' \code{maxDepth} alignments are also added to the value for \code{maxDepth} #' at index \code{maxDepth + 1}. #' diff --git a/man/getBaseCoverageForBam.Rd b/man/getBaseCoverageForBam.Rd index b8c6dca..7885f93 100644 --- a/man/getBaseCoverageForBam.Rd +++ b/man/getBaseCoverageForBam.Rd @@ -38,9 +38,9 @@ coverage in the soft-clipped ends and around indels.} decompressing BAM records.} } \value{ -A numeric vector of length \code{maxDepth + 1}, with values at +A named numeric vector of length \code{maxDepth + 1}, with values at index \code{i} giving the number of positions that were overlapped by -exactly \code{i} alignments. Positions overlapped by more than +exactly \code{i-1} alignments. Positions overlapped by more than \code{maxDepth} alignments are also added to the value for \code{maxDepth} at index \code{maxDepth + 1}. } diff --git a/src/utils.cpp b/src/utils.cpp index 3c28aa3..4c9e768 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1250,9 +1250,9 @@ static int resolve_regions_to_spans(const std::vector ®ionsvect, //' getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "full") //' getBaseCoverageForBam(modbamfile, "chr1:6000000-7000000", 12L, "simple") //' -//' @return A numeric vector of length \code{maxDepth + 1}, with values at +//' @return A named numeric vector of length \code{maxDepth + 1}, with values at //' index \code{i} giving the number of positions that were overlapped by -//' exactly \code{i} alignments. Positions overlapped by more than +//' exactly \code{i-1} alignments. Positions overlapped by more than //' \code{maxDepth} alignments are also added to the value for \code{maxDepth} //' at index \code{maxDepth + 1}. //' @@ -1282,6 +1282,9 @@ Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, sam_hdr_t *inbamhdr = NULL; uint64_t *hist = (uint64_t*)calloc(maxDepth + 1, sizeof(uint64_t)); // hist[maxDepth] = overflow bin Rcpp::NumericVector histvect(maxDepth + 1); + Rcpp::IntegerVector histvectNamesInt = Rcpp::seq(0, maxDepth); + Rcpp::CharacterVector histvectNames = Rcpp::as(histvectNamesInt); + histvect.attr("names") = histvectNames; // set default regions if NULL if (regions.isNotNull()) { diff --git a/src/utils.h b/src/utils.h index de17f13..3b5424d 100644 --- a/src/utils.h +++ b/src/utils.h @@ -88,22 +88,6 @@ int check_bam_format(samFile *infile, char *buffer, int &buffer_len); -static void region_cov(samFile*, - hts_itr_t*, - bam1_t*, - hts_pos_t, - hts_pos_t, - int32_t*, - uint64_t*, - uint); -static void region_cov_simple(samFile*, - hts_itr_t*, - bam1_t*, - hts_pos_t, - hts_pos_t, - int32_t*, - uint64_t*, - uint); Rcpp::NumericVector getBaseCoverageForBam(const std::string, Rcpp::Nullable>, const uint, diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R index 0fcb376..d0aefec 100644 --- a/tests/testthat/test-utils.R +++ b/tests/testthat/test-utils.R @@ -403,7 +403,7 @@ test_that("getBaseCoverageForBam works", { # ground truth (verified against samtools depth -a) expect_identical(head(vNULL, 11), - c(2728206677, 3841, 811, 563, 973, 596, 675, 463, 420, 1066, 6366)) + setNames(c(2728206677, 3841, 811, 563, 973, 596, 675, 463, 420, 1066, 6366), 0:10)) expect_true(all(vNULL[12:201] == 0)) # max depth is 10 for 10 reads ## -- single contig (no coordinates) -------------------------------------- ## @@ -411,14 +411,14 @@ test_that("getBaseCoverageForBam works", { # the whole genome, but with zero-count restricted to chr1 length. vChr1 <- getBaseCoverageForBam(bamfile, regions = "chr1") expect_identical(head(vChr1, 11), - c(195138505, 3841, 811, 563, 973, 596, 675, 463, 420, 1066, 6366)) + setNames(c(195138505, 3841, 811, 563, 973, 596, 675, 463, 420, 1066, 6366), 0:10)) expect_true(all(vChr1[12:201] == 0)) expect_equal(sum(vChr1), 195154279) # contig length ## -- bounded region ------------------------------------------------------ ## vRegion <- getBaseCoverageForBam(bamfile, regions = "chr1:6940000-6960000") expect_identical(head(vRegion, 5), - c(18370, 841, 651, 139, 0)) + setNames(c(18370, 841, 651, 139, 0), 0:4)) expect_true(all(vRegion[6:201] == 0)) expect_equal(sum(vRegion), 20001) # total positions in region @@ -426,7 +426,7 @@ test_that("getBaseCoverageForBam works", { vClamp <- getBaseCoverageForBam( bamfile, regions = "chr1:6940000-999999999") expect_identical(head(vClamp, 4), - c(188212649, 841, 651, 139)) + setNames(c(188212649, 841, 651, 139), 0:3)) # clamped to chr1 length, so sum = 195154279 - 6939999 expect_equal(sum(vClamp), 195154279 - 6939999) @@ -445,7 +445,7 @@ test_that("getBaseCoverageForBam works", { bamfile, regions = "chr1:6940000-6960000", maxDepth = 2L) expect_length(vSmall, 3L) # index 2 = overflow: positions covered by depth 3 AND 4 - expect_identical(vSmall, c(18370, 841, 651 + 139)) + expect_identical(vSmall, setNames(c(18370, 841, 651 + 139), 0:2)) ## -- nThreads parameter does not affect result ---------------------------- ## v2t <- getBaseCoverageForBam( @@ -492,10 +492,10 @@ test_that("getBaseCoverageForBam method='simple' works", { bamfile, regions = "chr1:6940000-6960000", method = "full") vRegionSimple <- getBaseCoverageForBam( bamfile, regions = "chr1:6940000-6960000", method = "simple") - expect_identical(vRegionFull, c(18370, 841, 651, 139, - rep(0, 197))) - expect_identical(vRegionSimple, c(18369, 829, 660, 143, - rep(0, 197))) + expect_identical(vRegionFull, + setNames(c(18370, 841, 651, 139, rep(0, 197)), 0:200)) + expect_identical(vRegionSimple, + setNames(c(18369, 829, 660, 143, rep(0, 197)), 0:200)) # both cover all positions in the region expect_equal(sum(vRegionSimple), 20001) # simple does not create depth > full in this case (only soft-clip @@ -514,14 +514,14 @@ test_that("getBaseCoverageForBam method='simple' works", { expect_false(isTRUE(all.equal(vChr1Full, vChr1Simple))) # specific first 5 bins expect_identical(head(vChr1Simple, 5), - c(195138477, 3839, 823, 538, 991)) + setNames(c(195138477, 3839, 823, 538, 991), 0:4)) ## -- "simple" on maxDepth overflow bin ------------------------------------ ## vSmall <- getBaseCoverageForBam( bamfile, regions = "chr1:6940000-6960000", method = "simple", maxDepth = 2L) # simple's max depth is also 4; index 2 catches depth 3+4 bins - expect_identical(vSmall, c(18369, 829, 660 + 143)) + expect_identical(vSmall, setNames(c(18369, 829, 660 + 143), 0:2)) ## -- "simple" on a region fully outside any read span -------------------- ## # positions never covered by any read -> all positions go into the @@ -531,7 +531,7 @@ test_that("getBaseCoverageForBam method='simple' works", { vNoneSimple <- getBaseCoverageForBam( bamfile, regions = "chr1:10000000-10001000", method = "simple") expect_identical(vNoneFull, vNoneSimple) - expect_equal(vNoneSimple[1L], 1001) # all positions at zero depth + expect_equal(vNoneSimple[1L], c("0" = 1001)) # all positions at zero depth expect_true(all(vNoneSimple[-1L] == 0)) expect_equal(sum(vNoneSimple), 1001) # region length }) From 865124702e503c56872ef9e33fba6e2021c85a18 Mon Sep 17 00:00:00 2001 From: mbstadler Date: Fri, 25 Sep 2026 08:42:55 +0200 Subject: [PATCH 09/13] add nocov to error cases not supported by current test data Co-authored-by: Charlotte Soneson Co-authored-by: Michael Stadler --- src/utils.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/utils.cpp b/src/utils.cpp index 4c9e768..72592a7 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1107,7 +1107,7 @@ static void region_cov(samFile *fp, hts_itr_t *it, bam1_t *b, while (sam_itr_next(fp, it, b) >= 0) { if (b->core.flag & (BAM_FUNMAP | BAM_FSECONDARY | BAM_FSUPPLEMENTARY)) { - continue; + continue; // # nocov } hts_pos_t pos = b->core.pos; const uint32_t *cig = bam_get_cigar(b); @@ -1153,7 +1153,7 @@ static void region_cov_simple(samFile *fp, hts_itr_t *it, bam1_t *b, while (sam_itr_next(fp, it, b) >= 0) { if (b->core.flag & (BAM_FUNMAP | BAM_FSECONDARY | BAM_FSUPPLEMENTARY)) { - continue; + continue; // # nocov } hts_pos_t s_rel = b->core.pos - beg; hts_pos_t e_rel = bam_endpos(b) - beg; @@ -1319,24 +1319,24 @@ Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, hts_pos_t end = std::get<2>(sp); hts_pos_t len = end - beg; if (len == 0) { - continue; + continue; // # nocov } int32_t *diff = (int32_t*)malloc((len + 1) * sizeof(int32_t)); if (diff == NULL) { - had_error = true; + had_error = true; // # nocov start snprintf(buffer, buffer_len, "Failed to allocate memory for coverage vector\n"); - goto end; + goto end; // # nocov end } iter = sam_itr_queryi(idx, tid, beg, end); if (iter == NULL) { - free(diff); + free(diff); // # nocov start diff = NULL; had_error = true; snprintf(buffer, buffer_len, "Failed to create region iterator\n"); - goto end; + goto end; // # nocov end } if (method == "full") { From 3e87b4de303876461124016db955e46e2122ba2b Mon Sep 17 00:00:00 2001 From: Charlotte Soneson Date: Fri, 25 Sep 2026 09:00:50 +0200 Subject: [PATCH 10/13] Use uint32_t instead of uint --- src/RcppExports.cpp | 4 ++-- src/utils.cpp | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 88d376f..f30f1d4 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -259,14 +259,14 @@ BEGIN_RCPP END_RCPP } // getBaseCoverageForBam -Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions, const uint maxDepth, const std::string method, int nThreads); +Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions, const uint32_t maxDepth, const std::string method, int nThreads); RcppExport SEXP _SingleMoleculeGenomicsIO_getBaseCoverageForBam(SEXP bamfileSEXP, SEXP regionsSEXP, SEXP maxDepthSEXP, SEXP methodSEXP, SEXP nThreadsSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const std::string >::type bamfile(bamfileSEXP); Rcpp::traits::input_parameter< Rcpp::Nullable> >::type regions(regionsSEXP); - Rcpp::traits::input_parameter< const uint >::type maxDepth(maxDepthSEXP); + Rcpp::traits::input_parameter< const uint32_t >::type maxDepth(maxDepthSEXP); Rcpp::traits::input_parameter< const std::string >::type method(methodSEXP); Rcpp::traits::input_parameter< int >::type nThreads(nThreadsSEXP); rcpp_result_gen = Rcpp::wrap(getBaseCoverageForBam(bamfile, regions, maxDepth, method, nThreads)); diff --git a/src/utils.cpp b/src/utils.cpp index 99c4b0b..baea80e 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1202,7 +1202,7 @@ int check_bam_format(samFile *infile, //' @keywords internal static void region_cov(samFile *fp, hts_itr_t *it, bam1_t *b, hts_pos_t beg, hts_pos_t end, int32_t *diff, - uint64_t *hist, uint maxd) { + uint64_t *hist, uint32_t maxd) { hts_pos_t len = end - beg; memset(diff, 0, (len + 1) * sizeof(int32_t)); @@ -1248,7 +1248,7 @@ static void region_cov(samFile *fp, hts_itr_t *it, bam1_t *b, //' @keywords internal static void region_cov_simple(samFile *fp, hts_itr_t *it, bam1_t *b, hts_pos_t beg, hts_pos_t end, int32_t *diff, - uint64_t *hist, uint maxd) { + uint64_t *hist, uint32_t maxd) { hts_pos_t len = end - beg; memset(diff, 0, (len + 1) * sizeof(int32_t)); @@ -1360,7 +1360,7 @@ static int resolve_regions_to_spans(const std::vector ®ionsvect, // [[Rcpp::export]] Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions = R_NilValue, - const uint maxDepth = 200, + const uint32_t maxDepth = 200, const std::string method = "full", int nThreads = 3) { From e3464cf03cada395cc5f73c69693a1225281d8eb Mon Sep 17 00:00:00 2001 From: Charlotte Soneson Date: Fri, 25 Sep 2026 09:07:36 +0200 Subject: [PATCH 11/13] Update also utils.h --- src/utils.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.h b/src/utils.h index 3b5424d..1208e93 100644 --- a/src/utils.h +++ b/src/utils.h @@ -90,6 +90,6 @@ int check_bam_format(samFile *infile, Rcpp::NumericVector getBaseCoverageForBam(const std::string, Rcpp::Nullable>, - const uint, + const uint32_t, const std::string, int); From c5538da715d649ab4069ce2ba1f2778923d66e4b Mon Sep 17 00:00:00 2001 From: Charlotte Soneson Date: Fri, 25 Sep 2026 09:19:25 +0200 Subject: [PATCH 12/13] Export getBaseCoverageForBam --- NAMESPACE | 1 + R/RcppExports.R | 1 + src/utils.cpp | 1 + 3 files changed, 3 insertions(+) diff --git a/NAMESPACE b/NAMESPACE index dcc7557..68bfb1b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -18,6 +18,7 @@ export(filterPositions) export(filterReads) export(filterReadsModBam) export(flattenReadLevelAssay) +export(getBaseCoverageForBam) export(getReadLevelAssayNames) export(getReadLevelColDataNames) export(getReadNamesBySample) diff --git a/R/RcppExports.R b/R/RcppExports.R index 937c6b5..540b2d2 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -858,6 +858,7 @@ complement <- function(n) { #' \code{maxDepth} alignments are also added to the value for \code{maxDepth} #' at index \code{maxDepth + 1}. #' +#' @export getBaseCoverageForBam <- function(bamfile, regions = NULL, maxDepth = 200L, method = "full", nThreads = 3L) { .Call(`_SingleMoleculeGenomicsIO_getBaseCoverageForBam`, bamfile, regions, maxDepth, method, nThreads) } diff --git a/src/utils.cpp b/src/utils.cpp index baea80e..cd496d8 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -1357,6 +1357,7 @@ static int resolve_regions_to_spans(const std::vector ®ionsvect, //' \code{maxDepth} alignments are also added to the value for \code{maxDepth} //' at index \code{maxDepth + 1}. //' +//' @export // [[Rcpp::export]] Rcpp::NumericVector getBaseCoverageForBam(const std::string bamfile, Rcpp::Nullable> regions = R_NilValue, From 5d8bee93efe5e455f8162b00d1eb19014fd3de18 Mon Sep 17 00:00:00 2001 From: Charlotte Soneson Date: Fri, 25 Sep 2026 09:39:21 +0200 Subject: [PATCH 13/13] Bump version, update NEWS --- DESCRIPTION | 2 +- NEWS.md | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index f8fd9f8..d1bb933 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: SingleMoleculeGenomicsIO Title: Tools for reading single molecule genomics data -Version: 0.1.7 +Version: 0.1.8 Authors@R: c( person("Michael", "Stadler", , "michael.stadler@fmi.ch", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-2269-4934")), diff --git a/NEWS.md b/NEWS.md index 337ca79..1b3aaae 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# SingleMoleculeGenomicsIO 0.1.8 + +* Add `getBaseCoverageForBam` and `countAlignmentsInBam` functions + # SingleMoleculeGenomicsIO 0.1.7 * Default to sample reads from all chromosomes in the bam files if no explicit sequence name(s) are given