R/simulate_bursts.R

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

#' Simulate the increase in copy numbers of sequences
#' @param sequence_counts The distribution of copy numbers of the sequences
#' @param num_jumps The number of time steps this simulation goes across
#' @param burst_probability The probability of a burst event occurring
#' @param burst_choice A choice of how to select the sequence to burst
#' "parameters" - takes the probabilities assigned to sequence_counts
#' "equal"      - gives each family the same weightage
#' "accumulate" - weighs a sequence in proportion to its copy number
#' @return A data frame of sequence counts with respect to time
#' @export
simulate_bursts <- function(sequence_counts, num_jumps, burst_threshhold, burst_choice)
{
  num_sequences <- length(sequence_counts$distribution)

  # the sequence counts with respect to time
  data <- data.frame(matrix(0, ncol=num_sequences+1, nrow=num_jumps))
  colnames(data) <- c("time", paste(as.character(1:num_sequences), "count"))

  for(t in 1:num_jumps)
  {
    data[t, ] <- c(t, sequence_counts$distribution)
    if(to_burst(burst_threshhold))
    {
      if(burst_choice == "equal")
      {
        i <- select_sequence(rep(1, num_sequences))
      }
      else if(burst_choice == "accumulate")
      {
        i <- select_sequence(sequence_counts$distribution)
      }
      else if(burst_choice == "parameters")
      {
        i <- select_sequence(sequence_counts$burst_probs)
      }
      else
      {
        stop("invalid selection for burst choice")
      }
      sequence_counts$distribution[i] <-
        sequence_counts$distribution[i] +
        burst_sequences(sequence_counts$burst_means[i])
    }
  }
  return(data)
}