compare_stats_global.R

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

library(rcombinator)
devtools::reload(".")

full_pipeline <- function(sequence_length,
                          timestep, # in years
                          num_jumps,
                          burst_threshold,
                          burst_mean,
                          flip_probability,
                          max_copy_number,
                          dist_method,
                          recombination)
{
    x <- get_random_sequence(sequence_length, 2)
    filename <- paste(sequence_length, timestep, num_jumps, burst_threshold,
                      max_copy_number, dist_method, sep="_")
    title <- paste0("Sequence distances over time with\n sequence length ", sequence_length,
                    ",\n time jumps of ", timestep, " years",
                    ",\n for ", num_jumps, " jumps",
                    ",\n with burst probability ", burst_threshold,
                    ",\n with maximum copy number ", max_copy_number,
                    ",\n using distance metric ", dist_method)

    if(recombination)
    {
        summary <- simulate_recombination(x,
                                          num_jumps=num_jumps,
                                          transition_matrix=get_k80_matrix(timestep),
                                          burst_threshold=burst_threshold,
                                          burst_mean=burst_mean,
                                          flip_probability=flip_probability,
                                          max_copy_number=max_copy_number,
                                          dist_method=dist_method,
                                          options="summary")
        filename <- paste0(filename,"_recomb")
        title <- paste0(title, " WITH recombination\n")
    }
    else
    {
        summary <- simulate_without_recombination(x,
                                          num_jumps=num_jumps,
                                          transition_matrix=get_k80_matrix(timestep),
                                          burst_threshold=burst_threshold,
                                          burst_mean=burst_mean,
                                          max_copy_number=max_copy_number,
                                          dist_method=dist_method,
                                          options="summary")
        filename <- paste0(filename,"_norecomb")
        title <- paste0(title, " WITHOUT recombination\n")
    }
    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")

    filename <- paste0(filename, ".png")
    p <- plot_summary_stats_complex(summary, title, timejump = timestep)
    ggplot2::ggsave(filename, width=12, height=14, plot=p)
}

full_pipeline(500, timestep=10^4, num_jumps=800, burst_threshold=0.9, burst_mean=1,
              flip_probability = 0.3, max_copy_number=50, dist_method="ape", recombination=TRUE)
full_pipeline(500, timestep=10^4, num_jumps=800, burst_threshold=0.9, burst_mean=1,
              flip_probability = 0.3, max_copy_number=50, dist_method="ape", recombination=FALSE)

full_pipeline(500, timestep=10^5, num_jumps=800, burst_threshold=0.9, burst_mean=1,
              flip_probability = 0.3, max_copy_number=50, dist_method="ape", recombination=TRUE)
full_pipeline(500, timestep=10^5, num_jumps=800, burst_threshold=0.9, burst_mean=1,
              flip_probability = 0.3, max_copy_number=50, dist_method="ape", recombination=FALSE)

full_pipeline(500, timestep=10^6, num_jumps=800, burst_threshold=0.9, burst_mean=1,
              flip_probability = 0.3, max_copy_number=50, dist_method="ape", recombination=TRUE)
full_pipeline(500, timestep=10^6, num_jumps=800, burst_threshold=0.9, burst_mean=1,
              flip_probability = 0.3, max_copy_number=50, dist_method="ape", recombination=FALSE)

full_pipeline(5000, timestep=10^4, num_jumps=800, burst_threshold=0.9, burst_mean=1,
              flip_probability = 0.3, max_copy_number=50, dist_method="ape", recombination=TRUE)
full_pipeline(5000, timestep=10^4, num_jumps=800, burst_threshold=0.9, burst_mean=1,
              flip_probability = 0.3, max_copy_number=50, dist_method="ape", recombination=FALSE)

full_pipeline(5000, timestep=10^5, num_jumps=800, burst_threshold=0.9, burst_mean=1,
              flip_probability = 0.3, max_copy_number=50, dist_method="ape", recombination=TRUE)
full_pipeline(5000, timestep=10^5, num_jumps=800, burst_threshold=0.9, burst_mean=1,
              flip_probability = 0.3, max_copy_number=50, dist_method="ape", recombination=FALSE)

full_pipeline(5000, timestep=10^6, num_jumps=800, burst_threshold=0.9, burst_mean=1,
              flip_probability = 0.3, max_copy_number=50, dist_method="ape", recombination=TRUE)
full_pipeline(5000, timestep=10^6, num_jumps=800, burst_threshold=0.9, burst_mean=1,
              flip_probability = 0.3, max_copy_number=50, dist_method="ape", recombination=FALSE)