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

【问题标题】：Counting phylogenetic tree topologies in R计算 R 中的系统发育树拓扑
【发布时间】：2020-06-26 15:59:46
【问题描述】：

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

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

我将拥有来自 15 种可能拓扑的 1000 棵树。将每种拓扑的树数制表的最简单方法是什么（即计数总和为 1000）。

【问题讨论】：

标签： r ape-phylo

【解决方案1】：

小树

对于较小的树（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)。

然后可以使用 user12728748 建议的边缘矩阵来比较树。

【讨论】：

【解决方案2】：

我不是这些树的专家，但在这里也许你可以对边缘矩阵进行排序，然后计算唯一的矩阵。这在本示例中在一定程度上有效，但是您可以翻转树的三种对称情况仍然单独计算（例如，第 8 和第 9 种情况，其中 ab 和 cd 与 cd 和 ab 是一个节点的子节点）。在下面的结果表中，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 年 6 月 26 日创建}

【讨论】：