@@ -11,13 +11,10 @@ BSmooth(BSseq,

         h = 1000,

         maxGap = 10^8,

         keep.se = FALSE,

-        verbose = TRUE,

         BPPARAM = bpparam(),

-        BACKEND = getRealizationBackend(),

-        ...,

-        parallelBy = c("sample", "chromosome"),

-        mc.preschedule = FALSE,

-        mc.cores = 1)

+        chunkdim = NULL,

+        level = NULL,

+        verbose = getOption("verbose"))

 }

 \arguments{

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

@@ -25,17 +22,9 @@ BSmooth(BSseq,

   \item{h}{The minimum smoothing window, in bases.}

   \item{maxGap}{The maximum gap between two methylation loci, before the smoothing is broken across the gap.

   The default smoothes each chromosome separately.}

-  \item{parallelBy}{\strong{Deprecated}.

-    See section 'Parallelization and progress monitoring' for further details.}

-  \item{mc.preschedule}{\strong{Deprecated}.

-    See section 'Parallelization and progress monitoring' for further details.}

-  \item{mc.cores}{\strong{Deprecated}.

-    See section 'Parallelization and progress monitoring' for further details.}

   \item{keep.se}{Should the estimated standard errors from the smoothing

     algorithm be kept. This will make the return object roughly 30

     percent bigger and is currently not be used for anything in \pkg{bsseq}.}

-  \item{verbose}{\strong{Deprecated}.

-    See section, 'Parallelization and progress monitoring' for further details.}

   \item{BPPARAM}{An optional \linkS4class{BiocParallelParam} instance

     determining the parallel back-end to be used during evaluation. Currently

     supported are \linkS4class{SerialParam} (Unix, Mac, Windows),

@@ -44,13 +33,18 @@ BSmooth(BSseq,

     \linkS4class{BatchJobsParam} (Unix, Mac, Windows, only with the in-memory

     realization backend). See sections 'Parallelization and progress

     monitoring' and 'Realization backends' for further details.}

-  \item{BACKEND}{\code{NULL} or a single string specifying the name of the

-    realization backend. When the backend is set to \code{NULL}, the

-    \code{coef} and \code{se.coef} assays are realized in memory as ordinary

-    matrices, otherwise these are realized with the given \code{BACKEND}.}

-  \item{...}{Additional arguments passed to the \linkS4class{RealizationSink}

-    constructor.

-    See section 'Realization backends' for further details.}

+  \item{chunkdim}{\strong{Only applicable if \code{BACKEND == "HDF5Array"}.}

+    The dimensions of the chunks to use for writing the data to

+    disk. By default, \code{\link[HDF5Array]{getHDF5DumpChunkDim}()} using the

+    dimensions of the returned \linkS4class{BSseq} object will be used. See

+    \code{?\link[HDF5Array]{getHDF5DumpChunkDim}} for more information.}

+  \item{level}{\strong{Only applicable if \code{BACKEND == "HDF5Array"}.}

+    The compression level to use for writing the data to disk. By

+    default, \code{\link[HDF5Array]{getHDF5DumpCompressionLevel}()} will be

+    used. See \code{?\link[HDF5Array]{getHDF5DumpCompressionLevel}} for more

+    information.}

+  \item{verbose}{A \code{logical(1)} indicating whether progress messages

+    should be printed (default \code{TRUE}).}

 }

 \details{

   \code{ns} and \code{h} are passed to the \code{locfit} function. The

@@ -141,14 +141,15 @@ slots for representing smoothed data. This class is an extension of

     \item{\code{strandCollapse(BSseq, shift = TRUE)}}{

       This function operates on a \code{BSseq} objects which has

-      stranded loci (ie. loci where the strand is one of \sQuote{+} pr

-      \sQuote{-}).  It will collapse the methylation and coverage

-      information across the two strands, into one position.  The

-      argument \code{shift} indicates whether the positions for the loci

-      on the reverse strand should be shifted one (ie. the positions for

-      these loci are the positions of the \sQuote{G} in the

-      \sQuote{CpG}; this is the case for Bismark output for example.}

+      stranded loci (i.e. loci where the strand is one of \sQuote{+} or

+      \sQuote{-}). It will collapse the methylation and coverage

+      information across the two strands, unstranding the loci in the process

+      and potentially re-ordering them.

+      The argument \code{shift} indicates whether the positions for the loci

+      on the reverse strand should be shifted one (i.e. the positions for

+      these loci are the positions of the \sQuote{G} in the

+      \sQuote{CpG}; this is the case for Bismark output for example).}

   }

 }

 \section{Coercion}{

@@ -183,18 +184,17 @@ M <- matrix(1:9, 3,3)

 colnames(M) <- c("A1", "A2", "A3")

 BStest <- BSseq(pos = 1:3, chr = c("chr1", "chr2", "chr1"), M = M, Cov = M + 2)

 chrSelectBSseq(BStest, seqnames = "chr1", order = TRUE)

-collapseBSseq(BStest, columns = c("A1" = "A", "A2" = "A", "A3" = "B"))

+collapseBSseq(BStest, group = c("A", "A", "B"))

 #-------------------------------------------------------------------------------

 # An example using a HDF5-backed BSseq object

 #

-library(HDF5Array)

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

-hdf5_BStest <- BSseq(pos = 1:3,

-                     chr = c("chr1", "chr2", "chr1"),

-                     M = hdf5_M,

-                     Cov = hdf5_M + 2)

+hdf5_BStest <- realize(BStest, "HDF5Array")

 chrSelectBSseq(hdf5_BStest, seqnames = "chr1", order = TRUE)

-collapseBSseq(hdf5_BStest, columns = c("A1" = "A", "A2" = "A", "A3" = "B"))

+collapseBSseq(

+    BSseq = hdf5_BStest,

+    group = c("A", "A", "B"),

+    BACKEND = "HDF5Array",

+    type = "integer")

 }

 \keyword{classes}

@@ -13,12 +13,12 @@

              rmZeroCov = FALSE,

              strandCollapse = TRUE,

              BPPARAM = bpparam(),

-             BACKEND = getRealizationBackend()

-             dir = tempdir(),

+             BACKEND = getRealizationBackend(),

+             dir = tempfile("BSseq"),

              replace = FALSE,

              chunkdim = NULL,

              level = NULL,

-             verbose = TRUE)

+             verbose = getOption("verbose"))

 }

 \arguments{

   \item{files}{The path to the files created by running Bismark's methylation

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

     (Unix, Mac, and Windows, limited to single-machine clusters), and

     \linkS4class{BatchJobsParam} (Unix, Mac, Windows, only with the in-memory

     realization backend).

-    See sections 'Parallelization and progress  monitoring' and 'Realization backends' for further details.}

+    See sections 'Parallelization and progress  monitoring' and 'Realization

+    backends' for further details.}

   \item{BACKEND}{\code{NULL} or a single string specifying the name of the

     realization backend.

     When the backend is set to \code{NULL}, the \code{M} and \code{Cov} assays

@@ -80,13 +81,36 @@

     should be printed (default \code{TRUE}).}

 }

+\section{File formats}{

+  The format of each file is automatically detected using the internal function \code{bsseq:::.guessBismarkFileType()}.

+  Files ending in \code{.gz}, \code{.bz2}, \code{.xz}, or \code{.zip} will be automatically decompressed to \code{\link[base]{tempdir}()}.

+  \subsection{Supported file formats}{

+    Bismark's 'genome wide cytosine report' (\url{https://github.com/FelixKrueger/Bismark/tree/master/Docs#the-genome-wide-cytosine-report-optional-is-tab-delimited-in-the-following-format-1-based-coords}) and 'coverage' (\url{https://github.com/FelixKrueger/Bismark/tree/master/Docs#the-coverage-output-looks-like-this-tab-delimited-1-based-genomic-coords}) formats are both supported.

+    If setting \code{loci = NULL}, then we strongly recommend using the 'genome wide cytosine report' output format because this includes strand information for each locus.

+     The 'coverage' output does not contain strand information and so the \code{\link[BiocGenerics]{strand}} of the returned \linkS4class{BSseq} object will be set to \code{*} unless stranded \code{loci} are supplied.

+  }

+

+  \subsection{Unsupported file formats}{

+     Neither the 'bedGraph' output format (\url{https://github.com/FelixKrueger/Bismark/tree/master/Docs#the-bedgraph-output-optional-looks-like-this-tab-delimited-0-based-start-1-based-end-coords}) nor the 'bismark_methylation_extractor' output format (\url{https://github.com/FelixKrueger/Bismark/tree/master/Docs#the-bismark_methylation_extractor-output-is-in-the-form-tab-delimited-1-based-coords}) are supported.

+     The former does not include the required counts of methylated and unmethylated reads hile the is an intermediate file containing read-level, rather than locus-level, data on methylation.

+  }

+  \subsection{One-based vs. zero-based genomic co-ordinates}{

+    The genomic co-ordinates of the Bismark output files may be zero-based or one-based depending on whether the \code{--zero_based} argument was used when running Bismark's methylation extractor.

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

+    \pkg{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.

+

+    Since Bioconductor packages typically use one-based co-ordinates, we strongly recommend that your Bismark files are also one-based.

+  }

+}

+

 \section{Efficient use of \code{read.bismark()}}{

   We recommend the following to achieve fast and efficient importing of Bismark files:

   \itemize{

     \item Specify the set of methylation loci via the \code{loci} argument.

+    \item Use Bismark files in the 'coverage' output format.

     \item Leave \code{rmZeroCov = FALSE}.

-    \item Use a \code{BPPARAM} with a large number of workers (cores).

+    \item Use a \code{BPPARAM} with a moderate number of workers (cores).

     \item Use \code{BACKEND = "HDF5Array"}.

   }

@@ -98,24 +122,28 @@

     \strong{If you are unsure whether the below-described shortcuts apply to your data, leave \code{loci = NULL} and let \code{read.bismark()} identify the set of candidate loci from \code{files}.}

-    If all \code{files} are 'genome wide cytosine reports' for samples aligned to the same reference genome, then all \code{files} contain the exact same set of methylation loci.

+    You may wish to use the \code{\link{findLoci}()} function to find all methylation loci of interest in your reference genome (e.g., all CpGs) and then pass the result via the \code{loci} argument.

+

+    Alternatively, if all \code{files} are 'genome wide cytosine reports' for samples aligned to the same reference genome, then all \code{files} contain the exact same set of methylation loci.

     In this case, you may wish to first construct \code{loci} using the internal function \code{bsseq:::.readBismarkAsFWGRanges()} applied to a single file, e.g., \code{loci = bsseq:::.readBismarkAsFWGRanges(files[1], rmZeroCov, strandCollapse)}.

+  }

-    % TODO: Uncomment once findCytosines is exported and documented.

-    % Alternatively, you may wish to use the \code{\link{findCytosines}()} function to construct all methylation loci of interest (e.g., all CpGs) in your reference genome and then pass this via the \code{loci} argument.

+  \subsection{Using the 'coverage' Bismark files}{

+    It will generally be faster to parse Bismark files in the 'coverage' output format than those in the 'genome wide cytosine report' format This is because the former only includes loci with non-zero coverage and so the file size is often considerably smaller, particularly for shallowly sequenced samples (e.g., those from single-cell bisulfite sequencing).

   }

   \subsection{Leaving \code{rmZeroCov = FALSE}}{

     If you set \code{rmZeroCov = TRUE}, then \code{read.bismark()} must first parse all the \code{files} to identify which loci have zero coverage in all samples and then filter these out from the set of candidate loci.

-    \strong{This will happen regardless of whether \code{loci} is \code{NULL} or is supplied as a \linkS4class{GenomicRanges} of candidate loci.}

+    \strong{This will happen even if you supply \code{loci} with a \linkS4class{GenomicRanges} of candidate loci.}

     Furthermore, any coverage-based filtering of methylation loci is best left until you have constructed your final \linkS4class{BSseq} object.

     In our experience, the final \linkS4class{BSseq} object is often the product of combining multiple \linkS4class{BSseq} objects, each constructed with a separate call to \code{read.bismark()}.

-    In such cases, it is premature to use \code{rmZeroCov = TRUE} when running each \code{read.bismark()} and combining these objects will regretably often lead to an inefficiently constructed final \linkS4class{BSseq} object.

+    In such cases, it is premature to use \code{rmZeroCov = TRUE} when running each \code{read.bismark()}; regretably, combining these objects will often then lead to an inefficiently stored \linkS4class{BSseq} object.

   }

-  \subsection{Using a \code{BPPARAM} with a large number of workers (cores)}{

-    Each file can be processed on its own, so you can process in parallel as many files as you have workers. We have good experience using as many as 40 cores via \code{BPPARAM = MulticoreParam(workers = 40)}.}

+  \subsection{Using a \code{BPPARAM} with a moderate number of workers (cores)}{

+    Each file can be processed on its own, so you can process in parallel as many files as you have workers. However, if using the HDF5Array backend, then writing to the HDF5 file cannot be performed in parallel and so becomes the bottleneck. Nonetheless, by using a moderate number of workers (2 - 10), we can ensure there is processed data available to write to disk as soon as the current write is completed.

+  }

   \subsection{Using \code{BACKEND = "HDF5Array"}}{

     By using the HDF5Array realization backend from \pkg{HDF5Array}, we reduce the amount of data that is kept in-memory at any one time.

@@ -126,28 +154,6 @@

     }

 }

-\section{File formats}{

-  The format of each file is automatically detected using the internal function \code{bsseq:::.guessBismarkFileType()}.

-  Files ending in \code{.gz}, \code{.bz2}, \code{.xz}, or \code{.zip} will be automatically uncompressed.

-  \subsection{Supported file formats}{

-     We strongly recommend that you use the 'genome wide cytosine report' output format

-     (\url{https://github.com/FelixKrueger/Bismark/tree/master/Docs#the-genome-wide-cytosine-report-optional-is-tab-delimited-in-the-following-format-1-based-coords}) as this includes trand information for each locus.

-     The 'coverage' output format (\url{https://github.com/FelixKrueger/Bismark/tree/master/Docs#the-coverage-output-looks-like-this-tab-delimited-1-based-genomic-coords}) is also supported, although all loci from these files will have their \code{\link[BiocGenerics]{strand}} set to \code{*} in the returned \linkS4class{BSseq} object because this file format does not include strand information.

-  }

-

-  \subsection{Unsupported file formats}{

-     Neither the 'bedGraph' output format (\url{https://github.com/FelixKrueger/Bismark/tree/master/Docs#the-bedgraph-output-optional-looks-like-this-tab-delimited-0-based-start-1-based-end-coords}) nor the 'bismark_methylation_extractor' output format (\url{https://github.com/FelixKrueger/Bismark/tree/master/Docs#the-bismark_methylation_extractor-output-is-in-the-form-tab-delimited-1-based-coords}) are supported.

-     The former does not include the required counts of methylated and unmethylated reads hile the is an intermediate file containing read-level, rather than locus-level, data on methylation.

-  }

-  \subsection{One-based vs. zero-based genomic co-ordinates}{

-    The genomic co-ordinates of the Bismark output files may be zero-based or one-based depending on whether the \code{--zero_based} argument was used when running Bismark's methylation extractor.

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

-    \pkg{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.

-

-    Since Bioconductor packages typically use one-based co-ordinates, we strongly recommend that your Bismark files are also one-based.

-  }

-}

-

 \section{Realization backends}{

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

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

@@ -211,10 +217,12 @@

 }

 \seealso{

-  \code{\link{read.bsmooth}()} for parsing output from the BSmooth

+  \itemize{

+    \item \code{\link{read.bsmooth}()} for parsing output from the BSmooth

   alignment suite.

-  \code{\link{read.umtab}()} for parsing legacy (old) formats from the BSmooth alignment suite.

-  \code{\link{collapseBSseq}()} for collapsing (aggregating) data from sample's with multiple Bismark methylation extractor files (e.g., technical replicates).

+    \item \code{\link{read.umtab}()} for parsing legacy (old) formats from the BSmooth alignment suite.

+    \item \code{\link{collapseBSseq}()} for collapsing (aggregating) data from sample's with multiple Bismark methylation extractor files (e.g., technical replicates).

+  }

 }

 \examples{

@@ -252,3 +260,5 @@

 \author{

   Peter Hickey \email{peter.hickey@gmail.com}

 }

+

+% TODO: Mention and link to findLoci().