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

 Package: bsseq

-Version: 0.9.9

+Version: 0.9.10

 Title: Analyze, manage and store bisulfite sequencing data

 Description: Tools for analyzing and visualizing bisulfite sequencing data

 Author: Kasper Daniel Hansen

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

         stopifnot(min(group2) >= 1 & max(group2) <= ncol(BSseq))

     } else stop("problems with argument 'group2'")

     stopifnot(length(intersect(group1, group2)) == 0)

+    stopifnot(length(group1) > 0)

+    stopifnot(length(group2) > 0)

+    stopifnot(length(group1) + length(group2) >= 3)

     if(estimate.var == "paired")

         stopifnot(length(group1) == length(group2))

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

     if(is.null(maxGap))

         maxGap <- BSseq@parameters$maxGap

+    if(is.null(maxGap))

+        stop("need to set argument 'maxGap'")

     if(verbose) cat("[BSmooth.tstat] preprocessing ... ")

     ptime1 <- proc.time()

@@ -65,20 +65,21 @@ plot.BSseqTstat <- function(x, y, ...) {

     }

 }

-getStats <- function(BSseqTstat, regions = NULL, column = "tstat.corrected") {

+getStats <- function(BSseqTstat, regions = NULL, stat = "tstat.corrected") {

     stopifnot(is(BSseqTstat, "BSseqTstat"))

-    stopifnot(column %in% colnames(BSseqTstat@stats))

-    if(class(regions) == "data.frame")

-        regions <- data.frame2GRanges(regions)

     if(is.null(regions))

         return(BSseqTstat@stats)

+    if(class(regions) == "data.frame")

+        regions <- data.frame2GRanges(regions)

+    stopifnot(stat %in% colnames(BSseqTstat@stats))

+    stopifnot(length(stat) == 1)

     stopifnot(is(regions, "GenomicRanges"))

     ov <- findOverlaps(BSseqTstat, regions)

     ov.sp <- split(queryHits(ov), subjectHits(ov))

     getRegionStats <- function(idx) {

         mat <- BSseqTstat@stats[idx,, drop=FALSE]

-        areaStat <- sum(mat[, column])

-        maxStat <- max(mat[, column])

+        areaStat <- sum(mat[, stat])

+        maxStat <- max(mat[, stat])

         c(areaStat, maxStat)

     }

     stats <- matrix(NA, ncol = 2, nrow = length(regions))

@@ -86,7 +87,7 @@ getStats <- function(BSseqTstat, regions = NULL, column = "tstat.corrected") {

     tmp <- lapply(ov.sp, getRegionStats)

     stats[as.integer(names(tmp)),] <- do.call(rbind, tmp)

     out <- as.data.frame(stats)

-    if(! column %in% c("tstat.corrected", "tstat"))

+    if(! stat %in% c("tstat.corrected", "tstat"))

         return(out)

     getRegionStats_ttest <- function(idx) {

         mat <- BSseqTstat@stats[idx,, drop=FALSE]

@@ -14,12 +14,14 @@ dmrFinder <- function(BSseqTstat, cutoff = NULL, qcutoff = c(0.025, 0.975),

     positions <- start(BSseqTstat)

     regions <- bsseq:::regionFinder3(direction, chr = chrs, pos = positions,

                                      maxGap = maxGap, verbose = subverbose)

+    if(is.null(regions$down) && is.null(regions$up))

+        return(NULL)

     if(verbose) cat("[dmrFinder] creating dmr data.frame\n")

     regions <- do.call(rbind, regions)

     rownames(regions) <- NULL

     regions$width <- regions$end - regions$start + 1

     regions$invdensity <- regions$width / regions$n

-    stats <- getStats(BSseqTstat, regions, column = stat)

+    stats <- getStats(BSseqTstat, regions, stat = stat)

     regions <- cbind(regions, stats)

     if(stat %in% c("tstat.corrected", "tstat")) {

         regions$direction <- ifelse(regions$meanDiff > 0, "hyper", "hypo")

@@ -21,6 +21,7 @@

     \item Allowed user specification of maxk, with a default of 10,000

     in BSmooth.

     \item Many bugfixes made necessary by the new class representation.

+    \item Better argument checking in BSmooth.tstat.

   }

 }