R/recombine_sequences.R

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

#' Takes two sequences and returns a recombination of the two sequences
#'
#' @param sequences The two sequences that are going to recombine
#' @param mean_num_events The mean number of recombination events that occur,
#' from a Poisson distribution
#' @return A new sequence, generated by recombination between two sequences
#' @export
#'
recombine_sequences <- function(sequences, mean_num_events)
{
  cur_index <- sample.int(2, size=1) # which sequence we are currently reading from
  new_sequence <- sequences[[cur_index]]

  # how many recombination events are happening
  num_recombination <- rpois(n=1, lambda=mean_num_events)

  if (num_recombination > 0) {
    n <- length(sequences[[cur_index]])
    other_index <- 3 - cur_index # the other sequence to borrow from

    posn_of_recomb <- sort(sample(2:n, size=num_recombination, replace = FALSE))
    # if we have an odd number of recombinations, make sure you recombine until the end
    if (num_recombination %% 2 != 0) {
      posn_of_recomb <- c(posn_of_recomb, n+1)
    }

    # actually recombine the sequences
    i <- 1
    while(i < length(posn_of_recomb)) {
      new_sequence[posn_of_recomb[i]:(posn_of_recomb[i+1]-1)] <-
        sequences[[other_index]][posn_of_recomb[i]:(posn_of_recomb[i+1]-1)]
      i <- i+2
    }
  }
  return(new_sequence)
}