R/retrotransposon_evolution.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
#' @param timestep The size of the timesteps of jumps
#' @param final_time The time to simulate until
#' @param initial_sequences The original sequences to begin with
#' @param mutation_function How to apply point mutations
#' @param copy_numbers The copy numbers for the sequences over time
#' @param recombination_function How to recombine sequences
#' @param dist_function How to recombine sequences
#' @return A list of sequences that have evolved over time with both point
#' mutation and recombination
#' @export
retrotransposon_evolution <- function(
                               timestep,
                               final_time,
                               initial_sequences,
                               mutation_function,
                               copy_numbers,
                               recombination_function,
                               dist_function,
                               num_seq_dumps,
                               num_summary_dumps,
                               num_pairmatrix_dumps)
{
  time <- copy_numbers$time
  num_jumps <- length(time)

  sequences <- initial_sequences
  initial_sequence <- sequences[[1]]

  seq_dump <- create_initial_seq_dump(time, sequences, num_seq_dumps)
  seq_dump_indices <- seq_dump$seq_dump_indices
  seq_dump_data <- seq_dump$seq_dump_data
  seq_dump_data_index <- seq_dump$seq_dump_data_index

  summary_dump <- create_initial_summary_dump(time, dist_function, copy_numbers,
                                              sequences, initial_sequence, num_summary_dumps)
  summary_dump_indices <- summary_dump$summary_dump_indices
  summary_dump_data <- summary_dump$summary_dump_data
  summary_dump_data_index <- summary_dump$summary_dump_data_index

  pairmatrix_dump <- create_initial_pairmatrix_dump(time, dist_function,
                                                    sequences, num_pairmatrix_dumps)
  pairmatrix_dump_indices <- pairmatrix_dump$pairmatrix_dump_indices
  pairmatrix_dump_data <- pairmatrix_dump$pairmatrix_dump_data
  pairmatrix_dump_data_index <- pairmatrix_dump$pairmatrix_dump_data_index

  new_sequences <- get_random_sequence(length(initial_sequence), 500)
  for (iter in 2:num_jumps) {
    cat(paste(iter, "/", num_jumps, "\n", sep=""))

    sequences <- mutation_function(sequences, timestep)
    new_seq_indices <- copy_numbers$sequences[[iter]]

    if (length(new_seq_indices) == 0) {
      print("Reached 0 active sequences")
      break
    }

    to_recombine <- duplicated(new_seq_indices)
    for (seq_index in 1:length(new_seq_indices)) {

      if (!to_recombine[seq_index]) {
        new_sequences[[seq_index]] <- sequences[[new_seq_indices[seq_index]]]
      }
      else {
        new_sequences[[seq_index]] <- recombination_function(
                                        sequences[[new_seq_indices[seq_index]]],
                                        sequences[[sample(new_seq_indices, size=1)]]
                                      )
      }
    }
    sequences <- new_sequences[1:length(new_seq_indices)]

    if (seq_dump_indices[iter]) {
      seq_dump_data$sequences[[seq_dump_data_index]] <- sequences
      seq_dump_data_index <- seq_dump_data_index + 1
    }
    if (pairmatrix_dump_indices[iter]) {
      pairmatrix_dump_data$pairmatrix[[pairmatrix_dump_data_index]] <- dist_function(sequences)
      pairmatrix_dump_data_index <- pairmatrix_dump_data_index + 1
    }
    if (summary_dump_indices[iter]) {
      summary_dump_data_init_dists <- sapply(sequences,
                                             function(x)
                                               dist_function(list(x, initial_sequence))[1])
      quantile_init_stats <- stats::quantile(summary_dump_data_init_dists, na.rm = TRUE)
      summary_dump_data$avg_init_dist[summary_dump_data_index] <- mean(summary_dump_data_init_dists, na.rm = TRUE)
      summary_dump_data$min_init_dist[summary_dump_data_index] <- quantile_init_stats["0%"]
      summary_dump_data$q25_init_dist[summary_dump_data_index] <- quantile_init_stats["25%"]
      summary_dump_data$med_init_dist[summary_dump_data_index] <- quantile_init_stats["50%"]
      summary_dump_data$q75_init_dist[summary_dump_data_index] <- quantile_init_stats["75%"]
      summary_dump_data$max_init_dist[summary_dump_data_index] <- quantile_init_stats["100%"]

      summary_dump_data_index <- summary_dump_data_index + 1
    }
  }

  return(list(initial_sequence=initial_sequence,
              seq_dump_data=seq_dump_data,
              summary_dump_data=summary_dump_data,
              pairmatrix_dump_data=pairmatrix_dump_data))
}