Having problems with a stepping stone model and get printed frequencies

[Héctor Alessandro López Hérnandez]

I want to retrieve the information of an allele "m2" and calculate their frequency in each subpopulation (COUNT x COUNT in a stepping stone model). The code is the next one:

initialize() {

defineConstant("COUNT",50);

initializeSex("A");

initializeMutationRate(1e-7);

initializeMutationType("m1", 0.5, "f", 0.0);

initializeMutationType("m2", 0.5, "f", 0.01);

initializeGenomicElementType("g1", m1, 1.0);

initializeGenomicElement(g1, 0, 99999);

initializeRecombinationRate(1e-7);

}

1 early() {

for (i in 0:(COUNT-1))

sim.addSubpop(i, 500);

subpops = sim.subpopulations;

for (i in 1:(COUNT-1)) subpops[i].setMigrationRates(i-1, 0.2);

}

1 late() {

target=sample(p1.genomes,50);

target.addNewDrawnMutation(m2,rdunif(1,0,sim.chromosome.lastPosition));

}

1:1000 late() {

for (i in 1:(COUNT)) {

sweep = sim.mutationsOfType(m2);

if (size(sweep) > 0) {

for (j in 1:(COUNT)) {

freq = sim.mutationFrequencies(get(paste0("p", j)), sweep);

catn(sim.cycle + ": " + freq + " p" + j);

}

}

}

}

1000 late() { sim.outputFixedMutations();}  

What I am trying in the last part of the code is explore the allele frecuency and obtain the thick+ the frequency of the allele + the population from which was retrieve.

Thanks

AL

[Ben Haller]

Hi Héctor!

I'm not sure in what way the algorithm you want to execute is COUNT x COUNT, rather than just COUNT.  You have two nested `for` loops, for `i` and `j`, but you do not seem to use the loop index variable `i` at all in the code:

1:1000 late() {

for (i in 1:(COUNT)) {

sweep = sim.mutationsOfType(m2);

if (size(sweep) > 0) {

for (j in 1:(COUNT)) {

freq = sim.mutationFrequencies(get(paste0("p", j)), sweep);

catn(sim.cycle + ": " + freq + " p" + j);

}

}

}

}

So I don't understand what the purpose of the outer loop in `i` is.  It looks to me like you could just write:

1:1000 late() {

    sweep = sim.mutationsOfType(m2);
    if (size(sweep) > 0) {

        for (j in seqLen(COUNT)) {
            freq = sim.mutationFrequencies(j, sweep);

            catn(sim.cycle + ": " + freq + " p" + j);
        }
    }
}

Note the use of `seqLen(COUNT)`; this is the best way to generate a sequence from 0 to COUNT-1 (which seems to me to be what you want, not 1 to COUNT).  You could use `seqLen(COUNT)` elsewhere in your code too, such as in your 1 early() event.

Also note that mutationFrequencies() can take a subpopulation identifier (`j` here); you don't need to construct a string using `paste0()` the way you were doing.

Finally, note that you could keep a reference to the sweep mutation with defineConstant(), rather than looking it up each time with sim.mutationsOfType(m2).  This is not very important, but might make your code a little more readable.  See section 9.9 for an example of this technique.

Apologies if I've misunderstood what you're trying to do -- particularly the purpose of the outer `for` loop in `i`.  If I did misunderstand, please clarify the intention.  Happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University


Héctor Alessandro López Hérnandez wrote on 5/14/24 2:35 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 on the web visit https://groups.google.com/d/msgid/slim-discuss/b31b78c7-336b-4751-b6e2-40dc16ec62aen%40googlegroups.com.

[Héctor Alessandro López Hérnandez]

Hi Ben

Thank you very much for your answer.

COUNT actually cames from the recipe #5.3.1 "A linear stepping stone model" which, I guess, is the number of subpopulations. What I am working on is a stepping stone model with a size of NxN where N is the number of subpopulations (being "COUNT" our "N") and adding a beneficial mutation in one subpopulation (p1), experimenting with rate migration to see if this mutation lost or fixed and finally tracking their changes in each subpopulation and in each generation.

I will check your cited example and thank you again! 

Héctor Alessandro López Hernández

Correo alternativo: Comunidad UNAM

"Per aspera, ad astra"

[Ben Haller]

Hi Héctor,

Sorry for the slow reply; life is a bit busy right now.  Finished teaching a SLiM workshop in Copenhagen, then relocated to Italy for the next workshop I will teach.  :->

Yes, what I find confusing about your code and your question is that you say you want to model an NxN subpopulation matrix, but your code only creates N subpopulations in a linear format – because you started from 5.3.1 and have not yet adapted it to be NxN instead of just N.  So then your output algorithm, where you seem to want to loop over N in both i and j, is kind of putting the cart before the horse, since you only have N subpopulations.  Anyway, I'd recommend that you look at sections 5.3.3 and 5.3.4, which show making an NxN matrix of subpopulations, which seems to be what you want.

Good luck and happy modeling!

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

Héctor Alessandro López Hérnandez wrote on 5/15/24 4:30 AM: