r - 比较复杂结构列表

Question

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

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 对象，如果它们相同则返回 TRUE。请参阅docs 。函数count()接受一个项目和一个列表，并使用 all.equal.phylo() 返回该项目在列表中出现的次数。 .

问题是函数 count()大部分时间返回 0。这应该是不可能的，因为列表 unique_multiphylo是 beast_output_rooted 的子列表，这意味着count()至少应返回 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)