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

 Package: bsseq

-Version: 0.9.7

+Version: 0.9.8

 Title: Analyze, manage and store bisulfite sequencing data

 Description: Tools for analyzing and visualizing bisulfite sequencing data

 Author: Kasper Daniel Hansen

@@ -9,7 +9,7 @@ Depends: R (>= 2.15), methods, BiocGenerics, IRanges, GenomicRanges,

 Imports: scales, stats, graphics, Biobase, locfit

 Suggests: RUnit, bsseqData

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

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

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

         gof_stats.R plotting.R fisher.R fixes.R

 License: Artistic-2.0

 LazyData: yes

@@ -148,128 +148,3 @@ BSmooth <- function(BSseq, ns = 70, h = 1000, maxGap = 10^8, parallelBy = c("sam

 }

-BSmooth.tstat <- function(BSseq, group1, group2, estimate.var = c("same", "paired", "group2"),

-                          local.correct = TRUE, maxGap = NULL, qSd = 0.75, k = 101, mc.cores = 1, verbose = TRUE){

-    smoothSd <- function(Sds, k) {

-        k0 <- floor(k/2)

-        if(all(is.na(Sds))) return(Sds)

-        thresSD <- pmax(Sds, quantile(Sds, qSd, na.rm = TRUE), na.rm = TRUE)

-        addSD <- rep(median(Sds, na.rm = TRUE), k0)

-        sSds <- as.vector(runmean(Rle(c(addSD, thresSD, addSD)), k = k))

-        sSds

-    }

-    compute.correction <- function(idx, qSd = 0.75) {

-        xx <- start(BSseq)[idx]

-        yy <- tstat[idx]

-        suppressWarnings({

-            drange <- diff(range(xx, na.rm = TRUE))

-        })

-        if(drange <= 25000)

-            return(yy)

-        tstat.function <- approxfun(xx, yy)

-        xx.reg <- seq(from = min(xx), to = max(xx), by = 2000)

-        yy.reg <- tstat.function(xx.reg)

-        fit <- locfit(yy.reg ~ lp(xx.reg, h = 25000, deg = 2, nn = 0),

-                      family = "huber", maxk = 50000) 

-        correction <- predict(fit, newdata = data.frame(xx.reg = xx))

-        yy - correction 

-    }

-

-    estimate.var <- match.arg(estimate.var)

-    stopifnot(is(BSseq, "BSseq"))

-    stopifnot(hasBeenSmoothed(BSseq))

-    if(is.numeric(group1)) {

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

-        group1 <- sampleNames(BSseq)[group1]

-    }

-    if(is.numeric(group2)) {

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

-        group2 <- sampleNames(BSseq)[group2]

-    }

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

-    if(estimate.var == "paired")

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

-    stopifnot(all(c(group1, group2) %in% sampleNames(BSseq)))

-

-    if(any(rowSums(getCoverage(BSseq)[, c(group1, group2)]) == 0))

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

-    

-    if(is.null(maxGap))

-        maxGap <- BSseq@parameters$maxGap

-    

-    if(verbose) cat("preprocessing ... ")

-    stime <- system.time({

-        clusterIdx <- makeClusters(BSseq, maxGap = maxGap, mc.cores = mc.cores)

-    })[3]

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

-        

-    if(verbose) cat("computing stats within groups ... ")

-    stime <- system.time({

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

-                         confint = FALSE)[, c(group1, group2)]

-        group1.means <- rowMeans(allPs[, group1, drop = FALSE], na.rm = TRUE)

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

-        })[3]                               

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

-    

-    if(verbose) cat("computing stats across groups ... ")

-    stime <- system.time({

-        switch(estimate.var,

-               "group2" = {

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

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

-                       smoothSd(rawSds[idx], k = k)

-                   }, mc.cores = mc.cores))

-                   scale <- sqrt(1/length(group1) + 1/length(group2))

-                   tstat.sd <- smoothSds * scale

-               },

-               "same" = {

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

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

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

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

-                       smoothSd(rawSds[idx], k = k)

-                   }, mc.cores = mc.cores))

-                   scale <- sqrt(1/length(group1) + 1/length(group2))

-                   tstat.sd <- smoothSds * scale

-               },

-               "paired" = {

-                   rawSds <- rowSds(allPs[, group1, drop = FALSE] - allPs[, group2, drop = FALSE])

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

-                       smoothSd(rawSds[idx], k = k)

-                   }, mc.cores = mc.cores))

-                   scale <- sqrt(1/length(group1))

-                   tstat.sd <- smoothSds * scale

-               })

-        tstat <- (group1.means - group2.means) / tstat.sd

-        is.na(tstat)[tstat.sd == 0] <- TRUE

-        if(local.correct) {

-            tstat.corrected <- do.call(c, mclapply(clusterIdx,

-                                                   compute.correction, qSd = qSd,

-                                                   mc.cores = mc.cores))

-        }

-    })[3]

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

-    

-    if(local.correct) {

-        stats <- cbind(rawSds, tstat.sd, group2.means, group1.means,

-                       tstat, tstat.corrected)

-        colnames(stats) <- c("rawSds", "tstat.sd",

-                             "group2.means", "group1.means", "tstat",

-                             "tstat.corrected")

- 

-    } else {

-        stats <- cbind(rawSds, tstat.sd, group2.means, group1.means,

-                       tstat)

-        colnames(stats) <- c("rawSds", "tstat.sd",

-                             "group2.means", "group1.means", "tstat")

-    }

-    

-    parameters <- c(BSseq@parameters,

-                    list(tstatText = sprintf("BSmooth.tstat (local.correct = %s, maxGap = %d)",

-                         local.correct, maxGap),

-                         group1 = group1, group2 = group2, k = k, qSd = qSd,

-                         local.correct = local.correct, maxGap = maxGap))

-    out <- BSseqTstat(gr = granges(BSseq), stats = stats, parameters = parameters)

-    out

-}

new file mode 100644

@@ -0,0 +1,130 @@

+BSmooth.tstat <- function(BSseq, group1, group2, estimate.var = c("same", "paired", "group2"),

+                          local.correct = TRUE, maxGap = NULL, qSd = 0.75, k = 101, mc.cores = 1, verbose = TRUE){

+    smoothSd <- function(Sds, k) {

+        k0 <- floor(k/2)

+        if(all(is.na(Sds))) return(Sds)

+        thresSD <- pmax(Sds, quantile(Sds, qSd, na.rm = TRUE), na.rm = TRUE)

+        addSD <- rep(median(Sds, na.rm = TRUE), k0)

+        sSds <- as.vector(runmean(Rle(c(addSD, thresSD, addSD)), k = k))

+        sSds

+    }

+    compute.correction <- function(idx, qSd = 0.75) {

+        xx <- start(BSseq)[idx]

+        yy <- tstat[idx]

+        suppressWarnings({

+            drange <- diff(range(xx, na.rm = TRUE))

+        })

+        if(drange <= 25000)

+            return(yy)

+        tstat.function <- approxfun(xx, yy)

+        xx.reg <- seq(from = min(xx), to = max(xx), by = 2000)

+        yy.reg <- tstat.function(xx.reg)

+        fit <- locfit(yy.reg ~ lp(xx.reg, h = 25000, deg = 2, nn = 0),

+                      family = "huber", maxk = 50000) 

+        correction <- predict(fit, newdata = data.frame(xx.reg = xx))

+        yy - correction 

+    }

+

+    estimate.var <- match.arg(estimate.var)

+    stopifnot(is(BSseq, "BSseq"))

+    stopifnot(hasBeenSmoothed(BSseq))

+    if(is.character(group1)) {

+        stopifnot(all(group1 %in% sampleNames(BSseq)))

+        group1 <- match(group1, sampleNames(BSseq))

+    }

+    if(is.numeric(group1)) {

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

+    } else stop("problems with argument 'group1'")

+    if(is.character(group2)) {

+        stopifnot(all(group2 %in% sampleNames(BSseq)))

+        group2 <- match(group2, sampleNames(BSseq))

+    }    

+    if(is.numeric(group2)) {

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

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

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

+    if(estimate.var == "paired")

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

+    

+    if(any(rowSums(getCoverage(BSseq)[, c(group1, group2)]) == 0))

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

+    

+    if(is.null(maxGap))

+        maxGap <- BSseq@parameters$maxGap

+    

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

+    stime <- system.time({

+        clusterIdx <- makeClusters(BSseq, maxGap = maxGap, mc.cores = mc.cores)

+    })[3]

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

+        

+    if(verbose) cat("[BSmooth.tstat] computing stats within groups ... ")

+    stime <- system.time({

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

+                         confint = FALSE)[, c(group1, group2)]

+        group1.means <- rowMeans(allPs[, group1, drop = FALSE], na.rm = TRUE)

+        group2.means <- rowMeans(allPs[, group2, drop = FALSE], na.rm = TRUE)

+        })[3]                               

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

+    

+    if(verbose) cat("[BSmooth.tstat] computing stats across groups ... ")

+    stime <- system.time({

+        switch(estimate.var,

+               "group2" = {

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

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

+                       smoothSd(rawSds[idx], k = k)

+                   }, mc.cores = mc.cores))

+                   scale <- sqrt(1/length(group1) + 1/length(group2))

+                   tstat.sd <- smoothSds * scale

+               },

+               "same" = {

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

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

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

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

+                       smoothSd(rawSds[idx], k = k)

+                   }, mc.cores = mc.cores))

+                   scale <- sqrt(1/length(group1) + 1/length(group2))

+                   tstat.sd <- smoothSds * scale

+               },

+               "paired" = {

+                   rawSds <- rowSds(allPs[, group1, drop = FALSE] - allPs[, group2, drop = FALSE])

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

+                       smoothSd(rawSds[idx], k = k)

+                   }, mc.cores = mc.cores))

+                   scale <- sqrt(1/length(group1))

+                   tstat.sd <- smoothSds * scale

+               })

+        tstat <- (group1.means - group2.means) / tstat.sd

+        is.na(tstat)[tstat.sd == 0] <- TRUE

+        if(local.correct) {

+            tstat.corrected <- do.call(c, mclapply(clusterIdx,

+                                                   compute.correction, qSd = qSd,

+                                                   mc.cores = mc.cores))

+        }

+    })[3]

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

+    

+    if(local.correct) {

+        stats <- cbind(rawSds, tstat.sd, group2.means, group1.means,

+                       tstat, tstat.corrected)

+        colnames(stats) <- c("rawSds", "tstat.sd",

+                             "group2.means", "group1.means", "tstat",

+                             "tstat.corrected")

+ 

+    } else {

+        stats <- cbind(rawSds, tstat.sd, group2.means, group1.means,

+                       tstat)

+        colnames(stats) <- c("rawSds", "tstat.sd",

+                             "group2.means", "group1.means", "tstat")

+    }

+    

+    parameters <- c(BSseq@parameters,

+                    list(tstatText = sprintf("BSmooth.tstat (local.correct = %s, maxGap = %d)",

+                         local.correct, maxGap),

+                         group1 = group1, group2 = group2, k = k, qSd = qSd,

+                         local.correct = local.correct, maxGap = maxGap))

+    out <- BSseqTstat(gr = granges(BSseq), stats = stats, parameters = parameters)

+    out

+}

@@ -37,3 +37,69 @@ BSseqTstat <- function(gr = NULL, stats = NULL, parameters = NULL) {

     out

 }

+summary.BSseqTstat <- function(object, ...) {

+    quant <- quantile(getStats(object)[, "tstat.corrected"],

+                      prob = c(0.0001, 0.001, 0.01, 0.5, 0.99, 0.999, 0.9999))

+    quant <- t(t(quant))

+    colnames(quant) <- "quantiles"

+    out <- list(quantiles = quant)

+    class(out) <- "summary.BSseqTstat"

+    out

+}

+

+print.summary.BSseqTstat <- function(x, ...) {

+    print(as.matrix(x$quantiles))

+}

+

+plot.BSseqTstat <- function(x, y, ...) {

+    tstat <- getStats(x)[, "tstat"]

+    plot(density(tstat), xlim = c(-10,10), col = "blue", main = "")

+    if("tstat.corrected" %in% colnames(getStats(x))) {

+        tstat.cor <- getStats(x)[, "tstat.corrected"]

+        lines(density(tstat.cor), col = "black")

+        legend("topleft", legend = c("uncorrected", "corrected"), lty = c(1,1),

+               col = c("blue", "black"))

+    } else {

+        legend("topleft", legend = c("uncorrected"), lty = 1,

+               col = c("blue"))

+    }

+}

+

+getStats <- function(BSseqTstat, regions = NULL, column = "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)

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

+        c(areaStat, maxStat)

+    }

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

+    colnames(stats) <- c("areaStat", "maxStat")

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

+        return(out)

+    getRegionStats_ttest <- function(idx) {

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

+        group1.mean <- mean(mat[, "group1.means"])

+        group2.mean <- mean(mat[, "group2.means"])

+        meanDiff <- mean(mat[, "group1.means"] - mat[, "group2.means"])

+        tstat.sd <- mean(mat[, "tstat.sd"])

+        c(meanDiff, group1.mean, group2.mean, tstat.sd)

+    }

+    stats_ttest<- matrix(NA, ncol = 4, nrow = length(regions))

+    colnames(stats_ttest) <- c("meanDiff", "group1.mean", "group2.mean", "tstat.sd")

+    tmp <- lapply(ov.sp, getRegionStats)

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

+    out <- cbind(out, as.data.frame(stats_ttest))

+    out

+}

@@ -1,23 +1,22 @@

 collapseBSseq <- function(BSseq, columns) {

     ## columns is a vector of new names, names(columns) is sampleNames or empty

+    stopifnot(is.character(columns))

     if(is.null(names(columns)) && length(columns) != ncol(BSseq))

         stop("if `columns' does not have names, it needs to be of the same length as `BSseq` has columns (samples)")

     if(!is.null(names(columns)) && !all(names(columns) %in% sampleNames(BSseq)))

         stop("if `columns` has names, they need to be sampleNames(BSseq)")

-    if(is.null(names(columns))) {

-        sp <- split(1:ncol(BSseq), columns)

-    } else {

-        sp <- split(names(columns), columns)

-    }

+    if(is.null(names(columns)))

+        columns.idx <- 1:ncol(BSseq)

+    else

+        columns.idx <- match(names(columns), sampleNames(BSseq))

+    sp <- split(columns.idx, columns)

     M <- do.call(cbind, lapply(sp, function(ss) {

         rowSums(getBSseq(BSseq, "M")[, ss, drop = FALSE])

     }))

-    colnames(M) <- names(sp)

     Cov <- do.call(cbind, lapply(sp, function(ss) {

         rowSums(getBSseq(BSseq, "M")[, ss, drop = FALSE])

     }))

-    colnames(Cov) <- names(sp)

-    BSseq(gr = getBSseq(BSseq, "gr"), M = M, Cov = Cov)

+    BSseq(gr = getBSseq(BSseq, "gr"), M = M, Cov = Cov, sampleNames = names(sp))

 }

 chrSelectBSseq <- function(BSseq, seqnames = NULL, order = FALSE) {

@@ -192,59 +191,4 @@ getCoverage <- function(BSseq, regions = NULL, type = c("Cov", "M"),

     return(outMatrix)

 }

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

-    stopifnot(is(BSseqTstat, "BSseqTstat"))

-    column <- match.arg(column)

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

-        regions <- data.frame2GRanges(regions)

-    if(is.null(regions))

-        return(BSseqTstat@stats)

-    stopifnot(is(regions, "GenomicRanges"))

-    ov <- findOverlaps(BSseqTstat, regions)

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

-    getRegionStats <- function(idx) {

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

-        meanDiff <- mean(mat[, "group1.means"] - mat[, "group2.means"])

-        group1.mean <- mean(mat[, "group1.means"])

-        group2.mean <- mean(mat[, "group2.means"])

-        tstat.sd <- mean(mat[, "tstat.sd"])

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

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

-        c(meanDiff, group1.mean, group2.mean, tstat.sd, areaStat, maxStat)

-    }

-    stats <- lapply(ov.sp, getRegionStats)

-    out <- matrix(NA, ncol = 6, nrow = length(regions))

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

-    colnames(out) <- c("meanDiff", "group1.mean", "group2.mean", "tstat.sd", "areaStat", "maxStat")

-    out <- as.data.frame(out)

-    out

-}

-

-

-summary.BSseqTstat <- function(object, ...) {

-    quant <- quantile(getStats(object)[, "tstat.corrected"],

-                      prob = c(0.0001, 0.001, 0.01, 0.5, 0.99, 0.999, 0.9999))

-    quant <- t(t(quant))

-    colnames(quant) <- "quantiles"

-    out <- list(quantiles = quant)

-    class(out) <- "summary.BSseqTstat"

-    out

-}

-print.summary.BSseqTstat <- function(x, ...) {

-    print(as.matrix(x$quantiles))

-}

-

-plot.BSseqTstat <- function(x, y, ...) {

-    tstat <- getStats(x)[, "tstat"]

-    plot(density(tstat), xlim = c(-10,10), col = "blue", main = "")

-    if("tstat.corrected" %in% colnames(getStats(x))) {

-        tstat.cor <- getStats(x)[, "tstat.corrected"]

-        lines(density(tstat.cor), col = "black")

-        legend("topleft", legend = c("uncorrected", "corrected"), lty = c(1,1),

-               col = c("blue", "black"))

-    } else {

-        legend("topleft", legend = c("uncorrected"), lty = 1,

-               col = c("blue"))

-    }

-}

@@ -1,24 +1,29 @@

 dmrFinder <- function(BSseqTstat, cutoff = NULL, qcutoff = c(0.025, 0.975),

-                             maxGap = 300, column = c("tstat.corrected", "tstat"), verbose = TRUE) {

-    column <- match.arg(column)

-    if(! column %in% colnames(BSseqTstat@stats))

-        stop("'column' needs to be a column of 'BSseqTstat@stats'")

+                             maxGap = 300, stat = "tstat.corrected", verbose = TRUE) {

+    subverbose <- max(as.integer(verbose) - 1L, 0L)

+    if(! stat %in% colnames(BSseqTstat@stats))

+        stop("'stat' needs to be a column of 'BSseqTstat@stats'")

     if(is.null(cutoff))

-        cutoff <- quantile(BSseqTstat@stats[, column], qcutoff)

-    direction <- as.integer(BSseqTstat@stats[, column] >= cutoff[2])

-    direction[BSseqTstat@stats[, column] <= cutoff[1]] <- -1L

+        cutoff <- quantile(BSseqTstat@stats[, stat], qcutoff)

+    if(length(cutoff) == 1)

+        cutoff <- c(-cutoff, cutoff)

+    direction <- as.integer(BSseqTstat@stats[, stat] >= cutoff[2])

+    direction[BSseqTstat@stats[, stat] <= cutoff[1]] <- -1L

     direction[is.na(direction)] <- 0L

     chrs <- as.character(seqnames(BSseqTstat@gr))

     positions <- start(BSseqTstat)

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

-    if(verbose) cat("creating dmr data.frame\n")

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

+                                     maxGap = maxGap, verbose = subverbose)

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

     regions <- do.call(rbind, regions)

     rownames(regions) <- NULL

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

-    regions <- cbind(regions, stats)

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

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

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

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

+    regions <- cbind(regions, stats)

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

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

+    }

     regions <- regions[order(abs(regions$areaStat), decreasing = TRUE),]

     regions

 }

@@ -46,11 +51,11 @@ clusterMaker <- function(chr, pos, order.it=TRUE, maxGap=300){

 regionFinder3 <- function(x, chr, positions, keep, maxGap = 300, verbose = TRUE) {

     getSegments3 <- function(x, cid, verbose = TRUE){

-        if(verbose) cat("segmenting\n")

+        if(verbose) cat("[regionFinder3] segmenting\n")

         segments <- cumsum(c(1, diff(x) != 0)) +

             cumsum(c(1, diff(cid) != 0))

         names(segments) <- NULL

-        if(verbose) cat("splitting\n")

+        if(verbose) cat("[regionFinder3] splitting\n")

         out <- list(up = split(which(x > 0), segments[x > 0]),

                     down = split(which(x < 0), segments[x < 0]),

                     nochange = split(which(x == 0), segments[x == 0]))

@@ -2,16 +2,21 @@ fisherTests <- function(BSseq, group1, group2, lookup = NULL,

                         returnLookup = TRUE, mc.cores = 1,

                         verbose = TRUE) {

     stopifnot(is(BSseq, "BSseq"))

-    if(is.integer(group1)) {

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

-        group1 <- sampleNames(BSseq)[group1]

+    if(is.character(group1)) {

+        stopifnot(all(group1 %in% sampleNames(BSseq)))

+        group1 <- match(group1, sampleNames(BSseq))

     }

-    if(is.integer(group2)) {

+    if(is.numeric(group1)) {

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

+    } else stop("problems with argument 'group1'")

+    if(is.character(group2)) {

+        stopifnot(all(group2 %in% sampleNames(BSseq)))

+        group2 <- match(group2, sampleNames(BSseq))

+    }    

+    if(is.numeric(group2)) {

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

-        group2 <- sampleNames(BSseq)[group2]

-    }

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

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

-    sampleNames <- c(group1, group2)

     if(is.null(lookup)) {

         oldkeys <- ""

     } else {

@@ -19,7 +24,7 @@ fisherTests <- function(BSseq, group1, group2, lookup = NULL,

                   setequal(colnames(lookup), c("p.value", "log2OR")))

         oldkeys <- rownames(lookup)

     }

-    if(verbose) cat(" setup ... ")

+    if(verbose) cat("[fisherTests] setup ... ")

     stime <- system.time({

         M1 <- rowSums(getCoverage(BSseq, type = "M")[, group1, drop = FALSE])

         M2 <- rowSums(getCoverage(BSseq, type = "M")[, group2, drop = FALSE])

@@ -34,16 +39,16 @@ fisherTests <- function(BSseq, group1, group2, lookup = NULL,

         expand2[,"U1"] <- expand[,"Cov1"] - expand[,"M1"]

         expand2[,"U2"] <- expand[,"Cov2"] - expand[,"M2"]

     })[3]

-    if(verbose) cat(round(stime,1), "secs\n")

-    if(verbose) cat(sprintf(" performing %d tests (using %d cores) ... ",

+    if(verbose) cat(round(stime,1), "done in secs\n")

+    if(verbose) cat(sprintf("[fisherTests] performing %d tests (using %d cores) ... ",

                             nrow(expand2), mc.cores))

     stime <- system.time({

         tests <- mclapply(1:nrow(expand2), function(ii) {

             unlist(fisher.test(matrix(expand2[ii,], ncol = 2))[c("p.value", "estimate")])

         }, mc.cores = mc.cores)

     })[3]

-    if(verbose) cat(round(stime,1), "secs\n")

-    if(verbose) cat(" finalizing ... ")

+    if(verbose) cat(round(stime,1), "done in secs\n")

+    if(verbose) cat("[fisherTests] finalizing ... ")

     stime <- system.time({

         tests <- do.call(rbind, tests)

         tests[,2] <- log2(tests[,2])

@@ -55,7 +60,7 @@ fisherTests <- function(BSseq, group1, group2, lookup = NULL,

     })[3]

-    if(verbose) cat(round(stime,1), "secs\n")

+    if(verbose) cat(round(stime,1), "done in secs\n")

     if(returnLookup)

         return(list(results = results, lookup = lookup))

     else

@@ -14,10 +14,12 @@ plotAnnoTrack <- function(gr, annoTrack) {

 }

 plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addRegions = NULL,

-                            annoTrack = NULL, col = NULL, lty = NULL, lwd = NULL, BSseqTstat = NULL,

+                            annoTrack = NULL, col = NULL, lty = NULL, lwd = NULL, 

+                            BSseqTstat = NULL, stat = "tstat.corrected", stat.col = "black",

+                            stat.lwd = 1, stat.lty = 1, stat.ylim = c(-8,8),

                             mainWithWidth = TRUE, regionCol = alpha("red", 0.1), addTicks = TRUE,

                             addPoints = FALSE, pointsMinCov = 5, highlightMain = FALSE, verbose = TRUE) {

-    cat("preprocessing ...")

+    cat("[plotManyRegions] preprocessing ...")

     if(!is.null(regions)) {

         if(is(regions, "data.frame"))

             gr <- data.frame2GRanges(regions, keepColumns = FALSE)

@@ -39,9 +41,11 @@ plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addReg

         main <- rep(main, length = length(gr))

     cat("done\n")

     for(ii in seq(along = gr)) {

-        if(verbose) cat(sprintf("plotting region %d (out of %d)\n", ii, nrow(regions)))

+        if(verbose) cat(sprintf("[plotManyRegions]   plotting region %d (out of %d)\n", ii, nrow(regions)))

         plotRegion(BSseq = BSseq, region = regions[ii,], extend = extend,

                    col = col, lty = lty, lwd = lwd, main = main[ii], BSseqTstat = BSseqTstat,

+                   stat = stat, stat.col = stat.col, stat.lwd = stat.lwd,

+                   stat.lty = stat.lty, stat.ylim = stat.ylim,

                    addRegions = addRegions, regionCol = regionCol, mainWithWidth = mainWithWidth,

                    annoTrack = annoTrack, addTicks = addTicks, addPoints = addPoints,

                    pointsMinCov = pointsMinCov, highlightMain = highlightMain)

@@ -173,30 +177,33 @@ plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addReg

                         regionCol = regionCol, highlightMain = highlightMain)

     if(addPoints) {

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

-            abline(v = positions[rawPs[, samp] > 0.1], col = "grey80", lty = 1)

+        sapply(1:ncol(BSseq), function(sampIdx) {

+            abline(v = positions[rawPs[, sampIdx] > 0.1], col = "grey80", lty = 1)

         })

     } # This adds vertical grey lines so we can see where points are plotted

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

-        .bsPlotLines(positions, smoothPs[, samp], col = colEtc$col[samp],

-                     lty = colEtc$lty[samp], lwd = colEtc$lwd[samp],

+    sapply(1:ncol(BSseq), function(sampIdx) {

+        .bsPlotLines(positions, smoothPs[, sampIdx], col = colEtc$col[sampIdx],

+                     lty = colEtc$lty[sampIdx], lwd = colEtc$lwd[sampIdx],

                      plotRange = c(start(gr), end(gr)))

     })

     if(addPoints) {

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

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

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

+        sapply(1:ncol(BSseq), function(sampIdx) {

+            .bsPlotPoints(positions, rawPs[, sampIdx], coverage[, sampIdx],

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

         })

     }

 }

 plotRegion <- function(BSseq, region = NULL, extend = 0, main = "", addRegions = NULL, annoTrack = NULL,

-                          col = NULL, lty = NULL, lwd = NULL, BSseqTstat = NULL, mainWithWidth = TRUE,

-                          regionCol = alpha("red", 0.1), addTicks = TRUE, addPoints = FALSE,

-                          pointsMinCov = 5, highlightMain = FALSE) {

+                       col = NULL, lty = NULL, lwd = NULL,

+                       BSseqTstat = NULL, stat = "tstat.corrected", stat.col = "black",

+                       stat.lwd = 1, stat.lty = 1, stat.ylim = c(-8,8),

+                       mainWithWidth = TRUE,

+                       regionCol = alpha("red", 0.1), addTicks = TRUE, addPoints = FALSE,

+                       pointsMinCov = 5, highlightMain = FALSE) {

     opar <- par(mar = c(0,4.1,0,0), oma = c(5,0,4,2), mfrow = c(1,1))

     on.exit(par(opar))

@@ -214,16 +221,14 @@ plotRegion <- function(BSseq, region = NULL, extend = 0, main = "", addRegions =

     if(!is.null(BSseqTstat)) {

         if(!is.null(BSseqTstat))

             BSseqTstat <- subsetByOverlaps(BSseqTstat, gr)

-        plot(positions[1], 0.5, type = "n", xaxt = "n", yaxt = "n",

-             ylim = c(-8,8), xlim = plotRange, xlab = "", ylab = "t-stat")

+        plot(start(gr), 0.5, type = "n", xaxt = "n", yaxt = "n",

+             ylim = stat.ylim, xlim = c(start(gr), end(gr)), xlab = "", ylab = "t-stat")

         axis(side = 2, at = c(-5,0,5))

         abline(h = 0, col = "grey60")

-        bsseq:::.bsPlotLines(start(BSseqTstat), BSseqTstat@stats[, "tstat"],

-                             lty = 1, plotRange = plotRange, col = "red", lwd = 1)

-        bsseq:::.bsPlotLines(start(BSseqTstat), BSseqTstat@stats[, "tstat.corrected"],

-                             lty = 2, plotRange = plotRange, col = "red", lwd = 1)

-        bsseq:::.bsPlotLines(start(BSseqTstat), 100*BSseqTstat@stats[, "tstat.sd"],

-                             lty = 2, plotRange = plotRange, col = "blue", lwd = 1)

+        mapply(function(stat, col, lty, lwd) {

+            bsseq:::.bsPlotLines(start(BSseqTstat), BSseqTstat@stats[, stat],

+                                 lty = lty, plotRange = c(start(gr), end(gr)), col = col, lwd = lwd)

+        }, stat = stat, col = stat.col, lty = stat.lty, lwd = stat.lwd)

     }

     if(!is.null(annoTrack))

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

     names(idxes) <- sampleNames

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

         if (verbose) {

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

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

         }

         stime <- system.time({

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

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

                          sampleNames = sampleNames[ii], rmZeroCov = rmZeroCov)

         })[3]

         if (verbose) {

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

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

         }

         out  

     })

     if (verbose) {

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

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

     }

     stime <- system.time({

         allOut <- combineList(allOut)

     })[3]

     if (verbose) {

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

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

     }

     allOut

 }

@@ -79,7 +79,7 @@ read.umtab2.chr <- function(files, sampleNames = NULL,

     if(!keepFilt) {

         what0 <- what0[!grepl("^filt", what0)]

     }

-    if(verbose) cat("reading", files[1], "\n")

+    if(verbose) cat("[read.umtab2.chr] reading", files[1], "\n")

     scanPars <- list(sep = "", quote = "", quiet = TRUE, skip = 1,

                      what = what0, na.strings = c("NA", "?"))

     intab <- do.call(scan, c(scanPars, file = files[1]))

@@ -158,7 +158,7 @@ read.umtab2.chr2 <- function(files, sampleNames = NULL,

     what0[["ref"]] <- character(0)

     if(stranded)

         what0[["strand"]] <- character(0)

-    if(verbose) cat("parsing all samples first\n")

+    if(verbose) cat("[read.umtab2.chr] parsing all samples first\n")

     scanPars <- list(sep = "", quote = "", quiet = TRUE, skip = 1,

                      what = what0, na.strings = c("NA", "?"))

     ## First we get the locations of everything in all files

@@ -202,7 +202,7 @@ read.umtab2.chr2 <- function(files, sampleNames = NULL,

     stopifnot(all(c(keepM, keepU) %in% c(allMnames, allUnames)))

     for(ii in seq_along(files)) {

-        if(verbose) cat("reading", files[ii], "\n")

+        if(verbose) cat("[read.umtab2.chr] reading", files[ii], "\n")

         intab <- do.call(scan, c(scanPars, file = files[ii]))

         if(length(intab[[1]]) == 0)

             next

@@ -290,7 +290,7 @@ read.umtab.chr <- function(files, sampleNames = NULL,

         return(list())

     line1 <- strsplit(line1, "\t")[[1]]

     stopifnot(all(line1 == columnHeaders))

-    if(verbose) cat("reading", files[1], "\n")

+    if(verbose) cat("[read.umtab.chr] reading", files[1], "\n")

     scanPars <- list(sep = "", quote = "", quiet = TRUE, skip = 1,

                   what = what0, na.strings = c("NA", "?"))

     intab <- do.call(scan, c(scanPars, file = files[1]))

@@ -318,7 +318,7 @@ read.umtab.chr <- function(files, sampleNames = NULL,

     Ucy[,1] <- intab[["Ucy"]]

     csums[,1] <- as.integer(sapply(intab[c(allMnames, allUnames)], sum))

     for(ii in seq_along(files[-1]) + 1) {

-        if(verbose) cat("reading", files[ii], "\n")

+        if(verbose) cat("[read.umtab2] reading", files[ii], "\n")

         intab <- do.call(scan, c(scanPars, file = files[ii]))

         if(length(intab[[1]]) == 0)

             next

@@ -373,7 +373,7 @@ read.bsmoothDirRaw <- function(dir, seqnames = NULL, keepCycle = FALSE, keepFilt

         what0[grep("^filt", names(what0))] <- replicate(length(grep("^filt", names(what0))), NULL)

     outList <- lapply(allChrFiles, function(thisfile) {

         if(verbose)

-            cat(sprintf("Reading '%s'\n", thisfile))

+            cat(sprintf("[read.bsmoothDirRaw] Reading '%s'\n", thisfile))

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

             con <- gzfile(thisfile)

         else

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

     idxes <- seq_along(dirs)

     names(idxes) <- sampleNames

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

-        if(verbose) cat(sprintf("Reading dir '%s' ... ", dirs[ii]))

+        if(verbose) cat(sprintf("[read.bsmooth] Reading dir '%s' ... ", dirs[ii]))

         stime <- system.time({

             if(returnRaw) {

                 out <- read.bsmoothDirRaw(dir = dirs[ii], seqnames = seqnames, keepCycle = TRUE,

@@ -447,15 +447,15 @@ read.bsmooth <- function(dirs, sampleNames = NULL, seqnames = NULL, returnRaw =

                 out <- sampleRawToBSseq(raw, qualityCutoff = qualityCutoff, rmZeroCov, sampleName = sampleNames[ii])

             }

         })[3]

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

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

         out

     })

     if(!returnRaw) {

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

+        if(verbose) cat(sprintf("[read.bsmooth] Joining samples ... "))

         stime <- system.time({

             allOut <- combineList(allOut)

         })[3]

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

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

     }

     allOut

 }

@@ -464,9 +464,9 @@ parsingPipeline <- function(dirs, qualityCutoff = 20, outDir, seqnames = NULL,

                             subdir = "ev_bt2_cpg_tab", timing = FALSE) {

     if(!all(file.exists(dirs)))

         stop("not all directories in 'dirs' exists.")

-    cat("Parsing all files.\n")

+    cat("[parsingPipeline] Parsing all files.\n")

     oneDir <- function(dir) {

-        cat("  dir", basename(dir), ": ")

+        cat("[parsingPipeline]  dir", basename(dir), ": ")

         base <- basename(dir)

         cat("parsing, ")

         stime <- system.time({

@@ -476,7 +476,7 @@ parsingPipeline <- function(dirs, qualityCutoff = 20, outDir, seqnames = NULL,

             assign(paste0(base, ".raw"), raw)

         })[3]

         if(timing) {

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

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

             print(gc())

         }

         cat("saving, ")

@@ -489,7 +489,7 @@ parsingPipeline <- function(dirs, qualityCutoff = 20, outDir, seqnames = NULL,

             seqlevels(bsseq)[seqlevels(bsseq) == "chrgi|9626243|ref|NC_001416.1|"] <- "chrLambda"

         })[3]

         if(timing) {

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

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

             print(gc())

         }

         cat("ordering, ")

@@ -10,14 +10,17 @@

     \item validObject(BSseq) has been extended to also check for

     sampleNames consistency.

     \item Fixed a bug related to validity checking.

+    \item Increased maxk from 10,000 to 50,000 in calls to locfit, to

+    allowing fitting the model on genomes with unusally long chromosomes

+    (Thanks to Brian Herb for reporting).

     \item The class representation for class 'BSseq' has changed

     completely. The new class is build on 'SummarizedExperiment' from

     GenomicRanges instead of 'hasGRanges'.  Use 'x <- updateObject(x)' to

     update serialized (saved) objects.

     \item Fixing a problem in orderBSseq related to chromosome names.

-    \item Increased maxk from 10,000 to 50,000 in calls to locfit, to

-    allowing fitting the model on genomes with unusally long chromosomes

-    (Thanks to Brian Herb for reporting).

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

+    in BSmooth.

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

   }

 }

@@ -20,8 +20,8 @@ test_BSseq <- function() {

     checkEquals(dim(BStest), c(3,3))

     checkEquals(nrow(BStest), 3)

     checkEquals(ncol(BStest), 3)

-    checkEquals(getCoverage(BStest, type = "M"), M)

-    checkEquals(getCoverage(BStest, type = "Cov"), M + 2)

+    checkEquals(getCoverage(BStest, type = "M"), unname(M))

+    checkEquals(getCoverage(BStest, type = "Cov"), unname(M + 2))

     checkEquals(sampleNames(BStest), colnames(M))

     BStest2 <- BSseq(pos = 3:1, chr = rep("chr1", 3), M = M[3:1,],

@@ -204,8 +204,8 @@ we will only keep CpGs where at least 2 cancer samples and at least 2 normal sam

 line breaks in Sweave)

 <<keepLoci>>=

 BS.cov <- getCoverage(BS.cancer.ex.fit)

-keepLoci.ex <- which(rowSums(BS.cov[, c("C1", "C2", "C3")] >= 2) >= 2 &

-                  rowSums(BS.cov[, c("N1", "N2", "N3")] >= 2) >= 2)

+keepLoci.ex <- which(rowSums(BS.cov[, BS.cancer.ex$Type == "cancer"] >= 2) >= 2 &

+                     rowSums(BS.cov[, BS.cancer.ex$Type == "normal"] >= 2) >= 2)

 length(keepLoci.ex)

 BS.cancer.ex.fit <- BS.cancer.ex.fit[keepLoci.ex,]

 @