Commit a8753c6f authored by Jannes Höpke's avatar Jannes Höpke

Upload New File

parent 8abcb101
###### R Script for paper "Phylogenetic and morphometric analysis of Plantago section
###### Coronopus (Plantaginaceae)" by Jannes Hoepke, Ladislav Mucina and Dirk C. Albach
###### Date: 11 August 2015
###### Date of revision: 25 February 2019
###### Script is based on ...
# 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)
# ... but was modified by Hoepke, J.
# This script was created to plot the support values (summarised by SumTrees.py) on the "best"
# (= most likely) tree from GARLI. For doing so, please copy supplementary Appendix files
# S5 to S23 in one folder (S5-S11 are for ITS, S12-S17 for trnLF, and S18-S23 for cpDNA).
################################# Initial set-up #######################################
### Loading packages
library(ape)
library(geiger)
library(phytools)
### Print information about used package versions (not necessarily matching with citations in the
# publication but it shows the program versions during my last successful run during revision
# of this file)
sink("Appendix_S10_Hoepke_et_al_tree_plotting_RSessionInfo.txt")
print(sessionInfo())
sink()
### Create (sub-)directories, e.g. for storing plots
system("mkdir Plots")
########################### Make preparations ##################################
### 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("ITS_best.tre")->bestTree #FigTree output (rooted)
read.nexus("ITS_mb_consensus.tre")->bayesTree #SumTrees.py output
read.nexus("ITS_garli_consensus.tre")->bootTree #SumTrees.py output
### Setting plot options
opar<-par()
par(mar=c(0,0,0,0),xpd=T)
########################### Root all trees ##################################
### Reroot bestTree
plot(bestTree,cex=tcex)
### Adjust bestTree via "rotate"
nodelabels()
rrbestTree<-rotate(bestTree, 61, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 62, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 63, polytom = c(2,3))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 72, polytom = c(1,3))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 73, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 74, polytom = c(1,3))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 79, polytom = c(3,13))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 89, polytom = c(3,13))
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()
i<-67
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, 68, polytom = c(1,2))
plot(rrbootTree,show.node.label=T,cex=tcex)
nodelabels()
rrbootTree<-rotate(rrbootTree, 72, polytom = c(1,2))
plot(rrbootTree,show.node.label=T,cex=tcex)
nodelabels()
rrbootTree<-rotate(rrbootTree, 74, polytom = c(1,2))
plot(rrbootTree,show.node.label=T,cex=tcex)
nodelabels()
rrbootTree<-rotate(rrbootTree, 75, polytom = c(18,19))
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<-67
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, 68, polytom = c(1,2))
plot(rrbayesTree,show.node.label=T,cex=tcex)
nodelabels()
rrbayesTree<-rotate(rrbayesTree, 75, polytom = c(1,2))
plot(rrbayesTree,show.node.label=T,cex=tcex)
nodelabels()
rrbayesTree<-rotate(rrbayesTree, 80, polytom = c(1,13))
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)
########################### Prepare comparisons of trees ##################################
### 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("ITS_bestTree.tiff", width = 20, height = 25, units = 'cm', res = 300)
par(mar=c(2,0,0,9),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(33)]<-""
bootNodes2<-rep("",length(bootNodes1))
bootNodes2[33]<-"100"
bayesNodes1<-bayesList
bayesNodes1[c(33)]<-""
bayesNodes2<-rep("",length(bayesList))
bayesNodes2[33]<-"1.00"
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.3, -0.3), frame="n", cex=scex, font=1)
nodelabels(bayesNodes2, adj=c(2.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