QC for 1855 - #1858
Conversation
|
Tick the box to add this pull request to the merge queue (same as
|
|
Hm. We need to disentangle this with current ongoing changes, seems like. I'll have a look. |
More specificially: for instance,
where the number on the left is what @silastittes has calculated and on the right is what's in @srgibbon17's PR. And gee, what you're doing sounds exactly like what I thought @srgibbon17 was doing? @srgibbon17, shouldn't just doing this and multiplying by (1-sigma)/(1+sigma) do the trick? (BTW, you should put this code in comments in the test file, @silastittes!) (Also perhaps add a test for the recombination rate at a particular point on some chromosome?) |
|
Oh, I see: @srgibbon17 has and @silastittes has Which one is right? |
|
I think What @silastittes has with |
|
Separate from that, I think the estimated mean recombination rates might be subtly different. @silastittes built Specifically, I calculated the average rates using One motivation for calculating the mean recombination rates this way is that the chromosome wide rate should match the result from calculating the mean recombination rate from the genetic map. So, for example, this test should pass with how I calculated things (since the same 0 recombination rate extension is applied in both cases) |
|
That makes sense. We should cite the Thinking about how to do the mean rate calculation for the test, and adding a local rate test now. |
No need to cite Caicedo. I read Nordborg et al, and everything we need is pretty clearly explained there. We don't need to cite someone else for telling us how to read Nordborg. And: can we make sure we're doing this right? (Which means carefully explaining things?) @silastittes it looks like you're using Also note that Nordborg says that
Our time units ought to be generations, so I think we ought not to use rescaled units. But I think that @srgibbon17 is pointing out that we're just using an Ne that includes the effect of selfing, so one generation in our sims probably actually corresponds to |
|
@silastittes can you also provide a rationale for whichever mean rate calculation you think we should use? |
|
So two issues now
I haven't thought about this yet, but want to acknowledge we still need to figure that out. Happy to help, but was in the weeds on the second issue ...
I made a little function to extend the genetic map when the chromosome is longer. When I apply it to the current stdpopsim release it reproduces the mean recombination rates that maps = ["/sietch_colab/stittes/tmp/salome2012_liftover_maps/TAIR10_lifted_chr1.txt",
"/sietch_colab/stittes/tmp/salome2012_liftover_maps/TAIR10_lifted_chr2.txt",
"/sietch_colab/stittes/tmp/salome2012_liftover_maps/TAIR10_lifted_chr3.txt",
"/sietch_colab/stittes/tmp/salome2012_liftover_maps/TAIR10_lifted_chr4.txt",
"/sietch_colab/stittes/tmp/salome2012_liftover_maps/TAIR10_lifted_chr5.txt",]
def rate_map_zero_adjusted(hapmap_file: str, chrom: str, spp: str = "AraTha"):
map_ex = msprime.RateMap.read_hapmap(hapmap_file)
species = stdpopsim.get_species(spp)
chrom_length = species.genome.get_chromosome(chrom).length
rate = map_ex.rate.copy()
#extend to end of chromosome and remake map
if map_ex.right[-1] < chrom_length:
new_position = np.append(map_ex.position, chrom_length)
new_rate = np.append(np.nan, rate)
new_rate_map = msprime.RateMap(position = new_position, rate = new_rate)
return new_rate_map
og_maps = {str(i+1): rate_map_zero_adjusted(f, str(i+1)) for i,f in enumerate(maps)}
#Can we recreate the old means?
for chr in ["1", "2", "3", "4", "5"]:
# with the current release installed, not this PR
aratha = stdpopsim.get_species("AraTha")
contig_cat = aratha.get_contig(chr, genetic_map = "SalomeAveraged_TAIR10")
old_map = contig_cat.recombination_map
print(np.isclose(old_map.mean_rate, og_maps[chr].mean_rate)) #all trueIs the above code correct, and/or is there an easier way in the context of a test? I used this function (which I believe works) for the rate values in the test and updated the selfing correction, but our values still differ. ============================================================================ short test summary info =============================================================================
FAILED tests/test_AraTha.py::TestGenome::test_recombination_selfing_correction[3] - AssertionError: assert 1.8481460212041185e-08 == 2.0991601361354922e-09
FAILED tests/test_AraTha.py::TestGenome::test_recombination_selfing_correction[1] - AssertionError: assert 1.7819993539790216e-08 == 2.0449377650848874e-09
FAILED tests/test_AraTha.py::TestGenome::test_recombination_selfing_correction[2] - AssertionError: assert 1.9247415602980858e-08 == 2.223812153186693e-09
FAILED tests/test_AraTha.py::TestGenome::test_recombination_selfing_correction[5] - AssertionError: assert 1.837704068834849e-08 == 2.1089090575158218e-09
FAILED tests/test_AraTha.py::TestGenome::test_recombination_selfing_correction[4] - AssertionError: assert 2.0915452375673585e-08 == 2.4002257476526107e-09
FAILED tests/test_AraTha.py::TestGenome::test_selfing_correction_local - assert np.float64(2.2387625242718484e-09) == 1.950902806431916e-08
========================================================================== 6 failed, 32 passed in 5.38s ==========================================================================The current values seem to small based on a quick sanity check for chromosome 1 Does the correction get applied twice somehow? Am I'm missing something else?
If I understand correctly, I agree it makes sense to use the mean rate that stdpopsim generates as @srgibbon17 suggested. If we don't there will be a slight mismatch between the chromosome-wide recombination rate and the value generated when using the Added a local rate test as you suggested too @petrelharp. |
|
I'm sorry, I should have realized this is bumping into #1832 as well. How different are the mean rates calculated in the two ways? Also: I'm surprised that adding |
The naive msprime means are: {f"{str(i+1)}": msprime.RateMap.read_hapmap(f).mean_rate for i,f in enumerate(maps)}
{'1': np.float64(3.450876386829853e-08),
'2': np.float64(3.712478235184855e-08),
'3': np.float64(3.675733646147423e-08),
'4': np.float64(4.3132238886295825e-08),
'5': np.float64(3.533591716785435e-08)}versus {'1': np.float64(3.510538727338673e-08),
'2': np.float64(3.791740873787229e-08),
'3': np.float64(3.6408476617721134e-08),
'4': np.float64(4.120344118007697e-08),
'5': np.float64(3.620277015604653e-08)}
I think it works because the the chromosomes all end in 0.0, so adding nan shifts the rates to the next interval. Agree it's a bit hacky. Here's chromosome 5 from the catalogue versus with the little function ... left | right | mid | span | rate
-- | -- | -- | -- | --
0 | 27344 | 13672 | 27344 | nan
27344 | 271380 | 149362 | 244036 | 8.7e-08
271380 | 508092 | 389736 | 236712 | 8.3e-08
508092 | 778237 | 643164.5 | 270145 | 7.5e-08
778237 | 1166715 | 972476 | 388478 | 6.4e-08
1166715 | 1384731 | 1275723 | 218016 | 4.4e-08
1384731 | 1603468 | 1494099.5 | 218737 | 3.9e-08
1603468 | 1917138 | 1760303 | 313670 | 3.8e-08
1917138 | 2229414 | 2073276 | 312276 | 3.5e-08
2229414 | 2287469 | 2258441.5 | 58055 | 2.5e-08
⋯ | ⋯ | ⋯ | ⋯ | ⋯
23779396 | 24052883 | 23916139.5 | 273487 | 3.2e-08
24052883 | 24514178 | 24283530.5 | 461295 | 3.2e-08
24514178 | 24739811 | 24626994.5 | 225633 | 3.2e-08
24739811 | 25595063 | 25167437 | 855252 | 3.1e-08
25595063 | 25666426 | 25630744.5 | 71363 | 4.8e-08
25666426 | 25882447 | 25774436.5 | 216021 | 5.1e-08
25882447 | 26022890 | 25952668.5 | 140443 | 5.7e-08
26022890 | 26337325 | 26180107.5 | 314435 | 5.6e-08
26337325 | 26766420 | 26551872.5 | 429095 | 4.2e-08
26766420 | 26975502 | 26870961 | 209082 | 0
...
left | right | mid | span | rate
-- | -- | -- | -- | --
0 | 27344 | 13672 | 27344 | nan
27344 | 271380 | 149362 | 244036 | 8.7e-08
271380 | 508092 | 389736 | 236712 | 8.3e-08
508092 | 778237 | 643164.5 | 270145 | 7.5e-08
778237 | 1166715 | 972476 | 388478 | 6.4e-08
1166715 | 1384731 | 1275723 | 218016 | 4.4e-08
1384731 | 1603468 | 1494099.5 | 218737 | 3.9e-08
1603468 | 1917138 | 1760303 | 313670 | 3.8e-08
1917138 | 2229414 | 2073276 | 312276 | 3.5e-08
2229414 | 2287469 | 2258441.5 | 58055 | 2.5e-08
⋯ | ⋯ | ⋯ | ⋯ | ⋯
23779396 | 24052883 | 23916139.5 | 273487 | 3.2e-08
24052883 | 24514178 | 24283530.5 | 461295 | 3.2e-08
24514178 | 24739811 | 24626994.5 | 225633 | 3.2e-08
24739811 | 25595063 | 25167437 | 855252 | 3.1e-08
25595063 | 25666426 | 25630744.5 | 71363 | 4.8e-08
25666426 | 25882447 | 25774436.5 | 216021 | 5.1e-08
25882447 | 26022890 | 25952668.5 | 140443 | 5.7e-08
26022890 | 26337325 | 26180107.5 | 314435 | 5.6e-08
26337325 | 26766420 | 26551872.5 | 429095 | 4.2e-08
26766420 | 26975502 | 26870961 | 209082 | 0
|
Wait, this sounds just plain wrong (we don't want to shift rates to the next interval?) but also there is no difference between the two things you posted? |
|
Did some digging and figured out why my weird append of nan at the beginning was working. Here's the docstring from msprime for that argument: The function in the comments of the PR us updates and now looks like this: def rate_map_zero_adjusted(hapmap_file: str, chrom: str, spp: str = "AraTha"):
map_ex = msprime.RateMap.read_hapmap(hapmap_file, rate_col=2)
species = stdpopsim.get_species(spp)
chrom_length = species.genome.get_chromosome(chrom).length
rate = map_ex.rate.copy()
#extend to end of chromosome and remake map
if map_ex.right[-1] < chrom_length:
new_position = np.append(map_ex.position, chrom_length)
new_rate = np.append(rate, 0)
new_rate_map = msprime.RateMap(position = new_position, rate = new_rate)
return new_rate_mapWithout applying the selfing correction, I get the same mean_rate values as we get in the release version of The tests are still failing for unclear reasons. Do you get the same rate values (right hand side) on your local install @srgibbon17? FAILED tests/test_AraTha.py::TestGenome::test_recombination_selfing_correction[1] - AssertionError: assert 1.781967764837082e-08 == 2.0449377650848874e-09
FAILED tests/test_AraTha.py::TestGenome::test_recombination_selfing_correction[2] - AssertionError: assert 1.937839693555237e-08 == 2.223812153186693e-09
FAILED tests/test_AraTha.py::TestGenome::test_recombination_selfing_correction[5] - AssertionError: assert 1.837712630491788e-08 == 2.1089090575158218e-09
FAILED tests/test_AraTha.py::TestGenome::test_recombination_selfing_correction[4] - AssertionError: assert 2.091567275873253e-08 == 2.4002257476526107e-09
FAILED tests/test_AraTha.py::TestGenome::test_recombination_selfing_correction[3] - AssertionError: assert 1.8292173775123127e-08 == 2.0991601361354922e-09
FAILED tests/test_AraTha.py::TestGenome::test_selfing_correction_local - assert np.float64(2.2387625242718484e-09) == 1.950902806431916e-08 |
|
Hold on - I don't think that explains things? The weird appending was happening after the call to As for why the tests aren't passing, one of you has |
Doh! Thank you. |
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #1858 +/- ##
=======================================
Coverage 99.82% 99.82%
=======================================
Files 143 143
Lines 5030 5033 +3
Branches 515 515
=======================================
+ Hits 5021 5024 +3
Misses 6 6
Partials 3 3 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
Turns out the reason that the "insert a nan into the rates" code above produced the same thing as stdpopsim was two wrongs making a right: see #1859. |
|
Going back to the questions above: this is a bit of a problem since the right thing to do here is to solve #857; so what I think we'd like to do now is to do the least-bad solution and clearly document what we do (unless someone wants to deal with #857 now). Let's see:
So, we have the choice of either agreeing with the mutation rate (per generation) or with the Ne values (per adjusted-generation). Note that there is only a difference of My vote is for sticking with "time units is generations", which is what @srgibbon17 has implemented (lower recomb by a factor of |
|
I've been out for a couple days, so apologies for the slow reply. I suppose the "least-bad/best" way to calculate the mean rates might depend on how #1832 gets resolved? But, whether the mean rate includes an extension of zero or NaN on the flanking edges of the chromosome doesn't seem to affect the mean rates too much (e.g., for chromosome 1, the mean rate without adding flanking edges is Re: how to pick the least-bad adjustment for the recombination rate? I agree with what @petrelharp outlined above. We have to pick between either the recombination rate or Ne agreeing with the time scale (but can't quite have both) and this is the closest to getting us to recombination rates being in "generations" which seems to be the least-bad option, in my opinion. Also, I hadn't taken a close look at #857, but I will now! |
|
Ah, sorry - I wasn't talking about the question about including the flanking zeros or not. What you're doing there seems good. |
|
Hey @srgibbon17, I'm working on addressing #1859. TL;DR I introduced a bug in the cM column when lifting over the coordinates to TAIR10. I noticed my mean rates are slightly different now that I'm using {'1': np.float64(3.550101369046314e-08),
'2': np.float64(3.7066078164228056e-08),
'3': np.float64(3.638229478106127e-08),
'4': np.float64(4.2697122435425374e-08),
'5': np.float64(3.66415642846763e-08)}Since our mean rate values matched before this, I'm guessing we both did the same thing previously and did read_hapmap without the |
|
Note that you are effectively editing his PR in, @silastittes, since your PR will go onto his. So I suggest just making that change here, after @srgibbon17 confirms. |
|
Hey @silastittes, that seems good to me. From Peter's comment, it seems like you'll make the necessary changes and edit the existing rates, so I don't need to edit my PR directly? But, if not, it's very simple to edit the mean rates myself |
|
Ok cool! Just made the updates. The tests passed locally cause I have the new map cached, but makes sense the local recombination rate will fail until we update the genetic map on AWS. |
| # # extend to end of chromosome and remake map | ||
| # if map_ex.right[-1] < chrom_length: | ||
| # new_position = np.append(map_ex.position, chrom_length) | ||
| # new_rate = np.append(rate, 0) |
There was a problem hiding this comment.
I think this should append nan not 0, since we don't know the recombination rate here (we don't think it's zero).
| # new_rate = np.append(rate, 0) | |
| # new_rate = np.append(rate, np.nan) |
But, won't that be equivalent to just computing the mean rate directly without extending the maps? (If so that's good, let's do that?)
|
say, @silastittes @srgibbon17 can we get this finished off? |
|
Yes! I was waiting for the updated genetic map to be on AWS, but would you prefer we move forward without that updated first? |
|
Ah sorry should we be bugging @andrewkern then? |
|
There is, however, the unresolved comment above? Namely: I don't think that we should be appending 0 to the maps, I think we should treat the rate on those regions as unknown, in which case I think we don't need to extend the maps at all to compute the mean rates. And, I haven't seen any discussion or justification of the choice you did make there? |
|
I agree that leaving it as unknown makes more sense. I guess the rational is the current set up here matches what we do in genetic maps. Would it be at all confusing if users run |
|
I think the current behavior there is a bug: #1832. |
|
Sounds good. |
|
Sorry for all the nitpicking! Thanks! |
|
Mind squashing these commits down, and then we'll merge? |
AraTha recombination rate (take 2) pre-commit style updates comments and long_description edits Update stdpopsim/catalog/AraTha/species.py Co-authored-by: Peter Ralph <petrel.harp@gmail.com> edit comments and add/edit citations Update stdpopsim/catalog/AraTha/genetic_maps.py Co-authored-by: Peter Ralph <petrel.harp@gmail.com> Update stdpopsim/catalog/AraTha/species.py Co-authored-by: Peter Ralph <petrel.harp@gmail.com> added citation to Platt et al for recombination rate update changelog update changelog AraTha first pass QC for 1855 local test and revisions for mean test corrected to genetic map mean rate chromosome length adjustment selfing correction bug flake8 black and flake8 edits for AraTha genetic map bug fix local test updated removed chromomsome extension for mean rate



QC for #1855
Tests fail!
Here's how I built the dictionary within the test. I downloaded the salome2012 genetic map
Made a list of the chromosome files and calculated the mean rate for each
Wrong map file? Missing step? Something else?
Do we need an additional test besides the chromosome means?