Сравнить список сложных конструкций - PullRequest
Новая
       5

Сравнить список сложных конструкций

0 голосов
/ 04 февраля 2019

У меня есть два списка сложных структур (каждый список представляет собой объект 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() принимает два филологических объекта и возвращает ИСТИНА, если они совпадают.См. документы .Функция 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))
}

Однако, похоже, что он отлично работает с этими смоделированными данными ...

Ответы [ 2 ]

0 голосов
/ 13 февраля 2019

Существует более короткое решение вашей проблемы: 

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

, поскольку unique_multiphylo содержит атрибут с информацией, которую вы ищете (см. ?unique.multiPhylo).

0 голосов
/ 07 февраля 2019

Мне наконец удалось получить правильные результаты.В функции all.equal.phylo() мне нужно было установить для параметра use.edge.length значение FALSE, чтобы сравнивались только топологии филогенетических деревьев.

Вот мой код:

(я изменил имена нескольких переменных, чтобы было понятнее, что я пытался сделать) 

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)
...