#Tree loding and manipulations

#1) load library

library(ape)

#2) change working directory to the want that you will use for input/output files

setwd("/Users/jcsantos/Desktop/R_class_winter_2015_home/0_Bio559R_course_final_files/week_6/data")

#3) We can read parenthetical format:

#File 'CYTB_best_tree_garli.newick. Notice the quotations and the ; at the end.

CYTB_best_tree_garli_vector <- '(((AF020227:1.0E-8,AF020230:0.00134935):0.16805302,AJ278511:0.23257378):0.07250545,(((AB218960:0.66999926,(AB218884:0.3772031,(JF718358:0.99125309,AB218883:0.47313352):0.08306359):0.04297558):0.07585159,(((EF653285:0.24354808,GQ272812:0.36068284):0.17224242,KC853797:0.59306894):0.05749033,(((KJ504993:0.21865681,KJ505798:0.27311239):0.12195138,HQ332321:0.44616643):0.23555534,(DQ316873:0.46648512,EU443141:0.52343904):0.07951119):0.16423278):0.05731512):0.03191975,(EU129908:0.96285016,EU116519:0.66866705):0.10784251):0.32820119,(HQ141260:0.29554602,(FJ535926:0.02431339,FJ535913:0.00650359):0.20829292):0.05159107);'
tree <- read.tree(text = CYTB_best_tree_garli_vector)

#4) Alternatively, you can read trees in the R-environment using 'read.tree' function

cytb_tree_garli <-read.tree(file = "CYTB_best_tree_garli.newick")
cat(readLines("CYTB_best_tree_garli.newick"), sep = "\n") # read the contents of the file

#5) We can plot both trees side by side

par(mfrow=c(1,2))
plot(tree, type = "phylogram", edge.width = 2)
plot(cytb_tree_garli, type = "phylogram", edge.width = 2)

#6) We can get visually identify the nodes of our tree by plotting

str(tree)

#7) We can explore each component:

tree$edge # Provides a matrix of node and tip numbers. It is numbered as follows: 
          # tips: 1 to the total number of ‘n’ tips (extant species)
          #  nodes: ‘n+1’ to the total number of nodes  

tree$Node #Indicates that the number of internal nodes in this case NULL (tree is unrooted)
tree$tip.label #names of taxa at the tips
tree$edge.length #branch length values


#8) Our tree is unrooted, so we can check this tree conformation and reroot this tree based on specific taxa

is.rooted(tree) #[1] FALSE
outgroup <- c("JF718358", "AB218883", "AB218884", "AB218960") # I choose this clade 
rooted_tree <- root(tree, outgroup, resolve.root=TRUE)
plot(rooted_tree)
is.rooted(rooted_tree) #[1] TRUE

#9) We can get visually identify the nodes of our tree by plotting

plot(rooted_tree, type = "phylogram", label.offset = 0.05, edge.width = 2)
nodelabels()# labeled in lightblue 
tiplabels()# labeled in yellow

#10) We can also can plot different formats based on ‘type’

par(mfrow=c(1,3)) #this will create a graphics space: 1 row, 3 columns
plot(rooted_tree, type = "phylogram", edge.width = 2, cex = 1, label.offset = 0.05)
add.scale.bar(cex = 1, font = 2, col = "red")
plot(rooted_tree, type = "cladogram", edge.width = 2, cex = 1, label.offset = 0.05)
plot(rooted_tree, type = "fan", edge.width = 2, cex = 1, label.offset = 0.05)

#11) We can do a series of changes on our tree based on the branch (edges) lengths:

#Collapse short branches:

rooted_tree_collapsed <- di2multi(rooted_tree, tol = 0.1) #tol provide the branch length min value

par(mfrow=c(1,2))
plot(rooted_tree, type = "phylogram", edge.width = 2)
add.scale.bar(cex = 1, font = 2, col = "red")
plot(rooted_tree_collapsed)
add.scale.bar(cex = 1, font = 2, col = "red")

#12) We can do a series of changes on our tree based on the branch (edges) lengths:

#Change branch lengths:

rooted_tree_bl_1 <- compute.brlen(rooted_tree, 1)  # all branches equal to 1
rooted_tree_bl_square <- compute.brlen(rooted_tree, power=2)  # all branches squared

rooted_tree_bl_abs_log_10 <-rooted_tree
rooted_tree_bl_abs_log_10$edge.length <- abs(log10(rooted_tree_bl_abs_log_10$edge.length+1))

par(mfrow=c(1,4))
plot(rooted_tree, type = "phylogram", edge.width = 2)
add.scale.bar(cex = 1, font = 2, col = "red")
plot(rooted_tree_bl_1)
add.scale.bar(cex = 1, font = 2, col = "red")
plot(rooted_tree_bl_square)
add.scale.bar(cex = 1, font = 2, col = "red")
plot(rooted_tree_bl_abs_log_10)
add.scale.bar(cex = 1, font = 2, col = "red")


#13) We can do a series of changes on our tree based on the tips:

#Drop taxa:

rooted_tree_reduced <- drop.tip(rooted_tree, tip = c("AB218884", "AB218883", "JF718358"))

par(mfrow=c(1,2))
plot(rooted_tree, type = "phylogram", edge.width = 2)
add.scale.bar(cex = 1, font = 2, col = "red")
plot(rooted_tree_reduced)
add.scale.bar(cex = 1, font = 2, col = "red")

#14) We can do a series of changes on our tree based on the tips:

#Identify internal nodes (MRCA-- most recent common ancestor):

getMRCA(rooted_tree, tip = c("AB218884", "AB218883", "JF718358"))# ID MRCA node 38

#15) We can do a series of changes on our tree appearance. For some of  these, we will need other handy R package 

install.packages("phytools")
library("phytools")

# Rotate Nodes:
rooted_tree_rotated <- rotate(rooted_tree, node = 22)

par(mfrow=c(1,2))
plot(rooted_tree_rotated, type = "phylogram", edge.width = 2)
add.scale.bar(cex = 1, font = 2, col = "red")
plot(rooted_tree, edge.width = 2)
add.scale.bar(cex = 1, font = 2, col = "red")

#16) Add tip to a tree:

# select place to add the new tip: 

which(rooted_tree$tip.label=="EU443141") # [1] 15 number of this tip

# get a list of tips and their branch lengths

setNames(rooted_tree$edge.length[sapply(1:length(rooted_tree$tip.label),function(x,y) which(y==x),y=rooted_tree$edge[,2])],rooted_tree$tip.label)

# EU443141
# 0.52343904

rooted_tree_plus_1_tip <- bind.tip(rooted_tree, tip.label = "New_tip", edge.length= 0.52343904, where=15, position=0.5)

par(mfrow=c(1,2))
plot(rooted_tree, type = "phylogram", edge.width = 2)
edgelabels(round(rooted_tree$edge.length,4)) #round values to 4 decimals
add.scale.bar(cex = 1, font = 2, col = "red")
plot(rooted_tree_plus_1_tip, edge.width = 2)
edgelabels(round(rooted_tree_plus_1_tip$edge.length,4))
add.scale.bar(cex = 1, font = 2, col = "red")

#17) Extract clade:

getMRCA(rooted_tree, tip = c("EU443141", "KJ504993"))# ID MRCA node is 33
Clade_interest <- extract.clade(rooted_tree, node=33, root.edge = 0)

par(mfrow=c(1,2))
plot(rooted_tree, type = "phylogram", edge.width = 2)
add.scale.bar(cex = 1, font = 2, col = "red")
plot(Clade_interest, edge.width = 2)
add.scale.bar(cex = 1, font = 2, col = "red") 

#18) Paste two trees:

tree_1 <- rooted_tree
tree_1$tip.label[1]<-"NA" #this node is for “AF020227”
tree_2 <- rooted_tree
tree_2$root.edge<-0.1
double_tree <- paste.tree(tree_1, tree_2)
 
par(mfrow=c(1,2))
plot(rooted_tree, type = "phylogram", edge.width = 2)
add.scale.bar(cex = 1, font = 2, col = "red")
plot(double_tree, edge.width = 2)
add.scale.bar(cex = 1, font = 2, col = "red") 

#19) We can save our tree in the to main formats. 

write.tree(rooted_tree, file = "rooted_tree.newick") # newick format 
cat(readLines("rooted_tree.newick"), sep = "\n")

writeNexus(rooted_tree, "rooted_tree.nexus") #nexus format 
cat(readLines("rooted_tree.nexus"), sep = "\n")