Tengo dos listas de estructuras complejas (cada lista es un objeto multiPhylo que contiene árboles filogenéticos) y me gustaría saber cuántas veces aparece cada elemento del primero en el segundo. Bastante sencillo, pero por algunas razones mi código devuelve resultados incorrectos.

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

La función all.equal.phylo() toma dos objetos phylo y devuelve VERDADERO si son iguales. Consulte documentos. La función count() toma un elemento y una lista y devuelve el número de veces que el elemento aparece en la lista usando all.equal.phylo().

El problema es que la función count() devuelve 0 la mayor parte del tiempo. Esto no debería ser posible ya que la lista unique_multiphylo es una sublista de beast_output_rooted, lo que significa que count() debería devolver al menos 1.

¿Qué pasa con mi código? ¿Y cómo puedo corregirlo? Muchas gracias por la ayuda!

EDITAR: aquí hay un ejemplo reproducible:

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

Sin embargo, parece funcionar perfectamente bien con estos datos simulados ...

2
Justine Vandendorpe 4 feb. 2019 a las 23:17

2 respuestas

La mejor respuesta

Hay una solución más breve para su problema:

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

Ya que unique_multiphylo contiene un atributo con la información que busca (consulte ?unique.multiPhylo).

1
klash 13 feb. 2019 a las 19:43

Finalmente logré obtener los resultados adecuados. En la función all.equal.phylo(), necesitaba establecer el parámetro use.edge.length en FALSE para que solo se comparan las topologías de los árboles filogenéticos.

Aquí está mi código:

(Cambié los nombres de un par de variables para aclarar lo que estaba tratando de hacer)

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)
2
Justine Vandendorpe 7 feb. 2019 a las 10:13