###################################################################################################################################################
###	R script for consensus trees, a posteriori timescaling, and phenogram																		###
### Written by V. Fischer 2020																													###
###	Associated paper:																															###
###	Fischer et al. A new elasmosaurid plesiosaurian from the Early Cretaceous of Russia marks an early attempt at neck elongation. ZJLS XXDOIXX	###
###################################################################################################################################################


#================#
#=== Packages ===#
#================#

library(ape)
library(paleotree)
library(strap)
library(phytools)
library(ggplot2)
library(psych)


#=========================#
#=== Working directory ===#
#=========================#

setwd("~/")


#===========================================================#
#=== Import trees and compute + timescale consensus tree ===#
#===========================================================#

###	Get phylogenetic data and drop tip according to the taxa passing the completeness threshold	###
trees_IW <- read.nexus("SUPP_matrix_trees_IW.nex")
trees_EW <- read.nexus("SUPP_matrix_trees_EW.nex")
ranges <- read.table("SUPP_ranges.txt", header=TRUE)
consensus_IW <- timePaleoPhy(consensus(trees_IW), timeData= ranges, type= "equal", vartime = 3)
consensus_IW_50 <- timePaleoPhy(consensus(trees_IW,p=0.5), timeData= ranges, type= "equal", vartime = 3)
consensus_EW <- timePaleoPhy(consensus(trees_EW), timeData= ranges, type= "equal", vartime = 3)
consensus_EW_50 <- timePaleoPhy(consensus(trees_EW,p=0.5), timeData= ranges, type= "equal", vartime = 3)


###	Timescaled figure ###
pdf("FigS1_consensus_IW_holotype.pdf",width=9, height=11)
geoscalePhylo(ladderize(consensus_IW), ranges, cex.ts=0.5, cex.tip=0.3)
#geoscalePhylo(ladderize(consensus_IW_50), ranges, cex.ts=0.5, cex.tip=0.3) #to get the 50% majority rule instead
dev.off()

pdf("FigS2_consensus_EW_holotype.pdf",width=9, height=11)
geoscalePhylo(ladderize(consensus_EW), ranges, cex.ts=0.5, cex.tip=0.3)
#geoscalePhylo(ladderize(consensus_EW_50), ranges, cex.ts=0.5, cex.tip=0.3) #to get the 50% majority rule instead
dev.off()



#=====================#
#=== Cervical data ===#
#=====================#


###	Get randomly selected most parsimonious tree and timescale it	###
random_tree <- timePaleoPhy(trees_IW[[sample(1:length(trees_IW),1)]], timeData=cbind(rowMeans(ranges),rowMeans(ranges)-0.1), type= "mbl", vartime = 5) #cbind(rowMeans(ranges),rowMeans(ranges)-0.1) is used to get the data point in the middle of the range for the phenogram, otherwise, it is put by default at the FAD by the phenogram function


###	Get and tidy anterior cervical data	###

full_cervical_data <- read.csv("SUPP_cervical_data.csv",header=TRUE,row.names=1)[,3:11]

#choose which dataset (add a '#' in front unwanted one)
cervical_data <- full_cervical_data[,1:4] ; choose <- "c10" #for cervical 10 dataset
#cervical_data <- full_cervical_data[,c(1,6:8)] ; choose <- "longest" #for relatively longest cervical dataset


cervical_data$L_H <- cervical_data[,3]/cervical_data[,2]
cervical_data$H_W <- cervical_data[,2]/cervical_data[,4]
L_H_phylo <- cervical_data[rownames(cervical_data) %in% random_tree$tip.label,"L_H"]
names(L_H_phylo) <- rownames(cervical_data[rownames(cervical_data) %in% random_tree$tip.label,])
cervical_data_supp <- cervical_data[!(rownames(cervical_data) %in% random_tree$tip.label),]

#taxa with similar values than Jucha:
subset(cervical_data, L_H>1.32 & L_H<1.39)


###	Create dataset for histograms
histo_data <- data.frame(L_H=c(cervical_data[cervical_data$Time>100.5,"L_H"],na.omit(cervical_data[cervical_data$Time<100.5,"L_H"])),bin=c(rep("EK",length(cervical_data[cervical_data$Time>100.5,"L_H"])),rep("LK",length(na.omit(cervical_data[cervical_data$Time<100.5,"L_H"])))))
histo_data <- histo_data[-c(1,2),] #remove non-elasmosaurids
range(histo_data[histo_data$bin=="EK",1]) #Early Cretaceous range
range(histo_data[histo_data$bin=="LK",1]) #Late Cretaceous range

# Histogram figure
ggplot(histo_data, aes(x=L_H, fill=bin,color=bin)) +
  geom_density(alpha=0.4, position = 'identity') +
  scale_colour_manual(values = c("#94cc79","#afd46c"))+
  scale_fill_manual(values = c("#94cc79","#afd46c"))+
  labs(title="Distribution of cervical elongation values",x="Length/height ratio")+
  theme(legend.position="none",panel.grid.major = element_blank(), panel.grid.minor = element_blank())
#ggsave(paste("Fig10_S2_histo_LH_",choose,".pdf",sep=""),device="pdf",units="mm",width=180,height=70)


###	Phenograms

#prune randomly selected tree to match taxon selection
random_tree_pruned <- drop.tip(random_tree,random_tree$tip.label[!(random_tree$tip.label %in% rownames(cervical_data))])

#compute ancestral values for nodes using maximum likelihood
anc_L_H <-  fastAnc(random_tree_pruned, L_H_phylo)

#Extract data from the phenogram function to add temporal information on the plots
pheno_data <- phenogram(random_tree_pruned,c(L_H_phylo,anc_L_H),spread.labels=FALSE) #!! if spread.labels is set on TRUE (default) then the x data is that of the label position, not the end point of the branch!!
root <- rowMeans(ranges)["Cardiocorax_mukulu"]+pheno_data["Cardiocorax_mukulu",1]  #the general value from which temporal should be substracted to be placed in the correct position on the graph; this is a weird behaviour of the phenogram function
supp_to_plot <- cbind(root-cervical_data_supp[!is.na(cervical_data_supp[,2]),"Time"],cervical_data_supp[!is.na(cervical_data_supp[,2]),"L_H"])
rownames(supp_to_plot) <- rownames(cervical_data_supp[!is.na(cervical_data_supp$H),]) ; colnames(supp_to_plot) <- c("x","y")

#Phenogram figure

#pdf(paste("Fig10-S2_phenogram_",choose,".pdf",sep=""), width=8, height=7)
phenogram(random_tree_pruned,c(L_H_phylo,anc_L_H),ylim=c(min(L_H_phylo)-0.18,max(L_H_phylo)),ylab="Anterior cervical elongation",fsize = 0.7,ftype="i",spread.labels=TRUE)
rect(xleft=root-145,ybottom=min(L_H_phylo)-0.05,xright=root-139.8,ytop=max(L_H_phylo)+0.05,col="#cccccc50",border="NA") #Berriasian
rect(xleft=root-132.9,ybottom=min(L_H_phylo)-0.05,xright=root-129.4,ytop=max(L_H_phylo)+0.05,col="#cccccc50",border="NA") #Hauterivian
rect(xleft=root-125,ybottom=min(L_H_phylo)-0.05,xright=root-113,ytop=max(L_H_phylo)+0.05,col="#cccccc50",border="NA") #Aptian
rect(xleft=root-100.5,ybottom=min(L_H_phylo)-0.05,xright=root-93.9,ytop=max(L_H_phylo)+0.05,col="#cccccc50",border="NA") #Cenomanian
rect(xleft=root-89.8,ybottom=min(L_H_phylo)-0.05,xright=root-86.3,ytop=max(L_H_phylo)+0.05,col="#cccccc50",border="NA") #Coniacian
rect(xleft=root-83.6,ybottom=min(L_H_phylo)-0.05,xright=root-72.1,ytop=max(L_H_phylo)+0.05,col="#cccccc50",border="NA") #Campanian
rect(xleft=root-145,ybottom=min(L_H_phylo)-0.15,xright=root-100.5,ytop=min(L_H_phylo)-0.05,col="#94cc79",border="NA") #Early K rectangle
rect(xleft=root-100.5,ybottom=min(L_H_phylo)-0.15,xright=root-66,ytop=min(L_H_phylo)-0.05,col="#afd46c",border="NA") #Late K rectangle
text(x=c(root-122.74,root-83.25),y=rep(min(L_H_phylo)-0.1,2),labels=c("Early Cretaceous","Late Cretaceous"),cex=0.7) 
text(x=c(root-145,root-100.5,root-66),y=rep(min(L_H_phylo)-0.17,3),labels=c(145,100.5,66),cex=0.7)
points(pheno_data,pch=20,col=c(rep("#01abe9",2),rep("#ff952c",3),rep("#FF0D68",2),rep("#ff952c",7),rep("#23CC8F",2),rep("#ff952c",4),rep("#23CC8F",2))) #greenish "#23CC8F" = non-euelasmos – # orange "#ff952c" = euelasmos – pink "#FF0D68" = aristonectines
points(supp_to_plot,pch=20,col="grey")
text(x=c(rep(0,4)),y=c(1.67,1.64,1.61,1.58,1.55),labels=c("leptocleidians","early elasmosaurids","euelasmosauridans","aristonectine euelasmosauridans","indeterminate elasmosaurids"),col=c("#01abe9","#23CC8F","#ff952c","#FF0D68","grey"),pos=4,cex=0.7) #color legend
#dev.off()


#============================================================#
#=== Modified neck data from Soul & Benson 2017 Evolution ===#
#============================================================#
SB2017 <- read.csv("SUPP_S&B2017_long_necked.csv",na="NA",row.names=1)

cor.test(SB2017$Cervical,SB2017$contrib, method="pearson")

pairs.panels(SB2017,method="pearson",pch=21,lm=TRUE,hist.col="grey",density=FALSE,ellipses=FALSE,stars=TRUE)



