How can I plot a tree (and squirrels) in R?

问题:

Here is my tree:

tree = data.frame(branchID = c(1,11,12,111,112,1121,1122), length = c(32, 21, 19, 5, 12, 6, 2))

> tree
  branchID length
1        1     32
2       11     21
3       12     19
4      111      5
5      112     12
6     1121      6
7     1122      2

This tree is in 2D and is made of branches. Each branch has an ID. 1 is the trunk. Then the trunk bifurcate into two branches, 11 on the left and 12 on the right. 11 bifurcates as well in the branches called 111 (going toward the left) and 112 (going toward the right). etc.. Each branch has a certain length.

On this tree there are squirrels:

squirrels = data.frame(branchID = c(1,11,1121,11,111), PositionOnBranch = c(23, 12, 4, 2, 1), name=c("FluffyTail", "Ginger", "NutCracker", "SuperSquirrel", "ChipnDale"))

> squirrels
  branchID PositionOnBranch          name
1        1               23    FluffyTail
2       11               12        Ginger
3     1121                4    NutCracker
4       11                2 SuperSquirrel
5      111                1     ChipnDale

Each squirrel is found on a specific branch. For example the FluffyTail is on the trunk at position 23 (the total length of the trunk being 32). ChipnDale is on the branch 111 at position 1 (the total length of the branch 111 is 5). The position is taken relatively to the lower extremity of the branch.

How can I plot my tree and my squirrels?

回答1:

I put a bit more thought/time into this, and have packaged up some horticultural functions in package trees, here.

With trees, you can:

• generate a random tree design (a random seed, so to speak) with seed();
• sow the seed to give rise to a magnificent tree with germinate();
• add randomly-located leaves (or squirrels) with foliate();
• add squirrels (for example) to specified locations with squirrels(); and
• prune() the tree.

---

# Install the package and set the RNG state
devtools::install_github('johnbaums/trees')
set.seed(1)

Let's fertilise a seed and grow a tree

# Create a tree seed    
s <- seed(70, 10, min.branch.length=0, max.branch.length=4,
          min.trunk.height=5, max.trunk.height=8)

head(s, 10)

#       branch    length
# 1          0 6.3039785
# 2          L 2.8500587
# 3         LL 1.5999775
# 4        LLL 1.3014086
# 5       LLLL 3.0283486
# 6      LLLLL 0.8107690
# 7     LLLLLR 2.8444849
# 8    LLLLLRL 0.4867677
# 9   LLLLLRLR 0.9819541
# 10 LLLLLRLRR 0.5732175

# Germinate the seed
g <- germinate(s, col='peachpuff4')

And add some leaves

leafygreens <- colorRampPalette(paste0('darkolivegreen', c('', 1:4)))(100)
foliate(g, 5000, 4, pch=24:25, col=NA, cex=1.5, bg=paste0(leafygreens, '30'))

Or some squirrels

plot(g, col='peachpuff4')
squirrels(g, 
          branches=c("LLLLRRRL", "LRLRR", "LRRLRLLL", "LRRRLL", "RLLLLLR", 
                     "RLLRL", "RLLRRLRR", "RRRLLRL", "RRRLLRR", "RRRRLR"),
          pos=c(0.22, 0.77, 0.16, 0.12, 0.71, 0.23, 0.18, 0.61, 0.8, 2.71),
          pch=20, cex=2.5)

Plotting @Remi.b's tree and squirrels

g <- germinate(list(trunk.height=32, 
                   branches=c(1, 2, 11, 12, 121, 122),
                   lengths=c(21, 19, 5, 12, 6, 2)), 
              left='1', right='2', angle=40)

xy <- squirrels(g, c(0, 1, 121, 1, 11), pos=c(23, 12, 4, 2, 1), 
               left='1', right='2', pch=21, bg='white', cex=3, lwd=2)
text(xy$x, xy$y, labels=seq_len(nrow(xy)), font=2)
legend('bottomleft', bty='n',
      legend=paste(seq_len(nrow(xy)), 
                   c('FluffyTail', 'Ginger', 'NutCracker', 'SuperSquirrel', 
                     'ChipnDale'), sep='. '))

---

EDIT:

Following @baptiste's hot tip about @ScottChamberlain's rphylopic package, it's time to upgrade those dots to squirrels (though they may resemble coffee beans).

library(rphylopic)
s <- seed(50, 10, min.branch.length=0, max.branch.length=5,
          min.trunk.height=5, max.trunk.height=8)
g <- germinate(s, trunk.width=15, col='peachpuff4')
leafygreens <- colorRampPalette(paste0('darkolivegreen', c('', 1:4)))(100)
foliate(g, 2000, 4, pch=24:25, col=NA, cex=1.2, bg=paste0(leafygreens, '50'))
xy <- foliate(g, 2, 2, 4, xy=TRUE, plot=FALSE)

# snazzy drop shadow
add_phylopic_base(
    image_data("5ebe5f2c-2407-4245-a8fe-397466bb06da", size = "64")[[1]], 
    1, xy$x, xy$y, ysize = 2.3, col='black')
add_phylopic_base(
    image_data("5ebe5f2c-2407-4245-a8fe-397466bb06da", size = "64")[[1]], 
    1, xy$x, xy$y, ysize = 2, col='darkorange3')

回答2:

I probably over-thought this, but... squirrels.

get.coords <- function(a, d, x0, y0) {
  a <- ifelse(a <= 90, 90 - a, 450 - a)
  data.frame(x = x0 + d * cos(a / 180 * pi), 
             y = y0+ d * sin(a / 180 * pi))
}


tree$angle <- sapply(gsub(2, '+45', gsub(1, '-45', tree$branchID)), 
                     function(x) eval(parse(text=x)))
tree$tipy <- tree$tipx <- tree$basey <- tree$basex <- NA

for(i in seq_len(nrow(tree))) {
  if(tree$branchID[i] == 0) {
    tree$basex[i] <- tree$basey[i] <- tree$tipx[i] <- 0
    tree$tipy[i] <- tree$length[i]
    next
  } else if(tree$branchID[i] %in% 1:2) {
    parent <- 0
  } else {
    parent <- substr(tree$branchID[i], 1, nchar(tree$branchID[i])-1)
  }
  tree$basex[i] <- tree$tipx[which(tree$branchID==parent)]
  tree$basey[i] <- tree$tipy[which(tree$branchID==parent)]
  tip <- get.coords(tree$angle[i], tree$length[i], tree$basex[i], tree$basey[i])
  tree$tipx[i] <- tip[, 1]
  tree$tipy[i] <- tip[, 2]
}  

squirrels$nesty <- squirrels$nestx <- NA
for (i in seq_len(nrow(squirrels))) {
  b <- tree[tree$branchID == squirrels$branchID[i], ]
  nest <- get.coords(b$angle, squirrels$PositionOnBranch[i], b$basex, b$basey)
  squirrels$nestx[i] <- nest[1]
  squirrels$nesty[i] <- nest[2]
}

And now we plot.

plot.new()
plot.window(xlim=range(tree$basex, tree$tipx), 
            ylim=range(tree$basey, tree$tipy), asp=1)
with(tree, segments(basex, basey, tipx, tipy, lwd=pmax(10/nchar(branchID), 1)))
points(squirrels[, c('nestx', 'nesty')], pch=21, cex=3, bg='white', lwd=2)
text(squirrels[, c('nestx', 'nesty')], labels=seq_len(nrow(squirrels)), font=2)
legend('bottomleft', legend=paste(seq_len(nrow(squirrels)), squirrels$name), bty='n')

And for kicks we will simulate a bigger tree (and put some apples on it like in Farmville):

twigs <- replicate(50, paste(rbinom(5, 1, 0.5) + 1, collapse=''))
branches <- sort(unique(c(sapply(twigs, function(x) sapply(seq_len(nchar(x)), function(y) substr(x, 1, y))))))
tree <- data.frame(branchID=c(0, branches), length=c(30, sample(10, length(branches), TRUE)), 
                   stringsAsFactors=FALSE)


tree$angle <- sapply(gsub(2, '+45', gsub(1, '-45', tree$branchID)), 
                     function(x) eval(parse(text=x)))
tree$tipy <- tree$tipx <- tree$basey <- tree$basex <- NA

for(i in seq_len(nrow(tree))) {
  if(tree$branchID[i] == 0) {
    tree$basex[i] <- tree$basey[i] <- tree$tipx[i] <- 0
    tree$tipy[i] <- tree$length[i]
    next
  } else if(tree$branchID[i] %in% 1:2) {
    parent <- 0
  } else {
    parent <- substr(tree$branchID[i], 1, nchar(tree$branchID[i])-1)
  }
  tree$basex[i] <- tree$tipx[which(tree$branchID==parent)]
  tree$basey[i] <- tree$tipy[which(tree$branchID==parent)]
  tip <- get.coords(tree$angle[i], tree$length[i], tree$basex[i], tree$basey[i])
  tree$tipx[i] <- tip[, 1]
  tree$tipy[i] <- tip[, 2]
}  

plot.new()
plot.window(xlim=range(tree$basex, tree$tipx), 
            ylim=range(tree$basey, tree$tipy), asp=1)
par(mar=c(0, 0, 0, 0))
with(tree, segments(basex, basey, tipx, tipy, lwd=pmax(20/nchar(branchID), 1)))

apple_branches <- sample(branches, 10)
sapply(apple_branches, function(x) {
  b <- tree[tree$branchID == x, ]
  apples <- get.coords(b$angle, runif(sample(2, 1), 0, b$length), b$basex, b$basey)
  points(apples, pch=20, col='tomato2', cex=2)
})

回答3:

Well, you could convert your data to define a "tree" as defined by the ape package. Here's a function that can convert your data.frame to the correct format.

library(ape)

to.tree <- function(dd) {
    dd$parent <- dd$branchID %/% 10

    root <- subset(dd, parent==0)
    dd <- subset(dd, parent!=0)

    ids <- unique(c(dd$parent, dd$branchID))
    tip <- !(ids %in% dd$parent)
    lvl <- ids[order(!tip, ids)]
    edg <- sapply(dd[,c("parent","branchID")], 
        function(x) as.numeric(factor(x, levels=lvl)))

    x<-list(
        edge=edg,
        edge.length=dd$length,
        tip.label=head(lvl, sum(tip)),
        node.label=tail(lvl, length(tip)-sum(tip)),
        Nnode = length(tip)-sum(tip),
        root.edge=root$length[1]
    )
    class(x)<-"phylo"
    reorder(x)    
}

Then we can plot it somewhat easily

xx <- to.tree(tree)
plot(xx, show.node.label=TRUE, root.edge=TRUE)

Now, if we want to add the squirrel information, we need to know where each branch is located. I'm going to borrow getphylo_x and getphylo_y from this answer. Then I can run

sx<-Vectorize(getphylo_x, "node")(xx, as.character(squirrels$branchID)) -
    tree$length[match(squirrels$branchID, tree$branchID)] +
    squirrels$PositionOnBranch
sy<-Vectorize(getphylo_y, "node")(xx, as.character(squirrels$branchID))

points(sx,sy)
text(sx,sy, squirrels$name, pos=3)

to add the squirrel information to the plot. The final result is

It's not perfect but it's not a bad start.

回答4:

The reshaping of this might take a while, but this is broadly possible. E.g., rejigging your data representation so it looks like:

library(igraph)
dat <- read.table(text="1 1n2
1n2 1.1
1n2 1.2
1.1 1.1.1
1.1 1.1.2
1.1.2 1.1.2.1
1.1.2 1.1.2.2",header=FALSE)

g <- graph.data.frame(dat)
tkplot(g)

And manually moving the tree parts around in tkplot, you can get:

Doing this automatically is a whole different story admittedly.

回答5:

A version that supports trees with more than two branches. A bit work is required to convert to a data.tree structure, and to add the squirrels to it. But once you're there, the plotting is straight forward.

df <- data.frame(branchID = c(1,11,12,13, 14, 111,112,1121,1122), length = c(32, 21, 12, 8, 19, 5, 12, 6, 2))
squirrels <- data.frame(branchID = c(1,11,1121,11,111), PositionOnBranch = c(23, 12, 4, 2, 1), squirrel=c("FluffyTail", "Ginger", "NutCracker", "SuperSquirrel", "ChipnDale"), stringsAsFactors = FALSE)

library(magrittr)

#derive pathString from branchID, so we can convert it to data.tree structure
df$branchID %>%
  as.character %>%
  sapply(function(x) strsplit(x, split = "")) %>%
  sapply(function(x) paste(x, collapse = "/")) ->
  df$pathString

df$type <- "branch"

library(data.tree)

tree <- FromDataFrameTable(df)

#climb, little squirrels!
for (i in 1:nrow(squirrels)) {
  squirrels[i, 'branchID'] %>%
    as.character %>%
    strsplit(split = "") %>%
    extract2(1) %>%
    extract(-1) -> path
  if (length(path) > 0) branch <- tree$Climb(path)
  else branch <- tree
  #actually, we add the squirrels as branches to our tree
  #What a symbiotic coexistence!
  #advantage: Our SetCoordinates can be re-used as is
  #disadvantage: may be confusing, and it requires us
  #to do some filtering later
  branch$AddChild(squirrels[i, 'squirrel'],
                 length = squirrels[i, 'PositionOnBranch'],
                 type = "squirrel")
}



SetCoordinates <- function(node, branch) {
  if (branch$isRoot) {
    node$x0 <- 0
    node$y0 <- 0
  } else {
    node$x0 <- branch$parent$x1
    node$y0 <- branch$parent$y1
  }

  #let's hope our squirrels didn't flunk in trigonometry ;-)
  angle <- branch$position / (sum(Get(branch$siblings, "type") == "branch") + 2)
  x <- - node$length * cospi(angle)
  y <- sqrt(node$length^2 - x^2)
  node$x1 <- node$x0 + x
  node$y1 <- node$y0 + y
}

#let it grow!
tree$Do(function(node) {
          SetCoordinates(node, node)
          node$lwd <- 10 * (node$root$height - node$level + 1) / node$root$height
        }, filterFun = function(node) node$type == "branch")
tree$Do(function(node) SetCoordinates(node, node$parent), filterFun = function(node) node$type == "squirrel")

Looking at the data:

print(tree, "type", "length", "x0", "y0", "x1", "y1")

This prints like so:

                    levelName     type length        x0       y0         x1       y1
1  1                            branch     32   0.00000  0.00000   0.000000 32.00000
2   ¦--1                        branch     21   0.00000 32.00000 -16.989357 44.34349
3   ¦   ¦--1                    branch      5 -16.98936 44.34349 -19.489357 48.67362
4   ¦   ¦   °--ChipnDale      squirrel      1 -16.98936 44.34349 -17.489357 45.20952
5   ¦   ¦--2                    branch     12 -16.98936 44.34349 -10.989357 54.73580
6   ¦   ¦   ¦--1                branch      6 -10.98936 54.73580 -13.989357 59.93195
7   ¦   ¦   ¦   °--NutCracker squirrel      4 -10.98936 54.73580 -12.989357 58.19990
8   ¦   ¦   °--2                branch      2 -10.98936 54.73580  -9.989357 56.46785
9   ¦   ¦--Ginger             squirrel     12   0.00000 32.00000  -9.708204 39.05342
10  ¦   °--SuperSquirrel      squirrel      2   0.00000 32.00000  -1.618034 33.17557
11  ¦--2                        branch     12   0.00000 32.00000  -3.708204 43.41268
12  ¦--3                        branch      8   0.00000 32.00000   2.472136 39.60845
13  ¦--4                        branch     19   0.00000 32.00000  15.371323 43.16792
14  °--FluffyTail             squirrel     23   0.00000  0.00000   0.000000 23.00000

Once we're here, plotting is also easy:

plot(c(min(tree$Get("x0")), max(tree$Get("x1"))),
     c(min(tree$Get("y0")), max(tree$Get("y1"))),
     type='n', asp=1, axes=FALSE, xlab='', ylab='')

tree$Do(function(node) segments(node$x0, node$y0, node$x1, node$y1, lwd = node$lwd),
        filterFun = function(node) node$type == "branch")

tree$Do(function(node) {
          points(node$x1, node$y1, lwd = 8, col = "saddlebrown")
          text(node$x1, node$y1, labels = node$name, pos = 2, cex = 0.7)
        },
        filterFun = function(node) node$type == "squirrel")