Mercurial > repos > jeremyjliu > region_motif_enrichment

comparison region_motif_lib/plotting.r @ 30:da0bc6dead0d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Restore deleted source files
author jeremyjliu
date Mon, 05 Jan 2015 22:28:37 -0500
parents 19d2cffb8db3
children
comparison
equal deleted inserted replaced
29:c78ebf00379d 30:da0bc6dead0d
1 library(graphics, quietly=TRUE)
2
3 plot.verbose=F
4 name.cleaner<-function(...,sep="",replace="_") {
5 plot.name=gsub(" ",replace,paste(...,sep=sep))
6 plot.name=gsub("/",replace,plot.name)
7 plot.name=gsub(",",replace,plot.name)
8 plot.name=gsub("'",replace,plot.name)
9 plot.name=gsub("\\+","plus",plot.name)
10 plot.name=gsub("\\(","",plot.name)
11 plot.name=gsub("\\)","",plot.name)
12 return(plot.name)
13 }
14 plot.namer <- function(..., date=0, fig.dir=0, format="png",sep="") {
15 plot.name=name.cleaner(...,sep=sep)
16 if(date==0) date=gsub("-","",as.character(Sys.Date()))
17 if(fig.dir==0) fig.dir="/Users/alver/allplots"
18 plot.name=paste(fig.dir,"/",date,plot.name,".",format,sep="")
19 if(plot.verbose) cat(" saving figure: ",plot.name,"\n")
20 return(plot.name)
21 }
22
23 plot.scatter <- function(x,y=NULL,f=0.1,same=FALSE,n.points=-1,draw.lowess=FALSE,write.r=TRUE,cex.r=1,col=NULL,col.line=NULL,lwd.line=1,
24 draw.loess=FALSE,span=0.5,bandwidth=bandwidth,draw.prof=FALSE,xlog=FALSE,ylog=FALSE,cor.method="pearson",log="",ind=NULL,
25 draw.spread=FALSE,...) {
26
27 ## if col is the same length as x, use col for each point matching x.
28 ## if col is the same length as ind, use col for each point matching x[ind].
29 ## else use densCols function based on col.
30 ## if col is null, densCols is used with bluetone for first plot and redtone for same=T.
31
32 #print(length(x))
33 #print(length(y))
34
35 xy <- xy.coords(x, y)
36 x=xy$x
37 y=xy$y
38
39 output=list()
40 col.use = col
41
42 if(!is.null(ind)) {
43 if(length(col.use)==length(x)) {
44 col.use=col.use[ind]
45 }
46 x=x[ind]
47 y=y[ind]
48 }
49
50 if(length(col.use)!=length(x)) {
51 col.use=rep(NA,length(x))
52 }
53
54
55 take=which(is.finite(x) & is.finite(y))
56 x=x[take]
57 y=y[take]
58 col.use=col.use[take]
59
60 if(grepl("x",log)) xlog=TRUE
61 if(grepl("y",log)) ylog=TRUE
62 if(xlog) log="x"
63 if(ylog) log=paste(log,"y",sep="")
64
65 if(draw.lowess) {
66 lo = lowess(x,y,f)
67 output$lowess=lo
68 }
69 if(draw.loess | draw.spread) {
70 px=x;py=y
71 if(xlog) px=log(x)
72 if(ylog) py=log(y)
73 ind = which(is.finite(px+py))
74 px=px[ind]
75 py=py[ind]
76 lo = loess(py ~ px,span=span,iterations=5)
77 lo.y=as.numeric(lo$fitted)
78 lo.x=as.numeric(lo$x)
79 if(draw.spread) lo.sd = loess((lo.y-py)^2 ~ lo.x,span=span*1.5,iterations=5)
80 if(xlog) lo.x=exp(lo.x)
81 if(ylog) lo.y=exp(lo.y)
82 lo =data.frame(x=lo.x,y=lo.y)
83 if(draw.spread) {
84 lo.sd=lo.sd$fitted
85 if(ylog) lo.sd=lo.sd*lo.y*lo.y
86 lo$sd=sqrt(pmax(0,lo.sd))
87 }
88 lo=unique(lo)
89 lo = lo[order(lo$x),]
90 output$loess=lo
91 }
92
93 if(draw.prof) {
94 px=x;py=y
95 if(xlog) px=log(x)
96 p=prof(px,py,50)
97 if(xlog) p$x=exp(p$x)
98 output$prof=p
99 }
100
101 r=cor(x,y,method=cor.method)
102 output$cor=r
103 output$cor.method=cor.method
104
105 len=length(x)
106 if(n.points>0 & n.points<len) {
107 take=sample(1:len,n.points)
108 x=x[take]
109 y=y[take]
110 col.use=col.use[take]
111 }
112
113 if(xlog) {
114 ind = which(x>0)
115 x=x[ind]
116 y=y[ind]
117 col.use=col.use[ind]
118 }
119 xcol=x
120 if(xlog) xcol=log(xcol)
121 if(ylog) {
122 ind = which(y>0)
123 x=x[ind]
124 xcol=xcol[ind]
125 y=y[ind]
126 col.use=col.use[ind]
127 }
128 ycol=y
129 if(ylog) ycol=log(ycol)
130
131 if(is.null(col)) {
132 if(!same) {
133 col=colorRampPalette(blues9[-(1:3)])
134 } else {
135 col=colorRampPalette(c("lightpink","red","darkred"))
136 }
137 }
138 if(!is.na(col.use[1])) {
139 col=col.use
140 } else {
141 col= suppressPackageStartupMessages(densCols(xcol,ycol,col =col,bandwidth=bandwidth,nbin=500))
142 }
143 if(!same) {
144 plot(x,y,col=col,log=log,...)
145 } else {
146 points(x,y,col=col,...)
147 }
148
149 if(is.null(col.line)) {
150 col.line="darkblue"
151 if(same) col.line="darkred"
152 }
153 if(draw.lowess | draw.loess) lines(lo,col=col.line,lwd=lwd.line)
154 if(draw.spread) {
155 lines(lo$x,lo$y+lo$sd,col=col.line,lwd=lwd.line)
156 lines(lo$x,lo$y-lo$sd,col=col.line,lwd=lwd.line)
157 }
158 if(draw.prof) {
159 points(p)
160 plot.prof(p)
161 }
162 if(write.r & !same) mtext(paste("r=",round(r,3),sep=""),cex=cex.r)
163 return(invisible(output))
164 }
165
166 #color.int=c(144,586,465,257,490,100,74,24)
167 #coli=1
168 #cols = integer()
169 colramp.bwr = vector()
170 colramp.byr = vector()
171 colramp.bw = vector()
172 colramp.bw2 = vector()
173
174 plot.save=F
175
176 setup.plotting <- function() {
177 pdf.options(useDingbats = FALSE)
178 # cols<<-colors()[color.int]
179 # cols<<-rep(cols,100)
180 colramp.bwr <<- colorRampPalette(c("blue","white","red"),space="Lab")(100);
181 colramp.byr <<- colorRampPalette(c("blue","yellow","red"),space="Lab")(100);
182 colramp.bw <<- colorRampPalette(c("white","black"),space="Lab")(100)
183 colramp.bw2 <<- colorRampPalette(c("grey92","grey5"),space="Lab")(100)
184 }
185
186
187 plot.cluster <- function(x,k, max.points.cl=-1, image.sep=-1,col=NULL, reorder=FALSE) {
188 x[which(is.na(x))]=0
189 if(reorder) {
190 o=hclust(dist(t(x)))$order
191 x=x[,o]
192 }
193 if(image.sep<0) {
194 if(max.points.cl>0) {
195 image.sep=ceiling(0.2*max.points.cl)
196 } else {
197 image.sep=ceiling(0.2 * nrow(x) / nrow(k$centers))
198 }
199 }
200
201 distances<-dist(k$centers)
202 hcl=hclust(distances)
203
204 adjust.branch.sep <-function(ddr,lengths) {
205 assign.branch.sep <- function(d,i.leaf) {
206 if(is.leaf(d)) {
207 attr(d,"members")<-lengths[i.leaf]
208 i.leaf=i.leaf+1
209 output=list(d=d,i.leaf=i.leaf)
210 return(output)
211 }
212 else{
213 input=assign.branch.sep(d[[1]],i.leaf)
214 i.leaf=input$i.leaf
215 d[[1]]=input$d
216
217 input=assign.branch.sep(d[[2]],i.leaf)
218 i.leaf=input$i.leaf
219 d[[2]]=input$d
220
221 attr(d,"members")<-attr(d[[1]],"members")+attr(d[[2]],"members")
222 output=list(d=d,i.leaf=i.leaf)
223 return(output)
224 }
225 }
226 ddr<-as.dendrogram(ddr)
227 ddr=assign.branch.sep(ddr,1)$d
228 return(ddr)
229 }
230
231 n.points.actual=k$size
232 if(max.points.cl>0) {
233 k$size[which(k$size>max.points.cl)] = max.points.cl
234 }
235
236 ddr<-adjust.branch.sep(hcl,k$size[hcl$order]+image.sep)
237 centers=length(hcl$order)
238
239 n.points=sum(k$size)
240 n.dims=ncol(x)
241 z=matrix(numeric((n.points+(centers-1)*image.sep)*n.dims),ncol=n.dims)
242
243
244 last.row=0
245 cluster.y.pos=numeric(centers)
246 for(i.c in hcl$order) {
247 n.p=k$size[i.c]
248 z[last.row+1:n.p,] = x[which(k$cluster==i.c)[1:n.p],]
249 cluster.y.pos[i.c]=last.row+n.p/2
250 last.row=last.row+n.p+image.sep
251 }
252
253 zlim=c(0,max(z))
254 if(min(z)<0) {
255 m=max(c(z,-z))
256 zlim=c(-m,m)
257 }
258 if(is.null(col)) {
259 if(min(z)>=0) {
260 col= colramp.bw
261 } else {
262 col= colorRampPalette(c("blue","yellow","red"),space="Lab")(100);
263 }
264 }
265 x.pl=seq1(n.dims+1)-0.5
266 y.pl=seq1(nrow(z)+1)-0.5
267 l <- layout(matrix(1:2,ncol=2),widths=c(1,5))
268 par(mar = c(6,0.5,6,0))
269 my.plot.dendrogram(ddr,horiz=T,axes=F,yaxs="i",xaxs="i",leaflab="none",center=T,lwd=10)
270 par(mar = c(6,0.1,6,2.1))
271 image(x=x.pl,y=y.pl,z=t(z),zlim=zlim,axes=FALSE,xlab="",col=col)
272 mtext("cluster",side=4,adj=1.1)
273 mtext("points",side=4,adj=1.1,line=1)
274 mtext(seq1(centers),side=4,at=cluster.y.pos)
275 mtext(n.points.actual,side=4,at=cluster.y.pos,line=1)
276
277 if(!is.null(dimnames(x)[[2]])) {
278 mtext(dimnames(x)[[2]],side=1,at=seq1(n.dims),las=2)
279 }
280 }
281
282 plot.cluster2 <- function(k, n.clusters=-1, n.clusters.per.panel=4, cols=c("black","red","blue","darkgreen","orange"),f=0,xshift=0,...) {
283 if(n.clusters<=0) n.clusters=nrow(k$centers)
284
285 n.elements=as.numeric(unlist(lapply(seq1(n.clusters), function(cl) length(which(abs(k$cluster)==cl)))))
286
287 distances<-dist(k$centers)
288 n.panels = ceiling(n.clusters/n.clusters.per.panel)
289 n.rows=ceiling(sqrt(n.panels))
290 n.cols=ceiling(n.panels/n.rows)
291 n.panels.layout=n.rows*n.cols
292
293 layout(matrix(seq1(n.panels.layout),nrow=n.rows,byrow=TRUE))
294
295 min=min(k$centers)
296 max=max(k$centers)
297
298 if(f>0) {
299 for(i.cluster in seq1(n.clusters)) {
300 k$centers[i.cluster,]=lowess(k$centers[i.cluster,],f=f)$y
301 }
302 }
303
304 ## hcl=hclust(distances)
305 hcl=list()
306 hcl$order=1:n.clusters
307
308 for(i.cluster in seq1(n.clusters)) {
309 if(i.cluster %% n.clusters.per.panel == 1 ) {
310 clusters.of.panel=i.cluster:(i.cluster+n.clusters.per.panel-1)
311 clusters.of.panel=clusters.of.panel[which(clusters.of.panel<=n.clusters)]
312 clusters.of.panel=hcl$order[clusters.of.panel]
313 plot(c(0,length(k$centers[1,]))+xshift,c(min,max),type="n",...)
314 mtext(paste(clusters.of.panel," (",n.elements[clusters.of.panel],")",sep=""),line=length(clusters.of.panel)-seq1(length(clusters.of.panel)),col=cols[seq1(length(clusters.of.panel)) %% n.clusters.per.panel+1] )
315 }
316 # lines(k$centers[hcl$order[i.cluster],],col=cols[i.cluster %% n.clusters.per.panel+1])
317 lines(seq1(length(k$centers[1,]))+xshift,k$centers[hcl$order[i.cluster],],col=cols[i.cluster %% n.clusters.per.panel+1])
318 }
319 }
320
321 my.colors <- function(n) {
322 few.colors=c("black","red","blue","green3","mediumorchid3","gold2","darkcyan","sienna2")
323 if(n<=length(few.colors)) return(few.colors [seq1(n)])
324 col=integer(n)
325 n.families=7
326 n.members=ceiling(n/n.families)
327 for(i in seq1(n)) {
328 member=ceiling(i/n.families)
329 ratio=(member-1)/(n.members-1)
330 c2=0+0.8*ratio
331 if(member %% 2 == 1) ratio=-ratio
332 c1=0.8-0.2*ratio
333 c3=0.75-0.2*ratio
334 if(i %% n.families == 1) {col[i]=rgb(c2,c2,c2)}
335 if(i %% n.families == 2) {col[i]=rgb(c1,c1/2,c1/2)}
336 if(i %% n.families == 3) {col[i]=rgb(c1/2,0.9*c1,c1/2)}
337 if(i %% n.families == 4) {col[i]=rgb(c1/2,c1/2,c1)}
338 if(i %% n.families == 5) {col[i]=rgb(c3,c3,c3/2)}
339 if(i %% n.families == 6) {col[i]=rgb(c3,c3/2,c3)}
340 if(i %% n.families == 0) {col[i]=rgb(c3/2,c3,c3)}
341 }
342 return(col)
343 }
344
345 plot.my.colors <-function(n) {
346 x11()
347 col=my.colors(n)
348 plot(x=c(0,n),y=c(0,1),type="n")
349 segments(seq1(n)-1,runif(n),seq1(n),runif(n),col=col)
350 }
351
352
353 plot.colors <-function() {
354 x11(width=10,height=10)
355 plot(c(0,26),c(0,26),type="n")
356 c=colors()
357 n=length(c)
358 i=1:n
359 x=i%%26
360 y=floor(i/26)
361 rect(x,y,x+1,y+1,col=c[i],border=c[i])
362 text(x+0.5,y+0.5,i)
363 }
364
365
366 adjust.branch.sep <-function(ddr,lengths) {
367 assign.branch.sep <- function(d,i.leaf) {
368 if(is.leaf(d)) {
369 attr(d,"members")<-lengths[i.leaf]
370 i.leaf=i.leaf+1
371 output=list(d=d,i.leaf=i.leaf)
372 return(output)
373 }
374 else{
375 input=assign.branch.sep(d[[1]],i.leaf)
376 i.leaf=input$i.leaf
377 d[[1]]=input$d
378
379 input=assign.branch.sep(d[[2]],i.leaf)
380 i.leaf=input$i.leaf
381 d[[2]]=input$d
382
383 attr(d,"members")<-attr(d[[1]],"members")+attr(d[[2]],"members")
384 output=list(d=d,i.leaf=i.leaf)
385 return(output)
386 }
387 }
388 ddr<-as.dendrogram(ddr)
389 ddr=assign.branch.sep(ddr,1)$d
390 return(ddr)
391 }
392 t.dhcol <- function(dr,h,cols=c(1)) {
393 # check child heights
394 if(attr(dr[[1]],"height")<h) {
395 # color
396 ecol <- cols[coli];
397 coli <<- coli+1;
398 dr[[1]] <- dendrapply(dr[[1]],function(e) { attr(e,"edgePar") <- list(col=ecol); e});
399 attr(dr[[1]],"edgePar") <- list(col=ecol,p.border=NA,p.col=NA,t.col=1,t.cex=1.3);
400 } else {
401 dr[[1]] <- t.dhcol(dr[[1]],h,cols);
402 }
403
404 if(attr(dr[[2]],"height")<h) {
405 # color
406 ecol <- cols[coli];
407 coli <<- coli+1;
408 dr[[2]] <- dendrapply(dr[[2]],function(e) { attr(e,"edgePar") <- list(col=ecol); e});
409 attr(dr[[2]],"edgePar") <- list(col=ecol,p.border=NA,p.col=NA,t.col=1,t.cex=1.3);
410 } else {
411 dr[[2]] <- t.dhcol(dr[[2]],h,cols);
412 }
413 return(dr);
414 }
415
416
417
418 ### The rest is PeterK's my.plot.dendogram
419
420 ## FIXME: need larger par("mar")[1] or [4] for longish labels !
421 ## {probably don't change, just print a warning ..}
422 my.plot.dendrogram <-
423 function (x, type = c("rectangle", "triangle"), center = FALSE,
424 edge.root = is.leaf(x) || !is.null(attr(x, "edgetext")),
425 nodePar = NULL, edgePar = list(),
426 leaflab = c("perpendicular", "textlike", "none"), dLeaf = NULL,
427 xlab = "", ylab = "", xaxt="n", yaxt="s",
428 horiz = FALSE, frame.plot = FALSE, ...)
429 {
430 type <- match.arg(type)
431 leaflab <- match.arg(leaflab)
432 hgt <- attr(x, "height")
433 if (edge.root && is.logical(edge.root))
434 edge.root <- 0.0625 * if(is.leaf(x)) 1 else hgt
435 mem.x <- .my.memberDend(x)
436 yTop <- hgt + edge.root
437 if(center) { x1 <- 0.5 ; x2 <- mem.x + 0.5 }
438 else { x1 <- 1 ; x2 <- mem.x }
439 xlim <- c(x1 - 1/2, x2 + 1/2)
440 ylim <- c(0, yTop)
441 if (horiz) {## swap and reverse direction on `x':
442 xl <- xlim; xlim <- rev(ylim); ylim <- xl
443 tmp <- xaxt; xaxt <- yaxt; yaxt <- tmp
444 }
445 plot(0, xlim = xlim, ylim = ylim, type = "n", xlab = xlab, ylab = ylab,
446 xaxt = xaxt, yaxt = yaxt, frame.plot = frame.plot, ...)
447 if(is.null(dLeaf))
448 dLeaf <- .75*(if(horiz) strwidth("w") else strheight("x"))
449
450 if (edge.root) {
451 ### FIXME: the first edge + edgetext is drawn here, all others in plotNode()
452 ### ----- maybe use trick with adding a single parent node to the top ?
453 x0 <- my.plotNodeLimit(x1, x2, x, center)$x
454 if (horiz)
455 segments(hgt, x0, yTop, x0)
456 else segments(x0, hgt, x0, yTop)
457 if (!is.null(et <- attr(x, "edgetext"))) {
458 my <- mean(hgt, yTop)
459 if (horiz)
460 text(my, x0, et)
461 else text(x0, my, et)
462 }
463 }
464 my.plotNode(x1, x2, x, type = type, center = center, leaflab = leaflab,
465 dLeaf = dLeaf, nodePar = nodePar, edgePar = edgePar, horiz = horiz)
466 }
467
468 ### the work horse: plot node (if pch) and lines to all children
469 my.plotNode <-
470 function(x1, x2, subtree, type, center, leaflab, dLeaf,
471 nodePar, edgePar, horiz = FALSE)
472 {
473 inner <- !is.leaf(subtree) && x1 != x2
474 yTop <- attr(subtree, "height")
475 bx <- my.plotNodeLimit(x1, x2, subtree, center)
476 xTop <- bx$x
477 usrpar <- par("usr");
478
479 ## handle node specific parameters in "nodePar":
480 hasP <- !is.null(nPar <- attr(subtree, "nodePar"))
481 if(!hasP) nPar <- nodePar
482
483 if(getOption("verbose")) {
484 cat(if(inner)"inner node" else "leaf", ":")
485 if(!is.null(nPar)) { cat(" with node pars\n"); str(nPar) }
486 cat(if(inner)paste(" height", formatC(yTop),"; "),
487 "(x1,x2)= (",formatC(x1,wid=4),",",formatC(x2,wid=4),")",
488 "--> xTop=", formatC(xTop, wid=8),"\n", sep="")
489 }
490
491 Xtract <- function(nam, L, default, indx)
492 rep(if(nam %in% names(L)) L[[nam]] else default,
493 length.out = indx)[indx]
494 asTxt <- function(x) # to allow 'plotmath' labels:
495 if(is.character(x) || is.expression(x) || is.null(x)) x else as.character(x)
496
497 i <- if(inner || hasP) 1 else 2 # only 1 node specific par
498
499 if(!is.null(nPar)) { ## draw this node
500 pch <- Xtract("pch", nPar, default = 1:2, i)
501 cex <- Xtract("cex", nPar, default = c(1,1), i)
502 col <- Xtract("col", nPar, default = par("col"), i)
503 bg <- Xtract("bg", nPar, default = par("bg"), i)
504 points(if (horiz) cbind(yTop, xTop) else cbind(xTop, yTop),
505 pch = pch, bg = bg, col = col, cex = cex)
506 }
507
508 if(leaflab == "textlike")
509 p.col <- Xtract("p.col", nPar, default = "white", i)
510 lab.col <- Xtract("lab.col", nPar, default = par("col"), i)
511 lab.cex <- Xtract("lab.cex", nPar, default = c(1,1), i)
512 lab.font <- Xtract("lab.font", nPar, default = par("font"), i)
513 if (is.leaf(subtree)) {
514 ## label leaf
515 if (leaflab == "perpendicular") { # somewhat like plot.hclust
516 if(horiz) {
517 X <- yTop + dLeaf * lab.cex
518 Y <- xTop; srt <- 0; adj <- c(0, 0.5)
519 }
520 else {
521 Y <- yTop - dLeaf * lab.cex
522 X <- xTop; srt <- 90; adj <- 1
523 }
524 nodeText <- asTxt(attr(subtree,"label"))
525 text(X, Y, nodeText, xpd = TRUE, srt = srt, adj = adj,
526 cex = lab.cex, col = lab.col, font = lab.font)
527 }
528 }
529 else if (inner) {
530 segmentsHV <- function(x0, y0, x1, y1) {
531 if (horiz)
532 segments(y0, x0, y1, x1, col = col, lty = lty, lwd = lwd)
533 else segments(x0, y0, x1, y1, col = col, lty = lty, lwd = lwd)
534 }
535 for (k in 1:length(subtree)) {
536 child <- subtree[[k]]
537 ## draw lines to the children and draw them recursively
538 yBot <- attr(child, "height")
539 if (getOption("verbose")) cat("ch.", k, "@ h=", yBot, "; ")
540 if (is.null(yBot))
541 yBot <- 0
542 xBot <-
543 if (center) mean(bx$limit[k:(k + 1)])
544 else bx$limit[k] + .my.midDend(child)
545
546 hasE <- !is.null(ePar <- attr(child, "edgePar"))
547 if (!hasE)
548 ePar <- edgePar
549 i <- if (!is.leaf(child) || hasE) 1 else 2
550 ## define line attributes for segmentsHV():
551 col <- Xtract("col", ePar, default = par("col"), i)
552 lty <- Xtract("lty", ePar, default = par("lty"), i)
553 lwd <- Xtract("lwd", ePar, default = par("lwd"), i)
554 if (type == "triangle") {
555 segmentsHV(xTop, yTop, xBot, yBot)
556 }
557 else { # rectangle
558 segmentsHV(xTop,yTop, xBot,yTop)# h
559 segmentsHV(xBot,yTop, xBot,yBot)# v
560 }
561 vln <- NULL
562 if (is.leaf(child) && leaflab == "textlike") {
563 nodeText <- asTxt(attr(child,"label"))
564 if(getOption("verbose"))
565 cat('-- with "label"',format(nodeText))
566 hln <- 0.6 * strwidth(nodeText, cex = lab.cex)/2
567 vln <- 1.5 * strheight(nodeText, cex = lab.cex)/2
568 rect(xBot - hln, yBot,
569 xBot + hln, yBot + 2 * vln, col = p.col)
570 text(xBot, yBot + vln, nodeText, xpd = TRUE,
571 cex = lab.cex, col = lab.col, font = lab.font)
572 }
573 if (!is.null(attr(child, "edgetext"))) {
574 edgeText <- asTxt(attr(child, "edgetext"))
575 if(getOption("verbose"))
576 cat('-- with "edgetext"',format(edgeText))
577 if (!is.null(vln)) {
578 mx <-
579 if(type == "triangle")
580 (xTop+ xBot+ ((xTop - xBot)/(yTop - yBot)) * vln)/2
581 else xBot
582 my <- (yTop + yBot + 2 * vln)/2
583 }
584 else {
585 mx <- if(type == "triangle") (xTop + xBot)/2 else xBot
586 my <- (yTop + yBot)/2
587 }
588 ## Both for "triangle" and "rectangle" : Diamond + Text
589
590 p.col <- Xtract("p.col", ePar, default = "white", i)
591 p.border <- Xtract("p.border", ePar, default = par("fg"), i)
592 ## edge label pars: defaults from the segments pars
593 p.lwd <- Xtract("p.lwd", ePar, default = lwd, i)
594 p.lty <- Xtract("p.lty", ePar, default = lty, i)
595 t.col <- Xtract("t.col", ePar, default = col, i)
596 t.cex <- Xtract("t.cex", ePar, default = 1, i)
597 t.font<- Xtract("t.font",ePar, default= par("font"), i)
598 t.shift <- Xtract("t.shift", ePar, default = 0.01, i)
599
600 vlm <- strheight(c(edgeText,"h"), cex = t.cex)/2
601 hlm <- strwidth (c(edgeText,"m"), cex = t.cex)/2
602 hl3 <- c(hlm[1], hlm[1] + hlm[2], hlm[1])
603 #polygon(mx+ c(-hl3, hl3), my + sum(vlm)*c(-1:1,1:-1),
604 # col = p.col, border= p.border, lty = p.lty, lwd = p.lwd)
605 #text(mx, my, edgeText, cex = t.cex, col = t.col, font = t.font)
606 if(horiz) {
607 text(my, mx+t.shift*abs(usrpar[3]-usrpar[4]), edgeText, cex = t.cex, col = t.col, font = t.font)
608 } else {
609 text(mx+t.shift*abs(usrpar[2]-usrpar[1]), my, edgeText, cex = t.cex, col = t.col, font = t.font)
610 }
611 }
612 my.plotNode(bx$limit[k], bx$limit[k + 1], subtree = child,
613 type, center, leaflab, dLeaf, nodePar, edgePar, horiz)
614 }
615 }
616 }
617
618 my.plotNodeLimit <- function(x1, x2, subtree, center)
619 {
620 ## get the left borders limit[k] of all children k=1..K, and
621 ## the handle point `x' for the edge connecting to the parent.
622 inner <- !is.leaf(subtree) && x1 != x2
623 if(inner) {
624 K <- length(subtree)
625 mTop <- .my.memberDend(subtree)
626 limit <- integer(K)
627 xx1 <- x1
628 for(k in 1:K) {
629 m <- .my.memberDend(subtree[[k]])
630 ##if(is.null(m)) m <- 1
631 xx1 <- xx1 + (if(center) (x2-x1) * m/mTop else m)
632 limit[k] <- xx1
633 }
634 limit <- c(x1, limit)
635 } else { ## leaf
636 limit <- c(x1, x2)
637 }
638 mid <- attr(subtree, "midpoint")
639 center <- center || (inner && !is.numeric(mid))
640 x <- if(center) mean(c(x1,x2)) else x1 + (if(inner) mid else 0)
641 list(x = x, limit = limit)
642 }
643
644 .my.memberDend <- function(x) {
645 r <- attr(x,"x.member")
646 if(is.null(r)) {
647 r <- attr(x,"members")
648 if(is.null(r)) r <- 1:1
649 }
650 r
651 }
652
653 .my.midDend <- function(x)
654 if(is.null(mp <- attr(x, "midpoint"))) 0 else mp
655
656
657 ## original Andy Liaw; modified RG, MM :
658 my.heatmap <- function (x, Rowv=NULL, Colv=if(symm)"Rowv" else NULL,
659 distfun = dist, hclustfun = hclust,
660 reorderfun = function(d,w) reorder(d,w),
661 add.expr, symm = FALSE, revC = identical(Colv, "Rowv"),
662 scale = c("row", "column", "none"), na.rm=TRUE,
663 margins = c(5, 5), ColSideColors, RowSideColors,
664 cexRow = 0.2 + 1/log10(nr), cexCol = 0.2 + 1/log10(nc),
665 labRow = NULL, labCol = NULL, main = NULL, xlab = NULL, ylab = NULL,
666 keep.dendro = FALSE,
667 verbose = getOption("verbose"), imageSize=4, imageVSize=imageSize,imageHSize=imageSize,lasCol=2, lasRow=2, respect=F, ...)
668 {
669 scale <- if(symm && missing(scale)) "none" else match.arg(scale)
670 if(length(di <- dim(x)) != 2 || !is.numeric(x))
671 stop("'x' must be a numeric matrix")
672 nr <- di[1]
673 nc <- di[2]
674 if(nr <= 1 || nc <= 1)
675 stop("'x' must have at least 2 rows and 2 columns")
676 if(!is.numeric(margins) || length(margins) != 2)
677 stop("'margins' must be a numeric vector of length 2")
678
679 doRdend <- !identical(Rowv,NA)
680 doCdend <- !identical(Colv,NA)
681 ## by default order by row/col means
682 if(is.null(Rowv)) Rowv <- rowMeans(x, na.rm = na.rm)
683 if(is.null(Colv)) Colv <- colMeans(x, na.rm = na.rm)
684
685 ## get the dendrograms and reordering indices
686
687 if(doRdend) {
688 if(inherits(Rowv, "dendrogram"))
689 ddr <- Rowv
690 else {
691 hcr <- hclustfun(distfun(x))
692 ddr <- as.dendrogram(hcr)
693 if(!is.logical(Rowv) || Rowv)
694 ddr <- reorderfun(ddr, Rowv)
695 }
696 if(nr != length(rowInd <- order.dendrogram(ddr)))
697 stop("row dendrogram ordering gave index of wrong length")
698 }
699 else rowInd <- 1:nr
700
701 if(doCdend) {
702 if(inherits(Colv, "dendrogram"))
703 ddc <- Colv
704 else if(identical(Colv, "Rowv")) {
705 if(nr != nc)
706 stop('Colv = "Rowv" but nrow(x) != ncol(x)')
707 ddc <- ddr
708 }
709 else {
710 hcc <- hclustfun(distfun(if(symm)x else t(x)))
711 ddc <- as.dendrogram(hcc)
712 if(!is.logical(Colv) || Colv)
713 ddc <- reorderfun(ddc, Colv)
714 }
715 if(nc != length(colInd <- order.dendrogram(ddc)))
716 stop("column dendrogram ordering gave index of wrong length")
717 }
718 else colInd <- 1:nc
719
720 ## reorder x
721 x <- x[rowInd, colInd]
722
723 labRow <-
724 if(is.null(labRow))
725 if(is.null(rownames(x))) (1:nr)[rowInd] else rownames(x)
726 else labRow[rowInd]
727 labCol <-
728 if(is.null(labCol))
729 if(is.null(colnames(x))) (1:nc)[colInd] else colnames(x)
730 else labCol[colInd]
731
732 if(scale == "row") {
733 x <- sweep(x, 1, rowMeans(x, na.rm = na.rm))
734 sx <- apply(x, 1, sd, na.rm = na.rm)
735 x <- sweep(x, 1, sx, "/")
736 }
737 else if(scale == "column") {
738 x <- sweep(x, 2, colMeans(x, na.rm = na.rm))
739 sx <- apply(x, 2, sd, na.rm = na.rm)
740 x <- sweep(x, 2, sx, "/")
741 }
742
743 ## Calculate the plot layout
744 lmat <- rbind(c(NA, 3), 2:1)
745 lwid <- c(if(doRdend) 1 else 0.05, imageHSize)
746 lhei <- c((if(doCdend) 1 else 0.05) + if(!is.null(main)) 0.2 else 0, imageVSize)
747 if(!missing(ColSideColors)) { ## add middle row to layout
748 if(!is.character(ColSideColors) || length(ColSideColors) != nc)
749 stop("'ColSideColors' must be a character vector of length ncol(x)")
750 lmat <- rbind(lmat[1,]+1, c(NA,1), lmat[2,]+1)
751 lhei <- c(lhei[1], 0.2, lhei[2])
752 }
753 if(!missing(RowSideColors)) { ## add middle column to layout
754 if(!is.character(RowSideColors) || length(RowSideColors) != nr)
755 stop("'RowSideColors' must be a character vector of length nrow(x)")
756 lmat <- cbind(lmat[,1]+1, c(rep(NA, nrow(lmat)-1), 1), lmat[,2]+1)
757 lwid <- c(lwid[1], 0.2, lwid[2])
758 }
759 lmat[is.na(lmat)] <- 0
760 if(verbose) {
761 cat("layout: widths = ", lwid, ", heights = ", lhei,"; lmat=\n")
762 print(lmat)
763 }
764
765 ## Graphics `output' -----------------------
766
767 op <- par(no.readonly = TRUE)
768 on.exit(par(op))
769 layout(lmat, widths = lwid, heights = lhei, respect = respect)
770 ## draw the side bars
771 if(!missing(RowSideColors)) {
772 par(mar = c(margins[1],0, 0,0.5))
773 image(rbind(1:nr), col = RowSideColors[rowInd], axes = FALSE)
774 }
775 if(!missing(ColSideColors)) {
776 par(mar = c(0.5,0, 0,margins[2]))
777 image(cbind(1:nc), col = ColSideColors[colInd], axes = FALSE)
778 }
779 ## draw the main carpet
780 par(mar = c(margins[1], 0, 0, margins[2]))
781 if(!symm || scale != "none")
782 x <- t(x)
783 if(revC) { # x columns reversed
784 iy <- nr:1
785 ddr <- rev(ddr)
786 x <- x[,iy]
787 } else iy <- 1:nr
788
789 image(1:nc, 1:nr, x, xlim = 0.5+ c(0, nc), ylim = 0.5+ c(0, nr),
790 axes = FALSE, xlab = "", ylab = "", ...)
791 axis(1, 1:nc, labels= labCol, las= lasCol, line= -0.5, tick= 0, cex.axis= cexCol)
792 if(!is.null(xlab)) mtext(xlab, side = 1, line = margins[1] - 1.25)
793 axis(4, iy, labels= labRow, las= lasRow, line= -0.5, tick= 0, cex.axis= cexRow)
794 if(!is.null(ylab)) mtext(ylab, side = 4, line = margins[2] - 1.25,las=lasRow)
795 if (!missing(add.expr))
796 eval(substitute(add.expr))
797
798 ## the two dendrograms :
799 par(mar = c(margins[1], 0, 0, 0))
800 if(doRdend)
801 my.plot.dendrogram(ddr, horiz = TRUE, axes = FALSE, yaxs = "i", leaflab = "none")
802 else frame()
803
804 par(mar = c(0, 0, if(!is.null(main)) 1 else 0, margins[2]))
805 if(doCdend)
806 my.plot.dendrogram(ddc, axes = FALSE, xaxs = "i", leaflab = "none")
807 else if(!is.null(main)) frame()
808
809 ## title
810 if(!is.null(main)) title(main, cex.main = 1.5*op[["cex.main"]])
811
812 invisible(list(rowInd = rowInd, colInd = colInd,
813 Rowv = if(keep.dendro && doRdend) ddr,
814 Colv = if(keep.dendro && doCdend) ddc ))
815 }