Simulation of a DNA site

In this notebook, we’ll see how to simulate DNA evolution along a branch of arbitrary length, under the Jukes and Cantor (1969) model.

Simulating along a branch

# Useful function for drawing a DNA state, for instance at the root
drawNewState <- function (){
    return(sample(x=c("A","C","G","T"), size=1))
}
# Function to simulate along a branch using the Gillespie algorithm under the Jukes Cantor model (1969)
# In this case, we consider that there can be a mutation towards the same state.
# So the sum of the rates is 0.25*4=1.0
simulateAlongBranch <- function(startingState, branchLength) {
    l <- 0.0
    current <- startingState
    while (l < branchLength) {
        l <- l + rexp(rate=1.0, n=1) # 1.0 is the sum of the rates of all possible events
        if (l < branchLength) {
            current <- drawNewState()
        }
    }
    return (current)
}

This way we can simulate along a branch. Let’s do it several times, for different branch lengths. First we implement a function to do that:

simulateSitesAndCollectFrequencies <- function (startingState, branchLength, numSites) {
  simulatedStates = c()
  for (i in 1:numSites) {
    simulatedStates = c(simulatedStates, simulateAlongBranch(startingState, branchLength))
  }  
  return(table(simulatedStates)/numSites)
}

Then we do the simulations. Branch length 0.1:

simulateSitesAndCollectFrequencies("A", 0.1, 1000)

simulatedStates
    A     C     G     T 
0.931 0.028 0.018 0.023 

Then we do the simulations. Branch length 1:

simulateSitesAndCollectFrequencies("A", 1, 1000)

simulatedStates
    A     C     G     T 
0.503 0.191 0.156 0.150 

Then we do the simulations. Branch length 10:

simulateSitesAndCollectFrequencies("A", 10, 1000)

simulatedStates
    A     C     G     T 
0.255 0.239 0.255 0.251 

Do you find those proportions expected or surprising?