#download and install both "strap" and "paleotree"
install.packages("strap")
install.packages("paleotree")
library (strap)
library (paleotree)

#set working directory
setwd("D:/DSR/Fossil Thylacinidae/Phylogeny")

#load FAD/LAD date ranges in a *.txt file, following format:
  # FAD	LAD	
  # Djarthia_murgonensis	54.7	54.6
  # Perameles_nasuta	0.0	0.0
  # Antechinus_flavipes	0.0	0.0
  # Badjcinus_turnbulli	27.8	23.0
  # Barinya_wangala	23.03	11.62
  # Dasycercus_cristicauda	0.0	0.0
  # Dasyuroides_achilpatna	2.9	2.5
  # Dasyurus_dunmalli	4.2	3.6
  # Dasyurus_hallucatus	0.0	0.0
  # Dasyurus_maculatus	0.0	0.0
  # Muribacinus_gadiyuli	15.1	12.9
  # Mutpuracinus_archibaldi	15.1	12.9
  # Ngamalacinus_timmulvaneyi	18.5	16.2
  # Nimbacinus_dicksoni	17.0	12.0
  # Phascogale_tapoatafa	0.0	0.0
  # Sminthopsis_floravillensis	2.58	0.1
  # Sminthopsis_murina	0.0	0.0
  # Thylacinus_cynocephalus	0.0	0.0
  # Thylacinus_macknessi	18.5	16.2
  # Thylacinus_megiriani	7.5	4.5
  # Thylacinus_potens	8.5	5.5
  # Thylacinus_yorkellus	5.3	3.6
  # Tyarrpecinus_rothi	8.5	5.5
  # Wabulacinus_ridei	18.5	17.0
Thylacinidae.ages <- read.table("DatesTotal.txt", header = T)

#load tree in Newick format - NB reconstruct this tree as one line, no spaces
#(Djarthia_murgonensis,(Perameles_nasuta,((Mutpuracinus_archibaldi,Barinya_wangala,
#(Sminthopsis_murina,Dasyuroides_achilpatna,(Antechinus_flavipes,Phascogale_tapoatafa),
#(Sminthopsis_floravillensis,Dasyurus_dunmalli,(Dasycercus_cristicauda,
#(Dasyurus_hallucatus,Dasyurus_maculatus))))),(Badjcinus_turnbulli,
#((Nimbacinus_dicksoni,(Muribacinus_gadiyuli,Ngamalacinus_timmulvaneyi)),
#(Tyarrpecinus_rothi,(Wabulacinus_ridei,(Thylacinus_macknessi,(Thylacinus_potens,
#(Thylacinus_yorkellus,Thylacinus_cynocephalus,Thylacinus_megiriani))))))))));
Thylacinidae.tree <- read.tree("Thyla.tre")

#time scale tree using the "equal" method (modified from Brusatte et al 2008, rlength = 1 to 
#eliminate zero-length branches, add.terminal to include stratigraphic range at tips
Thylacinidae_calibrate.tree <- DatePhylo(Thylacinidae.tree, 
  Thylacinidae.ages, 
  rlen = 1, 
  method = "equal", 
  add.terminal = TRUE)

#plot the time-scaled tree against the chronostratigraphic scale
geoscalePhylo(ladderize(Thylacinidae_calibrate.tree, right=FALSE), 
  Thylacinidae.ages, 
  width = 3, 
  cex.age = 0.75, 
  cex.ts = 0.5, 
  cex.tip = 1)

#node ages can be found by using the dateNodes function in the package paleotree
#remember to use your time-scaled tree
#first values (1-11 here) are LAD of the tips, latter values (here 12-20) 
#are the nodes starting at the root and here 
#ending at the cynocephalus + megiriani + yorkellus polytomy
dateNodes(Thylacinidae_calibrate.tree)