python trait tests

[petrelharp]

Here's python code to test for trait consistency. Currently, I'm getting failures only for these two:

FAILED test_trait_consistency.py::TestTraits::test_traits_consistency[recipe_WF_adds.slim] - assert False
FAILED test_trait_consistency.py::TestTraits::test_traits_consistency[recipe_WF_no_substitutions.slim] - assert False

This is confusing to me - I should have more! - because I seem not to be triggering #671? And, why is "WF no substitutions" failing but not "nonWF" (which has no substitutions)?

There's probably some dumb errors, is why. But there's a lot going on here.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 76.77%. Comparing base (0a0d915) to head (95337db).

Additional details and impacted files

@@              Coverage Diff               @@
##           multitrait     #672      +/-   ##
==============================================
- Coverage       76.79%   76.77%   -0.03%     
==============================================
  Files             117      117              
  Lines           78866    78866              
  Branches        14320    14316       -4     
==============================================
- Hits            60565    60549      -16     
- Misses          18301    18317      +16     

see 10 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[petrelharp]

Okay, I know (part of?) what's going on:

Options for accumulateBaseline and convertToSubstitution:

1. accumulate=T, convert=T: mutations that fix count towards the trait
2. accumulate=T, convert=F: mutations that fix count towards the trait
3. accumulate=F, convert=T: mutations that fix do not count towards the trait
4. accumulate=F, convert=F: mutations that fix count towards the trait

In the current code, I was including fixed mutations in trait computation if accumulate=T. But if we're not converting mutations, this is not correct.

And, from the tree sequence I have no way of knowing whether (a) a given mutation was converted to a substitution, or (b) whether a given mutation type has convertToSubstitution turned on or not.

For testing purposes, I could record the convert choice in metadata and so do these tests.

I guess this means that you could save out a tree sequence with accumulate=F, then reload it, and get different trait values, but only if you change the convertToSubstitution property when you reload? I haven't tested that. But, that seems like not a problem.

I'm not sure if the inability to compute traits - in all situations, without extra information from the script - from the genotype information in the tree sequence is a problem. I suppose not.

[petrelharp]

Hah! I got the logic right for convert/accumulate. Still one failing, and that's the one where I add mutations: recipe_WF_adds.slim.

[petrelharp]

Okay, there's something funny with addDrawnMutation and substitutions. In the recipe_WF_adds.slim script (reproduced here):

initialize()
{
    setSeed(23);
    initializeSLiMOptions(keepPedigrees=T);
    initializeTreeSeq(timeUnit="generations");
	source("init.slim");
	 initializeSex();
    defineConstant("L", 100);
    initializeChromosome(1, L, "A", "A");
    initializeMutationRate(1e-2);
    initializeMutationType("m1", 0.5, "f", 0.01);
    initializeMutationType("m2", 0.5, "f", 0.01);
   // m2.convertToSubstitution = F;
    initializeGenomicElementType("g1", m1, 1.0);
    initializeGenomicElement(g1, 0, L-1);
    initializeRecombinationRate(1e-2);
}

1 early() { 
    sim.addSubpop("p1", 10);
}
1: late() {
    for (chrom in sim.chromosomes) {
        haps = p1.individuals.haplosomesForChromosomes(chrom, includeNulls=F);
        if (length(haps) > 0) {
          sample(haps, 1 + asInteger(length(haps)/2)).addNewDrawnMutation(m2, rdunif(10, 0, L-1));
        }
    }
}

30 late() {
    saveTreeSeqTraits();
    catn("Done.");
    sim.simulationFinished();
}

if I uncomment the line that makes m2 not convert to substitutions (and m2 are the added mutations) then this passes.

This is vaguely ringing a bell? Do addDrawnMutations interact with baseline offset differently or something?

[bhaller]

Okay, I know (part of?) what's going on:

Options for accumulateBaseline and convertToSubstitution:

1. accumulate=T, convert=T: mutations that fix count towards the trait
2. accumulate=T, convert=F: mutations that fix count towards the trait
3. accumulate=F, convert=T: mutations that fix do not count towards the trait
4. accumulate=F, convert=F: mutations that fix count towards the trait

For (2) and (4), it's moot because if convert=F, no mutations get turned into substitutions. They stay mutations, and they continue to count towards the trait just like other mutations. Not sure if we're on the same page there or not.

In the current code, I was including fixed mutations in trait computation if accumulate=T. But if we're not converting mutations, this is not correct.

So here I'm puzzled. As I commented above, if convert=F, mutations stay mutations, and continue to count; the value of accumulate then doesn't even matter, since substitution doesn't happen. But you say here "I was including fixed mutations in trait computation... this is not correct." In what way is it not correct?

And, from the tree sequence I have no way of knowing whether (a) a given mutation was converted to a substitution, or (b) whether a given mutation type has convertToSubstitution turned on or not.

For testing purposes, I could record the convert choice in metadata and so do these tests.

I'm confused about this, too. Why do you even need to know whether a mutation was converted to a substitution or not? On the python side, they're all just mutations, and they all count; there is no distinction between mutations and substitutions on the python side. I thought you weren't using the "substitution offset" on the python side at all.

I guess this means that you could save out a tree sequence with accumulate=F, then reload it, and get different trait values, but only if you change the convertToSubstitution property when you reload? I haven't tested that. But, that seems like not a problem.

Yes, to get correct results it is required that the simulation be configured in the same way at save and load. SLiM doesn't check every possible way that that could be violated (because then we'd have to write everything out to metadata, just to be able to do that check), but it is required, and if the user goes off-piste the consequences are Somebody Else's Problem.

I'm not sure if the inability to compute traits - in all situations, without extra information from the script - from the genotype information in the tree sequence is a problem. I suppose not.

I am puzzled as to why this is not possible; I'm not following your logic above. Clarify please?

[bhaller]

By the way, you asked somewhere (Slack?) why case 3 is allowed, since it would appear to produce incorrect trait values; a mutation gets substituted, and its effects just go away. The reason is basically for backward compatibility. This is precisely how things worked before SLiM 6 when simulating quantitative traits, since there was no built-in support for them, just sumOfMutationsOfType(). Back in that paradigm, it was recommended to always set convertToSubstitution=F if you were modeling quantitative traits, so that this wouldn't bite you. But, some clever people, like Vítor Sudbrack, realized that if they added the effects from substitutions into their trait values they could still use convertToSubstitution=T and things would work correctly. So there's code out there that does exactly that. The default trait (which, again, is there for backward compatibility) has accumulate=F; and old models that use sumOfMutationsOfType() with the default trait will continue to work; and if they also set convertToSubstitution=T and do their own substitution effect accumulation, that, too, will continue to work. The goal was not to break those existing models. And I think you can/should assume that if the user is in case 3, they are in fact doing their own substitution effect accumulation; just assume that, because if they're not, their model is totally broken anyway. Maybe issue a warning, if you want to.

But pondering this, probably that combination should only be legal when using the default trait...? In models that use initializeTrait(), and are thus in the new paradigm, probably having convert=T but accumulate=F should be an error. Or at least a strongly worded warning; maybe SLiM shouldn't make such a strong assumption that it knows better than the user. But a strongly worded warning seems right (and again, I think on the Python side you can/should assume that they are doing their own substitution effect accumulation in that case). Do you agree?

So maybe this case 3 is the source of a lot of your woes above? Sorry about that, if so!

[bhaller]

Okay, there's something funny with addDrawnMutation and substitutions. In the recipe_WF_adds.slim script (reproduced here):
...
if I uncomment the line that makes m2 not convert to substitutions (and m2 are the added mutations) then this passes.

This is vaguely ringing a bell? Do addDrawnMutations interact with baseline offset differently or something?

No, I don't think so...? So, hmm. Looking at your code:

1: late() {
    for (chrom in sim.chromosomes) {
        haps = p1.individuals.haplosomesForChromosomes(chrom, includeNulls=F);
        if (length(haps) > 0) {
          sample(haps, 1 + asInteger(length(haps)/2)).addNewDrawnMutation(m2, rdunif(10, 0, L-1));
        }
    }
}

Parsing what this code is doing... It loops over the chromosomes and, for each chromosome, gets the vector of all haplosomes. Then it calls sample() to sample from that vector of haplosomes, selecting about half of them. And then it calls addNewDrawnMutation() on each of those haplosomes (in a vectorized method call) and asks each one to add 10 new drawn m2 mutations at 10 random positions that might or might not be the same.

So I guess the question is, what's special about this code? What is it doing that pushes a weird button, either in SLiM or in your Python code? My first guess is that the problem occurs when two of the draws from rdunif(10, 0, L-1) are actually the same position. I'm not sure, without delving into the code for addNewDrawnMutation () to check (and it is surprisingly complicated code!), what exactly happens in that case. In principle, I guess one should be added and then the other should be added, each invoking the stacking policy that is in effect (which looks to just be the default of "s", stack), and the derived state of the second one should include the first one. Should they be recorded as occurring at exactly the same time, or does tskit require their times to be different such that we ought to have some tiny fractional time offset? Anyhow, I'd suggest trying sample(0:(L-1), 10, replace=F) to draw 10 unique positions and see if that fixes the problem; and if it does, then the positional collision is in some way causal.

Note that the SLiM core no longer produces this situation, ever; it uniques the mutation positions that it draws, specifically to avoid this kind of situation because it was giving us problems! I don't remember the details of the problems it caused, but it seems likely to me that the same problems are arising here. If so, then perhaps addNew[Drawn]Mutation() should check for non-unique positions and throw an error...? But the user could still add mut A and then, in a separate call, add mut B on top of it, so maybe that check is not sufficient; maybe when adding each mutation we'd need to check "is any existing mutation at this position recorded as having now as its time of origin" and throw an error on that?

Anyhow, let's first figure out if positional collision is indeed the problem, then figure out why that's a problem, then figure out what to do about it. :->

[petrelharp]

But pondering this, probably that combination should only be legal when using the default trait...?

I agree about the strongly worded warning!

[petrelharp]

In the current code, I was including fixed mutations in trait computation if accumulate=T. But if we're not converting mutations, this is not correct.

So here I'm puzzled. As I commented above, if convert=F, mutations stay mutations, and continue to count; the value of accumulate then doesn't even matter, since substitution doesn't happen. But you say here "I was including fixed mutations in trait computation... this is not correct." In what way is it not correct?

LOL sorry. My first pass here was "if a mutation is at frequency 1.0 then only include it if accumulate=T". This assumes that convert=T. In fact, I need to not include the effects of a mutation that I see at frequency 1.0 only if convert=T and we are not accumulating. In other words, in case (3).

[petrelharp]

Nope, uniquing the mutation location (or just applying one mutation per hap) doesn't fix the problem.

[petrelharp]

It's being very tricky to get this into a part of parameter space where either (a) it's totally minimal, or (b) the behavior doesn't depend on seed. But, this is pretty good (note it's only adding a single m2 mutation in the first generation):

initialize()
{
    setSeed(37);
    initializeSLiMOptions(keepPedigrees=T);
    initializeTreeSeq(timeUnit="generations");
	 initializeSex();
    defineConstant("L", 1);
    initializeChromosome(1, L, "A", "A");
    initializeMutationRate(0.99);
    initializeMutationType("m1", 0.5, "f", 0.01);
    initializeMutationType("m2", 0.5, "f", 0.01);
    initializeGenomicElementType("g1", m1, 1.0);
    initializeGenomicElement(g1, 0, L-1);
    initializeRecombinationRate(1e-2);
}

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

1 late() {
   hap = sample(p1.haplosomes, 1);
   loc = sample(0:(L-1), 1);
   hap.addNewDrawnMutation(m2, loc);
}

30 late() {
    catn("Done.");
    sim.demandPhenotypes(NULL);
    sim.treeSeqOutput("out.trees");
    sim.simulationFinished();
}

This has traits that differ from python. Changing the seed to 38 does not. What's the difference between those? In both those seeds, (a) there's both an m1 and an m2 mutation at the same position in a single individual; and (b) the m2 mutations are no longer around at the end of the simulation.

The trait values in metadata agree with what SLiMgui reports on the fly.

Making the number of generations smaller (tends to?) make the problem go away: for instance, at 20 total generations, the problem only occurs about half the time; here it occurs most of the time. A similar thing happens with the number of individuals.

Here's the phenotypes: first column in python, second is as reported by SLiM in metadata, third is the difference:

[[1.34207301 1.31562853 0.02644447]
 [1.33539603 1.30908312 0.0263129 ]
 [1.34160802 1.31517271 0.02643531]
 [1.34160802 1.31517271 0.02643531]
 [1.34207301 1.31562853 0.02644447]
 [1.33539603 1.30908312 0.0263129 ]
 [1.33493335 1.30862956 0.02630379]]

This is typical: across seeds, the last column is different, but all the values agree with each other in the first two decimals. (e.g., 0.026 in this seed, 0.015 in another).

[petrelharp]

Okay, the problem is that SLiM thinks that these individuals have fewer mutations than python does! Now, which ones?

[petrelharp]

And, the mutations in the first individual that python is counting that SLiM is not are 13, 13, 20, 20, 40, the first ones! These are fixed.

For reference, here's the python code I used to investigate (note that in this model, fitness is 1.01^(num_mutations/2)):

import numpy as np

node_muts = np.zeros(ts.num_samples)
for v in ts.variants(isolated_as_missing=False):
    for j, g in enumerate(v.genotypes):
        a = v.alleles[g]
        if a != "":
            node_muts[j] += len(a.split(","))

ind_muts = np.zeros(ts.num_individuals)
for j, k in enumerate(ts.samples()):
    n = ts.node(k)
    if n.individual >= 0:
        ind_muts[n.individual] += node_muts[j]
print(node_muts)
print(ind_muts)
print(1.01**(ind_muts/2))
print(ts.tables.individuals.metadata_vector(["per_trait", 0, "phenotype"]))

v = next(ts.variants(samples=ts.individual(0).nodes, isolated_as_missing=False))
m = list(map(int, (",".join(v.alleles[g] for g in v.genotypes)).split(",")))
m.sort()
m

And, here's the mutations carried by the first individual according to the tree sequence:

[13, 13, 20, 20, 40, 44, 49, 58, 63, 72, 81, 88, 95, 103, 108, 112, 115, 129, 131, 145, 146, 156, 160, 170, 171, 183, 187, 195, 206, 212, 213, 227, 231, 243, 246, 252, 262, 265, 266, 279, 281, 294, 302, 310, 315, 320, 321, 334, 340, 349, 358, 367, 368, 379, 384, 393, 394, 401, 402]

[petrelharp]

OKAY THE PROBLEM IS IN PYTHON

You know what my code thinks the frequency of SLiM mutation 13 is? 0.9999999999999997. (Thus, not fixed.)

[petrelharp]

Sorry for all the noise there! That was very annoying. The problem was that I was adding up the frequencies of all the variants that carried mutation 13 to get the frequency of SLiM mutation 13; however, floating-point error. Changing to adding up the counts fixed the problem.

[bhaller]

Yeah, can't count frequencies due to roundoff. I had that bug in SLiM somewhere, at one point. Excellent detective work!

[petrelharp]

Now I'm trying to get it to use development pyslim...

[petrelharp]

what the heck, I swear this passed the tests before I squashed those commits