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

 Package: bsseq

-Version: 0.7.4

+Version: 0.7.5

 Title: Analyze, manage and store bisulfite sequencing data

 Description: Tools for analyzing and visualizing bisulfite sequencing data

 Author: Kasper Daniel Hansen

@@ -251,15 +251,32 @@ combineList <- function(x, ...) {

     stopifnot(all(sapply(x, class) == "BSseq"))

     gr <- getBSseq(x[[1]], "gr")

     trans <- getBSseq(x[[1]], "trans")

-    ok <- sapply(x[-1], function(xx) {

-        identical(gr, getBSseq(xx, "gr")) &&

-            identical(trans, getBSseq(xx, "trans"))

+    sameTrans <- sapply(x[-1], function(xx) {

+        identical(trans, getBSseq(xx, "trans"))

     })

-    if(!all(ok))

-        stop("all elements of '...' in combineList needs to have the same granges and trans")

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

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

-    if(any(!sapply(x, hasBeenSmoothed))) {

+    if(!all(sameTrans))

+        stop("all elements of '...' in combineList needs to have the same trans")

+    sameGr <- sapply(x[-1], function(xx) {

+        identical(trans, getBSseq(xx, "gr"))

+    })

+    if(all(sameGr)) {

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

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

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

+        colnames(M) <- sampleNames

+        Cov <- M

+        idxes <- split(seq(along = sampleNames), rep(seq(along = x), times = sapply(x, ncol)))

+        for(ii in seq(along = idxes)) {

+            ov <- findOverlaps(gr, granges(x[[ii]]))

+            idx <- idxes[[ii]]

+            M[queryHits(ov), idx] <- getBSseq(x[[ii]], "M")[subjectHits(ov),]

+            Cov[queryHits(ov), idx] <- getBSseq(x[[ii]], "Cov")[subjectHits(ov),]

+        }

+    }

+    if(any(!sapply(x, hasBeenSmoothed)) || !(all(sameGr))) {

         coef <- NULL

         se.coef <- NULL

         trans <- NULL

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

     lines(xx, yy, col = col, lty = lty, lwd = lwd)

 }

-.bsPlotPoints <- function(x, y, z, col) {

+.bsPlotPoints <- function(x, y, z, col, pointsMinCov) {

     points(x[z>pointsMinCov], y[z>pointsMinCov], col = col, pch = 16, cex = 0.5)

 }

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

 .plotSmoothData <- function(BSseq, region, extend, addRegions, col, lty, lwd, regionCol,

                             addTicks, addPoints, pointsMinCov, highlightMain) {

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

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

     BSseq <- subsetByOverlaps(BSseq, gr)

     ## Extract basic information

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

     if(addPoints) {

         sapply(sampleNames(BSseq), function(samp) {

             .bsPlotPoints(positions, rawPs[, samp], coverage[, samp],

-                          col = colEtc$col[samp])

+                          col = colEtc$col[samp], pointsMinCov = pointsMinCov)

         })

     }

 }

@@ -1,69 +1,64 @@

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

-  ## Argument checking

-  if (anyDuplicated(files)){

-    stop("duplicate entries in 'files'")

-  }

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

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

-  }

-  ## Process each file

-  idxes <- seq_along(files)

-  names(idxes) <- sampleNames

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

-    if (verbose) {

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

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

+    ## Argument checking

+    if (anyDuplicated(files)){

+        stop("duplicate entries in 'files'")

     }

-    stime <- system.time({

-      if (returnRaw) {

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

-      } else{

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

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

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

-        elementMetadata(raw) <- NULL

-        out <- BSseq(gr = raw, M = M, Cov = Cov, sampleNames = sampleNames[ii], rmZeroCov = rmZeroCov)

-      }

-    })[3]

-    if (verbose) {

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

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

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

     }

-    out  

-  })

-  if (!returnRaw) {

+    ## Process each file

+    idxes <- seq_along(files)

+    names(idxes) <- sampleNames

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

+        if (verbose) {

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

+        }

+        stime <- system.time({

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

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

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

+            elementMetadata(raw) <- NULL

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

+                         sampleNames = sampleNames[ii], rmZeroCov = rmZeroCov)

+        })[3]

+        if (verbose) {

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

+        }

+        out  

+    })

     if (verbose) {

-      cat(sprintf("Joining samples ... "))

+        cat(sprintf("Joining samples ... "))

     }

     stime <- system.time({

-      allOut <- Reduce("combine", allOut)

+        allOut <- combineList(allOut)

     })[3]

     if (verbose) {

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

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

     }

-  }

-  allOut

+    allOut

 }

 read.bismarkFileRaw <- function(thisfile, verbose = TRUE){

-  ## Set up the 'what' argument for scan()

-  columnHeaders <- c("chr", "start", "end", "mPerc", "mCount", "uCount")

-  what0 <- replicate(length(columnHeaders), character(0))

-  names(what0) <- columnHeaders

-  int <- which(columnHeaders %in% c("start", "end", "mCount", "uCount"))

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

-  dbl <- which(columnHeaders %in% "mPerc")

-  what0[dbl] <- replicate(length(dbl), double(0))

-  ## Read in the file

-  if (grepl("\\.gz$", thisfile)) 

-    con <- gzfile(thisfile)

-  else 

-    con <- file(thisfile, open = "r")

-  out <- scan(file = con, what = what0, sep = "\t", quote = "", na.strings = "NA", quiet = TRUE)

-  close(con)

-  ## Create GRanges instance from 'out'

-  gr <- GRanges(seqnames = out[["chr"]], ranges = IRanges(start = out[["start"]], width = 1))

-  out[["chr"]] <- out[["start"]] <- out[["end"]] <- NULL

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

-  df <- DataFrame(out)

-  elementMetadata(gr) <- df

-  gr

+    ## Set up the 'what' argument for scan()

+    columnHeaders <- c("chr", "start", "end", "mPerc", "mCount", "uCount")

+    what0 <- replicate(length(columnHeaders), character(0))

+    names(what0) <- columnHeaders

+    int <- which(columnHeaders %in% c("start", "end", "mCount", "uCount"))

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

+    null <- which(columnHeaders %in% "mPerc")

+    what0[null] <- replicate(length(null), NULL)

+    ## Read in the file

+    if (grepl("\\.gz$", thisfile)) 

+        con <- gzfile(thisfile)

+    else 

+        con <- file(thisfile, open = "r")

+    out <- scan(file = con, what = what0, sep = "\t", quote = "", na.strings = "NA", quiet = TRUE)

+    close(con)

+    ## Create GRanges instance from 'out'

+    gr <- GRanges(seqnames = out[["chr"]], ranges = IRanges(start = out[["start"]], width = 1))

+    out[["chr"]] <- out[["start"]] <- out[["end"]] <- NULL

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

+    df <- DataFrame(out)

+    elementMetadata(gr) <- df

+    gr

 }

@@ -401,20 +401,17 @@ read.bsmoothDirRaw <- function(dir, seqnames = NULL, keepCycle = FALSE, keepFilt

 sampleRawToBSseq <- function(gr, qualityCutoff = 20, sampleName = NULL, rmZeroCov = FALSE) {

-    numberQualsGreaterThan.old <- function(char) {

-        quals <- as.integer(charToRaw(char)) - 33L

-        sum(quals >= qualityCutoff)

-    }

     numberQualsGreaterThan <- function(cvec) {

         onestring <- paste(cvec, collapse = "")

         greater <- (as.integer(charToRaw(onestring)) - 33L >= qualityCutoff)

-        tapply(greater, rep(1:length(cvec), times = nchar(cvec)), sum)

+        out <- tapply(greater, rep(1:length(cvec), times = nchar(cvec)), sum)

+        out

     }

     strToCov <- function(vec) {

         Cov <- rep(0, length(vec))

         wh <- which(! vec %in% c("", "0"))

-        ## Cov[wh] <- sapply(vec[wh], numberQualsGreaterThan.old)

-        Cov[wh] <- numberQualsGreaterThan(vec[wh])

+        if(length(wh) > 0)

+            Cov[wh] <- numberQualsGreaterThan(vec[wh])

         Cov

     }

     M <- matrix(strToCov(elementMetadata(gr)[, "Mstr"]), ncol = 1)

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

     if(!returnRaw) {

         if(verbose) cat(sprintf("Joining samples ... "))

         stime <- system.time({

-            ## FIXME: we can do much better than this ... need to write combineList

-            allOut <- Reduce("combine", allOut)

+            allOut <- combineList(allOut)

         })[3]

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

     }

@@ -4,8 +4,15 @@

 \section{Version 0.7.x}{

   \itemize{

-    \item Exposed combineList as a faster alternative to Reduce(combine,

-    list).

+    \item Removed the returnRaw argument to read.bismark() as it was

+    unnecessary (Bismark output files does not have additional

+    information beyond M and Cov and genomic positions, unlike BSmooth).

+    \item Moved the Bismark example data from data to inst/extdata.

+    \item combineList() now deals with the case where the list of BSseq

+    objects have different genomic locations.  This speeds up

+    read.bismark() substantially.

+    \item Exposed combineList() as a faster alternative to

+    Reduce(combine, list).

     \item Updated the code for the plotting routines (plotRegion).  This

     should not have an impact on user-visible code.

     \item Added read.bismark() function to parse output from the Bismark

similarity index 100%

rename from data/CpG_context_test_data.fastq_bismark.bedGraph

rename to inst/extdata/CpG_context_test_data.fastq_bismark.bedGraph

@@ -7,51 +7,66 @@

   Parsing output from the Bismark alignment suite.

 }

 \usage{

-  read.bismark(files, sampleNames, returnRaw = FALSE, rmZeroCov = FALSE, verbose = TRUE)

+  read.bismark(files, sampleNames, rmZeroCov = FALSE, verbose = TRUE)

 }

 \arguments{

-  \item{files}{Input files. Each sample is in a different file. Input files are created by running Bismark's \code{methylation_extractor}; see Note for details.} 

+  \item{files}{Input files. Each sample is in a different file. Input

+    files are created by running Bismark's \code{methylation_extractor};

+    see Note for details.}  

   \item{sampleNames}{sample names, based on the order of \code{files}.}

-  \item{returnRaw}{Should the function return the complete information in the output files?}

-  \item{rmZeroCov}{Should methylation loci that have zero coverage in all samples be removed. This will result in a much smaller object if the data originates from (targeted) capture bisulfite sequencing.} 

+  \item{rmZeroCov}{Should methylation loci that have zero coverage in

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

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

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

 }

 \note{

   Input files can either be gzipped or not.

-  Input files are created by running Bismark's \code{methylation_extractor} and \code{genome_methylation_bismark2bedGraph_v4.pl} scripts over the Bismark alignment file. For example (run from the command line):

+  Input files are created by running Bismark's

+  \code{methylation_extractor} and

+  \code{genome_methylation_bismark2bedGraph_v4.pl} scripts over the

+  Bismark alignment file. For example (run from the command line): 

-  \code{methylation_extractor -s --comprehensive test_data.fastq_bismark.sam}

+  \code{methylation_extractor -s --comprehensive

+  test_data.fastq_bismark.sam}

-  \code{genome_methylation_bismark2bedGraph_v4.pl --counts CpG_context_test_data.fastq_bismark.txt > CpG_context_test_data.fastq_bismark.bedGraph}

+  \code{genome_methylation_bismark2bedGraph_v4.pl --counts

+  CpG_context_test_data.fastq_bismark.txt >

+  CpG_context_test_data.fastq_bismark.bedGraph} 

-  The \code{--comprehensive} argument to \code{methylation_extractor} and the \code{--counts} argument to \code{genome_methylation_bismark2bedGraph_v4.pl} are required.

+  The \code{--comprehensive} argument to \code{methylation_extractor}

+  and the \code{--counts} argument to

+  \code{genome_methylation_bismark2bedGraph_v4.pl} are required. 

-  In this example, the file \code{CpG_context_test_data.fastq_bismark.bedGraph} is then the input file to \code{read.bismark}. 

+  In this example, the file

+  \code{CpG_context_test_data.fastq_bismark.bedGraph} is then the input

+  file to \code{read.bismark}.  

-  See \url{http://rpubs.com/PeteHaitch/readBismark} for a worked example using Bismark and \code{read.bismark}. Please consult the Bismark website for full details of these scripts and the latest versions (\url{http://www.bioinformatics.babraham.ac.uk/projects/download.html#bismark})

+  See \url{http://rpubs.com/PeteHaitch/readBismark} for a worked example

+  using Bismark and \code{read.bismark}. Please consult the Bismark

+  website for full details of these scripts and the latest versions

+  (\url{http://www.bioinformatics.babraham.ac.uk/projects/download.html#bismark}) 

 }

 \value{

-  Either an object of class \code{"BSseq"} (if \code{returnRaw = FALSE})

-  or a list of \code{"GRanges"} which each component coming from a

-  file. 

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

 }

 \seealso{

-  \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 collapse

-  (merging or summing) the data in two different directories.

+  \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 collapse  (merging or summing) the

+  data in two different directories. 

 }

 \examples{

   \dontrun{

-  bismarkBedGraph <- system.file('data/CpG_context_test_data.fastq_bismark.bedGraph', package = 'bsseq')

-  bismarkBSseq <- read.bismark(files = bismarkBedGraph, sampleNames = 'test_data', returnRaw = FALSE, rmZeroCov = FALSE, verbose = TRUE)

+  bismarkBedGraph <- system.file("extdata/CpG_context_test_data.fastq_bismark.bedGraph", package = 'bsseq')

+  bismarkBSseq <- read.bismark(files = bismarkBedGraph, sampleNames = "test_data", rmZeroCov = FALSE, verbose = TRUE)

   bismarkBSseq

   }

 }

 \author{

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

-

 }