Simulated FST value ~1/4 of the expected theoretical value
Boki Lai
Hello everyone,
I'm currently using SLiM 5 to simulate population divergence under various migration rates. For each parameter set, I ran 100 replicate simulations and calculated the average FST values using calcFST(). When I'm comparing the resulting FST values to the theoretical expectations based on Wright’s island model. However, I’ve consistently found that the FST values produced by SLiM are approximately one-fourth of the expected values.
For example, with two populations of 10,000 diploid individuals and a symmetric migration rate of 0.001, the expected FST (using the formula FST = 1 / (4Nm + 1)) is ~0.024; however, SLiM consistently reports a value of around 0.006.
I would like to know if this is expected behavior in SLiM, or if I’m missing something important in the manual?
Bless,
Ariel
Ben Haller
Step one would be for you to supply your complete SLiM script here, in a form that can be run directly in SLiMgui (i.e., not requiring any command-line arguments to be supplied). Without that, it is hard to guess what you might be doing wrong. See https://doi.org/10.1371/journal.pcbi.1010540.
Beyond that, there are several steps you can take to investigate this yourself:
- Output your simulation state into, say, VCF or MS, and use other existing tools to measure the FST in your model, to compare to the value you're getting from SLiM.
- Code up the same simulation in a different evolutionary simulation framework, to see whether you get the same or a different result.
- Look at the code for calcFST(), available in the Eidos terminal with functionSource("calcFST"), to see whether it is doing what you expect or not, as there are different ways of measuring FST; the documentation might also be informative on that count.
If you pursue any of those steps, you might report your findings back here. But as I said, step one is to show us what you're actually doing by supplying your script.
Cheers,
-B.
Benjamin C. Haller
Messer Lab
Cornell University
Boki Lai wrote on 8/4/25 5:16 AM:
--
SLiM forward genetic simulation: http://messerlab.org/slim/
---
You received this message because you are subscribed to the Google Groups "slim-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to slim-discuss...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/8a4bda88-3394-4631-a9b9-ac00a7447436n%40googlegroups.com.
Boki Lai
Below is my SLiM script. The first 20,000 generations are intended to allow FST values to stabilize under different migration rates, and after generation 20,000, I simulate alternating bottlenecks. The simulation tests three pairs of populations, each with a different migration rate (5e-4, 1e-4, and 1e-3).
Bless,
Ben Haller
Took a quick look at your model. My initial impression is that the runtime of the model doesn't look nearly long enough for it to equilibrate; I would guess that it is still building up genetic diversity, and that the FST values are still rising. Do you have reason to believe otherwise?
As I wrote before, it would be helpful if you were to measure the FST using a different method, such as a VCF or MS dump from the simulation, to confirm that the calcFST() method is producing reasonable results. I had a couple of other suggestions too. Good luck!
Cheers,
-B.
Benjamin C. Haller
Messer Lab
Cornell University
To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/6605c805-6e1a-4f63-a9a2-32e49e29d53cn%40googlegroups.com.
Jose Luis Mijangos
Boki Lai
Blessings,
Michael Kopp
initialize() {
initializeMutationRate(1e-7);
initializeMutationType("m1", 0.5, "f", 0.0);
initializeGenomicElementType("g1", m1, 1.0);
initializeGenomicElement(g1, 0, 999999);
initializeRecombinationRate(1e-8);
}
1 early() {
sim.addSubpop("p1", 1000);
sim.addSubpop("p2", 1000);
p1.setMigrationRates(p2, 0.001);
p2.setMigrationRates(p1, 0.001);
log = community.createLogFile("sim_log.txt", logInterval=10);
log.addCycle();
log.addCustomColumn("FST", "calcFST(p1.haplosomes, p2.haplosomes);");
}
20000 late() { }
With N = 1000 and m = 0.001, the expected Fst = 1/(1 + 4Nm) = 0.2, but the value recorded in the log-file fluctuates around a value of approximately 0.5. It would be nice to have an explanation for this discrepancy. Is my interpretation of parameters wrong?
Best wishes,
Michael.
Gregor Gorjanc
Gregor
--
SLiM forward genetic simulation: http://messerlab.org/slim/
---
You received this message because you are subscribed to the Google Groups "slim-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to slim-discuss...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/f71ecbb8-c3de-48a4-915e-e2b58a3e6231n%40googlegroups.com.
Miguel de Navascués
Note that you are using the expected Fst for the INFINITE island model.
In the finite island model, expected Fst is lower and depends on the
number of demes, if I remember right.
Best,
Miguel
> However, I’ve consistently found that the FST values produced by
> SLiM are approximately one-fourth of the expected values.
>
> For example, with two populations of 10,000 diploid individuals and
> a symmetric migration rate of 0.001, the expected FST (using the
> formula FST = 1 / (4Nm + 1)) is ~0.024; however, SLiM consistently
> reports a value of around 0.006.
>
> I would like to know if this is expected behavior in SLiM, or if I’m
> missing something important in the manual?
>
>
> Bless,
> Ariel
>
> SLiM forward genetic simulation: http://messerlab.org/slim/ <http://
> messerlab.org/slim/>
> You received this message because you are subscribed to the Google
> Groups "slim-discuss" group.
> To unsubscribe from this group and stop receiving emails from it, send
> discuss+u...@googlegroups.com>.
> discuss/f71ecbb8-c3de-48a4-915e-e2b58a3e6231n%40googlegroups.com
> e2b58a3e6231n%40googlegroups.com?utm_medium=email&utm_source=footer>.
--
Miguel de Navascués
Researcher at INRAE
Board Member of PCI Evolutionary Biology & PC Journal
Adress:
UMR CBGP, INRAE
Centre de Biologie pour la Gestion des Populations
755 avenue du campus Agropolis, CS30016
34988 Montferrier-sur-Lez cedex (France)
miguel.navascues AT inrae.fr
@Miguel_de...@mastodon.online
https://mnavascues.github.io/Website/
https://umr-cbgp.fr/membre/de-navascues/
Ben Haller
Fst = 1 / (1 + 4Nem(d/(d-1))^2)
That comes out to 0.05882353 for the parameters used here. I added a bit of plotting and calculation code to the model (shown after my signature), and SLiM produces a mean of 0.0588 for the run I did – right on the money. (Of course there's a fair bit of stochasticity, so this run being so close to the expected value is probably a bit lucky!) A plot generated by SLiMgui with this code is attached.
Cheers,
-B.
initialize()
{
initializeMutationRate(1e-7);
initializeMutationType("m1", 0.5, "f", 0.0);
initializeGenomicElementType("g1", m1, 1.0);
initializeGenomicElement(g1, 0, 999999);
initializeRecombinationRate(1e-8);
}
1 late()
{
sim.addSubpop("p1", 1000);
sim.addSubpop("p2", 1000);
p1.setMigrationRates(p2, 0.001);
p2.setMigrationRates(p1, 0.001);
defineConstant("LOG", community.createLogFile("sim_log.txt"));
LOG.addCycle();
LOG.addCustomColumn("FST", "calcFST(p1.haplosomes, p2.haplosomes);");
}
late()
{
if (sim.cycle % 10 == 0)
{
LOG.logRow();
cycles = slimgui.logFileData(LOG, "cycle");
fst = slimgui.logFileData(LOG, "FST");
plot = slimgui.createPlot("FST",
xrange=c(1, community.estimatedLastTick()), yrange=c(0, 1),
xlab="cycle", ylab="FST", width=1000, height=300);
plot.lines(cycles, fst);
if (sim.cycle > 10000)
{
equilibrated_cycles = (cycles >= 10000);
mean_fst = mean(fst[equilibrated_cycles]);
plot.abline(h=mean_fst, color="black");
plot.text(50000, mean_fst + 0.03, labels=format("%.4f", mean_fst));
}
}
}
50000 late()
{
Michael Kopp
initialize() {
defineConstant("d", 3); // number of demes
defineConstant("N", 1000); // local population size
defineConstant("m", 0.001); // migration rate
defineConstant("TMAX", 20000); // number of generations
initializeMutationRate(1e-7);
initializeMutationType("m1", 0.5, "f", 0.0);
initializeGenomicElementType("g1", m1, 1.0);
initializeGenomicElement(g1, 0, 99999);
initializeRecombinationRate(1e-8);
}
1 early() {
for(i in 1:d){
sim.addSubpop(paste0("p" + i), N);
}
for(pi in sim.subpopulations)
for(pj in sim.subpopulations) {
if (pi != pj) pi.setMigrationRates(pj, m/(d-1));
}
log = community.createLogFile("fst_log.txt", logInterval=10);
log.addCycle();
log.addCustomColumn("FST", "calcFST(p1.haplosomes, p2.haplosomes);");
catn("expected Fst = " + 1/(1 + 4 * N * m * (d/(d-1))^2));
}
TMAX late() {
fst = readCSV("fst_log.txt", colTypes = "ff");
fst = fst.asMatrix();
plot = slimgui.createPlot("Fst between demes 1 and 2",
xrange=c(0,TMAX),
xlab="Generation", ylab="Fst");
plot.lines(fst[,0], fst[,1], "black"); // fst over time
x = length(fst[,1]);
catn("Mean Fst = " + mean(fst[asInteger(x/2):(x-1),1])); // second half of simulation
}
Best wishes,
Michael.
Ben Haller
The equation being used here is just an approximation; it is not exact. I imagine you're bumping against the limitations of that approximation. As you add more demes the equilibration time might also increase, so 20000 might not be sufficient (which is one reason why I increased the max time to 100000 in my version of your code). I think there's probably literature on all this, and an exact formula might be available somewhere, but that's beyond my knowledge; if you really want to delve that deeply I'd suggest doing a literature search. It might also be helpful to compare against results from a different simulator. Good luck!
Cheers,
-B.
Benjamin C. Haller
Messer Lab
Cornell University
To unsubscribe from this group and stop receiving emails from it, send an email to slim-discuss...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/01ded0d1-343c-44c2-b0c5-712daab256f3n%40googlegroups.com.