R/find_parameter_balance_copy_number_variation.R

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

#' To find a balance between various parameters that results
#' in linear growth of active transposons
#' @export
find_paramater_balance_copy_number_variation <- function(num_samples = 25,
                                                         num_jumps = 500,
                                                         init_num_seq = 20,
                                                         burst_threshold_min = -3,
                                                         burst_threshold_max = -1,
                                                         burst_mean = 1,
                                                         inactivity_rate_min = -4,
                                                         inactivity_rate_max = -1,
                                                         max_copy_number = 10^6,
                                                         max_active_copy_number = 10^3,
                                                         to_display = TRUE)
{
  burst_threshold_v <- rep(10^seq(burst_threshold_min, burst_threshold_max,
                                  length.out=num_samples),
                           times=num_samples)
  inactivity_rate_v <- rep(10^seq(inactivity_rate_min, inactivity_rate_max,
                                  length.out=num_samples),
                           each=num_samples)
  results <- character(length = num_samples*num_samples)
  percentage_actives <- rep(NA_real_, num_samples*num_samples)
  final_copy_numbers <- rep(NA_integer_, num_samples*num_samples)

  for (counter in 1:(num_samples*num_samples))
  {
    burst_threshold <- burst_threshold_v[counter]
    inactivity_rate <- inactivity_rate_v[counter]

    x <- copy_number_growth(init_num_seq, num_jumps, burst_threshold,
                            burst_mean, inactivity_rate,
                            max_copy_number, max_active_copy_number,
                            to_display = FALSE)
    results[counter] <- x$result
    percentage_actives[counter] <- x$percentage_active
    final_copy_numbers[counter] <- x$final_copy_number

  }

  data <- data.frame(burst_threshold = burst_threshold_v,
                     inactivity_rate = inactivity_rate_v,
                     result = results,
                     percentage_active = percentage_actives,
                     final_copy_number = final_copy_numbers)

  if (to_display)
  {
    cols <- scales::brewer_pal(palette="Set1")(9)
    p <- ggplot2::ggplot(data = data,
                         ggplot2::aes(x=inactivity_rate,
                                      y=burst_threshold)) +
         ggplot2::geom_point(ggplot2::aes(shape=result,
                                          colour=percentage_active,
                                          size=final_copy_number)) +
         ggplot2::scale_x_log10() +
         ggplot2::scale_y_log10() +
         ggplot2::scale_colour_gradientn(colours=cols) +
         ggplot2::scale_size_continuous(range=c(2, 4), trans="log10") +
         ggplot2::labs(x="Probability of an active sequence becoming inactive",
                       y="Probability of an active sequence bursting",
                       shape="Result of simulation",
                       colour="% of active TEs",
                       size="Final number of TEs") +
         ggplot2::theme(
           axis.title  = ggplot2::element_text(size = 18, face = "bold"),
           axis.text.y = ggplot2::element_text(size = 18, face = "bold"),
           axis.text.x = ggplot2::element_text(size = 18, face = "bold")
         ) +
         ggplot2::theme(
           legend.title = ggplot2::element_text(size = 18),
           legend.text = ggplot2::element_text(size = 18)
         )


    print(p)
  }
  return(p)
}