Browse code

fixes to plotRegions/plotManyRegions and plot.BSseqTstat

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

khansen authored on 23/05/2013 01:08:15
Showing 5 changed files

... ...
@@ -1,5 +1,5 @@
1 1
 Package: bsseq
2
-Version: 0.9.0
2
+Version: 0.9.1
3 3
 Title: Analyze, manage and store bisulfite sequencing data
4 4
 Description: Tools for analyzing and visualizing bisulfite sequencing data
5 5
 Author: Kasper Daniel Hansen
... ...
@@ -237,9 +237,14 @@ print.summary.BSseqTstat <- function(x, ...) {
237 237
 
238 238
 plot.BSseqTstat <- function(x, y, ...) {
239 239
     tstat <- getStats(x)[, "tstat"]
240
-    tstat.cor <- getStats(x)[, "tstat.corrected"]
241 240
     plot(density(tstat), xlim = c(-10,10), col = "blue", main = "")
242
-    lines(density(tstat.cor), col = "black")
243
-    legend("topleft", legend = c("uncorrected", "corrected"), lty = c(1,1),
244
-           col = c("blue", "black"))
241
+    if("tstat.corrected" %in% colnames(getStats(x))) {
242
+        tstat.cor <- getStats(x)[, "tstat.corrected"]
243
+        lines(density(tstat.cor), col = "black")
244
+        legend("topleft", legend = c("uncorrected", "corrected"), lty = c(1,1),
245
+               col = c("blue", "black"))
246
+    } else {
247
+        legend("topleft", legend = c("uncorrected"), lty = 1,
248
+               col = c("blue"))
249
+    }
245 250
 }
... ...
@@ -23,22 +23,25 @@ dmrFinder <- function(BSseqTstat, cutoff = NULL, qcutoff = c(0.025, 0.975),
23 23
     regions
24 24
 }
25 25
 
26
-clusterMaker <- function(chr,pos,order.it=TRUE,maxGap=300){
27
-  nonaIndex=which(!is.na(chr) & !is.na(pos))
28
-  Indexes=split(nonaIndex,chr[nonaIndex])
29
-  clusterIDs=rep(NA,length(chr))
30
-  LAST=0
31
-  for(i in seq(along=Indexes)){
32
-    Index=Indexes[[i]]
33
-    x=pos[Index]
34
-    if(order.it){ Index=Index[order(x)];x=pos[Index] }
35
-    y=as.numeric(diff(x)>maxGap)
36
-    z=cumsum(c(1,y))
37
-    clusterIDs[Index]=z+LAST
38
-    LAST=max(z)+LAST
39
-  }
40
-  clusterIDs
41
-}  
26
+clusterMaker <- function(chr, pos, order.it=TRUE, maxGap=300){
27
+    nonaIndex <- which(!is.na(chr) & !is.na(pos))
28
+    Indexes <- split(nonaIndex, chr[nonaIndex])
29
+    clusterIDs <- rep(NA, length(chr))
30
+    LAST <- 0
31
+    for(i in seq(along = Indexes)){
32
+        Index <- Indexes[[i]]
33
+        x <- pos[Index]
34
+        if(order.it){
35
+            Index <- Index[order(x)]
36
+            x <- pos[Index]
37
+        }
38
+        y <- as.numeric(diff(x) > maxGap)
39
+        z <- cumsum(c(1, y))
40
+        clusterIDs[Index] <- z + LAST
41
+        LAST <- max(z) + LAST
42
+    }
43
+    clusterIDs
44
+}
42 45
 
43 46
 
44 47
 regionFinder3 <- function(x, chr, positions, keep, maxGap = 300, verbose = TRUE) {
... ...
@@ -13,10 +13,10 @@ plotAnnoTrack <- function(gr, annoTrack) {
13 13
         })
14 14
 }
15 15
 
16
-plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addRegions = NULL, annoTrack = NULL, 
17
-                            col = NULL, lty = NULL, lwd = NULL, BSseqTstat = NULL, mainWithWidth = TRUE,
18
-                            regionCol = alpha("red", 0.1), addTicks = TRUE, addPoints = FALSE,
19
-                            pointsMinCov = 5, highlightMain = FALSE, verbose = TRUE) {
16
+plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addRegions = NULL,
17
+                            annoTrack = NULL, col = NULL, lty = NULL, lwd = NULL, BSseqTstat = NULL,
18
+                            mainWithWidth = TRUE, regionCol = alpha("red", 0.1), addTicks = TRUE,
19
+                            addPoints = FALSE, pointsMinCov = 5, highlightMain = FALSE, verbose = TRUE) {
20 20
     cat("preprocessing ...")
21 21
     if(!is.null(regions)) {
22 22
         if(is(regions, "data.frame"))
... ...
@@ -69,6 +69,7 @@ plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addReg
69 69
     ## regions <- pintersect(region, rep(gr, length(regions)))
70 70
     ## regions <- regions[width(regions) == 0]
71 71
     regions <- subsetByOverlaps(regions, gr)
72
+    regions <- pintersect(regions, rep(gr, length(regions)))
72 73
     if(length(regions) == 0)
73 74
         return(NULL)
74 75
     rect(xleft = start(regions), xright = end(regions), ybottom = ylim[1],
... ...
@@ -115,7 +116,7 @@ plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addReg
115 116
                            format(plotRange[2], big.mark = ",", scientific = FALSE))
116 117
     if(mainWithWidth) {
117 118
         regionWidth <- sprintf("width = %s, extended = %s", 
118
-                               format(width(gr) - 2*extend, big.mark = ",", scientific = FALSE),
119
+                               format(width(gr), big.mark = ",", scientific = FALSE),
119 120
                                format(extend, big.mark = ",", scientific = FALSE))
120 121
         regionCoord <- sprintf("%s (%s)", regionCoord, regionWidth)
121 122
     }
... ...
@@ -2,6 +2,14 @@
2 2
 \title{bsseq news}
3 3
 \encoding{UTF-8}
4 4
 
5
+\section{Version 0.7.x}{
6
+  \itemize{
7
+    \item Fixed a problem with "width" in the title of bsseq plots.
8
+    \item plot.BSseqTstat now allows for BSseqTstat objects computed
9
+    without correction.
10
+  }
11
+}
12
+
5 13
 \section{Version 0.7.x}{
6 14
   \itemize{
7 15
     \item Removed the returnRaw argument to read.bismark() as it was