Moran_simulation_clx
clx
2025-01-31
#moran_simulation function
#input: population table of individuals and fitness
moran_simulation <- function(population){
#settings and parameters
#population size
N <- nrow(population)
#initial time
time <- 0
sim_pop<-population$id
events <- data.frame(time = numeric(),
reproducer = numeric(),
replaced = numeric())
cat("Population:", sim_pop, "Time:", time, "\n")
while (length(unique(sim_pop)) > 1) {
#sum of the fitness
total_rate <- sum(population$fitness[sim_pop])
#probability to be choose
rates <- population$fitness[sim_pop] / total_rate
#waiting time
wait_time <- rexp(1, rate = total_rate)
time <- time + wait_time
#choose one individual to reproduce
reproducer_index <- sample(1:N, 1, prob = rates)
#randomly choose an individual that be replaced
to_replace_index <- sample(1:N, 1)
sim_pop[to_replace_index]<-sim_pop[reproducer_index]
cat("Population:", sim_pop, "Time:", time, "\n")
# record replacement relationship and time point
events <- rbind( events, data.frame(time = wait_time, reproducer = reproducer_index, replaced = to_replace_index) )
}
return(events)
}generate tree format from the history
build_coalescent_tree <- function(events, population) {
N <- nrow(population)
nodes <- as.character(population$id)
branch_lengths <- rep(0, N)
ancestors <- nodes
# go through the events backwards
for (i in nrow(events):1) {
event <- events[i, ]
time <- event$time
#branch increasing
branch_lengths <- branch_lengths + time
#renew ancestors of individuals by the events
reproducer<-event$reproducer
ancestors[which(ancestors == event$replaced)] <- reproducer
#emerge nodes that have same ancestor
if (sum(duplicated(ancestors)) > 0) {
#if two individuals have same ancestor
#find the individuals that have same ancestor
rm_index <- which(duplicated(ancestors) | duplicated(ancestors, fromLast = TRUE))
a <- rm_index[1]
b <- rm_index[2]
#generate a new node that is the combination of these two
new_node <- paste0("(", nodes[a], ":", branch_lengths[a], ",", nodes[b], ":", branch_lengths[b], ")")
#remove coalescent nodes
nodes <- nodes[-rm_index]
branch_lengths <- branch_lengths[-rm_index]
ancestors <- ancestors[-rm_index]
#add new node
nodes <- c(nodes, new_node)
branch_lengths <- c(branch_lengths, 0)
ancestors <- c(ancestors, reproducer)
}
}
return(paste0(nodes, ";"))
}n<-6
id<-c(1:n)
fitness<- c(1.1, 1.2, 1.3, 1, 1, 1)
population <- data.frame(id = id, fitness = fitness)
#simulate the population history
events <- moran_simulation(population)## Population: 1 2 3 4 5 6 Time: 0
## Population: 1 4 3 4 5 6 Time: 0.0744159
## Population: 1 4 3 6 5 6 Time: 0.2090928
## Population: 1 1 3 6 5 6 Time: 0.3266537
## Population: 1 1 6 6 5 6 Time: 0.4289079
## Population: 1 1 6 6 1 6 Time: 0.4721206
## Population: 1 1 6 6 1 6 Time: 0.6445721
## Population: 1 1 6 6 1 6 Time: 0.8298878
## Population: 1 1 6 6 1 6 Time: 0.9320347
## Population: 1 1 1 6 1 6 Time: 1.108167
## Population: 1 1 1 6 1 6 Time: 1.151871
## Population: 6 1 1 6 1 6 Time: 1.294143
## Population: 6 1 1 6 1 1 Time: 1.411419
## Population: 6 1 1 6 1 1 Time: 1.568244
## Population: 6 1 1 6 1 1 Time: 1.915598
## Population: 6 6 1 6 1 1 Time: 2.304489
## Population: 1 6 1 6 1 1 Time: 2.615534
## Population: 1 6 1 1 1 1 Time: 2.754292
## Population: 1 6 6 1 1 1 Time: 2.84451
## Population: 1 6 6 1 1 1 Time: 2.967296
## Population: 1 6 6 1 1 1 Time: 3.075922
## Population: 1 6 6 1 1 1 Time: 3.173771
## Population: 1 6 6 1 1 6 Time: 3.450212
## Population: 6 6 6 1 1 6 Time: 3.646756
## Population: 6 6 6 1 1 6 Time: 3.664243
## Population: 6 6 6 1 1 1 Time: 3.838538
## Population: 6 6 6 1 1 6 Time: 4.172361
## Population: 6 6 6 6 1 6 Time: 4.531478
## Population: 6 6 6 6 1 6 Time: 4.837657
## Population: 6 6 6 6 6 6 Time: 4.879973 #events <- events[events$reproducer != events$replaced, ]
events## time reproducer replaced
## 1 0.07441590 4 2
## 2 0.13467694 6 4
## 3 0.11756082 1 2
## 4 0.10225428 4 3
## 5 0.04321265 1 5
## 6 0.17245154 5 5
## 7 0.18531564 2 1
## 8 0.10214696 4 6
## 9 0.17613235 1 3
## 10 0.04370376 1 2
## 11 0.14227244 4 1
## 12 0.11727548 2 6
## 13 0.15682500 5 2
## 14 0.34735415 6 6
## 15 0.38889124 1 2
## 16 0.31104505 6 1
## 17 0.13875821 5 4
## 18 0.09021800 2 3
## 19 0.12278555 5 5
## 20 0.10862643 3 3
## 21 0.09784854 6 5
## 22 0.27644136 2 6
## 23 0.19654392 3 1
## 24 0.01748657 2 3
## 25 0.17429538 4 6
## 26 0.33382307 2 6
## 27 0.35911657 3 4
## 28 0.30617880 2 1
## 29 0.04231614 2 5#draw tree from the history events
library(ape)## Warning: 程辑包'ape'是用R版本4.2.3 来建造的 newick_tree <- build_coalescent_tree(events, population)
#cat( newick_tree)
plot(read.tree(text = newick_tree), type = "phylogram", show.tip.label = TRUE, edge.width = 2)
axisPhylo()
main function of simulation
library(ape)
main_moran_simu<-function(population){
##main:
#simulate the population history
events <- moran_simulation(population)
events <- events[events$reproducer != events$replaced, ]
#events
#draw tree from the history events
newick_tree <- build_coalescent_tree(events, population)
#cat( newick_tree)
plot(read.tree(text = newick_tree), type = "phylogram", show.tip.label = TRUE, edge.width = 2)
axisPhylo()
}an example
n<-6
id<-c(1:n)
fitness<- c(1.1, 1.2, 1.3, 1, 1, 1)
population <- data.frame(id = id, fitness = fitness)
main_moran_simu(population)## Population: 1 2 3 4 5 6 Time: 0
## Population: 1 3 3 4 5 6 Time: 0.4667031
## Population: 1 3 3 4 5 3 Time: 0.6115389
## Population: 1 3 3 3 5 3 Time: 0.748682
## Population: 1 3 3 3 1 3 Time: 0.8133516
## Population: 1 1 3 3 1 3 Time: 0.8167021
## Population: 1 3 3 3 1 3 Time: 0.9607275
## Population: 1 3 3 3 1 3 Time: 1.208978
## Population: 1 3 3 3 1 3 Time: 1.536519
## Population: 1 3 3 3 1 3 Time: 1.567063
## Population: 1 3 3 1 1 3 Time: 1.613481
## Population: 1 3 3 1 3 3 Time: 1.655239
## Population: 1 1 3 1 3 3 Time: 1.72249
## Population: 1 1 3 1 3 3 Time: 2.016673
## Population: 1 1 3 1 3 3 Time: 2.425347
## Population: 1 1 3 1 1 3 Time: 2.546316
## Population: 1 1 3 1 1 3 Time: 2.70798
## Population: 1 1 1 1 1 3 Time: 3.364686
## Population: 1 1 1 1 1 3 Time: 3.508249
## Population: 1 1 1 1 1 1 Time: 3.856845
W-F model:
random_coalescent <- function(N) {
# poopulation initialization
# individuals and the corresponding time till coalescent
population <- as.character(paste("ind", c(1:N)))
pop_t <- rep(0, N)
pop_frame <- data.frame(population, pop_t, stringsAsFactors = FALSE)
# coalescent times from N to 2 individuals
t <- sapply(c(2:N), function(x) rexp(1, 1/choose(x,2)))
# coalescent process
for (i in 1:(N-1)) {
#renew waiting time for tips
pop_frame$pop_t <- pop_frame$pop_t + t[i]
#randomly choose 2 tips
co_pair <- sample(1:length(pop_frame[,1]), 2, replace = FALSE)
co_ind <- pop_frame[co_pair, ]
pop_frame <- pop_frame[-co_pair, ]
# construct new tip by the choosen 2 tips
new_branch <- paste(
"(", co_ind[1, 1], ":", as.numeric(co_ind[1, 2]), ",",
co_ind[2, 1], ":", as.numeric(co_ind[2, 2]), ")"
)
# renew the tips dataframe
pop_frame <- rbind(pop_frame, data.frame(population = new_branch, pop_t = 0, stringsAsFactors = FALSE))
}
# return final tree in newick format
return(paste0(pop_frame$population[1], ";"))
}library(ape)
# an example for 10 individuals
plot.phylo(read.tree(text = random_coalescent(6)))
axisPhylo()