Counting fixed mutations in subpopulations?

[David Rinker]

Hi,

I have been running some simulations (SLiMulations?) where a founding population gets split into two subpopulations which then evolve indepepdently.

I'm interested in seeing which mutations fix in each subpopulation. The option outputFixedMutations() only summarizes which mutations are no longer segregating in any populaiton, so it is not exactly what I want. 

Wondering if there is an easy way to report fixed mutations on a per-population basis?

Thanks fo any help

-David

[Ben Haller]

  Hi!  Well, as you say, mutations that have fixed in only one subpopulation but not in the other are not considered “fixed” in SLiM, so they won’t appear in the output generated by outputFixedMutations().  To be “fixed” in SLiM, a mutation has to be fixed across the entire population.  (Unless, that is, you intervene in SLiM’s normal mutation processing by calling removeMutations() yourself; if you do that, you’re on your own :->).

  If you want to find/count mutations that have fixed in one subpopulation but not in the other, the mutationFrequencies() and/or mutationCounts() methods of SLiMSim should be useful.  You can pass in the subpopulation you’re interested in, and get frequencies/counts for all mutations for that subpopulation only.  A frequency of 1.0 would then mean “fixed in the specified subpopulation”.  And then you can dump whatever information you’re interested in, regarding those mutations, to the output yourself, or perhaps use outputMutations() to do so.

  Let me know if you have further questions!  :->

Cheers,

-B.

Benjamin C. Haller

Messer Lab

Cornell University

[David Rinker]

Thank you Ben.  Very helpful and exactly what I needed to know

[David Rinker]

Hello again.

I'm still struggling to get the output that I want for my sub-populations--basically what I would like to get out is a summary (per-population) of mutation frequencies along with each mutation's ID...

What I've tried:

1)  sim.mutationFrequencies() outputs a list of frequencies but I can't tell how to associate each of them back to a specific mutation

2) I have also tried simply setting convertToSubstitution=F  figuring I could just parse outputFull()  for what I need.  However, this does not do what I expected it to and it appears that fixed sites are still removed and only included in the separate outputFixedMutations().

Any further hand-holding here would be appreciated! :)

  

On Tuesday, February 27, 2018 at 6:24:17 AM UTC-6, Ben Haller wrote:

[Ben Haller]

  Hi David.  OK, so.  Regarding your first question, on the output of information about mutations, I think what you want to do is:

1. Get the vector of mutations you're interested in (from sim.mutations, or whatever other way you like); let's call that vector "muts"

2. Call sim.mutationFrequencies(pX, muts) for each subpopulation pX that you're interested in, and output the results of each call somewhere

3. Call muts.id to get the mutation identifiers, and output those somewhere

  Now you have the frequencies and the mutation identifiers, and they correspond one-to-one, so you can match them up – this frequency is for this mutation ID.  And that lets you also cross-correlate the information with the output of outputFull() or outputFixedMutations(), if you wish, using the mutation identifiers as a lookup key.

  Regarding your second question, on convertToSubstitution=F, if you're not seeing any effect in your model, that might perhaps be because the mutations you're looking at have not, in fact, fixed?  Remember that mutations have to fix population-wide, not just in a single subpopulation, before they are considered "fixed" by SLiM.  If they never reach population-level fixation, the value of that flag will make no difference.

  So, just as example, if we look at a simplified version of the model we were discussing in the other thread:

initialize() {

   initializeMutationRate(1e-7);

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

   initializeGenomicElementType("g1", m1, 1.0);

   initializeGenomicElement(g1, 0, 99999);

   initializeRecombinationRate(1e-8);

}

1 { sim.addSubpop("p1", 10000); }

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

  This printed out a list of about 100 fixed mutations at the end, in the run I did.  If I add a call to set convertToSubstitution to F for m1, changing it to:

initialize() {

   initializeMutationRate(1e-7);

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

   m1.convertToSubstitution = F;

   initializeGenomicElementType("g1", m1, 1.0);

   initializeGenomicElement(g1, 0, 99999);

   initializeRecombinationRate(1e-8);

}

1 { sim.addSubpop("p1", 10000); }

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

  Now this model prints out no fixed mutations at the end, because SLiM has been told not to detect/substitute fixed m1 mutations.  It also runs more slowly (since all mutations are kept as segregating in the population), and it looks a bit different in SLiM since more and more mutations build up at a frequency of 1 in the chromosome view.  If you wanted to output information on these mutations, you would use outputFull() (or your own custom output code, of course).

  So, if you've got a model in which this flag does not seem to be working properly, please supply the full code for that model and I'll see if I can clarify what's going on.  I hope this helps!

Cheers,

-B.

Benjamin C. Haller

Messer Lab

Cornell University

[David Rinker]

Thank you.  Both solutions are working for me.

Regarding my second question, I simply didn't realize that convertToSubstitution required specifying the mutation type.  I was going by the SLiM Reference Sheet which doesn't mention this, so wrongly assumed that not specifying mutations type would apply to all.  Guess not :)

Ok, looks like I've got what I need right now.

Thanks again!

[David Rinker]

One question on generation output

What's the secret to passing mutationFrequencies() through sapply() (ie, in a line-by-line manner)?  

So far, everything I've tried results in the entire vector being written into each line of my output file.

[Ben Haller]

Ah, I guess you mean in the generation cycle flowchart?  It is shown elsewhere on the reference sheet as being a property of MutationType, but I can see how the flowchart's mention of it could be confusing.  I'll fix that, somehow.

Cheers,

-B.

[Ben Haller]

  The mutationFrequencies() method is passed a vector of mutations, and returns a vector of frequencies, so I'm not sure what you're trying to do with sapply().  Why don't you contact me off-list about this; doesn't sound like it's likely to be of general interest, so let's take it off-list.  :->

Cheers,

-B.