Add recombination rate to stdpopsim.DemographicModel

[gregorgorjanc]

Short description of what you'd like added:

Add recombination rate specification in demographic models in addition to in species

Rationale and considerations:

We can specify recombination rate in the species file, but demographic models for that species could have assumed different recombination rate (and mutation rate) in estimation. We can already specify mutation rate separately for species and demographic models, so it would be good to offer the same for recombination rate.

This came about in QC of the BosTau #579, where the difference between the most recent recombination estimate used for species definition and a bit older recombination rate estimate for demographic model are not too different, so the issue is not huge there, but could be much bigger in other species.

How can I help?

[gregorgorjanc]

@gregorgorjanc look at how demography-specific mutation rate was handled in #839

[gregorgorjanc]

I am looking into this since many demographic models are inferred/calibrated by assuming a specific recombination rate, maybe even a specific genetic map, which will invariably differ from the species recombination rate and genetic map we have implemented in stdpopsim.

I looked across the codebase, but found it tricky, hence writing this summary to see if I am following the structure correctly.

We have:

• Species class with:

  • genetic_maps

• Genome class with:

  • mean_recombination_rate property function (length-weighted mean recombination rate across all chromosomes - crossovers, not DSB)
  • mean_mutation_rate property function (length-weighted mean mutation rate across all chromosomes)

• Chromosome class with:

  • recombination_rate attribute (rate of crossovers, when genetic_map is not used)
  • mutation_rate attribute
  • there was a helpful docstring Note that although the recombination rate for a Chromosome is the rate of crossovers, the recombination map for a Contig gives the rate of both crossovers and gene conversion, so will differ if the gene conversion fraction is nonzero.

• Contig class with:

  • genetic_map attribute (undocumented? TODO) (I assume map of crossover rates)
  • recombination_map attribute (The recombination map for the region, that gives the rates of double-stranded breaks. Note that if gene conversion fraction is nonzero, then this map can have a larger rate than the corresponding chromosome's recombination rate, since the chromosome's rate describes only crossovers)
  • mutation_rate attribute
  • recombination_rate is passed to some functions to build recombination_map

• DemographicModel class with:

  • mutation_rate attribute (if not provided, species default is used)

The docs and code interchange between "recombination map" and "genetic map", which I understand as they are synonyms, but Contig class distinguishes between them. Based on the above summary, this is because of the possibility of gene conversion, so "recombination map" in stdpopsim is "genetic map" combined with gene conversion, right? Should we then be a bit clearer in the docs to avoid this confusion.

Potential TODOs:

• Document genetic_map attribute in Contig class
• Clarify/consolidate recombination/genetic map language in the docs, possibly also code (use genetic map to refer to crossover rate map and recombination map to refer to crossover rate and gene conversion map combined?)
• Expand the DemographicModel class with:
  • recombination_rate (to specify which recombination rate was used in constructing demographic model, if not provided, species default is used)
  • genetic_map (to specify which genetic map was used in constructing demographic model [assuming some did?], if provided, overrules recombination_rate)
• Pass DemographicModel.recombination_rate or DemographicModel.genetic_map to engines (by feeding them into Contig)
• Expand tutorial docs on how to pass
  • DemographicModel.recombination_rate to Contig for work with Python API
  • DemographicModel.genetic_map to Contig for work with Python API

Any suggestions or comments on this plan?

[gregorgorjanc]

I have implemented the above (except genetic map in a demographic model - not sure how important that is - maybe for future).

Using that PR code I get this behaviour:

> stdpopsim BosTau \
  --chromosome 25 \
  --right 10000 \
  --output BosTau.trees \
  pop_0:1
...
Mean recombination rate: 9.26e-09
Mean mutation rate: 1.2e-08
...

> stdpopsim BosTau \
  --chromosome 25 \
  --right 10000 \
  --demographic-model HolsteinFriesian_1M13 \
  --output BosTau.trees \
  Holstein_Friesian:1
...
Mean recombination rate: 9.4e-09
Mean mutation rate: 1e-08
...

So, we now have a cleaner / more correct way to define demographic models as they are published.

Python API would work like this:

import stdpopsim
species = stdpopsim.get_species("BosTau")
contig = species.get_contig('25')
model = species.get_demographic_model("HolsteinFriesian_1M13")

print("contig mutation rate:", contig.mutation_rate)
# contig mutation rate: 1.2e-08
print("model mutation rate:", model.mutation_rate)
# model mutation rate: 9.4e-09
print("contig recombination rate:", f"{contig.recombination_map.mean_rate:.3}")
# contig recombination rate: 9.26e-09
print("model recombination rate:", model.recombination_rate)
# model recombination rate: 1e-08

contig = species.get_contig("chr25", mutation_rate=model.mutation_rate, recombination_rate=model.recombination_rate)
print(contig.mutation_rate == model.mutation_rate)
# True
print(contig.recombination_map.mean_rate == model.recombination_rate)
# True