# Time Calibrated Tree Code
install.packages("paleotree")
install.packages(c("ape","phangorn"))
install.packages("devtools")
devtools::install_github("santiagosnchez/rBt")


library(diversitree)
library(phytools)
library(paleotree)
library(ggtree)
library(ape)
library(phangorn)
library(rBt)

library(geiger) #A package for "macroevolutionary simulation and estimating parameters related to diversification"

library(phangorn) 
library(convevol)
library(ellipse)
library(cluster)


non_calibrated_tree_ML <- read.tree("AllGenes.fas.treefile")
ladderize(non_calibrated_tree_ML)
plot(non_calibrated_tree_ML, cex = 0.5)
tiplabels(cex = 0.5)

rootposition <- 0.5 * non_calibrated_tree_ML$edge.length[which(non_calibrated_tree_ML$edge[ ,2] == 70)]

non_calibrated_tree_ML.rooted <- reroot(non_calibrated_tree_ML, 126, position = rootposition) 

plot(non_calibrated_tree_ML.rooted, cex = 0.5)
nodelabels(cex = 0.3, node= 313, bg = NULL)

#*Formatting your node age data (advanced!)
  
#  getMRCA(tree,c(“taxonA”,”taxonB”))

#If you want to double-check nodes, the opposite of getMRCA is getDescendants(tree, node number)

#calibration node 1 - Stem lineage of Balistidae, dating the MRCA of Abalistes and Cantherhines.
getMRCA(non_calibrated_tree_ML.rooted,c("Abalistes_stellatus","Cantherhines_pullus"))
# node 229
# Absolute age offset = 37.2 Ma 
# 95% soft upper bound = 125 Ma 
# Log(Mean) = 2.907
# Log(StDev) = 0.8
calib1 <- getDescendants(non_calibrated_tree_ML.rooted, 229)

#calibration node 2 - Stem lineage Diodontidae, dating the MRCA of Diodontidae (Diodon and Chilomycterus) and Tetraodon.
getMRCA(non_calibrated_tree_ML.rooted,c("Diodon_holocanthus","Chilomycterus_schoepfii", "Tetraodon_miurus"))
# node 449
# Absolute age offset = 50 Ma 
# 95% soft upper bound = 125 Ma 
# Log(Mean) = 2.75 
# Log(StDev) = 0.8
calib2 <- getDescendants(non_calibrated_tree_ML.rooted, 449)



##calibration node 3 - Stem lineage Siganidae, dating the MRCA of Siganus and Scatophagidae (Scatophagus and Selenotoca).
getMRCA(non_calibrated_tree_ML.rooted,c("Siganus_spinus","Scatophagus_argus", "Selenotoca_multifasciata"))
# node 245
# Absolute age offset = 55.8 Ma 
# 95% soft upper bound = 125 Ma 
# Log(Mean) = 2.669
# Log(StDev) = 0.8
calib3 <- getDescendants(non_calibrated_tree_ML.rooted, 245)

#calibration node 4 - Stem lineage Echeneidae, dating the MRCA of Echeneis, Coryphaena, and Rachycentron.
getMRCA(non_calibrated_tree_ML.rooted,c("Echeneis_naucrates","Coryphaena_hippurus", "Rachycentron_canadum"))
# node 298
# Absolute age offset = 30.1 Ma 
# 95% soft upper bound = 125 Ma 
# Log(Mean) = 2.985
# Log(StDev) = 0.8
calib4 <- getDescendants(non_calibrated_tree_ML.rooted, 298)



#calibration node 5 - Stem lineage Carangidae, dating the MRCA of Carangidae (Caranx, Seriola, and Trachinotus), Echeneis, Coryphaena, and Rachycentron.
getMRCA(non_calibrated_tree_ML.rooted,c("Caranx_crysos","Seriola_dumerili", "Trachinotus_carolinus", "Echeneis_naucrates","Coryphaena_hippurus", "Rachycentron_canadum"))
# node 295
# Absolute age offset = 55.8 Ma 
# 95% soft upper bound = 125 Ma 
# Log(Mean) = 2.669
# Log(StDev) = 0.8
calib5 <- getDescendants(non_calibrated_tree_ML.rooted, 295)



#calibration node 6 - Stem lineage Zenopsis, dating the MRCA of Zenopsis and Zeus.
getMRCA(non_calibrated_tree_ML.rooted,c("Zenopsis_conchifer","Zeus_faber"))
# node 331
# Absolute age offset = 32 Ma 
# 95% soft upper bound = 70.6 Ma 
# Log(Mean) = 2.0852
# Log(StDev) = 0.8
calib6 <- getDescendants(non_calibrated_tree_ML.rooted, 331)


#calibration node 7 - Stem lineage Zeiformes, dating the MRCA of Zeiformes (Cyttopsis, Zenopsis, and Zeus), Gadiformes (Lota and Coryphaenoides), and Stylephorus chordatus.
getMRCA(non_calibrated_tree_ML.rooted,c("Zenopsis_conchifer","Zeus_faber", "Cyttopsis_rosea", "Lota_lota", "Coryphaenoides_rupestris", "Stylephorus_chordatus"))
# node 329
# Absolute age offset = 70.6 Ma 
# 95% soft upper bound = 125 Ma 
# Log(Mean) = 2.4278 
# Log(StDev) = 0.8
calib7 <- getDescendants(non_calibrated_tree_ML.rooted, 329)


#calibration node 8 - Stem lineage Lampridae dating the MRCA of Lampris, Regalecus, and Trachipterus.
getMRCA(non_calibrated_tree_ML.rooted,c("Lampris_guttatus","Regalecus_russelii", "Trachipterus_arcticus"))
# node 416
# Absolute age offset = 55.8 Ma 
# 95% soft upper bound = 125 Ma 
# Log(Mean) = 2.669
# Log(StDev) = 0.8
calib8 <- getDescendants(non_calibrated_tree_ML.rooted, 416)


#calibration node 9 -  Stem lineage Aphredoderidae, dating the MRCA of Aphredoderus and Chologaster.
getMRCA(non_calibrated_tree_ML.rooted,c("Aphredoderus_sayanus","Chologaster_cornuta"))
# node 335
# Absolute age offset = 34.1 Ma 
# 95% soft upper bound = 57 Ma 
# Log(Mean) = 1.563 
# Log(StDev) = 0.8
calib9 <- getDescendants(non_calibrated_tree_ML.rooted, 335)


#calibration node 10 - Stem lineage Percopsidae, dating the MRCA of Percopsiformes (Percopsis, Aphredoderus, and Chologaster)
getMRCA(non_calibrated_tree_ML.rooted,c("Percopsis_omiscomaycus","Aphredoderus_sayanus", "Chologaster_cornuta"))
# 334
# Absolute age offset = 57 Ma 
# 95% soft upper bound = 93.6 Ma 
# Log(Mean) = 2.032
# Log(StDev) = 0.8
calib10 <- getDescendants(non_calibrated_tree_ML.rooted, 334)


#calibration node 11 - Stem lineage Polymixiiformes, dating the MRCA of Polymixia and Percopsiformes (Percopsis, Aphredoderus, and Chologaster).
getMRCA(non_calibrated_tree_ML.rooted,c("Percopsis_omiscomaycus","Aphredoderus_sayanus", "Chologaster_cornuta", "Polymixia_japonica"))
# node 327
# Absolute age offset = 93.6 Ma 
# 95% soft upper bound = 220 Ma
# Log(Mean) = 1.88 
# Log(StDev) = 0.8
calib11 <- getDescendants(non_calibrated_tree_ML.rooted, 327)


#calibration node 12 - Stem lineage Alepisauridae, dating the MRCA of Alepisaurus, and Omosudis.
getMRCA(non_calibrated_tree_ML.rooted,c("Alepisaurus_ferox","Omosudis_lowii"))
# node 408
# Absolute age offset = 100 Ma
# 95% soft upper bound = 150 Ma 
# Log(Mean) = 2.344
# Log(StDev) = 0.8
calib12 <- getDescendants(non_calibrated_tree_ML.rooted, 408)


#calibration node 13 - Stem lineage Myctophiformes, dating the MRCA of Myctophiformes and Acanthomorpha. (USED Scopelengys and Abalistes).
getMRCA(non_calibrated_tree_ML.rooted,c("Scopelengys_tristis","Abalistes_stellatus"))
# node 336
# Absolute age offset = 96 Ma 
# 95% soft upper bound = 150 Ma 
# Log(Mean) = 2.421
# Log(StDev) = 0.8
calib13 <- getDescendants(non_calibrated_tree_ML.rooted, 336)


#calibration node 14 - Stem lineage Salmoninae, dating the MRCA of Thymallus, Salvelinus, and Coregonus.
getMRCA(non_calibrated_tree_ML.rooted,c("Thymallus_brevirostris","Salvelinus_alpinus", "Coregonus_clupeaformis"))
# node 345
# Absolute age offset = 51.8 Ma 
# 95% soft upper bound = 150 Ma 
# Log(Mean) = 3.019
# Log(StDev) = 0.8
calib14 <- getDescendants(non_calibrated_tree_ML.rooted, 345)


#calibration node 15 - Stem lineage Chanidae, dating the MRCA of Chanos and Cromeria.
getMRCA(non_calibrated_tree_ML.rooted,c("Chanos_chanos","Cromeria_nilotica"))
# node 376
# Absolute age offset = 133.9 Ma 
# 95% soft upper bound = 150 Ma 
# Log(Mean) = 1.26
# Log(StDev) = 0.8
calib15 <- getDescendants(non_calibrated_tree_ML.rooted, 376)


#calibration node 16 - Stem lineage Notopteridae, dating the MRCA of Gnathonemus, and Chitala.
getMRCA(non_calibrated_tree_ML.rooted,c("Gnathonemus_petersii","Chitala_chitala"))
# node 369
# Absolute age offset = 99.6 Ma
# 95% soft upper bound = 220 Ma 
# Log(Mean) = 3.223
# Log(StDev) = 0.8
calib16 <- getDescendants(non_calibrated_tree_ML.rooted, 369)


#calibration node 17 - Stem lineage Teleostei, dating the MRCA of Teleostei and Amia. (USED Abalistes_stellatus).
getMRCA(non_calibrated_tree_ML.rooted,c("Amia_calva","Abalistes_stellatus"))
# node 370
# Absolute age offset = 237 Ma
# 95% soft upper bound = 325.5 Ma
# Log(Mean) = 2.915
# Log(StDev) = 0.8
calib17 <- getDescendants(non_calibrated_tree_ML.rooted, 370)


#calibration node 18 - Stem Lineage Euteleostei, dating the MRCA of Euteleostei, Alepocephaloidei, Clupeiformes, and Ostariophysi. (USED lepidogalaxias and Abalistes_stellatus)
getMRCA(non_calibrated_tree_ML.rooted,c("Lepidogalaxias_salamandroides","Abalistes_stellatus"))
# node 351
# Absolute age offset = 150 Ma 
# 95% soft upper bound = 220 Ma 
# Log(Mean) = 2.345
# Log(StDev) = 0.8
calib18 <- getDescendants(non_calibrated_tree_ML.rooted, 351)


#calibration node 19 (ROOT NODE) -  Crown-group Actinopterygii (used polypterus and Abalistes_stellatus)
getMRCA(non_calibrated_tree_ML.rooted,c("Polypterus_ornatipinnis","Abalistes_stellatus"))
# node 228
# Absolute age offset = 382.5 Ma 
# 95% soft upper bound = 438 Ma 
# Log(Mean) = 2.17
# Log(StDev) = 0.8
calib19 <- getDescendants(non_calibrated_tree_ML.rooted, 228)


#calibration node 20 - Crown-group Actinopteri. (used Polyodon and Abalistes_stellatus)
getMRCA(non_calibrated_tree_ML.rooted,c("Polyodon_spathula","Abalistes_stellatus"))
# node 453
# Absolute age offset = 325.5 Ma 
# 95% soft upper bound = 382.5 Ma 
# Log(Mean) = 2.475
# Log(StDev) = 0.8
calib20 <- getDescendants(non_calibrated_tree_ML.rooted, 453)



n <- max(length(calib19))
length(calib1) <- n    
length(calib2) <- n 
length(calib3) <- n 
length(calib4) <- n 
length(calib5) <- n 
length(calib6) <- n 
length(calib7) <- n 
length(calib8) <- n 
length(calib9) <- n 
length(calib10) <- n 
length(calib11) <- n 
length(calib12) <- n 
length(calib13) <- n 
length(calib14) <- n 
length(calib15) <- n 
length(calib16) <- n 
length(calib17) <- n 
length(calib18) <- n 
length(calib20) <- n 


Calibration_taxa_groups <- cbind(calib1, calib2, calib3, calib4, calib5, calib6, calib7, calib8, calib9, calib10, calib11, calib12, calib13, calib14, calib15, calib16, calib17, calib18, calib19, calib20)

Calibration_taxa_groups <- as.data.frame(Calibration_taxa_groups)

write.csv(Calibration_taxa_groups, file = "Calibration_taxa_groups.csv")


#Build calibration data fram in order of calibration number (1-20)
calib <- data.frame(age.min = c(37.2,50,55.8,30.1,55.8,32,70.6,55.8,34.1,57,93.6,100,96,51.8,133.9,99.6,237,150,382.5,325.5), age.max = c(125,125,125,125,125,70.6,125,125,57,93.6,220,150,150,150,150,220,325.5,220,438,382.5), node = c(229,449,245,298,295,331,329,416,335,334,327,408,336,345,376,369,370,351,228,453))

#Transform to ultrametric
tr.dated = chronos(non_calibrated_tree_ML.rooted, calibration=calib)

#TEST ULTRAMETRIC

# check the root age
branching.times(tr.dated)[1]

# then plot
plot(tr.dated, cex=0.3)
simpleAxisPhylo()

write.tree(tr.dated, file="AllGenesUltrametric.tre")

#Get taxa and tip numbers
non_calibrated_tree_ML.rooted$tip.label

