Commit 82c373c2 authored by Jannes Höpke's avatar Jannes Höpke

Upload New File

parent 057ccae5
### Script by:
### GLOR, Rich (2008). "R Tip: Labeling Trees w/ Posterior Probability and Bootstrap Support".
### http://treethinkers.blogspot.de/2008/10/labeling-trees-posterior-probability.html
### (access date: 4th August 2015)
### Modified by:
### HÖPKE, Jannes
#Set working directory
setwd("~/Documents/Uni Oldenburg/6. Semester/Bachelorarbeitsmodul/Phylogenetische Bäume in R/comb_2")
#Fist we need to open some necessary libraries
library(ape)
library(geiger)
#Some plot options
#define cex-size for tiplabel
tcex<-0.7
#define cex-size for support
scex<-0.6
#define cex-size for tree line
lcex<-1.5
#Read in trees
read.nexus("comb_best.tre")->bestTree #FigTree output (rooted)
read.nexus("comb_mb_consensus.tre")->bayesTree #SumTrees.py output
read.nexus("comb_garli_consensus.tre")->bootTree #SumTrees.py output
#Setting plot options
opar<-par()
par(mar=c(0,0,0,0),xpd=T)
#Reroot bestTree
plot(bestTree,cex=tcex)
library(phytools)
#Adjust bestTree via "rotate"
nodelabels()
rrbestTree<-rotate(bestTree, 64, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 88, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 89, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 77, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 78, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 80, polytom = c(1,3))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 69, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 70, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 96, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
#Create new labels for bestTree
bestLabels<-read.csv("bestTree_tips.csv", header=T, sep=";")
str(bestLabels)
attach(bestLabels)
bestLabels<-data.frame(lapply(bestLabels, as.character), stringsAsFactors=FALSE)
#converts all colums from a dataframe into characters
detach(bestLabels)
attach(bestLabels)
str(bestLabels)
bestTips <- mixedFontLabel(new1," ",new2,", ",new3," ",new4, italic = 1:3, parenthesis = 7, sep="")
#Reroot bootTree
plot(bootTree,cex=tcex)
nodelabels()
tiplabels()
library(phytools)
i<-70
rbootTree <- reroot(bootTree, i, position = 0.5 * bootTree$edge.length[which(bootTree$edge[,2] == i)])
plot(rbootTree,show.node.label=T,cex=tcex)
#Adjust bootTree via "rotate"
nodelabels()
rrbootTree<-rotate(rbootTree, 76, polytom = c(2,3))
plot(rrbootTree,show.node.label=T,cex=tcex)
nodelabels()
rrbootTree<-rotate(rrbootTree, 76, polytom = c(1,2))
plot(rrbootTree,show.node.label=T,cex=tcex)
nodelabels()
rrbootTree<-rotate(rrbootTree, 77, polytom = c(1,2))
plot(rrbootTree,show.node.label=T,cex=tcex)
nodelabels()
rrbootTree<-rotate(rrbootTree, 78, polytom = c(4,6))
plot(rrbootTree,show.node.label=T,cex=tcex)
#Delete "Root" labeling
rrbootTree$node.label[which(rrbootTree$node.label=="Root")]<-""
plot(rrbootTree,show.node.label=T,cex=tcex)
#Reroot bayesTree
plot(bayesTree, show.node.label=T,cex=tcex)
nodelabels()
i<-70
rbayesTree <- reroot(bayesTree, i, position = 0.5 * bayesTree$edge.length[which(bayesTree$edge[,2] == i)])
plot(rbayesTree,show.node.label = T,cex=tcex)
#Adjust bootTree via "rotate"
nodelabels()
rrbayesTree<-rotate(rbayesTree, 77, polytom = c(1,2))
plot(rrbayesTree,show.node.label=T,cex=tcex)
#Delete "Root" labeling
rrbayesTree$node.label[which(rrbayesTree$node.label=="Root")]<-""
plot(rrbayesTree,show.node.label=T,cex=tcex)
#The getAllSubTrees function below is a necessary subfunction that atomizes a tree into
#each individual subclade and was provided compliments of Luke Harmon.
getAllSubtrees<-function(phy, minSize=2)
{
res<-list()
count=1
ntip<-length(phy$tip.label)
for(i in 1:phy$Nnode)
{
l<-tips(phy, ntip+i)
bt<-match(phy$tip.label, l)
if(sum(is.na(bt))==0)
{
st<-phy
}
else st<-drop.tip(phy, phy$tip.label[is.na(bt)])
if(length(st$tip.label)>=minSize)
{
res[[count]]<-st
count<-count+1
}
}
res
}
#The plotBootBayes function below plots both posterior probability and bootstrap values
#on each node of the consensus tree obtained from your Maximum-Likelihood analysis.
#Bootstrap values will appear above and to the left of the node they support, whereas
#Bayesian posterior probabilies will appear below and to the left of the node.
tiff("comb_bestTree.tiff", width = 20, height = 25, units = 'cm', res = 300)
par(mar=c(2,0,0,8),xpd=T)
getAllSubtrees(rrbayesTree)->bayesSub
getAllSubtrees(rrbootTree)->bootSub
getAllSubtrees(rrbestTree)->bestSub
bayesList<-matrix("",Nnode(rrbestTree),1)
bootList<-matrix("",Nnode(rrbestTree),1)
#The commands below compare all the subclades in the Bayes tree to all the subclades
#in the bootstrap tree, and vice versa, and identifies all those clades that are
#identical.
for(i in 1:Nnode(rrbestTree)) {
for(j in 1:Nnode(rrbayesTree)) {
match(bestSub[[i]]$tip.label[order(bestSub[[i]]$tip.label)], bayesSub[[j]]$tip.label[order(bayesSub[[j]]$tip.label)])->shared
match(bayesSub[[j]]$tip.label[order(bayesSub[[j]]$tip.label)], bestSub[[i]]$tip.label[order(bestSub[[i]]$tip.label)])->shared2
if(sum(is.na(c(shared,shared2)))==0) {
rrbayesTree$node.label[j]->bayesList[i]
}}}
for(i in 1:Nnode(rrbestTree)) {
for(j in 1:Nnode(rrbootTree)) {
match(bestSub[[i]]$tip.label[order(bestSub[[i]]$tip.label)], bootSub[[j]]$tip.label[order(bootSub[[j]]$tip.label)])->shared
match(bootSub[[j]]$tip.label[order(bootSub[[j]]$tip.label)], bestSub[[i]]$tip.label[order(bestSub[[i]]$tip.label)])->shared2
if(sum(is.na(c(shared,shared2)))==0) {
rrbootTree$node.label[j]->bootList[i]
}}}
plot(rrbestTree, cex=tcex, edge.width =lcex, show.tip.label=F) #Plots your Maximum-Likelihood consensus tree
tiplabels(bestTips, adj = -0.02, frame = "none", col = "black", bg="white",cex=tcex)
bootNodes1<-bootList
bootNodes1[c(20,21)]<-""
bootNodes2<-rep("",length(bootNodes1))
bootNodes2[20]<-"85"
bootNodes3<-rep("",length(bootNodes1))
bootNodes3[21]<-"61"
bayesNodes1<-bayesList
bayesNodes1[c(20,21)]<-""
bayesNodes2<-rep("",length(bayesList))
bayesNodes2[20]<-"1.00"
bayesNodes3<-rep("",length(bayesList))
bayesNodes3[21]<-"0.92"
nodelabels(bootNodes1, adj=c(1.2, -0.3), frame="n", cex=scex, font=1)
nodelabels(bayesNodes1, adj=c(1.12, 1.3), frame="n", cex=scex, font=1)
nodelabels(bootNodes2, adj=c(2.5, -0.3), frame="n", cex=scex, font=1)
nodelabels(bayesNodes2, adj=c(1.8, 1.3), frame="n", cex=scex, font=1)
nodelabels(bootNodes3, adj=c(1, -0.3), frame="n", cex=scex, font=1)
nodelabels(bayesNodes3, adj=c(1, 1.3), frame="n", cex=scex, font=1)
add.scale.bar(x = 0, y = -1, length = NULL, cex=tcex)
dev.off()
# --> graphical output!
par(opar)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment