# Supplementary file S4 for  
# Phylogenetic relationships of Neogene hamsters (Mammalia, Rodentia, Cricetinae) revealed under Bayesian inference and maximum parsimony
# Moritz Dirnberger1, Pablo Peláez-Campomanes2, Raquel López-Antoñanzas1
# 1ISEM, Univ Montpellier, CNRS, IRD, Montpellier, France
# 2Departamento de Paleobiología, Museo Nacional de Ciencias Naturales-CSIC, Madrid, Spain


	######## Identifying rogue taxa based on MrBayes output ########

# load libraries
library("TreeTools")
library("Rogue")

filename <- "infile"

# for the case of four runs
run1 <- ape::read.nexus(paste0(filename,".nex.run1.t"))
run2 <- ape::read.nexus(paste0(filename,".nex.run2.t"))
run3 <- ape::read.nexus(paste0(filename,".nex.run3.t"))
run4 <- ape::read.nexus(paste0(filename,".nex.run4.t"))

burninfrac <- 0.3
burningens <- as.integer(length(run1)*burninfrac):length(run1)
trees <- c(run1[burningens], run2[burningens], run3[burningens], run4[burningens])

plenary <- Consensus(trees, p = 0.5)

# plot instability of leaves
par(mar = rep(0, 4), cex = 0.85)
plot(plenary, tip.color = ColByStability(trees))
SpectrumLegend(0.06, 0.06, legend = c("Stable", "Unstable"),
               palette = hcl.colors(131, 'inferno')[1:101])

# identify rogue taxa
rogues <- RogueTaxa(trees, mreOptimization = TRUE)
rogueTaxa <- rogues$taxon[-1]

# compare consensus trees
windows(width = 20, height = 15)
par(mar = rep(0, 4), mfrow = c(1, 2))

reduced <- ConsensusWithout(trees, rogueTaxa, p = 0.5)

plot(plenary, tip.color = ifelse(plenary$tip.label %in% rogueTaxa, 2, 1))
LabelSplits(plenary, round(SplitFrequency(plenary, trees)/length(trees), 2), frame = "none")

plot(reduced)
LabelSplits(reduced, round(SplitFrequency(reduced, trees)/length(trees), 2), frame = "none")




	######## Assessing strategraphic fit of resulting topologies ########

# load libraries
library(strap)

# prepare .txt file with FAD and LAD of each taxon

tree<-read.nexus("CricetinaeMB_IGR_fossiltip_tree.t.con.tre")
ages<-read.table("Cricetinae_Ages.txt")

X <- StratPhyloCongruence(tree, ages, hard=TRUE, randomly.sample.ages=FALSE, fix.topology=TRUE, fix.outgroup=TRUE, outgroup="Eucricetodon_wangae")
write.table(X$input.tree.results, "File.csv", sep=",")




	######## Ancestral state reconstructions based on MrBayes output ########

# load libraries
library(ape)
library(phangorn)
library(phytools)

# delete spaces between character states in the file

characters <- read.nexus.data("infile.nex")
tree <- read.nexus("CricetinaeMB_IGR_fossiltip_fossiltree.t.con.tre")

# convert character data to data.frame
charMat <- t(data.frame(characters))
write.csv(charMat, "character_matrix.csv")
charMat <- read.csv("character_matrix.csv", row.names=1)
charMat <- charMat[tree$tip.label,]

# some functions needed later
addVal <- function(score, Len=4)	{
	out <- rep(0, Len)

	if (score == "?")	{
		out <- rep(1/Len, Len)
	}
	else if (grepl("/", score, fixed = TRUE))	{
		Splits <- as.numeric(strsplit(score, "/")[[1]])+1
		out[Splits] <- 0.5
	}
	else	{
		Spot <- as.numeric(score) + 1
		out[Spot] <- 1
	}
	return(out)
}
prepColumn <- function(trait)	{
	Levels <- unique(trait)
	if (class(Levels[1]) != "integer")	{
		# Break character scorings apart to get split scores (eg 0/1)
		Lvls <- sapply(Levels, strsplit, split="")
		# Render them numeric
		nLvls <- na.omit(as.numeric(unlist(Lvls)))
		# Find the biggest
		Len <- max(nLvls) + 1
	}
	else	{
		Len <- max(Levels) + 1
	}
	Mat <- t(sapply(trait, addVal, Len=Len))
	return(Mat)	
}

## Set colors to plot
Colors <- c('#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f','#ff7f00','#cab2d6','#6a3d9a')
names(Colors) <- 0:9


# Map trait X
Trait <- 1
toMap <- prepColumn(charMat[,Trait])
rownames(toMap) <- tree$tip.label
colnames(toMap) <- 0:(ncol(toMap)-1)
	
  ## test for which model is best-fitting
ER <- fitMk(tree, toMap, model="ER", nsim=100, pi="estimated")
SYM <- fitMk(tree, toMap, model="SYM", nsim=100, pi="estimated")
ARD <- fitMk(tree, toMap, model="ARD", nsim=100, pi="estimated")
anova(ER, SYM, ARD)	
	
  ## choose model		
Model <- "SYM"
	
sMap <- make.simmap(tree, toMap, model= Model, nsim=100, pi="estimated", Q="mcmc")

pdf(file = paste("stochasticMapofTrait_", Trait, Model, ".pdf", sep=""))
plot(describe.simmap(sMap), cex=0.5, colors=Colors, fsize=0.5)
add.simmap.legend(colors=Colors, x=0, y=38, vertical=F, prompt=F) # adjust x/y to look nice
title(main = paste("stochastic map of trait ", Trait, sep=""), line = -1, outer = TRUE)
dev.off()