1. bsseq
  2. R
  3. BSseqTstat_class.R
Browse code

Revert to ordinary matrix for BSseqStat and BSseqTstat

Peter Hickey authored on 30/09/2018 01:08:13
Showing 1 changed files
R/BSseqTstat_class.R
History View file @ 7563138
... ... 
@@ -1,5 +1,5 @@
1 1 
 setClass("BSseqTstat", contains = "hasGRanges",
2 
-         representation(stats = "DelayedMatrix",
2 
+         representation(stats = "matrix",
3 3 
                         parameters = "list")
4 4 
          )
5 5 
 setValidity("BSseqTstat", function(object) {
... ... 
@@ -17,11 +17,6 @@ setMethod("show", signature(object = "BSseqTstat"),
17 17 
               cat(" ", object@parameters$smoothText, "\n")
18 18 
               cat("with parameters\n")
19 19 
               cat(" ", object@parameters$tstatText, "\n")
20 
-              if (.isHDF5ArrayBacked(object@stats)) {
21 
-                  cat("'stats' slot is HDF5Array-backed\n")
22 
-              } else {
23 
-                  cat("'stats' slot is in-memory\n")
24 
-              }
25 20 
           })
26 21
27 22 
 setMethod("[", "BSseqTstat", function(x, i, ...) {
... ... 
@@ -37,7 +32,7 @@ setMethod("[", "BSseqTstat", function(x, i, ...) {
37 32 
 BSseqTstat <- function(gr = NULL, stats = NULL, parameters = NULL) {
38 33 
     out <- new("BSseqTstat")
39 34 
     out@gr <- gr
40 
-    out@stats <- .DelayedMatrix(stats)
35 
+    out@stats <- stats
41 36 
     out@parameters <- parameters
42 37 
     out
43 38 
 }
... ... 
@@ -75,7 +70,7 @@ plot.BSseqTstat <- function(x, y, ...) {
75 70 
 setMethod("updateObject", "BSseqTstat",
76 71 
           function(object, ...) {
77 72 
               stats <- object@stats
78 
-              stats <- .DelayedMatrix(stats)
73 
+              stats <- stats
79 74 
               object@stats <- stats
80 75 
               object
81 76 
           }
Browse code

bsseq now uses DelayedMatrix objects from the DelayedArray package for all matrix-like data

Peter Hickey authored on 07/04/2017 17:42:29
Showing 1 changed files
R/BSseqTstat_class.R
History View file @ 95e9059
... ... 
@@ -1,5 +1,5 @@
1 
-setClass("BSseqTstat", contains = "hasGRanges",
2 
-         representation(stats = "matrix",
1 
+setClass("BSseqTstat", contains = "hasGRanges",
2 
+         representation(stats = "DelayedMatrix",
3 3 
                         parameters = "list")
4 4 
          )
5 5 
 setValidity("BSseqTstat", function(object) {
... ... 
@@ -17,6 +17,11 @@ setMethod("show", signature(object = "BSseqTstat"),
17 17 
               cat(" ", object@parameters$smoothText, "\n")
18 18 
               cat("with parameters\n")
19 19 
               cat(" ", object@parameters$tstatText, "\n")
20 
+              if (.isHDF5ArrayBacked(object@stats)) {
21 
+                  cat("'stats' slot is HDF5Array-backed\n")
22 
+              } else {
23 
+                  cat("'stats' slot is in-memory\n")
24 
+              }
20 25 
           })
21 26
22 27 
 setMethod("[", "BSseqTstat", function(x, i, ...) {
... ... 
@@ -32,14 +37,14 @@ setMethod("[", "BSseqTstat", function(x, i, ...) {
32 37 
 BSseqTstat <- function(gr = NULL, stats = NULL, parameters = NULL) {
33 38 
     out <- new("BSseqTstat")
34 39 
     out@gr <- gr
35 
-    out@stats <- stats
40 
+    out@stats <- .DelayedMatrix(stats)
36 41 
     out@parameters <- parameters
37 42 
     out
38 43 
 }
39 44
40 45 
 summary.BSseqTstat <- function(object, ...) {
41 
-    quant <- quantile(getStats(object)[, "tstat.corrected"],
42 
-                      prob = c(0.0001, 0.001, 0.01, 0.5, 0.99, 0.999, 0.9999))
46 
+    quant <- .quantile(getStats(object)[, "tstat.corrected"],
47 
+                       prob = c(0.0001, 0.001, 0.01, 0.5, 0.99, 0.999, 0.9999))
43 48 
     quant <- t(t(quant))
44 49 
     colnames(quant) <- "quantiles"
45 50 
     out <- list(quantiles = quant)
... ... 
@@ -53,9 +58,11 @@ print.summary.BSseqTstat <- function(x, ...) {
53 58
54 59 
 plot.BSseqTstat <- function(x, y, ...) {
55 60 
     tstat <- getStats(x)[, "tstat"]
61 
+    tstat <- as.array(tstat)
56 62 
     plot(density(tstat), xlim = c(-10,10), col = "blue", main = "")
57 63 
     if("tstat.corrected" %in% colnames(getStats(x))) {
58 64 
         tstat.cor <- getStats(x)[, "tstat.corrected"]
65 
+        tstat.cor <- as.array(tstat.cor)
59 66 
         lines(density(tstat.cor), col = "black")
60 67 
         legend("topleft", legend = c("uncorrected", "corrected"), lty = c(1,1),
61 68 
                col = c("blue", "black"))
... ... 
@@ -65,3 +72,11 @@ plot.BSseqTstat <- function(x, y, ...) {
65 72 
     }
66 73 
 }
67 74
75 
+setMethod("updateObject", "BSseqTstat",
76 
+          function(object, ...) {
77 
+              stats <- object@stats
78 
+              stats <- .DelayedMatrix(stats)
79 
+              object@stats <- stats
80 
+              object
81 
+          }
82 
+)
Browse code

fixed conflict

From: Kasper Daniel Hansen <kasperdanielhansen@gmail.com>

git-svn-id: https://hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/bsseq@110074 bc3139a8-67e5-0310-9ffc-ced21a209358

khansen authored on 29/10/2015 19:12:43
Showing 1 changed files
R/BSseqTstat_class.R
History View file @ 65141d1
... ... 
@@ -65,42 +65,3 @@ plot.BSseqTstat <- function(x, y, ...) {
65 65 
     }
66 66 
 }
67 67
68 
-getStats <- function(BSseqTstat, regions = NULL, stat = "tstat.corrected") {
69 
-    stopifnot(is(BSseqTstat, "BSseqTstat"))
70 
-    if(is.null(regions))
71 
-        return(BSseqTstat@stats)
72 
-    if(class(regions) == "data.frame")
73 
-        regions <- data.frame2GRanges(regions)
74 
-    stopifnot(stat %in% colnames(BSseqTstat@stats))
75 
-    stopifnot(length(stat) == 1)
76 
-    stopifnot(is(regions, "GenomicRanges"))
77 
-    ov <- findOverlaps(BSseqTstat, regions)
78 
-    ov.sp <- split(queryHits(ov), subjectHits(ov))
79 
-    getRegionStats <- function(idx) {
80 
-        mat <- BSseqTstat@stats[idx,, drop=FALSE]
81 
-        areaStat <- sum(mat[, stat])
82 
-        maxStat <- max(mat[, stat])
83 
-        c(areaStat, maxStat)
84 
-    }
85 
-    stats <- matrix(NA, ncol = 2, nrow = length(regions))
86 
-    colnames(stats) <- c("areaStat", "maxStat")
87 
-    tmp <- lapply(ov.sp, getRegionStats)
88 
-    stats[as.integer(names(tmp)),] <- do.call(rbind, tmp)
89 
-    out <- as.data.frame(stats)
90 
-    if(! stat %in% c("tstat.corrected", "tstat"))
91 
-        return(out)
92 
-    getRegionStats_ttest <- function(idx) {
93 
-        mat <- BSseqTstat@stats[idx,, drop=FALSE]
94 
-        group1.mean <- mean(mat[, "group1.means"])
95 
-        group2.mean <- mean(mat[, "group2.means"])
96 
-        meanDiff <- mean(mat[, "group1.means"] - mat[, "group2.means"])
97 
-        tstat.sd <- mean(mat[, "tstat.sd"])
98 
-        c(meanDiff, group1.mean, group2.mean, tstat.sd)
99 
-    }
100 
-    stats_ttest<- matrix(NA, ncol = 4, nrow = length(regions))
101 
-    colnames(stats_ttest) <- c("meanDiff", "group1.mean", "group2.mean", "tstat.sd")
102 
-    tmp <- lapply(ov.sp, getRegionStats_ttest)
103 
-    stats_ttest[as.integer(names(tmp)),] <- do.call(rbind, tmp)
104 
-    out <- cbind(out, as.data.frame(stats_ttest))
105 
-    out
106 
-}
Browse code

better argument cheking in various functions

git-svn-id: https://hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/bsseq@80445 bc3139a8-67e5-0310-9ffc-ced21a209358

khansen authored on 14/09/2013 15:06:17
Showing 1 changed files
R/BSseqTstat_class.R
History View file @ c5bee10
... ... 
@@ -65,20 +65,21 @@ plot.BSseqTstat <- function(x, y, ...) {
65 65 
     }
66 66 
 }
67 67
68 
-getStats <- function(BSseqTstat, regions = NULL, column = "tstat.corrected") {
68 
+getStats <- function(BSseqTstat, regions = NULL, stat = "tstat.corrected") {
69 69 
     stopifnot(is(BSseqTstat, "BSseqTstat"))
70 
-    stopifnot(column %in% colnames(BSseqTstat@stats))
71 
-    if(class(regions) == "data.frame")
72 
-        regions <- data.frame2GRanges(regions)
73 70 
     if(is.null(regions))
74 71 
         return(BSseqTstat@stats)
72 
+    if(class(regions) == "data.frame")
73 
+        regions <- data.frame2GRanges(regions)
74 
+    stopifnot(stat %in% colnames(BSseqTstat@stats))
75 
+    stopifnot(length(stat) == 1)
75 76 
     stopifnot(is(regions, "GenomicRanges"))
76 77 
     ov <- findOverlaps(BSseqTstat, regions)
77 78 
     ov.sp <- split(queryHits(ov), subjectHits(ov))
78 79 
     getRegionStats <- function(idx) {
79 80 
         mat <- BSseqTstat@stats[idx,, drop=FALSE]
80 
-        areaStat <- sum(mat[, column])
81 
-        maxStat <- max(mat[, column])
81 
+        areaStat <- sum(mat[, stat])
82 
+        maxStat <- max(mat[, stat])
82 83 
         c(areaStat, maxStat)
83 84 
     }
84 85 
     stats <- matrix(NA, ncol = 2, nrow = length(regions))
... ... 
@@ -86,7 +87,7 @@ getStats <- function(BSseqTstat, regions = NULL, column = "tstat.corrected") {
86 87 
     tmp <- lapply(ov.sp, getRegionStats)
87 88 
     stats[as.integer(names(tmp)),] <- do.call(rbind, tmp)
88 89 
     out <- as.data.frame(stats)
89 
-    if(! column %in% c("tstat.corrected", "tstat"))
90 
+    if(! stat %in% c("tstat.corrected", "tstat"))
90 91 
         return(out)
91 92 
     getRegionStats_ttest <- function(idx) {
92 93 
         mat <- BSseqTstat@stats[idx,, drop=FALSE]
Browse code

important bugfix to getStats (affecting dmrFinder)

git-svn-id: https://hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/bsseq@78749 bc3139a8-67e5-0310-9ffc-ced21a209358

khansen authored on 22/07/2013 19:43:14
Showing 1 changed files
R/BSseqTstat_class.R
History View file @ f4870c2
... ... 
@@ -98,7 +98,7 @@ getStats <- function(BSseqTstat, regions = NULL, column = "tstat.corrected") {
98 98 
     }
99 99 
     stats_ttest<- matrix(NA, ncol = 4, nrow = length(regions))
100 100 
     colnames(stats_ttest) <- c("meanDiff", "group1.mean", "group2.mean", "tstat.sd")
101 
-    tmp <- lapply(ov.sp, getRegionStats)
101 
+    tmp <- lapply(ov.sp, getRegionStats_ttest)
102 102 
     stats_ttest[as.integer(names(tmp)),] <- do.call(rbind, tmp)
103 103 
     out <- cbind(out, as.data.frame(stats_ttest))
104 104 
     out
Browse code

various bugfixes

git-svn-id: https://hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/bsseq@78387 bc3139a8-67e5-0310-9ffc-ced21a209358

khansen authored on 12/07/2013 20:53:31
Showing 1 changed files
R/BSseqTstat_class.R
History View file @ 1e3d27d
... ... 
@@ -37,3 +37,69 @@ BSseqTstat <- function(gr = NULL, stats = NULL, parameters = NULL) {
37 37 
     out
38 38 
 }
39 39
40 
+summary.BSseqTstat <- function(object, ...) {
41 
+    quant <- quantile(getStats(object)[, "tstat.corrected"],
42 
+                      prob = c(0.0001, 0.001, 0.01, 0.5, 0.99, 0.999, 0.9999))
43 
+    quant <- t(t(quant))
44 
+    colnames(quant) <- "quantiles"
45 
+    out <- list(quantiles = quant)
46 
+    class(out) <- "summary.BSseqTstat"
47 
+    out
48 
+}
49 
+
50 
+print.summary.BSseqTstat <- function(x, ...) {
51 
+    print(as.matrix(x$quantiles))
52 
+}
53 
+
54 
+plot.BSseqTstat <- function(x, y, ...) {
55 
+    tstat <- getStats(x)[, "tstat"]
56 
+    plot(density(tstat), xlim = c(-10,10), col = "blue", main = "")
57 
+    if("tstat.corrected" %in% colnames(getStats(x))) {
58 
+        tstat.cor <- getStats(x)[, "tstat.corrected"]
59 
+        lines(density(tstat.cor), col = "black")
60 
+        legend("topleft", legend = c("uncorrected", "corrected"), lty = c(1,1),
61 
+               col = c("blue", "black"))
62 
+    } else {
63 
+        legend("topleft", legend = c("uncorrected"), lty = 1,
64 
+               col = c("blue"))
65 
+    }
66 
+}
67 
+
68 
+getStats <- function(BSseqTstat, regions = NULL, column = "tstat.corrected") {
69 
+    stopifnot(is(BSseqTstat, "BSseqTstat"))
70 
+    stopifnot(column %in% colnames(BSseqTstat@stats))
71 
+    if(class(regions) == "data.frame")
72 
+        regions <- data.frame2GRanges(regions)
73 
+    if(is.null(regions))
74 
+        return(BSseqTstat@stats)
75 
+    stopifnot(is(regions, "GenomicRanges"))
76 
+    ov <- findOverlaps(BSseqTstat, regions)
77 
+    ov.sp <- split(queryHits(ov), subjectHits(ov))
78 
+    getRegionStats <- function(idx) {
79 
+        mat <- BSseqTstat@stats[idx,, drop=FALSE]
80 
+        areaStat <- sum(mat[, column])
81 
+        maxStat <- max(mat[, column])
82 
+        c(areaStat, maxStat)
83 
+    }
84 
+    stats <- matrix(NA, ncol = 2, nrow = length(regions))
85 
+    colnames(stats) <- c("areaStat", "maxStat")
86 
+    tmp <- lapply(ov.sp, getRegionStats)
87 
+    stats[as.integer(names(tmp)),] <- do.call(rbind, tmp)
88 
+    out <- as.data.frame(stats)
89 
+    if(! column %in% c("tstat.corrected", "tstat"))
90 
+        return(out)
91 
+    getRegionStats_ttest <- function(idx) {
92 
+        mat <- BSseqTstat@stats[idx,, drop=FALSE]
93 
+        group1.mean <- mean(mat[, "group1.means"])
94 
+        group2.mean <- mean(mat[, "group2.means"])
95 
+        meanDiff <- mean(mat[, "group1.means"] - mat[, "group2.means"])
96 
+        tstat.sd <- mean(mat[, "tstat.sd"])
97 
+        c(meanDiff, group1.mean, group2.mean, tstat.sd)
98 
+    }
99 
+    stats_ttest<- matrix(NA, ncol = 4, nrow = length(regions))
100 
+    colnames(stats_ttest) <- c("meanDiff", "group1.mean", "group2.mean", "tstat.sd")
101 
+    tmp <- lapply(ov.sp, getRegionStats)
102 
+    stats_ttest[as.integer(names(tmp)),] <- do.call(rbind, tmp)
103 
+    out <- cbind(out, as.data.frame(stats_ttest))
104 
+    out
105 
+}
Browse code

moved to a SummarizedExperiment backend for class 'BSseq'

git-svn-id: https://hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/bsseq@78143 bc3139a8-67e5-0310-9ffc-ced21a209358

khansen authored on 06/07/2013 21:26:40
Showing 1 changed files
R/BSseqTstat_class.R
History View file @ 8e0ebd6
1 1 
new file mode 100644
... ... 
@@ -0,0 +1,39 @@
1 
+setClass("BSseqTstat", contains = "hasGRanges",
2 
+         representation(stats = "matrix",
3 
+                        parameters = "list")
4 
+         )
5 
+setValidity("BSseqTstat", function(object) {
6 
+    msg <- NULL
7 
+    if(length(object@gr) != nrow(object@stats))
8 
+        msg <- c(msg, "length of 'gr' is different from the number of rows of 'stats'")
9 
+    if(is.null(msg)) TRUE else msg
10 
+})
11 
+
12 
+setMethod("show", signature(object = "BSseqTstat"),
13 
+          function(object) {
14 
+              cat("An object of type 'BSseqTstat' with\n")
15 
+              cat(" ", length(object), "methylation loci\n")
16 
+              cat("based on smoothed data:\n")
17 
+              cat(" ", object@parameters$smoothText, "\n")
18 
+              cat("with parameters\n")
19 
+              cat(" ", object@parameters$tstatText, "\n")
20 
+          })
21 
+
22 
+setMethod("[", "BSseqTstat", function(x, i, ...) {
23 
+    if(missing(i))
24 
+        stop("need [i] for subsetting")
25 
+    if(missing(i))
26 
+        return(x)
27 
+    x@gr <- x@gr[i]
28 
+    x@stats <- x@stats[i,, drop = FALSE]
29 
+    x
30 
+})
31 
+
32 
+BSseqTstat <- function(gr = NULL, stats = NULL, parameters = NULL) {
33 
+    out <- new("BSseqTstat")
34 
+    out@gr <- gr
35 
+    out@stats <- stats
36 
+    out@parameters <- parameters
37 
+    out
38 
+}
39 
+