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

 Package: bsseq

-Version: 1.15.3

+Version: 1.15.4

 Encoding: UTF-8

 Title: Analyze, manage and store bisulfite sequencing data

 Description: A collection of tools for analyzing and visualizing bisulfite

@@ -26,7 +26,7 @@ Imports:

     data.table,

     S4Vectors,

     R.utils (>= 2.0.0),

-    matrixStats (>= 0.50.0),

+    DelayedMatrixStats (>= 1.1.12),

     permute,

     limma,

     DelayedArray (>= 0.5.27),

@@ -16,8 +16,8 @@ importFrom(graphics, "abline", "axis", "layout", "legend", "lines",

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

 import(parallel)

 importFrom(locfit, "locfit", "lp")

-importFrom(matrixStats, "rowSds", "rowVars", "rowMaxs", "rowMins", "colMeans2",

-           "colSums2")

+importFrom(DelayedMatrixStats, "rowSds", "rowVars", "colMeans2", "colSums2",

+           "rowSums2", "rowMeans2")

 import(IRanges)

 import(GenomicRanges)

 import(SummarizedExperiment)

@@ -36,7 +36,8 @@ importFrom(R.utils, "isGzipped", "gunzip")

 import(limma)

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

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

-importFrom(DelayedArray, "DelayedArray", "plogis", "pmin2", "pmax2")

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

+           "rowMins")

 ##

 ## Exporting

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

         ##       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)) {

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

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

         }

     }

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

     ptime1 <- proc.time()

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

                      confint = FALSE)

-    group1.means <- rowMeans(allPs[, group1, drop = FALSE], na.rm = TRUE)

-    group2.means <- rowMeans(allPs[, group2, drop = FALSE], na.rm = TRUE)

+    group1.means <- rowMeans2(allPs[, group1, drop = FALSE], na.rm = TRUE)

+    group2.means <- rowMeans2(allPs[, group2, drop = FALSE], na.rm = TRUE)

     ptime2 <- proc.time()

     stime <- (ptime2 - ptime1)[3]

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

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

         }

     }

     if(rmZeroCov) {

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

+        wh <- which(rowSums2(Cov) == 0)

         if(length(wh) > 0) {

             gr <- gr[-wh]

             M <- M[-wh,,drop = FALSE]

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

             return(getBSseq(BSseq, type = type))

         }

         if (what == "perRegionTotal") {

-            return(colSums(getBSseq(BSseq, type = type)))

+            return(colSums2(getBSseq(BSseq, type = type)))

         }

         if (what == "perRegionAverage") {

-            return(colMeans(getBSseq(BSseq, type = type)))

+            return(colMeans2(getBSseq(BSseq, type = type)))

         }

     }

     if (class(regions) == "data.frame") {

@@ -44,15 +44,15 @@

     is(x@seed, "matrix")

 }

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

+# NOTE: Equivalent to rowSums2(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) {

+.delayed_rowSums2 <- 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

+# NOTE: Equivalent to colSums2(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) {

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

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

 }

@@ -63,11 +63,11 @@

     if (MARGIN == 1) {

         if (is.null(BACKEND)) {

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

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

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

             }))

         } else {

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

-                .delayed_rowSums(x, j)

+                .delayed_rowSums2(x, j)

             }))

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

             colnames(collapsed_x) <- names(sp)

@@ -75,11 +75,11 @@

     } else if (MARGIN == 2) {

         if (is.null(BACKEND)) {

             collapsed_x <- do.call(rbind, lapply(sp, function(i) {

-                colSums(x[i, , drop = FALSE])

+                colSums2(x[i, , drop = FALSE])

             }))

         } else {

             collapsed_x <- do.call(rbind, lapply(sp, function(i) {

-                .delayed_colSums(x, i)

+                .delayed_colSums2(x, i)

             }))

             # NOTE: Need to manually add rownames when using this method

             rownames(collapsed_x) <- names(sp)

@@ -12,7 +12,7 @@ fisherTests <- function(BSseq, group1, group2, lookup = NULL,

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

@@ -26,10 +26,10 @@ fisherTests <- function(BSseq, group1, group2, lookup = NULL,

     }

     if(verbose) cat("[fisherTests] setup ... ")

     ptime1 <- proc.time()

-    M1 <- rowSums(getCoverage(BSseq, type = "M")[, group1, drop = FALSE])

-    M2 <- rowSums(getCoverage(BSseq, type = "M")[, group2, drop = FALSE])

-    Cov1 <- rowSums(getCoverage(BSseq, type = "Cov")[, group1, drop = FALSE])

-    Cov2 <- rowSums(getCoverage(BSseq, type = "Cov")[, group2, drop = FALSE])

+    M1 <- rowSums2(getCoverage(BSseq, type = "M")[, group1, drop = FALSE])

+    M2 <- rowSums2(getCoverage(BSseq, type = "M")[, group2, drop = FALSE])

+    Cov1 <- rowSums2(getCoverage(BSseq, type = "Cov")[, group1, drop = FALSE])

+    Cov2 <- rowSums2(getCoverage(BSseq, type = "Cov")[, group2, drop = FALSE])

     keys <- paste(Cov1, Cov2, M1, M2, sep = "_")

     newkeys <- setdiff(keys, oldkeys)

     expand <- matrix(as.integer(do.call(rbind, strsplit(newkeys, "_"))), ncol = 4)

@@ -1,8 +1,8 @@

 poissonGoodnessOfFit <- function(BSseq, nQuantiles = 10^5) {

     Cov <- getBSseq(BSseq, type = "Cov")

-    lambda <- rowMeans(Cov)

+    lambda <- rowMeans2(Cov)

     out <- list(chisq = NULL, df = ncol(Cov) - 1)

-    out$chisq <- rowSums((Cov - lambda)^2 / lambda)

+    out$chisq <- rowSums2((Cov - lambda)^2 / lambda)

     if(nQuantiles > length(out$chisq))

         nQuantiles <- length(out$chisq)

     probs <- seq(from = 0, to = 1, length.out = nQuantiles)

@@ -22,14 +22,14 @@ binomialGoodnessOfFit <- function(BSseq, method = c("MLE"), nQuantiles = 10^5) {

     M <- getCoverage(BSseq, type = "M", what = "perBase")

     Cov <- getCoverage(BSseq, type = "Cov", what = "perBase")

     if(method == "MLE") # This is the MLE under iid

-        p <- rowSums(M) / rowSums(Cov) 

+        p <- rowSums2(M) / rowSums2(Cov)

     ## p <- rowMeans(M / Cov, na.rm = TRUE)

     out <- list(chisq = NULL, quantiles = NULL, df = ncol(Cov) - 1)

     ## This gets tricky.  For the binomial model, we need Cov > 0, but

     ## this implies that a given location may have less then nCol(BSseq)

     ## number of observations, and hence the degrees of freedom will

-    ## vary.  

-    out$chisq <- rowSums((M - Cov*p)^2 / sqrt(Cov * p * (1-p)))

+    ## vary.

+    out$chisq <- rowSums2((M - Cov*p)^2 / sqrt(Cov * p * (1-p)))

     probs <- seq(from = 0, to = 1, length.out = nQuantiles)

     ## Next line will remove all locations where not all samples have coverage

     out$quantiles <- quantile(out$chisq, prob = probs, na.rm = TRUE)

@@ -11,7 +11,7 @@ read.umtab2 <- function(dirs, sampleNames = NULL, rmZeroCov = FALSE,

                                    readCycle = readCycle, keepFilt = keepFilt,

                                    verbose = verbose)

         if(rmZeroCov && length(chrRead) != 0){

-            wh <- which(rowSums(chrRead$U + chrRead$M) > 0)

+            wh <- which(rowSums2(chrRead$U + chrRead$M) > 0)

             chrRead$M <- chrRead$M[wh,]

             chrRead$U <- chrRead$U[wh,]

             if(readCycle) {

@@ -62,7 +62,7 @@ read.umtab2 <- function(dirs, sampleNames = NULL, rmZeroCov = FALSE,

 }

-read.umtab2.chr <- function(files, sampleNames = NULL, 

+read.umtab2.chr <- function(files, sampleNames = NULL,

                             keepM, keepU, readCycle = FALSE, keepFilt = FALSE,

                             verbose = TRUE) {

     columnHeaders <- strsplit(readLines(files[1], n = 1), "\t")[[1]]

@@ -148,7 +148,7 @@ read.umtab2.chr <- function(files, sampleNames = NULL,

                 csums = csums))

 }

-read.umtab2.chr2 <- function(files, sampleNames = NULL, 

+read.umtab2.chr2 <- function(files, sampleNames = NULL,

                              keepM, keepU, verbose = TRUE) {

     columnHeaders <- strsplit(readLines(files[1], n = 1), "\t")[[1]]

     if("strand" %in% columnHeaders)

@@ -187,7 +187,7 @@ read.umtab2.chr2 <- function(files, sampleNames = NULL,

     }, grLoc)

     nPos <- length(grLoc)

     nSamples <- length(files)

-    

+

     M <- matrix(0L, nrow = nPos, ncol = nSamples)

     colnames(M)  <- sampleNames

     U <- Mcy <- Ucy <- M

@@ -248,7 +248,7 @@ read.umtab <- function(dirs, sampleNames = NULL, rmZeroCov = FALSE,

                                   keepM = keepM, keepU = keepU, sampleNames = sampleNames,

                                   verbose = verbose)

         if(rmZeroCov && length(chrRead) != 0){

-            wh <- which(rowSums(chrRead$U + chrRead$M) > 0)

+            wh <- which(rowSums2(chrRead$U + chrRead$M) > 0)

             chrRead$M <- chrRead$M[wh,]

             chrRead$U <- chrRead$U[wh,]

             chrRead$Mcy <- chrRead$Mcy[wh,]

@@ -369,7 +369,7 @@ read.bsmoothDirRaw <- function(dir, seqnames = NULL, keepCycle = FALSE, keepFilt

     what0[int] <- replicate(length(int), integer(0))

     if(!keepCycle)

         what0[c("Mcy", "Ucy")] <- replicate(2, NULL)

-    if(!keepFilt) 

+    if(!keepFilt)

         what0[grep("^filt", names(what0))] <- replicate(length(grep("^filt", names(what0))), NULL)

     outList <- lapply(allChrFiles, function(thisfile) {

         if(verbose)

@@ -448,7 +448,7 @@ read.bsmooth <- function(dirs, sampleNames = NULL, seqnames = NULL, returnRaw =

         }

         ptime2 <- proc.time()

         stime <- (ptime2 - ptime1)[3]

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

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

         out

     })

     if(!returnRaw) {