Stdpopsim: AraTha recombination rate is out by a factor of 10

[grahamgower]

_Note: this bug has been fixed. The issue is described here for posterity._

The Arabidopsis species definition uses recombination_rate=8.1e-9. Above the species definition is the following comment:

# recombination value from Huber et al 2014 MBE
# rho=200/Mb, assume Ne=124,000, rho=2*Ne*r

Indeed, the values in the comment appear to accurately reflect the Huber et al. publication. However, I calculate r = rho/Ne/2 = (200/1e6)/124000/2 = 8.064516129032259e-10.

[grahamgower]

Fixed in #527.

[andrewkern]

just seeing this-- this doesn't make sense to me as the reported recombination rates for Arabidopsis are on the order of ~4cM/Mb. e.g., https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3313057/

[grahamgower]

So we're still out by a factor of 2? And this isn't a haploid versus diploid thing? Should we just use the mean rates as calculated from the Salome map?

$ cat aratha_rr.py
import stdpopsim

aratha = stdpopsim.get_species("AraTha")

for chrom in aratha.genome.chromosomes:
    contig = aratha.get_contig(chrom.id, genetic_map="SalomeAveraged_TAIR7")
    rr = contig.recombination_map.mean_recombination_rate
    print(chrom.id, rr)

$ python aratha_rr.py
chr1 3.511099648074139e-08
chr2 3.8060257952620545e-08
chr3 3.598006781811692e-08
chr4 4.0758147612457255e-08
chr5 3.642943789368259e-08

[grahamgower]

Oh, shit, that's a factor of 50!

[grahamgower]

From Huber et al. 2014:

In the simulations, we set recombination rate in A. thaliana to a value of 2.5 cM/Mb, which together with the ancestral population size of 124,000 results in an effective population recombination rate ρ of about 200 (2Ner) per 1 Mb, assuming an inbreeding rate of 97% (Platt et al. 2010). We checked whether this amount of recombination leads to similar LD levels as in the data by estimating the effective recombination rate ρ for both model and data, in southern and northern Sweden separately (fig. 3). ρ was estimated in 1-Mb windows by filtering SNPs for allele frequency of at least 5% and fitting a curve according to (Hill and Weir 1988) to the plot of r2 versus distance in bp. The single parameter of this curve is ρ and is estimated by a nonlinear least squares approach in R.

Salomé calculated recombinations from F2 crosses. So maybe the difference is caused by adjusting for the quoted 97% inbreeding rate?

[andrewkern]

seems there is a lot of variation in recombination rate in these plants! every paper i look at has a different number.... Graham are you thinking that inbreeding affects Ne and so that could be why Huber et al's number is strange?

[apragsdale]

I'm not sure how the inferences were made, but at least for simulating self-compatible populations there is a straight-forward rescaling of mutation and recombination parameters based on the selfing rate (https://www.genetics.org/content/154/2/923.short). Slim can directly allow any selfing rate, I imagine, but for the msprime engine this could be a potential feature for a future release to allow selfing?

[grahamgower]

Actually, I was thinking that recombination during selfing is just not observable (so the lower rate would match the empirical data). In any event, it's not clear what Huber et al.'s adjustment was.

[jeromekelleher]

Do we agree that the currently released version accurately gives the rate estimated by Huber et al?

[andrewkern]

right. that is the germane question and i think the answer is yes. i think.

[jeromekelleher]

I think we can close this issue?