夜雨聆风
R包神器 | ape (三) 树的写和读,根和外群
在ape包中，将树写出、再读入之后，树对象”phylo”的序列标签的顺序会发生变化 (具体变化见下面代码的文字注释)。
这种变化反而可能提示哪些序列可作为根(Root)或外群 (OutGroup)。当然，你也可以根据进化树的信息或相关文献，来自己指定根与外群。
ape中有关根和外群的函数：
root函数 – 相对于指定的外群，或在节点 (Node)中指定的节点，对系统发育树重新定根(Reroot)。
unroot函数 – 在系统发育树中去掉根，或者不变 (若已经是无根树)。
扫码、联系客服老师报名，领取资料、上手分析
读取aln.fasta文件
library(ape)
# 读取序列
ex.dna.all <- ape::read.dna("./result/pairsnp/core.Temporal.aln.fasta", format = "fasta", as.character=F)
# 选项format："fasta", "sequential", "clustal", "interleaved"(序列名含空格时), 或任何明确的缩写
ex.dna = ex.dna.all将树写至文件
# 由序列的遗传距离生成1个邻接树
trw <- bionj(dist.dna(ex.dna, model = "raw"))
str(trw)## List of 4
## $ edge : int [1:41, 1:2] 23 25 26 37 37 26 25 27 28 32 ...
## $ edge.length: num [1:41] 6.78e-05 1.76e-04 1.95e-02 7.50e-01 1.42e-02 ...
## $ tip.label : chr [1:22] "ERR2245277" "ERR2245279" "ERR2245288" "ERR2245293" ...
## $ Nnode : int 20
## - attr(*, "class")= chr "phylo"
## - attr(*, "order")= chr "cladewise"# phy - an object of class "phylo" or "multiPhylo"
trw##
## Phylogenetic tree with 22 tips and 20 internal nodes.
##
## Tip labels:
## ERR2245277, ERR2245279, ERR2245288, ERR2245293, ERR2245294, ERR2245295, ...
##
## Unrooted; includes branch lengths.# Phylogenetic tree with 22 tips and 20 internal nodes.
# Tip labels:
# ERR2245277, ERR2245279, ERR2245288, ERR2245293, ERR2245294, ERR2245295, ...
# Unrooted; includes branch lengths.
plot(trw)
注意：此时，树的tip.label的顺序与树的来源(e.g.距离矩阵)中序列的顺序相同。
但是，后续把树写出来、再读入后，tip.label的顺序会发生变化 (进化树图的倒序，或者Newick文件的正序?)
trw$tip.label## [1] "ERR2245277" "ERR2245279" "ERR2245288" "ERR2245293" "ERR2245294"
## [6] "ERR2245295" "ERR2245344" "ERR2245378" "ERR2245380" "ERR2245386"
## [11] "ERR2245388" "ERR2245393" "ERR2245400" "ERR2245401" "ERR2245402"
## [16] "ERR2245406" "ERR2245411" "ERR2245418" "ERR2245420" "ERR2245425"
## [21] "ERR2512377" "Reference"# "ERR2245277" "ERR2245279" "ERR2245288" ... "ERR2245425" "ERR2512377" "Reference"
labels(ex.dna)## [1] "ERR2245277" "ERR2245279" "ERR2245288" "ERR2245293" "ERR2245294"
## [6] "ERR2245295" "ERR2245344" "ERR2245378" "ERR2245380" "ERR2245386"
## [11] "ERR2245388" "ERR2245393" "ERR2245400" "ERR2245401" "ERR2245402"
## [16] "ERR2245406" "ERR2245411" "ERR2245418" "ERR2245420" "ERR2245425"
## [21] "ERR2512377" "Reference"# "ERR2245277" "ERR2245279" "ERR2245288" ... "ERR2245425" "ERR2512377" "Reference"
# 将树写入文件 - Write Tree File in Parenthetic Format
write.tree(trw, file = "./tree.txt", append = FALSE, digits = 10, tree.names = FALSE)读取树
注意：把树写到文件、再读进来之后，tip.label的顺序发生了变化，排到前面的可能与其它成员距离较远 (提示其可作为根或外群？)
# 从文件读取树，获得类(Class)对象‘phylo’
tree.f <- read.tree('./tree.txt')
# This stores the phylogenetic tree is in an object named tree1 of class ‘phylo’
tree.f##
## Phylogenetic tree with 22 tips and 20 internal nodes.
##
## Tip labels:
## Reference, ERR2512377, ERR2245277, ERR2245425, ERR2245418, ERR2245411, ...
##
## Unrooted; includes branch lengths.# Phylogenetic tree with 22 tips and 20 internal nodes. (内部节点)
# Tip labels:
# Reference, ERR2512377, ERR2245277, ERR2245425, ERR2245418, ERR2245411, ...
# Unrooted; includes branch lengths. (无根，含枝长，无自举值？)
tree.f$tip.label## [1] "Reference" "ERR2512377" "ERR2245277" "ERR2245425" "ERR2245418"
## [6] "ERR2245411" "ERR2245344" "ERR2245420" "ERR2245401" "ERR2245406"
## [11] "ERR2245380" "ERR2245393" "ERR2245288" "ERR2245293" "ERR2245402"
## [16] "ERR2245386" "ERR2245294" "ERR2245400" "ERR2245378" "ERR2245388"
## [21] "ERR2245295" "ERR2245279"# "Reference" "ERR2512377" "ERR2245277" ... "ERR2245388" "ERR2245295" "ERR2245279"
plot(tree.f)
定义树的根(Root)或外群(OutGroup) 留意代码中的文字注释
# root/unroot函数：需要输入的对象的类是"phylo"或"multiPhylo"
# 常用参数：
# phy - (树变量) an object of class "phylo" or "multiPhylo".
# outgroup - (数字或字符，来指定外群) mode numeric or character specifying the new outgroup.
# node - (可选，Node的编号) alternatively, a node number where to root the tree.
# resolve.root - (是否解离新产生的根，成为二分叉节点) resolve the new root as a bifurcating node.
# interactive - (交互式选择节点) select the node by clicking on the tree which must be plotted.
# edgelabel - 是否将节点的标签视为边的标签，并切换，以便在使用drawSupportOnEdges时将其与正确的边相关联，见：Czech et al. 2016.
# is.rooted函数：逻辑判断
is.rooted(tree.f)## [1] FALSEtree.f$tip.label## [1] "Reference" "ERR2512377" "ERR2245277" "ERR2245425" "ERR2245418"
## [6] "ERR2245411" "ERR2245344" "ERR2245420" "ERR2245401" "ERR2245406"
## [11] "ERR2245380" "ERR2245393" "ERR2245288" "ERR2245293" "ERR2245402"
## [16] "ERR2245386" "ERR2245294" "ERR2245400" "ERR2245378" "ERR2245388"
## [21] "ERR2245295" "ERR2245279"tree.f.r = root(tree.f, outgroup = "Reference")
is.rooted(tree.f.r)## [1] FALSE# [1] FALSE
# 在重定根之前，统一转为了无根树 (Unrooted first before rerooting)
tree.f.r$tip.label## [1] "Reference" "ERR2512377" "ERR2245277" "ERR2245425" "ERR2245418"
## [6] "ERR2245411" "ERR2245344" "ERR2245420" "ERR2245401" "ERR2245406"
## [11] "ERR2245380" "ERR2245393" "ERR2245288" "ERR2245293" "ERR2245402"
## [16] "ERR2245386" "ERR2245294" "ERR2245400" "ERR2245378" "ERR2245388"
## [21] "ERR2245295" "ERR2245279"# 此时已经含有外群标记，删除这个外群后，即为有根树。因为删除之前，含有"基部三分叉(Basal trichotomy)"
is.rooted(drop.tip(tree.f.r, "Reference"))## [1] TRUE# 所以，将三分叉解离为二分叉，也能获得有根树
is.rooted(multi2di(tree.f.r))## [1] TRUE# 保留"基部三分叉(Basal trichotomy)"，但强制转为有根树
tree.f.r.t=tree.f.r
tree.f.r.t$root.edge <- 0
is.rooted(tree.f.r.t)## [1] TRUEplot(tree.f.r.t)
plot(tree.f.r)
# 二者图形一致。内容可能不一致？
# 也可以直接获得有根树 (加上选项`r = TRUE`)
tree.f.rr = root(tree.f, outgroup = "Reference", r = TRUE)
is.rooted(tree.f.rr)## [1] TRUEtree.f.rr$tip.label## [1] "Reference" "ERR2512377" "ERR2245277" "ERR2245425" "ERR2245418"
## [6] "ERR2245411" "ERR2245344" "ERR2245420" "ERR2245401" "ERR2245406"
## [11] "ERR2245380" "ERR2245393" "ERR2245288" "ERR2245293" "ERR2245402"
## [16] "ERR2245386" "ERR2245294" "ERR2245400" "ERR2245378" "ERR2245388"
## [21] "ERR2245295" "ERR2245279"tree.f.rr##
## Phylogenetic tree with 22 tips and 21 internal nodes.
##
## Tip labels:
## Reference, ERR2512377, ERR2245277, ERR2245425, ERR2245418, ERR2245411, ...
##
## Rooted; includes branch lengths.# Phylogenetic tree with 22 tips and 21 internal nodes.
# Tip labels:
# Reference, ERR2512377, ERR2245277, ERR2245425, ERR2245418, ERR2245411, ...
# Rooted; includes branch lengths.
plot(tree.f.rr) # 此时图中序列名的顺序，与tip.label变得不一致
同时定义2个及以上的Tips为根 (这2个Tips的进化关系/拓扑结构不会变)
tree.f.rr2 = root(tree.f, outgroup = c("Reference", "ERR2512377"), r = TRUE)
is.rooted(tree.f.rr2)## [1] TRUEtree.f.rr2$tip.label## [1] "Reference" "ERR2512377" "ERR2245277" "ERR2245425" "ERR2245418"
## [6] "ERR2245411" "ERR2245344" "ERR2245420" "ERR2245401" "ERR2245406"
## [11] "ERR2245380" "ERR2245393" "ERR2245288" "ERR2245293" "ERR2245402"
## [16] "ERR2245386" "ERR2245294" "ERR2245400" "ERR2245378" "ERR2245388"
## [21] "ERR2245295" "ERR2245279"tree.f.rr2##
## Phylogenetic tree with 22 tips and 21 internal nodes.
##
## Tip labels:
## Reference, ERR2512377, ERR2245277, ERR2245425, ERR2245418, ERR2245411, ...
##
## Rooted; includes branch lengths.# Phylogenetic tree with 22 tips and 21 internal nodes.
# Tip labels:
# Reference, ERR2512377, ERR2245277, ERR2245425, ERR2245418, ERR2245411, ...
# Rooted; includes branch lengths.
plot(tree.f.rr2)
上面所有的代码中，需要注意的地方：
小 节
参数outgroup可以是字符、数字。在第1种情况下，需给出了新外群顶端/Tips的标签；在第2种情况下，给出了向量phy$tip.label中标签的编号。
如果外群的长度为1 (即单个字符或数字)，则使用该顶端/Tip下的节点作为新的根，来重新对进化树定根。
如果外群的长度为2或2以上 (即指定2个或以上的序列为外群)，则使用内群(Ingroup)的最近共同祖先 (Most recent common ancestor, MRCA) 作为新根。注意，树在被重新定根之前显示为无根，因此如果outgroup已经是outgroup (运行过root函数，已标记了外群，但此时仍是无根树，因为含有基部三分叉，见上面的截图)，那么返回的树与原始树不同（参见示例代码）。如果外群不是单源的(Monophyletic，即树上有2个相同的序列名、且被指定为外群？)，则运行失败并提示错误。
如果设置了resolve.root=TRUE(即解离三分叉为二分叉)，root会在内群的MRCA的下方添加一个长度为0的分枝。
ape认为有根或无根的原理是：如果只有2个分枝连接到根，或者有root.edge元素 (即tree$root.edge不为NULL)，则认为树是有根的；在所有其它情况下，is.rooted函数都返回FALSE。
如果指定的节点已经是树的当前的根，则resolve.root=TRUE与node=，二者一起使用会产生错误。这是因为：在没有显式的/明确的(Explicit)外群的情况下，解离无根树中的节点时存在歧义。如果这个节点不是当前的根，则通过将节点右侧的分支/Clade (当树以默认方式绘制时)视为内群来任意地解决这个分歧。可参阅此处的详细说明：https://www.mail-archive.com/r-sig-phylo@r-project.org/msg03805.html
扫码、联系客服老师报名，领取资料、上手分析
