[Phylobase-commits] r258 - in branches/pdcgsoc: R misc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Aug 6 21:33:06 CEST 2008
Author: pdc
Date: 2008-08-06 21:33:06 +0200 (Wed, 06 Aug 2008)
New Revision: 258
Added:
branches/pdcgsoc/R/treePlot.R
Removed:
branches/pdcgsoc/misc/temp.R
Log:
move tree plotting functions into the R directory so they will be loaded with the package
Added: branches/pdcgsoc/R/treePlot.R
===================================================================
--- branches/pdcgsoc/R/treePlot.R (rev 0)
+++ branches/pdcgsoc/R/treePlot.R 2008-08-06 19:33:06 UTC (rev 258)
@@ -0,0 +1,432 @@
+require(phylobase)
+require(grid)
+require(lattice)
+treePlot <- function(phy,
+ type = c('phylogram', 'cladogram'),
+ show.tip.label = TRUE,
+ show.node.label = FALSE,
+ tip.order = NULL,
+ plot.data = is(phy, 'phylo4d'),
+ rot = 0,
+ tip.plot.fun = 'bubbles',
+ edge.color = 'black', ## TODO colors for branhes and nodes seperately?
+ node.color = 'black',
+ tip.color = 'black',
+ edge.width = 1, ## TODO currently only one width is allowed allow many?
+ ...
+ )
+{
+ width <- height <- 0.9
+ type <- match.arg(type)
+ phy.orig <- phy
+ Nedges <- nrow(phy at edge)
+ Ntips <- length(phy at tip.label)
+ if(is.null(edgeLength(phy)) || type == 'cladogram') {
+ # TODO there should be an abstraction for assigning branch lengths
+ phy at edge.length <- rep(1, nrow(phy at edge))
+ }
+ xxyy <- phyloXXYY(phy, tip.order)
+ phy <- xxyy$phy
+ # TODO this is pointless no? simply returns 1:Ntips
+ tindex <- phy at edge[phy at edge[, 2] <= Ntips, 2]
+ if(type == 'cladogram') {
+ xxyy$xx[phy at edge[, 2] <= Ntips] <- 1
+ }
+
+ # TODO cladogram methods incorrect
+ # TODO abstract, make ultrametric? good algorithms for this?
+ grid.newpage()
+ ## because we may reoder the tip, we need to update the phy objec
+
+ if(!plot.data) {
+ phyplotlayout <- grid.layout(nrow = 1, ncol = 1)
+ # TODO for very long plots, alternative margin setting useful
+ pushViewport(viewport(width = width, height = height,
+ layout = phyplotlayout,
+ name = 'phyplotlayout', angle = -rot))
+ pushViewport(viewport(layout.pos.col = 1, layout.pos.row = 1))
+ tree.plot(xxyy, show.tip.label, show.node.label,
+ edge.color, node.color, tip.color,
+ edge.width, rot)
+ upViewport()
+ upViewport()
+ # TODO should return something useful
+ return(invisible())
+ }
+
+ if(plot.data) {
+ if(is.function(tip.plot.fun)) {
+ # to keep the tip plots square and non-overlapping create a
+ # Ntip * Ntip + 1 grid. The last column is respected which scales
+ # the width of the data col with the height of the window
+ # this might be excessive with a performance penalty
+ datalayout <- grid.layout(nrow = Ntips, ncol = Ntips + 1,
+ # width = unit.c(unit(rep(1, Ntips), rep('null', Ntips)), unit(1/Ntips, 'npc')),
+ respect = matrix(c(rep(0, Ntips * Ntips), rep(0, Ntips - 1), 1), nrow = Ntips)
+ )
+ # TODO this is done multiple times,
+ pushViewport(viewport(width = width, height = height,
+ layout = datalayout,
+ name = 'datalayout', angle = -rot))
+ pushViewport(viewport(layout.pos.col = 1:Ntips))
+ tree.plot(xxyy, show.tip.label, show.node.label,
+ edge.color, node.color, tip.color,
+ edge.width, rot)
+ upViewport()
+
+ pushViewport(viewport(
+ yscale = c(-0.5/Ntips, 1 + 0.5/Ntips),
+ layout.pos.col = Ntips + 1,
+ name = 'data_plots'))
+ # grid.rect(gp = gpar(col = 2))
+ ## TODO should plots float at tips, or only along edge?
+ for(i in xxyy$yy[which(phy at edge[, 2] <= Ntips)]) {
+ pushViewport(viewport(
+ y = i,
+ x = 0.5,
+ height = unit(1/Ntips, 'npc'), # snpc keeps the viewports sq
+ width = unit(1, 'snpc'),
+ name = paste('data_plot', i),
+ just = "center"
+ ))
+ grid.rect()
+ tip.plot.fun()
+ upViewport()
+ }
+ upViewport()
+ upViewport()
+ } else {
+ # use phylobubbles as default
+ dlabwdth <- max(stringWidth(colnames(phy at tip.data)))
+ phyplotlayout <- grid.layout(nrow = 2, ncol = 2,
+ heights = unit.c(unit(1, 'null', NULL), dlabwdth),
+ widths = unit(c(1, 1), c('null', 'null'), list(NULL, NULL))
+ )
+ pushViewport(viewport(width = width, height = height,
+ layout = phyplotlayout,
+ name = 'phyplotlayout', angle = -rot))
+ pushViewport(viewport(layout.pos.row = 1:2, layout.pos.col = 2,
+ height = unit(1, 'npc', NULL) +
+ convertUnit(dlabwdth, 'npc'),
+ default.units = 'native'))
+ phylobubbles(xxyy, ...)
+ upViewport()
+ pushViewport(viewport(layout.pos.col = 1, layout.pos.row = 1))
+ tree.plot(xxyy, show.tip.label, show.node.label,
+ edge.color, node.color, tip.color,
+ edge.width, rot)
+ upViewport()
+ upViewport()
+ }
+ }
+}
+
+tree.plot <- function(xxyy, show.tip.label, show.node.label, edge.color,
+ node.color, tip.color, edge.width, rot)
+{
+ # TODO switch to phylobase abstractions
+ phy <- xxyy$phy
+ Nedges <- nrow(phy at edge)
+ Ntips <- length(phy at tip.label)
+ tindex <- phy at edge[phy at edge[, 2] <= Ntips, 2]
+ eindex <- match(phy at edge[,2], xxyy$phy.orig at edge[,2])
+ segs <- segs(XXYY = xxyy)
+
+ ## TODO check that colors are valid?
+ ## TODO edge colors are required to be in the order of edge matrix
+ if(length(edge.color) != Nedges) {
+ edge.color <- rep(edge.color, length.out = Nedges)
+ }
+ edge.color <- edge.color[eindex]
+
+ ## TODO check that colors are valid?
+ nindex <- sort(eindex[phy at edge[, 2] > Ntips], index.return = TRUE)$ix
+ if(length(node.color) != length(nindex)) {
+ node.color <- rep(node.color, length.out = length(nindex))
+ }
+ node.color <- node.color[nindex]
+
+ if(show.tip.label) {
+ labw <- max(stringWidth(phy at tip.label))
+ treelayout <- grid.layout(nrow = 1, ncol = 2,
+ widths = unit.c(unit(1, 'null', NULL), labw + unit(0.02, 'npc'))
+ )
+ tindex <- phy at edge[phy at edge[, 2] <= Ntips, 2]
+ if(length(tip.color) != Ntips) {
+ tip.color <- rep(tip.color, length.out = Ntips)
+ }
+ } else {
+ treelayout <- grid.layout(nrow = 1, ncol = 1)
+ }
+ # grid.show.layout(treelayout)
+ pushViewport(viewport(
+ x = 0.5, y = 0.5,
+ width = 1, height = 1,
+ # rotataion set here
+ layout = treelayout, name = 'treelayout', angle = -rot))
+ pushViewport(viewport(
+ layout = treelayout, layout.pos.col = 1,
+ name = 'tree'))
+ vseg <- grid.segments( # draws vertical lines
+ x0 = segs$v0x, y0 = segs$v0y,
+ x1 = segs$v1x, y1 = segs$v1y,
+ name = "vert", gp = gpar(col = node.color, lwd = 1))
+ hseg <- grid.segments( # draws horizontal lines
+ x0 = segs$h0x, y0 = segs$h0y,
+ x1 = segs$h1x, y1 = segs$h1y,
+ name = "horz", gp = gpar(col = edge.color, lwd = 1))
+ upViewport()
+ if(show.tip.label) {
+ pushViewport(viewport(
+ layout = treelayout, layout.pos.col = 1,
+ ))
+ labtext <- grid.text(
+ phy at tip.label[tindex],
+ x = xxyy$xx[phy at edge[, 2] %in% tindex] + 0.02,
+ ## TODO yuck!!
+ y = xxyy$yy[phy at edge[, 2] %in% tindex],
+ default.units = 'npc',
+ rot = rot, just = 'left', gp = gpar(col = tip.color[tindex])
+ )
+ upViewport()
+ }
+ # TODO probably want to be able to adjust the location of these guys
+ if(show.node.label) {
+ pushViewport(viewport(
+ layout = treelayout, layout.pos.col = 1,
+ ))
+ labtext <- grid.text(
+ phy at node.label[nindex],
+ x = xxyy$xx[phy at edge[, 2] > Ntips],
+ ## TODO yuck!!
+ y = xxyy$yy[phy at edge[, 2] > Ntips],
+ default.units = 'npc',
+ rot = rot, just = 'left', gp = gpar(col = node.color[nindex])
+ )
+ upViewport()
+ }
+ upViewport()
+ # grobTree(vseg, hseg, labtext)
+}
+
+
+phyloXXYY <- function(phy, tip.order = NULL) {
+ ## initalize the output
+ Nedges <- nrow(phy at edge)
+ phy.orig <- phy
+ Ntips <- length(phy at tip.label)
+ xxyy = list(
+ yy = rep(NA, Nedges),
+ xx = numeric(Nedges),
+ ## record the order that nodes are visited in
+ traverse = NULL)
+ if(is.null(edgeLength(phy))) {
+ # TODO there should be an abstraction for assigning branch lengths
+ stop('Phylogeny has no branch lengths, cannot calculate x coordinates')
+ }
+
+ # TODO tip ordering should be dealt with at a higher level
+ # if(!is.null(tip.order)) {
+ ## TODO do we need to acount for line weight when plotting close to edges?
+ # yy[which(phy at edge[, 2] == tip.order)] <- seq(
+ ## TODO perhaps we want to use match here?
+ ## 0, 1, length.out = Ntips)
+ # } else {
+ ## reoder the phylo and assign even y spacing to the tips
+ phy <- reorder(phy, 'pruningwise')
+ xxyy$yy[phy at edge[, 2] <= Ntips] <- seq(
+ 0, 1, length.out = Ntips
+ )
+ # }
+
+ ## a recurvise preorder traversal
+ ## node -- initalized to be root, is the starting point for the traversal
+ ## phy -- the phylogeny
+ ## xxyy -- the list initalized below that holds the output
+ ## prevx -- the sum of ancestral branch lengths
+ calc.node.xy <- function(node, phy, xxyy, prevx = 0) {
+ ## if node == root node, and there is no root edge set get descendants
+ ## and set index to NULL index is used for indexing output
+ if(any(phy at edge[, 2] == node) == FALSE) {
+ decdex <- which(phy at edge[, 1] == node)
+ index <- NULL
+ ## if root node start at x = 0
+ newx <- 0
+ } else {
+ ## non-root node behavior
+ ## get row in edge matrix corresponding to node, get descendants
+ index <- which(phy at edge[, 2] == node)
+ decdex <- which(phy at edge[, 1] == phy at edge[index, 2])
+ ## non-root node x location
+ newx <- xxyy$xx[index] <- phy at edge.length[index] + prevx
+ ## if the x value is already set we are at a tip and we return
+ if(!is.na(xxyy$yy[index])) { return(xxyy) }
+ }
+ for(i in phy at edge[decdex, 2]) {
+ ## for each decendant call the function again
+ xxyy <- calc.node.xy(i, phy, xxyy, newx)
+ }
+ if(!is.null(index)) {
+ ## set y value by averaging the decendants
+ xxyy$yy[index] <- mean(xxyy$yy[decdex])
+ }
+ ## TODO performance improvement here? rely on above ordering?
+ ## keep track of the nodes and order we visited them
+ xxyy$traverse <- c(xxyy$traverse, phy at edge[decdex, 2])
+ xxyy
+ }
+ ## call function for the first time
+ xxyy <- calc.node.xy(Ntips + 1, phy, xxyy)
+ ## scale the x values
+ xxyy$xx <- xxyy$xx / max(xxyy$xx)
+ # TODO return an index vector instead of a second phy object
+ c(xxyy, phy = list(phy), phy.orig = list(phy.orig))
+}
+
+segs <- function(XXYY) {
+ phy <- XXYY$phy
+ treelen <- rep(NA, nrow(phy at edge) + 1)
+ segs <- list(v0x = treelen, v0y = treelen, v1x = treelen, v1y = treelen,
+ h0x = treelen, h0y = treelen, h1x = treelen, h1y = treelen)
+ troot <- length(phy at tip.label) + 1
+
+ get.coor <- function(node, segs) {
+ if(any(phy at edge[, 2] == node) == FALSE) {
+ decdex <- which(phy at edge[, 1] == node)
+ index <- length(treelen)
+ segs$v0x[index] <- segs$v1x[index] <- 0
+
+ segs$h0y[index] <- segs$h1y[index] <- NA
+ segs$h0x[index] <- segs$h1x[index] <- NA
+ segs$h0x[decdex] <- 0
+ } else {
+ index <- which(phy at edge[, 2] == node)
+ if(!any(phy at edge[, 1] == node)) {
+ return(segs)
+ }
+ decdex <- which(phy at edge[, 1] == phy at edge[index, 2])
+ segs$v0x[index] <- segs$v1x[index] <- XXYY$xx[index]
+ segs$h0x[decdex] <- XXYY$xx[index]
+ }
+ segs$h1x[decdex] <- XXYY$xx[decdex]
+ segs$h0y[decdex] <- segs$h1y[decdex] <- XXYY$yy[decdex]
+
+ segs$v0y[index] <- min(XXYY$yy[decdex])
+ segs$v1y[index] <- max(XXYY$yy[decdex])
+
+ for(i in phy at edge[decdex, 2]) {
+ segs <- get.coor(i, segs)
+ }
+ segs
+ }
+ get.coor(troot, segs)
+}
+
+phylobubbles <- function(XXYY, square = FALSE) {
+ phy <- XXYY$phy
+
+ tys <- XXYY$yy[phy at edge[, 2] <= nTips(phy)]
+
+ traits <- phy at tip.data
+
+ maxr <- ifelse(ncol(traits) > nTips(phy), 1/ncol(traits), 1/nTips(phy))
+
+ tnames <- names(traits)
+ traits <- scale(traits)
+ traits <- apply(traits, 2, function(x) maxr * x / max(abs(x), na.rm = T))
+ names(traits) <- tnames
+
+ if(ncol(traits) == 1) {
+ xpos <- 0.5
+ } else {
+ xpos <- seq(0+maxr, 1-maxr, length.out = ncol(traits))
+ }
+ tys <- tys # * (1 - (2 * maxr)) + maxr
+ xrep <- rep(xpos, each = length(tys))
+ ccol <- ifelse(traits < 0, 'black', 'white')
+ nays <- tys[apply(traits, 1, function(x) any(is.na(x)))]
+ naxs <- xpos[apply(traits, 2, function(x) any(is.na(x)))]
+ traits[is.na(traits)] <- 0
+
+ datalabwidth <- max(stringWidth(colnames(traits)))
+ tiplabwidth <- max(stringWidth(phy at tip.label))
+
+ bublayout <- grid.layout(nrow = 2, ncol = 2,
+ widths = unit.c(unit(1, 'null', NULL), tiplabwidth),
+ heights = unit.c(unit(1, 'null', NULL), datalabwidth))
+ pushViewport(viewport(
+ x = 0.5, y = 0.5,
+ width = 1, height = 1,
+ layout = bublayout, name = 'bublayout'))
+ pushViewport(viewport(
+ name = 'bubble_plots',
+ layout = bublayout,
+ layout.pos.col = 1,
+ layout.pos.row = 1
+ ))
+ grid.segments(x0 = 0, x1 = 1, y0 = tys, y1 = tys, gp = gpar(col = 'grey'))
+ grid.segments(x0 = xpos, x1 = xpos, y0 = 0, y1 = 1, gp = gpar(col = 'grey'))
+ grid.text('x', naxs, nays)
+ if(square) {
+ # to keep the squares square, yet resize nicely use the square npc
+ sqedge <- unit(unlist(traits), 'snpc')
+ grid.rect(x = xrep, y = tys,
+ width = sqedge,
+ height = sqedge,
+ gp=gpar(fill = ccol))
+ } else {
+ grid.circle(xrep, tys, r = unlist(traits), gp = gpar(fill = ccol))
+ }
+ popViewport()
+ pushViewport(viewport(
+ name = 'bubble_tip_labels',
+ layout = bublayout,
+ layout.pos.col = 2,
+ layout.pos.row = 1
+ ))
+ grid.text(phy at tip.label, 0.2, tys, just = 'left')
+ popViewport()
+ pushViewport(viewport(
+ name = 'bubble_data_labels',
+ layout = bublayout,
+ layout.pos.col = 1,
+ layout.pos.row = 2
+ ))
+ grid.text(colnames(traits), xpos, .8, rot = 90, just = 'right')
+ popViewport()
+
+ popViewport()
+}
+
+## How do we translate this info into a plot?
+## Test code
+# out <- phyloXXYY(foo <- as(rcoal(3), 'phylo4'))
+# data(geospiza)
+# foo <- phyloXXYY(geospiza)
+# phylobubbles(foo)
+## TODO true arbitary functions with data from associated data frames
+
+# ff <- function() {grid.lines(1:10/10, runif(10))}
+#
+# p1 <- treePlot(
+# geospiza,
+# # show.tip.label = FALSE,
+# show.node.label = TRUE,
+# edge.color = rainbow(nrow(geospiza at edge)),
+# # node.color = rainbow(nrow(geospiza at edge)),
+# # plot.data = FALSE,
+# # tip.plot.fun = ff,
+# tip.color = c('red', 'black', 'blue'),
+# square = FALSE
+# )
+#
+# treeWpoly <- as(read.tree(text = '((a,b,c),d);'), 'phylo4')
+# print(phyloXXYY(treeWpoly))
+# treePlot(treeWpoly, type = "cladogram")
+
+# n <- 10
+# tree1 <- as(rtree(n), 'phylo4')
+# tree1 at tip.label <- replicate(n, paste(sample(LETTERS, sample(2:20, 1)), collapse = ""))
+# treePlot(tree1)
+
+
Deleted: branches/pdcgsoc/misc/temp.R
===================================================================
--- branches/pdcgsoc/misc/temp.R 2008-08-06 19:26:41 UTC (rev 257)
+++ branches/pdcgsoc/misc/temp.R 2008-08-06 19:33:06 UTC (rev 258)
@@ -1,432 +0,0 @@
-require(phylobase)
-require(grid)
-require(lattice)
-treePlot <- function(phy,
- type = c('phylogram', 'cladogram'),
- show.tip.label = TRUE,
- show.node.label = FALSE,
- tip.order = NULL,
- plot.data = is(phy, 'phylo4d'),
- rot = 0,
- tip.plot.fun = 'bubbles',
- edge.color = 'black', ## TODO colors for branhes and nodes seperately?
- node.color = 'black',
- tip.color = 'black',
- edge.width = 1, ## TODO currently only one width is allowed allow many?
- ...
- )
-{
- width <- height <- 0.9
- type <- match.arg(type)
- phy.orig <- phy
- Nedges <- nrow(phy at edge)
- Ntips <- length(phy at tip.label)
- if(is.null(edgeLength(phy)) || type == 'cladogram') {
- # TODO there should be an abstraction for assigning branch lengths
- phy at edge.length <- rep(1, nrow(phy at edge))
- }
- xxyy <- phyloXXYY(phy, tip.order)
- phy <- xxyy$phy
- # TODO this is pointless no? simply returns 1:Ntips
- tindex <- phy at edge[phy at edge[, 2] <= Ntips, 2]
- if(type == 'cladogram') {
- xxyy$xx[phy at edge[, 2] <= Ntips] <- 1
- }
-
- # TODO cladogram methods incorrect
- # TODO abstract, make ultrametric? good algorithms for this?
- grid.newpage()
- ## because we may reoder the tip, we need to update the phy objec
-
- if(!plot.data) {
- phyplotlayout <- grid.layout(nrow = 1, ncol = 1)
- # TODO for very long plots, alternative margin setting useful
- pushViewport(viewport(width = width, height = height,
- layout = phyplotlayout,
- name = 'phyplotlayout', angle = -rot))
- pushViewport(viewport(layout.pos.col = 1, layout.pos.row = 1))
- tree.plot(xxyy, show.tip.label, show.node.label,
- edge.color, node.color, tip.color,
- edge.width, rot)
- upViewport()
- upViewport()
- # TODO should return something useful
- return(invisible())
- }
-
- if(plot.data) {
- if(is.function(tip.plot.fun)) {
- # to keep the tip plots square and non-overlapping create a
- # Ntip * Ntip + 1 grid. The last column is respected which scales
- # the width of the data col with the height of the window
- # this might be excessive with a performance penalty
- datalayout <- grid.layout(nrow = Ntips, ncol = Ntips + 1,
- # width = unit.c(unit(rep(1, Ntips), rep('null', Ntips)), unit(1/Ntips, 'npc')),
- respect = matrix(c(rep(0, Ntips * Ntips), rep(0, Ntips - 1), 1), nrow = Ntips)
- )
- # TODO this is done multiple times,
- pushViewport(viewport(width = width, height = height,
- layout = datalayout,
- name = 'datalayout', angle = -rot))
- pushViewport(viewport(layout.pos.col = 1:Ntips))
- tree.plot(xxyy, show.tip.label, show.node.label,
- edge.color, node.color, tip.color,
- edge.width, rot)
- upViewport()
-
- pushViewport(viewport(
- yscale = c(-0.5/Ntips, 1 + 0.5/Ntips),
- layout.pos.col = Ntips + 1,
- name = 'data_plots'))
- # grid.rect(gp = gpar(col = 2))
- ## TODO should plots float at tips, or only along edge?
- for(i in xxyy$yy[which(phy at edge[, 2] <= Ntips)]) {
- pushViewport(viewport(
- y = i,
- x = 0.5,
- height = unit(1/Ntips, 'npc'), # snpc keeps the viewports sq
- width = unit(1, 'snpc'),
- name = paste('data_plot', i),
- just = "center"
- ))
- grid.rect()
- tip.plot.fun()
- upViewport()
- }
- upViewport()
- upViewport()
- } else {
- # use phylobubbles as default
- dlabwdth <- max(stringWidth(colnames(phy at tip.data)))
- phyplotlayout <- grid.layout(nrow = 2, ncol = 2,
- heights = unit.c(unit(1, 'null', NULL), dlabwdth),
- widths = unit(c(1, 1), c('null', 'null'), list(NULL, NULL))
- )
- pushViewport(viewport(width = width, height = height,
- layout = phyplotlayout,
- name = 'phyplotlayout', angle = -rot))
- pushViewport(viewport(layout.pos.row = 1:2, layout.pos.col = 2,
- height = unit(1, 'npc', NULL) +
- convertUnit(dlabwdth, 'npc'),
- default.units = 'native'))
- phylobubbles(xxyy, ...)
- upViewport()
- pushViewport(viewport(layout.pos.col = 1, layout.pos.row = 1))
- tree.plot(xxyy, show.tip.label, show.node.label,
- edge.color, node.color, tip.color,
- edge.width, rot)
- upViewport()
- upViewport()
- }
- }
-}
-
-tree.plot <- function(xxyy, show.tip.label, show.node.label, edge.color,
- node.color, tip.color, edge.width, rot)
-{
- # TODO switch to phylobase abstractions
- phy <- xxyy$phy
- Nedges <- nrow(phy at edge)
- Ntips <- length(phy at tip.label)
- tindex <- phy at edge[phy at edge[, 2] <= Ntips, 2]
- eindex <- match(phy at edge[,2], xxyy$phy.orig at edge[,2])
- segs <- segs(XXYY = xxyy)
-
- ## TODO check that colors are valid?
- ## TODO edge colors are required to be in the order of edge matrix
- if(length(edge.color) != Nedges) {
- edge.color <- rep(edge.color, length.out = Nedges)
- }
- edge.color <- edge.color[eindex]
-
- ## TODO check that colors are valid?
- nindex <- sort(eindex[phy at edge[, 2] > Ntips], index.return = TRUE)$ix
- if(length(node.color) != length(nindex)) {
- node.color <- rep(node.color, length.out = length(nindex))
- }
- node.color <- node.color[nindex]
-
- if(show.tip.label) {
- labw <- max(stringWidth(phy at tip.label))
- treelayout <- grid.layout(nrow = 1, ncol = 2,
- widths = unit.c(unit(1, 'null', NULL), labw + unit(0.02, 'npc'))
- )
- tindex <- phy at edge[phy at edge[, 2] <= Ntips, 2]
- if(length(tip.color) != Ntips) {
- tip.color <- rep(tip.color, length.out = Ntips)
- }
- } else {
- treelayout <- grid.layout(nrow = 1, ncol = 1)
- }
- # grid.show.layout(treelayout)
- pushViewport(viewport(
- x = 0.5, y = 0.5,
- width = 1, height = 1,
- # rotataion set here
- layout = treelayout, name = 'treelayout', angle = -rot))
- pushViewport(viewport(
- layout = treelayout, layout.pos.col = 1,
- name = 'tree'))
- vseg <- grid.segments( # draws vertical lines
- x0 = segs$v0x, y0 = segs$v0y,
- x1 = segs$v1x, y1 = segs$v1y,
- name = "vert", gp = gpar(col = node.color, lwd = 1))
- hseg <- grid.segments( # draws horizontal lines
- x0 = segs$h0x, y0 = segs$h0y,
- x1 = segs$h1x, y1 = segs$h1y,
- name = "horz", gp = gpar(col = edge.color, lwd = 1))
- upViewport()
- if(show.tip.label) {
- pushViewport(viewport(
- layout = treelayout, layout.pos.col = 1,
- ))
- labtext <- grid.text(
- phy at tip.label[tindex],
- x = xxyy$xx[phy at edge[, 2] %in% tindex] + 0.02,
- ## TODO yuck!!
- y = xxyy$yy[phy at edge[, 2] %in% tindex],
- default.units = 'npc',
- rot = rot, just = 'left', gp = gpar(col = tip.color[tindex])
- )
- upViewport()
- }
- # TODO probably want to be able to adjust the location of these guys
- if(show.node.label) {
- pushViewport(viewport(
- layout = treelayout, layout.pos.col = 1,
- ))
- labtext <- grid.text(
- phy at node.label[nindex],
- x = xxyy$xx[phy at edge[, 2] > Ntips],
- ## TODO yuck!!
- y = xxyy$yy[phy at edge[, 2] > Ntips],
- default.units = 'npc',
- rot = rot, just = 'left', gp = gpar(col = node.color[nindex])
- )
- upViewport()
- }
- upViewport()
- # grobTree(vseg, hseg, labtext)
-}
-
-
-phyloXXYY <- function(phy, tip.order = NULL) {
- ## initalize the output
- Nedges <- nrow(phy at edge)
- phy.orig <- phy
- Ntips <- length(phy at tip.label)
- xxyy = list(
- yy = rep(NA, Nedges),
- xx = numeric(Nedges),
- ## record the order that nodes are visited in
- traverse = NULL)
- if(is.null(edgeLength(phy))) {
- # TODO there should be an abstraction for assigning branch lengths
- stop('Phylogeny has no branch lengths, cannot calculate x coordinates')
- }
-
- # TODO tip ordering should be dealt with at a higher level
- # if(!is.null(tip.order)) {
- ## TODO do we need to acount for line weight when plotting close to edges?
- # yy[which(phy at edge[, 2] == tip.order)] <- seq(
- ## TODO perhaps we want to use match here?
- ## 0, 1, length.out = Ntips)
- # } else {
- ## reoder the phylo and assign even y spacing to the tips
- phy <- reorder(phy, 'pruningwise')
- xxyy$yy[phy at edge[, 2] <= Ntips] <- seq(
- 0, 1, length.out = Ntips
- )
- # }
-
- ## a recurvise preorder traversal
- ## node -- initalized to be root, is the starting point for the traversal
- ## phy -- the phylogeny
- ## xxyy -- the list initalized below that holds the output
- ## prevx -- the sum of ancestral branch lengths
- calc.node.xy <- function(node, phy, xxyy, prevx = 0) {
- ## if node == root node, and there is no root edge set get descendants
- ## and set index to NULL index is used for indexing output
- if(any(phy at edge[, 2] == node) == FALSE) {
- decdex <- which(phy at edge[, 1] == node)
- index <- NULL
- ## if root node start at x = 0
- newx <- 0
- } else {
- ## non-root node behavior
- ## get row in edge matrix corresponding to node, get descendants
- index <- which(phy at edge[, 2] == node)
- decdex <- which(phy at edge[, 1] == phy at edge[index, 2])
- ## non-root node x location
- newx <- xxyy$xx[index] <- phy at edge.length[index] + prevx
- ## if the x value is already set we are at a tip and we return
- if(!is.na(xxyy$yy[index])) { return(xxyy) }
- }
- for(i in phy at edge[decdex, 2]) {
- ## for each decendant call the function again
- xxyy <- calc.node.xy(i, phy, xxyy, newx)
- }
- if(!is.null(index)) {
- ## set y value by averaging the decendants
- xxyy$yy[index] <- mean(xxyy$yy[decdex])
- }
- ## TODO performance improvement here? rely on above ordering?
- ## keep track of the nodes and order we visited them
- xxyy$traverse <- c(xxyy$traverse, phy at edge[decdex, 2])
- xxyy
- }
- ## call function for the first time
- xxyy <- calc.node.xy(Ntips + 1, phy, xxyy)
- ## scale the x values
- xxyy$xx <- xxyy$xx / max(xxyy$xx)
- # TODO return an index vector instead of a second phy object
- c(xxyy, phy = list(phy), phy.orig = list(phy.orig))
-}
-
-segs <- function(XXYY) {
- phy <- XXYY$phy
- treelen <- rep(NA, nrow(phy at edge) + 1)
- segs <- list(v0x = treelen, v0y = treelen, v1x = treelen, v1y = treelen,
- h0x = treelen, h0y = treelen, h1x = treelen, h1y = treelen)
- troot <- length(phy at tip.label) + 1
-
- get.coor <- function(node, segs) {
- if(any(phy at edge[, 2] == node) == FALSE) {
- decdex <- which(phy at edge[, 1] == node)
- index <- length(treelen)
- segs$v0x[index] <- segs$v1x[index] <- 0
-
- segs$h0y[index] <- segs$h1y[index] <- NA
- segs$h0x[index] <- segs$h1x[index] <- NA
- segs$h0x[decdex] <- 0
- } else {
- index <- which(phy at edge[, 2] == node)
- if(!any(phy at edge[, 1] == node)) {
- return(segs)
- }
- decdex <- which(phy at edge[, 1] == phy at edge[index, 2])
- segs$v0x[index] <- segs$v1x[index] <- XXYY$xx[index]
- segs$h0x[decdex] <- XXYY$xx[index]
- }
- segs$h1x[decdex] <- XXYY$xx[decdex]
- segs$h0y[decdex] <- segs$h1y[decdex] <- XXYY$yy[decdex]
-
- segs$v0y[index] <- min(XXYY$yy[decdex])
- segs$v1y[index] <- max(XXYY$yy[decdex])
-
- for(i in phy at edge[decdex, 2]) {
- segs <- get.coor(i, segs)
- }
- segs
- }
- get.coor(troot, segs)
-}
-
-phylobubbles <- function(XXYY, square = FALSE) {
- phy <- XXYY$phy
-
- tys <- XXYY$yy[phy at edge[, 2] <= nTips(phy)]
-
- traits <- phy at tip.data
-
- maxr <- ifelse(ncol(traits) > nTips(phy), 1/ncol(traits), 1/nTips(phy))
-
- tnames <- names(traits)
- traits <- scale(traits)
- traits <- apply(traits, 2, function(x) maxr * x / max(abs(x), na.rm = T))
- names(traits) <- tnames
-
- if(ncol(traits) == 1) {
- xpos <- 0.5
- } else {
- xpos <- seq(0+maxr, 1-maxr, length.out = ncol(traits))
- }
- tys <- tys # * (1 - (2 * maxr)) + maxr
- xrep <- rep(xpos, each = length(tys))
- ccol <- ifelse(traits < 0, 'black', 'white')
- nays <- tys[apply(traits, 1, function(x) any(is.na(x)))]
- naxs <- xpos[apply(traits, 2, function(x) any(is.na(x)))]
- traits[is.na(traits)] <- 0
-
- datalabwidth <- max(stringWidth(colnames(traits)))
- tiplabwidth <- max(stringWidth(phy at tip.label))
-
- bublayout <- grid.layout(nrow = 2, ncol = 2,
- widths = unit.c(unit(1, 'null', NULL), tiplabwidth),
- heights = unit.c(unit(1, 'null', NULL), datalabwidth))
- pushViewport(viewport(
- x = 0.5, y = 0.5,
- width = 1, height = 1,
- layout = bublayout, name = 'bublayout'))
- pushViewport(viewport(
- name = 'bubble_plots',
- layout = bublayout,
- layout.pos.col = 1,
- layout.pos.row = 1
- ))
- grid.segments(x0 = 0, x1 = 1, y0 = tys, y1 = tys, gp = gpar(col = 'grey'))
- grid.segments(x0 = xpos, x1 = xpos, y0 = 0, y1 = 1, gp = gpar(col = 'grey'))
- grid.text('x', naxs, nays)
- if(square) {
- # to keep the squares square, yet resize nicely use the square npc
- sqedge <- unit(unlist(traits), 'snpc')
- grid.rect(x = xrep, y = tys,
- width = sqedge,
- height = sqedge,
- gp=gpar(fill = ccol))
- } else {
- grid.circle(xrep, tys, r = unlist(traits), gp = gpar(fill = ccol))
- }
- popViewport()
- pushViewport(viewport(
- name = 'bubble_tip_labels',
- layout = bublayout,
- layout.pos.col = 2,
- layout.pos.row = 1
- ))
- grid.text(phy at tip.label, 0.2, tys, just = 'left')
- popViewport()
- pushViewport(viewport(
- name = 'bubble_data_labels',
- layout = bublayout,
- layout.pos.col = 1,
- layout.pos.row = 2
- ))
- grid.text(colnames(traits), xpos, .8, rot = 90, just = 'right')
- popViewport()
-
- popViewport()
-}
-
-## How do we translate this info into a plot?
-## Test code
-# out <- phyloXXYY(foo <- as(rcoal(3), 'phylo4'))
-data(geospiza)
-# foo <- phyloXXYY(geospiza)
-# phylobubbles(foo)
-## TODO true arbitary functions with data from associated data frames
-
-ff <- function() {grid.lines(1:10/10, runif(10))}
-
-p1 <- treePlot(
- geospiza,
- # show.tip.label = FALSE,
- show.node.label = TRUE,
- edge.color = rainbow(nrow(geospiza at edge)),
- # node.color = rainbow(nrow(geospiza at edge)),
- # plot.data = FALSE,
- # tip.plot.fun = ff,
- tip.color = c('red', 'black', 'blue'),
- square = FALSE
-)
-
-treeWpoly <- as(read.tree(text = '((a,b,c),d);'), 'phylo4')
-# print(phyloXXYY(treeWpoly))
-# treePlot(treeWpoly, type = "cladogram")
-
-# n <- 10
-# tree1 <- as(rtree(n), 'phylo4')
-# tree1 at tip.label <- replicate(n, paste(sample(LETTERS, sample(2:20, 1)), collapse = ""))
-# treePlot(tree1)
-
-
More information about the Phylobase-commits
mailing list