Repository: popsim-consortium/stdpopsim (language: python)

Issue:
QC for Caenorhabditis elegans
**PR for new species:** #904

Further details in issue #834.

If you volunteer to QC this species, please use the checklist below.
While this list is intended to be comprehensive, it may not be exhaustive.
Where relevant, the QC reviewer should identify that parameter values match
those given in the linked citation(s).

The QC reviewer should start a pull request that fills out the test stubs
with independently obtained values. The reviewer may look at the python code
for rationale provided in comments, but should ignore the actual code
as much as possible - comments in the code should give enough information
that it's obvious how to get the correct value from the provided references.
(In particular, we shouldn't copy-paste the value from the code into the test!)

- [ ] Recombination rate.
  - This might be genome-wide, or per-chromosome. Both are fine.
  - Check there's a comment describing where it came from, and/or how calculated.
  - From a publication? Check the value(s) match the publication.
  - Calculated somehow? Average over a recombination map? Redo the calculation.
- [ ] Mutation rate.
  - This might be genome-wide, or per-chromosome. Both are fine.
  - Check there's a comment describing where it came from, and/or how calculated.
  - From a publication? Check the value(s) match the publication.
  - Calculated somehow? Redo the calculation.
- [ ] Recombination map (if present).
  - Does it match the assembly? Liftover is fine, if clearly stated.
  - Is the description/long_description a good summary of how the map was created?
- [ ] Population size.
- [ ] Generation time.

For each citation, check:
- Doi link.
- Is publication a preprint? Is there a peer-reviewed publication instead?
- Is the year correct.
- Is the author correct (spelling, accents/ligatures/etc.)

Citations are required for:
- [ ] Genome reference assembly.
- [ ] Mutation rate.
- [ ] Recombination rate.
- [ ] Recombination map(s) (if relevant).
- [ ] Population size.
- [ ] Generation time.

The final PR should:
- [ ] fill out the test stubs
- [ ] delete the `pytest.mark.skip` lines that make the tests not run
- [ ] make sure they pass, talking to the original author to figure out discrepancies

C. elegans genetic map
A follow-up PR to #904.

Note: The genetic map has only 3 columns.


Discussion / hints:
Hi @grahamgower, I'm going to do the QC for CaeEle
Hey @attrna ! I have some questions about the citations and values for C. elegans.

- **Generation time**. The number I gleaned from the [citation](https://doi.org/10.7554/eLife.05849") was 3 days: "newly hatched individuals pass through four larval stages (L1, L2, L3, L4) and reach the adult stage after 3 days." It looks like you were going with the 150gens/year = 2.43 days? I thought that value was reported less confidently in the paper. 
- **Recombination rate**. I read in the [citation](https://doi.org/10.1371/journal.pgen.1000419) that each chromosome has a genetic length of 50cM. Shouldn't we use that to calculate rec per chromosome? 
- **Mutation rate**. It looks like the mtDNA mutation rate isn't being implemented. We should specify the mtDNA estimate from the second citation right?
- **Population size**. There were two citations([#1](https://doi.org/10.1016/j.cub.2005.06.022) and [#2](https://doi.org/10.1534/genetics.105.048207)). I was leaning towards the second one, because it's representative of all of Europe? Either way we should choose one of the two, right?


Hi @attrna. Are you still interested in having this work merged into stdpopsim? It looks like there's a bit of polishing needed, but it looks like you've already done most of the hard work.
Hi @grahamgower! Yes, I do. I'm so sorry it took so long!
Please could you help me with that? Should I do rebasing or something else too as it's been a while since I did the PR?


Thank you!
Anastasia


Yes, you'll definitely need to rebase, because there have been a lot of changes in the past 12 months. It would be great if you could also address the other small comments I made last year. There are lots of instructions in the [dev docs](https://popsim-consortium.github.io/stdpopsim-docs/latest/development.html) for rebasing, but this can be a pain even for experienced git users. So if you get tangled just let me know and I can have a go at it.
Hi @grahamgower! I made all the small changes you suggested above, but I'm not sure what commits to squash and where exactly to rebase (how to pick the base for rebasing). I guess I should not commit the changes before rebasing, is it right? Please could you help me with that? 
Hi @attrna. You will need to commit your changes before rebasing. Check the dev docs for additional info: https://popsim-consortium.github.io/stdpopsim-docs/latest/development.html#rebasing

* After committing, fetch the latest changes from the stdpopsim upstream repository. `git fetch upstream`, where "upstream" is the name of the upstream repository. This will probably be called "upstream" but you can double check by looking for the popsim-consortium/stdpopsim repository in the output of `git remote -v`.
* Make a copy of your working branch (`git branch celegansgenmap-copy`), so if you get stuck you can reset your working branch to match the copy (`git reset --hard celegansgenmap-copy`). You can also push your working branch to your github fork as a backup (it should also appear here in your pull request as usual).
* rebase onto the "upstream/main" branch (where "upstream" is once again the name for the upstream repository). `git rebase -i upstream/main`.
* Don't worry about squashing anything right now, just see if the rebase works. You can rebase again later and squash things.
* If you get conflicts during rebasing, you'll need to manually edit the conflicting files to fix things. Make sure to carefully read the messages that git gives you, and calling `git status` between edits can be helpful.
* If all goes well, you can try rebasing again and squashing commits. "One idea, one commit" is a good rule of thumb. The main idea for your pull request is to add the genetic map for CaeEle, so probably one commit is fine.
Hi @grahamgower! Thank you so much for your help! I did it, probably I messed something up because it doesn't let me push. Should I do git pull before all that?

```
anastasia@Anastasias-MacBook-Pro:CaeEle$ git push
To https://github.com/attrna/stdpopsim
 ! [rejected]        celegansgenmap -> celegansgenmap (non-fast-forward)
error: failed to push some refs to 'https://github.com/attrna/stdpopsim'
hint: Updates were rejected because the tip of your current branch is behind
hint: its remote counterpart. Integrate the remote changes (e.g.
hint: 'git pull ...') before pushing again.
hint: See the 'Note about fast-forwards' in 'git push --help' for details  
```
It seems that every time I'm doing "git rebase -i upstream/main" it finds some conflicts and I don't get how to resolve them, because it's basically the old version on the file.
# [Codecov](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/910?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=popsim-consortium) Report
> Merging [#910](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/910?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=popsim-consortium) (2764d95) into [main](https://codecov.io/gh/popsim-consortium/stdpopsim/commit/446062384d6c482277d21b032fa7f748713389c8?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=popsim-consortium) (4460623) will **increase** coverage by `0.00%`.
> The diff coverage is `100.00%`.

> :exclamation: Current head 2764d95 differs from pull request most recent head 1eb7b82. Consider uploading reports for the commit 1eb7b82 to get more accurate results

```diff
@@           Coverage Diff           @@
##             main     #910   +/-   ##
=======================================
  Coverage   99.56%   99.56%           
=======================================
  Files         102      103    +1     
  Lines        3243     3248    +5     
  Branches      415      415           
=======================================
+ Hits         3229     3234    +5     
  Misses          6        6           
  Partials        8        8           
```


| [Impacted Files](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/910?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=popsim-consortium) | Coverage Δ | |
|---|---|---|
| [stdpopsim/catalog/CaeEle/\_\_init\_\_.py](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/910/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=popsim-consortium#diff-c3RkcG9wc2ltL2NhdGFsb2cvQ2FlRWxlL19faW5pdF9fLnB5) | `100.00% <100.00%> (ø)` | |
| [stdpopsim/catalog/CaeEle/genetic\_maps.py](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/910/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=popsim-consortium#diff-c3RkcG9wc2ltL2NhdGFsb2cvQ2FlRWxlL2dlbmV0aWNfbWFwcy5weQ==) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/910?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=popsim-consortium).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=popsim-consortium)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/910?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=popsim-consortium). Last update [4460623...1eb7b82](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/910?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=popsim-consortium). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=popsim-consortium).

Thanks @attrna. I rebased for you and *force pushed* the branch back to the github repository. I also had issues with the rebase because of the same contents being present multiple times, so I'm not sure what happened there. After rebasing, one usually has to `git push -f` to force push the branch.
Thank you so much, @grahamgower!!  

Fix patch:
diff --git a/stdpopsim/catalog/CaeEle/genetic_maps.py b/stdpopsim/catalog/CaeEle/genetic_maps.py
index b7557869b..d9b97e36a 100644
--- a/stdpopsim/catalog/CaeEle/genetic_maps.py
+++ b/stdpopsim/catalog/CaeEle/genetic_maps.py
@@ -18,10 +18,8 @@
         first recombination breakpoint observed in the RIAILs.
         The genetic map corresponds to the assembly ce11 (GCA_000002985.3).
         """,
-    url=(
-        "https://stdpopsim.s3.us-west-2.amazonaws.com/genetic_maps/"
-        "CaeEle/rockman2009_maps.tgz"
-    ),
+    url="https://stdpopsim.s3.us-west-2.amazonaws.com/genetic_maps/CaeEle/"
+    "rockman2009_maps.tgz",
     sha256="ef0efa0aec3aa8fcd9800830d83d80bb2525788cf7b28105b7292072fedad8fb",
     file_pattern="genetic_map/C.elegans.Rockman.Kruglyak.2009.{id}.hapmap.txt",
     citations=[
diff --git a/stdpopsim/catalog/CaeEle/species.py b/stdpopsim/catalog/CaeEle/species.py
index 1dfbc7ddc..679941fb1 100644
--- a/stdpopsim/catalog/CaeEle/species.py
+++ b/stdpopsim/catalog/CaeEle/species.py
@@ -55,18 +55,44 @@
     reasons={stdpopsim.CiteReason.POP_SIZE},
 )
 
-# Mean chromosomal rates, calculated from the Rockman and Kruglyak 2009 map.
-# !!! The recombination rate is multiplied by 0.001 because the outcrosssing rate is 0.1%
+_SivasundarAndHey2003 = stdpopsim.Citation(
+    doi="https://doi.org/10.1093/genetics/163.1.147",
+    year=2003,
+    author="Sivasundar & Hey",
+    reasons={stdpopsim.CiteReason.POP_SIZE},
+)
+
+_SivasundarAndHey2005 = stdpopsim.Citation(
+    doi="https://doi.org/10.1016/j.cub.2005.08.034",
+    year=2005,
+    author="Sivasundar & Hey",
+    reasons={stdpopsim.CiteReason.POP_SIZE},
+)
+
+# Mean chromosomal rates, calculated from the
+# Rockman and Kruglyak 2009 map. !!! The recombination
+# rate is multiplied by 0.001 because the outcrosssing
+# rate is 0.1%
 _recombination_rate_data = {
-    "I": 3.330038e-11,
-    "II": 3.999342e-11,
-    "III": 4.484974e-11,
-    "IV": 2.417689e-11,
-    "V": 2.722476e-11,
-    "X": 3.447911e-11,
+    "I": 3.1216265402124167e-11,
+    "II": 3.529290802315087e-11,
+    "III": 3.906598767640363e-11,
+    "IV": 2.712098077556377e-11,
+    "V": 2.4705737572511805e-11,
+    "X": 2.9472374817864404e-11,
     "MtDNA": 0,
 }
 
+genomic_dna_mu = 1.84e-9
+_mutation_rate_data = {
+    "I": genomic_dna_mu,
+    "II": genomic_dna_mu,
+    "III": genomic_dna_mu,
+    "IV": genomic_dna_mu,
+    "V": genomic_dna_mu,
+    "X": genomic_dna_mu,
+    "MtDNA": 1.05e-7,
+}
 
 _chromosomes = []
 for name, data in genome_data.data["chromosomes"].items():
@@ -75,14 +101,15 @@
             id=name,
             length=data["length"],
             synonyms=data["synonyms"],
-            mutation_rate=1.84e-9,  # _Konrad et al. de-nove mutation rate,
-            # it's not uniform and it's much better to use a mutation map.
+            mutation_rate=_mutation_rate_data[name],
+            # _Konrad et al. de-nove mutation rate,
+            # it's not uniform and it's much better to use
+            # a mutation map.
             # mutation_rate=_mutation_rate_data[name],
             recombination_rate=_recombination_rate_data[name],
         )
     )
 
-
 _genome = stdpopsim.Genome(
     chromosomes=_chromosomes,
     assembly_name=genome_data.data["assembly_name"],
@@ -97,13 +124,16 @@
 
 _species = stdpopsim.Species(
     id="CaeEle",
-    ensembl_id="",
+    ensembl_id="caenorhabditis_elegans",
     name="Caenorhabditis elegans",
     common_name="C. elegans",
     genome=_genome,
-    generation_time=0.01,  # the generation time in the lab ~150
-    # generation per year (0.00666), it should be less in the wild
-    population_size=10000,
+    generation_time=0.01,  # One of the estimates reported
+    # by the paper is ~150 generations per year (0.00666).
+    # We settled on 0.01, because generations are expected
+    # to be fewer in nature due to dauer larvae stage.
+    population_size=10000,  # Lots of estimates available,
+    # they are all over the place. Settling on 10,000.
     #    selfing_rate,=0.999,
     citations=[
         _FrezalAndFelix2015.because(stdpopsim.CiteReason.GEN_TIME),


Test patch:
diff --git a/tests/test_CaeEle.py b/tests/test_CaeEle.py
index e6dfa97fd..356fd84da 100644
--- a/tests/test_CaeEle.py
+++ b/tests/test_CaeEle.py
@@ -8,8 +8,8 @@ class TestSpeciesData(test_species.SpeciesTestBase):
 
     species = stdpopsim.get_species("CaeEle")
 
-    #   def test_ensembl_id(self):
-    #        assert self.species.ensembl_id == ""
+    def test_ensembl_id(self):
+        assert self.species.ensembl_id == "caenorhabditis_elegans"
 
     def test_name(self):
         assert self.species.name == "Caenorhabditis elegans"
@@ -21,31 +21,86 @@ def test_common_name(self):
     # independently referring to the citations provided in the
     # species definition, filling in the appropriate values
     # and deleting the pytest "skip" annotations.
-    @pytest.mark.skip("Population size QC not done yet")
     def test_qc_population_size(self):
-        assert self.species.population_size == -1
+        assert self.species.population_size == 10000
 
-    @pytest.mark.skip("Generation time QC not done yet")
     def test_qc_generation_time(self):
-        assert self.species.generation_time == -1
+        assert self.species.generation_time == 0.01
 
 
 class TestGenomeData(test_species.GenomeTestBase):
 
     genome = stdpopsim.get_species("CaeEle").genome
+    mu = 1.84e-9
+
+    # # downloading genetic map, for calculating total cM per chromosome
+    # comm = "wget http://sesame.uoregon.edu/~ateterina/rockman2009_maps.tgz"
+    # os.system(comm)
+    # try:
+    #     comm = "tar -zxvf rockman2009_maps.tgz"
+    #     os.system(comm)
+    # except:
+    #     pass
+    # my_recs = {}
+    # genome = stdpopsim.get_species("CaeEle").genome
+    # chromlist = ["I", "II", "III", "IV", "V", "X"]
+    # for chrom in chromlist:
+    #     bps, rates, cM = [], [], 0
+    #     with open(
+    #         "genetic_map/C.elegans.Rockman.Kruglyak.2009." + chrom + ".hapmap.txt"
+    #     ) as infile:
+    #         infile.readline()  # header
+    #         newline = infile.readline().strip().split()
+    #         bp = float(newline[1])
+    #         rate = float(newline[2])
+    #         bps.append(bp)
+    #         rates.append(rate)
+    #         for line in infile:
+    #             newline = line.strip().split()
+    #             bp = float(newline[1])
+    #             rate = float(newline[2])
+    #             bps.append(bp)
+    #             rates.append(rate)
+    #             previous_index = bps.index(bp) - 1
+    #             length = bp - bps[previous_index]
+    #             cM += rates[previous_index] * length
+    #     total_rate = cM / genome.get_chromosome(chrom).length
+    #     total_rate /= 100  # cM to rate
+    #     total_rate /= 1000000  # Mb to bp
+    #     total_rate /= 1000  # 0.1% outcrossing
+    #     my_recs[chrom] = total_rate
+    #     print(my_recs)
+
+    my_recs = {
+        "I": 3.1216265402124167e-11,
+        "II": 3.529290802315087e-11,
+        "III": 3.906598767640363e-11,
+        "IV": 2.712098077556377e-11,
+        "V": 2.4705737572511805e-11,
+        "X": 2.9472374817864404e-11,
+        "MtDNA": 0,
+    }
 
-    @pytest.mark.skip("Recombination rate QC not done yet")
     @pytest.mark.parametrize(
         ["name", "rate"],
-        {"I": -1, "II": -1, "III": -1, "IV": -1, "V": -1, "X": -1, "MtDNA": -1}.items(),
+        my_recs.items(),
     )
     def test_recombination_rate(self, name, rate):
-        assert pytest.approx(rate, self.genome.get_chromosome(name).recombination_rate)
+        assert rate == pytest.approx(
+            self.genome.get_chromosome(name).recombination_rate
+        )
 
-    @pytest.mark.skip("Mutation rate QC not done yet")
     @pytest.mark.parametrize(
         ["name", "rate"],
-        {"I": -1, "II": -1, "III": -1, "IV": -1, "V": -1, "X": -1, "MtDNA": -1}.items(),
+        {
+            "I": mu,
+            "II": mu,
+            "III": mu,
+            "IV": mu,
+            "V": mu,
+            "X": mu,
+            "MtDNA": 1.05e-7,
+        }.items(),
     )
     def test_mutation_rate(self, name, rate):
-        assert pytest.approx(rate, self.genome.get_chromosome(name).mutation_rate)
+        assert rate == pytest.approx(self.genome.get_chromosome(name).mutation_rate)