Simulation of a DNA site for likelihood computation

In this notebook, we’ll see how to simulate DNA evolution along a branch of arbitrary length, under the Jukes and Cantor (1969) model. Then we will compute the likelihood of a phylogenetic tree using such simulations.

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.926 0.024 0.029 0.021 

Then we do the simulations. Branch length 1:

simulateSitesAndCollectFrequencies("A", 1, 1000)

simulatedStates
    A     C     G     T 
0.549 0.145 0.148 0.158 

Then we do the simulations. Branch length 10:

simulateSitesAndCollectFrequencies("A", 10, 1000)

simulatedStates
    A     C     G     T 
0.254 0.217 0.261 0.268 

Do you find those proportions expected or surprising?

Computing the likelihood of a tree using simulations

# Compute the probability of a tree
computeTreeProbabilityBySimulating <- function(bl_a1, bl_c, bl_a2, bl_y) {
    Nsuccess = 0
    Nreps = 10000
    iter = 0
    while (iter < Nreps) {
        x <- drawNewState()
        y <- simulateAlongBranch(x, bl_y)
        A1 <- simulateAlongBranch(y, bl_a1)
        C <- simulateAlongBranch(y, bl_c)
        A2 <- simulateAlongBranch(x, bl_a2)
        if (A1=="A" && C=="C" && A2=="A") {
            Nsuccess <- Nsuccess + 1
        }
        iter <- iter + 1
    }
    return (Nsuccess/Nreps)
}
p = computeTreeProbabilityBySimulating(0.1, 0.1, 0.4, 0.3)
print(paste ("Likelihood of ((A:0.1, C:0.1):0.3, A:0.4); for pattern {A,C,A}: ", p, sep=""))

[1] "Likelihood of ((A:0.1, C:0.1):0.3, A:0.4); for pattern {A,C,A}: 0.0039"

print(paste("Log-likelihood of ((A:0.1, C:0.1):0.3, A:0.4); for pattern {A,C,A}: ", log(p), sep=""))

[1] "Log-likelihood of ((A:0.1, C:0.1):0.3, A:0.4); for pattern {A,C,A}: -5.54677872584654"

Comparison to likelihood computed with usual software

In common software, substitutions to the same state are not considered, and therefore branch lengths need to be rescaled. Indeed, branch lengths correspond to expected numbers of substitutions: in our case, we expect 4/3 times more substitutions than in common software. Therefore, to compare the likelihoods, we have to run common software on branch lengths that are 3/4 those of the tree, i.e.: ((A1:0.075, C:0.075):0.225, A2:0.3); When I do that, I get : loglk=-5.480185. So the Gillespie sampling approach is not too far from the Felsenstein algorithm exact approach.