#首先導(dǎo)入ape
if(!require(ape)){
install.packages("ape") #'ape' will be installed if needed
library(ape)
}
tree使用 “phylo” 類
我們使用 rtree() function 生成一個樹，參數(shù) n= 30個 tips (可以理解為末端節(jié)點)
phy <- rtree(n=30) #n specifies the number of tips we want.
plot(phy)
查看樹是由哪些數(shù)據(jù)構(gòu)成的，我們先來理解下“phylo”類
class(phy) #查看phy類型
## [1] "phylo"
phy #包含30個tips 和29個內(nèi)部節(jié)點（可以理解為樹的分叉點）
## Phylogenetic tree with 30 tips and 29 internal nodes.
##
## Tip labels:
## t13, t28, t8, t22, t19, t16, ...
##
## Rooted; includes branch lengths.
str(phy) #其實phy類似于list類型
## List of 4
## $ edge : int [1:58, 1:2] 31 32 33 34 35 36 37 37 38 38 ...
## $ tip.label : chr [1:30] "t13" "t28" "t8" "t22" ...
## $ edge.length: num [1:58] 0.2235 0.4492 0.0535 0.9588 0.2777 ...
## $ Nnode : int 29
## - attr(*, "class")= chr "phylo"
## - attr(*, "order")= chr "cladewise"
依次來看一下各數(shù)據(jù)結(jié)構(gòu)及代表什么
head(phy$edge) #一個兩列的數(shù)據(jù)框
## [,1] [,2]
## [1,] 31 32
## [2,] 32 33
## [3,] 33 34
## ..
phy$tip.label #樹末端的label
## [1] "t13" "t28" "t8" "t22" "t19" "t16" "t7" "t11" "t4" "t24" "t17"
## [12] "t26" "t15" "t1" "t12" "t30" "t27" "t6" "t29" "t3" "t20" "t25"
## [23] "t10" "t21" "t23" "t14" "t18" "t9" "t5" "t2"
phy$edge.length #樹枝的長度
## [1] 0.223515275 0.449214761 0.053517260 0.958796252 0.277726818
## ..
## [46] 0.339903616 0.173343718 0.993939807 0.176320243 0.087404924
## [51] 0.113624431 0.108095456 0.150655877 0.317285047 0.572281965
## [56] 0.019839419 0.529848145 0.998013315
phy$Nnode #中間節(jié)點的個數(shù)
## [1] 29
Newick格式的tree文件
我們構(gòu)建一個tips為 A-F的樹，使用 read.tree()讀取
mini.phy <- read.tree(text = "((((A,B), C), (D,E)),F);") #defining a tree in bracket form
plot(mini.phy, cex=2)
這個tree只包含3個數(shù)據(jù)，因為我們沒提供edge的長度
str(mini.phy)
## List of 3
## $ edge : int [1:10, 1:2] 7 8 9 10 10 9 8 11 11 7 ...
## $ tip.label: chr [1:6] "A" "B" "C" "D" ...
## $ Nnode : int 5
## - attr(*, "class")= chr "phylo"
## - attr(*, "order")= chr "cladewise"
我們看看這些邊（edges）是怎么表示的
mini.phy$edge
## [,1] [,2]
## [1,] 7 8
## [2,] 8 9
## [3,] 9 10
## [4,] 10 1
## [5,] 10 2
## [6,] 9 3
## [7,] 8 11
## [8,] 11 4
## [9,] 11 5
## [10,] 7 6
我們看到一個矩陣有10行和2列。這個矩陣表示節(jié)點和下一節(jié)點組合，定義了樹的每一個分支。下面圖中標(biāo)出每個節(jié)點的數(shù)字
plot(mini.phy, label.offset=0.2)
nodelabels() #add node numbers
tiplabels() #add tip numbers
使用tiplabels()給tip加上不同的顏色
mycol<-c("blue", "blue", "blue", "red", "red", "red") #顏色向量
plot(mini.phy, adj=0, label.offset=0.75) #adjust if necessary/desired!
tiplabels(pch=21, col="black", adj=1, bg=mycol, cex=2) #ditto! try different values!
交換樹枝的位置
有時旋轉(zhuǎn)節(jié)點是很有必要的，使用rotate()function。繼續(xù)使用mini.phy數(shù)據(jù)
plot(mini.phy, label.offset=0.2) #same as before
nodelabels() #add node numbers
tiplabels() #add tip numbers
交換分支clades (D, E) and (A, B, C)? 圖中可以看出他們的最近共同祖先是node 8.
rot.phy <- rotate(mini.phy, node=8)
plot(rot.phy, label.offset=0.2) #still the same as before
nodelabels()
tiplabels()
找到指定分支下的所有tips，可以使用geiger包 的the tips() function得到分支下所有端點
library(geiger)
cladeABC <- tips(rot.phy, node=9) # node 9 defines the clade composed of (A, B, C)
cladeABC
## [1] "A" "B" "C"
另外還有個修剪枝的函數(shù)drop.tip()，下面是剪切掉cladeABC分支
pruned.phy <- drop.tip(rot.phy, cladeABC)
plot(pruned.phy, label.offset=0.2)
nodelabels()
tiplabels() #did it work?
It may also be useful to select all branches in a specific group of tips. This is implemented in theapepackage; the which.edge() function finds all edges in a specified group. For example, let’s identify all branches of the cladeABC group as defined above.
cladeABCbranches <- which.edge(rot.phy, cladeABC)
#cladeABC was defined earlier, using the tips() function
cladeABCbranches #this should be a numerical vector containing 6, 7, 8, 9
## [1] 6 7 8 9
And as we can see, rows 6-9 of the
edge.length matrix represent the expected branches. Let’s first plot the tree, and then look at the edge matrix for cross-checking.
plot(rot.phy, label.offset=0.2) #we are sticking to the same plot settings
nodelabels()
tiplabels()
rot.phy$edge
## [,1] [,2]
## [1,] 7 8
## [2,] 8 11
## [3,] 11 4
## [4,] 11 5
## [5,] 8 9
## [6,] 9 10
## [7,] 10 1
## [8,] 10 2
## [9,] 9 3
## [10,] 7 6
Here would be a good opportunity to show you how to assign different branch colors. For example, how can one emphasize the branches of the clade formed by A, B, and C? We first create a color vector, only consisting of grey colors. Then we’ll assign black to all branches of clade ABC.
clcolr <- rep("darkgrey", dim(rot.phy$edge)[1]) #defining a color vector; look up rep() and dim() functions!
clcolr[cladeABCbranches] <- "red" #specifying color of branches in clade (A, B, C)
lot(rot.phy, edge.color=clcolr, edge.width=2) #also making the branches a little thicker
Let’s conclude this section with one last exercise: combining trees. Assume you have two different phylogenies, with two different sets of taxa (no overlap). Another assumption is that you have knowledge how the trees may fit together. Then the bind.tree() function of ape package can help. The function takes in two phylo objects. The position of where the trees are bound is defined by tip- or node number within the first tree. Note that you can also specify the “root” as binding position.
#the first tree
tree1 <- rtree(n=10)
tree1$tip.label <- rep("apples", 10) #renaming labels to enhance readability
plot(tree1, label.offset=0.1); nodelabels(); tiplabels()
#the second tree
tree2 <- rtree(n=10)
tree2$tip.label <- rep("oranges", 10)
plot(tree2)
#merging trees
combined.tree <- bind.tree(tree1, tree2, where=5) #combining trees at position "5"
plot(combined.tree)
Try at least two different merging positions. For each example, plot the two original trees and the merged version in one figure.
Fully resolved and polytomous trees
All tree examples so far were fully resolved, i.e. each tree was fully binary. It’s very easy to access visually for small trees, but one can also do this more formally:
is.binary.tree(mini.phy)
## [1] TRUE
Let’s make a tree that is not fully resolved.
poly.phy <- read.tree(text = "(((A,B,C),(D,E)),F);") #A, B, and C form a polytomy
plot(poly.phy)
Is this tree binary?
is.binary.tree(poly.phy) # "FALSE"!
## [1] FALSE
OK. Many comparative methods require fully resolved trees. But what to do if that’s not the case? The multi2di() function can resolve polytomies, either randomly or in the order in which they appear in the tree. The default setting is to resolve polytomies randomly.
resolved.phy <- multi2di(poly.phy)
is.binary.tree(resolved.phy)
## [1] TRUE
plot(resolved.phy) # visual inspection.
If you repeat the above lines a few times you will see the effect of randomly resolving the polytomy.
Please plot two randomly resolved trees next to each other!
畫圖參數(shù)設(shè)置
通過修改參數(shù)，改變tree的樣式
phy <- rtree(n=10) # n specifies the number of tips we want.
plot(phy) #default orientation is right-handed
plot(phy, direction="upwards", # "rightwards" (default), "leftwards", and "downwards"
show.tip.label=FALSE, #顯示tip.label
#cex=.3 #tip.lable 字體大小
edge.width=2, #枝寬度
edge.color="blue") # 樹枝顏色 ；藍色
There are many other options… type ?plot.phylo in the console to read more. One last useful command I’d like to introduce is the ladderize() function. It reorganizes the tree structure, normally yielding much more readable trees.
ladderized.phy <- ladderize(phy)
plot(ladderized.phy, direction="upwards",
show.tip.label=FALSE, edge.width=2, edge.color="blue")
