@@ -3,3 +3,4 @@

 .RData$

 auto/

 *~

+.Rhistory

@@ -1,19 +1,19 @@

 Package: bsseq

-Version: 1.9.2

+Version: 1.11.1

+Encoding: UTF-8

 Title: Analyze, manage and store bisulfite sequencing data

 Description: A collection of tools for analyzing and visualizing bisulfite

     sequencing data.

 Authors@R: c(person(c("Kasper", "Daniel"), "Hansen", role = c("aut", "cre"),

     email = "kasperdanielhansen@gmail.com"),

-    person("Peter", "Hickey", role = "ctb", email = "peter.hickey@gmail.com"))

+    person("Peter", "Hickey", role = "aut", email = "peter.hickey@gmail.com"))

 Depends:

     R (>= 3.3),

     methods,

     BiocGenerics,

     GenomicRanges (>= 1.23.7),

     SummarizedExperiment (>= 0.1.1),

-    parallel,

-    limma

+    parallel

 Imports:

     IRanges (>= 2.5.17),

     GenomeInfoDb,

@@ -27,19 +27,22 @@ Imports:

     S4Vectors,

     R.utils (>= 2.0.0),

     matrixStats (>= 0.50.0),

-    permute

+    permute,

+    limma,

+    DelayedArray,

+    HDF5Array

 Suggests:

     RUnit,

     bsseqData,

     BiocStyle,

-    knitr

+    knitr,

 Collate:

+    utils.R

     hasGRanges.R

     BSseq_class.R

     BSseqTstat_class.R

     BSseq_utils.R

     combine.R

-    utils.R

     read.bsmooth.R

     read.bismark.R

     BSmooth.R

@@ -52,6 +55,9 @@ Collate:

     BSmooth.fstat.R

     BSseqStat_class.R

     getStats.R

+    hdf5_utils.R

+    combine_utils.R

+    DelayedArray_utils.R

 License: Artistic-2.0

 VignetteBuilder: knitr

 URL: https://github.com/kasperdanielhansen/bsseq

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

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

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

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

-           "qnorm", "qqplot", "qunif", "cov2cor")

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

+           "plogis")

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

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

 import(parallel)

@@ -33,6 +34,9 @@ importFrom(data.table, "fread", "setnames")

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

 import(limma)

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

+importClassesFrom(DelayedArray, "DelayedArray", "DelayedMatrix")

+importFrom(DelayedArray, "DelayedArray", "plogis", "pmin2", "pmax2",

+           "setRealizeBackend", "type")

 ##

 ## Exporting

@@ -17,7 +17,7 @@ makeClusters <- function(hasGRanges, maxGap = 10^8) {

     names(clusterIdx) <- NULL

     clusterIdx

 }

-    

+

 BSmooth <- function(BSseq, ns = 70, h = 1000, maxGap = 10^8, parallelBy = c("sample", "chromosome"),

                     mc.preschedule = FALSE, mc.cores = 1, keep.se = FALSE, verbose = TRUE) {

@@ -29,8 +29,8 @@ BSmooth <- function(BSseq, ns = 70, h = 1000, maxGap = 10^8, parallelBy = c("sam

         if(verbose >= 2)

             cat(sprintf("[BSmooth]   smoothing start: sample:%s, chr:%s, nLoci:%s\n",

                         this_sample_chr[1], this_sample_chr[2], length(idxes)))

-        Cov <- getCoverage(BSseq, type = "Cov")[idxes, sampleIdx]

-        M <- getCoverage(BSseq, type = "M")[idxes, sampleIdx]

+        Cov <- getCoverage(BSseq, type = "Cov")[idxes, sampleIdx, drop = FALSE]

+        M <- getCoverage(BSseq, type = "M")[idxes, sampleIdx, drop = FALSE]

         pos <- start(BSseq)[idxes]

         stopifnot(all(diff(pos) > 0))

         wh <- which(Cov != 0)

@@ -45,8 +45,9 @@ BSmooth <- function(BSseq, ns = 70, h = 1000, maxGap = 10^8, parallelBy = c("sam

                         trans = NULL, h = h, nn = nn))

         }

         sdata <- data.frame(pos = pos[wh],

-                            M = pmin(pmax(M[wh], 0.01), Cov[wh] - 0.01),

-                            Cov = Cov[wh])

+                            M = as.array(pmin2(pmax2(M[wh, ], 0.01),

+                                               Cov[wh, ] - 0.01)),

+                            Cov = Cov[wh, ])

         fit <- locfit(M ~ lp(pos, nn = nn, h = h), data = sdata,

                       weights = Cov, family = "binomial", maxk = 10000)

         pp <- preplot(fit, where = "data", band = "local",

@@ -119,11 +120,15 @@ BSmooth <- function(BSseq, ns = 70, h = 1000, maxGap = 10^8, parallelBy = c("sam

         coef <- do.call(rbind, lapply(out, function(xx) xx$coef))

         se.coef <- do.call(rbind, lapply(out, function(xx) xx$se.coef))

     })

+    coef <- .DelayedMatrix(coef)

+    if (!is.null(se.coef)) {

+        se.coef <- .DelayedMatrix(se.coef)

+    }

     ptime.outer2 <- proc.time()

     stime.outer <- (ptime.outer2 - ptime.outer1)[3]

     if(verbose)

         cat(sprintf("[BSmooth] smoothing done in %.1f sec\n", stime.outer))

-    

+

     rownames(coef) <- NULL

     colnames(coef) <- sampleNames(BSseq)

     if(!is.null(se.coef)) {

@@ -135,16 +140,8 @@ BSmooth <- function(BSseq, ns = 70, h = 1000, maxGap = 10^8, parallelBy = c("sam

         assay(BSseq, "coef") <- coef

     if(!is.null(se.coef))

         assay(BSseq, "se.coef") <- se.coef

-    mytrans <- function(x) {

-        y <- x

-        ix <- which(x < 0)

-        ix2 <- which(x > 0)

-        y[ix] <- exp(x[ix])/(1 + exp(x[ix]))

-        y[ix2] <- 1/(1 + exp(-x[ix2]))

-        y

-    }

-    environment(mytrans) <- baseenv()

-    BSseq@trans <- mytrans

+    BSseq@trans <- plogis

+

     parameters <- list(smoothText = sprintf("BSmooth (ns = %d, h = %d, maxGap = %d)", ns, h, maxGap),

                        ns = ns, h = h, maxGap = maxGap)

     BSseq@parameters <- parameters

@@ -9,6 +9,7 @@ BSmooth.fstat <- function(BSseq, design, contrasts, verbose = TRUE){

     ptime1 <- proc.time()

     allPs <- getMeth(BSseq, type = "smooth", what = "perBase",

                      confint = FALSE)

+    allPs <- as.array(allPs)

     fit <- lmFit(allPs, design)

     fitC <- contrasts.fit(fit, contrasts)

     ## Need

@@ -17,7 +18,7 @@ BSmooth.fstat <- function(BSseq, design, contrasts, verbose = TRUE){

     ##   tstats <- fitC$coefficients / fitC$stdev.unscaled / fitC$sigma

     ##   rawSds <- fitC$sigma

     ##   cor.coefficients <- cov2cor(fitC$cov.coefficients)

-    rawSds <- fitC$sigma

+    rawSds <- as.matrix(fitC$sigma)

     cor.coefficients <- cov2cor(fitC$cov.coefficients)

     rawTstats <- fitC$coefficients / fitC$stdev.unscaled / fitC$sigma

     names(dimnames(rawTstats)) <- NULL

@@ -31,7 +32,8 @@ BSmooth.fstat <- function(BSseq, design, contrasts, verbose = TRUE){

                   cor.coefficients = cor.coefficients,

                   rawTstats = rawTstats)

     out <- BSseqStat(gr = granges(BSseq),

-                     stats = stats, parameters = parameters)

+                     stats = stats,

+                     parameters = parameters)

     out

 }

@@ -57,8 +59,12 @@ smoothSds <- function(BSseqStat, k = 101, qSd = 0.75, mc.cores = 1,

     if(verbose) cat(sprintf("done in %.1f sec\n", stime))

     smoothSds <- do.call("c",

                          mclapply(clusterIdx, function(idx) {

-                             smoothSd(getStats(BSseqStat, what = "rawSds")[idx], k = k, qSd = qSd)

+                             rawSds <- getStats(BSseqStat,

+                                                what = "rawSds")[idx, ]

+                             rawSds <- as.array(rawSds)

+                             smoothSd(rawSds, k = k, qSd = qSd)

                          }, mc.cores = mc.cores))

+    smoothSds <- .DelayedMatrix(as.matrix(smoothSds))

     if("smoothSds" %in% names(getStats(BSseqStat)))

         BSseqStat@stats[["smoothSds"]] <- smoothSds

     else

@@ -68,21 +74,30 @@ smoothSds <- function(BSseqStat, k = 101, qSd = 0.75, mc.cores = 1,

 # Quieten R CMD check

 globalVariables("tstat")

+

+# NOTE: Realises in memory a matrix with nrow = length(BSseqStat) and

+#       ncol = length(coef)

 computeStat <- function(BSseqStat, coef = NULL) {

     stopifnot(is(BSseqStat, "BSseqStat"))

     if(is.null(coef)) {

         coef <- 1:ncol(getStats(BSseqStat, what = "rawTstats"))

     }

-    tstats <- getStats(BSseqStat, what = "rawTstats")[, coef, drop = FALSE]

-    tstats <- tstats * getStats(BSseqStat, what = "rawSds") /

+    raw_tstats <- getStats(BSseqStat, what = "rawTstats")[, coef, drop = FALSE]

+    scaled_sds <- getStats(BSseqStat, what = "rawSds") /

         getStats(BSseqStat, what = "smoothSds")

+    scaled_sds_matrix <- do.call(cbind, replicate(ncol(raw_tstats), scaled_sds))

+    tstats <- raw_tstats * scaled_sds_matrix

     if(length(coef) > 1) {

-        cor.coefficients <- getStats(BSseqStat, what = "cor.coefficients")[coef,coef]

-        stat <- as.numeric(classifyTestsF(tstats, cor.coefficients,

-                                          fstat.only = TRUE))

+        cor.coefficients <- getStats(BSseqStat,

+                                     what = "cor.coefficients")[coef,coef]

+        # NOTE: classifyTestsF() calls as.matrix(tstats) and so realises this

+        #       array

+        stat <- .DelayedMatrix(as.matrix(classifyTestsF(tstats,

+                                                        cor.coefficients,

+                                                        fstat.only = TRUE)))

         stat.type <- "fstat"

     } else {

-        stat <- as.numeric(tstats)

+        stat <- .DelayedMatrix(tstats)

         stat.type <- "tstat"

     }

     if("stat" %in% names(getStats(BSseqStat))) {

@@ -114,9 +129,9 @@ localCorrectStat <- function(BSseqStat, threshold = c(-15,15), mc.cores = 1, ver

         xx.reg <- seq(from = min(xx), to = max(xx), by = 2000)

         yy.reg <- tstat.function(xx.reg)

         fit <- locfit(yy.reg ~ lp(xx.reg, h = 25000, deg = 2, nn = 0),

-                      family = "huber", maxk = 50000) 

+                      family = "huber", maxk = 50000)

         correction <- predict(fit, newdata = data.frame(xx.reg = xx))

-        yy - correction 

+        yy - correction

     }

     maxGap <- BSseqStat$parameters$maxGap

     if(verbose) cat("[BSmooth.tstat] preprocessing ... ")

@@ -20,9 +20,9 @@ BSmooth.tstat <- function(BSseq, group1, group2, estimate.var = c("same", "paire

         xx.reg <- seq(from = min(xx), to = max(xx), by = 2000)

         yy.reg <- tstat.function(xx.reg)

         fit <- locfit(yy.reg ~ lp(xx.reg, h = 25000, deg = 2, nn = 0),

-                      family = "huber", maxk = 50000) 

+                      family = "huber", maxk = 50000)

         correction <- predict(fit, newdata = data.frame(xx.reg = xx))

-        yy - correction 

+        yy - correction

     }

     estimate.var <- match.arg(estimate.var)

@@ -38,7 +38,7 @@ BSmooth.tstat <- function(BSseq, group1, group2, estimate.var = c("same", "paire

     if(is.character(group2)) {

         stopifnot(all(group2 %in% sampleNames(BSseq)))

         group2 <- match(group2, sampleNames(BSseq))

-    }    

+    }

     if(is.numeric(group2)) {

         stopifnot(min(group2) >= 1 & max(group2) <= ncol(BSseq))

     } else stop("problems with argument 'group2'")

@@ -48,22 +48,29 @@ BSmooth.tstat <- function(BSseq, group1, group2, estimate.var = c("same", "paire

     stopifnot(length(group1) + length(group2) >= 3)

     if(estimate.var == "paired")

         stopifnot(length(group1) == length(group2))

-    

-    if(any(rowSums(getCoverage(BSseq)[, c(group1, group2)]) == 0))

-        warning("Computing t-statistics at locations where there is no data; consider subsetting the 'BSseq' object first")

-    

+

+    if (!.isHDF5ArrayBacked(getCoverage(BSseq))) {

+        ## TODO: This check can be super expensive if Cov is a HDF5Matrix, so

+        ##       we skip it for the time being. Really noticeable when Cov

+        ##       needs row subsetted/reordered (i.e. has Cov@index[[1]] exists)

+        ##

+        if(any(rowSums(getCoverage(BSseq)[, c(group1, group2)]) == 0)) {

+            warning("Computing t-statistics at locations where there is no data; consider subsetting the 'BSseq' object first")

+        }

+    }

+

     if(is.null(maxGap))

         maxGap <- BSseq@parameters$maxGap

     if(is.null(maxGap))

         stop("need to set argument 'maxGap'")

-    

+

     if(verbose) cat("[BSmooth.tstat] preprocessing ... ")

     ptime1 <- proc.time()

     clusterIdx <- makeClusters(BSseq, maxGap = maxGap)

     ptime2 <- proc.time()

     stime <- (ptime2 - ptime1)[3]

     if(verbose) cat(sprintf("done in %.1f sec\n", stime))

-        

+

     if(verbose) cat("[BSmooth.tstat] computing stats within groups ... ")

     ptime1 <- proc.time()

     allPs <- getMeth(BSseq, type = "smooth", what = "perBase",

@@ -73,12 +80,12 @@ BSmooth.tstat <- function(BSseq, group1, group2, estimate.var = c("same", "paire

     ptime2 <- proc.time()

     stime <- (ptime2 - ptime1)[3]

     if(verbose) cat(sprintf("done in %.1f sec\n", stime))

-    

+

     if(verbose) cat("[BSmooth.tstat] computing stats across groups ... ")

     ptime1 <- proc.time()

     switch(estimate.var,

            "group2" = {

-               rawSds <- rowSds(allPs, cols = group2, na.rm = TRUE)

+               rawSds <- .rowSds(x = allPs, cols = group2, na.rm = TRUE)

                smoothSds <- do.call(c, mclapply(clusterIdx, function(idx) {

                    smoothSd(rawSds[idx], k = k)

                }, mc.cores = mc.cores))

@@ -86,9 +93,9 @@ BSmooth.tstat <- function(BSseq, group1, group2, estimate.var = c("same", "paire

                tstat.sd <- smoothSds * scale

            },

            "same" = {

-               rawSds <- sqrt( ((length(group1) - 1) * rowVars(allPs, cols = group1) +

-                                (length(group2) - 1) * rowVars(allPs, cols = group2)) /

-                              (length(group1) + length(group2) - 2))

+               rawSds <- sqrt( ((length(group1) - 1) * .rowVars(x = allPs, cols = group1) +

+                                    (length(group2) - 1) * .rowVars(x = allPs, cols = group2)) /

+                                   (length(group1) + length(group2) - 2))

                smoothSds <- do.call(c, mclapply(clusterIdx, function(idx) {

                    smoothSd(rawSds[idx], k = k)

                }, mc.cores = mc.cores))

@@ -96,7 +103,7 @@ BSmooth.tstat <- function(BSseq, group1, group2, estimate.var = c("same", "paire

                tstat.sd <- smoothSds * scale

            },

            "paired" = {

-               rawSds <- rowSds(allPs[, group1, drop = FALSE] - allPs[, group2, drop = FALSE])

+               rawSds <- .rowSds(x = allPs[, group1, drop = FALSE] - allPs[, group2, drop = FALSE])

                smoothSds <- do.call(c, mclapply(clusterIdx, function(idx) {

                    smoothSd(rawSds[idx], k = k)

                }, mc.cores = mc.cores))

@@ -113,24 +120,24 @@ BSmooth.tstat <- function(BSseq, group1, group2, estimate.var = c("same", "paire

     ptime2 <- proc.time()

     stime <- (ptime2 - ptime1)[3]

     if(verbose) cat(sprintf("done in %.1f sec\n", stime))

-    

+

     if(local.correct) {

         stats <- cbind(rawSds, tstat.sd, group2.means, group1.means,

                        tstat, tstat.corrected)

         colnames(stats) <- c("rawSds", "tstat.sd",

                              "group2.means", "group1.means", "tstat",

                              "tstat.corrected")

- 

+

     } else {

         stats <- cbind(rawSds, tstat.sd, group2.means, group1.means,

                        tstat)

         colnames(stats) <- c("rawSds", "tstat.sd",

                              "group2.means", "group1.means", "tstat")

     }

-    

+

     parameters <- c(BSseq@parameters,

                     list(tstatText = sprintf("BSmooth.tstat (local.correct = %s, maxGap = %d)",

-                         local.correct, maxGap),

+                                             local.correct, maxGap),

                          group1 = group1, group2 = group2, k = k, qSd = qSd,

                          local.correct = local.correct, maxGap = maxGap))

     out <- BSseqTstat(gr = granges(BSseq), stats = stats, parameters = parameters)

@@ -1,4 +1,4 @@

-setClass("BSseqStat", contains = "hasGRanges", 

+setClass("BSseqStat", contains = "hasGRanges",

          representation(stats = "list",

                         parameters = "list")

          )

@@ -8,17 +8,12 @@ setValidity("BSseqStat", function(object) {

     if(is.null(names(object@stats)) || any(names(object@stats) == "") ||

        anyDuplicated(names(object@stats)))

         msg <- validMsg(msg, "the 'stats' list needs to be named with unique names.")

-    for(name in c("rawSds", "smoothsSds", "stat")) {

-        if(name %in% names(object@stats) && !is.vector(object@stats[[name]]))

-            msg <- validMsg(msg, sprintf("component '%s' of slot 'stats' have to be a vector", name))

-        if(name %in% names(object@stats) && length(object@stats[[name]]) != length(object@gr))

-            msg <- validMsg(msg, sprintf("component '%s' of slot 'stats' have to have the same length as slot 'gr'", name))

-    }

-    for(name in c("rawTstats")) {

-        if(name %in% names(object@stats) && !is.matrix(object@stats[[name]]))

-            msg <- validMsg(msg, sprintf("component '%s' of slot 'stats' have to be a matrix", name))

-        if(name %in% names(object@stats) && nrow(object@stats[[name]]) != length(object@gr))

-            msg <- validMsg(msg, sprintf("component '%s' of slot 'stats' have to have the same number of rows as slot 'gr' is long", name))

+    for(name in c("rawSds", "smoothsSds", "stat", "rawTstats")) {

+        if(name %in% names(object@stats) &&

+           !is(object@stats[[name]], "DelayedMatrix"))

+            msg <- validMsg(msg, sprintf("component '%s' of slot 'stats' has to be a DelayedMatrix object", name))

+        if(name %in% names(object@stats) && isTRUE(nrow(object@stats[[name]]) != length(object@gr)))

+            msg <- validMsg(msg, sprintf("component '%s' of slot 'stats' has to have the same number of rows as slot 'gr' is long", name))

     }

     if(is.null(msg)) TRUE else msg

 })

@@ -29,8 +24,17 @@ setMethod("show", signature(object = "BSseqStat"),

               cat(" ", length(object), "methylation loci\n")

               cat("based on smoothed data:\n")

               cat(" ", object@parameters$smoothText, "\n")

-              cat("with parameters\n")

-              cat(" ", object@parameters$tstatText, "\n")

+              not_delayed_matrices <- c("cor.coefficients", "stat.type")

+              delayed_matrices <- setdiff(names(object@stats),

+                                          not_delayed_matrices)

+              is_HDF5Array_backed <- vapply(object@stats[delayed_matrices],

+                                            .isHDF5ArrayBacked,

+                                            logical(1L))

+              if (any(is_HDF5Array_backed)) {

+                  cat("'stats' slot is HDF5Array-backed\n")

+              } else {

+                  cat("'stats' slot is in-memory\n")

+              }

           })

 setMethod("[", "BSseqStat", function(x, i, ...) {

@@ -42,14 +46,11 @@ setMethod("[", "BSseqStat", function(x, i, ...) {

     statnames <- names(x@stats)

     names(statnames) <- statnames

     x@stats <- lapply(statnames, function(nam) {

-        if(nam %in% c("rawTstats", "modelCoefficients")) {

-            stopifnot(is.matrix(x@stats[[nam]]))

+        if(nam %in% c("rawTstats", "modelCoefficients", "rawSds", "smoothSds",

+                      "stat")) {

+            stopifnot(is(x@stats[[nam]], "DelayedMatrix"))

             return(x@stats[[nam]][i,,drop=FALSE])

         }

-        if(nam %in% c("rawSds", "smoothSds", "stat")) {

-            stopifnot(is.vector(x@stats[[nam]]))

-            return(x@stats[[nam]][i])

-        }

         x@stats[[nam]]

     })

     x

@@ -58,11 +59,26 @@ setMethod("[", "BSseqStat", function(x, i, ...) {

 BSseqStat <- function(gr = NULL, stats = NULL, parameters = NULL) {

     out <- new("BSseqStat")

     out@gr <- gr

+    not_delayed_matrices <- c("cor.coefficients", "stat.type")

+    delayed_matrices <- setdiff(names(stats), not_delayed_matrices)

+    stats[delayed_matrices] <- endoapply(stats[delayed_matrices],

+                                         .DelayedMatrix)

     out@stats <- stats

     out@parameters <- parameters

     out

 }

+setMethod("updateObject", "BSseqStat",

+          function(object, ...) {

+              not_delayed_matrices <- c("cor.coefficients", "stat.type")

+              delayed_matrices <- setdiff(names(stats), not_delayed_matrices)

+              stats <- object@stats

+              stats[delayed_matrices] <- endoapply(stats[delayed_matrices],

+                                                   .DelayedMatrix)

+              object@stats <- stats

+              object

+          }

+)

 ## summary.BSseqStat <- function(object, ...) {

 ##     quant <- quantile(getStats(object)[, "tstat.corrected"],

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

-setClass("BSseqTstat", contains = "hasGRanges", 

-         representation(stats = "matrix",

+setClass("BSseqTstat", contains = "hasGRanges",

+         representation(stats = "DelayedMatrix",

                         parameters = "list")

          )

 setValidity("BSseqTstat", function(object) {

@@ -17,6 +17,11 @@ setMethod("show", signature(object = "BSseqTstat"),

               cat(" ", object@parameters$smoothText, "\n")

               cat("with parameters\n")

               cat(" ", object@parameters$tstatText, "\n")

+              if (.isHDF5ArrayBacked(object@stats)) {

+                  cat("'stats' slot is HDF5Array-backed\n")

+              } else {

+                  cat("'stats' slot is in-memory\n")

+              }

           })

 setMethod("[", "BSseqTstat", function(x, i, ...) {

@@ -32,14 +37,14 @@ setMethod("[", "BSseqTstat", function(x, i, ...) {

 BSseqTstat <- function(gr = NULL, stats = NULL, parameters = NULL) {

     out <- new("BSseqTstat")

     out@gr <- gr

-    out@stats <- stats

+    out@stats <- .DelayedMatrix(stats)

     out@parameters <- parameters

     out

 }

 summary.BSseqTstat <- function(object, ...) {

-    quant <- quantile(getStats(object)[, "tstat.corrected"],

-                      prob = c(0.0001, 0.001, 0.01, 0.5, 0.99, 0.999, 0.9999))

+    quant <- .quantile(getStats(object)[, "tstat.corrected"],

+                       prob = c(0.0001, 0.001, 0.01, 0.5, 0.99, 0.999, 0.9999))

     quant <- t(t(quant))

     colnames(quant) <- "quantiles"

     out <- list(quantiles = quant)

@@ -53,9 +58,11 @@ print.summary.BSseqTstat <- function(x, ...) {

 plot.BSseqTstat <- function(x, y, ...) {

     tstat <- getStats(x)[, "tstat"]

+    tstat <- as.array(tstat)

     plot(density(tstat), xlim = c(-10,10), col = "blue", main = "")

     if("tstat.corrected" %in% colnames(getStats(x))) {

         tstat.cor <- getStats(x)[, "tstat.corrected"]

+        tstat.cor <- as.array(tstat.cor)

         lines(density(tstat.cor), col = "black")

         legend("topleft", legend = c("uncorrected", "corrected"), lty = c(1,1),

                col = c("blue", "black"))

@@ -65,3 +72,11 @@ plot.BSseqTstat <- function(x, y, ...) {

     }

 }

+setMethod("updateObject", "BSseqTstat",

+          function(object, ...) {

+              stats <- object@stats

+              stats <- .DelayedMatrix(stats)

+              object@stats <- stats

+              object

+          }

+)

@@ -1,22 +1,30 @@

-setClass("BSseq", contains = "RangedSummarizedExperiment", 

+setClass("BSseq", contains = "RangedSummarizedExperiment",

          representation(trans = "function",

                         parameters = "list"))

 setValidity("BSseq", function(object) {

     msg <- validMsg(NULL, .checkAssayNames(object, c("Cov", "M")))

+    msg <- validMsg(msg, .checkAssayClasses(object,

+                                            c("Cov", "M", "coef", "se.coef")))

     if(class(rowRanges(object)) != "GRanges")

         msg <- validMsg(msg, sprintf("object of class '%s' needs to have a 'GRanges' in slot 'rowRanges'", class(object)))

     ## benchmarking shows that min(assay()) < 0 is faster than any(assay() < 0) if it is false

     if(is.null(colnames(object)))

         msg <- validMsg(msg, "colnames (aka sampleNames) need to be set")

-    if(min(assay(object, "M")) < 0)

-        msg <- validMsg(msg, "the 'M' assay has negative entries")

-    if(min(assay(object, "Cov")) < 0)

-        msg <- validMsg(msg, "the 'Cov' assay has negative entries")

-    if(max(assay(object, "M") - assay(object, "Cov")) > 0.5)

-        msg <- validMsg(msg, "the 'M' assay has at least one entry bigger than the 'Cov' assay")

-    if(!is.null(rownames(assay(object, "M"))) ||

-       !is.null(rownames(assay(object, "Cov"))) ||

+    M <- assay(object, "M", withDimnames = FALSE)

+    Cov <- assay(object, "Cov", withDimnames = FALSE)

+    if (!.isHDF5ArrayBacked(M) && !.isHDF5ArrayBacked(Cov)) {

+        ## TODO: This check is super expensive if M or Cov is a HDF5Matrix, so

+        ##       we skip it for the time being

+        if(min(assay(object, "M", withDimnames = FALSE)) < 0)

+            msg <- validMsg(msg, "the 'M' assay has negative entries")

+        if(min(assay(object, "Cov", withDimnames = FALSE)) < 0)

+            msg <- validMsg(msg, "the 'Cov' assay has negative entries")

+        if(max(assay(object, "M", withDimnames = FALSE) -

+               assay(object, "Cov", withDimnames = FALSE)) > 0.5)

+            msg <- validMsg(msg, "the 'M' assay has at least one entry bigger than the 'Cov' assay")

+    }

+    if(!is.null(rownames(M)) || !is.null(rownames(Cov)) ||

        ("coef" %in% assayNames(object) && !is.null(rownames(assay(object, "coef")))) ||

        ("se.coef" %in% assayNames(object) && !is.null(rownames(assay(object, "se.coef")))))

         warning("unnecessary rownames in object")

@@ -34,6 +42,11 @@ setMethod("show", signature(object = "BSseq"),

               } else {

                   cat("has not been smoothed\n")

               }

+              if (.isHDF5ArrayBacked(object)) {

+                  cat("Some assays are HDF5Array-backed\n")

+              } else {

+                  cat("All assays are in-memory\n")

+              }

           })

 setMethod("pData", "BSseq", function(object) {

@@ -93,7 +106,7 @@ getBSseq <- function(BSseq, type = c("Cov", "M", "gr", "coef", "se.coef", "trans

         return(BSseq@parameters)

     if(type == "gr")

         return(BSseq@rowRanges)

-    

+

 }

 BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

@@ -111,10 +124,14 @@ BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

         stop("'gr' needs to have widths of 1")

     if(is.null(M) || is.null(Cov))

         stop("Need M and Cov")

-    if(!is.matrix(M))

-        stop("'M' needs to be a matrix")

-    if(!is.matrix(Cov))

-        stop("'Cov' needs to be a matrix")

+    M <- .DelayedMatrix(M)

+    Cov <- .DelayedMatrix(Cov)

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

+        coef <- .DelayedMatrix(coef)

+    }

+    if (!is.null(se.coef)) {

+        se.coef <- .DelayedMatrix(se.coef)

+    }

     if(length(gr) != nrow(M) ||

        length(gr) != nrow(Cov) ||

        ncol(Cov) != ncol(M))

@@ -139,9 +156,15 @@ BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

     if(length(unique(sampleNames)) != ncol(M))

         stop("sampleNames need to be unique and of the right length.")

     ## check that 0 <= M <= Cov and remove positions with Cov = 0

-    if(any(M < 0) || any(M > Cov) || any(is.na(M)) || any(is.na(Cov)) ||

-       any(is.infinite(Cov)))

-        stop("'M' and 'Cov' may not contain NA or infinite values and 0 <= M <= Cov")

+    if (!.isHDF5ArrayBacked(M) && !.isHDF5ArrayBacked(Cov)) {

+        ## TODO: This check is super expensive if M or Cov is a HDF5Matrix, so

+        ##       we skip it for the time being

+        if (any(M < 0) || any(M > Cov) || any(is.na(M)) || any(is.na(Cov)) ||

+            any(is.infinite(Cov))) {

+            stop("'M' and 'Cov' may not contain NA or infinite values and ",

+                 "0 <= M <= Cov")

+        }

+    }

     if(rmZeroCov) {

         wh <- which(rowSums(Cov) == 0)

         if(length(wh) > 0) {

@@ -156,7 +179,7 @@ BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

         mm <- as.matrix(findOverlaps(grR, gr))

         mm <- mm[order(mm[,1]),]

         if(length(grR) == length(gr)) {

-            ## only re-ordering is necessary 

+            ## only re-ordering is necessary

             gr <- grR

             M <- M[mm[,2],,drop = FALSE]

             Cov <- Cov[mm[,2],,drop = FALSE]

@@ -171,8 +194,31 @@ BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

             gr <- grR

             sp <- split(mm[,2], mm[,1])[as.character(1:length(grR))]

             names(sp) <- NULL

-            M <- do.call(rbind, lapply(sp, function(ii) colSums(M[ii,, drop = FALSE])))

-            Cov <- do.call(rbind, lapply(sp, function(ii) colSums(Cov[ii,, drop = FALSE])))

+            # TODO: .collapseDelayedMatrix() always return numeric; it may be

+            #       worth coercing M and Cov to integer DelayedMatrix objects,

+            #       which would halve storage requirements and impose some more

+            #       structure on the BSseq class (via new validity method

+            #       checks)

+            # NOTE: Tries to be smart about how collapsed DelayedMatrix should

+            #       be realized

+            if (.isHDF5ArrayBacked(M)) {

+                M_BACKEND <- "HDF5Array"

+            } else {

+                M_BACKEND <- NULL

+            }

+            M <- .collapseDelayedMatrix(x = M,

+                                        sp = sp,

+                                        MARGIN = 2,

+                                        BACKEND = M_BACKEND)

+            if (.isHDF5ArrayBacked(Cov)) {

+                Cov_BACKEND <- "HDF5Array"

+            } else {

+                Cov_BACKEND <- NULL

+            }

+            Cov <- .collapseDelayedMatrix(x = Cov,

+                                          sp = sp,

+                                          MARGIN = 2,

+                                          BACKEND = Cov_BACKEND)

         }

     }

     if(is.null(colnames(M)) || any(sampleNames != colnames(M)))

@@ -180,8 +226,8 @@ BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

     if(is.null(colnames(Cov)) || any(sampleNames != colnames(Cov)))

         colnames(Cov) <- sampleNames

     if(!is.null(coef)) {

-        if(!is.matrix(coef) ||

-           nrow(coef) != nrow(M) ||

+

+        if(nrow(coef) != nrow(M) ||

            ncol(coef) != ncol(M))

             stop("'coef' does not have the right dimensions")

         if(is.null(colnames(coef)) || any(sampleNames != colnames(coef)))

@@ -190,8 +236,7 @@ BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

             rownames(coef) <- NULL

     }

     if(!is.null(se.coef)) {

-        if(!is.matrix(se.coef) ||

-           nrow(se.coef) != nrow(M) ||

+        if(nrow(se.coef) != nrow(M) ||

            ncol(se.coef) != ncol(M))

             stop("'se.coef' does not have the right dimensions")

         if(is.null(colnames(se.coef)) || any(sampleNames != colnames(se.coef)))

@@ -213,19 +258,30 @@ BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

     BSseq

 }

-

 setMethod("updateObject", "BSseq",

           function(object, ...) {

-               if(.hasSlot(object, "assays")) {

-                   # call method for RangedSummarizedExperiment objects

-                   callNextMethod()

-               } else {

-                   BSseq(gr = object@gr, M = object@M, Cov = object@Cov,

-                         coef = object@coef, se.coef = object@se.coef,

-                         trans = object@trans, parameters = object@parameters,

-                         pData = object@phenoData@data)

-               }

-           })

+              # NOTE: identical() is too strong

+              if (isTRUE(all.equal(getBSseq(object, "trans"), .oldTrans))) {

+                  object@trans <- plogis

+              }

+              if (.hasSlot(object, "assays")) {

+                  assays(object) <- endoapply(

+                      assays(object, withDimnames = FALSE), function(assay) {

+                          if (is.null(assay)) {

+                              return(assay)

+                          } else {

+                              .DelayedMatrix(assay)

+                          }

+                      })

+                  # call method for RangedSummarizedExperiment objects

+                  callNextMethod()

+              } else {

+                  BSseq(gr = object@gr, M = object@M, Cov = object@Cov,

+                        coef = object@coef, se.coef = object@se.coef,

+                        trans = object@trans, parameters = object@parameters,

+                        pData = object@phenoData@data)

+              }

+          })

 strandCollapse <- function(BSseq, shift = TRUE) {

@@ -285,7 +341,7 @@ strandCollapse <- function(BSseq, shift = TRUE) {

 ##                ## Based on the delta method

 ##                se.d <- sqrt((data$se.coef[, sample1] * p1 * (1-p1))^2 +

 ##                             (data$se.coef[, sample2] * p2 * (1-p2))^2)

-               

+

 ##                lower <- d - z * se.d

 ##                upper <- d + z * se.d

 ##                out <- data.frame(d = d, lower = lower, upper = upper)

@@ -10,12 +10,33 @@ collapseBSseq <- function(BSseq, columns) {

     else

         columns.idx <- match(names(columns), sampleNames(BSseq))

     sp <- split(columns.idx, columns)

-    M <- do.call(cbind, lapply(sp, function(ss) {

-        rowSums(getBSseq(BSseq, "M")[, ss, drop = FALSE])

-    }))

-    Cov <- do.call(cbind, lapply(sp, function(ss) {

-        rowSums(getBSseq(BSseq, "Cov")[, ss, drop = FALSE])

-    }))

+    # TODO: .collapseDelayedMatrix() always return numeric; it may be

+    #       worth coercing M and Cov to integer DelayedMatrix objects,

+    #       which would halve storage requirements and impose some more

+    #       structure on the BSseq class (via new validity method

+    #       checks)

+    # NOTE: Tries to be smart about how collapsed DelayedMatrix should

+    #       be realized

+    M <- getBSseq(BSseq, "M")

+    if (.isHDF5ArrayBacked(M)) {

+        M_BACKEND <- "HDF5Array"

+    } else {

+        M_BACKEND <- NULL

+    }

+    M <- .collapseDelayedMatrix(x = M,

+                                sp = sp,

+                                MARGIN = 1,

+                                BACKEND = M_BACKEND)

+    Cov <- getBSseq(BSseq, "Cov")

+    if (.isHDF5ArrayBacked(Cov)) {

+        Cov_BACKEND <- "HDF5Array"

+    } else {

+        Cov_BACKEND <- NULL

+    }

+    Cov <- .collapseDelayedMatrix(x = Cov,

+                                  sp = sp,

+                                  MARGIN = 1,

+                                  BACKEND = Cov_BACKEND)

     BSseq(gr = getBSseq(BSseq, "gr"), M = M, Cov = Cov, sampleNames = names(sp))

 }

@@ -108,11 +129,10 @@ getMeth <- function(BSseq, regions = NULL, type = c("smooth", "raw"),

         outMatrix <- matrix(NA, ncol = ncol(BSseq), nrow = length(regions))

         colnames(outMatrix) <- sampleNames(BSseq)

         outMatrix[as.integer(rownames(out)),] <- out

-        return(outMatrix)

+        .DelayedMatrix(outMatrix)

     }

 }

-

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

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

     stopifnot(is(BSseq, "BSseq"))

@@ -175,7 +195,7 @@ getCoverage <- function(BSseq, regions = NULL, type = c("Cov", "M"),

     outMatrix <- matrix(NA, ncol = ncol(BSseq), nrow = length(regions))

     colnames(outMatrix) <- sampleNames(BSseq)

     outMatrix[as.integer(rownames(out)),] <- out

-    return(outMatrix)

+    .DelayedMatrix(outMatrix)

 }

new file mode 100644

@@ -0,0 +1,91 @@

+# Functions/methods that would be good to have in DelayedArray

+

+.rowVars <- function(x, rows = NULL, cols = NULL, ...) {

+    if (is(x, "DelayedArray")) {

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

+            x <- x[rows, ]

+        }

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

+            x <- x[, cols]

+        }

+        row_vars <- rowVars(as.array(x), ...)

+    } else {

+        row_vars <- rowVars(x, rows = rows, cols = cols, ...)

+    }

+    row_vars

+}

+

+.rowSds <- function(x, rows = NULL, cols = NULL, ...) {

+    row_vars <- .rowVars(x, rows = rows, cols = cols, ...)

+    sqrt(row_vars)

+}

+

+.quantile <- function(x, ...) {

+    if (is(x, "DelayedArray")) {

+        x <- as.array(x)

+    }

+    quantile(x, ...)

+}

+

+.DelayedMatrix <- function(x) {

+    x_name <- deparse(substitute(x))

+    X <- try(DelayedArray(x), silent = TRUE)

+    if (is(X, "try-error")) {

+        stop("Could not construct DelayedMatrix from '", x_name, "'",

+             call. = FALSE)

+    }

+    if (!is(X, "DelayedMatrix")) {

+        stop("'", x_name, "' must be matrix-like", call. = FALSE)

+    }

+    X

+}

+

+.isSimpleDelayedMatrix <- function(x) {

+    is(x@seed, "matrix")

+}

+

+# NOTE: Equivalent to rowSums(x[, j, drop = FALSE]) but does it using a

+#       delayed operation and always returns a nrow(x) x 1 DelayedMatrix

+.delayed_rowSums <- function(x, j) {

+    Reduce(`+`, lapply(j, function(jj) x[, jj, drop = FALSE]))

+}

+

+# NOTE: Equivalent to colSums(x[i, , drop = FALSE]) but does it using a

+#       delayed operation and always returns a 1 x ncol(x) DelayedMatrix

+.delayed_colSums <- function(x, i) {

+    Reduce(`+`, lapply(i, function(ii) x[ii, , drop = FALSE]))

+}

+

+# MARGIN = 1: collapse using rowSums

+# MARGIN = 2: collapse using colSums

+.collapseDelayedMatrix <- function(x, sp, MARGIN, BACKEND = NULL) {

+    stopifnot(is(x, "DelayedMatrix"))

+    if (MARGIN == 1) {

+        if (is.null(BACKEND)) {

+            collapsed_x <- do.call(cbind, lapply(sp, function(j) {

+                rowSums(x[, j, drop = FALSE])

+            }))

+        } else {

+            collapsed_x <- do.call(cbind, lapply(sp, function(j) {

+                .delayed_rowSums(x, j)

+            }))

+            # NOTE: Need to manually add colnames when using this method

+            colnames(collapsed_x) <- names(sp)