cpp/analyse_output/analyse_output.R

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

# To analyse the output produced by a cpp simulation

# Global settings
HOME_DIR <- "~/SysBioProject/rcombinator/cpp/"
setwd(HOME_DIR)
source('analyse_output/header.R')
library(cowplot)

OUTPUT_LOC = "./unified_output/"
SAVE_LOC = "../../report/images/"

# Create filename from parameters
create_filename <- function(
                     timestep, final_time, num_runs,
                     init_num_seq, num_recomb, recomb_means,
                     num_sensitive, p_sensitive,
                     burst_mean, burst_p,
                     max_total_copies, max_active_copies,
                     sim_type)
{
  if (sim_type == NO_BURSTS)
  {
    filename = paste(timestep, final_time, num_runs,
                     init_num_seq,
                     sep="_")
  }
  else if (sim_type == WITHOUT_FLAGS)
  {
    if (burst_p == 1) {str_burst_p <- paste0(burst_p, ".0") }
    else { str_burst_p <- paste0(burst_p) }
    filename = paste(timestep, final_time, num_runs,
                     burst_mean, str_burst_p,
                     max_total_copies,
                     sep="_")
  }
  else if (sim_type == WITH_FLAGS)
  {
    if (burst_p == 1) {str_burst_p <- paste0(burst_p, ".0") }
    else { str_burst_p <- paste0(burst_p) }
    filename = paste(timestep, final_time, num_runs,
                     num_sensitive, p_sensitive,
                     burst_mean, str_burst_p,
                     sep="_")
  }
  else
  {
    print("NOT A VALID SIM TYPE")
  }
  filename = paste0(sim_type, "_", filename)
  return(filename)
}

# PARAMETER SCAN
source("./../R/copy_number_growth.R")
source("./../R/find_parameter_balance_copy_number_variation.R")
p <- find_paramater_balance_copy_number_variation()
save_plot(paste0(SAVE_LOC, "parameter_scan.png"),
          p, base_height = 12, base_width = 18)

# GENE OUTPUT BASIC POINT MUTATION MODEL
sim_type         <- NO_BURSTS
timestep         <- 1
final_time       <- 500
num_runs         <- 20
init_num_seq_v   <- c(10, 100, 1000)

source('analyse_output/header.R')
for(init_num_seq in init_num_seq_v)
{
  gene_f <- create_filename(
                        sim_type = sim_type,
                        timestep = timestep,
                        final_time = final_time,
                        num_runs = num_runs,
                        init_num_seq = init_num_seq
                        )
  gene <- get_data(gene_f, sim_type, OUTPUT_LOC)
  gene_p <-no_bursts_plot(gene)
  save_plot(paste0(SAVE_LOC,"gene_over_time_",init_num_seq,".png"),
            gene_p, base_height = 10, base_width = 13)
  # init in blue
  # pair in red
}

# FLAGGING NO BURSTS
sim_type         <- WITH_FLAGS
timestep         <- 1
final_time       <- 500
num_runs         <- 10
init_num_seq     <- 100
num_sensitive_v  <- c(2, 10, 50)
p_sensitive_v    <- c(1.0, 0.5, 0.1)

source('analyse_output/header.R')
for(num_sensitive in num_sensitive_v)
{
for(p_sensitive in p_sensitive_v)
{
  flags_no_bursts_f <- paste("flags_no_bursts",
                             timestep, final_time, num_runs,
                             init_num_seq, num_sensitive, p_sensitive,
                    sep="_")
  if (p_sensitive == 1)
  {
    flags_no_bursts_f <- paste0(flags_no_bursts_f, ".0")
  }
  flags_no_bursts <- get_data(flags_no_bursts_f, sim_type, OUTPUT_LOC)
  flags_no_bursts_p <- flags_no_bursts_plot_all_runs(flags_no_bursts)
  save_plot(paste0(SAVE_LOC,"flags_no_bursts_",init_num_seq,"_",num_sensitive,
                   "_",p_sensitive,".png"),
            flags_no_bursts_p, base_height = 10, base_width = 13)
}
}

# comparing without flags
sim_type           <- WITHOUT_FLAGS
timestep           <- 1
final_time         <- 500
num_runs           <- 20
#                       burst mean    burst p     max total cop
burst_mean_v       <- c(1, 1)     #c(0.25,  1,  4,    1,   1,   1)
burst_p_v          <- c(1.0, 0.1) #c( 0.1,0.1,0.1, 0.001, 1.0, 0.1)
max_total_copies_v <- c(100, 10)  #c( 100,100,100,  100, 100,  10)

source('analyse_output/header.R')
for (i in 1:length(max_total_copies_v))
{
  burst_mean <- burst_mean_v[i]
  burst_p    <- burst_p_v[i]
  max_total_copies <- max_total_copies_v[i]

  without_f <- create_filename(sim_type = sim_type,
                        timestep = timestep, final_time = final_time, num_runs = num_runs,
                        burst_mean = burst_mean, burst_p = burst_p,
                        max_total_copies = max_total_copies
                        )
  gene_f <- create_filename(
                        sim_type = NO_BURSTS,
                        timestep = timestep,
                        final_time = final_time,
                        num_runs = num_runs,
                        init_num_seq = max_total_copies,
                        )
  without <- get_data(without_f, sim_type, OUTPUT_LOC)
  gene <- get_data(gene_f, NO_BURSTS, OUTPUT_LOC)
  without_p_init <- without_plot_init(without, gene)
  without_p_pair <- without_plot_pair(without, gene)
  save_plot(paste0(SAVE_LOC,"without_init_",
                   burst_mean, "_", burst_p, "_", max_total_copies, ".png"),
            without_p_init, base_height = 10, base_width = 12)
  save_plot(paste0(SAVE_LOC,"without_pair_",
                   burst_mean, "_", burst_p, "_", max_total_copies, ".png"),
            without_p_pair, base_height = 10, base_width = 12)
}

# comparing with flags
sim_type         <- WITH_FLAGS
timestep         <- 1
final_time       <- 500
num_runs         <- 20
num_sensitive_v  <- c(10 ,  10)
p_sensitive_v    <- c(0.5, 0.5)
burst_mean_v     <- c(1.0,0.25)
burst_p_v        <- c(0.1, 0.1)
max_time_v       <- c(500,150)
init_num_seq     <- 100

for (i in 1:length(p_sensitive_v))
{
  num_sensitive <- num_sensitive_v[i]
  p_sensitive <- p_sensitive_v[i]
  burst_mean <- burst_mean_v[i]
  burst_p <- burst_p_v[i]
  max_time <- max_time_v[i]
  source('analyse_output/header.R')
  with_f <- create_filename( sim_type = sim_type,
                        timestep = timestep, final_time = final_time, num_runs = num_runs,
                        num_sensitive = num_sensitive, p_sensitive = p_sensitive,
                        burst_mean = burst_mean, burst_p = burst_p
                        )
  gene_f <- create_filename( sim_type = NO_BURSTS,
                        timestep = timestep, final_time = final_time, num_runs = num_runs,
                        init_num_seq = init_num_seq
                        )
  with <- get_data(with_f, sim_type, OUTPUT_LOC)
  gene <- get_data(gene_f, NO_BURSTS, OUTPUT_LOC)
  with_p_init <- with_plot_init(with, gene, max_time)
  with_p_pair <- with_plot_pair(with, gene, max_time)
  save_plot(paste0(SAVE_LOC,"with_init_", num_sensitive, "_", p_sensitive, "_",
                   burst_mean, "_", burst_p, ".png"),
            with_p_init, base_height = 10, base_width = 12)
  save_plot(paste0(SAVE_LOC,"with_pair_", num_sensitive, "_", p_sensitive, "_",
                   burst_mean, "_", burst_p, ".png"),
            with_p_pair, base_height = 10, base_width = 12)
}

#################################################
################### END #########################
#################################################