Add a fossil tip to a tree

Question

I am having a hard time trying to add several fossil tips to a 100 phylogenetic time trees.

Here is a drawing of what I want to do:

I would like to do it for my 100 trees contained in a single nexus file.

Thank you in advance.

EDIT:

Okay, I figured out how to do it. Here is the code for the example:

library(phytools)

tree <- read.tree(text = '(A:20,(B1:5,B2:5):15);')
plot(tree)

nodelabels()
edgelabels()

case1 <- bind.tip(tree, tip.label="Fossil tip 1", where=which(tree$tip.label=="A"), position = 10, edge.length=10)
plot(case1)

case2 <- bind.tip(tree, tip.label="Fossil tip 2", where=5, position=0.5*tree$edge.length[2], edge.length=0.5*tree$edge.length[2]+tree$edge.length[3])
plot(case2)

case3 <- bind.tip(tree, tip.label="Fossil tip 3", where=5, position=0.5*tree$edge.length[2], edge.length=(0.5*tree$edge.length[2]+tree$edge.length[3])-10)
plot(case3)

And here is the code fo my data:

library(phytools)

mammif <- read.nexus(file.choose()) #100 phyogenetic time trees
tree <- mammif$tree_7736 #for the first tree
plot(tree)

nodelabels()
edgelabels()

case1 <- bind.tip(tree, tip.label="Fossil tip 1", where=which(tree$tip.label=="Tupaia_belangeri"), position = 30, edge.length=30)
plot(case1)

case2 <- bind.tip(tree, tip.label="Fossil tip 2", where=36, position=0.5*tree$edge.length[21], edge.length=0.5*tree$edge.length[21]+tree$edge.length[22])
plot(case2)

case3 <- bind.tip(tree, tip.label="Fossil tip 3", where=33, position=0.5*tree$edge.length[13], edge.length=(0.5*tree$edge.length[13]+tree$edge.length[14])-10)
plot(case3)

The code works just fine for the first tree. I would like to loop it for the rest of the trees but the thing is that the nodes and the edges can change between the trees. I think it would be easier to designate the taxa of a clade rather than a node and the edge leading to a taxon rather than an edge number. Any idea?

Answer 1

Here is a general approach to find the most recent common ancestor (MRCA) of some tips and bind a (fossil) branch to that node. I am going to show this with some random trees, but write back if you have trouble implementing with your real data.

library(ape)
library(phytools)

First we generate some random ultrametric trees.

set.seed(1)
random_trees <-
  replicate(n = 4,
            expr = rcoal(n = 5, rooted = TRUE),
            simplify = FALSE)

# plot
par(mfrow=c(2,2))
lapply(random_trees, plot.phylo)

Then we add a root edge with some arbitrary length. This is required for phytools::bind.tip (which is a wrapper around ape::bind.tree). Depending on your real trees, you might not need to worry about the root branch length.

random_trees <-
  lapply(seq_along(random_trees), function(i) {
    random_trees[[i]]$root.edge <- 0.2;
    return(random_trees[[i]])})

lapply(random_trees, "[[", "root.edge")

Then, lets asumme we want to add a fossil node on the branch representing the MRCA of tips t1 and t2. We can cycle over the trees to find the wanted node (MRCA of t1 and t2). The topologies here can be different, all we want is that target node.

target_taxa <- paste0("t", 1:2)

target_mrcas <-
  lapply(random_trees, function(x)
    getMRCA(x, tip = target_taxa))

In a final pass, we cycle over the trees and bind our fossil to the desired branch

trees_with_fossil <-
  lapply(seq_along(random_trees), function(i)
    phytools::bind.tip(
      tree = random_trees[[i]],
      tip.label = "fossil",
      edge.length = 0.1,
      where = target_mrcas[[i]],
      position = 0.02
    ))

Plot:

lapply(trees_with_fossil, plot.phylo)