R/full_pipeline_comparison.R

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

#' Full pipeline for simulation of retrotransposon evolution
#' @export
full_pipeline_comparison <- function(
                              sequence_length, num_initial_sequences,
                              timestep, final_time,
                              burst_probability, burst_mean, max_copy_number,
                              recombination_mean,
                              dist_method,
                              num_seq_dumps = 0, num_summary_dumps = 0, num_pairmatrix_dumps = 0)
{
  initial_sequences <- get_random_sequence(sequence_length, num_initial_sequences)
  filename <- paste(sequence_length, num_initial_sequences,
                    timestep, final_time,
                    burst_probability, burst_mean, max_copy_number,
                    recombination_mean,
                    dist_method, sep="_")
  filename_norecomb <- paste0(filename, "_norecomb.RData")
  filename_recomb <- paste0(filename, "_recomb.RData")

  if (dist_method == "ape") { dist_function <-  stringdist_ape }
  else if (dist_method == "levenshtein") { dist_function <-  stringdist_levenshtein }
  else if (dist_method == "biostrings") { dist_function <-  stringdist_biostrings }
  else if (dist_method == "biostrings_percentage") {
    dist_function <- function(x) stringdist_biostrings(x, similarity_percentage = TRUE)
  }
  else { stop("Choose a correct distance method") }

  copy_numbers <- copy_number_evolution(timestep,
                                        final_time,
                                        length(initial_sequences),
                                        burst_function=burst_copies,
                                        burst_function_args=
                                          list(burst_mean=burst_mean,
                                               burst_probability=burst_probability),
                                        kill_function=kill_copies,
                                        kill_function_args=
                                          list(max_copy_number=max_copy_number))
    data_norecomb <- retrotransposon_evolution(
                timestep, final_time,
                initial_sequences,
                mutation_function =
                  function(s, t) point_mutation_sequences(s, get_k80_matrix, t),
                copy_numbers,
                recombination_function =
                  function(a, b) recombine_sequences(list(a, b), 0),
                dist_function,
                num_seq_dumps,
                num_summary_dumps,
                num_pairmatrix_dumps)

    data_recomb <- retrotransposon_evolution(
                timestep, final_time,
                initial_sequences,
                mutation_function =
                  function(s, t) point_mutation_sequences(s, get_k80_matrix, t),
                copy_numbers,
                recombination_function =
                  function(a, b) recombine_sequences(list(a, b), recombination_mean),
                dist_function,
                num_seq_dumps,
                num_summary_dumps,
                num_pairmatrix_dumps)

    save(data_norecomb, file=filename_norecomb)
    save(data_recomb, file=filename_recomb)
    return(list(norecomb=filename_norecomb,
                recomb=filename_recomb))
}