#Loading a tree and performing and ancestral character state reconstruction
library(diversitree)
library(phytools)
library(paleotree)
library(ggtree)
library(geiger)

#Load data

Fluorescence_data <- read.csv("BiofluorChar_SpToFamNames.csv")
summary(Fluorescence_data)

Fluorescence_tree_FamNames <- read.tree("100mil_trim_FamNames.tre")
plot.phylo(Fluorescence_tree_FamNames)
Fluorescence_tree_FamNames <- ladderize(Fluorescence_tree_FamNames) 

#--------ANCESTRAL CHARACTER STATE RECONSTRUCTION

#----RedOnly Data-----

# Using the heterodonty of each specimen
Fluorescence_Red_YN <- Fluorescence_data[,4]
Fluorescence_Red_YN

Fluorescence_Red_YN <- setNames(Fluorescence_Red_YN, Fluorescence_data$Family)
Fluorescence_Red_YN
str(Fluorescence_Red_YN)

Fluorescence_Red_YN_Tree_ER <- ace(Fluorescence_Red_YN, Fluorescence_tree_FamNames, model = "ER", type = "discrete", CI = TRUE)
Fluorescence_Red_YN_Tree_ER
round(Fluorescence_Red_YN_Tree_ER$lik.anc,3)

#----GreenOnly Data-----
Fluorescence_Green <- Fluorescence_data[,5]
Fluorescence_Green

Fluorescence_Green <- setNames(Fluorescence_Green, Fluorescence_data$Family) 
str(Fluorescence_Green)

Fluorescence_Green_Tree <- ace(Fluorescence_Green, Fluorescence_tree_FamNames, model = "ER", type = "discrete", CI = TRUE)
Fluorescence_Green_Tree
round(Fluorescence_Green_Tree$lik.anc,3)

#----Biofluor Absense/Pres Data-----
Fluorescence_Biofluor <- Fluorescence_data[,9]
Fluorescence_Biofluor

Fluorescence_Biofluor <- setNames(Fluorescence_Biofluor, Fluorescence_data$Family) 
str(Fluorescence_Biofluor)

Fluorescence_Biofluor_Tree <- ace(Fluorescence_Biofluor, Fluorescence_tree_FamNames, model = "ER", type = "discrete", CI = TRUE)
Fluorescence_Biofluor_Tree
round(Fluorescence_Biofluor_Tree$lik.anc,3)

#----Red and Green Data-----
Fluorescence_RandG <- Fluorescence_data[,8]
Fluorescence_RandG

Fluorescence_RandG <- setNames(Fluorescence_RandG, Fluorescence_data$BINO) 
str(Fluorescence_RandG)

Fluorescence_RandG_Tree <- ace(Fluorescence_RandG, Fluorescence_tree_FamNames, model = "ER", type = "discrete", CI = TRUE)
Fluorescence_RandG_Tree
round(Fluorescence_RandG_Tree$lik.anc,3)

#--------PLOT------
quartz(width=15, height = 15)

plot.phylo(Fluorescence_tree_FamNames, type = "fan", cex = .3, lwd= 0.4, rotate.tree = 180, 
           use.edge.length = F, label.offset = 0.04, node.width = 0.2, node.depth = 2)

nodelabels(pie = Fluorescence_Red_YN_Tree$lik.anc, cex = 0.11, piecol = c('white','red'))
nodelabels(pie = Fluorescence_Green_Tree$lik.anc, cex = 0.11, piecol = c('white','green'))
nodelabels(pie = Fluorescence_Biofluor_Tree$lik.anc, cex = 0.11, piecol = c('white', 'skyblue'))
nodelabels(pie = Fluorescence_RandG_Tree$lik.anc, cex = 0.11, piecol = c('white', 'darkblue'))

