offsets and Y-like chromosomes

[petrelharp]

Consider this simulation of a simple quantitative trait on a Y chromosome.

initialize() {
   setSeed(35);
   initializeSLiMOptions(keepPedigrees=T);
	
	initializeSex();
	
	quant1T = initializeTrait("quant1T", "a", 0.0, 0.0, 0.0, directFitnessEffect=F);
	initializeMutationType("m1", 0.4, "e", 0.05);
	initializeGenomicElementType("g1", m1, c(1));
	
	symbols = c("Y");
	types = c("Y");
	names = c("Y");
	ids = seqAlong(names);
	lengths = rdunif(size(ids), 1e3, 2e3);
	
	for (id in ids, symbol in symbols, length in lengths, type in types, name in names)
	{
		initializeChromosome(id, length, type, symbol, name);
		initializeMutationRate(1e-3);
		initializeRecombinationRate(1e-4);
		initializeGenomicElement(g1);
	}
}

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

1: late() {
    sim.demandPhenotype(sim.subpopulations);
    inds = p1.individuals;
    f = inds[inds.sex == "F"];
    catn(sum(sim.substitutions.effectSize));
    catn(f.quant1T);
}

10 late() {
    sim.demandPhenotype(sim.subpopulations);
    sim.simulationFinished();
}

This is printing the phenotypes of the females, whch should all be 0.0, having no offsets. However, I get the following, which the first number being the sum of all substitutions effects, and the vector being female phenotypes:

0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.211879
0.423758 0.423758 0.423758 0.423758 0.423758 0.423758 0.423758 0.423758 0.423758 0.423758
0.211879
0.423758 0.423758 0.423758 0.423758 0.423758 0.423758 0.423758 0.423758 0.423758 0.423758
0.211879
0.423758 0.423758 0.423758 0.423758 0.423758 0.423758 0.423758 0.423758 0.423758 0.423758

As you can see, the baselineOffset is applying to females as well, but AFAICT it should not?

This is another wrinkle with the baselineOffset. I haven't thought hard about this, but won't this be solved by having more than one baseline offset? If there was one baselineOffset per chromosome, then individuals not carrying that chromosome could simply not include that term. (And tangentially, if each chromosome had its own diploidBaselineOffset and hemizygousBaselineOffset then that would remove the restriction that hemizygous dominance coefficients by 1.0, I think?)

Otherwise, I think we just can't use baseline accumulation with Y-like chromosomes?

ps. I caught this since my calculation of traits on the tree sequence didn't match SLiM's values.

[bhaller]

Ah, interesting. Yes, this slipped through the cracks; good catch! (Comprehensive unit tests FTW once again! I've got unit tests in SLiM that should have caught this as well, actually, so one question is why those unit tests didn't raise for this. So, redundancy in test regimes FTW, too! :->)

The design is that every individual gets the baseline offset, and I don't want to change that; I'm not sure what it would even mean for it to be per-chromosome, since situations like the Y can arise for type A or H chromosomes as well, in situations like haplodiploidy. Having to loop through chromosomes and check for null haplosomes in each individual just to calculate the baseline offset for individuals would be a bit time-consuming, and also complicated to document and explain; and calculation of trait values is a core bottleneck that needs to be as fast as possible. Having a variable number of baseline offsets (depending on the number of chromosomes) would also make lots of other things annoying/inefficient; the baseline offsets would have to be kept in a malloced vector (extra pointer dereference to access them), it would make for more complication in the metadata schema, etc. So I don't think more than one baseline offset is the solution.

Essentially the same problem arose already in the context of hemizygosity. If the hemizygous dominance coefficient for a mutation type in use is not 1.0, baseline accumulation will raise an error if a substitution occurs because it would not work correctly for the hemizygous case (since, again, every individual gets the baseline offset). Here's an example of that:

initialize() {
	initializeSex();
	initializeTrait("trait1", "a");
	initializeMutationType("m1", 0.5, "f", 0.0);
	m1.setDefaultHemizygousDominanceForTrait(NULL, 0.8);
	initializeGenomicElementType("g1", m1, 1.0);
	
	initializeChromosome(1, 1e5, type="X");
	initializeGenomicElement(g1, 0, 99999);
	initializeMutationRate(1e-7);
	initializeRecombinationRate(1e-8);
}
1 early() {
	sim.addSubpop("p1", 500);
}
20000 late() { }

If you run this model, SLiM doesn't raise an error right away (after initialize()) because it doesn't yet know that the configuration will lead to a problem in practice. But when a mutation actually fixes and gets substituted, an error gets raised:

ERROR (Species::DoBaselineAccumulationForSubstitution): baseline accumulation cannot be enabled for trait 'trait1', because a substitution has a hemizygous dominance coefficient other than 1.0 for that trait. The effect of the changed baseline offset would therefore not match the effect of the original mutation, making baseline accumulation invalid. Either (1) hemizygous dominance coefficients must be 1.0 for the trait for all mutations, (2) baseline accumulation must be turned off for the trait, or (3) substitution must be disabled, with convertToSubstitution=F, for all mutation types where the hemizygous dominance coefficient is not 1.0 for the trait.

(I notice that this error can be made more narrow than it is; since the mutations in this model are all neutral, substitutions of those mutations would not affect the baseline offset. It ought to check that the effect size for the trait in question is actually non-zero, and raise only in that case. So I'll fix that in passing here, too.)

This mechanism should just be extended to the Y chromosome and all other similar cases, I think. Specifically, if a mutation fixes and gets substituted, and baseline accumulation is enabled for a given trait, and the mutation is non-neutral for that trait, THEN IF the mutation is on a chromosome that is known to be absent in some individuals an error should be raised. That would include the Y, as well as chromosome types W, FL, and ML.

The loose end is chromosome types A and H in models like haplodiploidy, etc.; they are allowed to be hemizygous or completely absent to support such models, but SLiM has no way of knowing whether that might potentially happen at some future point in time. But that – having all null haplosomes for a type A or H chromosome – can happen only as a result of addRecombinant() or addMultiRecombinant(); that's the only way to generate individuals with unusual haplosome configurations like that. (Well, addCloned() would also, but only by cloning an individual previously made with those other methods.) So I can (a) add a flag, per chromosome, that says "baseline accumulation occurred for a non-neutral substitution on this chromosome", that gets set when that substitution occurs, which is easy to check for, and also (b) add a flag, per chromosome, that says "this chromosome has been seen to be completely absent (all null haplosomes) in an individual", and set that flag in addRecombinant() and addMultiRecombinant() if that arises, which is easy to check for. If flags (a) and (b) both get set, then there is a potential violation of this requirement, and an error gets raised similarly to the hemizygous case that I described above. A design along these lines should be straightforward, and should add close to zero overhead.

Finally, note that this additional requirement would not be particularly onerous for the user to handle. If they want baseline accumulation in their model, and they also want a Y chromosome, then they just need to (a) use a separate mutation type for mutations on the Y, and (b) turn off substitution for that mutation type with convertToSubstitution=F. Substitution then won't occur on the Y, and so baseline accumulation won't occur for mutations on the Y, and so they won't get an error. As for the error message in the hemizygous case, shown above, the error message for the new case you've discovered could also spell out the options for avoiding the error. They won't get the speedup of substitution and baseline accumulation for their Y mutations; but that's not the end of the world. The Y is usually short anyway :->.

I'll do this as soon as I'm done with #536. Short on time right now since my parents are visiting, but today is a stay-at-home day so perhaps I'll make some progress.

[petrelharp]

That sounds fine, but two notes:

• If you had a femaleBaselineOffset and maleBaselineOffset then that would suffice for everything except for the scripts that make an individual with less chromosomes (eg the haplodiploid script)?
• These offsets wouldn't necessarily need to be user-visible: you could keep the user-visible (single) baselineOffset meaning the thing provided by the user when initializing the trait; they can recover the baseline offset by something like sum(sim.substitutions.effectSize) anyhow?

I should move/copy those tests over to the treerec/tests.

[petrelharp]

Another note: I assume this error will not be raised for the default trait, and so aren't there situations with a Y chromosome for which the trait value in individual metadata won't match the value you'd calculate from the tree sequence? This is harmless, because fitness is relative, of course - I'm just mentioning it in case it makes you think of an actual problem.

[petrelharp]

Hm, baseline accumulation is off for the default trait, so that's moot.

[bhaller]

Right, baseline accumulation is off for the default trait for exactly this sort of reason – maximum backward compatibility etc.

Yes, a male/female case would cover a lot, but it sounds more complicated to me, both conceptually and in the code – especially since it doesn't solve the whole problem, so the other solution that I proposed above would still have to also be in the mix. I prefer to keep the stake that there is a single baseline offset, and it gets added to everyone. Anything else and I have to rewrite section 1.7 and chapter 13, redesign the flowchart diagram that shows how trait calculations work, etc. (Because I do think this would have to end up being user-visible, because the baseline offset itself is user-visible and read/write, etc. The definition you propose, "meaning the thing provided by the user when initializing the trait", is different from the definition now.)

[petrelharp]

Okay - it's worth pondering a bit more, I think, because these caveats create extra complexity and pitfalls? For instance:

Finally, note that this additional requirement would not be particularly onerous for the user to handle. If they want baseline accumulation in their model, and they also want a Y chromosome, then they just need to (a) use a separate mutation type for mutations on the Y, and (b) turn off substitution for that mutation type with convertToSubstitution=F.

I guess I'm not sure I'd describe this as "not particularly onerous" because: (a) users will get a easy-to-deal with error that says "turn off substitution for this mutation type"... but then they have to go figure out how to do that and also probably try to figure out what that the ramifications of that are; and (b) having a different mutation type for the same biological type of mutation is going to create complexities downstream (and also upstream - I'm thinking of stdpopsim), so a more comon fix will probably be just to turn off substitution - but I guess I don't have a good sense of what the performance implications of that are?

[bhaller]

Okay - it's worth pondering a bit more, I think, because these caveats create extra complexity and pitfalls?

OK. So you're proposing a separate male/female baseline offset to take care of the differences for all of the sex-dependent chromosome types? For the haplodiploid case and other similar things where the user is using add[Multi]Recombinant() to do unusual things, that would still be handled as I describe above, and the user would have to figure out what to do about the error message they get. And for the case of a hemizygous dominance coefficient other than 1.0, similarly they would get an error and have to figure out what to do. Those two concessions are necessary because we haven't figured out a better way to handle those cases; they just seem kind of incompatible with baseline accumulation.

I'm OK with this, since you feel strongly, and I do see the benefits of it. It does deal with the 95% case, and it'd be nice to be able to keep substitution turned on for that case and not have to have a separate mutation type and so forth. It feels like the best route forward, if we do this, is to make public and user-visible the distinction that we now have in the top-level metadata, between the user-defined baseline offset and the additional baseline offset that comes from substitutions. Maybe there is utility for both of these to have separate versions for the two sexes; it makes good biological sense for males and females to have different baseline offsets, for cases like sexual dimorphism in body size for example. (Presumably such dimorphism is in some way governed by effects of sex-linked mutations or some such, but the point of baseline offset is to let the user just say "Here's the baseline when no mutations are present" without having to explicitly model the genetics underlying that baseline.)

So Trait would change to have:

baselineOffset <-> float$
baselineOffsetM <-> float$
baselineOffsetF <-> float$

These would provide the user-defined baseline offsets for hermaphrodites, males, and females respectively. Accessing a sex-dependent property in a non-sex model, or vice versa, would raise. These are read/write; the user is free to change the baseline offset if they wish, to reflect perhaps an environmental change that shifts the baseline.

And Trait would also have:

substitutionOffset => float$
substitutionOffsetM => float$
substitutionOffsetF => float$

These would similarly provide the substitution offsets, which are now conceptually different from the baseline offset. (Probably the name of "baseline accumulation" changes to "substitution accumulation".) They are read-only. In a fresh simulation they start at zero (no substitutions); whenever a substitution is made (including when reading in a .trees file, if the various conditions for substitution are met, including convertToSubstitution being true for the mutation type in question, and substitution accumulation being enabled for the focal trait) they accumulate the effect of the mutation into that property. They get written out by SLiM into top-level metadata, but do not get used by SLiM on read; they get recalculated instead, as per the design we have worked out now.

It's a bit clunky having three properties for each of these, but I don't think there's a better solution really, and in practice it shouldn't matter much; people will rarely actually interact with there properties in their script beyond setting the initial value of the baseline offset in initializeTrait(). Speaking of which, that should probably be removed from initializeTrait(), since it isn't clear how to make the API nice given that there might be one (hermaphroditic) value or two (sex-linked) values; better to just let the user set the baseline offsets they want with the properties above, after the initializeTrait() call?

What else? The diagram about trait calculation will change to show both the baseline offset and the substitution offset feeding in. All calls to initializeTrait() everywhere will change to not take a baseline offset. All doc related to baseline accumulation will shift to the new terminology/paradigm. The top-level JSON tree-seq metadata will change to persist all of this; for the property values that don't apply to the model (hermaphroditic values in a sexual model, and vice versa) the value in the metadata will be NAN. All trait calculation code (in SLiM and pyslim) will shift over to using the correct version of the baseline offset and the substitution offset.

This is a LOT of change, so let's make sure we're on the same page and this design will make us both happy. What are your thoughts on the proposal above?

[petrelharp]

And for the case of a hemizygous dominance coefficient other than 1.0, similarly they would get an error and have to figure out what to do.

Note that this would be (usually) dealt by the male/female thing.

[petrelharp]

Probably the name of "baseline accumulation" changes to "substitution accumulation".

I like that!

that should probably be removed from initializeTrait()

That seems okay.

[petrelharp]

Okay, big-picture thoughts: I wouldn't say I feel strongly? And I don't have as good an idea of the amount of work as you do. But my feeling is that this would clean things up enough (particularly if you agree it also takes care of the hemizygous coeffienct not being 1.0) that I think it's worth a day or so? Not worth a week?

[bhaller]

I think it is the right path, the more I ponder it. And hmm, I think it does take care of the hemizygous dominance not being 1.0, but only for the case where the hemizygosity is the result of being sex-linked; it takes care of it for the X and the Z, in other words. But that is helpful, for sure. So it takes care of all of the issues with baseline accumulation that are related to sex and hemizygosity. And it doesn't take care of any of those problems for cases related to haplodiploidy, etc.; those will be checked for and produce a raise if there is a problem. I think I will implement this when I have time. Thanks for the suggestion, I am convinced!

[petrelharp]

Note that H and Y also will work, I think.

And gee, naively, wouldn't even the haplodiploid example in the docs work with this proposal if the chromosomes were all set to be X chromosomes?