####
# Produce dated tree plot, LTT plot, and MEDUSA (shift of diversification rate)
###

# Produced by Sven BUERKI, October 2025

###
# Load packages
###
# Install packages if needed
# install.packages(c("ape", "phytools"))
# install.packages("strap")  # if not installed yet

library(ape)
library(phytools)
library(strap)
library(geiger)

###
# Dated tree plot
###

data(timescales)  # load the geological timescales data

# 1. Read treePL tree (NEXUS format)
tree_file <- "turr_dated_1_v18.tre"
tree <- read.nexus(tree_file)
# Ladderize tree
tree <- ladderize(tree, right = F)

pdf("Figures/Fig3_Dated_tree_Turraea.pdf")
# Add geological timescale below tree
# phytools geoscalePlot expects an ultrametric tree, make sure the tree is ultrametric:
if (!is.ultrametric(tree)) {
  message("Tree is not ultrametric, forcing ultrametric with chronos()")
  tree <- chronos(tree)
}

# Plot geological timescale with default Phanerozoic scale
# Add margin to bottom for geological timescale axis
par(mar = c(5, 4, 4, 2) + 0.1)
tree_height <- max(node.depth.edgelength(tree))
tree$root.time <- tree_height
strap::geoscalePhylo(tree, units = c("Epoch"), boxes = "Epoch", width = 1.5, show.tip.label = TRUE)
dev.off()

###
# LTT plot
###
pdf("Figures/LTT_plot.pdf")
ltt.plot(tree, log = "y", col = "black", lwd = 2, 
         main = "Lineages-Through-Time Plot",
         ylab = "Number of Lineages (log scale)",
         xlab = "Time (Ma)")
dev.off()

###
# MEDUSA
###
res <- medusa(tree)
summary(res)

# Model.ID Shift.Node Cut.At Model Ln.Lik.part         r epsilon     r.low    r.high
# 1        1        120   node  yule   -96.80832 0.0403947      NA 0.0260406 0.0591962
# 2        2        129   stem  yule   -295.7554 0.1169110      NA 0.0948735 0.1421514
#   One diversification rate shift detected
# 
# Background rate (r₁): 0.04
# 
# Shifted clade rate (r₂): 0.12
# 
# Interpretation: A lineage stemming from node 129 experienced a roughly threefold increase in net diversification rate relative to the rest of the tree.