Chrcov

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")),

NAMESPACE

@@ -18,6 +18,7 @@ export(filterPositions)
 export(filterReads)
 export(filterReadsModBam)
 export(flattenReadLevelAssay)
+export(getBaseCoverageForBam)
 export(getReadLevelAssayNames)
 export(getReadLevelColDataNames)
 export(getReadNamesBySample)

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

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)
+}
+

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<std::vector<std::string>> 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<std::vector<std::string>> >::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}
 };