R/get_k80_matrix.R

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

#' Chooses a transition matrix according the Kimura-80 model
#'
#' Returns a table T where Probability(X becoming Y) = T[X, Y]
#' @param time_jump The number of years we jump by for our transition matrix
#' To be valid, we need it to be less than 9*10^7
#' @param K80_ALPHA probability of transition per year
#' @param K80_BETA  probability of transversion per year
#' @return A table of numerical values that are the transition matrix for the K-80 model
#' @export

get_k80_matrix <- function(time_jump,
                           K80_ALPHA = 10^(-8),
                           K80_BETA  = 10^(-9))
{
  if (time_jump > 9*(10^7)) {
    stop("Time jump too large for the Kimura 80 matrix")
  }

  k80_alpha_effective = K80_ALPHA * time_jump
  k80_beta_effective  = K80_BETA  * time_jump
  mat <- matrix(1:16, nrow = 4, ncol = 4)
  dimnames(mat) <- list(NUCLEOTIDES, NUCLEOTIDES)
  mat <- as.table(mat)

  temp <- c(rep(k80_alpha_effective, 2), rep(k80_beta_effective, 2))
  mat["T", ] <- temp
  mat["C", ] <- temp
  mat["A", ] <- rev(temp)
  mat["G", ] <- rev(temp)
  diag(mat) <- 1-(k80_alpha_effective + 2*k80_beta_effective)

  return(mat)
}