@@ -1,33 +1,12 @@

 language: r

-sudo: required

-

-# Choose your operating system(s) (default: linux)

-os:

-  - linux

-  - osx

-

-# Bioconductor-specific stuff. If targetting bioc-devel then need to use

-# "bioc_use_devel: true", otherwise remove or comment-out this line.

-# TODO: Should --timings be added to r_check_args to mimic BioC build machines?

+r:

+  - devel

+sudo: false

+cache: packages

 bioc_required: true

-bioc_use_devel: true

-# On CRAN all warnings are treated as errors; this is not true on Bioconductor.

 warnings_are_errors: false

-

-# Need this if also using covr to get code coverage

+# Used for code coverage testing with covr and codecov

 r_github_packages:

   - jimhester/covr

 after_success:

   - Rscript -e 'covr::codecov()'

-

-# bsseq-specific stuff

-r_build_args: "--no-build-vignettes --no-manual"

-# TODO: Travis is reporting an error when checking the vignette, hence the

-# inclusion of --no-vignettes in r_check_args

-r_check_args: "--no-vignettes"

-# Must 'manually' install packages in SUGGESTS in order to run full tests

-r_binary_packages:

-  - RUnit

-bioc_packages:

-  - bsseqData

-  - BiocStyle

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

 Package: bsseq

-Version: 1.7.7

+Version: 1.7.13

 Title: Analyze, manage and store bisulfite sequencing data

 Description: A collection of tools for analyzing and visualizing bisulfite

     sequencing data.

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

 Depends:

-    R (>= 2.15),

+    R (>= 3.3),

     methods,

     BiocGenerics,

     GenomicRanges (>= 1.23.7),

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

     data.table,

     S4Vectors,

     R.utils (>= 2.0.0),

-    matrixStats,

+    matrixStats (>= 0.50.0),

     permute

 Suggests:

     RUnit,

@@ -68,7 +68,7 @@ export("BSseq", "getMeth", "getCoverage", "getBSseq", "getStats",

        "read.umtab", "read.umtab2", "read.bsmooth", "read.bismark",

        "poissonGoodnessOfFit", "binomialGoodnessOfFit",

        "data.frame2GRanges", "BSseqTstat",

-       "BSseqStat", "BSmooth.fstat", "smoothSds", "computeStat")

+       "BSseqStat")

 S3method("print", "chisqGoodnessOfFit")

 S3method("plot", "chisqGoodnessOfFit")

@@ -1,10 +1,10 @@

-BSmooth.fstat <- function(BSseq, design, contrasts, returnModelCoefficients = FALSE, verbose = TRUE){

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

     stopifnot(is(BSseq, "BSseq"))

     stopifnot(hasBeenSmoothed(BSseq))

-        

+

     ## if(any(rowSums(getCoverage(BSseq)[, unlist(groups)]) == 0))

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

-    

+

     if(verbose) cat("[BSmooth.fstat] fitting linear models ... ")

     ptime1 <- proc.time()

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

@@ -30,9 +30,6 @@ BSmooth.fstat <- function(BSseq, design, contrasts, returnModelCoefficients = FA

     stats <- list(rawSds = rawSds,

                   cor.coefficients = cor.coefficients,

                   rawTstats = rawTstats)

-    if(returnModelCoefficients) {

-        stats$modelCoefficients <- fit$coefficients

-    }

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

                      stats = stats, parameters = parameters)

     out

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

     BSseqStat

 }

-

+# Quieten R CMD check

+globalVariables("tstat")

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

     stopifnot(is(BSseqStat, "BSseqStat"))

     if(is.null(coef)) {

@@ -140,11 +138,10 @@ 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)

+                           contrasts = contrasts)

     bstat <- smoothSds(bstat)

     bstat <- computeStat(bstat, coef = coef)

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

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

     idxMatrix <- permuteAll(nperm, design)

     nullDist <- getNullDistribution_BSmooth.fstat(BSseq = BSseq,

                                                   idxMatrix = idxMatrix,

@@ -164,3 +161,38 @@ fstat.pipeline <- function(BSseq, design, contrasts, cutoff, fac, nperm = 1000,

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

 }

+fstat.comparisons.pipeline <- function(BSseq, design, contrasts, cutoff, fac,

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

+                                       verbose = TRUE) {

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

+                           design = design,

+                           contrasts = contrasts,

+                           verbose = verbose)

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

+    # NOTE: Want to keep the bstat object corresponding to the original fstat

+    #       and not that from the last t-tsts in the following lapply()

+    bstat_f <- bstat

+    if (verbose) {

+        cat(paste0("[fstat.comparisons.pipeline] making ", ncol(contrasts),

+                   " comparisons ... \n"))

+    }

+    l <- lapply(seq_len(ncol(contrasts)), function(coef) {

+        if (verbose) {

+            cat(paste0("[fstat.comparisons.pipeline] ",

+                       colnames(contrasts)[coef], "\n"))

+        }

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

+        dmrs <- dmrFinder(bstat, cutoff = cutoff, maxGap = maxGap.dmr,

+                          verbose = verbose)

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

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

+        }

+        dmrs

+    })

+    names(l) <- colnames(contrasts)

+    list(bstat = bstat, dmrs = l)

+}

+

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

     ptime1 <- proc.time()

     switch(estimate.var,

            "group2" = {

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

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

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

                    smoothSd(rawSds[idx], k = k)

                }, mc.cores = mc.cores))

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

                tstat.sd <- smoothSds * scale

            },

            "same" = {

-               rawSds <- sqrt( ((length(group1) - 1) * rowVars(allPs[, group1, drop = FALSE]) +

-                                (length(group2) - 1) * rowVars(allPs[, group2, drop = FALSE])) /

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

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

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

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

                    smoothSd(rawSds[idx], k = k)

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

                          rmZeroCov = FALSE,

                          strandCollapse = TRUE,

                          fileType = c("cov", "oldBedGraph", "cytosineReport"),

+                         mc.cores = 1,

                          verbose = TRUE) {

     ## Argument checking

     if (anyDuplicated(files)) {

@@ -15,10 +16,18 @@ read.bismark <- function(files,

     if (verbose) {

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

     }

+    if (strandCollapse && (fileType == "cov" || fileType == "oldBedGraph")) {

+        stop("'strandCollapse' must be 'FALSE' if 'fileType' is '", fileType, "'")

+    }

+    ## SummarizedExperiment validator makes calls to identical() that will fail

+    ## due to how sampleNames are propagated sometimes with and without names().

+    ## To make life simpler, we simply strip the names()

+    sampleNames <- unname(sampleNames)

+

     ## Process each file

     idxes <- seq_along(files)

     names(idxes) <- sampleNames

-    allOut <- lapply(idxes, function(ii) {

+    allOut <- mclapply(idxes, function(ii) {

         if (verbose) {

             cat(sprintf("[read.bismark] Reading file '%s' ... ", files[ii]))

         }

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

                                                  rmZeroCov = rmZeroCov,

                                                  keepContext = FALSE)

         }

-        if (strandCollapse && !all(runValue(strand(out)) == "*")) {

+        if (strandCollapse) {

             out <- strandCollapse(out)

         }

         ptime2 <- proc.time()

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

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

         }

         out

-    })

+    }, mc.cores = mc.cores)

     if (verbose) {

         cat(sprintf("[read.bismark] Joining samples ... "))

@@ -121,5 +130,6 @@ read.bismarkCytosineReportRaw <- function(thisfile,

     ## Create BSseq instance from 'out'

     BSseq(gr = gr, sampleNames = thisSampleName,

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

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

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

+          rmZeroCov = rmZeroCov)

 }

@@ -9,6 +9,15 @@ test_read.bismark_coverage <- function() {

                           fileType = "cov",

                           verbose = FALSE)

     checkTrue(is(bsseq, "BSseq"))

+    # Regression check

+    # Should not fail if sampleNames have names()

+    bsseq <- read.bismark(files = infile,

+                          sampleNames = c(test = "test_data"),

+                          rmZeroCov = FALSE,

+                          strandCollapse = FALSE,

+                          fileType = "cov",

+                          verbose = FALSE)

+    checkTrue(is(bsseq, "BSseq"))

     # Should also work because the "cov" fileType is the same as the

     # "oldBedGraph" fileType

new file mode 100644

@@ -0,0 +1,62 @@

+\name{BSmooth.fstat}

+\alias{BSmooth.fstat}

+\title{

+    Compute F-statistics based on smoothed whole-genome bisulfite

+    sequencing data.

+}

+\description{

+    Compute F-statistics based on smoothed whole-genome bisulfite

+    sequencing data.

+}

+\usage{

+BSmooth.fstat(BSseq, design, contrasts, verbose = TRUE)

+}

+\arguments{

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

+    \item{design}{The design matrix of the bisulfite-sequencing experiment,

+        with rows corresponding to arrays and columns to coefficients to be

+        estimated.}

+    \item{contrasts}{Numeric matrix with rows corresponding to columns in

+        \code{design} and columns containing contrasts. May be a vector if

+        there is only one contrast.}

+    \item{verbose}{Should the function be verbose?}

+}

+\details{

+    \strong{TODO}

+}

+\value{

+    An object of class \linkS4class{BSseqStat}.

+}

+\author{

+    Kasper Daniel Hansen \email{khansen@jhsph.edu}

+}

+\seealso{

+    \code{\link{BSmooth}} for the input object and

+    \linkS4class{BSseq} for its class.

+    \linkS4class{BSseqTstat} describes the return class. This

+    function is likely to be followed by the use of

+    \code{\link{dmrFinder}}.

+}

+\examples{

+  \dontrun{

+    if(require(bsseqData)) {

+        data(keepLoci.ex)

+        data(BS.cancer.ex.fit)

+        BS.cancer.ex.fit <- updateObject(BS.cancer.ex.fit)

+        ## Remember to subset the BSseq object, see vignette for explanation

+        ## TODO: Kind of a forced example

+        design <- model.matrix(~0 + BS.cancer.ex.fit$Type)

+        colnames(design) <- gsub("BS\\\.cancer\\\.ex\\\.fit\\\$Type", "",

+                                 colnames(design))

+        contrasts <- makeContrasts(

+            cancer_vs_normal = cancer - normal,

+            levels = design

+        )

+        BS.stat <- BSmooth.fstat(BS.cancer.ex.fit[keepLoci.ex,],

+                                  design,

+                                  contrasts)

+        BS.stat

+    }

+  }

+}

+\keyword{internal}

new file mode 100644

@@ -0,0 +1,64 @@

+\name{BSseqStat-class}

+\Rdversion{1.1}

+\docType{class}

+\alias{BSseqStat}

+\alias{BSseqStat-class}

+\alias{[,BSseqStat-method}

+\alias{[,BSseqStat,ANY-method}

+\alias{show,BSseqStat-method}

+\title{Class BSseqStat}

+\description{

+    A class for representing statistics for smoothed whole-genome

+    bisulfite sequencing data.

+}

+\usage{

+    BSseqStat(gr = NULL, stats = NULL, parameters = NULL)

+}

+\arguments{

+    \item{gr}{The genomic locations as an object of class \code{GRanges}.}

+    \item{stats}{The statistics, as a matrix.}

+    \item{parameters}{A list of parameters.}

+}

+\section{Objects from the Class}{

+    Objects can be created by calls of the form \code{BSseqStat(...)}.

+    However, usually objects are returned by \code{BSmooth.fstat(...)} and

+    not constructed by the user.

+}

+\section{Slots}{

+    \describe{

+        \item{\code{stats}:}{This is a matrix with columns representing

+            various statistics for methylation loci along the genome.}

+        \item{\code{parameters}:}{Object of class \code{list}.  A list of

+            parameters representing how the statistics were computed.}

+        \item{\code{gr}:}{Object of class \code{GRanges} giving genomic

+            locations.}

+    }

+}

+\section{Extends}{

+    Class \code{\linkS4class{hasGRanges}}, directly.

+}

+\section{Methods}{

+    \describe{

+        \item{[}{The subsetting operator; one may only subset in one

+            dimension, corresponding to methylation loci.}

+        \item{show}{The show method.}

+    }

+}

+\section{Utilities}{

+    This class extends \code{hasGRanges} and therefore inherits a number

+    of useful \code{GRanges} methods that operate on the \code{gr} slot,

+    used for accessing and setting the genomic locations and also do

+    \code{subsetByOverlaps}.

+}

+\author{

+    Kasper Daniel Hansen \email{khansen@jhsph.edu}

+}

+\seealso{

+    \code{\linkS4class{hasGRanges}} for accessing the genomic locations.

+    \code{\link{BSmooth.fstat}} for a function

+    that returns objects of class \code{BSseqStat}, and \code{\link{smoothSds}},

+    \code{\link{computeStat}} and \code{\link{dmrFinder}}

+    for functions that operate based on these statistics.  Also see

+    the more specialised \code{\linkS4class{BSseqTstat}}.

+}

+\keyword{classes}

@@ -31,7 +31,7 @@ not constructed by the user..

     \item{\code{parameters}:}{Object of class \code{list}.  A list of

       parameters representing how the t-statistics were computed.}

     \item{\code{gr}:}{Object of class \code{GRanges} giving genomic

-      locations.} 

+      locations.}

   }

 }

 \section{Extends}{

@@ -48,7 +48,7 @@ Class \code{\linkS4class{hasGRanges}}, directly.

   This class extends \code{hasGRanges} and therefore inherits a number

   of useful \code{GRanges} methods that operate on the \code{gr} slot,

   used for accessing and setting the genomic locations and also do

-  \code{subsetByOverlaps}. 

+  \code{subsetByOverlaps}.

   }

 \author{

   Kasper Daniel Hansen \email{khansen@jhsph.edu}

@@ -56,7 +56,7 @@ Class \code{\linkS4class{hasGRanges}}, directly.

 \seealso{

   The package vignette(s).  \code{\linkS4class{hasGRanges}} for accessing

   the genomic locations.  \code{\link{BSmooth.tstat}} for a function

-  that objects of class \code{BSseqTstat}, and \code{\link{dmrFinder}}

+  that returns objects of class \code{BSseqTstat}, and \code{\link{dmrFinder}}

   for a function that computes DMRs based on the t-statistics.  Also see

   \code{\link[bsseqData]{BS.cancer.ex.tstat}} for an example of the

   class in the \pkg{bsseqData} package.

new file mode 100644

@@ -0,0 +1,66 @@

+\name{computeStat}

+\alias{computeStat}

+\title{

+    Compute a test statistic based on smoothed whole-genome bisulfite

+    sequencing data.

+ }

+\description{

+    Compute a test statistic based on smoothed whole-genome bisulfite

+    sequencing data.

+ }

+\usage{

+    computeStat(BSseqStat, coef = NULL)

+}

+

+\arguments{

+    \item{BSseqStat}{An object of class \code{BSseqStat}, typically an object

+        returned by \code{\link{smoothSds}(...)} and not constructed by

+        the user.}

+    \item{coef}{A vector indicating for which coefficients the statistic is to

+        be computed (\code{coef = NULL} corresponds to testing all

+        coefficients). If the length of the \code{coef} is 1 then the

+        corresponding t-statistic is computed, otherwise the corresponding

+        F-statistic is computed.}

+}

+\details{

+    \strong{TODO}

+}

+\value{

+  An object of class \linkS4class{BSseqStat}. More speciically, the input

+  \linkS4class{BSseqStat} object with the computed statistics added to the

+  \code{stats} slot (accessible with \code{\link{getStats}}).

+}

+\author{

+  Kasper Daniel Hansen \email{khansen@jhsph.edu}

+}

+

+\seealso{

+  \code{\link{smoothSds}} for the function to create the appropriate

+  \code{\linkS4class{BSseqStat}} input object.

+  \code{\linkS4class{BSseqStat}} also describes the return class.  This

+  function is likely to be followed by the use of \code{\link{dmrFinder}}.}

+\examples{

+  \dontrun{

+    if(require(bsseqData)) {

+        data(keepLoci.ex)

+        data(BS.cancer.ex.fit)

+        BS.cancer.ex.fit <- updateObject(BS.cancer.ex.fit)

+        ## Remember to subset the BSseq object, see vignette for explanation

+        ## TODO: Kind of a forced example

+        design <- model.matrix(~0 + BS.cancer.ex.fit$Type)

+        colnames(design) <- gsub("BS\\\.cancer\\\.ex\\\.fit\\\$Type", "",

+                                 colnames(design))

+        contrasts <- makeContrasts(

+            cancer_vs_normal = cancer - normal,

+            levels = design

+        )

+        BS.stat <- BSmooth.fstat(BS.cancer.ex.fit[keepLoci.ex,],

+                                  design,

+                                  contrasts)

+        BS.stat <- smoothSds(BS.stat)

+        BS.stat <- computeStat(BS.stat)

+        BS.stat

+    }

+  }

+}

+\keyword{internal}

@@ -11,6 +11,7 @@

              rmZeroCov = FALSE,

              strandCollapse = TRUE,

              fileType = c("cov", "oldBedGraph", "cytosineReport"),

+             mc.cores = 1,

              verbose = TRUE)

 }

 \arguments{

@@ -22,8 +23,13 @@

     all samples be removed. This will result in a much smaller object if

     the data originates from (targeted) capture bisulfite sequencing.}

   \item{strandCollapse}{Should strand-symmetric methylation loci, e.g., CpGs,

-    be collapsed across strands.}

+    be collapsed across strands. This option is only available if

+    \code{fileType = "cytosineReport"} since the other file types do not

+    contain the necessary strand information.}

   \item{fileType}{The format of the input file; see Note for details.}

+  \item{mc.cores}{The number of cores used.  Note that setting

+    \code{mc.cores} to a value greater than 1 is not  supported on MS

+    Windows, see the help page for \code{mclapply}.}

   \item{verbose}{Make the function verbose.}

 }

 \note{

@@ -56,7 +62,7 @@

   There is no standard file extension for this file. The "\code{C-context}" and

   "\code{trinucleotide context}" columns are not currently used by bsseq.

-  The following is an incomplete list of issues to be aware of when using

+  The following is a list of some issues to be aware of when using

   output from Bismark's methylation extractor:

   \itemize{

@@ -75,7 +81,9 @@

     Furthermore, the default co-ordinate system varies by version of Bismark.

     bsseq makes no assumptions about the basis of the genomic co-ordinates and

     it is left to the user to ensure that the appropriate basis is used in the

-    analysis of their data.

+    analysis of their data. Since Bioconductor packages and

+    \code{\linkS4class{GRanges}} use one-based co-ordinates, it

+    is recommended that your Bismark files are also one-based.

   }

 }

 \value{

@@ -96,7 +104,7 @@

   bismarkBSseq <- read.bismark(files = infile,

                                sampleNames = "test_data",

                                rmZeroCov = FALSE,

-                               strandCollapse = TRUE,

+                               strandCollapse = FALSE,

                                fileType = "cov",

                                verbose = TRUE)

   bismarkBSseq

new file mode 100644

@@ -0,0 +1,73 @@

+\name{smoothSds}

+\alias{smoothSds}

+\title{

+    Smooth the standard deviations using a thresholded running mean based on

+    smoothed whole-genome bisulfite sequencing data.

+}

+\description{

+    Smooth the standard deviations using a thresholded running mean based on

+    smoothed whole-genome bisulfite sequencing data.

+}

+\usage{

+    smoothSds(BSseqStat, k = 101, qSd = 0.75, mc.cores = 1, maxGap = 10^8,

+              verbose = TRUE)

+}

+

+\arguments{

+    \item{BSseqStat}{An object of class \code{BSseqStat}, typically an object

+        returned by \code{\link{BSmooth.fstat}(...)} and not constructed by

+        the user.}

+    \item{k}{A positive scalar, see details.}

+    \item{qSd}{A scalar between 0 and 1, see details.}

+    \item{mc.cores}{The number of cores used.  Note that setting

+        \code{mc.cores} to a value greater than 1 is not supported on MS

+        Windows, see the help page for \code{mclapply}.}

+    \item{maxGap}{A scalar greater than 0, see details.}

+    \item{verbose}{Should the function be verbose?}

+}

+\details{

+    The standard deviation estimates are smoothed using a running mean with a

+    width of \code{k} and thresholded using \code{qSd} which sets the minimum

+    standard deviation to be the \code{qSd}-quantile.

+}

+\value{

+    An object of class \linkS4class{BSseqStat}. More speciically, the input

+    \linkS4class{BSseqStat} object with the computed statistics added to the

+    \code{stats} slot (accessible with \code{\link{getStats}}).

+}

+\author{

+    Kasper Daniel Hansen \email{khansen@jhsph.edu}

+}

+

+\seealso{

+    \code{\link{BSmooth.fstat}} for the function to create the appropriate

+    \code{\linkS4class{BSseqStat}} input object.

+    \code{\linkS4class{BSseqStat}} also describes the return class.  This

+    function is likely to be followed by the use of \code{\link{computeStat}}.}

+\examples{

+  \dontrun{

+    if(require(bsseqData)) {

+        data(keepLoci.ex)

+        data(BS.cancer.ex.fit)

+        BS.cancer.ex.fit <- updateObject(BS.cancer.ex.fit)

+        ## Remember to subset the BSseq object, see vignette for explanation

+        ## TODO: Kind of a forced example

+        design <- model.matrix(~0 + BS.cancer.ex.fit$Type)

+        colnames(design) <- gsub("BS\\\.cancer\\\.ex\\\.fit\\\$Type", "",

+                                 colnames(design))

+        contrasts <- makeContrasts(

+            cancer_vs_normal = cancer - normal,

+            levels = design

+        )

+        BS.stat <- BSmooth.fstat(BS.cancer.ex.fit[keepLoci.ex,],

+                                 design,

+                                 contrasts)

+        BS.stat <- smoothSds(BS.stat)

+        ## Comparing the raw standard deviations to the smoothed standard

+        ## deviations

+        summary(getStats(BS.stat, what = "rawSds"))

+        summary(getStats(BS.stat, what = "smoothSds"))

+    }

+  }

+}

+\keyword{internal}

@@ -158,7 +158,7 @@ We also allow for easy computation of Fisher's exact test for each loci, using t

 Objects of class \Rcode{BSseq} contains a \Rcode{GRanges} object with the genomic locations.  This

 \Rcode{GRanges} object can be obtained by \Rcode{granges}.  A number of standard \Rcode{GRanges}

-methods works directly on the \R.code{BSseq} object, such as \Rcode{start}, \Rcode{end},

+methods works directly on the \Rcode{BSseq} object, such as \Rcode{start}, \Rcode{end},

 \Rcode{seqnames} (chromosome names) etc.

 These objects also contains a \Rcode{phenoData} object for sample pheno data.  Useful methods are

@@ -201,7 +201,7 @@ Even in comparisons where we do not observe these large-scale methylation differ

 Once t-statistics have been computed, we can compute differentially methylated regions (DMRs) by

 thresholding the t-statistics.  Here we use a cutoff of $4.6$, which was chosen by looking at the

-quantiles of the t-statistics (for the entire genome).

+quantiles of the t-statistics (for the entire genome)\footnote{See \url{https://support.bioconductor.org/p/78227/} for further discussion on the choice of cutoff. }.

 <<dmrs>>=

 dmrs0 <- dmrFinder(BS.cancer.ex.tstat, cutoff = c(-4.6, 4.6))

 dmrs <- subset(dmrs0, n >= 3 & abs(meanDiff) >= 0.1)