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

 Package: bsseq

-Version: 0.9.3

+Version: 0.9.4

 Title: Analyze, manage and store bisulfite sequencing data

 Description: Tools for analyzing and visualizing bisulfite sequencing data

 Author: Kasper Daniel Hansen

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

         parallel, matrixStats

 Imports: scales, stats, graphics, Biobase, locfit

 Suggests: RUnit, bsseqData

-Collate: hasGRanges.R BSseq_class.R BSseq_utils.R utils.R read.bsmooth.R

-        read.bismark.R BSmooth.R dmrFinder.R gof_stats.R plotting.R fisher.R fixes.R

+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 

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

 License: Artistic-2.0

 LazyData: yes

 biocViews: DNAMethylation

@@ -7,7 +7,7 @@ importFrom(BiocGenerics, "anyDuplicated", "cbind", "colnames",

            "combine", "density", "intersect", "lapply", "ncol",

            "nrow", "order", "paste", "pmax", "pmin", "rbind",

            "Reduce", "rep.int", "rownames", "sapply", "setdiff",

-           "strand", "strand<-", "union", "unique")

+           "strand", "strand<-", "union", "unique", "updateObject")

 importFrom(stats, "approxfun", "fisher.test", "ppoints",

            "predict", "preplot", "qchisq",

@@ -25,9 +25,8 @@ importFrom(scales, "alpha")

 importClassesFrom(Biobase, "AnnotatedDataFrame")

 importMethodsFrom(Biobase, "annotatedDataFrameFrom",

                   "pData", "pData<-",

-                  "phenoData", "phenoData<-",

                   "sampleNames", "sampleNames<-")

-

+importFrom(Biobase, "validMsg")

 ##

 ## Exporting

@@ -50,10 +49,9 @@ exportMethods("[", "show",

               "granges",

               "dim", "nrow", "ncol",

               "sampleNames", "sampleNames<-",

-              "phenoData", "phenoData<-",

               "pData", "pData<-",

               "findOverlaps", "subsetByOverlaps",

-              "combine")

+              "combine", "updateObject")

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

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

@@ -122,7 +122,11 @@ BSmooth <- function(BSseq, ns = 70, h = 1000, maxGap = 10^8, parallelBy = c("sam

         rownames(se.coef) <- NULL

         colnames(se.coef) <- sampleNames(BSseq)

     }

-    

+

+    if(!is.null(coef))

+        assay(BSseq, "coef") <- coef

+    if(!is.null(se.coef))

+        assay(BSseq, "se.coef") <- se.coef

     mytrans <- function(x) {

         y <- x

         ix <- which(x < 0)

@@ -133,8 +137,6 @@ BSmooth <- function(BSseq, ns = 70, h = 1000, maxGap = 10^8, parallelBy = c("sam

     }

     environment(mytrans) <- baseenv()

     BSseq@trans <- mytrans

-    BSseq@coef <- coef

-    BSseq@se.coef <- se.coef

     parameters <- list(smoothText = sprintf("BSmooth (ns = %d, h = %d, maxGap = %d)", ns, h, maxGap),

                        ns = ns, h = h, maxGap = maxGap)

     BSseq@parameters <- parameters

new file mode 100644

@@ -0,0 +1,39 @@

+setClass("BSseqTstat", contains = "hasGRanges", 

+         representation(stats = "matrix",

+                        parameters = "list")

+         )

+setValidity("BSseqTstat", function(object) {

+    msg <- NULL

+    if(length(object@gr) != nrow(object@stats))

+        msg <- c(msg, "length of 'gr' is different from the number of rows of 'stats'")

+    if(is.null(msg)) TRUE else msg

+})

+

+setMethod("show", signature(object = "BSseqTstat"),

+          function(object) {

+              cat("An object of type 'BSseqTstat' with\n")

+              cat(" ", length(object), "methylation loci\n")

+              cat("based on smoothed data:\n")

+              cat(" ", object@parameters$smoothText, "\n")

+              cat("with parameters\n")

+              cat(" ", object@parameters$tstatText, "\n")

+          })

+

+setMethod("[", "BSseqTstat", function(x, i, ...) {

+    if(missing(i))

+        stop("need [i] for subsetting")

+    if(missing(i))

+        return(x)

+    x@gr <- x@gr[i]

+    x@stats <- x@stats[i,, drop = FALSE]

+    x

+})

+

+BSseqTstat <- function(gr = NULL, stats = NULL, parameters = NULL) {

+    out <- new("BSseqTstat")

+    out@gr <- gr

+    out@stats <- stats

+    out@parameters <- parameters

+    out

+}

+

@@ -1,69 +1,25 @@

-setClassUnion("matrixOrNULL", c("matrix", "NULL"))

-

-setClass("BSseq", contains = "hasGRanges", 

-         representation(M = "matrix",

-                        Cov = "matrix",

-                        coef = "matrixOrNULL",

-                        se.coef = "matrixOrNULL",

-                        trans = "function",

-                        parameters = "list",

-                        phenoData = "AnnotatedDataFrame"))

-

-

-setClass("BSseqTstat", contains = "hasGRanges", 

-         representation(stats = "matrix",

-                        parameters = "list")

-         )

-         

-

-##

-## Simple accessor and replacement functions

-##

+setClass("BSseq", contains = "SummarizedExperiment", 

+         representation(trans = "function",

+                        parameters = "list"))

 setValidity("BSseq", function(object) {

-    msg <- NULL

-    nCpGs <- length(object@gr)

-    nSamples <- ncol(object@M)

-    if(nCpGs != nrow(object@M) ||

-       nCpGs != nrow(object@Cov))

-        msg <- c(msg, "length of the 'gr' slot is not equal to the nrows the 'M' or 'Cov' slots")

-    if(nSamples != ncol(object@Cov))

-        msg <- c(msg, "ncols of the 'M' slot is not equal to the ncols of the 'Cov' slot")

-    if(nSamples != nrow(object@phenoData))

-        msg <- c(msg, "the 'phenoData' slot does not match the dimensions of the 'M' slot")

-    if(any(object@M < 0))

-        msg <- c(msg, "the 'M' slot has negative entries")

-    if(any(object@Cov < 0))

-        msg <- c(msg, "the 'Cov' slot has negative entries")

-    if(any(object@M > object@Cov))

-        msg <- c(msg, "the 'M' slot has at least one entry bigger than the 'Cov' slot")

-    if(!is.null(object@coef) && (nCpGs != nrow(object@coef) ||

-                                 nSamples != ncol(object@coef)))

-        msg <- c(msg, "nrows and/or ncols of the 'coef' slot does not have the same dimensions as the 'M' slot")

-    if(!is.null(object@se.coef) && (nCpGs != nrow(object@se.coef) ||

-                                    nSamples != ncol(object@se.coef)))

-        msg <- c(msg, "nrows and/or ncols of the 'se.coef' slot does not have the same dimensions as the 'M' slot")

-    if(!is.null(rownames(object@M)) ||

-       !is.null(rownames(object@Cov)) ||

-       !is.null(rownames(object@coef)) ||

-       !is.null(rownames(object@se.coef))) warning("unnecessary rownames in object")

-    ## FIXME: check samplenames

-    if(any(colnames(object@M) != rownames(pData(object))) ||

-       any(colnames(object@M) != rownames(object@Cov)))

-        msg <- c(msg, "sample names are messed up: colnames of the M slot has to equal colnames of the Cov slot which has to equal rownames of the phenoData slot")

-    if(is.null(msg)) TRUE else msg

-})

-

-

-

-setValidity("BSseqTstat", function(object) {

-    msg <- NULL

-    if(length(object@gr) != nrow(object@stats))

-        msg <- c(msg, "length of 'gr' is different from the number of rows of 'stats'")

-    if(is.null(msg)) TRUE else msg

+    msg <- validMsg(NULL, .checkAssayNames(object, c("Cov", "M")))

+    if(class(rowData(object)) != "GRanges")

+        msg <- validMsg(msg, sprintf("object of class '%s' needs to have a 'GRanges' in slot 'rowData'", class(object)))

+    if(any(assay(object, "M") < 0))

+        msg <- validMsg(msg, "the 'M' assay has negative entries")

+    if(any(assay(object, "Cov") < 0))

+        msg <- validMsg(msg, "the 'Cov' assay has negative entries")

+    if(any(assay(object, "M") > assay(object, "Cov")))

+        msg <- validMsg(msg, "the 'M' assay has at least one entry bigger than the 'Cov' assay")

+    if(!is.null(rownames(assay(object, "M"))) ||

+       !is.null(rownames(assay(object, "Cov"))) ||

+       ("coef" %in% names(assays(object)) && !is.null(rownames(assay(object, "coef")))) ||

+       ("se.coef" %in% names(assays(object)) && !is.null(rownames(assay(object, "se.coef")))))

+        warning("unnecessary rownames in object")

+    if (is.null(msg)) TRUE else msg

 })

-

 setMethod("show", signature(object = "BSseq"),

           function(object) {

               cat("An object of type 'BSseq' with\n")

@@ -77,251 +33,68 @@ setMethod("show", signature(object = "BSseq"),

               }

           })

-setMethod("show", signature(object = "BSseqTstat"),

-          function(object) {

-              cat("An object of type 'BSseqTstat' with\n")

-              cat(" ", length(object), "methylation loci\n")

-              cat("based on smoothed data:\n")

-              cat(" ", object@parameters$smoothText, "\n")

-              cat("with parameters\n")

-              cat(" ", object@parameters$tstatText, "\n")

-          })

-

-

-##

-## eSet stuff

-##

-

-setMethod("dim", "BSseq", function(x) {

-    dim(x@M)

-})

-setMethod("nrow", "BSseq", function(x) {

-    nrow(x@M)

-})

-setMethod("ncol", "BSseq", function(x) {

-    ncol(x@M)

+setMethod("pData", "BSseq", function(object) {

+    object@colData

 })

-setMethod("phenoData", "BSseq", function(object) {

-    object@phenoData

-})

-setReplaceMethod("phenoData", signature = signature(

-                              object = "BSseq",

-                              value = "AnnotatedDataFrame"),

+setReplaceMethod("pData",

+                 signature = signature(

+                      object = "BSseq",

+                      value = "data.frame"),

                  function(object, value) {

-                     object@phenoData <- value

+                     object@colData <- as(value, "DataFrame")

                      object

                  })

-setMethod("pData", "BSseq", function(object) {

-    pData(phenoData(object))

-})

-setReplaceMethod("pData", signature = signature(

-                              object = "BSseq",

-                              value = "data.frame"),

+setReplaceMethod("pData",

+                 signature = signature(

+                      object = "BSseq",

+                      value = "DataFrame"),

                  function(object, value) {

-                     pd <- object@phenoData

-                     pData(pd) <- value

-                     object@phenoData <- pd

+                     object@colData <- value

                      object

                  })

 setMethod("sampleNames", "BSseq", function(object) {

-    sampleNames(phenoData(object))

+    colnames(object)

 })

 setReplaceMethod("sampleNames",

-                 signature(object = "BSseq", value = "ANY"),

+                 signature = signature(

+                     object = "BSseq",

+                     value = "ANY"),

                  function(object, value) {

-                     sampleNames(phenoData(object)) <- value

-                     colnames(object@M) <- value

-                     colnames(object@Cov) <- value

-                     if(!is.null(object@coef))

-                         colnames(object@coef) <- value

-                     if(!is.null(object@se.coef))

-                         colnames(object@se.coef) <- value

+                     colnames(object) <- value

                      object

                  })

-hasBeenSmoothed <- function(BSseq) {

-    !is.null(BSseq@coef)

-}

-

-##

-## Subsetting and combining

-##

-

-setMethod("[", "BSseq", function(x, i, j, ...) {

-    if(missing(drop))

-        drop <- FALSE

-    if(missing(i) && missing(j))

-        stop("need [i,j] for subsetting")

-    if(!missing(j))

-        x@phenoData <- phenoData(x)[j,, ..., drop = drop]

-    else

-        x@phenoData <- phenoData(x)

-    if(missing(i)) {

-        x@M <- getBSseq(x, "M")[, j, drop = FALSE]

-        x@Cov <- getBSseq(x, "Cov")[, j, drop = FALSE]

-        x@coef <- getBSseq(x, "coef")[, j, drop = FALSE]

-        x@se.coef <- getBSseq(x, "se.coef")[, j, drop = FALSE]

-        x@gr <- granges(x)

-    } else {

-        x@M <- getBSseq(x, "M")[i, j, drop = FALSE]

-        x@Cov <- getBSseq(x, "Cov")[i, j, drop = FALSE]

-        x@coef <- getBSseq(x, "coef")[i, j, drop = FALSE]

-        x@se.coef <- getBSseq(x, "se.coef")[i, j, drop = FALSE]

-        x@gr <- granges(x)[i]

-    }

-    x

-})

-

-setMethod("[", "BSseqTstat", function(x, i, ...) {

-    if(missing(i))

-        stop("need [i] for subsetting")

-    if(missing(i))

-        return(x)

-    x@gr <- x@gr[i]

-    x@stats <- x@stats[i,, drop = FALSE]

-    x

-})

-

-setMethod("combine", signature(x = "BSseq", y = "BSseq"), function(x, y, ...) {

-    ## All of this assumes that we are place x and y "next" to each other,

-    ##  ie. we are not combining the same set of samples sequenced at different times

-    if (class(x) != class(y))

-        stop(paste("objects must be the same class, but are ",

-                   class(x), ", ", class(y), sep=""))

-    if(hasBeenSmoothed(x) && hasBeenSmoothed(y) && !all.equal(x@trans, y@trans))

-        stop("'x' and 'y' need to be smoothed on the same scale")

-    phenoData <- combine(phenoData(x), phenoData(y))

-    if(identical(granges(x), granges(y))) {

-        gr <- granges(x)

-        M <- cbind(getBSseq(x, "M"), getBSseq(y, "M"))

-        Cov <- cbind(getBSseq(x, "Cov"), getBSseq(y, "Cov"))

-        if(!hasBeenSmoothed(x) || !hasBeenSmoothed(y)) {

-            coef <- NULL

-            se.coef <- NULL

-            trans <- NULL

-        } else {

-            coef <- cbind(getBSseq(x, "coef"), getBSseq(y, "coef"))

-            se.coef <- cbind(getBSseq(x, "se.coef"), getBSseq(y, "se.coef"))

-            trans <- getBSseq(x, "trans")

-        }

-    } else {

-        gr <- reduce(c(granges(x), granges(y)), min.gapwidth = 0L)

-        mm.x <- as.matrix(findOverlaps(gr, granges(x)))

-        mm.y <- as.matrix(findOverlaps(gr, granges(y)))

-        sampleNames <- c(sampleNames(x), sampleNames(y))

-        ## FIXME: there is no check that the two sampleNames are disjoint.

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

-        colnames(M) <- colnames(Cov) <- sampleNames

-        M[mm.x[,1], 1:ncol(x)] <- getBSseq(x, "M")[mm.x[,2],]

-        M[mm.y[,1], ncol(x) + 1:ncol(y)] <- getBSseq(y, "M")[mm.y[,2],]

-        Cov[mm.x[,1], 1:ncol(x)] <- getBSseq(x, "Cov")[mm.x[,2],]

-        Cov[mm.y[,1], ncol(x) + 1:ncol(y)] <- getBSseq(y, "Cov")[mm.y[,2],]

-        if(!hasBeenSmoothed(x) || !hasBeenSmoothed(y)) {

-            coef <- NULL

-            se.coef <- NULL

-            trans <- NULL

-        } else {

-            trans <- x@trans

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

-            colnames(coef) <- sampleNames(phenoData)

-            if(hasBeenSmoothed(x))

-                coef[mm.x[,1], 1:ncol(x)] <- getBSseq(x, "coef")[mm.x[,2],]

-            if(hasBeenSmoothed(y))

-                coef[mm.y[,1], ncol(x) + 1:ncol(y)] <- getBSseq(y, "coef")[mm.y[,2],]

-            if(is.null(getBSseq(x, "se.coef")) && is.null(getBSseq(x, "se.coef")))

-                se.coef <- NULL

-            else {

-                se.coef <- matrix(0, nrow = length(gr), ncol = length(sampleNames))

-                colnames(se.coef) <- sampleNames(phenoData)

-                if(!is.null(getBSseq(x, "se.coef")))

-                    se.coef[mm.x[,1], 1:ncol(x)] <- getBSseq(x, "se.coef")[mm.x[,2],]

-                if(!is.null(getBSseq(y, "se.coef")))

-                    se.coef[mm.y[,1], ncol(x) + 1:ncol(y)] <- getBSseq(y, "se.coef")[mm.y[,2],]

-            }

-        }

-    }

-    BSseq(gr = gr, M = M, Cov = Cov, coef = coef, se.coef = se.coef,

-          phenoData = phenoData, trans = trans, rmZeroCov = FALSE)

+setMethod("length", "BSseq", function(x) {

+    length(granges(x))

 })

-combineList <- function(x, ...) {

-    if(class(x) == "BSseq")

-        x <- list(x, ...)

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

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

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

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

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

-    })

-    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

-    } else {

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

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

-    }

-    phenoData <- Reduce(combine, lapply(x, phenoData))

-    BSseq(gr = gr, M = M, Cov = Cov, coef = coef, se.coef = se.coef,

-          phenoData = phenoData, trans = trans, rmZeroCov = FALSE)

+hasBeenSmoothed <- function(BSseq) {

+    "coef" %in% names(assays(BSseq))

 }

 getBSseq <- function(BSseq, type = c("Cov", "M", "gr", "coef", "se.coef", "trans", "parameters")) {

     type <- match.arg(type)

-    switch(type,

-           "Cov" = {

-               return(BSseq@Cov)

-           },

-           "M" = {

-               return(BSseq@M)

-           },

-           "gr" = {

-               return(BSseq@gr)

-           },

-           "coef" = {

-               return(BSseq@coef)

-           },

-           "se.coef" = {

-               return(BSseq@se.coef)

-           },

-           "trans" = {

-               return(BSseq@trans)

-           },

-           "parameters" = {

-               return(BSseq@parameters)

-           })

+    if(type %in% c("M", "Cov"))

+        return(assay(BSseq, type))

+    if(type %in% c("coef", "se.coef") && type %in% names(assays(BSseq)))

+        return(assay(BSseq, type))

+    if(type %in% c("coef", "se.coef"))

+       return(NULL)

+    if(type == "trans")

+        return(BSseq@trans)

+    if(type == "parameters")

+        return(BSseq@parameters)

+    if(type == "gr")

+        return(BSseq@rowData)

+    

 }

-

-

 BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

-                  trans = NULL, parameters = NULL, phenoData = NULL,

+                  trans = NULL, parameters = NULL, pData = NULL,

                   gr = NULL, pos = NULL, chr = NULL, sampleNames = NULL,

                   rmZeroCov = FALSE) {

     if(is.null(gr)) {

@@ -350,8 +123,10 @@ BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

     if(!is.null(names(gr)))

         names(gr) <- NULL

     ## deal with sampleNames

-    if(is.null(sampleNames) && !is.null(phenoData) && !is.null(sampleNames(phenoData)))

-        sampleNames <- sampleNames(phenoData)

+    if(!is(pData, "DataFrame"))

+        pData <- as(pData, "DataFrame")

+    if(is.null(sampleNames) && !is.null(pData) && !is.null(rownames(pData)))

+        sampleNames <- rownames(pData)

     if(is.null(sampleNames) && !is.null(colnames(M)))

         sampleNames <- colnames(M)

     if(is.null(sampleNames) && !is.null(colnames(Cov)))

@@ -401,31 +176,33 @@ BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

         colnames(M) <- sampleNames

     if(is.null(colnames(Cov)) || any(sampleNames != colnames(Cov)))

         colnames(Cov) <- sampleNames

-    if(is.null(phenoData))

-        phenoData <- annotatedDataFrameFrom(M, byrow = FALSE)

-    BSseq <- new("BSseq", gr = gr, M = M, Cov = Cov, phenoData = phenoData)

     if(!is.null(coef)) {

         if(!is.matrix(coef) ||

-           nrow(coef) != nrow(BSseq) ||

-           ncol(coef) != ncol(BSseq))

+           nrow(coef) != nrow(M) ||

+           ncol(coef) != ncol(M))

             stop("'coef' does not have the right dimensions")

         if(is.null(colnames(coef)) || any(sampleNames != colnames(coef)))

             colnames(coef) <- sampleNames

         if(!is.null(rownames(coef)))

             rownames(coef) <- NULL

-        BSseq@coef <- coef

     }

     if(!is.null(se.coef)) {

         if(!is.matrix(se.coef) ||

-           nrow(se.coef) != nrow(BSseq) ||

-           ncol(se.coef) != ncol(BSseq))

+           nrow(se.coef) != nrow(M) ||

+           ncol(se.coef) != ncol(M))

             stop("'se.coef' does not have the right dimensions")

         if(is.null(colnames(se.coef)) || any(sampleNames != colnames(se.coef)))

             colnames(se.coef) <- sampleNames

         if(!is.null(rownames(se.coef)))

             rownames(se.coef) <- NULL

-        BSseq@se.coef <- se.coef

     }

+    assays <- SimpleList(M = M, Cov = Cov, coef = coef, se.coef = se.coef)

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

+    if(is.null(pData) || all(dim(pData) == c(0,0)))

+        BSseq <- SummarizedExperiment(assays = assays, rowData = gr)

+    else

+        BSseq <- SummarizedExperiment(assays = assays, rowData = gr, colData = pData)

+    BSseq <- as(BSseq, "BSseq")

     if(is.function(trans))

         BSseq@trans <- trans

     if(is.list(parameters))

@@ -434,13 +211,15 @@ BSseq <- function(M = NULL, Cov = NULL, coef = NULL, se.coef = NULL,

 }

-BSseqTstat <- function(gr = NULL, stats = NULL, parameters = NULL) {

-    out <- new("BSseqTstat")

-    out@gr <- gr

-    out@stats <- stats

-    out@parameters <- parameters

-    out

-}

+setMethod("updateObject", "BSseq",

+          function(object, ...) {

+               if(!is(try(object@assays, silent = TRUE), "try-error"))

+                   return(object)

+               BSseq(gr = object@gr, M = object@M, Cov = object@Cov,

+                     coef = object@coef, se.coef = object@se.coef,

+                     trans = object@trans, parameters = object@parameters,

+                     pData = object@phenoData@data)

+           })

@@ -10,14 +10,14 @@ collapseBSseq <- function(BSseq, columns) {

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

     }

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

-        rowSums(BSseq@M[, ss, drop = FALSE])

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

     }))

     colnames(M) <- names(sp)

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

-        rowSums(BSseq@Cov[, ss, drop = FALSE])

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

     }))

     colnames(Cov) <- names(sp)

-    BSseq(gr = BSseq@gr, M = M, Cov = Cov)

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

 }

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

new file mode 100644

@@ -0,0 +1,106 @@

+setMethod("combine", signature(x = "BSseq", y = "BSseq"), function(x, y, ...) {

+    ## All of this assumes that we are place x and y "next" to each other,

+    ##  ie. we are not combining the same set of samples sequenced at different times

+    if (class(x) != class(y))

+        stop(paste("objects must be the same class, but are ",

+                   class(x), ", ", class(y), sep=""))

+    if(hasBeenSmoothed(x) && hasBeenSmoothed(y) && !all.equal(x@trans, y@trans))

+        stop("'x' and 'y' need to be smoothed on the same scale")

+    pData <- combine(as(pData(x), "data.frame"), as(pData(y), "data.frame"))

+    if(identical(granges(x), granges(y))) {

+        gr <- granges(x)

+        M <- cbind(getBSseq(x, "M"), getBSseq(y, "M"))

+        Cov <- cbind(getBSseq(x, "Cov"), getBSseq(y, "Cov"))

+        if(!hasBeenSmoothed(x) || !hasBeenSmoothed(y)) {

+            coef <- NULL

+            se.coef <- NULL

+            trans <- NULL

+        } else {

+            coef <- cbind(getBSseq(x, "coef"), getBSseq(y, "coef"))

+            se.coef <- cbind(getBSseq(x, "se.coef"), getBSseq(y, "se.coef"))

+            trans <- getBSseq(x, "trans")

+        }

+    } else {

+        gr <- reduce(c(granges(x), granges(y)), min.gapwidth = 0L)

+        mm.x <- as.matrix(findOverlaps(gr, granges(x)))

+        mm.y <- as.matrix(findOverlaps(gr, granges(y)))

+        sampleNames <- c(sampleNames(x), sampleNames(y))

+        ## FIXME: there is no check that the two sampleNames are disjoint.

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

+        colnames(M) <- colnames(Cov) <- sampleNames

+        M[mm.x[,1], 1:ncol(x)] <- getBSseq(x, "M")[mm.x[,2],]

+        M[mm.y[,1], ncol(x) + 1:ncol(y)] <- getBSseq(y, "M")[mm.y[,2],]

+        Cov[mm.x[,1], 1:ncol(x)] <- getBSseq(x, "Cov")[mm.x[,2],]

+        Cov[mm.y[,1], ncol(x) + 1:ncol(y)] <- getBSseq(y, "Cov")[mm.y[,2],]

+        if(!hasBeenSmoothed(x) || !hasBeenSmoothed(y)) {

+            coef <- NULL

+            se.coef <- NULL

+            trans <- NULL

+        } else {

+            trans <- x@trans

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

+            colnames(coef) <- rownames(pData)

+            if(hasBeenSmoothed(x))

+                coef[mm.x[,1], 1:ncol(x)] <- getBSseq(x, "coef")[mm.x[,2],]

+            if(hasBeenSmoothed(y))

+                coef[mm.y[,1], ncol(x) + 1:ncol(y)] <- getBSseq(y, "coef")[mm.y[,2],]

+            if(is.null(getBSseq(x, "se.coef")) && is.null(getBSseq(x, "se.coef")))

+                se.coef <- NULL

+            else {

+                se.coef <- matrix(0, nrow = length(gr), ncol = length(sampleNames))

+                colnames(se.coef) <- sampleNames(pData)

+                if(!is.null(getBSseq(x, "se.coef")))

+                    se.coef[mm.x[,1], 1:ncol(x)] <- getBSseq(x, "se.coef")[mm.x[,2],]

+                if(!is.null(getBSseq(y, "se.coef")))

+                    se.coef[mm.y[,1], ncol(x) + 1:ncol(y)] <- getBSseq(y, "se.coef")[mm.y[,2],]

+            }

+        }

+    }

+    BSseq(gr = gr, M = M, Cov = Cov, coef = coef, se.coef = se.coef,

+          pData = pData, trans = trans, rmZeroCov = FALSE)

+})

+

+combineList <- function(x, ...) {

+    if(class(x) == "BSseq")

+        x <- list(x, ...)

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

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

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

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

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

+    })

+    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

+    } else {

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

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

+    }

+    pData <- Reduce(combine, lapply(x, pData))

+    BSseq(gr = gr, M = M, Cov = Cov, coef = coef, se.coef = se.coef,

+          pData = pData, trans = trans, rmZeroCov = FALSE)

+}

+

@@ -26,3 +26,11 @@ data.frame2GRanges <- function(df, keepColumns = FALSE, ignoreStrand = FALSE) {

     gr

 }

+.checkAssayNames <- function(object, names) {

+    nms <- names(assays(object, withDimnames = FALSE))

+    if(!all(names %in% nms))

+        return(sprintf("object of class '%s' needs to have assay slots with names '%s'",

+                       class(object), paste0(names, collapse = ", ")))

+    else

+        NULL

+}

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

@@ -10,6 +10,10 @@

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

     sampleNames consistency.

     \item Fixed a bug related to validity checking.

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

   }

 }

@@ -13,7 +13,7 @@

 }

 \usage{data(BS.chr22)}

 \format{

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

+  An object of class \code{BSseq}.

 }

 \details{

   All coordinates are in hg18.

@@ -12,7 +12,7 @@ BSmooth(BSseq, ns = 70, h = 1000, maxGap = 10^8,

   mc.cores = 1, keep.se = FALSE, verbose = TRUE)

 }

 \arguments{

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

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

   \item{ns}{The minimum number of methylation loci in a smoothing window.}

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

   \item{maxGap}{The maximum gap between two methylation loci, before the

@@ -43,7 +43,7 @@ BSmooth(BSseq, ns = 70, h = 1000, maxGap = 10^8,

 }

 \value{

-  An object of class \code{"BSseq"}, containing smoothed values and

+  An object of class \code{BSseq}, containing smoothed values and

   optionally standard errors for these.

 }

 \author{

@@ -59,6 +59,7 @@ BSmooth(BSseq, ns = 70, h = 1000, maxGap = 10^8,

 }

 \examples{

 \dontrun{

+data(BS.chr22)

 BS.fit <- BSmooth(BS.chr22, verbose = TRUE)

 BS.fit

 }

@@ -14,7 +14,7 @@ BSmooth.tstat(BSseq, group1, group2,

   maxGap = NULL, qSd = 0.75, k = 101, mc.cores = 1, verbose = TRUE)

 }

 \arguments{

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

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

   \item{group1}{A vector of sample names or indexes for the \sQuote{treatment}

     group.} 

   \item{group2}{A vector of sample names or indexes for the \sQuote{control}

@@ -48,7 +48,7 @@ BSmooth.tstat(BSseq, group1, group2,

   Additional details in the reference.

 }

 \value{

-  An object of class \code{"BSseqTstat"}.

+  An object of class \code{BSseqTstat}.

 }

 \author{

   Kasper Daniel Hansen \email{khansen@jhsph.edu}

@@ -70,6 +70,7 @@ BSmooth.tstat(BSseq, group1, group2,

 if(require(bsseqData)) {

   data(keepLoci.ex)

   data(BS.cancer.ex.fit)

+  BS.cancer.ex.fit <- updateObject(BS.cancer.ex.fit)

   ## Remember to subset the BSseq object, see vignette for explanation

   BS.tstat <- BSmooth.tstat(BS.cancer.ex.fit[keepLoci.ex,],

                             group1 = c("C1", "C2", "C3"),

@@ -5,15 +5,13 @@

 \alias{[,BSseq-method}

 \alias{combine,BSseq,BSseq-method}

 \alias{combineList}

-\alias{dim,BSseq-method}

-\alias{ncol,BSseq-method}

-\alias{nrow,BSseq-method}

+\alias{length,BSseq-method}

 \alias{pData,BSseq-method}

-\alias{phenoData,BSseq-method}

-\alias{phenoData<-,BSseq,AnnotatedDataFrame-method}

 \alias{pData<-,BSseq,data.frame-method}

+\alias{pData<-,BSseq,DataFrame-method}

 \alias{sampleNames,BSseq-method}

 \alias{sampleNames<-,BSseq,ANY-method}

+\alias{updateObject,BSseq-method}

 \alias{show,BSseq-method}

 \alias{getBSseq}

 \alias{collapseBSseq}

@@ -31,30 +29,16 @@ An object from the class links together several pieces of information.

 (1) genomic locations stored as a \code{GRanges} object, a location by

 samples matrix of M values, a location by samples matrix of Cov

 (coverage) values and phenodata information.  In addition, there are

-slots for representing smoothed data.  Objects can be created by calls

-of the form \code{BSseq(...)}.   

+slots for representing smoothed data. This class is an extension of

+\code{SummarizedExperiment}.

 }

 \section{Slots}{

   \describe{

-    \item{\code{gr}:}{Object of class \code{"GRanges"} giving genomic

-      locations.}

-    \item{\code{M}:}{Object of class \code{"matrix"}.  This is a

-      location by sample matrix of the number of reads supporting

-      methylation.} 

-    \item{\code{Cov}:}{Object of class \code{"matrix"}.  This is a

-      location by sample matrix of the coverage.}

-    \item{\code{coef}:}{Object of class \code{"matrixOrNULL"}.  This is

-      an optional slot representing smoothed data.}

-    \item{\code{se.coef}:}{Object of class \code{"matrixOrNULL"}.  This

-      is an optional slot representing standard errors of the

-      smoothing.} 

-    \item{\code{trans}:}{Object of class \code{"function"}.  This

+    \item{\code{trans}:}{Object of class \code{function}.  This

       function transforms the \code{coef} slot from the scale the

       smoothing was done to the 0-1 methylation scale.}

-    \item{\code{parameters}:}{Object of class \code{"list"}.  A list of

+    \item{\code{parameters}:}{Object of class \code{list}.  A list of

       parameters representing for example how the data was smoothed.}

-    \item{\code{phenoData}:}{Object of class

-      \code{"AnnotatedDataFrame"}.  Sample information.}

   }

 }

 \section{Methods}{

@@ -62,17 +46,13 @@ of the form \code{BSseq(...)}.

     \item{[}{\code{signature(x = "BSseq")}: Subsetting by location

       (using integer indices) or sample (using integers or sample

       names).}

-    \item{dim}{The dimensions of the object (number of locations by

-      number of samples).}

-    \item{ncol}{The number of columns (equal to the number of samples).} 

-    \item{nrow}{The number of rows (equal to the number of genomic

-      locations).} 

+    \item{length}{Unlike for \code{SummarizedExperiment},

+      \code{length()} is the number of methylation loci (equal to

+      \code{length(granges(x))}).} 

     \item{sampleNames,sampleNames<-}{Sample names and its replacement

-      function for the object.}

-    \item{phenoData,phenoData<-}{Obtain and replace the \code{phenoData}

-      slot.}

+      function for the object.  This is an alias for \code{colnames}.}

     \item{pData,pData<-}{Obtain and replace the \code{pData} slot of the

-      \code{phenoData} slot.}

+      \code{phenoData} slot. This is an alias for \code{colData}.}

     \item{show}{The show method.}

     \item{combine}{This function combines two \code{BSSeq} objects.

       The genomic locations of the new object is the union of the

@@ -88,13 +68,13 @@ of the form \code{BSseq(...)}.

   \code{subsetByOverlaps}. 

   There are a number of almost methods-like functions for operating on

-  objects of class \code{"BSseq"}, including \code{getBSseq},

+  objects of class \code{BSseq}, including \code{getBSseq},

   \code{collapseBSseq}, and \code{orderBSseq}.  They are detailed below.

   \describe{