Appendix_S16_Hoepke_et_al_trnLF_plotting.r 6.47 KB
Newer Older
Jannes Höpke's avatar
Jannes Höpke committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188
### 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/trnLF_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("trnLF_best.tre")->bestTree #FigTree output (rooted)
read.nexus("trnLF_mb_consensus.tre")->bayesTree #SumTrees.py output
read.nexus("trnLF_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, 59, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 60, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 61, polytom = c(1,2))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 67, polytom = c(1,7))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 67, polytom = c(2,10))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 67, polytom = c(3,11))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 67, polytom = c(4,6))
plot(rrbestTree,show.node.label=F,cex=tcex)
nodelabels()
rrbestTree<-rotate(rrbestTree, 67, polytom = c(5,7))
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, 72, polytom = c(3,4))
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<-63
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, 64, polytom = c(1,2))
plot(rrbootTree,show.node.label=T,cex=tcex)
nodelabels()
rrbootTree<-rotate(rrbootTree, 72, polytom = c(5,9))
plot(rrbootTree,show.node.label=T,cex=tcex)
nodelabels()
rrbootTree<-rotate(rrbootTree, 66, polytom = c(1,8))
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<-63
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, 64, polytom = c(1,2))
plot(rrbayesTree,show.node.label=T,cex=tcex)
nodelabels()
rrbayesTree<-rotate(rrbayesTree, 73, polytom = c(1,3))
plot(rrbayesTree,show.node.label=T,cex=tcex)
nodelabels()
rrbayesTree<-rotate(rrbayesTree, 68, 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("trnLF_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)

nodelabels(bootList, adj=c(1.2, -0.3), frame="n", cex=scex, font=1) 
nodelabels(bayesList, adj=c(1.12, 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)