@@ -33,7 +33,8 @@ Imports:

     BiocParallel,

     readr,

     BSgenome,

-    Biostrings

+    Biostrings,

+    utils

 Suggests:

     testthat,

     bsseqData,

@@ -7,11 +7,11 @@ importFrom(BiocGenerics, "anyDuplicated", "cbind", "colnames",

            "combine", "density", "intersect", "lapply", "ncol",

            "nrow", "order", "paste", "pmax", "pmin", "rbind",

            "Reduce", "rep.int", "rownames", "sapply", "setdiff",

-           "strand", "strand<-", "union", "unique", "updateObject")

+           "strand", "strand<-", "union", "unique", "updateObject", "unstrand")

 importFrom(stats, "approxfun", "fisher.test", "ppoints",

            "predict", "preplot", "qchisq",

            "qnorm", "qqplot", "qunif", "cov2cor",

-           "plogis")

+           "plogis", "setNames")

 importFrom(graphics, "abline", "axis", "layout", "legend", "lines",

            "mtext", "par", "plot", "points", "polygon", "rect", "rug", "text")

 import(parallel)

@@ -27,15 +27,15 @@ importMethodsFrom(Biobase, "annotatedDataFrameFrom",

                   "pData", "pData<-",

                   "sampleNames", "sampleNames<-")

 importFrom(Biobase, "validMsg")

-importMethodsFrom(GenomeInfoDb, "seqlengths", "seqlengths<-", "seqinfo", "seqinfo<-",

-                  "seqnames", "seqnames<-", "seqlevels", "seqlevels<-")

+importMethodsFrom(GenomeInfoDb, "seqlengths", "seqlengths<-", "seqinfo",

+                  "seqinfo<-", "seqnames", "seqnames<-", "seqlevels",

+                  "seqlevels<-", "sortSeqlevels")

+importFrom(GenomeInfoDb, "Seqinfo")

 import(S4Vectors)

 importFrom(gtools, "combinations")

-# importFrom(data.table, "setnames", "setDT", "data.table", ":=", "setkey")

-# TODO: More careful use of imports to avoid clashes (e.g., S4Vectors::shift()

-#       and data.table::shift()).

-import(data.table)

-importFrom(R.utils, "isGzipped", "gunzip")

+# NOTE: data.table has some NAMESPACE clashes with functions in Bioconductor,

+#       e.g., shift(). If new ones are discovered, add them to this list.

+import(data.table, except = c(shift, first, second))

 import(limma)

 importFrom(permute, "shuffleSet", "how")

 import(DelayedArray)

@@ -43,12 +43,12 @@ import(BiocParallel)

 importFrom(readr, "cols", "cols_only", "col_character", "col_integer",

            "col_skip", "col_double", "col_factor", "read_tsv", "tokenizer_tsv")

 importFrom(Biostrings, "DNAString", "vmatchPattern", "reverseComplement")

+importFrom(utils, "read.delim")

 ##

 ## Exporting

 ##

-

 exportClasses("hasGRanges",

               "BSseq",

               "BSseqTstat",

@@ -226,7 +226,11 @@ BSmooth <- function(BSseq,

         se.coef_sink <- NULL

         sink_lock <- NULL

     } else if (BACKEND == "HDF5Array") {

-        coef_sink <- HDF5RealizationSink(

+        if (!requireNamespace("HDF5Array", quietly = TRUE)) {

+            stop("HDF5Array package required for HDF5Array backend",

+                 call. = FALSE)

+        }

+        coef_sink <- HDF5Array::HDF5RealizationSink(

             dim = dim(M),

             # NOTE: Never allow dimnames.

             dimnames = NULL,

@@ -237,7 +241,7 @@ BSmooth <- function(BSseq,

         sink_lock <- ipcid()

         on.exit(ipcremove(sink_lock), add = TRUE)

         if (keep.se) {

-            se.coef_sink <- HDF5RealizationSink(

+            se.coef_sink <- HDF5Array::HDF5RealizationSink(

                 dim = dim(M),

                 # NOTE: Never allow dimnames.

                 dimnames = NULL,

@@ -53,3 +53,78 @@

         stop("'type' = ", type, " is not supported")

     }

 }

+

+# Helper functions for setting up ArrayGrid instances --------------------------

+

+# NOTE: Copy of minfi:::colGrid()

+# TODO: Perhaps move this to DelayedMatrixStats?

+colGrid <- function(x) {

+    block_maxlen <- max(nrow(x),

+                        DelayedArray:::get_default_block_maxlength(type(x)))

+    spacings <- DelayedArray:::get_spacings_for_linear_capped_length_blocks(

+        refdim = dim(x),

+        block_maxlen = block_maxlen)

+    RegularArrayGrid(dim(x), spacings)

+}

+

+# NOTE: Copy of minfi:::rowGrid()

+# TODO: Perhaps move this to DelayedMatrixStats?

+rowGrid <- function(x) {

+    block_maxlen <- max(ncol(x),

+                        DelayedArray:::get_default_block_maxlength(type(x)))

+    spacings <- DelayedArray:::get_spacings_for_hypercube_capped_length_blocks(

+        refdim = dim(x),

+        block_maxlen = block_maxlen)

+    RegularArrayGrid(dim(x), spacings)

+}

+

+# Advanced block processing routines -------------------------------------------

+

+# NOTE: Copy of minfi:::blockApplyWithRealization()

+# TODO: Perhaps move this to DelayedMatrixStats?

+# NOTE: DelayedArray::blockApply() with the option to write the blocks to

+#       'sink'. Useful, for example, to apply a function across column-blocks

+#       of a DelayedMatrix, write these results to disk, and then wrap

+#       these in a DelayedMatrix.

+# TODO: See https://github.com/Bioconductor/DelayedArray/issues/10

+blockApplyWithRealization <- function(x, FUN, ..., sink = NULL, x_grid = NULL,

+                                      sink_grid = NULL, BPREDO = list(),

+                                      BPPARAM = bpparam()) {

+    FUN <- match.fun(FUN)

+

+    # Check conformable dots_grids and sinks_grids

+    x_grid <- DelayedArray:::.normarg_grid(x_grid, x)

+    sink_grid <- DelayedArray:::.normarg_grid(sink_grid, sink)

+    if (!identical(dim(x_grid), dim(sink_grid))) {

+        stop("non-conformable 'x_grid' and 'sink_grid'")

+    }

+

+    # Loop over blocks of `x` and write to `sink`

+    nblock <- length(x_grid)

+    bplapply(seq_len(nblock), function(b) {

+        if (DelayedArray:::get_verbose_block_processing()) {

+            message("Processing block ", b, "/", nblock, " ... ",

+                    appendLF = FALSE)

+        }

+        x_viewport <- x_grid[[b]]

+        sink_viewport <- sink_grid[[b]]

+        block <- DelayedArray:::extract_block(x, x_viewport)

+        if (!is.array(block)) {

+            block <- DelayedArray:::.as_array_or_matrix(block)

+        }

+        attr(block, "from_grid") <- x_grid

+        attr(block, "block_id") <- b

+        block_ans <- FUN(block, ...)

+        # NOTE: This is the only part different from DelayedArray::blockApply()

+        if (!is.null(sink)) {

+            write_block_to_sink(block_ans, sink, sink_viewport)

+            block_ans <- NULL

+        }

+        if (DelayedArray:::get_verbose_block_processing()) {

+            message("OK")

+        }

+    },

+    BPREDO = BPREDO,

+    BPPARAM = BPPARAM)

+}

+

@@ -231,11 +231,12 @@ setMethod("findOverlaps", c("FWGRanges", "FWGRanges"), .findOverlaps_FWGRanges)

             # strand since no longer valid.

             dt[strand == "-", start := start - 1L][, strand := NULL]

             # Aggregate counts at loci with the same 'seqnames' and 'start'.

-            dt <- dt[, .(M = sum(M), U = sum(U)), by = c("seqnames", "start")]

+            dt <- dt[,

+                     list(M = sum(M), U = sum(U)), by = c("seqnames", "start")]

         }

         # Identify loci with non-zero coverage then drop 'M' and 'U' as no

         # longer required.

-        dt <- dt[(M + U) > 0][, c("M", "U") := .(NULL, NULL)]

+        dt <- dt[(M + U) > 0][, c("M", "U") := list(NULL, NULL)]

     } else {

         dt <- .readBismarkAsDT(

             file = file,

@@ -42,7 +42,7 @@ setMethod(

                 dim = c(nrow(x), length(idx)),

                 dimnames = list(rownames(x), names(idx)),

                 type = "double")

-            x_grid <- minfi:::colGrid(x)

+            x_grid <- colGrid(x)

             sink_grid <- RegularArrayGrid(

                 refdim = dim(sink),

                 spacings = c(nrow(sink), ncol(sink) / length(x_grid)))

@@ -53,7 +53,7 @@ setMethod(

                 dimnames = list(names(idx), colnames(x)),

                 type = "double")

             on.exit(close(sink))

-            x_grid <- minfi:::rowGrid(x)

+            x_grid <- rowGrid(x)

             sink_grid <- RegularArrayGrid(

                 refdim = dim(sink),

                 spacings = c(nrow(sink) / length(x_grid), ncol(sink)))

@@ -62,7 +62,7 @@ setMethod(

         }

         # Loop over blocks of 'x' and write to 'sink'.

-        minfi:::blockApplyWithRealization(

+        blockApplyWithRealization(

             x = x,

             FUN = .collapseMatrixLike,

             idx = idx,

@@ -63,7 +63,7 @@

     on.exit(close(sink))

     # Set up ArrayGrid instance over columns of `sink`.

-    sink_grid <- minfi:::colGrid(sink)

+    sink_grid <- colGrid(sink)

     # Loop over column grid of 'sink', identify samples required for that

     # block, bring that data into memory, pass down to .combineList_matrix(),

@@ -184,7 +184,7 @@

         # strand since no longer valid.

         dt[strand == "-", start := start - 1L][, strand := NULL]

         # Aggregate counts at loci with the same 'seqnames' and 'start'.

-        dt <- dt[, .(M = sum(M), U = sum(U)), by = c("seqnames", "start")]

+        dt <- dt[, list(M = sum(M), U = sum(U)), by = c("seqnames", "start")]

     }

     # Construct FWGRanges ------------------------------------------------------

@@ -214,7 +214,7 @@

     M <- matrix(rep(0L, length(loci)), ncol = 1)

     Cov <- matrix(rep(0L, length(loci)), ncol = 1)

     M[subjectHits(ol)] <- dt[queryHits(ol), ][["M"]]

-    Cov[subjectHits(ol)] <- dt[queryHits(ol), .(Cov = (M + U))][["Cov"]]

+    Cov[subjectHits(ol)] <- dt[queryHits(ol), list(Cov = (M + U))][["Cov"]]

     # Return 'M' and 'Cov' or write them to the RealizationSink objects --------

@@ -241,9 +241,13 @@

         Cov_sink <- NULL

         sink_lock <- NULL

     } else if (BACKEND == "HDF5Array") {

+        if (!requireNamespace("HDF5Array", quietly = TRUE)) {

+            stop("HDF5Array package required for HDF5Array backend",

+                 call. = FALSE)

+        }

         # TODO: HDF5Array is only in suggests, so need to qualify the use of

         #       HDF5RealizationSink()

-        M_sink <- HDF5RealizationSink(

+        M_sink <- HDF5Array::HDF5RealizationSink(

             dim = c(ans_nrow, ans_ncol),

             # NOTE: Never allow dimnames.

             dimnames = NULL,

@@ -255,7 +259,7 @@

             # level = NULL,

             ...)

         on.exit(close(M_sink), add = TRUE)

-        Cov_sink <- HDF5RealizationSink(

+        Cov_sink <- HDF5Array::HDF5RealizationSink(

             dim = c(ans_nrow, ans_ncol),

             # NOTE: Never allow dimnames.

             dimnames = NULL,

@@ -294,7 +298,7 @@

     #       bar. Only SnowParam (and MulticoreParam by inheritance) have a

     #       bptasks<-() method.

     # TODO: Check that setting number of tasks doesn't affect things (e.g.,

-    #       the cost of transfering loci_dt to the workers may be substantial).

+    #       the cost of transfering loci to the workers may be substantial).

     if (is(BPPARAM, "SnowParam") && bpprogressbar(BPPARAM)) {

         bptasks(BPPARAM) <- length(grid)

     }

@@ -445,8 +449,7 @@ read.bismark <- function(files,

             #       backends. If the realization backend is NULL then an

             #       ordinary matrix is returned rather than a matrix-backed

             #       DelayedMatrix.

-            stop("The '", realization_backend, "' realization backend is ",

-                 "not supported.\n",

+            stop("The '", BACKEND, "' realization backend is not supported.\n",

                  "See help(\"BSmooth\") for details.",

                  call. = FALSE)

         }

@@ -70,7 +70,7 @@ BSmooth(BSseq,

   The choice of realization backend is controlled by the \code{BACKEND}

   argument, which defaults to the current value of

-  \code{DelayedArray::\link[DelayedArray]{getRealizationBackend()}}.

+  \code{DelayedArray::\link[DelayedArray]{getRealizationBackend}()}.

   \code{BSmooth} supports the following realization backends:

@@ -82,9 +82,8 @@ BSseq(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

   #

   library(HDF5Array)

-  # See ?HDF5Array::writeHDF5Array for details

-  hdf5_M <- writeHDF5Array(M)

-  hdf5_Cov <- writeHDF5Array(Cov)

+  hdf5_M <- realize(M, "HDF5Array")

+  hdf5_Cov <- realize(Cov, "HDF5Array")

   hdf5_BS1 <- BSseq(chr = c("chr1", "chr2", "chr1"),

                     pos = c(1, 2, 3),

                     M = hdf5_M,

@@ -8,7 +8,8 @@

 }

 \usage{

 getCoverage(BSseq, regions = NULL, type = c("Cov", "M"),

-  what = c("perBase", "perRegionAverage", "perRegionTotal"))

+  what = c("perBase", "perRegionAverage", "perRegionTotal"),

+  withDimnames = TRUE)

 }

 \arguments{

   \item{BSseq}{An object of class \code{BSseq}.}

@@ -17,6 +18,10 @@ getCoverage(BSseq, regions = NULL, type = c("Cov", "M"),

   \item{type}{This returns either coverage or the total

     evidence for methylation at the loci.}

   \item{what}{The type of return object, see details.}

+  \item{withDimnames}{A \code{logical(1)}, indicating whether dimnames should

+    be applied to extracted coverage elements. Setting

+    \code{withDimnames = FALSE} increases the speed and memory efficiency with

+    which coverage is extracted.}

 }

 \value{

   \strong{NOTE:} The return type of \code{getCoverage} varies depending on its

@@ -8,7 +8,8 @@

 }

 \usage{

 getMeth(BSseq, regions = NULL, type = c("smooth", "raw"),

-  what = c("perBase", "perRegion"), confint = FALSE, alpha = 0.95)

+  what = c("perBase", "perRegion"), confint = FALSE, alpha = 0.95,

+  withDimnames = TRUE)

 }

 \arguments{

   \item{BSseq}{An object of class \code{BSseq}.}

@@ -20,6 +21,10 @@ getMeth(BSseq, regions = NULL, type = c("smooth", "raw"),

     methylation estimates (see below).  This is only supported if

     \code{what} is equal to \code{perBase}.}

   \item{alpha}{alpha value for the confidence interval.}

+  \item{withDimnames}{A \code{logical(1)}, indicating whether dimnames should

+    be applied to extracted coverage elements. Setting

+    \code{withDimnames = FALSE} increases the speed and memory efficiency with

+    which coverage is extracted.}

 }

 \note{

   A \code{BSseq} object needs to be smoothed by the function

@@ -136,7 +136,7 @@

 \section{Realization backends}{

   The \code{read.bismark()} function creates a \linkS4class{BSseq} object with two assays, \code{M} and \code{Cov}.

   The choice of \emph{realization backend} controls whether these assays are stored in-memory as an ordinary \link[base]{matrix} or on-disk as a \linkS4class{HDF5Array}, for example.

-  The choice of realization backend is controlled by the \code{BACKEND} argument, which defaults to the current value of \code{DelayedArray::\link[DelayedArray]{getRealizationBackend()}}.

+  The choice of realization backend is controlled by the \code{BACKEND} argument, which defaults to the current value of \code{DelayedArray::\link[DelayedArray]{getRealizationBackend}()}.

   \code{read.bismark()} supports the following realization backends:

@@ -1,114 +1,5 @@

 context("combine()")

-# TODO: .subassignRowsDelayedMatrix() is no longer needed in bsseq; test can be

-#       removed.

-# test_that(".subassignRowsDelayedMatrix()", {

-#     nrow <- 1000L

-#     ncol <- 10L

-#     x <- realize(matrix(seq_len(nrow * ncol), ncol = ncol), "HDF5Array")

-#     x_i <- seq(1L, 2L * nrow, 2L)

-#     y <- realize(matrix(seq(-1L, -nrow * ncol, -1L), ncol = 10), "HDF5Array")

-#     y_i <- seq(2L, nrow, 2L)

-#

-#     z1 <- bsseq:::.subassignRowsDelayedMatrix(x = x,

-#                                               i = x_i,

-#                                               nrow = 2L * nrow,

-#                                               fill = NA_integer_,

-#                                               BACKEND = NULL)

-#     z2 <- bsseq:::.subassignRowsDelayedMatrix(x = x,

-#                                               i = x_i,

-#                                               nrow = 2L * nrow,

-#                                               fill = NA_integer_,

-#                                               BACKEND = "HDF5Array",

-#                                               by_row = FALSE)

-#     z3 <- bsseq:::.subassignRowsDelayedMatrix(x = x,

-#                                               i = x_i,

-#                                               nrow = 2L * nrow,

-#                                               fill = NA_integer_,

-#                                               BACKEND = "HDF5Array",

-#                                               by_row = TRUE)

-#     expect_identical(as.array(z1), as.array(z2))

-#     expect_identical(as.array(z1), as.array(z3))

-# })

-

-# TODO: .combineListOfDelayedMatrixObjects() is no longer needed in bsseq; test

-#        can be removed.

-# test_that(".combineListOfDelayedMatrixObjects()", {

-#     nrow <- 10

-#     ncol <- 4

-#     x <- matrix(seq_len(nrow),

-#                 ncol = ncol / 2,

-#                 dimnames = list(NULL, letters[1:2]))

-#     y <- matrix(100L + seq_len(nrow),

-#                 ncol = ncol / 2,

-#                 dimnames = list(NULL, letters[3:4]))

-#     x_i <- seq(1, nrow, ncol / 2)

-#     y_i <- seq(2, nrow, ncol / 2)

-#     fill <- NA_integer_

-#

-#     # The expected output

-#     z <- matrix(fill,

-#                 nrow = nrow,

-#                 ncol = ncol,

-#                 dimnames = list(NULL, letters[seq_len(ncol)]))

-#     # NOTE: as.array(x) is a no-op if x is a matrix and realises a

-#     #       DelayedMtrix in memory

-#     z[x_i, seq(1, ncol(x))] <- x

-#     z[y_i, seq(ncol(x) + 1, ncol(x) + ncol(y))] <- y

-#

-#     # # Test with in-memory DelayedMatrix objects

-#     X <- bsseq:::.DelayedMatrix(x)

-#     Y <- bsseq:::.DelayedMatrix(y)

-#

-#     Z <- bsseq:::.combineListOfDelayedMatrixObjects(

-#         X = list(X, Y),

-#         I = list(x_i, y_i),

-#         nrow = nrow,

-#         ncol = ncol,

-#         dimnames = list(NULL, c(colnames(X), colnames(Y))),

-#         fill = fill,

-#         BACKEND = NULL)

-#     expect_identical(z, as.array(Z))

-#     expect_true(!bsseq:::.isHDF5ArrayBacked(Z))

-#

-#

-#     # Test with HDF5Array-backed DelayedMatrix objects

-#     hdf5_X <- realize(X, BACKEND = "HDF5Array")

-#     hdf5_Y <- realize(Y, BACKEND = "HDF5Array")

-#

-#     hdf5_Z <- bsseq:::.combineListOfDelayedMatrixObjects(

-#         X = list(hdf5_X, hdf5_Y),

-#         I = list(x_i, y_i),

-#         nrow = nrow,

-#         ncol = ncol,

-#         dimnames = list(NULL, c(colnames(hdf5_X), colnames(hdf5_Y))),

-#         fill = fill,

-#         BACKEND = "HDF5Array")

-#     expect_identical(z, as.array(hdf5_Z))

-#     expect_true(bsseq:::.isHDF5ArrayBacked(hdf5_Z))

-# })

-

-checkBSseqAssaysIdentical <- function(x, y) {

-    stopifnot(is(x, "BSseq") && is(y, "BSseq"))

-    assay_names <- c("M", "Cov", "coef", "se.coef")

-    check_identical <- vapply(assay_names, function(an) {

-        if (!is.null(getBSseq(x, an))) {

-            identical(as.array(getBSseq(x, an)), as.array(getBSseq(y, an)))

-        } else {

-            identical(getBSseq(x, an), getBSseq(y, an))

-        }

-    }, logical(1L))

-    checkTrue(all(check_identical))

-}

-

-checkBSseqIdentical <- function(x, y) {

-    checkTrue(identical(rowRanges(x), rowRanges(y)) &&

-                  identical(getBSseq(x, "trans"), getBSseq(y, "trans")) &&

-                  identical(getBSseq(x, "parameters"),

-                            getBSseq(y, "parameters")) &&

-                  checkBSseqAssaysIdentical(x, y))

-}

-

 test_that("combine()", {

     bsseq_fit <- BSmooth(bsseq_test)

     BSSEQ_TEST <- realize(bsseq_test, "HDF5Array")

@@ -168,6 +59,7 @@ test_that("combineList()", {

     Z <- combineList(BSSEQ_FIT[, 1], BSSEQ_FIT[, 2])

     expect_equivalent_SE(Z, BSSEQ_FIT)

 })

+

 test_that(

     "Test bug fix reported in https://github.com/hansenlab/bsseq/pull/54/", {

         M <- matrix(0:8, 3, 3)