Browse code

reorganization of plotting code

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

Kasper D. Hansen authored on 12/11/2012 02:08:38
Showing1 changed files

... ...
@@ -48,12 +48,193 @@ plotManyRegions <- function(BSseq, regions = NULL, extend = 0, main = "", addReg
48 48
     }
49 49
 }
50 50
 
51
+.bsPlotLines <- function(x, y, col, lty, lwd, plotRange) {
52
+    if(sum(!is.na(y)) <= 1)
53
+        return(NULL)
54
+    xx <- seq(from = plotRange[1], to = plotRange[2], length.out = 500)
55
+    yy <- approxfun(x, y)(xx)
56
+    lines(xx, yy, col = col, lty = lty, lwd = lwd)
57
+}
58
+
59
+.bsPlotPoints <- function(x, y, z, col) {
60
+    points(x[z>pointsMinCov], y[z>pointsMinCov], col = col, pch = 16, cex = 0.5)
61
+}
62
+
63
+.bsHighlightRegions <- function(regions, gr, ylim, regionCol, highlightMain) {
64
+    if(is.data.frame(regions))
65
+        regions <- data.frame2Granges(regions)
66
+    if(highlightMain)
67
+        regions <- c(regions, gr)
68
+    if(is.null(regions)) return(NULL)
69
+    ## regions <- pintersect(region, rep(gr, length(regions)))
70
+    ## regions <- regions[width(regions) == 0]
71
+    regions <- subsetByOverlaps(regions, gr)
72
+    if(length(regions) == 0)
73
+        return(NULL)
74
+    rect(xleft = start(regions), xright = end(regions), ybottom = ylim[1],
75
+         ytop = ylim[2], col = regionCol, border = NA)
76
+}
77
+
78
+.bsGetCol <- function(object, col, lty, lwd) {
79
+    ## Assumes that object has pData and sampleNames methods
80
+    if(is.null(col) & "col" %in% names(pData(object)))
81
+        col <- pData(object)[["col"]]
82
+    else
83
+        col <- rep("black", ncol(pData(object)))
84
+    if(is.null(names(col)))
85
+        names(col) <- sampleNames(object)
86
+    
87
+    if(is.null(lty) & "lty" %in% names(pData(object)))
88
+        lty <- pData(object)[["lty"]]
89
+    else
90
+        lty <- rep(1, ncol(pData(object)))
91
+    if(is.null(names(lty)))
92
+        names(lty) <- sampleNames(object)
93
+    
94
+    if(is.null(lwd) & "lwd" %in% names(pData(object)))
95
+        lwd <- pData(object)[["lwd"]]
96
+    else
97
+        lwd <- rep(1, ncol(object))
98
+    if(is.null(names(lwd)))
99
+        names(lwd) <- sampleNames(object)
100
+                   
101
+    return(list(col = col, lty = lty, lwd = lwd))
102
+}
51 103
 
104
+.plotSmoothData <- function(BSseq, region, extend, addRegions, col, lty, lwd, regionCol,
105
+                            addTicks, addPoints, pointsMinCov, highlightMain) {
106
+    if(is.data.frame(region))
107
+        region <- data.frame2GRanges(region)
108
+
109
+    if(!is.null(region)) {
110
+        if(is(region, "data.frame"))
111
+            gr <- data.frame2GRanges(region, keepColumns = FALSE)
112
+        else
113
+            gr <- region
114
+        if(!is(gr, "GRanges") || length(gr) != 1)
115
+            stop("'region' needs to be either a 'data.frame' (with a single row) or a 'GRanges' (with a single element)")
116
+    } else {
117
+        gr <- GRanges(seqnames = seqnames(BSseq)[1],
118
+                      ranges = IRanges(start = min(start(BSseq)),
119
+                      end = max(start(BSseq))))
120
+    }
121
+    gr <- resize(gr, width = 2*extend + width(gr), fix = "center")
122
+    BSseq <- subsetByOverlaps(BSseq, gr)
123
+    
124
+    ## Extract basic information
125
+    sampleNames <- sampleNames(BSseq)
126
+    names(sampleNames) <- sampleNames
127
+    positions <- start(BSseq)
128
+    smoothPs <- getMeth(BSseq, type = "smooth")
129
+    rawPs <- getMeth(BSseq, type = "raw")
130
+    coverage <- getCoverage(BSseq)
131
+    plotRange
132
+    
133
+    ## get col, lwd, lty
134
+    colEtc <- bsseq:::.bsGetCol(object = BSseq, col = col, lty = lty, lwd = lwd)
135
+    
136
+    ## The actual plotting
137
+    plot(positions[1], 0.5, type = "n", xaxt = "n", yaxt = "n",
138
+         ylim = c(0,1), xlim = c(start(gr), end(gr)), xlab = "", ylab = "Methylation")
139
+    axis(side = 2, at = c(0.2, 0.5, 0.8))
140
+    if(addTicks)
141
+        rug(positions)
142
+
143
+    .bsHighlightRegions(regions = regions, gr = gr, ylim = c(0,1),
144
+                        regionCol = regionCol, highlightMain = highlightMain)
145
+    
146
+    if(addPoints) {
147
+        sapply(sampleNames(BSseq), function(samp) {
148
+            abline(v = positions[rawPs[, samp] > 0.1], col = "grey80", lty = 1)
149
+        })
150
+    } # This adds vertical grey lines so we can see where points are plotted
151
+
152
+    sapply(sampleNames(BSseq), function(samp) {
153
+        .bsPlotLines(positions, smoothPs[, samp], col = colEtc$col[samp],
154
+                     lty = colEtc$lty[samp], lwd = colEtc$lwd[samp],
155
+                     plotRange = c(start(gr), end(gr)))
156
+    })
157
+
158
+    if(addPoints) {
159
+        sapply(sampleNames(BSseq), function(samp) {
160
+            .bsPlotPoints(positions, rawPs[, samp], coverage[, samp],
161
+                          col = colEtc$col[samp])
162
+        })
163
+    }
164
+}
165
+
166
+
167
+newPlotRegion <- function(BSseq, region = NULL, extend = 0, main = "", addRegions = NULL, annoTrack = NULL,
168
+                          col = NULL, lty = NULL, lwd = NULL, BSseqTstat = NULL, mainWithWidth = TRUE,
169
+                          regionCol = alpha("red", 0.1), addTicks = TRUE, addPoints = FALSE,
170
+                          pointsMinCov = 5, highlightMain = FALSE) {
171
+    
172
+    opar <- par(mar = c(0,4.1,0,0), oma = c(5,0,4,2), mfrow = c(1,1))
173
+    on.exit(par(opar))
174
+    if(is.null(BSseqTstat))
175
+        layout(matrix(1:2, ncol = 1), heights = c(2,1))
176
+    else
177
+        layout(matrix(1:3, ncol = 1), heights = c(2,2,1))
178
+
179
+    bsseq:::.plotSmoothData(BSseq = BSseq, region = region, extend = extend, addRegions = addRegions,
180
+                            col = col, lty = lty, lwd = lwd, regionCol = regionCol,
181
+                            addTicks = addTicks, addPoints = addPoints,
182
+                            pointsMinCov = pointsMinCov, highlightMain = highlightMain)
183
+
184
+    if(!is.null(BSseqTstat)) {
185
+        if(!is.null(BSseqTstat))
186
+            BSseqTstat <- subsetByOverlaps(BSseqTstat, gr)
187
+        plot(positions[1], 0.5, type = "n", xaxt = "n", yaxt = "n",
188
+             ylim = c(-8,8), xlim = plotRange, xlab = "", ylab = "t-stat")
189
+        axis(side = 2, at = c(-5,0,5))
190
+        abline(h = 0, col = "grey60")
191
+        bsseq:::.bsPlotLines(start(BSseqTstat), BSseqTstat@stats[, "tstat"],
192
+                             lty = 1, plotRange = plotRange, col = "red", lwd = 1)
193
+        bsseq:::.bsPlotLines(start(BSseqTstat), BSseqTstat@stats[, "tstat.corrected"],
194
+                             lty = 2, plotRange = plotRange, col = "red", lwd = 1)
195
+        bsseq:::.bsPlotLines(start(BSseqTstat), 100*BSseqTstat@stats[, "tstat.sd"],
196
+                             lty = 2, plotRange = plotRange, col = "blue", lwd = 1)
197
+    }
198
+    
199
+    if(!is.null(annoTrack))
200
+        bsseq:::plotAnnoTrack(gr, annoTrack)
201
+
202
+    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
+    
52 210
 
53 211
 plotRegion <- function(BSseq, region = NULL, extend = 0, main = "", addRegions = NULL, annoTrack = NULL,
54 212
                        col = NULL, lty = NULL, lwd = NULL, BSseqTstat = NULL, mainWithWidth = TRUE,
55 213
                        regionCol = alpha("red", 0.1), addTicks = TRUE, addPoints = FALSE,
56 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
+    }
57 238
     plotRects <- function(ylim) {
58 239
         if(!is.null(addRegions))
59 240
             rect(xleft = addRegions$start, xright = addRegions$end, ybottom = ylim[1],
... ...
@@ -117,29 +298,11 @@ plotRegion <- function(BSseq, region = NULL, extend = 0, main = "", addRegions =
117 298
     if(!is.null(BSseqTstat))
118 299
         BSseqTstat <- subsetByOverlaps(BSseqTstat, gr)
119 300
     positions <- start(BSseq)
120
-    if(length(positions) == 0)
121
-        stop("No overlap between BSseq data and region")
122
-
123
-    ## This next section manages the plot title
124
-    regionCoord <- sprintf("%s: %s - %s", plotChr, 
125
-                               format(plotRange[1], big.mark = ",", scientific = FALSE),
126
-                               format(plotRange[2], big.mark = ",", scientific = FALSE))
127
-    regionWidth <- sprintf("width = %s, extended = %s", 
128
-                           format(origWidth, big.mark = ",", scientific = FALSE),
129
-                           format(extend, big.mark = ",", scientific = FALSE))
130
-    if(mainWithWidth)
131
-        regionCoord <- sprintf("%s (%s)", regionCoord, regionWidth)
132
-    
133
-    if(is.null(main)) {
134
-        main <- ""
135
-    } else {
136
-        if(main != "") {
137
-            main <- sprintf("%s\n%s", main, regionCoord)
138
-        } else {
139
-            main <- regionCoord
140
-        }
301
+    if(length(positions) == 0) {
302
+        warning("No overlap between BSseq data and region")
303
+        return(NULL)
141 304
     }
142
-
305
+        
143 306
     ## Now for some plotting
144 307
     opar <- par(mar = c(0,4.1,0,0), oma = c(5,0,4,2), mfrow = c(1,1))
145 308
     on.exit(par(opar))
... ...
@@ -219,8 +382,11 @@ plotRegion <- function(BSseq, region = NULL, extend = 0, main = "", addRegions =
219 382
 
220 383
     if(!is.null(annoTrack))
221 384
         bsseq:::plotAnnoTrack(gr, annoTrack)
222
-    
223
-    mtext(side = 3, text = main, outer = TRUE, cex = 1)
385
+
386
+    if(!is.null(main)) {
387
+        main <- makePlotTitle(gr = gr, extend = extend, main = main, mainWithWidth = mainWithWidth)
388
+        mtext(side = 3, text = main, outer = TRUE, cex = 1)
389
+    }
224 390
     return(invisible(NULL))
225 391
 }
226 392