Browse code

update to the plotting functions

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

Kasper D. Hansen authored on 12/11/2012 14:46:59
Showing 1 changed files

... ...
@@ -62,7 +62,7 @@ plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addReg
62 62
 
63 63
 .bsHighlightRegions <- function(regions, gr, ylim, regionCol, highlightMain) {
64 64
     if(is.data.frame(regions))
65
-        regions <- data.frame2Granges(regions)
65
+        regions <- data.frame2GRanges(regions)
66 66
     if(highlightMain)
67 67
         regions <- c(regions, gr)
68 68
     if(is.null(regions)) return(NULL)
... ...
@@ -77,30 +77,57 @@ plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addReg
77 77
 
78 78
 .bsGetCol <- function(object, col, lty, lwd) {
79 79
     ## Assumes that object has pData and sampleNames methods
80
-    if(is.null(col) & "col" %in% names(pData(object)))
80
+    if(is.null(col) && "col" %in% names(pData(object)))
81 81
         col <- pData(object)[["col"]]
82 82
     else
83
-        col <- rep("black", ncol(pData(object)))
83
+        col <- rep("black", nrow(pData(object)))
84 84
     if(is.null(names(col)))
85 85
         names(col) <- sampleNames(object)
86 86
     
87
-    if(is.null(lty) & "lty" %in% names(pData(object)))
87
+    if(is.null(lty) && "lty" %in% names(pData(object)))
88 88
         lty <- pData(object)[["lty"]]
89 89
     else
90
-        lty <- rep(1, ncol(pData(object)))
90
+        lty <- rep(1, nrow(pData(object)))
91 91
     if(is.null(names(lty)))
92 92
         names(lty) <- sampleNames(object)
93 93
     
94
-    if(is.null(lwd) & "lwd" %in% names(pData(object)))
94
+    if(is.null(lwd) && "lwd" %in% names(pData(object)))
95 95
         lwd <- pData(object)[["lwd"]]
96 96
     else
97
-        lwd <- rep(1, ncol(object))
97
+        lwd <- rep(1, nrow(pData(object)))
98 98
     if(is.null(names(lwd)))
99 99
         names(lwd) <- sampleNames(object)
100 100
                    
101 101
     return(list(col = col, lty = lty, lwd = lwd))
102 102
 }
103 103
 
104
+.bsPlotTitle <- function(gr, extend, main, mainWithWidth) {
105
+    if(is.data.frame(gr))
106
+        gr <- data.frame2GRanges(gr)
107
+    if(length(gr) > 1) {
108
+        warning("plotTitle: gr has more than one element")
109
+        gr <- gr[1]
110
+    }
111
+    plotChr <- as.character(seqnames(gr))
112
+    plotRange <- c(start(gr), end(gr))
113
+    regionCoord <- sprintf("%s: %s - %s", plotChr, 
114
+                           format(plotRange[1], big.mark = ",", scientific = FALSE),
115
+                           format(plotRange[2], big.mark = ",", scientific = FALSE))
116
+    if(mainWithWidth) {
117
+        regionWidth <- sprintf("width = %s, extended = %s", 
118
+                               format(width(gr) - 2*extend, big.mark = ",", scientific = FALSE),
119
+                               format(extend, big.mark = ",", scientific = FALSE))
120
+        regionCoord <- sprintf("%s (%s)", regionCoord, regionWidth)
121
+    }
122
+    if(main != "") {
123
+        main <- sprintf("%s\n%s", main, regionCoord)
124
+    } else {
125
+        main <- regionCoord
126
+    }
127
+    main
128
+}
129
+
130
+
104 131
 .plotSmoothData <- function(BSseq, region, extend, addRegions, col, lty, lwd, regionCol,
105 132
                             addTicks, addPoints, pointsMinCov, highlightMain) {
106 133
     if(is.data.frame(region))
... ...
@@ -128,8 +155,7 @@ plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addReg
128 155
     smoothPs <- getMeth(BSseq, type = "smooth")
129 156
     rawPs <- getMeth(BSseq, type = "raw")
130 157
     coverage <- getCoverage(BSseq)
131
-    plotRange
132
-    
158
+        
133 159
     ## get col, lwd, lty
134 160
     colEtc <- bsseq:::.bsGetCol(object = BSseq, col = col, lty = lty, lwd = lwd)
135 161
     
... ...
@@ -140,7 +166,7 @@ plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addReg
140 166
     if(addTicks)
141 167
         rug(positions)
142 168
 
143
-    .bsHighlightRegions(regions = regions, gr = gr, ylim = c(0,1),
169
+    .bsHighlightRegions(regions = addRegions, gr = gr, ylim = c(0,1),
144 170
                         regionCol = regionCol, highlightMain = highlightMain)
145 171
     
146 172
     if(addPoints) {
... ...
@@ -164,7 +190,7 @@ plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addReg
164 190
 }
165 191
 
166 192
 
167
-newPlotRegion <- function(BSseq, region = NULL, extend = 0, main = "", addRegions = NULL, annoTrack = NULL,
193
+plotRegion <- function(BSseq, region = NULL, extend = 0, main = "", addRegions = NULL, annoTrack = NULL,
168 194
                           col = NULL, lty = NULL, lwd = NULL, BSseqTstat = NULL, mainWithWidth = TRUE,
169 195
                           regionCol = alpha("red", 0.1), addTicks = TRUE, addPoints = FALSE,
170 196
                           pointsMinCov = 5, highlightMain = FALSE) {
... ...
@@ -200,198 +226,14 @@ newPlotRegion <- function(BSseq, region = NULL, extend = 0, main = "", addRegion
200 226
         bsseq:::plotAnnoTrack(gr, annoTrack)
201 227
 
202 228
     if(!is.null(main)) {
203
-        main <- makePlotTitle(gr = gr, extend = extend, main = main, mainWithWidth = mainWithWidth)
204
-        mtext(side = 3, text = main, outer = TRUE, cex = 1)
205
-    }
206
-    return(invisible(NULL))
207
-}
208
-
209
-    
210
-
211
-plotRegion <- function(BSseq, region = NULL, extend = 0, main = "", addRegions = NULL, annoTrack = NULL,
212
-                       col = NULL, lty = NULL, lwd = NULL, BSseqTstat = NULL, mainWithWidth = TRUE,
213
-                       regionCol = alpha("red", 0.1), addTicks = TRUE, addPoints = FALSE,
214
-                       pointsMinCov = 5, highlightMain = FALSE) {
215
-    makeTitle <- function(gr, extend, main, mainWithWidth) {
216
-        if(length(gr) > 1) {
217
-            warning("plotTitle: gr has more than one element")
218
-            gr <- gr[1]
219
-        }
220
-        plotChr <- as.character(seqnames(gr))
221
-        plotRange <- c(start(gr), end(gr))
222
-        regionCoord <- sprintf("%s: %s - %s", plotChr, 
223
-                               format(plotRange[1], big.mark = ",", scientific = FALSE),
224
-                               format(plotRange[2], big.mark = ",", scientific = FALSE))
225
-        if(mainWithWidth) {
226
-            regionWidth <- sprintf("width = %s, extended = %s", 
227
-                                   format(width(gr) - 2*extend, big.mark = ",", scientific = FALSE),
228
-                                   format(extend, big.mark = ",", scientific = FALSE))
229
-            regionCoord <- sprintf("%s (%s)", regionCoord, regionWidth)
230
-        }
231
-        if(main != "") {
232
-            main <- sprintf("%s\n%s", main, regionCoord)
233
-        } else {
234
-            main <- regionCoord
235
-        }
236
-        main
237
-    }
238
-    plotRects <- function(ylim) {
239
-        if(!is.null(addRegions))
240
-            rect(xleft = addRegions$start, xright = addRegions$end, ybottom = ylim[1],
241
-                 ytop = ylim[2], col = regionCol, border = NA)
242
-    }
243
-    restrictRegions <- function(regions, plotRange, plotChr) {
244
-        if(is.null(regions)) return(NULL)
245
-        regions <- regions[regions$chr == plotChr &
246
-                           ((regions$start >= plotRange[1] &
247
-                             regions$start <= plotRange[2]) |
248
-                           (regions$end >= plotRange[1] &
249
-                            regions$end <= plotRange[2])),, drop = FALSE]
250
-        if(nrow(regions) == 0)
251
-            regions <- NULL
252
-        regions
253
-    }
254
-    plotLines <- function(x, y, lty, col, lwd, plotRange) {
255
-        if(sum(!is.na(y)) <= 1)
256
-            return(NULL)
257
-        xx <- seq(from = plotRange[1], to = plotRange[2], length.out = 2000)
258
-        yy <- approxfun(x, y)(xx)
259
-        lines(xx, yy, col = col, lty = lty, lwd = lwd)
260
-    }
261
-    plotPoints <- function(x, y, z, col) {
262
-        points(x[z>pointsMinCov], y[z>pointsMinCov], col = col, pch = 16, cex = 0.5)
263
-        ## sample, label = "", col) {
264
-        ## plot(positions[1], 0.5, type = "n", xaxt = "n", yaxt = "n",
265
-        ##      ylim = c(0,1), xlim = plotRange, xlab = "", ylab = "")
266
-        ## plotRects(c(0,1))
267
-        ## rawp <- getP(BSseq, sample = sample, type = "raw", addPositions = TRUE, addConfint = TRUE)
268
-        ## cols <- rep(alpha("black", 0.5), nrow(rawp))
269
-        ## segments(x0 = rawp$pos, y0 = rawp$lower, y1 = rawp$upper, col = cols)
270
-        ## points(positions, rawp$p, col = cols)
271
-        ## if(nrow(BSseq$coef) > 0) {
272
-        ##     fitp <- getP(BSseq, sample = sample, type = "fit", addPosition = TRUE, addConfint = FALSE) 
273
-        ##     lines(fitp$pos, fitp$p, col = col,, lty = 1)
274
-        ##     ## lines(fitp$pos, fitp$lower, col = col, lty = 2)
275
-        ##     ## lines(fitp$pos, fitp$upper, col = col, lty = 2)
276
-        ##     text(plotRange[1], 0.1, labels = label)
277
-        ## }
278
-    }
279
-    
280
-    ## First we create a basic GRanges which will be the plotting region
281
-    if(!is.null(region)) {
282
-        if(is(region, "data.frame"))
283
-            gr <- data.frame2GRanges(region, keepColumns = FALSE)
284
-        else
285
-            gr <- region
286
-        if(!is(gr, "GRanges") || length(gr) != 1)
287
-            stop("'region' needs to be either a 'data.frame' (with a single row) or a 'GRanges' (with a single element)")
288
-    } else {
289
-        gr <- GRanges(seqnames = seqnames(BSseq)[1],
290
-                      ranges = IRanges(start = min(start(BSseq)),
291
-                      end = max(start(BSseq))))
292
-    }
293
-    origWidth <- width(gr)
294
-    gr <- resize(gr, width = 2*extend + width(gr), fix = "center")
295
-    plotRange <- c(start(gr), end(gr))
296
-    plotChr <- as.character(seqnames(gr))[1]
297
-    BSseq <- subsetByOverlaps(BSseq, gr)
298
-    if(!is.null(BSseqTstat))
299
-        BSseqTstat <- subsetByOverlaps(BSseqTstat, gr)
300
-    positions <- start(BSseq)
301
-    if(length(positions) == 0) {
302
-        warning("No overlap between BSseq data and region")
303
-        return(NULL)
304
-    }
305
-        
306
-    ## Now for some plotting
307
-    opar <- par(mar = c(0,4.1,0,0), oma = c(5,0,4,2), mfrow = c(1,1))
308
-    on.exit(par(opar))
309
-    if(is.null(BSseqTstat))
310
-        layout(matrix(1:2, ncol = 1), heights = c(2,1))
311
-    else
312
-        layout(matrix(1:3, ncol = 1), heights = c(2,2,1))
313
-       
314
-    sampleNames <- sampleNames(BSseq)
315
-    names(sampleNames) <- sampleNames
316
-    plot(positions[1], 0.5, type = "n", xaxt = "n", yaxt = "n",
317
-         ylim = c(0,1), xlim = plotRange, xlab = "", ylab = "Methylation")
318
-    axis(side = 2, at = c(0.2, 0.5, 0.8))
319
-    if(addTicks)
320
-        rug(positions)
321
-
322
-    addRegions <- restrictRegions(addRegions, plotRange = plotRange, plotChr = plotChr)
323
-    if(highlightMain)
324
-        addRegions <- rbind(region[, c("chr", "start", "end")],
325
-                            addRegions[, c("chr", "start", "end")])
326
-    if(!is.null(addRegions))
327
-        plotRects(c(0,1))
328
-
329
-    smoothPs <- getMeth(BSseq, type = "smooth")
330
-    rawPs <- getMeth(BSseq, type = "raw")
331
-    coverage <- getCoverage(BSseq)
332
-    
333
-    if(is.null(col) & "col" %in% names(pData(BSseq)))
334
-        col <- pData(BSseq)[["col"]]
335
-    else
336
-        col <- rep("black", ncol(BSseq))
337
-    if(is.null(names(col)))
338
-        names(col) <- sampleNames(BSseq)
339
-    if(is.null(lty) & "lty" %in% names(pData(BSseq)))
340
-        lty <- pData(BSseq)[["lty"]]
341
-    else
342
-        lty <- rep(1, ncol(BSseq))
343
-    if(is.null(names(lty)))
344
-        names(lty) <- sampleNames(BSseq)
345
-    if(is.null(lwd) & "lwd" %in% names(pData(BSseq)))
346
-        lwd <- pData(BSseq)[["lwd"]]
347
-    else
348
-        lwd <- rep(1, ncol(BSseq))
349
-    if(is.null(names(lwd)))
350
-        names(lwd) <- sampleNames(BSseq)
351
-
352
-    if(addPoints) {
353
-        sapply(sampleNames(BSseq), function(samp) {
354
-            abline(v = positions[rawPs[, samp] > 0.1], col = "grey80", lty = 1)
355
-        })
356
-    }
357
-
358
-    sapply(sampleNames(BSseq), function(samp) {
359
-        plotLines(positions, smoothPs[, samp], col = col[samp],
360
-                  lty = lty[samp], lwd = lwd[samp], plotRange = plotRange)
361
-    })
362
-
363
-    if(addPoints) {
364
-        sapply(sampleNames(BSseq), function(samp) {
365
-            plotPoints(positions, rawPs[, samp], coverage[, samp], col = col[samp])
366
-        })
367
-    }
368
-    
369
-    if(!is.null(BSseqTstat)) {
370
-        plot(positions[1], 0.5, type = "n", xaxt = "n", yaxt = "n",
371
-             ylim = c(-8,8), xlim = plotRange, xlab = "", ylab = "t-stat")
372
-        axis(side = 2, at = c(-5,0,5))
373
-        abline(h = 0, col = "grey60")
374
-        plotLines(start(BSseqTstat), BSseqTstat@stats[, "tstat"],
375
-                  lty = 1, plotRange = plotRange, col = "red", lwd = 1)
376
-        plotLines(start(BSseqTstat), BSseqTstat@stats[, "tstat.corrected"],
377
-                  lty = 2, plotRange = plotRange, col = "red", lwd = 1)
378
-        plotLines(start(BSseqTstat), 100*BSseqTstat@stats[, "tstat.sd"],
379
-                  lty = 2, plotRange = plotRange, col = "blue", lwd = 1)
380
-    }
381
-    
382
-
383
-    if(!is.null(annoTrack))
384
-        bsseq:::plotAnnoTrack(gr, annoTrack)
385
-
386
-    if(!is.null(main)) {
387
-        main <- makePlotTitle(gr = gr, extend = extend, main = main, mainWithWidth = mainWithWidth)
229
+        main <- bsseq:::.bsPlotTitle(gr = region, extend = extend, main = main,
230
+                                     mainWithWidth = mainWithWidth)
388 231
         mtext(side = 3, text = main, outer = TRUE, cex = 1)
389 232
     }
390 233
     return(invisible(NULL))
391 234
 }
392 235
 
393
-
394
-
236
+ 
395 237
 ##     plotP <- function(sample, label = "", col) {
396 238
 ##         plot(positions[1], 0.5, type = "n", xaxt = "n", yaxt = "n",
397 239
 ##              ylim = c(0,1), xlim = plotRange, xlab = "", ylab = "")