计算 R 中的系统发育树拓扑

Question

给定 R 中的 multiPhylo 对象，计算重复拓扑数的最简单方法是什么。例如，如果我从 4 tip 拓扑的所有 15 种可能分辨率中随机抽样：

library(ape)
library(phytools)
m <- do.call(c, lapply(1:1000, function(x) multi2di(starTree(c('a','b','c','d')))))

我将从 15 种可能的拓扑中得到 1000 棵树。将每个拓扑的树数制成表格的最简单方法是什么（即计数总和将为 1000）。

Answer 1

我不是那些树的专家，但在这里也许你可以对边矩阵进行排序，然后计算唯一的边矩阵。这在一定程度上适用于此示例，但您可以翻转树的三个对称情况仍然单独计算（例如第 8 和第 9 个，其中 ab 和 cd 与 cd 和 ab 是一个节点的子节点）。在下面的结果 table 中，N 将是计数，而 first 将是此拓扑在 m 中的第一次出现。您可以通过 for(i in res$first) plot(m[[i]]) 绘制独特的树并检查那些对称的情况。

library(ape)
library(phytools)
#> Loading required package: maps
library(data.table)
set.seed(123)
m <- do.call(c, lapply(1:1000, function(x) multi2di(starTree(c('a','b','c','d')))))
edges <- data.table(do.call(rbind, lapply(m, function(x) unlist(data.table(x$edge, key=c("V1", "V2"))))))
res <- edges[,.(.N, first=head(.I, 1)), by = names(edges)]
res
#>     V11 V12 V13 V14 V15 V16 V21 V22 V23 V24 V25 V26  N first
#>  1:   5   5   6   6   7   7   2   6   1   7   3   4 53     1
#>  2:   5   5   6   6   7   7   4   6   2   7   1   3 65     2
#>  3:   5   5   6   6   7   7   4   6   3   7   1   2 60     3
#>  4:   5   5   6   6   7   7   2   6   4   7   1   3 56     4
#>  5:   5   5   6   6   7   7   6   7   1   4   2   3 63     5
#>  6:   5   5   6   6   7   7   6   7   2   3   1   4 55     6
#>  7:   5   5   6   6   7   7   1   6   3   7   2   4 66     7
#>  8:   5   5   6   6   7   7   6   7   3   4   1   2 55     8
#>  9:   5   5   6   6   7   7   6   7   1   2   3   4 40    15
#> 10:   5   5   6   6   7   7   6   7   1   3   2   4 66    16
#> 11:   5   5   6   6   7   7   3   6   2   7   1   4 52    20
#> 12:   5   5   6   6   7   7   1   6   4   7   2   3 46    24
#> 13:   5   5   6   6   7   7   3   6   1   7   2   4 44    27
#> 14:   5   5   6   6   7   7   1   6   2   7   3   4 54    28
#> 15:   5   5   6   6   7   7   2   6   3   7   1   4 60    29
#> 16:   5   5   6   6   7   7   4   6   1   7   2   3 53    32
#> 17:   5   5   6   6   7   7   6   7   2   4   1   3 63    39
#> 18:   5   5   6   6   7   7   3   6   4   7   1   2 49    43

^{由 reprex package (v0.3.0)}

于 2020-06-26 创建

Answer 2

小树

对于较小的树（< ~20 个叶子），您可以使用“TreeTools”包将每个树拓扑转换为唯一的整数：

library('TreeTools')
library('phytools')
m <- do.call(c, lapply(1:1000, function(x) multi2di(starTree(c('a','b','c','d')))))

# Tabulate unique topologies
table(vapply64(m, as.TreeNumber, 1))

您可以使用

绘制每个编号的拓扑

topologyToPlot <- 2
plot(as.phylo(topologyToPlot, nTip = 4))

大树

对于较大的树，您可以确保具有等效拓扑结构的树在 R 中表示相同：

（如有必要）确保树的提示内部表示是一致的 m <- RenumberTips(m, m[[1]]).
reordering 树的内部边和节点编号使用 m <- Preorder(m).

然后可以按照的建议使用边矩阵比较树。

计算 R 中的系统发育树拓扑

Counting phylogenetic tree topologies in R

r

ape-phylo

小树

大树