|
31
|
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 }
|
