I have to do this analysis as suggested by reviewers for the manuscript. I have to ancestral state reconstruction.
Here is the path to the folder where the analysis is done and files and codes are saved: /uufs/chpc.utah.edu/common/home/u6007910/projects/lmelissa_hostAdaptation/ancestral_state/
I copied the tree with 0 migration events from the Lmelissa_treemix folder. The file is the output file from Treemix.
library(ape)
library(phytools)
#read in the tree
mytree = read.tree("out_treemix0.treeout")
#check the tree
mytree
#get tree objects
str(mytree)
class(mytree) #phylo
#find the edge lengths or branch length of the tree
mytree$edge.length
#assign each object of the tree to some vector
branches <- mytree$edge
species <- mytree$tip.label
populations <- mytree$tip.label
brlength <- mytree$edge.length
nodes <- mytree$Nnode #27
#designate character state (Medicago(1) or native(0) for each tip label (population)). Note the outgroup are given 0 for native plant
char.tree <- c(0,0,1,0,1,1,1,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,1,1,0,1,1)
length(char.tree) #28 for 28 populations
names(char.tree) <- species
#since the branch lengths are lesser than 0, these lengths need to be modified to run the ancestral state reconstruction as I was getting an error. To do this, the recommended is using the compute.brlen() funtion in library(ape)
##now this is the new tree with modified branch lengths
newtree <- compute.brlen(mytree)
## use this tree to do further analysis for ancestral state recontruction
#using ace function to do the ancestral state reconstruction to get posterior probabilites
#contructing marginal likelihoods using all rates equal model
er.tree <- ace(char.tree, newtree, type ="discrete", model="ER")
#using simmap to get the tree
mytree.anc <- make.simmap(newtree, x = char.tree, model="ARD")
mytree.mapsum <- describe.simmap(mytree.anc, plot=FALSE)
#plotting the tree with posterior probabilies
cols <-setNames(palette()[1:length(unique(char.tree))], sort(unique(char.tree)))
#using the estimated posterior probabilities from er model
round(er.tree$lik.anc, 3)
plotTree(newtree)
nodelabels(node=1:newtree$Nnode+Ntip(newlength.brlen), pie = er.tree$lik.anc, piecol=cols, cex=0.5)
tiplabels(pie=to.matrix(char.tree, sort(unique(char.tree))), piecol=cols, cex=0.3)
add.simmap.legend(colors=cols, prompt=FALSE, x=0.9*par()$usr[1], y = -max(nodeHeights(newlength.brlen)), fsize=0.8)
mytree.mapsum <- describe.simmap(mytree.anc, plot=FALSE)
mytree.mapsum