It would be nice to be able to run simulations using different mutation or recombination rates than the defaults or options available for each species. Currently, running a simulation such as
engine = stdpopsim.get_engine("msprime")
ts = engine.simulate(model, contig, samples, msprime_model="dtwf")
works perfectly fine. But if I want to, say, double the mutation rate, I can't quite figure out how to do so. Setting contig.mutation_rate = my_mut_rate gives me a FrozenInstanceError, and the only solution I've come up with so far is to take the output tree sequence, remove all mutations, and then re-mutate the tree sequence with my preferred mutation rate.
In general, it would be nice to make it easier to run simulations with mutation rates and recombination maps that aren't supported as defaults within stdpopsim, that instead the user locally defines, and we should then also document how to do so.
Probably an easier way around this is make a new contig object. E.g.
new_contig = stdpopsim.Contig(
mutation_rate=my_mut_rate,
recombination_map=old_contig.recombination_map,
genetic_map=old_contig.genetic_map)
But maybe we should change stdpopsim.Contig to not be frozen?
Oh, and almost nothing in python is really immutable. You can mutate frozen attr classes using this hack:
object.__setattr__(contig, "mutation_rate", my_mut_rate)
But don't tell anyone I suggested that. :P
But maybe we should change stdpopsim.Contig to not be frozen?
Yeah, no good reason for us to make this frozen. We should just turn this off.
Once this gets sorted out, someone should document it. @apragsdale, since you're doing it already, want to throw a simple example together?
Yes, I'm happy to @petrelharp. Though my current workaround doesn't do as Graham suggested, making these attributes non-frozen seems like a good idea.
This might be a good chance to also wrap in #501, where we want to allow simulations under the mutation rate that the model was inferred with. I think this shouldn't be hard to put together, because in the original papers for a given model, they should have cited the mutation rate that they used in their inference.
Logistically, would this mean that the contig and demographic models would need to be able to communicate? Could imagine something like, "if demographic model has a set mutation rate, use that rate, and if not use the species' default rate." And if the user passes some argument mutation_rate = other_mut_rate, then the simulation would override the default rate and use the input rate instead?
I think just taking out frozen=True from the contig class is a good first step @apragsdale - the other stuff is complicated and I'd like to push it out until after the next release.
There are currently three classes marked with frozen=true (Chromosome, Contig, and Engine). I can't really think of any reason they need to be frozen, so maybe they should all be unfrozen?
Ok, @jeromekelleher, that sounds reasonable. That can be pushed off until after this next release, and we can just make this minor adjustment instead. I'll put together a PR that does that and also adds to the documentation how to simulate with your own mutation and recombination rates/map.
All of these attrs classes have been unfrozen now, btw.
@hyanwong added docs on this in #632 .
Most helpful comment
Probably an easier way around this is make a new contig object. E.g.
But maybe we should change stdpopsim.Contig to not be frozen?