比较复杂结构列表

Question

我有两个复杂结构列表（每个列表都是一个包含系统发育树的 multiPhylo 对象），我想找出第一个的每个元素在第二个中出现了多少次。非常简单，但由于某些原因，我的代码 return 的结果是错误的。

library(devtools)
install_github('santiagosnchez/rBt')
library(rBt)

beast_output <- read.annot.beast('strict_BD_starBEAST_logcomb.species.trees')
beast_output_rooted <- root(beast_output, c('taxon_A', 'taxon_B')) # length == 20,000
unique_multiphylo <- unique.multiPhylo(beast_output_rooted) # length == 130

count <- function(item, list) {
  total = 0
  for (i in 1:length(list)) {
    if (all.equal.phylo(item, list[[i]])) {
      total = total + 1
    }
  }
  return(total)
}

result <- data.frame(un_tr = rep(0, 130), count = rep(0, 130))
for (i in 1:length(unique_multiphylo)) {
  result[i, ] <- c(i, count(unique_multiphylo[[i]], beast_output_rooted))
}

函数 all.equal.phylo() 接受两个 phylo 对象，如果它们相同，则 return 为 TRUE。参见 docs。函数 count() 接受一个项目和一个列表，并且 returns 使用 all.equal.phylo().

计算项目在列表中出现的次数

问题是函数 count() return 大多数时候都是 0。这应该是不可能的，因为列表 unique_multiphylo 是 beast_output_rooted 的子列表，这意味着 count() 应该至少 return 1.

我的代码有什么问题？我该如何纠正它？非常感谢您的帮助！

编辑：这是一个可重现的例子：

install.packages('ape')
library(ape)

set.seed(42)
trees <- lapply(rep(c(10, 25, 50, 100), 3), rtree) # length == 12
class(trees) <- 'multiPhylo'
unique_multiphylo <- unique.multiPhylo(trees) # length == 12

count <- function(item, list) {
  total = 0
  for (i in 1:length(list)) {
    if (all.equal.phylo(item, list[[i]])) {
      total = total + 1
    }
  }
  return(total)
}

result <- data.frame(un_tr = rep(0, 12), count = rep(0, 12))
for (i in 1:length(unique_multiphylo)) {
  result[i, ] <- c(i, count(unique_multiphylo[[i]], trees))
}

但是，这些模拟数据似乎工作得很好...

Answer 1

我终于得到了正确的结果。在函数 all.equal.phylo() 中，我需要将参数 use.edge.length 设置为 FALSE，以便只比较系统发育树的 topologies。

这是我的代码：

（我更改了几个变量的名称以使其更清楚我正在尝试做什么）

install.packages('devtools')
library(devtools)
install_github('santiagosnchez/rBt')
library(rBt)

beast_output <- read.annot.beast('beast_output.trees')
beast_output_rooted <- root.multiPhylo(beast_output, c('taxon_A', 'taxon_B'))
unique_topologies <- unique.multiPhylo(beast_output_rooted)

count <- function(item, list) {
  total = 0
  for (i in 1:length(list)) {
    if (all.equal.phylo(item, list[[i]], use.edge.length = FALSE)) {
      total = total + 1
    }
  }
  return(total)
}

result <- data.frame(unique_topology = rep(0, length(unique_topologies)),
                     count = rep(0, length(unique_topologies)))
for (i in 1:length(unique_topologies)) {
  result[i, ] <- c(i, count(unique_topologies[[i]], beast_output_rooted))
}

result$percentage <- ((result$count/length(beast_output_rooted))*100)

Answer 2

您的问题有更短的解决方案：

table( attr(unique_multiphylo, "old.index") )

因为 unique_multiphylo 包含一个包含您所需要的信息的属性（参见 ?unique.multiPhylo）。

比较复杂结构列表

Compare list of complex structures

algorithm

comparison

r

phylogeny