R/simulate_recombination_summary.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
#' This function does not return the sequences over time, but only the summary stats
#' of the distances
#' @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_recombination_summary <- function(sequences, num_jumps, transition_matrix,
                                         burst_threshold, burst_mean, flip_probability,
                                         max_copy_number, dist_method)
{
  data_seq <- data.frame(time = rep(NA_integer_, 1),
                         sequences = rep(NA, 1))
  data_seq[1, ] <- c(1, list(list(sequences)))

  data <- get_summary_stats_timestep(pairwise_distances_timestamp(data_seq,
                                                                  1,
                                                                  dist_method),
                                      1)

  for(t in 2:num_jumps)
  {
    cat(paste(t, "/", num_jumps, "\n", sep=""))
    # point mutations
    sequences <- mutate_sequences(sequences, transition_matrix)
    num_sequences <- length(sequences)

    # check and see which members have a burst event
    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, flip_probability)
      {
        seq1 <- sequences[[select_sequence(burst_numbers)]]
        seq2 <- sequences[[select_sequence(burst_numbers)]]
        new_sequence <- recombine(list(seq1, seq2), flip_probability)
        return(new_sequence)
      }
      # generate a set of new sequences for each burst event
      new_sequences <-  purrr::map((num_sequences+1):new_num_sequences,
                                   create_new_sequence,
                                   sequences,
                                   burst_numbers, flip_probability)
      burst_mean <- update_burst_mean(burst_mean, num_sequences, new_num_sequences)
      sequences <- c(sequences, new_sequences)
      sequences <- kill_sequences(sequences, max_copy_number)
    }
    new_row <- c(t, list(list(sequences)))
    data_seq[1, ] <- new_row
    data_new <- pairwise_distances_bursts(data_seq, t, dist_method)
    data_new <- get_summary_stats_timestep(data_new, t)
    data <- rbind(data, data_new)
  }
  return(data)
}