compare_stats_families.R

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

library(rcombinator)
devtools::reload(".")
setwd("./output_for_comparison/")

full_pipeline <- function(sequence_length,
                          timestep, # in years
                          final_time,
                          burst_probability,
                          burst_mean,
                          max_copy_number,
                          recombination_rate,
                          dist_method,
                          num_matrix_dumps,
                          num_seq_dumps)
{
    x <- get_random_sequence(sequence_length, 2)
    filename <- paste(sequence_length, timestep, final_time, burst_probability,
                      burst_mean, max_copy_number, recombination_rate,
                      dist_method, sep="_")
    title <- paste0("Sequence distances over time with",
                    "\nsequence length ", sequence_length,
                    ",\ntime jumps of ", timestep, " years",
                    ",\nfor ", final_time, " years",
                    ",\nwith burst probability ", burst_probability,
                    ",\nwith burst mean ", burst_mean,
                    ",\nwith maximum copy number ", max_copy_number,
                    ",\nand recombination rate ", recombination_rate,
                    ",\nusing distance metric ", dist_method)

    data <- simulate_burst_evolution(x,
                                   timestep,
                                   final_time,
                                   get_k80_matrix(timestep),
                                   burst_probability,
                                   burst_mean,
                                   max_copy_number,
                                   recombination_rate,
                                   dist_method,
                                   num_matrix_dumps,
                                   num_seq_dumps)
    save(data, file = paste0(filename, ".RData"))
    summary <- data[[2]]
    title <- paste(title, "\nFinal mean, median, quant25, quant75:",
                  tail(summary$mean_dist, 1),
                  tail(summary$median_dist, 1),
                  tail(summary$quantiles_dist, 1)[[1]]["25%"],
                  tail(summary$quantiles_dist, 1)[[1]]["75%"],
                  "\n")
    p <- plot_summary_stats_complex(summary, title, timestep)
    ggplot2::ggsave(paste0(filename, ".png"), width=12, height=14, plot=p)
}

sequence_length <- 5000
timestep_v <- c(10^6, 10^5)
final_time <- 500*10^6
burst_probability_v <- c(0.01, 0.1, 1)
burst_mean <- 1
max_copy_number <- 5
recombination_rate <- 0.3
dist_method <- "ape"
num_matrix_dumps <- 100
num_seq_dumps <- 0

for(timestep in timestep_v)
{
    for(burst_probability in burst_probability_v)
    {
        full_pipeline(sequence_length,
                      timestep, final_time,
                      burst_probability, burst_mean,
                      max_copy_number,
                      recombination_rate,
                      dist_method,
                      num_matrix_dumps,
                      num_seq_dumps)

        full_pipeline(sequence_length,
                      timestep, final_time,
                      burst_probability, burst_mean,
                      max_copy_number,
                      0,
                      dist_method,
                      num_matrix_dumps,
                      num_seq_dumps)
    }
}