@@ -1,14 +1,14 @@

 Package: bsseq

-Version: 1.5.2

-Title: Analyze, manage and store bisulfite sequencing data.

-Description: Tools for analyzing and visualizing bisulfite sequencing data

+Version: 1.5.3

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

-Depends: R (>= 2.15), methods, BiocGenerics, S4Vectors, IRanges (>= 2.1.10),

+Depends: R (>= 2.15), methods, BiocGenerics, IRanges (>= 2.1.10),

         GenomicRanges (>= 1.19.6), SummarizedExperiment (>= 0.1.1), parallel,

         matrixStats, GenomeInfoDb

-Imports: scales, stats, graphics, Biobase, locfit, gtools

-Suggests: RUnit, bsseqData

+Imports: scales, stats, graphics, Biobase, locfit, gtools, data.table, S4Vectors

+Suggests: RUnit, bsseqData, BiocStyle

 Collate: hasGRanges.R BSseq_class.R BSseqTstat_class.R BSseq_utils.R combine.R 

         utils.R read.bsmooth.R read.bismark.R BSmooth.R BSmooth.tstat.R dmrFinder.R

         gof_stats.R plotting.R fisher.R permutations.R

@@ -29,6 +29,7 @@ importMethodsFrom(GenomeInfoDb, "seqlengths", "seqlengths<-", "seqinfo", "seqinf

                   "seqnames", "seqnames<-", "seqlevels", "seqlevels<-")

 import(S4Vectors)

 importFrom(gtools, "combinations")

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

 ##

 ## Exporting

@@ -58,7 +59,7 @@ exportMethods("[", "show",

 export("BSseq", "getMeth", "getCoverage", "getBSseq", "getStats",

        "collapseBSseq", "orderBSseq", "hasBeenSmoothed", "chrSelectBSseq",

        "BSmooth", "BSmooth.tstat", "dmrFinder", "fisherTests",

-       "combineList",

+       "combineList", "strandCollapse",

        "plotRegion", "plotManyRegions", 

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

        "poissonGoodnessOfFit", "binomialGoodnessOfFit",

@@ -228,8 +228,24 @@ setMethod("updateObject", "BSseq",

            })

-

-                  

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

+    if(all(runValue(strand(BSseq)) == "*")) {

+        warning("All loci are unstranded; nothing to collapse")

+        return(BSseq)

+    }

+    if(!(all(runValue(strand(BSseq)) %in% c("+", "-"))))

+        stop("'BSseq' object has a mix of stranded and unstranded loci.")

+    BS.forward <- BSseq[strand(BSseq) == "+"]

+    strand(BS.forward) <- "*"

+    BS.reverse <- BSseq[strand(BSseq) == "-"]

+    strand(BS.reverse) <- "*"

+    if(shift)

+        rowRanges(BS.reverse) <- shift(granges(BS.reverse), -1L)

+    sampleNames(BS.reverse) <- paste0(sampleNames(BS.reverse), "_REVERSE")

+    BS.comb <- combine(BS.forward, BS.reverse)

+    newBSseq <- collapseBSseq(BS.comb, columns = rep(sampleNames(BSseq), 2))

+    newBSseq

+}

 ## getD <- function(data, sample1, sample2, type = c("raw", "fit"),

 ##                  addPositions = FALSE, alpha = 0.95) {

@@ -41,7 +41,8 @@ getNullDmrs <- function(BSseq, idxMatrix1, idxMatrix2, estimate.var = "same",

 subsetDmrs <- function(xx) {

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

         return(NULL)

-    out <- subset(xx, n >= 3 & abs(meanDiff) > 0.1 & invdensity <= 300)

+    out <- xx[ xx[,"n"] >= 3 & abs(xx[, "meanDiff"]) > 0.1 &

+                  xx[, "invdensity"] <= 300, ]

     if(nrow(out) == 0)

         return(NULL)

     out

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

-read.bismark <- function(files, sampleNames, rmZeroCov = FALSE, verbose = TRUE){

+read.bismark <- function(files, sampleNames, rmZeroCov = FALSE, strandCollapse = TRUE,

+                         fileType = c("cytosineReport", "oldBedGraph"), verbose = TRUE){

     ## Argument checking

     if (anyDuplicated(files)){

         stop("duplicate entries in 'files'")

@@ -6,6 +7,7 @@ read.bismark <- function(files, sampleNames, rmZeroCov = FALSE, verbose = TRUE){

     if (length(sampleNames) != length(files) | anyDuplicated(sampleNames)){

         stop("argument 'sampleNames' has to have the same length as argument 'files', without duplicate entries")

     }

+    fileType <- match.arg(fileType)

     ## Process each file

     idxes <- seq_along(files)

     names(idxes) <- sampleNames

@@ -14,12 +16,19 @@ read.bismark <- function(files, sampleNames, rmZeroCov = FALSE, verbose = TRUE){

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

         }

         ptime1 <- proc.time()

-        raw <- read.bismarkFileRaw(thisfile = files[ii])

+        if(fileType == "oldBedGraph") {

+            raw <- read.bismarkFileRaw(thisfile = files[ii])

+        }

+        if(fileType == "cytosineReport") {

+            raw <- read.bismarkCytosineRaw(thisfile = files[ii], keepContext = FALSE)

+        }

         M <- matrix(mcols(raw)[, "mCount"], ncol = 1)

         Cov <- M + mcols(raw)[, "uCount"]

         mcols(raw) <- NULL

         out <- BSseq(gr = raw, M = M, Cov = Cov,

                      sampleNames = sampleNames[ii], rmZeroCov = rmZeroCov)

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

+            out <- strandCollapse(out)

         ptime2 <- proc.time()

         stime <- (ptime2 - ptime1)[3]

         if (verbose) {

@@ -64,3 +73,26 @@ read.bismarkFileRaw <- function(thisfile, verbose = TRUE){

     mcols(gr) <- df

     gr

 }

+

+read.bismarkCytosineRaw <- function(thisfile, keepContext = FALSE) {

+    out <- fread(thisfile)

+    if(length(out) != 6 && length(out) != 7)

+        stop("unknown file format")

+    if(length(out) == 6) {

+        setnames(out, c("chr", "start", "end", "methPerc", "mCount", "uCount"))

+        gr <- GRanges(seqnames = out[["chr"]],

+                      ranges = IRanges(start = out[["start"]], width = 1),

+                      mCount = out[["mCount"]], uCount = out[["uCount"]])

+    }

+    if(length(out) == 7) {

+        setnames(out, c("chr", "start", "strand", "mCount", "uCount", "type", "context"))

+        gr <- GRanges(seqnames = out[["chr"]], strand = out[["strand"]],

+                      ranges = IRanges(start = out[["start"]], width = 1),

+                      mCount = out[["mCount"]], uCount = out[["uCount"]])

+        if(keepContext) {

+            gr$type <- out[["type"]]

+            gr$context <- out[["context"]]

+        }                      

+    }

+    gr

+}

Binary files a/data/BS.chr22.rda and b/data/BS.chr22.rda differ

@@ -2,9 +2,18 @@

 \title{bsseq news}

 \encoding{UTF-8}

+\section{Version 1.5.x}{

+  \itemize{

+    \item new function strandCollapse for collapsing forward and reverse

+    strand data to be unstranded.

+    \item Updated read.bismark to support the cytosine report files;

+    both formats are support.

+  }

+}

+

 \section{Version 0.11.x}{

   \itemize{

-    \item Converted to using Authors@R in the DESCRIPTIOn file.

+    \item Converted to using Authors@R in the DESCRIPTION file.

     \item plotRegion did not respect the col/lty/lwd argument if given

     explicitely as opposed to through pData().  Reported by Karen

     Conneely <kconnee@emory.edu>.

@@ -17,6 +17,7 @@

 \alias{show,BSseq-method}

 \alias{getBSseq}

 \alias{collapseBSseq}

+\alias{strandCollapse}

 \alias{orderBSseq}

 \alias{chrSelectBSseq}

 \alias{hasBeenSmoothed}

@@ -128,6 +129,18 @@ slots for representing smoothed data. This class is an extension of

       methylation loci.  The input is a list, with each component an object

       of class \code{BSseq}.  The (slower) alternative is to use

       \code{Reduce(combine, list)}. }

+  

+    \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.}

+

   }

 }

 \section{Coercion}{

@@ -2,35 +2,18 @@

 %\VignetteDepends{bsseq}

 %\VignettePackage{bsseq}

 \documentclass[12pt]{article}

-<<options,echo=FALSE,results=hide>>=

-options(width=70)

-@ 

-\SweaveOpts{eps=FALSE,echo=TRUE}

-\usepackage{times}

-\usepackage{color,hyperref}

-\usepackage{fullpage}

-\usepackage[numbers]{natbib}

-\definecolor{darkblue}{rgb}{0.0,0.0,0.75}

-\hypersetup{colorlinks,breaklinks,

-            linkcolor=darkblue,urlcolor=darkblue,

-            anchorcolor=darkblue,citecolor=darkblue}

-

-\newcommand{\Rcode}[1]{{\texttt{#1}}}

-\newcommand{\Rpackage}[1]{{\texttt{#1}}}

-\newcommand{\software}[1]{\textsf{#1}}

-\newcommand{\R}{\software{R}}

+<<style-Sweave, eval=TRUE, echo=FALSE, results=tex>>=

+BiocStyle::latex()

+@

 \newcommand{\bold}[1]{\textbf{#1}} % Fix for R bibentry

-

 \title{The bsseq user's guide}

-\author{Kasper Daniel Hansen \\ \texttt{khansen@jhsph.edu}}

-\date{Modified: October 13, 2012.  Compiled: \today}

+\author{Kasper Daniel Hansen \\ \texttt{kasperdanielhansen@gmail.com}}

+\date{Modified: June 10, 2015.  Compiled: \today}

 \begin{document}

-\setlength{\parskip}{1\baselineskip}

-\setlength{\parindent}{0pt}

 \maketitle

-\section{Capabilities}

+\section*{Capabilities}

 The packages required for this vignette are

@@ -112,7 +95,7 @@ If you use this package, please cite our BSmooth paper:

 print(citation("bsseq"), style = "latex")

 @ 

-\section{Overview}

+\section*{Overview}

 The package assumes that the following data has been extracted from alignments:

 \begin{enumerate}

@@ -170,7 +153,7 @@ be of marginal interest to most users.  See the man page for \Rcode{binomialGood

 We also allow for easy computation of Fisher's exact test for each loci, using the function

 \Rcode{fisherTests}. 

-\section{Using objects of class BSseq}

+\section*{Using objects of class BSseq}

 \subsubsection*{Basic operations}

@@ -301,7 +284,7 @@ getMeth(BS.chr22, regions, type = "raw")

 getMeth(BS.chr22, regions[2], type = "raw", what = "perBase")

 @ 

-\section{Reading data}

+\section*{Reading data}

 \subsubsection*{Alignment output from the BSmooth alignment suite}

@@ -380,7 +363,6 @@ Fisher's exact tests may be efficiently computed using the function \Rcode{fishe

 Binomial and poisson goodness of fit tests statistics may be computed using

 \Rcode{binomialGoodnessOfFit} and \Rcode{poissonGoodnessOfFit}.

-\bibliographystyle{unsrturl}

 \bibliography{bsseq}

@@ -3,30 +3,13 @@

 %\VignetteDepends{bsseqData}

 %\VignettePackage{bsseq}

 \documentclass[12pt]{article}

-<<echo=FALSE,results=hide>>=

-options(width=70)

-@ 

-\SweaveOpts{eps=FALSE,echo=TRUE,eval=TRUE}

-\usepackage{times}

-\usepackage{color,hyperref}

-\usepackage{fullpage}

-\usepackage[numbers]{natbib}

-\usepackage{parskip}

-\definecolor{darkblue}{rgb}{0.0,0.0,0.75}

-\hypersetup{colorlinks,breaklinks,

-            linkcolor=darkblue,urlcolor=darkblue,

-            anchorcolor=darkblue,citecolor=darkblue}

-

-\newcommand{\Rcode}[1]{{\texttt{#1}}}

-\newcommand{\Rclass}[1]{{\texttt{"#1"}}}

-\newcommand{\Rpackage}[1]{{\texttt{#1}}}

-\newcommand{\software}[1]{\textsf{#1}}

-\newcommand{\R}{\software{R}}

+<<style-Sweave, eval=TRUE, echo=FALSE, results=tex>>=

+BiocStyle::latex()

+@

 \newcommand{\bold}[1]{\textbf{#1}} % Fix for R bibentry

-

 \title{Analyzing WGBS with the bsseq package}

-\author{Kasper Daniel Hansen\\ \texttt{khansen@jhsph.edu}}

-\date{Modified: October 13, 2012.  Compiled: \today}

+\author{Kasper Daniel Hansen\\ \texttt{kasperdanielhansen@gmail.com}}

+\date{Modified: June 10, 2015.  Compiled: \today}

 \begin{document}

 \maketitle

@@ -344,7 +327,6 @@ space and would add complications to the internal data structures.

 %% Backmatter

 %%

-\bibliographystyle{unsrturl}

 \bibliography{bsseq}

 \section*{SessionInfo}