R/copy_number_evolution.R

In sams25/rcombinator_old: Tools to Simulate the Evolution of Genetic Sequences

#' Simulates the copy number evolution for a set of transposons (very basic function)
#' @param timestep The size of the timesteps of jumps
#' @param final_time The time to simulate until
#' @param init_num_sequences The initial number of sequences we have
#' @param burst_function To convert from copy numbers before burst to copy numbers after burst
#' @param kill_function To saturate a copy number at some given value
#' @return A data frame with time and copy number for all timesteps
#' @export
copy_number_evolution <- function(
                               timestep,
                               final_time,
                               init_num_sequences,
                               burst_function,
                               burst_function_args,
                               kill_function,
                               kill_function_args)
{
  num_jumps = ceiling(final_time/timestep) + 1

  times <- timestep*(0:(num_jumps-1))
  num_copies <- data.frame(time = times)

  # the copy numbers at each time
  sequences <- rep(list(c()), num_jumps)
  cur_seqs <- 1:init_num_sequences
  cur_copy_nums <- rep(1, length(cur_seqs))

  sequences[[1]] <- cur_seqs

  for (iter in 2:num_jumps) {
    # burst and kill
    cur_copy_nums <- do.call(burst_function, c(list(copy_numbers=cur_copy_nums),
                                               burst_function_args))
    cur_copy_nums <- do.call(kill_function, c(list(copy_numbers=cur_copy_nums),
                                              kill_function_args))

    sequences[[iter]] <- rep(cur_seqs, cur_copy_nums)
    cur_seqs <- 1:sum(cur_copy_nums)
    cur_copy_nums <- rep(1, length(cur_seqs))
  }

  num_copies$sequences <- sequences
  return(num_copies)
}