R/simulate_without_recombination_expensive.R

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

#' Simulates the evolution of a family of sequences by taking into
#' consideration both point mutation and burst events and saves the sequeces
#' @param sequences The original sequences to begin with
#' @param num_jumps The number of time steps to simulate for
#' @param transition_matrix The transition matrix for the point mutations
#' @param burst_threshhold The probability of a burst even occurring
#' @param burst_mean The mean (for a Poisson distribution) increase in copy
#' number during a burst event
#' @return A list of sequences that have evolved over time with both point#
#' mutation and recombination
simulate_without_recombination_expensive <- function(sequences, num_jumps, transition_matrix,
                             burst_threshold, burst_mean)
{
  sequence_length <- length(sequences[[1]])
  data <- data.frame(time = rep(NA_integer_, num_jumps),
                     sequences = rep(NA, num_jumps))

  data[1, ] <- c(1, list(list(sequences)))
  for(t in 2:num_jumps)
  {
    sequences <- mutate_sequences(sequences, transition_matrix)
    num_sequences <- length(sequences)

    burst_numbers <- burst_sequences(burst_mean, num_sequences,
                                  to_burst(burst_threshold, num_sequences))

    new_num_sequences <- sum(burst_numbers)
    if(new_num_sequences > num_sequences)
    {
      create_new_sequence <- function(index, sequences, burst_numbers)
      {
        seq1 <- sequences[[select_sequence(burst_numbers)]]
        seq2 <- sequences[[select_sequence(burst_numbers)]]
        new_sequence <- recombine(list(seq1, seq2), 0)
        return(new_sequence)
      }
      new_sequences <-  purrr::map((num_sequences+1):new_num_sequences,
                                   create_new_sequence,
                                   sequences,
                                   burst_numbers)
      burst_mean <- update_burst_mean(burst_mean, num_sequences, new_num_sequences)
      sequences <- c(sequences, new_sequences)
    }
    new_row <- c(t, list(list(sequences)))
    data[t, ] <- new_row
  }
  return(data)
}