EricArcher
commented
4 years ago If I've understood the simulation you wanted to specify correctly, here is a script that uses the new fastsimcoal function format that does what you want:
rm(list = ls())
library(strataG)
popSize <- 10000
# mutation rates in loci
mutRates <- c(0.0001, 0.0002, 0.0005, 0.001)
# max no. of alleles in loci
maxRange <- c(3, 6, 9, 12, 15, 20)
# bottleneck sizes
botSize <- c(500, 50, 20)
#historical events
botTime <- 20 # since t0
botLength <- 20 # since the beginning of the bottleneck
splitTime <- 50 # since t0
numPops <- length(botSize)
# Setting the parameters for fastsimcoal
# set demes
demes <- do.call(
fscSettingsDemes,
lapply(1:numPops, function(i) fscDeme(popSize, popSize))
)
# Parameters of microsatellite loci
msat.params <- expand.grid(mu = mutRates, range = maxRange)
msats <- lapply(1:nrow(msat.params), function(i) {
fscBlock_microsat(
num.loci = 1,
mut.rate = msat.params$mu[i],
range.constraint = msat.params$range[i],
chromosome = i
)
})
genetics <- do.call(fscSettingsGenetics, msats)
# Parameters of coalescent analysis
source.sink <- combn(numPops, 2) - 1
after.bot <- lapply(1:ncol(source.sink), function(i) {
source.pop <- source.sink[2, i] + 1
fscEvent(
event.time = botTime,
source = source.sink[1, i],
sink = source.sink[2, i],
new.size = botSize[source.pop] / popSize
)
})
during.bot <- lapply(1:ncol(source.sink), function(i) {
source.pop <- source.sink[2, i] + 1
fscEvent(
event.time = botTime + botLength,
source = source.sink[1, i],
sink = source.sink[2, i],
new.size = popSize / botSize[source.pop]
)
})
before.bot <- lapply(1:ncol(source.sink), function(i) {
fscEvent(
event.time = botTime + botLength,
source = source.sink[1, i],
sink = source.sink[2, i]
)
})
events <- do.call(fscSettingsEvents, c(after.bot, during.bot, before.bot))
p <- fscWrite(
demes = demes,
genetics = genetics,
events = events
)
p <- fscRun(p, num.cores = 5)
konopinski
Hi, I guess this is a very hard time for you after releasing such a revolutionary version of the soft. I really appreciate your work. But I had to downgrade back to 2.0.2 because of the new method of introducing data to fastsimcoal. I managed to rewrite the script to produce input files automatically (24 microsatellite loci, 6 different allele nos. limits, 4 different mutation rates, and 4 populations with different histories) make some automated simulations, but did not manage to make fastsimcoal running. I know there's some problem possibly with the system, because I was not able to run it manually, but the par file from 2.0.2 works excellent. Below you will find a piece of script I used for simulations so that you can understand what automation I'm talking about. Maybe that will inspire you ;) I managed to rewrite most of the blocks, except for the new format event block. The main obstacle is that you cannot provide multiple data as I did before. In the end I found I can build normal data.frames and process them by the new functions so they get the desired classes and attributes, but building them is much more tricky. And after all that work I was very disappointed to see that fsc hangs. I will try to attach param files from the old and new versions, but I'm waiting for my simulations to finish, and I have to recreate all the loops I wrote for building new param file. Cheers Maciek Konopiński