R/kill_copies.R

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

#' Selects a set of sequences randomly given the max copy number
#' This simulates the 'death' of certain transposons
#' @param copy_numbers The copy numbers for each sequence
#' @param max_copy_number The number of sequences for the next generation
#' @return A new set of copy numbers after some sequences have died
#' @export
kill_copies <- function(copy_numbers, max_copy_number)
{
  # the sequences that we have
  seq_indices <- 1:length(copy_numbers)

  # calculating how many we have to remove
  total_copy_number <- sum(copy_numbers)
  effective_max <- min(rpois(n=1, lambda=max_copy_number),
                       total_copy_number)
  number_to_kill <- total_copy_number - effective_max

  if(number_to_kill > 0)
  {
    for(i in 1:number_to_kill)
    {
      killed_index <- sample(seq_indices, size=1,
                             prob = normalise_distribution(copy_numbers))
      copy_numbers[killed_index] <- copy_numbers[killed_index] - 1
    }
  }
  return(copy_numbers)
}