gregorgorjanc
commented
6 months ago I now looked into this.
First a bit of background.
We initially implemented Ne trends for cattle based on SNP array estimates and some basic cattle genome biology params. We got about ~50K sites per chromosome in total, which seemed enough.
Then we got a reviewer saying in one manuscript that we should be getting ~1M sites per chromosomes. To achieve that I implemented the Macleod et al. (2013) model as the best demographic model I have found. That bumped sites up qiute a bit, but it was still too low. I was at the same time reading a human genetics paper on de-novo mutation rate (~2.5e-8) and once I used that rate, we got the suggested ~1M sites.
As I was implementing the stdpopsim BosTaurus model, I was advised that we should not be changing model parameters freely because estimated model parameters depends on some inputs and that it is best to faithfully represent the published model with all inputs and estimates.
Currently we have these values:
- cattleChrSum = 3e9
- cattleChrBp = cattleChrSum / 30 # 1e8
- recRate = 1.00e-8
- genLen = recRate * cattleChrBp # 1.00
- mutRate = 2.5e-08 # from human literature
- eyeballed Ne from Macleod et al. (2013) figure 4B
- this gives me ~660K sites with 100 samples (the number grows with more samples)
I am proposing this change to be fully in line Macleod et al. (2013):
- cattleChrSum = 2.8e9 # ~7% smaller genome
- cattleChrBp = cattleChrSum / 30 # 0.93e8
- recRate = 9.26e-9 # from Macleod et al. (2013)
- genLen = recRate * cattleChrBp # 0.86
- mutRate = 1.20e-8 # from Macleod et al. (2013)
- exact Ne from Macleod et al. (2013) supplementary table S1
- this gives me ~240K sites with 100 samples (the number grows with more samples), which makes sense as we have reduced mutation rate and shortened the genome
It should also run a bit faster;) Will push a PR for this soon.
When implementing BosTau & Holstein demography between in
runMacs()I eyeballed numbers from figure 4 in Macleod et al. paper, while I later implemented the same in stdpopsim and figured out that there is a table with exact numbers in the paper appendix - see https://popsim-consortium.github.io/stdpopsim-docs/stable/catalog.html#sec_catalog_bostau_models_holsteinfriesian_1m13