R/copy_number_growth.R

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

#' To model copy number variation over time
#' @export
copy_number_growth <- function(init_num_seq, # the initial number of sequences (all active)
                               num_jumps,    # the number of generations to simulate for
                               burst_threshold, # the probability of a burst event
                               burst_mean, # how much to increase the copy number by
                               inactivity_probability, # the probability of inactivity
                               max_copy_number, # what is the maximum copy number to reach
                               max_active_copy_number, # what is the maximum number of active copies allowed
                               to_display = TRUE) # dummy parameter for display
{

  num_seq <- init_num_seq
  num_active <- init_num_seq

  total_pop  <- rep(NA_integer_, num_jumps)
  active_pop <- rep(NA_integer_, num_jumps)
  total_pop[1] <- num_seq
  active_pop[1] <- num_active

  for(t in 2:num_jumps)
  {
    # add mutations that decrease activity
    num_active <- num_active - rbinom(1, num_active, inactivity_probability)

    new_seq <- 0
    # allow burst events that increase the number of active elements
    for (j in 1:num_active)
    {
      if (runif(1) < burst_threshold)
      {
        new_seq <- new_seq + rpois(1, burst_mean)
      }
    }
    # add these back into the mix
    num_seq <- num_seq + new_seq
    num_active <- num_active + new_seq

    total_pop[t] <- num_seq
    active_pop[t] <- num_active

    if (num_seq > max_copy_number || num_active > max_active_copy_number ||
        num_active < 1)
    {
      break
    }
  }

  if (to_display)
  {
    plot(total_pop, type="l", col="red",
         xlab = "Population", ylab = "Time")
    lines(active_pop, col="blue")
    legend(x=1, y = max(total_pop, na.rm = TRUE),
           legend = c("Total", "Active"), col = c("red", "blue"),
           lty = c(1, 1), cex=1)
  }

  final_copy_number <- num_seq
  percentage_active <- (100*num_active)/num_seq

  if (any(total_pop > max_copy_number, na.rm=TRUE))
  {
    result <- "total pop exceeded"
  }
  else if (any(active_pop > max_active_copy_number, na.rm=TRUE))
  {
    result <- "active pop exceeded"
  }
  else if(0 %in% active_pop)
  {
    result <- "loss of activity"
  }
  else
  {
    obs_cor <- cor(total_pop,
                   seq(from=init_num_seq, to=max(total_pop, na.rm = TRUE),
                       length.out = length(total_pop)),
                   use = "complete.obs")
    if(is.na(obs_cor))
    {
      result <- "no growth"
    }
    else if(obs_cor > 0.9)
    {
      result <- "linear growth"
    }
    else
    {
      result <- "non-linear growth"
    }
  }
  return(list(result=result,
              percentage_active=percentage_active,
              final_copy_number=final_copy_number))
}