rloopsfor-loopphylogenyape-phylo

Selecting random phylogenetic trees from data set, runing analysis and trying to save results


I want to calculate the Moran index for several random phylogenetic trees from a data set.

As the result I want a list of: the values of Moran index calculation and the selected tree for that given calculation of Moran index (at each loop).

I need to know which tree was used to run Moran index.

I was able to get Moran values for the random trees, however I was not able to save in any way the random tree that were selected at each loop.

Since my data set is big, I created the example below:

library(ape)

x <- rmtree(20, 20)
names(x) <- paste("t", 1:10, sep = "")
var <- as.data.frame(matrix(rnorm(20),nrow=20))
nams <- data.frame(paste0("t",1:20,sep=""))
var_list <- cbind(nams,var)
colnames(var_list) <- c("names","variable")
var_list1 <- var_list$variable
names(var_list1) <- c(var_list[,1])

# Loop to select a random tree and run the Moran index

tr <- function(x, var_list1){
  resulist <- as.list(1:15) 
# treelist <- as.list(as.phylo(x[[i]])) #### not working
  for(i in 1:15){
    tutreer <- sample(x,size=1)[[1]]
    inr <- 1/cophenetic(tutreer)
    diag(inr) <- 0
    mran <- Moran.I(var_list1,inr,scaled = TRUE)
    mran
    
    resulist[[i]] <- list(obs=mran$observed,expect=mran$expected,
                            sd=mran$sd,pval=mran$p.value)
    # treelist[[i]] <- list(tutreer) #### Here, not working
    
  }
  return(resulist)
# return(treelist)
}

tr(x,var_list1)

Any ideas on how could I save the trees that was selected (and used to calculate the index) at each loop?


Solution

  • Why not just add it to the resulist[[i]] on each iteration?:

    resulist[[i]] <- list(obs=mran$observed,expect=mran$expected,
                              sd=mran$sd,pval=mran$p.value, tree=tutreer)
    

    Or, return the entire sampled tree, like this:

    tr <- function(x, var_list1){
      resulist <- as.list(1:15) 
      for(i in 1:15){
        sampled_tree <- sample(x,size=1)
        inr <- 1/cophenetic(sampled_tree[[1]])
        diag(inr) <- 0
        mran <- Moran.I(var_list1,inr,scaled = TRUE)
        resulist[[i]] <- c(mran, list(tree = sampled_tree))
      }
      return(resulist)
    }