Efficient simulation of pleiotropy

[Joshua Schraiber]

Hi,

I'm trying to do some simulations with a braindead model of pleiotropy (perhaps to be made more complicated later) where every mutation has an iid effect on n different traits, and then selection also acts in an independent way across traits. I think that with Gaussian selection, that means that the selection acting on an individual is proportional to exp(\sum_j x_{ij}^2/(2*Vs)) where x_{ij} is the the phenotype of the ith individual at the jth trait, and in an additive model x_{ij} = \sum_l \beta_{jl} *G_{il} where \beta_{jl} is the effect size of locus l on trait j and G_{il} is the genotype of individual i at locus l. So all I need to really do is calculate \sum_j x_{ij}^2 for the fitness.

However, I've tried a couple ways to set this up and they are all very slow. I think the fundamental issue is that I can't use sumOfMutationsOfType, so I have to do a loop myself. 

First, here's how I set up the causal mutations:

mutation(m2) {

effects = rnorm(n,0,sd=sqrt(sigma2));

mut.setValue("effect", effects);

return T;

}

Then here is one way I've tried to compute the G2 = \sum_j x_{ij}^2:

function (float)get_G2(o<Individual> inds) {

//get fixed effects

fixed_effects = float(n);

for (sub in sim.substitutions) {

fixed_effects = fixed_effects + sub.getValue("effect");

}

 res = float(length(inds));

 for (i in 0:(length(inds)-1)) {

  muts = inds[i].genomes.mutationsOfType(m2);

  effect = length(muts) ? muts.getValue("effect") else rep(0,n);

  effect = matrix(effect,nrow=n);

  G = apply(effect,0,"sum(applyValue);");

  G = G + 2*fixed_effects;

  G2 = G*G;

  res[i] = sum(G2);

 }

return res;

}

and the second way is 

function (float)get_G2(o<Individual> inds) {

//get fixed effects

fixed_effects = float(n);

for (sub in sim.substitutions) {

fixed_effects = fixed_effects + sub.getValue("effect");

}

res = float(length(inds));

for (i in 0:(length(inds)-1)) {

muts = inds[i].genomes.mutationsOfType(m2);

cur_G = 0;

for (mut in muts) {

cur_effects = mut.getValue("effect");

cur_G = cur_G + cur_effects;

}

cur_G = cur_G + 2*fixed_effects; //NB THE 2 DUE TO DIPLOIDY

G2 = cur_G*cur_G;

res[i] = sum(G2);

}

return res;

}

Both are pretty excruciating---presumably it's the sum over individuals that's the big problem here. Do you think there's any way to speed this up? 

Thanks!

PS: I'm using SLIM 4 since this is an update to some old scripts that were in SLIM 4, but I can port to SLIM 5 if you think that will help!

[Ben Haller]

Hi Joshua!

Yes, it's pretty painful to do this kind of thing without sumOfMutationsOfType().  I'm not sure that I have any solution that would obviously be more performant, although you could certainly play around with things like putting each of the n effects of a mutation into a separate setValue() property; that would then let you do something like sum(muts.getValue("EFFECT_3")) which might be faster than the games you're having to play now.  Beyond that, I have two comments:

- I am actually surprised to see that Eidos doesn't have the rowSums() and colSums() functions that R has; I thought I had added those by now.  It seems like that would probably help your algorithm somewhat.  I'd be happy to add them if you think they'd be useful to you; if so, please open a new issue on SLiM on GitHub, so that we can connect about it there.

- The next big project on my plate is adding intrinsic support for multiple traits, pleiotropy, etc., which should make this model of yours a cake walk.  If you'd like to hear more, and perhaps be an early tester of it (although I'm not at that stage yet), please contact me off-list at bha...@mac.com for more info.

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University


Joshua Schraiber wrote on 10/6/25 8:27 PM:

--
SLiM forward genetic simulation: http://messerlab.org/slim/
Before posting, please read http://benhaller.com/slim/bugs_and_questions.html
---
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/5e8ab68b-6b4a-4a49-aa38-e62f518bdfecn%40googlegroups.com.

[Ben Haller]

A quick followup on this thread.  I've just added rowSums() and colSums() functions to Eidos.  For rowSums() it is about 12x faster than using apply(); for colSums() it is about 8x faster.  (The difference is to the locality of memory access patterns, I would guess.)  So you might try using these, Joshua.  Your code will probably still be painfully slow, but perhaps a bit less painfully slow than before.  :->  See my previous reply for other recommendations.

This addition is now in the master branch of SLiM on GitHub, which you can build following the instructions in chapter 2 of the SLiM manual; they will be released in the next release of SLiM (which is presently unplanned).

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

[Joshua Schraiber]

Awesome, thanks for that! I'll try to get this installed soon. Really appreciate it.

You received this message because you are subscribed to a topic in the Google Groups "slim-discuss" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/slim-discuss/w2-Q4B5jxhU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to slim-discuss...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/b7648992-13a4-44c0-b813-cc9f44987a15n%40googlegroups.com.

[Joshua Schraiber]

A quick follow-up regarding downstream stuff from this model: it looks like saving the tree sequence *does not* save anything that's stored by setValue to the mutations? At least it doesn't look like it's living in the mutation objects in tskit. 

[Peter Ralph]

That's right. What is stored in the tree sequence is described here (and in the SLiM manual):

 https://tskit.dev/pyslim/docs/latest/metadata.html

and extending this to save arbitrary inforrmation that's attached to the various objects would not be efficient. To pass along information like this, put it in top-level metadata. Something like this (untested!):

sim.treeSeqOutput("out.trees", metadata=Dictionary(

   "mut_ids", sim.mutations.id,

   "effects", sim.mutations.getValue("effects"))

);

Then in python you get this by doing

mut_ids = ts.metadata['SLiM']['user_metadata']['mut_ids']

effects = ts.metadata['SLiM']['user_metadata']['effects']

and then you'd match them up to the mutations in the tree sequence using the mutation IDs that are stored in the derived state of each mutation.

--peter

---

From: slim-d...@googlegroups.com <slim-d...@googlegroups.com> on behalf of Joshua Schraiber <jgsch...@gmail.com>
Sent: Wednesday, October 8, 2025 4:37 PM
To: slim-discuss <slim-d...@googlegroups.com>
Subject: Re: Efficient simulation of pleiotropy

 

[Ben Haller]

That is correct, and documented.  If you have information that you want to persist, beyond what is saved by SLiM, you can assemble a metadata dictionary and pass it to treeSeqOutput().

Cheers,
-B.

Benjamin C. Haller
Messer Lab
Cornell University

Joshua Schraiber wrote on 10/8/25 7:37 PM:

To view this discussion visit https://groups.google.com/d/msgid/slim-discuss/d2041f4c-178a-4bb9-9265-bb6a33d54292n%40googlegroups.com.