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

Issue:
include gene conversion feature of msprime/slim into stdpopsim 
When we created stdpopsim, msprime did not have the gene conversion feature.
As gc is now possible in msprime as well as slim and gene conversion is a major factor of microbial evolution, it would be worth to include this as a possible parameter to the species models.

For E.Coli there might already be a commented part in the code that relates to the gene conversion of E.Coli.

Gene conversion can basically be handled like recombination. So if the parameter is not set, then gc is zero.
Otherwise, we have to specify a genome wide average for the gene conversion initiation rate and the mean tract length.

Add gene conversion parameters to E. coli. species model
When gene conversion can be handled as discussed in #846 we should add the corresponding parameters from Lapierre etl al. into the existing E.coli. species.


depends on #846


Discussion / hints:
I propose we add a ``gene_conversion_rate`` and ``gene_conversion_tract_length`` parameters to the Chromosome class, which default to zero. This can then be picked up in the contig and then passed to the simulation engine. 


Fix patch:
diff --git a/stdpopsim/catalog/EscCol/species.py b/stdpopsim/catalog/EscCol/species.py
index 430a81257..571b0a65c 100644
--- a/stdpopsim/catalog/EscCol/species.py
+++ b/stdpopsim/catalog/EscCol/species.py
@@ -33,6 +33,8 @@
             # from synonymous substitutions over 40,000 generations.
             mutation_rate=8.9e-11,
             recombination_rate=0.0,
+            gene_conversion_rate=8.9e-11,
+            gene_conversion_length=345,
         )
     )
 
diff --git a/stdpopsim/cli.py b/stdpopsim/cli.py
index 008fe42cb..165fda711 100644
--- a/stdpopsim/cli.py
+++ b/stdpopsim/cli.py
@@ -248,6 +248,14 @@ def get_species_help(species_id):
     species_text += (
         f"Recombination rate: {species.genome.mean_recombination_rate:.4g}\n"
     )
+    species_text += (
+        f"Gene conversion rate: {species.genome.mean_gene_conversion_rate:.4g}\n"
+    )
+    species_text += (
+        f"Gene conversion tract length: "
+        f"{species.genome.range_gene_conversion_lengths}\n"
+    )
+
     return species_text
 
 
diff --git a/stdpopsim/engines.py b/stdpopsim/engines.py
index a79f7a80a..50ca991e4 100644
--- a/stdpopsim/engines.py
+++ b/stdpopsim/engines.py
@@ -254,6 +254,8 @@ def simulate(
         ts = msprime.sim_ancestry(
             samples=samples,
             recombination_rate=contig.recombination_map,
+            gene_conversion_rate=contig.gene_conversion_rate,
+            gene_conversion_tract_length=contig.gene_conversion_length,
             demography=demographic_model.model,
             ploidy=2,
             random_seed=seeds[0],
diff --git a/stdpopsim/genomes.py b/stdpopsim/genomes.py
index 22c619b0b..1879adab8 100644
--- a/stdpopsim/genomes.py
+++ b/stdpopsim/genomes.py
@@ -23,7 +23,7 @@ class Genome:
     :vartype chromosomes: list
     :ivar citations: A list of :class:`.Citation` objects
         providing the source for the genome assembly,
-        mutation rate and recombination rate estimates.
+        mutation rate, recombination rate, and gene conversion rate estimates.
     :vartype citations: list
     :ivar length: The total length of the genome.
     :vartype length: int
@@ -93,6 +93,34 @@ def get_chromosome(self, id):
                 return chrom
         raise ValueError("Chromosome not found")
 
+    @property
+    def mean_gene_conversion_rate(self):
+        """
+        The length-weighted mean gene conversion rate across all chromosomes.
+        """
+        length = self.length
+        mean_gene_conversion_rate = 0
+        for chrom in self.chromosomes:
+            normalized_weight = chrom.length / length
+            if chrom.gene_conversion_rate is None:
+                rate = 0.0
+            else:
+                rate = chrom.gene_conversion_rate
+            cont = rate * normalized_weight
+            mean_gene_conversion_rate += cont
+        return mean_gene_conversion_rate
+
+    @property
+    def range_gene_conversion_lengths(self):
+        """
+        The range of gene conversion tract lengths across chromosomes.
+        """
+        prelim_lengths = [chrom.gene_conversion_length for chrom in self.chromosomes]
+        gc_tract_lengths = [0 if gctl is None else gctl for gctl in prelim_lengths]
+        min_gctl = min(gc_tract_lengths)
+        max_gctl = max(gc_tract_lengths)
+        return f"{min_gctl}" if min_gctl == max_gctl else f"{min_gctl} - {max_gctl}"
+
     @property
     def mean_recombination_rate(self):
         """
@@ -131,6 +159,10 @@ class Chromosome:
         chromosome.
     :ivar float recombination_rate: The recombination rate used when simulating
         this chromosome (if not using a genetic map).
+    :ivar float gene_conversion_rate: The gene conversion rate used when
+        simulating this chromosome.
+    :ivar float gene_conversion_length: The mean tract length of gene
+        conversion events used when simulating this chromosome.
     :ivar synonyms: List of synonyms that may be used when requesting this
         chromosome by ID, e.g. from the command line interface.
     :vartype synonyms: list of str
@@ -139,6 +171,8 @@ class Chromosome:
     id = attr.ib(type=str, kw_only=True)
     length = attr.ib(kw_only=True)
     recombination_rate = attr.ib(type=float, kw_only=True)
+    gene_conversion_rate = attr.ib(default=None, type=float, kw_only=True)
+    gene_conversion_length = attr.ib(default=None, type=float, kw_only=True)
     mutation_rate = attr.ib(type=float, kw_only=True)
     synonyms = attr.ib(factory=list, kw_only=True)
 
@@ -148,8 +182,8 @@ class Contig:
     """
     Class representing a contiguous region of genome that is to be simulated.
     This contains the information about mutation rates, distributions of
-    fitness effects (DFEs), and recombination rates that are needed to simulate
-    this region.
+    fitness effects (DFEs), gene conversion rates and recombination rates
+    that are needed to simulate this region.
 
     Information about targets of selection are contained in the ``dfe_list``
     and ``interval_list``. These must be of the same length, and the k-th DFE
@@ -157,6 +191,11 @@ class Contig:
 
     :ivar mutation_rate: The rate of mutation per base per generation.
     :vartype mutation_rate: float
+    :ivar gene_conversion_rate: The rate of gene conversion initiation per base per
+        generation.
+    :vartype gene_conversion_rate: float
+    :ivar gene_conversion_length: The mean tract length of gene conversions.
+    :vartype gene_conversion_length: float
     :ivar recombination_map: The recombination map for the region. See the
         :class:`msprime.RateMap` for details.
     :vartype recombination_map: msprime.RateMap
@@ -192,6 +231,8 @@ class Contig:
 
     recombination_map = attr.ib()
     mutation_rate = attr.ib(type=float)
+    gene_conversion_rate = attr.ib(default=None, type=float, kw_only=True)
+    gene_conversion_length = attr.ib(default=None, type=float, kw_only=True)
     genetic_map = attr.ib(default=None)
     inclusion_mask = attr.ib(default=None)
     exclusion_mask = attr.ib(default=None)
@@ -205,9 +246,21 @@ def __attrs_post_init__(self):
         )
 
     @staticmethod
-    def basic_contig(*, length, mutation_rate=0, recombination_rate=0):
+    def basic_contig(
+        *,
+        length,
+        mutation_rate=0,
+        recombination_rate=0,
+        gene_conversion_rate=None,
+        gene_conversion_length=None,
+    ):
         recomb_map = msprime.RateMap.uniform(length, recombination_rate)
-        return Contig(recombination_map=recomb_map, mutation_rate=mutation_rate)
+        return Contig(
+            recombination_map=recomb_map,
+            mutation_rate=mutation_rate,
+            gene_conversion_rate=gene_conversion_rate,
+            gene_conversion_length=gene_conversion_length,
+        )
 
     @staticmethod
     def species_contig(
@@ -218,6 +271,9 @@ def species_contig(
         length_multiplier=1,
         length=None,
         mutation_rate=None,
+        use_species_gene_conversion=False,
+        gene_conversion_rate=None,
+        gene_conversion_length=None,
         inclusion_mask=None,
         exclusion_mask=None,
     ):
@@ -234,6 +290,24 @@ def species_contig(
                     "https://github.com/popsim-consortium/stdpopsim/issues/406"
                 )
             )
+        if use_species_gene_conversion is True:
+            if gene_conversion_rate is not None or gene_conversion_length is not None:
+                raise ValueError(
+                    "Cannot use species gene conversion rates "
+                    "and set gene conversion rate or length"
+                )
+        if gene_conversion_rate is not None:
+            if gene_conversion_length is None:
+                raise ValueError(
+                    "Cannot set gene conversion rate without setting "
+                    "gene conversion length"
+                )
+        if gene_conversion_length is not None:
+            if gene_conversion_rate is None:
+                raise ValueError(
+                    "Cannot set gene conversion length without setting "
+                    "gene conversion rate"
+                )
         if chromosome is None:
             if genetic_map is not None:
                 raise ValueError("Cannot use genetic map with generic contig")
@@ -246,18 +320,32 @@ def species_contig(
             L_tot = 0
             r_tot = 0
             u_tot = 0
+            gc_tot = 0
+            gcl_tot = 0
             for chrom_data in species.genome.chromosomes:
                 if chrom_data.id.lower() not in non_autosomal_lower:
                     L_tot += chrom_data.length
                     r_tot += chrom_data.length * chrom_data.recombination_rate
                     u_tot += chrom_data.length * chrom_data.mutation_rate
+                    if chrom_data.gene_conversion_rate is None:
+                        local_gc_rate = 0
+                    else:
+                        local_gc_rate = chrom_data.gene_conversion_rate
+                        gcl_tot += chrom_data.length * chrom_data.gene_conversion_length
+                    gc_tot += chrom_data.length * local_gc_rate
             if mutation_rate is None:
                 mutation_rate = u_tot / L_tot
             r = r_tot / L_tot
+            if use_species_gene_conversion is True:
+                if gc_tot > 0 and gcl_tot > 0:
+                    gene_conversion_rate = gc_tot / L_tot
+                    gene_conversion_length = gcl_tot / L_tot
             contig = Contig.basic_contig(
                 length=length,
                 mutation_rate=mutation_rate,
                 recombination_rate=r,
+                gene_conversion_rate=gene_conversion_rate,
+                gene_conversion_length=gene_conversion_length,
             )
         else:
             if length is not None:
@@ -302,10 +390,16 @@ def species_contig(
             if mutation_rate is None:
                 mutation_rate = chrom.mutation_rate
 
+            if use_species_gene_conversion is True:
+                gene_conversion_rate = chrom.gene_conversion_rate
+                gene_conversion_length = chrom.gene_conversion_length
+
             contig = stdpopsim.Contig(
                 recombination_map=recomb_map,
                 mutation_rate=mutation_rate,
                 genetic_map=gm,
+                gene_conversion_rate=gene_conversion_rate,
+                gene_conversion_length=gene_conversion_length,
                 inclusion_mask=inclusion_intervals,
                 exclusion_mask=exclusion_intervals,
             )
@@ -439,12 +533,15 @@ def __str__(self):
         gmap = "None" if self.genetic_map is None else self.genetic_map.id
         s = (
             "Contig(length={:.2G}, recombination_rate={:.2G}, "
-            "mutation_rate={:.2G}, genetic_map={}, dfe_list={}, "
+            "mutation_rate={:.2G}, gene_conversion_rate={}, "
+            "gene_conversion_length{}, genetic_map={}, dfe_list={}, "
             "interval_list={})"
         ).format(
             self.length,
             self.recombination_map.mean_rate,
             self.mutation_rate,
+            self.gene_conversion_rate,
+            self.gene_conversion_length,
             gmap,
             self.dfe_list,
             self.interval_list,
diff --git a/stdpopsim/species.py b/stdpopsim/species.py
index 95e6f96c0..716b0ac1d 100644
--- a/stdpopsim/species.py
+++ b/stdpopsim/species.py
@@ -170,6 +170,9 @@ def get_contig(
         length_multiplier=1,
         length=None,
         mutation_rate=None,
+        use_species_gene_conversion=False,
+        gene_conversion_rate=None,
+        gene_conversion_length=None,
         inclusion_mask=None,
         exclusion_mask=None,
     ):
@@ -195,6 +198,19 @@ def get_contig(
             ``genetic_map`` argument.
         :param float mutation_rate: The per-base mutation rate. If none is given,
             the mutation rate defaults to the rate specified by species chromosomes.
+        :param bool use_species_gene_conversion: If set to True the parameters for gene
+            conversion of the species chromosome are used if available. For "generic"
+            contigs the gene conversion rate is given by the (mean) rate specified by
+            species chromosomes and the gene conversion length defaults to the mean
+            length specified by species chromosomes.
+        :param float gene_conversion_rate: The per-base gene conversion rate. Can only
+            be set together with gene_conversion_length and if
+            use_species_gene_conversion is False; If none is given, and
+            use_species_gene_conversion is True the gene conversion rate specified by
+            species chromosomes is used. If none and use_species_gene_conversion is
+            False gene conversion is disabled.
+        :param float gene_conversion_length: The mean gene conversion length.
+            May only be provided if gene_conversion_rate is also specified.
         :param inclusion_mask: If specified, simulated genomes are subset to only
             inlude regions given by the mask. The mask can be specified by the
             path and file name of a bed file or as a list or array of intervals
@@ -217,6 +233,9 @@ def get_contig(
             length_multiplier=length_multiplier,
             length=length,
             mutation_rate=mutation_rate,
+            use_species_gene_conversion=use_species_gene_conversion,
+            gene_conversion_rate=gene_conversion_rate,
+            gene_conversion_length=gene_conversion_length,
             inclusion_mask=inclusion_mask,
             exclusion_mask=exclusion_mask,
         )
diff --git a/validation.py b/validation.py
index 24d8a64fb..d1bf7a03b 100755
--- a/validation.py
+++ b/validation.py
@@ -42,6 +42,8 @@ def irradiate(contig, x=20):
     """
     return stdpopsim.Contig(
         recombination_map=contig.recombination_map,
+        gene_conversion_rate=contig.gene_conversion_rate,
+        gene_conversion_length=contig.gene_conversion_length,
         mutation_rate=x * contig.mutation_rate,
         genetic_map=contig.genetic_map,
     )


Test patch:
diff --git a/tests/test_EscCol.py b/tests/test_EscCol.py
index 4657a8fcb..3b8526826 100644
--- a/tests/test_EscCol.py
+++ b/tests/test_EscCol.py
@@ -26,3 +26,12 @@ class TestGenome(test_species.GenomeTestBase):
 
     def test_basic_attributes(self):
         assert len(self.genome.chromosomes) == 1
+
+    def test_mutation_rate(self):
+        assert self.genome.get_chromosome("Chromosome").mutation_rate == 8.9e-11
+
+    def test_gene_conversion_rate(self):
+        assert self.genome.get_chromosome("Chromosome").gene_conversion_rate == 8.9e-11
+
+    def test_gene_conversion_length(self):
+        assert self.genome.get_chromosome("Chromosome").gene_conversion_length == 345
diff --git a/tests/test_cli.py b/tests/test_cli.py
index 8c94dbdf1..812474f2d 100644
--- a/tests/test_cli.py
+++ b/tests/test_cli.py
@@ -254,7 +254,7 @@ def test_durvusula_2017_msmc(self):
         self.verify(cmd, num_samples=7)
 
     def test_lapierre_constant(self):
-        cmd = "EscCol -L 100 2"
+        cmd = "EscCol -L 1000 2"
         self.verify(cmd, num_samples=2)
 
 
@@ -985,10 +985,10 @@ def verify_noQC_warning(self, cmd):
                 capture_output(stdpopsim.cli.stdpopsim_main, cmd.split())
 
     def test_noQC_warning(self):
-        self.verify_noQC_warning("EscCol -d FakeModel -D 10 -L 10")
+        self.verify_noQC_warning("EscCol -d FakeModel -D 10 -L 1000")
 
     def test_noQC_warning_quiet(self):
-        self.verify_noQC_warning("-q EscCol -d FakeModel -D 10 -L 10")
+        self.verify_noQC_warning("-q EscCol -d FakeModel -D 10 -L 1000")
 
     def verify_noQC_citations_not_written(self, cmd, caplog):
         # Non-QCed models shouldn't be used in publications, so citations
@@ -1011,14 +1011,14 @@ def verify_noQC_citations_not_written(self, cmd, caplog):
     @pytest.mark.usefixtures("caplog")
     def test_noQC_citations_not_written(self, caplog):
         self.verify_noQC_citations_not_written(
-            "EscCol -d FakeModel -D 10 -L 10", caplog
+            "EscCol -d FakeModel -D 10 -L 1000", caplog
         )
 
     @pytest.mark.filterwarnings("ignore::stdpopsim.QCMissingWarning")
     @pytest.mark.usefixtures("caplog")
     def test_noQC_citations_not_written_verbose(self, caplog):
         self.verify_noQC_citations_not_written(
-            "-vv EscCol -d FakeModel -D 10 -L 10", caplog
+            "-vv EscCol -d FakeModel -D 10 -L 1000", caplog
         )
 
 
diff --git a/tests/test_genomes.py b/tests/test_genomes.py
index 7d441a2cc..862f8aae3 100644
--- a/tests/test_genomes.py
+++ b/tests/test_genomes.py
@@ -3,6 +3,7 @@
 """
 import numpy as np
 import pytest
+import copy
 import stdpopsim
 
 
@@ -25,6 +26,11 @@ def verify_dfe_breakpoints(self, contig):
                 k = np.searchsorted(breaks, left)
                 assert dfe_labels[k] == j
 
+    def test_default_gc(self):
+        contig = stdpopsim.Contig.basic_contig(length=42)
+        assert contig.gene_conversion_rate is None
+        assert contig.gene_conversion_length is None
+
     def test_default_dfe(self):
         contig = stdpopsim.Contig.basic_contig(length=100)
         assert len(contig.dfe_list) == 1
@@ -149,3 +155,41 @@ def test_is_neutral(self):
                 contig.add_dfe(intervals=np.array([[30, 40]]), DFE=dfes[1])
                 # exponential with mean zero doesn't count as neutral!
                 assert contig.is_neutral is (neutral and dist == "f")
+
+
+class TestGeneConversion(object):
+    species = stdpopsim.get_species("EscCol")
+
+    def test_modifiying_gene_conversion(self):
+        genome = self.species.get_contig(
+            chromosome="Chromosome", gene_conversion_rate=0.1, gene_conversion_length=1
+        )
+        assert genome.gene_conversion_rate == 0.1
+        assert genome.gene_conversion_length == 1
+
+    def test_use_species_gene_conversion(self):
+        genome = self.species.get_contig(
+            chromosome="Chromosome", use_species_gene_conversion=True
+        )
+        assert genome.gene_conversion_rate == 8.9e-11
+        assert genome.gene_conversion_length == 345
+
+    def test_unnamed_contig_with_gc(self):
+        genome = self.species.get_contig(
+            length=42, gene_conversion_rate=0.1, gene_conversion_length=1
+        )
+        assert genome.gene_conversion_rate == 0.1
+        assert genome.gene_conversion_length == 1
+
+    def test_use_species_gene_conversion_unnamed_contig(self):
+        genome = self.species.get_contig(length=42, use_species_gene_conversion=True)
+        assert genome.gene_conversion_rate == 8.9e-11
+        assert genome.gene_conversion_length == 345
+
+    def test_use_species_gene_conversion_unnamed_contig_undefined_gc(self):
+        mod_species = copy.deepcopy(self.species)
+        mod_species.genome.get_chromosome("Chromosome").gene_conversion_rate = None
+        mod_species.genome.get_chromosome("Chromosome").gene_conversion_length = None
+        genome = mod_species.get_contig(length=42, use_species_gene_conversion=True)
+        assert genome.gene_conversion_rate is None
+        assert genome.gene_conversion_length is None
diff --git a/tests/test_species.py b/tests/test_species.py
index c15d02b96..202457cc9 100644
--- a/tests/test_species.py
+++ b/tests/test_species.py
@@ -151,6 +151,15 @@ def test_generation_time_defined(self):
     def test_population_size_defined(self):
         assert self.species.population_size > 0
 
+    def test_default_gc(self):
+        for chrom in self.species.genome.chromosomes:
+            contig = self.species.get_contig(chrom.id)
+            assert contig.gene_conversion_rate is None
+            assert contig.gene_conversion_length is None
+            contig = self.species.get_contig(chrom.id, use_species_gene_conversion=True)
+            assert chrom.gene_conversion_rate == contig.gene_conversion_rate
+            assert chrom.gene_conversion_length == contig.gene_conversion_length
+
 
 class GenomeTestBase:
     """
@@ -253,38 +262,66 @@ def test_genetic_map(self):
         assert isinstance(contig.recombination_map, msprime.RateMap)
 
     def test_contig_options(self):
-        with pytest.raises(ValueError):
+        with pytest.raises(ValueError, match="Cannot use genetic map"):
             # cannot use genetic map with generic contig
             self.species.get_contig(genetic_map="ABC")
-        with pytest.raises(ValueError):
+        with pytest.raises(ValueError, match="Cannot use length multiplier"):
             # cannot use length multiplier with generic contig
             self.species.get_contig(length_multiplier=0.1)
-        with pytest.raises(ValueError):
+        with pytest.raises(ValueError, match="Must specify sequence length"):
             # must specify length with generic contig or give chromosome name
             self.species.get_contig()
-        with pytest.raises(ValueError):
+        with pytest.raises(ValueError, match="Cannot specify sequence length"):
             # cannot specify length with named chromosome
             self.species.get_contig("chr1", length=1e6)
-        with pytest.raises(ValueError):
-            # cannot specify mask for generic contig
+        with pytest.raises(ValueError, match="Cannot use mask"):
+            # cannot specify inclusion mask for generic contig
             self.species.get_contig(length=1e4, inclusion_mask=[(0, 100)])
-        with pytest.raises(ValueError):
-            # connot specify mask for generic contig
+        with pytest.raises(ValueError, match="Cannot use mask"):
+            # cannot specify exclusion mask for generic contig
             self.species.get_contig(length=1e4, exclusion_mask=[(0, 100)])
-        with pytest.raises(ValueError):
-            # cannot use length multiplier with mask
+        with pytest.raises(ValueError, match="Cannot use length multiplier"):
+            # cannot use length multiplier with inclusion mask
             self.species.get_contig(
                 "chr22", inclusion_mask=[(0, 100)], length_multiplier=0.1
             )
-        with pytest.raises(ValueError):
-            # cannot use length multiplier with mask
+        with pytest.raises(ValueError, match="Cannot use length multiplier"):
+            # cannot use length multiplier with exclusion mask
             self.species.get_contig(
                 "chr22", exclusion_mask=[(0, 100)], length_multiplier=0.1
             )
+        with pytest.raises(
+            ValueError, match="Cannot use species gene conversion rates"
+        ):
+            # cannot specify use_species_gene_conversion and custom gc rate
+            self.species.get_contig(
+                "chr1", use_species_gene_conversion=True, gene_conversion_rate=0.1
+            )
+        with pytest.raises(ValueError, match="without setting gene conversion length"):
+            # cannot specify custom gene conversion rate without gene conversion length
+            self.species.get_contig("chr1", gene_conversion_rate=0.1)
+        with pytest.raises(ValueError, match="without setting gene conversion rate"):
+            # cannot specify custom gene conversion length without gene conversion rate
+            self.species.get_contig("chr1", gene_conversion_length=1)
+
+    def test_use_species_gene_conversion(self):
+        contig = self.species.get_contig("chr22", use_species_gene_conversion=True)
+        if self.species.genome.get_chromosome("chr22").gene_conversion_rate is None:
+            assert contig.gene_conversion_rate is None
+            assert contig.gene_conversion_length is None
+        else:
+            assert (
+                contig.gene_conversion_rate
+                == self.species.genome.get_chromosome("chr22").gene_conversion_rate
+            )
+            assert (
+                contig.gene_conversion_length
+                == self.species.genome.get_chromosome("chr22").gene_conversion_length
+            )
 
     def test_generic_contig(self):
         L = 1e6
-        contig = self.species.get_contig(length=L)
+        contig = self.species.get_contig(length=L, use_species_gene_conversion=True)
         assert contig.recombination_map.sequence_length == L
 
         chrom_ids = np.arange(1, 23).astype("str")
@@ -299,6 +336,22 @@ def test_generic_contig(self):
             for c in self.species.genome.chromosomes
             if c.id in chrom_ids
         ]
+        gcs = [
+            c.gene_conversion_rate if c.gene_conversion_rate is not None else 0
+            for c in self.species.genome.chromosomes
+            if c.id in chrom_ids
+        ]
+        gcls = [
+            c.gene_conversion_length if c.gene_conversion_length is not None else 0
+            for c in self.species.genome.chromosomes
+            if c.id in chrom_ids
+        ]
 
         assert contig.mutation_rate == np.average(us, weights=Ls)
         assert contig.recombination_map.mean_rate == np.average(rs, weights=Ls)
+        if np.average(gcs, weights=Ls) == 0 or np.average(gcls, weights=Ls) == 0:
+            assert contig.gene_conversion_rate is None
+            assert contig.gene_conversion_length is None
+        else:
+            assert contig.gene_conversion_rate == np.average(gcs, weights=Ls)
+            assert contig.gene_conversion_length == np.average(gcls, weights=Ls)