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/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/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 diff --git a/R/RcppExports.R b/R/RcppExports.R index 0a25135..540b2d2 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -711,6 +711,20 @@ NULL #' @keywords internal NULL +#' Helper function for getBaseCoverageForBam +#' @noRd +#' @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. @@ -802,3 +816,50 @@ 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. 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 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 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 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-1} alignments. Positions overlapped by more than +#' \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/man/getBaseCoverageForBam.Rd b/man/getBaseCoverageForBam.Rd new file mode 100644 index 0000000..7885f93 --- /dev/null +++ b/man/getBaseCoverageForBam.Rd @@ -0,0 +1,66 @@ +% 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, + method = "full", + 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. 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.} + +\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.} +} +\value{ +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-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{ +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. +"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/RcppExports.cpp b/src/RcppExports.cpp index 623df6f..f30f1d4 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -258,6 +258,21 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// getBaseCoverageForBam +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 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)); + return rcpp_result_gen; +END_RCPP +} static const R_CallMethodDef CallEntries[] = { {"_SingleMoleculeGenomicsIO_estimateNoise", (DL_FUNC) &_SingleMoleculeGenomicsIO_estimateNoise, 4}, @@ -275,6 +290,7 @@ static const R_CallMethodDef CallEntries[] = { {"_SingleMoleculeGenomicsIO_countAlignmentsInBam", (DL_FUNC) &_SingleMoleculeGenomicsIO_countAlignmentsInBam, 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, 5}, {NULL, NULL, 0} }; diff --git a/src/utils.cpp b/src/utils.cpp index f4669ac..cd496d8 100644 --- a/src/utils.cpp +++ b/src/utils.cpp @@ -3,6 +3,8 @@ #include #include #include +#include +#include #include #include "utils.h" @@ -264,6 +266,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*. @@ -1192,3 +1195,283 @@ 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, + uint64_t *hist, uint32_t 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; // # nocov + } + 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]); + if (op == BAM_CMATCH || op == BAM_CEQUAL || op == BAM_CDIFF) { + // 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 += opl; + } else if (op == BAM_CDEL || op == BAM_CREF_SKIP) { + pos += opl; + } + // S / I / P do not advance the reference position + } + } + int32_t d = 0; + for (hts_pos_t i = 0; i < len; i++) { + d += diff[i]; + hist[d < maxd ? d : maxd]++; + } +} + +//' 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, uint32_t 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; // # nocov + } + 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 +// 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++) { + 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; + } + 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, "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); + if (end > hdr_len) end = hdr_len; + 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. 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 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 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 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-1} alignments. Positions overlapped by more than +//' \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, + const uint32_t 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; + 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; + 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()) { + 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; + } + + // 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; + } + + // 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; // # nocov + } + + int32_t *diff = (int32_t*)malloc((len + 1) * sizeof(int32_t)); + if (diff == NULL) { + had_error = true; // # nocov start + snprintf(buffer, buffer_len, + "Failed to allocate memory for coverage vector\n"); + goto end; // # nocov end + } + + iter = sam_itr_queryi(idx, tid, beg, end); + if (iter == NULL) { + free(diff); // # nocov start + diff = NULL; + had_error = true; + snprintf(buffer, buffer_len, "Failed to create region iterator\n"); + goto end; // # nocov end + } + + 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); + } + iter = NULL; + free(diff); + diff = NULL; + } + + // copy result into return vector before releasing `hist` + for (int i = 0; i <= maxDepth; i++) { + histvect[i] = (double) hist[i]; + } + + end: + // cleanup + if (hist) free(hist); + if (bamdata) bam_destroy1(bamdata); + if (idx) hts_idx_destroy(idx); + if (inbamhdr) sam_hdr_destroy(inbamhdr); + if (inbamfile) sam_close(inbamfile); + + 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..1208e93 100644 --- a/src/utils.h +++ b/src/utils.h @@ -87,3 +87,9 @@ int check_bam_format(samFile *infile, bool &had_error, char *buffer, int &buffer_len); + +Rcpp::NumericVector getBaseCoverageForBam(const std::string, + Rcpp::Nullable>, + const uint32_t, + const std::string, + int); diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R index b815246..41b4e21 100644 --- a/tests/testthat/test-utils.R +++ b/tests/testthat/test-utils.R @@ -420,3 +420,165 @@ 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"), + "Invalid region") + expect_error(getBaseCoverageForBam(bamfile, regions = "chr1:5-3"), + "Invalid region") + + ## -- 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), + 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) -------------------------------------- ## + # 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), + 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), + 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 + + ## -- region end beyond contig end is clamped ------------------------------ ## + vClamp <- getBaseCoverageForBam( + bamfile, regions = "chr1:6940000-999999999") + expect_identical(head(vClamp, 4), + setNames(c(188212649, 841, 651, 139), 0:3)) + # 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, setNames(c(18370, 841, 651 + 139), 0:2)) + + ## -- nThreads parameter does not affect result ---------------------------- ## + 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) + +}) + +## 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, + 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 + # 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), + 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, 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 + # 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], c("0" = 1001)) # all positions at zero depth + expect_true(all(vNoneSimple[-1L] == 0)) + expect_equal(sum(vNoneSimple), 1001) # region length +}) +