In jaishimuku/WrightLab_FBD: Expanded Functions for Use with FossilSim

foss_data <- read.table("mukunda_foss.log", sep = '\t',  head = TRUE)

n.ages <- function(tree){
  depth = ape::node.depth.edgelength(tree)
  node.ages = max(depth) - depth
  names(node.ages) <- 1:(tree$Nnode+length(tree$tip))
  # adding possible offset if tree fully extinct
  if(!is.null(tree$root.time)) node.ages = node.ages + tree$root.time - max(node.ages)
  return(node.ages)
}

sfbdrt <- function (n, numbsim, lambda, mu, psi, times, complete = FALSE) 
{
  if (length(lambda) != length(times)) 
    stop("Length mismatch between rate shift times and birth rates")
  if (length(mu) != length(times)) 
    stop("Length mismatch between rate shift times and death rates")
  if (length(psi) != length(times)) 
    stop("Length mismatch between rate shift times and sampling rates")
  trees = TreeSim::sim.rateshift.taxa(n, numbsim, lambda, 
                                      mu, rep(1, length(times)), times, complete = TRUE)
  for (i in 1:length(trees)) {
    t = trees[[i]]
    origin = max(n.ages(t)) + t$root.edge
    horizons = c(times, origin)
    f <- sim.fossils.intervals(tree = t, interval.ages = horizons, 
                               rates = psi)
    tree = SAtree.from.fossils(t, f)
    node.ages = n.ages(tree)
    if (complete == FALSE) {
      fossil.tips = is.extinct(tree, tol = 1e-06)
      sa.tips = tree$tip.label[tree$edge[, 2][(tree$edge[, 
                                                         2] %in% 1:length(tree$tip.label)) & (tree$edge.length == 
                                                                                                0)]]
      unsampled.tips = fossil.tips[!(fossil.tips %in% 
                                       sa.tips)]
      tree = ape::drop.tip(tree, unsampled.tips)
      node.ages = n.ages(tree)
    }
    trees[[i]] = tree
    trees[[i]]$root.edge = origin - max(node.ages)
    trees[[i]] = SAtree(trees[[i]], complete)
  }

  f$h <- (f$hmin + f$hmax)/2

  return(list(f, trees))
}

library(FossilSim)
library(tidyverse)
foss_sample <- function(n_extant, n_trees, data_frame, mu, lambda, psi, row, mers)
{

  max_num <- max(c(length(mu), length(lambda), length(psi)))

  if(length(mu) == 1) 
        {  
          mi <- rep(data_frame[row, mu], times = max_num)
      } else {
          mi <- c()
          for (i in mu) {
            mi <- append(mi, data_frame[row,i])
          }
      }

  if(length(lambda) == 1) 
        {  
          li <- rep(data_frame[row, lambda], times = max_num)
      } else {
          li <- c()
          for (i in lambda) {
            li <- append(li, data_frame[row,i])
          }
      }

  if(length(psi) == 1) 
        {  
          pi <- rep(data_frame[row, psi], times = max_num)
      } else {
          pi <- c()
          for (i in psi) {
            pi <- append(pi, data_frame[row,i])
          }
      }

  my_fuction <- sfbdrt(n_extant, n_trees, li, mi, pi, times = mers, complete = TRUE)
   return(my_fuction)
}

my_summary <- function(phylo_list, interval_ages)
{
  tree_origin <-tree.max(phylo_list[[2]][[1]])
  num_fossils <- count.fossils(phylo_list[[1]])
  num_fossils_binned <- count.fossils.binned(phylo_list[[1]], interval_ages)

  return(list(tree_origin, num_fossils, num_fossils_binned))
}

Example list generation by SFB.

my_tree <- foss_sample(666, 1, foss_data, "mu", "lambda", c("psi.1.", "psi.2.", "psi.3."), 120, c(0, 29, 70))
View(my_tree)
print(my_tree)

Time And Number

print("Origin Time")
tree.max(my_tree[[2]][[1]]) 
print("Number/s")
count.fossils(my_tree[[1]]) 

summary_table <- my_summary(my_tree, c(0,29,70))
print(summary_table)

jaishimuku/WrightLab_FBD documentation built on July 17, 2020, 12:05 a.m.