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

 Package: bsseq

-Version: 1.7.2

+Version: 1.7.6

 Title: Analyze, manage and store bisulfite sequencing data

 Description: A collection of tools for analyzing and visualizing bisulfite

     sequencing data.

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

     data.table,

     S4Vectors,

     R.utils (>= 2.0.0),

-    matrixStats

+    matrixStats,

+    permute

 Suggests:

     RUnit,

     bsseqData,

@@ -15,7 +15,7 @@ importFrom(graphics, "abline", "axis", "layout", "legend", "lines",

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

 import(parallel)

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

-importFrom(matrixStats, "rowSds", "rowVars")

+importFrom(matrixStats, "rowSds", "rowVars", "rowMaxs", "rowMins")

 import(IRanges)

 import(GenomicRanges)

 import(SummarizedExperiment)

@@ -32,6 +32,7 @@ importFrom(gtools, "combinations")

 importFrom(data.table, "fread", "setnames")

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

 import(limma)

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

 ##

 ## Exporting

@@ -63,7 +63,7 @@ smoothSds <- function(BSseqStat, k = 101, qSd = 0.75, mc.cores = 1,

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

                          }, mc.cores = mc.cores))

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

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

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

     else

         BSseqStat@stats <- c(getStats(BSseqStat), list(smoothSds = smoothSds))

     BSseqStat

@@ -73,7 +73,7 @@ smoothSds <- function(BSseqStat, k = 101, qSd = 0.75, mc.cores = 1,

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

     stopifnot(is(BSseqStat, "BSseqStat"))

     if(is.null(coef)) {

-        coef <- 1:ncol(BSseqStat$rawTstats)

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

     }

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

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

@@ -136,3 +136,31 @@ localCorrectStat <- function(BSseqStat, threshold = c(-15,15), mc.cores = 1, ver

     BSseqTstat

 }

+fstat.pipeline <- function(BSseq, design, contrasts, cutoff, fac, nperm = 1000,

+                           coef = NULL, maxGap.sd = 10 ^ 8, maxGap.dmr = 300,

+                           mc.cores = 1) {

+    bstat <- BSmooth.fstat(BSseq = BSseq, design = design,

+                           contrasts = contrasts,

+                           returnModelCoefficients = TRUE)

+    bstat <- smoothSds(bstat)

+    bstat <- computeStat(bstat, coef = coef)

+    dmrs <- dmrFinder(bstat, cutoff = cutoff)

+    idxMatrix <- permuteAll(nperm, design)

+    nullDist <- getNullDistribution_BSmooth.fstat(BSseq = BSseq,

+                                                  idxMatrix = idxMatrix,

+                                                  design = design,

+                                                  contrasts = contrasts,

+                                                  coef = coef,

+                                                  cutoff = cutoff,

+                                                  maxGap.sd = maxGap.sd,

+                                                  maxGap.dmr = maxGap.dmr,

+                                                  mc.cores = mc.cores)

+    fwer <- getFWER.fstat(null = c(list(dmrs), nullDist), type = "dmrs")

+    dmrs$fwer <- fwer

+    meth <- getMeth(BSseq, dmrs, what = "perRegion")

+    meth <- t(apply(meth, 1, function(xx) tapply(xx, fac, mean)))

+    dmrs <- cbind(dmrs, meth)

+    dmrs$maxDiff <- rowMaxs(meth) - rowMins(meth)

+    list(bstat = bstat, dmrs = dmrs, idxMatrix = idxMatrix, nullDist = nullDist)

+}

+

@@ -39,8 +39,10 @@ setMethod("[", "BSseqStat", function(x, i, ...) {

     if(missing(i))

         return(x)

     x@gr <- x@gr[i]

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

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

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

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

         }

@@ -79,7 +79,7 @@ combineList <- function(x, ...) {

         Cov <- do.call(cbind, lapply(x, function(xx) getBSseq(xx, "Cov")))

     } else {

         gr <- sort(reduce(do.call(c, unname(lapply(x, granges))), min.gapwidth = 0L))

-        sampleNames <- do.call(c, lapply(x, sampleNames))

+        sampleNames <- do.call(c, unname(lapply(x, sampleNames)))

         M <- matrix(0, ncol = length(sampleNames), nrow = length(gr))

         colnames(M) <- sampleNames

         Cov <- M

@@ -121,7 +121,7 @@ setMethod("findOverlaps",

           signature(query = "hasGRanges", subject = "GenomicRanges"),

           function (query, subject, maxgap = 0L, minoverlap = 1L,

                     type = c("any", "start", "end", "within", "equal"),

-                    select = c("all", "first"),

+                    select = c("all", "first", "last", "arbitrary"),

                     ignore.strand = FALSE, ...) {

               findOverlaps(query = granges(query), subject = subject,

                            maxgap = maxgap, minoverlap = minoverlap,

@@ -133,7 +133,7 @@ setMethod("findOverlaps",

           signature(query = "hasGRanges", subject = "hasGRanges"),

           function (query, subject, maxgap = 0L, minoverlap = 1L,

                     type = c("any", "start", "end", "within", "equal"),

-                    select = c("all", "first"),

+                    select = c("all", "first", "last", "arbitrary"),

                     ignore.strand = FALSE, ...) {

               findOverlaps(query = granges(query), subject = granges(subject),

                            maxgap = maxgap, minoverlap = minoverlap,

@@ -145,7 +145,7 @@ setMethod("findOverlaps",

           signature(query = "GenomicRanges", subject = "hasGRanges"),

           function (query, subject, maxgap = 0L, minoverlap = 1L,

                     type = c("any", "start", "end", "within", "equal"),

-                    select = c("all", "first"),

+                    select = c("all", "first", "last", "arbitrary"),

                     ignore.strand = FALSE, ...) {

               findOverlaps(query = query, subject = granges(subject),

                            maxgap = maxgap, minoverlap = minoverlap,

@@ -19,12 +19,22 @@ getNullDistribution_BSmooth.tstat <- function(BSseq, idxMatrix1, idxMatrix2,

     nullDist

 }

-getNullDistribution_BSmooth.fstat <- function(BSseq, idxMatrix, design, contrasts,

-                                              coef = NULL, cutoff, maxGap.sd, maxGap.dmr, mc.cores = 1) {

-    message(sprintf("[getNullDistribution_BSmooth.fstat] performing %d permutations\n", nrow(idxMatrix)))

-    nullDist <- mclapply(1:nrow(idxMatrix), function(ii) {

+getNullDistribution_BSmooth.fstat <- function(BSseq,

+                                              idxMatrix,

+                                              design, contrasts,

+                                              coef = NULL, cutoff,

+                                              maxGap.sd, maxGap.dmr,

+                                              mc.cores = 1) {

+

+    message(sprintf("[getNullDistribution_BSmooth.fstat] performing %d permutations\n",

+                    nrow(idxMatrix)))

+    # NOTE: Using mc.preschedule = TRUE

+    nullDist <- mclapply(seq_len(nrow(idxMatrix)), function(ii) {

         ptime1 <- proc.time()

-        bstat <- BSmooth.fstat(BSseq[, idxMatrix[ii,]], design = design, contrasts = contrasts)

+        # NOTE: More efficient to permute design matrix using idxMatrix[ii, ] than

+        #       permute the raw data with BSseq[, idxMatrix[ii, ]]

+        bstat <- BSmooth.fstat(BSseq, design = design[idxMatrix[ii, ], ],

+                               contrasts = contrasts)

         bstat <- smoothSds(bstat, maxGap = maxGap.sd)

         bstat <- computeStat(bstat, coef = coef)

         dmrs0 <- dmrFinder(bstat, cutoff = cutoff, maxGap = maxGap.dmr)

@@ -32,10 +42,18 @@ getNullDistribution_BSmooth.fstat <- function(BSseq, idxMatrix, design, contrast

         stime <- (ptime2 - ptime1)[3]

         message(sprintf("[getNullDistribution_BSmooth.fstat] completing permutation %d in %.1f sec\n", ii, stime))

         dmrs0

-    }, mc.cores = min(nrow(idxMatrix), mc.cores), mc.preschedule = FALSE)

+    }, mc.cores = min(nrow(idxMatrix), mc.cores), mc.preschedule = TRUE)

     nullDist

 }

+permuteAll <- function(nperm, design) {

+    message(sprintf("[permuteAll] performing %d unrestricted permutations of the design matrix\n", nperm))

+

+    CTRL <- how(nperm = nperm)

+    # NOTE: shuffleSet() returns a nperm * nrow(design) matrix of permutations

+    idxMatrix <- shuffleSet(n = nrow(design), control = CTRL)

+}

+

 getNullBlocks_BSmooth.tstat <- function(BSseq, idxMatrix1, idxMatrix2, estimate.var = "same",

                           mc.cores = 1) {

     getNullDistribution_BSmooth.tstat(BSseq = BSseq, idxMatrix1 = idxMatrix1,

@@ -54,9 +72,6 @@ getNullDmrs_BSmooth.tstat <- function(BSseq, idxMatrix1, idxMatrix2, estimate.va

                        mc.cores = mc.cores)

 }

-

-

-

 subsetDmrs <- function(xx) {

     if(is.null(xx) || is(xx, "try-error"))

         return(NULL)

@@ -81,26 +96,66 @@ getFWER <- function(null, type = "blocks") {

     null <- null[-1]

     null <- null[!sapply(null, is.null)]

     better <- sapply(1:nrow(reference), function(ii) {

-        meanDiff <- abs(reference$meanDiff[ii])

+        # meanDiff <- abs(reference$meanDiff[ii])

+        areaStat <- abs(reference$areaStat[ii])

+        width <- reference$width[ii]

+        n <- reference$n[ii]

+        if (type == "blocks") {

+            out <- sapply(null, function(nulldist) {

+                # any(abs(nulldist$meanDiff) >= meanDiff &

+                        # nulldist$width >= width)

+                any(abs(nulldist$areaStat) >= areaStat &

+                        nulldist$width >= width)

+            })

+        }

+        if (type == "dmrs") {

+            out <- sapply(null, function(nulldist) {

+                # any(abs(nulldist$meanDiff) >= meanDiff &

+                #     nulldist$n >= n)

+                any(abs(nulldist$areaStat) >= areaStat &

+                        nulldist$n >= n)

+            })

+        }

+        sum(out)

+    })

+    better

+}

+

+# NOTE: Identical to getFWER() except uses areaStat rather than meanDiff

+#       to compare regions.

+getFWER.fstat <- function(null, type = "blocks") {

+    reference <- null[[1]]

+    null <- null[-1]

+    null <- null[!sapply(null, is.null)]

+    # TODO: Will break if null == list(), which can occur in practice (although

+    #       rarely).

+    better <- sapply(1:nrow(reference), function(ii) {

+        # meanDiff <- abs(reference$meanDiff[ii])

+        areaStat <- abs(reference$areaStat[ii])

         width <- reference$width[ii]

         n <- reference$n[ii]

-        if(type == "blocks") {

+        if (type == "blocks") {

             out <- sapply(null, function(nulldist) {

-                any(abs(nulldist$meanDiff) >= meanDiff &

+                # any(abs(nulldist$meanDiff) >= meanDiff &

+                # nulldist$width >= width)

+                any(abs(nulldist$areaStat) >= areaStat &

                         nulldist$width >= width)

             })

         }

-        if(type == "dmrs") {

+        if (type == "dmrs") {

             out <- sapply(null, function(nulldist) {

-                any(abs(nulldist$meanDiff) >= meanDiff &

-                    nulldist$n >= n)

+                # any(abs(nulldist$meanDiff) >= meanDiff &

+                #     nulldist$n >= n)

+                any(abs(nulldist$areaStat) >= areaStat &

+                        nulldist$n >= n)

             })

         }

         sum(out)

     })

     better

 }

-       

+

+# TODO: Simplify makeIdxMatrix() by using permute package

 makeIdxMatrix <- function(group1, group2, testIsSymmetric = TRUE, includeUnbalanced = TRUE) {

     groupBoth <- c(group1, group2)

     idxMatrix1 <- NULL

@@ -15,7 +15,7 @@ plotAnnoTrack <- function(gr, annoTrack) {

 plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addRegions = NULL,

                             annoTrack = NULL, col = NULL, lty = NULL, lwd = NULL, 

-                            BSseqTstat = NULL, stat = "tstat.corrected", stat.col = "black",

+                            BSseqStat = NULL, stat = "tstat.corrected", stat.col = "black",

                             stat.lwd = 1, stat.lty = 1, stat.ylim = c(-8,8),

                             mainWithWidth = TRUE, regionCol = alpha("red", 0.1), addTicks = TRUE,

                             addPoints = FALSE, pointsMinCov = 5, highlightMain = FALSE, verbose = TRUE) {

@@ -32,8 +32,8 @@ plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addReg

     }

     gr <- resize(gr, width = 2*extend + width(gr), fix = "center")

     BSseq <- subsetByOverlaps(BSseq, gr)

-    if(!is.null(BSseqTstat))

-        BSseqTstat <- subsetByOverlaps(BSseqTstat, gr)

+    if(!is.null(BSseqStat))

+        BSseqStat <- subsetByOverlaps(BSseqStat, gr)

     if(length(start(BSseq)) == 0)

         stop("No overlap between BSseq data and regions")

@@ -43,7 +43,7 @@ plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addReg

     for(ii in seq(along = gr)) {

         if(verbose) cat(sprintf("[plotManyRegions]   plotting region %d (out of %d)\n", ii, nrow(regions)))

         plotRegion(BSseq = BSseq, region = regions[ii,], extend = extend,

-                   col = col, lty = lty, lwd = lwd, main = main[ii], BSseqTstat = BSseqTstat,

+                   col = col, lty = lty, lwd = lwd, main = main[ii], BSseqStat = BSseqStat,

                    stat = stat, stat.col = stat.col, stat.lwd = stat.lwd,

                    stat.lty = stat.lty, stat.ylim = stat.ylim,

                    addRegions = addRegions, regionCol = regionCol, mainWithWidth = mainWithWidth,

@@ -211,7 +211,7 @@ plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addReg

 plotRegion <- function(BSseq, region = NULL, extend = 0, main = "", addRegions = NULL, annoTrack = NULL,

                        col = NULL, lty = NULL, lwd = NULL,

-                       BSseqTstat = NULL, stat = "tstat.corrected", stat.col = "black",

+                       BSseqStat = NULL, stat = "tstat.corrected", stat.col = "black",

                        stat.lwd = 1, stat.lty = 1, stat.ylim = c(-8,8),

                        mainWithWidth = TRUE,

                        regionCol = alpha("red", 0.1), addTicks = TRUE, addPoints = FALSE,

@@ -219,7 +219,7 @@ plotRegion <- function(BSseq, region = NULL, extend = 0, main = "", addRegions =

     opar <- par(mar = c(0,4.1,0,0), oma = c(5,0,4,2), mfrow = c(1,1))

     on.exit(par(opar))

-    if(is.null(BSseqTstat))

+    if(is.null(BSseqStat))

         layout(matrix(1:2, ncol = 1), heights = c(2,1))

     else

         layout(matrix(1:3, ncol = 1), heights = c(2,2,1))

@@ -230,16 +230,25 @@ plotRegion <- function(BSseq, region = NULL, extend = 0, main = "", addRegions =

                     pointsMinCov = pointsMinCov, highlightMain = highlightMain)

     gr <- .bsGetGr(BSseq, region, extend)

-    if(!is.null(BSseqTstat)) {

-        BSseqTstat <- subsetByOverlaps(BSseqTstat, gr)

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

+        BSseqStat <- subsetByOverlaps(BSseqStat, gr)

+        if(is(BSseqStat, "BSseqTstat")) {

+            stat.values <- getStats(BSseqStat, what = stat)

+            stat.type <- "tstat"

+        }

+        if(is(BSseqStat, "BSseqStat")) {

+            stat.type <- getStats(BSseqStat, what = "stat.type")

+            if(stat.type == "tstat")

+                stat.values <- getStats(BSseqStat, what = "stat")

+            if(stat.type == "fstat")

+                stat.values <- sqrt(getStats(BSseqStat, what = "stat"))

+        }

         plot(start(gr), 0.5, type = "n", xaxt = "n", yaxt = "n",

-             ylim = stat.ylim, xlim = c(start(gr), end(gr)), xlab = "", ylab = "t-stat")

+             ylim = stat.ylim, xlim = c(start(gr), end(gr)), xlab = "", ylab = stat.type)

         axis(side = 2, at = c(-5,0,5))

         abline(h = 0, col = "grey60")

-        mapply(function(stat, col, lty, lwd) {

-            .bsPlotLines(start(BSseqTstat), BSseqTstat@stats[, stat],

-                         lty = lty, plotRange = c(start(gr), end(gr)), col = col, lwd = lwd)

-        }, stat = stat, col = stat.col, lty = stat.lty, lwd = stat.lwd)

+        .bsPlotLines(start(BSseqStat), stat.values, lty = stat.lty, col = stat.col, lwd = stat.lwd,

+                     plotRange = c(start(gr), end(gr)))

     }

     if(!is.null(annoTrack))

@@ -13,7 +13,7 @@ read.bismark <- function(files,

     }

     fileType <- match.arg(fileType)

     if (verbose) {

-        message(paste0("Assuming file type is", fileType))

+        message(paste0("Assuming file type is ", fileType))

     }

     ## Process each file

     idxes <- seq_along(files)

@@ -119,7 +119,7 @@ read.bismarkCytosineReportRaw <- function(thisfile,

                   strand = out[[3]])

     ## Create BSseq instance from 'out'

-    BSseq(gr = gr,

+    BSseq(gr = gr, sampleNames = thisSampleName,

           M = as.matrix(out[[4L]]),

           Cov = as.matrix(out[[4L]] + out[[5L]]))

 }

@@ -2,6 +2,14 @@

 \title{bsseq news}

 \encoding{UTF-8}

+\section{Version 1.7.x}{

+  \itemize{

+    \item Fixing an error with reading multiple samples aligned with

+    bismark in the format "cytosineReport".

+  }

+}

+

+

 \section{Version 1.5.x}{

   \itemize{

     \item new function strandCollapse for collapsing forward and reverse