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

Issue:
Support for gene conversion in SLiM engine
With #1106 the msprime engine supports gene conversion, and the rate/mean tract length parameters are stored as attributes of `stdpopsim.Contig`. It'd be nice to also allow gene conversion with the SLiM engine. This would require:

1. Having `slim_makescript` add a line with `initializeGeneConversion(...);` to the `initialize()` block of the SLiM script.
2. Reconciling differences in parameterization between the SLiM and msprime GC models. If I'm understanding the documentation correctly, msprime independently samples the locations of simple crossovers and GC events, using separately provided rates (a constant rate in the case of GC). Whereas SLiM draws breakpoints according to rates in a supplied recombination map, then randomly assigns each event to be a simple crossover or GC according to a supplied probability. So rates of GC events can vary across the sequence in SLiM, but not in msprime (currently). 
3. Support for heteroduplex mismatch repair -- there are two additional parameters in SLiM that control the proportion of GC tracts that are complex and the nucleotide bias resulting from mismatch repair within these tracts. These don't have an influence on ancestry but do influence stored mutations (?)

Relevant docs for [msprime](https://tskit.dev/msprime/docs/latest/ancestry.html#gene-conversion) and [SLiM](http://benhaller.com/slim/SLiM_Manual.pdf#page=32).


Discussion / hints:
Do you have input on how to set this up, @fbaumdicker?
Also TODO would be to add a bit to `slim_makescript` to ask SLiM what the GC rate is and record this in ts metadata when debugging is on, so that the tests can verify that GC is being done as requested; we can copy the way that we check recomb maps and DFEs.
With regards to SLiM's ability to model bias in nucleotide composition due to heteroduplex mismatch repair: currently, neutral mutations are added to SLiM-produced tree sequences in `stdpopsim._SLiMEngine._recap_and_rescale` using msprime. So, I think it'd be challenging to support mismatch repair bias because we'd have to do something complicated to keep track of mutations in complex gene conversion tracts?

For example, if there are nucleotide differences between the copy/non-copy strands within a complex gene conversion tract, then mismatches that are "repaired" to match the non-copy strand will be recorded as new mutations by SLiM. This requires knowing what mutations are present on the copy/non-copy strands during a particular complex gene conversion event. I don't think it's possible to know this, if neutral mutations are added post-SLiM-simulation on a possibly recapitated tree sequence.

So maybe we shouldn't worry about supporting nucleotide bias, and should instead focus on getting something as close as possible to the gene conversion model in msprime?
I agree - we can think about GC bias later.
OK, notes here:
- msprime's gc set-up is [described here](https://tskit.dev/pyslim/docs/latest/python_api.html#pyslim.generate_nucleotides)
- SLiM's gc set-up is described in section 1.5.6 of the manual

We have two parameters added in #1106, which are [passed straight to msprime.sim_ancestry](https://github.com/popsim-consortium/stdpopsim/blob/ec4fedd9859dc23b7b12013d69172d851a5d43d2/stdpopsim/engines.py#L257):
- `gene_conversion_rate`
- `gene_conversion_length`

Ok, so in msprime:
- lengths are Geometric with mean `gene_conversion_length`
- at a (Poisson) rate `gene_conversion_rate` (that must be constant)

Now, in SLiM, we specify:
- a recombation RateMap; let's say the rate at `x` is `r(x)`;
- lengths are Geometric *in each direction* (hence negative binomial) with a total mean of both directions of `meanLength`
- a proportion `nonCrossoverFraction` of DSBs are resolved as non-crossovers (so crossover rate is `r * (1 - nonCrossoverFraction)`)
- a proportion `simpleConversionFraction` of non-crossovers that resolve as "simple" (ie not involving heteroduplex repair)

So, we want `gene_conversion_length = meanLength`, and `simpleConversionFraction=1.0`. And, with msprime's (simple) recombination rate at `x` equal to `recombination_rate(x)`,
```
  recombination_rate(x) = r(x) ( 1- noncrossoverFraction )
  gene_conversion_rate = r(x) * noncrossoverFraction
```
Then we'd like to solve for `r(x)` and `noncrossoverFraction`.

However, these two models are fundamentally incompatible with non-constant rate maps, since:
- msprime only simulates a constant rate of GCs
- SLiM only simulates a rate of GCs that is a constant proportion of recombinations

Well, I'm glad we're noticing this before the release. The most obvious thing to do would be to set `noncrossoverFraction` in a given simulation so that the expected number of GCs is equal to `gene_conversion_rate * sequence_length`. But, maybe there's a better thing to do? Or, a better way to parameterize this?
> The most obvious thing to do would be to set noncrossoverFraction in a given simulation so that the expected number of GCs is equal to gene_conversion_rate * sequence_length. But, maybe there's a better thing to do? Or, a better way to parameterize this?

Another option I see is to change the way we store GC rates, effectively to match SLiM. This would make sense, as SLiM's set-up is closer to the biology (as I understand it, in diploid, sexual organisms at least). So, let's say that the crossover rate is `1e-8` and the crossover fraction is `0.1`; then we'd need to simulate in SLiM a recombination rate of `1e-7` and a noncrossover fraction of 0.9; in msprime we'd simulate with a recombination rate of `1e-8` and a gc rate of `0.9e-7`. This *almost* works, except in the case (eg of EscCol) where we don't want *any* crossovers; since the formula is `total_rate = crossover_rate / crossover_fraction`, if `crossover_rate = 0` (ie we set recombination rate to 0) then this doesn't work.

Note also that we already don't match SLiM because with GC turned on, "recombination rate" in SLiM means "rate of crossovers + noncrossovers", while in stdpopsim "recombination rate" means "rate of crossovers", and I don't think we want to change *that*, since if we did then just turning on GC would dramatically change what "recombination rate" means, counter-intuitively. (The situation is less counter-intuitive in SLiM because of the API, but in stdpopsim that one would affect the other is not so obvious.)

So, in summary, I see three ways forward:
1. match msprime's parameterization (what we currently do)
2. match SLiM's parameterization (I think we don't want to, because of the last paragraph above)
3. do a mix-and-match (doesn't work, as above)

... and, I've talked myself back to where we started out at. Thinking this through: what if someone wants to say that DSBs happen at rate proportional to a given recomb map, but then 90% of them resolve as GC events? If we require separate RateMaps for recomb and GC, then they'd be passing in two identical-up-to-scaling RateMaps. And, come to think of it, it would be more *realistic* for e.g. humans to specify just a `noncrossover_fraction` and pass that straight into SLiM instead of doing through the contortions to get SLiM to simulate from msprime's flat GC rate model (which we can't do anyhow).

Okay, so now here's

4. allow species to have *either* `gene_conversion_rate` *or* `gene_conversion_fraction`, since the first makes sense for EscCol and the second for HomSap. We'd have to figure out what to do in the cross-cases, but we can do that.

This seems right to me, but it makes things more complicated, which I don't like.

Unfortunately, I don't currently plan to include non-constant gene conversion rates in msprime, and I believe Ben recently mentioned that he also doesn't plan to allow the use of two independent recombination maps in slim.
Both would resolve the incompatibility here.

It might be worth looking at the solution @jeanrjc and @bhaller came up with when simulating bacterial recombination (= gene conversion) with slim in this paper: https://doi.org/10.24072/pcjournal.72
If we could use something similar, we would even get the same tract length distribution in msprime and slim.
At least for bacteria, i.e. if the recombination rate is 0 but the gc rate > 0, this might be the best solution.
However, it might also be possible to similarly include constant rate gene conversions in addition to a variable recombination map in eukaryotes.

If that doesn't solve the issue:
I like the simplicity of simply defining a `gene_conversion_fraction`, since this is how many people think about gc in eukaryotes. This would work well in the case of constant rate maps for slim and msprime.
But for non-constant recombination maps, it would not be possible to simulate them in this way in msprime. The only thing we could do is to set the constant gene conversion rate in msprime to a value such that the mean number of gene conversions is the same, with all the difficulties you describe.

The more I think about it, the more I conclude that the `gene_conversion_rate` is a species property, while the `gene_conversion_fraction` is more a recombination map property, since the gene conversion fraction estimates depend on the corresponding recombination estimates. 

So maybe we should distinguish between constant and non-constant recombination maps only when gc is enabled?

Constant rates --> use msprime as currently implemented, and do the same in slim with your formula above.
non-constant map --> use the constant gc rates in msprime and the fractions in slim, users need to specify how to (re)calculate the rates. We should avoid the scenario where turning on gc changes the crossover recombination rates by default.
> It might be worth looking at the solution @jeanrjc and @bhaller came up with when simulating bacterial recombination (= gene conversion) with slim in this paper: https://doi.org/10.24072/pcjournal.72 If we could use something similar, we would even get the same tract length distribution in msprime and slim. At least for bacteria, i.e. if the recombination rate is 0 but the gc rate > 0, this might be the best solution. However, it might also be possible to similarly include constant rate gene conversions in addition to a variable recombination map in eukaryotes.

Possible, yes.  Quite unfortunate to not be able to use the built-in gene conversion facility in SLiM, though, since it's really pretty sophisticated – even allowing for complex gene conversion tracts and heteroduplex mismatch repairs.  Abandoning that in favor of a scripted approach... well, it doesn't seem optimal.  But I get that you're in a bind here.

I'm not sure I understand all the issues at play enough to comment more than that.
I'm not worried about getting exactly the same model in msprime and SLiM - mostly, we'd like to get the same model in both because we have to make fewer decisions. However, we shouldn't fall back to a less realistic model just to get consistency across engines. (For instance, we're not restricting ourselves to neutral sims just because msprime won't do much selection.) And, it seems to me that SLiM is the more realistic one (?), so I'm not enthusiastic about doing something more complicated with a less realistic outcome.

> The more I think about it, the more I conclude that the gene_conversion_rate is a species property, while the gene_conversion_fraction is more a recombination map property, since the gene conversion fraction estimates depend on the corresponding recombination estimates.

I'm not sure what you mean? Are there recombination maps that have estimated non-constant gene conversion fractions? I thought we had at best a genome wide guess at gene conversion fraction per species?

Ah, I've thought of another alternative: remove `gene_conversion_rate`, and special-case bacterial: so, have a flag like `bacterial_recombination` that would be False for things that do crossing over, and these would have `gene_conversion_fraction`, while it would be True for things doing bacterial recombination, and for these "recombination rate" would be implemented as gene conversion rate (with no crossing over). So, EscCol would have `recombination_rate` equal to what's now in `gene_conversion_rate`, and `bacterial_recombination=True`.
> Ah, I've thought of another alternative: remove `gene_conversion_rate`, and special-case bacterial: so, have a flag like `bacterial_recombination` that would be False for things that do crossing over, and these would have `gene_conversion_fraction`, while it would be True for things doing bacterial recombination, and for these "recombination rate" would be implemented as gene conversion rate (with no crossing over). So, EscCol would have `recombination_rate` equal to what's now in `gene_conversion_rate`, and `bacterial_recombination=True`.

This sounds like an excellent proposal to me.  Biologically, horizontal gene transfer in bacteria is a totally different process than crossing over in diploids, and it makes sense to treat them separately.  If the gene conversion support in stdpopsim eventually goes further (with, for example, heteroduplex mismatch repair and GC-biased gene conversion), treating these two cases separately will make it much easier to expand the design.
I completely agree. Using the gene conversion parameters to mimic bacterial recombination should not limit the gene conversion features in diploids. So using a `bacterial_recombination` flag is the most natural thing to do and will at the same time avoid confusion.
We would still need the currently used `gene_conversion_tract_length` parameter in bacteria as well as diploids, which we might want to rename and split into two parameters? E.g. `bacterial_recombination_tract_length` and `gene_conversion_tract_length`.

> I'm not sure what you mean? Are there recombination maps that have estimated non-constant gene conversion fractions? I thought we had at best a genome wide guess at gene conversion fraction per species?

While there are some attempts to estimate the local gene conversion fractions along chromosomes, e.g. here: https://doi.org/10.1038/ejhg.2013.30 , that was not what I meant.
I am concerned that if we have different recombination maps in the catalog for the same species, we might need to also use different values for the `gene_conversion_fraction`, depending on the recombination map and the estimation method for the gene conversion parameters.
The same applies to the "gene_conversion_tract_length". This might also differ between different recombination maps.
Perhaps it would be best to have both, a default value for "gene_conversion_fraction", and a "gene_conversion_fraction" for each recombination map. This is basically what we currently do with the recombination rate itself anyway.



> We would still need the currently used gene_conversion_tract_length parameter in bacteria as well as diploids, which we might want to rename and split into two parameters? E.g. bacterial_recombination_tract_length and gene_conversion_tract_length.

I'd vote for leaving them named the same thing, since no species can have *both* properties, and they're implemented just the same. It'll be less confusing then adding an additional parameter.

>  Perhaps it would be best to have both, a default value for "gene_conversion_fraction", and a "gene_conversion_fraction" for each recombination map. This is basically what we currently do with the recombination rate itself anyway.

Is your idea to make things ready for the future in which recombination maps also estimate GC parameters/maps? Currently, are there any recombination maps that also have estimates of GC things (like, described in the same publication or something?)? We can always add the property in the future if it turns out to be important - but, is there any situation now where it'd be important to have?
Synthesizing the above discussion, to make things concrete (but with @fbaumdicker's suggestions still pending), here's a proposal:
- add a `bacterial_recombination` flag and a `gene_conversion_fraction` parameter (must be >= 0 and strictly < 1)
- if the flag is True then "recombination rate" means effectively GC rate, and the `gene_conversion_tract_length` must be defined; and `gene_conversion_fraction` must be not defined (ie None)
- furthermore if the flag is True then simulations will automatically turn on "gene conversion" (ie, bacterial recombination)
- if the flag is False then things work as currently implemented: by default there is no GC, and only if gene conversion is requested in the call to `simulate( )` is it turned on
- if gene conversion is turned on, then, in SLiM, the `gene_conversion_fraction` is passed straight to SLiM and the recombination rate map is multiplied by `1/(1-gene_conversion_fraction)` (so that the rate of crossing over is unchanged)
- if gene conversion is turned on in *msprime* then we work out what would be the average rate across the genome of gene conversion under SLiM's scheme (which would be `(total_map_length / sequence_length) * (gene_conversion_fraction / (1 - gene_conversion_fraction))`, I think), and add gene conversion at this rate
> > We would still need the currently used gene_conversion_tract_length parameter in bacteria as well as diploids, which we might want to rename and split into two parameters? E.g. bacterial_recombination_tract_length and gene_conversion_tract_length.
> 
> I'd vote for leaving them named the same thing, since no species can have _both_ properties, and they're implemented just the same. It'll be less confusing then adding an additional parameter.

I thought they were going to be implemented quite differently?  Bacterial recombination would be implemented using `addRecombinant()` as in section 16.13 in the SLiM manual, gene conversion would be implemented using SLiM's built-in gene conversion facilities.  I think of them as being completely different processes; gene conversion is gene conversion, bacterial recombination is horizontal gene transfer.  I personally find it a bit confusing for the two to be conflated; they are different biological processes and should be treated as such.  My two cents.  :->
> I thought they were going to be implemented quite differently?

That is true for forward-in-time simulations, but bacterial recombination has often been mimicked by gene conversion in coalescent-based backward simulations. So it depends a little bit on the user's background whether using the `gene_conversion_rate` parameter for bacterial recombination will be confusing or not.
> Is your idea to make things ready for the future in which recombination maps also estimate GC parameters/maps? Currently, are there any recombination maps that also have estimates of GC things (like, described in the same publication or something?)? We can always add the property in the future if it turns out to be important - but, is there any situation now where it'd be important to have?

I am not really an expert regarding recombination maps but thought that the current estimates of gene conversion fractions would depend on the assumption of which recombination map is assumed or co-estimated.
There are three related publications I know of. 
One very recent paper by Derek Setter and Konrad Lohse, where they estimate chromosome-specific recombination and gene conversion rates. https://doi.org/10.1093/genetics/iyac100
And one publication from 2007 by S Myers and Gill McVean on a method to jointly estimate crossover and non-crossover rates. https://doi.org/10.1534/genetics.107.078907
There is also a publication estimating two separate maps for crossing over and gene conversion in Drosophila:
https://doi.org/10.1371/journal.pgen.1002905

However, I agree that this would be easy to add later, so maybe it is best to focus on the species parameters for now. Users would also always be able to use custom values for the gene conversion fractions.


> Synthesizing the above discussion, to make things concrete (but with @fbaumdicker's suggestions still pending), here's a proposal: ...

I think the proposal is fine but we should probably carefully decide which mechanisms to entangle.
Personally, I like the bacterial recombination flag idea and using the recombination rate for the bacterial recombination initiation rate. I am also slightly in favor of using as few parameters as possible which would mean using the `gene_conversion_tract_length` also for the length of bacterial homologous recombination tracts.
> I thought they were going to be implemented quite differently? Bacterial recombination would be implemented using addRecombinant() as in section 16.13 in the SLiM manual, gene conversion would be implemented using SLiM's built-in gene conversion facilities. I

Hm: the main difference I see between 16.13 and naively just doing `crossoverFraction=1.0` is handling ploidy correctly, which would be real nice but a whole different thing. There's also a different distribution of horizontally exchanged tract lengths (uniform instead of negative binomial).

So - I think that something like 16.13 probably *is* the right thing to do; however, so as to not snowball changes here, I'd like to just keep it simple, as I describe above - it's not perfect, but SLiM will not be less realistic than msprime, anyhow.
What do you mean by "handling ploidy correctly"?  Are you guys planning to simulate bacteria with a diploid model?
> There's also a different distribution of horizontally exchanged tract lengths (uniform instead of negative binomial).

But this would be super easy to change using a uniformly distributed first breakpoint and a geometrically distributed length to the next breakpoint. By changing only these two lines of code in the manual in 16.14

```C
// draw two distinct locations; redraw if we get a duplicate
do breaks = rdunif(2, max=L-1);
while (breaks[0] == breaks[1]);
```
The manual even suggests to do so:
> The details of the breakpoint generation here might need to be modified in a more realistic
model. Here we draw the start and end positions of the transfer region independently, but perhaps
it would be better to draw the start location randomly and then draw a transfer length from a
geometric distribution or some other distribution

So here it would be possible to make it exactly as msprime, but honestly, I think that the negative binomial will not make a big difference compared to the geometrically distributed tract lengths. So for the moment it might not make a huge difference to simply use `crossoverFraction=1.0`, but in the long run it might be clever to have separate setups for bacterial_recombination True and False.
> What do you mean by "handling ploidy correctly"? Are you guys planning to simulate bacteria with a diploid model?

That is, in fact, currently what happens. PRs welcome! Personally, dealing with sex chromosomes correctly (which we currently don't) is higher on my list of things to do.

>  in the long run it might be clever to have separate setups for bacterial_recombination True and False.

Agreed! And, this will set us up for that.


> > What do you mean by "handling ploidy correctly"? Are you guys planning to simulate bacteria with a diploid model?
> 
> That is, in fact, currently what happens. PRs welcome! Personally, dealing with sex chromosomes correctly (which we currently don't) is higher on my list of things to do.

Do you mean simulating just the X, or just the Y?  Or do you mean simulating the X and the Y and one or more autosomes?
> Do you mean simulating just the X, or just the Y? Or do you mean simulating the X and the Y and one or more autosomes?

Right now if you ask to simulate the X or the Y (or the Z or W) you just get an autosomal simulation (and a warning saying that it's not right).
> > Do you mean simulating just the X, or just the Y? Or do you mean simulating the X and the Y and one or more autosomes?
> 
> Right now if you ask to simulate the X or the Y (or the Z or W) you just get an autosomal simulation (and a warning saying that it's not right).

Gotcha.  I was wondering because if you want to support multiple-chromosome simulations it would make life easier if you waited for SLiM to add support for that first.  :->  But if it's just simulating the X or the Y by itself, that is already supported, as you know, so that would hopefully be easier to add in stdpopsim.  :->
I'm going to go ahead with this proposal, but let me know if you're not happy with it!
I'm fine with whatever choice you folks settle on, just giving input.  :->
I'm making (final) edits to the **adding species paper**, and realized after discussion with @fbaumdicker that we need to incorporate the changes you mention here in the manuscript text. Currently, the section on **Updates to stdpopsim and the species catalog** includes a short paragraph on gene conversion specifying how the parameterization is set (using Franz's original setup). I'm not sure that we want to dive into the detail of the different cases because this will distract readers from the main point, and might also change in the future. Currently, the only concrete example we give is for bacteria (*E. Coli*). So, I thought that we'd mostly focus on parameterization for bacteria and archea, which is fairly simple  -- with the `bacterial_recombination` flag set and the `gene_conversion_tract_length` set to the appropriate value. We could also mention that gene conversion may be used in other species together with crossover recombination, and parameterization there is more complex (without diving into the details). That leads to another question: do we have an example in the catalog for a species where we model gene conversion and crossover recombination?

Fix patch:
diff --git a/docs/_ext/speciescatalog.py b/docs/_ext/speciescatalog.py
index 05ddaa963..0ad32005d 100644
--- a/docs/_ext/speciescatalog.py
+++ b/docs/_ext/speciescatalog.py
@@ -572,6 +572,31 @@ def run(self):
         # genomes:
         genome_section = nodes.section(ids=[f"sec_catalog_{species.id}_genome"])
         genome_section += nodes.title(text="Genome")
+        if species.genome.bacterial_recombination is True:
+            genome_section += self.make_field_list(
+                [
+                    (
+                        "Bacterial recombination with tract length range",
+                        species.genome.range_gene_conversion_lengths,
+                        None,
+                    )
+                ]
+            )
+        elif species.genome.mean_gene_conversion_fraction > 0:
+            genome_section += self.make_field_list(
+                [
+                    (
+                        "Mean gene conversion fraction",
+                        species.genome.mean_gene_conversion_fraction,
+                        None,
+                    ),
+                    (
+                        "Range gene conversion lengths",
+                        species.genome.range_gene_conversion_lengths,
+                        None,
+                    ),
+                ]
+            )
         genome_section += self.chromosomes_table(species)
         genome_section += nodes.paragraph(
             text="Mutation and recombination rates "
diff --git a/docs/api.rst b/docs/api.rst
index 89899794f..f22ec0cf9 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -78,6 +78,73 @@ Demographic Models
 .. autoclass:: stdpopsim.Population()
     :members:
 
+
+.. _sec_api_gene_conversion:
+
+*******************************************
+Gene conversion and bacterial recombination
+*******************************************
+
+Some species have estimates of the process of gene conversion.
+However, by default gene conversion does *not* happen, because
+(unfortunately) it can make simulations take much longer,
+and result in much bigger files.
+This is because the rate of gene conversion is usually much higher than
+the rate of crossing-over, so enabling gene conversion effectively
+increases the rate of recombination, which can have a strong effect on runtimes.
+To enable gene conversion for a species that has estimates,
+pass the `use_species_gene_conversion=True` argument to
+:meth:`Engine.simulate`.
+
+Gene conversion attributes are visible through the `gene_conversion_fraction`
+and `gene_conversion_length` properties of the species' `genome`.
+The "gene conversion fraction" gives the (average) fraction of
+double-stranded breaks that resolve as gene conversion
+(a.k.a non-crossover) events,
+and the "gene conversion length" is the mean length of the gene conversion tracts.
+So, if the rate of crossing over (which is often referred to as
+"recombination rate") is `r` and the gene conversion fraction is `f`,
+then the *total* rate of double-stranded breaks is `r / (1-f)`,
+and the gene conversion rate is `r * f / (1-f)`.
+
+A consequence of this is that
+the `recombination_map` attribute of a :class:`Contig` is the rate of
+double-stranded break initiation, not the rate of crossovers. The terminology
+conflicts somewhat with the `recombination_rate` property of
+a :class:`Chromosome`, which specifies the rate of crossovers. So, creating
+a contig with `use_species_gene_conversion=True` will result in a Contig
+with larger "recombination rates" than otherwise, because these rates
+include gene conversion as well as crossing over:
+
+.. code-block:: python
+
+    species = stdpopsim.get_species("DroMel")
+    contig = species.get_contig("2L")
+    contig.recombination_map.mean_rate
+    # 2.40462600791e-08
+    contig = species.get_contig("2L", use_species_gene_conversion=True)
+    contig.recombination_map.mean_rate
+    # 1.414485887005882e-07
+
+
+Different engines implement gene conversion in different ways: at the time
+of writing, msprime only allows a *constant* rate of gene conversion,
+while SLiM only allows gene conversion to be a fixed fraction of double-stranded
+break events. So, in SLiM we use the recombination map for the local
+rate of double-stranded breaks along the genome, while in msprime
+we specify a constant rate of gene conversion so that the average total number
+is the same as we would get from SLiM.
+
+In principle, **bacterial recombination** is a completely different issue,
+since bacterial recombination is by horizontal transfer of genomic segments.
+However, bacterial recombination is at present implemented using the engine's
+gene conversion mechanisms. So, bacterial species have the `bacterial_recombination`
+flag set, and need `gene_conversion_length` to be defined. (It could also be called horizontal
+transfer or homologous recombination segment length, but that would add another
+option.) For such species, gene conversion does not happen, so `gene_conversion_fraction`
+must not be set, and the "recombination rate" is the rate of bacterial recombination.
+
+
 .. _sec_api_dfes:
 
 ******************
diff --git a/stdpopsim/catalog/DroMel/species.py b/stdpopsim/catalog/DroMel/species.py
index 0b5e876ff..ce7522da9 100644
--- a/stdpopsim/catalog/DroMel/species.py
+++ b/stdpopsim/catalog/DroMel/species.py
@@ -33,6 +33,7 @@
     author="Comeron et al",
     doi="https://doi.org/10.1371/journal.pgen.1002905",
     year=2012,
+    reasons={stdpopsim.CiteReason.REC_RATE, stdpopsim.CiteReason.GENE_CONVERSION},
 )
 
 # Mean chromosomal rates, calculated by taking the
@@ -42,7 +43,7 @@
 # Chromosome 4 isn't in this map, so the average of
 # 2L, 2R, 3L and 3R weighted by their respective
 # chromosome lengths was used instead.
-_recombination_rate_data = {
+_recombination_rate = {
     "2L": 2.40462600791e-08,
     "2R": 2.23458641776e-08,
     "3L": 1.79660308862e-08,
@@ -54,27 +55,52 @@
 }
 
 
-_chromosomes = []
-for name, data in genome_data.data["chromosomes"].items():
-    _chromosomes.append(
-        stdpopsim.Chromosome(
-            id=name,
-            length=data["length"],
-            synonyms=data["synonyms"],
-            mutation_rate=5.49e-9,  # _SchriderEtAl de novo mutation rate
-            recombination_rate=_recombination_rate_data[name],
-        )
-    )
-
-_genome = stdpopsim.Genome(
-    chromosomes=_chromosomes,
-    assembly_name=genome_data.data["assembly_name"],
-    assembly_accession=genome_data.data["assembly_accession"],
+# Comeron et al:
+# - GC avg track length = 518bp
+# - 83% of DSBs are resolved as Gene Conversions
+# - avg rate of GC is 1.25E-07 / bp / female meiosis
+# (not sure if the latter estimate agrees)
+
+_gene_conversion_fraction = {
+    "2L": 0.83,
+    "2R": 0.83,
+    "3L": 0.83,
+    "3R": 0.83,
+    "4": 0,
+    "X": 0.83,
+    "Y": 0,
+    "mitochondrion_genome": 0,
+}
+
+_gene_conversion_length = {c: 518 for c in genome_data.data["chromosomes"]}
+
+_mutation_rate = {c: 5.49e-9 for c in genome_data.data["chromosomes"]}
+
+# _chromosomes = []
+# for name, data in genome_data.data["chromosomes"].items():
+#     _chromosomes.append(
+#         stdpopsim.Chromosome(
+#             id=name,
+#             length=data["length"],
+#             synonyms=data["synonyms"],
+#             mutation_rate=5.49e-9,  # _SchriderEtAl de novo mutation rate
+#             recombination_rate=_recombination_rate_data[name],
+#             gene_conversion_fraction=_gene_conversion_fraction_data[name],
+#             gene_conversion_length=_gene_conversion_length,
+#         )
+#     )
+
+_genome = stdpopsim.Genome.from_data(
+    genome_data.data,
+    recombination_rate=_recombination_rate,
+    mutation_rate=_mutation_rate,
+    gene_conversion_fraction=_gene_conversion_fraction,
+    gene_conversion_length=_gene_conversion_length,
     citations=[
         _SchriderEtAl.because(stdpopsim.CiteReason.MUT_RATE),
         _DosSantosEtAl,
         _HoskinsEtAl,
-        _ComeronEtAl.because(stdpopsim.CiteReason.REC_RATE),
+        _ComeronEtAl,
     ],
 )
 stdpopsim.utils.append_common_synonyms(_genome)
diff --git a/stdpopsim/catalog/EscCol/species.py b/stdpopsim/catalog/EscCol/species.py
index 571b0a65c..ff7749b69 100644
--- a/stdpopsim/catalog/EscCol/species.py
+++ b/stdpopsim/catalog/EscCol/species.py
@@ -22,6 +22,12 @@
     doi="https://doi.org/10.1126/science.277.5331.1453",
 )
 
+_didelot_et_al = stdpopsim.Citation(
+    author="Didelot et al.",
+    year=2012,
+    doi="https://doi.org/10.1186/1471-2164-13-256",
+)
+
 _chromosomes = []
 for name, data in genome_data.data["chromosomes"].items():
     _chromosomes.append(
@@ -32,22 +38,22 @@
             # Wielgoss et al. (2011) calculated for strain REL606,
             # 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,
+            recombination_rate=8.9e-11,
+            gene_conversion_length=542,
         )
     )
 
-# mean_conversion_rate=8.9e-11 # not implemented yet!
-# mean_conversion_length=542 # not implemented yet!
-
 _genome = stdpopsim.Genome(
     chromosomes=_chromosomes,
+    bacterial_recombination=True,
     assembly_name=genome_data.data["assembly_name"],
     assembly_accession=genome_data.data["assembly_accession"],
     citations=[
-        _wielgoss_et_al.because(stdpopsim.CiteReason.MUT_RATE),
+        _wielgoss_et_al.because(
+            {stdpopsim.CiteReason.MUT_RATE, stdpopsim.CiteReason.GENE_CONVERSION}
+        ),
         _blattner_et_al.because(stdpopsim.CiteReason.ASSEMBLY),
+        _didelot_et_al.because(stdpopsim.CiteReason.GENE_CONVERSION),
     ],
 )
 stdpopsim.utils.append_common_synonyms(_genome)
diff --git a/stdpopsim/catalog/StrAga/species.py b/stdpopsim/catalog/StrAga/species.py
index 8e7217fdc..67b34d0e4 100644
--- a/stdpopsim/catalog/StrAga/species.py
+++ b/stdpopsim/catalog/StrAga/species.py
@@ -16,9 +16,6 @@
 # mu =  0.56/(1e6*365)
 _mutation_rate = {"1": 1.53e-9}
 
-# no cross-over recombination in bacteria
-_recombination_rate = {"1": 0}
-
 # gene conversion rate (HGT rate):
 # hard to estimate precisely. What is often estimated is the r/m or rho/theta ratio.
 # in Oliveira et. al, PNAS, 2016, they estimate those rates on coregenome
@@ -35,7 +32,7 @@
 # that is already not recombining a lot.
 # Overall, as most of the estimation might be conservatives one can use a relative
 # rate rho/theta = 10%.
-# _gene_conversion_rate = {"1": _mutation_rate/10}
+_recombination_rate = {c: _mutation_rate[c] / 10 for c in _mutation_rate}
 
 # Mean gene conversion  tract length
 # In Brochet et al., 2008, PNAS, they estimate the size of transferred DNA experimentally
@@ -46,7 +43,7 @@
 # discussed in more detail in response to a reviewer, that can be seen there :
 # https://evolbiol.peercommunityin.org/articles/rec?id=356
 #      > Revision round #3 > Author's reply (pdf).
-# _mean_gene_conversion_tract_length = 120000
+_gene_conversion_length = {"1": 120000}
 
 
 _genome = stdpopsim.Genome.from_data(
@@ -57,6 +54,8 @@
         _DaCunha_et_al.because(stdpopsim.CiteReason.ASSEMBLY),
         _DaCunha_et_al.because(stdpopsim.CiteReason.MUT_RATE),
     ],
+    bacterial_recombination=True,
+    gene_conversion_length=_gene_conversion_length,
 )
 stdpopsim.utils.append_common_synonyms(_genome)
 
diff --git a/stdpopsim/citations.py b/stdpopsim/citations.py
index 6d614761d..ef994df14 100644
--- a/stdpopsim/citations.py
+++ b/stdpopsim/citations.py
@@ -20,6 +20,7 @@ class CiteReason:
     GEN_TIME = "generation time"
     MUT_RATE = "mutation rate"
     REC_RATE = "recombination rate"
+    GENE_CONVERSION = "gene conversion parameters"
     ASSEMBLY = "genome assembly"
     ANNOTATION = "genome annotation"
     DFE = "distribution of fitness effects"
diff --git a/stdpopsim/cli.py b/stdpopsim/cli.py
index 89706c2ac..2c95dfe7d 100644
--- a/stdpopsim/cli.py
+++ b/stdpopsim/cli.py
@@ -248,13 +248,22 @@ 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"
-    )
+    if species.genome.bacterial_recombination is True:
+        species_text += (
+            f"Bacterial recombination with tract length: "
+            f"{species.genome.range_gene_conversion_lengths}\n"
+        )
+    # TODO: there is no way to enable GC on the CLI
+    # elif ((species.genome.mean_gene_conversion_fraction is not None) and
+    #     (species.genome.mean_gene_conversion_fraction > 0.0)):
+    #     species_text += (
+    #         f"Gene conversion fraction: "
+    #         f"{species.genome.mean_gene_conversion_fraction:.4g}\n"
+    #     )
+    #     species_text += (
+    #         f"Gene conversion tract length: "
+    #         f"{species.genome.range_gene_conversion_lengths}\n"
+    #     )
 
     return species_text
 
@@ -883,6 +892,11 @@ def write_simulation_summary(
     dry_run_text += f"{indent}Mean recombination rate: {mean_recomb_rate}\n"
     dry_run_text += f"{indent}Mean mutation rate: {mut_rate}\n"
     dry_run_text += f"{indent}Genetic map: {gmap}\n"
+    if contig.bacterial_recombination is True:
+        dry_run_text += (
+            f"{indent}Bacterial recombination "
+            f"with tract length {contig.gene_conversion_length}\n"
+        )
     if dfe is not None:
         dry_run_text += f"{indent}DFE: {dfe}\n"
         dry_run_text += f"{indent}DFE applied to: {dfe_interval}\n"
diff --git a/stdpopsim/engines.py b/stdpopsim/engines.py
index 50ca991e4..bb1e8a30e 100644
--- a/stdpopsim/engines.py
+++ b/stdpopsim/engines.py
@@ -251,10 +251,30 @@ def simulate(
         rng = np.random.default_rng(seed)
         seeds = rng.integers(1, 2**31 - 1, size=2)
 
+        # if bacterial_recombination=True then "recombination_rate"
+        # is implemented by msprime as gene conversion rate
+        # (which must be constant)
+        recombination_map = contig.recombination_map
+        gc_rate = None
+        gc_frac = contig.gene_conversion_fraction
+        if contig.bacterial_recombination:
+            gc_rate = contig.recombination_map.mean_rate
+            recombination_map = msprime.RateMap(
+                position=[0.0, contig.length], rate=[0.0]
+            )
+        elif (gc_frac is not None) and (gc_frac > 0.0):
+            # the recombination map is a map of double-stranded breaks,
+            # so we need to separate out GC from crossovers
+            gc_rate = contig.recombination_map.mean_rate * gc_frac / (1 - gc_frac)
+            recombination_map = msprime.RateMap(
+                position=recombination_map.position,
+                rate=recombination_map.rate * (1 - gc_frac),
+            )
+
         ts = msprime.sim_ancestry(
             samples=samples,
-            recombination_rate=contig.recombination_map,
-            gene_conversion_rate=contig.gene_conversion_rate,
+            recombination_rate=recombination_map,
+            gene_conversion_rate=gc_rate,
             gene_conversion_tract_length=contig.gene_conversion_length,
             demography=demographic_model.model,
             ploidy=2,
diff --git a/stdpopsim/genomes.py b/stdpopsim/genomes.py
index 4780d2089..c770dbaf0 100644
--- a/stdpopsim/genomes.py
+++ b/stdpopsim/genomes.py
@@ -21,9 +21,13 @@ class Genome:
 
     :ivar chromosomes: A list of :class:`.Chromosome` objects.
     :vartype chromosomes: list
+    :ivar bacterial_recombination: Whether recombination is via horizontal gene
+        transfer (if this is True) or via crossover and possibly gene
+        conversion (if this is False). Default: False.
+    :vartype bacterial_recombination: bool
     :ivar citations: A list of :class:`.Citation` objects
         providing the source for the genome assembly,
-        mutation rate, recombination rate, and gene conversion rate estimates.
+        mutation rate, recombination rate, and gene conversion estimates.
     :vartype citations: list
     :ivar length: The total length of the genome.
     :vartype length: int
@@ -34,10 +38,20 @@ class Genome:
     chromosomes = attr.ib(factory=list)
     assembly_name = attr.ib(default=None, kw_only=True)
     assembly_accession = attr.ib(default=None, kw_only=True)
+    bacterial_recombination = attr.ib(type=bool, default=False, kw_only=True)
     citations = attr.ib(factory=list, kw_only=True)
 
     @staticmethod
-    def from_data(genome_data, *, recombination_rate, mutation_rate, citations):
+    def from_data(
+        genome_data,
+        *,
+        recombination_rate,
+        mutation_rate,
+        citations,
+        bacterial_recombination=False,
+        gene_conversion_fraction=None,
+        gene_conversion_length=None,
+    ):
         """
         Construct a Genome object from the specified dictionary of
         genome information from Ensembl, recombination_rate and
@@ -48,6 +62,12 @@ def from_data(genome_data, *, recombination_rate, mutation_rate, citations):
         chr_names = set(genome_data["chromosomes"].keys())
         assert set(recombination_rate.keys()) == chr_names
         assert set(mutation_rate.keys()) == chr_names
+        if gene_conversion_fraction is None:
+            gene_conversion_fraction = {k: None for k in chr_names}
+        assert set(gene_conversion_fraction.keys()) == chr_names
+        if gene_conversion_length is None:
+            gene_conversion_length = {k: None for k in chr_names}
+        assert set(gene_conversion_length.keys()) == chr_names
         chromosomes = []
         for name, data in genome_data["chromosomes"].items():
             chromosomes.append(
@@ -57,12 +77,15 @@ def from_data(genome_data, *, recombination_rate, mutation_rate, citations):
                     synonyms=data["synonyms"],
                     mutation_rate=mutation_rate[name],
                     recombination_rate=recombination_rate[name],
+                    gene_conversion_fraction=gene_conversion_fraction[name],
+                    gene_conversion_length=gene_conversion_length[name],
                 )
             )
         return Genome(
             chromosomes=chromosomes,
             assembly_name=genome_data["assembly_name"],
             assembly_accession=genome_data["assembly_accession"],
+            bacterial_recombination=bacterial_recombination,
             citations=citations,
         )
 
@@ -94,21 +117,25 @@ def get_chromosome(self, id):
         raise ValueError("Chromosome not found")
 
     @property
-    def mean_gene_conversion_rate(self):
+    def mean_gene_conversion_fraction(self):
         """
-        The length-weighted mean gene conversion rate across all chromosomes.
+        The genetic length-weighted mean gene conversion fraction 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
+        if self.bacterial_recombination is True:
+            return 0.0
+        else:
+            total_length = 0
+            total_gc_fraction = 0
+            for chrom in self.chromosomes:
+                length = chrom.length * chrom.recombination_rate
+                total_length += length
+                if chrom.gene_conversion_fraction is None:
+                    frac = 0.0
+                else:
+                    frac = chrom.gene_conversion_fraction
+                total_gc_fraction += frac * length
+            return total_gc_fraction / total_length
 
     @property
     def range_gene_conversion_lengths(self):
@@ -125,6 +152,7 @@ def range_gene_conversion_lengths(self):
     def mean_recombination_rate(self):
         """
         The length-weighted mean recombination rate across all chromosomes.
+        (Note that this is the rate of crossovers, not all double-stranded breaks.)
         """
         length = self.length
         mean_recombination_rate = 0
@@ -151,15 +179,18 @@ def mean_mutation_rate(self):
 @attr.s
 class Chromosome:
     """
-    Class representing a single chromosome for a species.
+    Class representing a single chromosome for a species. Note that although
+    the recombination rate for a Chromosome is the rate of crossovers,
+    the recombination map for a Contig gives the rate of both crossovers
+    and gene conversion, so will differ if the gene conversion fraction is nonzero.
 
     :ivar str ~.id: The string identifier for the chromosome.
     :ivar int length: The length of the chromosome.
     :ivar float mutation_rate: The mutation rate used when simulating this
         chromosome.
-    :ivar float recombination_rate: The recombination rate used when simulating
+    :ivar float recombination_rate: The rate of crossovers used when simulating
         this chromosome (if not using a genetic map).
-    :ivar float gene_conversion_rate: The gene conversion rate used when
+    :ivar float gene_conversion_fraction: The gene conversion fraction used when
         simulating this chromosome.
     :ivar float gene_conversion_length: The mean tract length of gene
         conversion events used when simulating this chromosome.
@@ -171,7 +202,7 @@ 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_fraction = 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)
@@ -194,14 +225,27 @@ 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.
+    :ivar recombination_map: The recombination map for the region, that gives
+        the rates of double-stranded breaks. Note that if gene conversion
+        fraction is nonzero, then this map can have a larger rate than the
+        corresponding chromosome's recombination rate, since the chromosome's
+        rate describes only crossovers. A :class:`msprime.RateMap` object.
     :vartype recombination_map: msprime.RateMap
+    :ivar bacterial_recombination: Whether the model of recombination is
+        by horizontal gene transfer or not (default is False, i.e.,
+        recombination is by crossovers and possibly gene conversion).
+    :ivar gene_conversion_fraction: The fraction of recombinations that resolve
+        as gene conversion events. The recombination map gives the rates of
+        *both* crossovers and gene conversions, with a fraction of gene
+        conversions given by this value. Defaults to None, i.e., no gene
+        conversion. Must be None or between 0 and 1. Must be None if
+        bacterial_recombination is True.
+    :vartype gene_conversion_fraction: float
+    :ivar gene_conversion_length: The mean tract length of gene conversions (if
+        bacterial_recombination is False), or horizontally tranferred segments (if
+        it is True). Must be present, and above 1 if either gene_conversion_fraction
+        is given or bacterial_recombination is True.
+    :vartype gene_conversion_length: float
     :ivar mask_intervals: Intervals to keep in simulated tree sequence, as a list
         of (left_position, right_position), such that intervals are non-overlapping
         and in ascending order. Should have shape Nx2, where N is the number of
@@ -223,10 +267,11 @@ class Contig:
     :vartype original_coordinates: tuple
 
     .. note::
-        To run stdpopsim simulations with alternative, user-specified mutation
-        or recombination rates, a new contig can be created based on an existing
-        one. For instance, the following will create a ``new_contig`` that,
-        when simulated, will have double the mutation rate of the ``old_contig``:
+        To run stdpopsim simulations with alternative, user-specified mutation,
+        recombination, or gene conversion rates, a new contig can be created
+        based on an existing one. For instance, the following will create a
+        ``new_contig`` that, when simulated, will have double the mutation rate
+        of the ``old_contig``:
 
         .. code-block:: python
 
@@ -239,7 +284,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)
+    bacterial_recombination = attr.ib(default=False, type=bool, kw_only=True)
+    gene_conversion_fraction = 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)
@@ -263,14 +309,45 @@ def basic_contig(
         length,
         mutation_rate=0,
         recombination_rate=0,
-        gene_conversion_rate=None,
+        bacterial_recombination=False,
+        gene_conversion_fraction=None,
         gene_conversion_length=None,
     ):
+        if bacterial_recombination:
+            if gene_conversion_fraction is not None:
+                raise ValueError(
+                    "Cannot set gene conversion fraction for bacterial recombination."
+                )
+            if gene_conversion_length is None:
+                raise ValueError(
+                    "Must set gene conversion length for bacterial recombination."
+                )
+        else:
+            if gene_conversion_fraction is not None:
+                if gene_conversion_length is None:
+                    raise ValueError(
+                        "Cannot set gene conversion fraction without setting "
+                        "gene conversion length"
+                    )
+            if gene_conversion_length is not None:
+                if gene_conversion_fraction is None:
+                    raise ValueError(
+                        "Cannot set gene conversion length without setting "
+                        "gene conversion fraction"
+                    )
+        if (gene_conversion_fraction is not None) and (
+            gene_conversion_fraction < 0 or gene_conversion_fraction > 1
+        ):
+            raise ValueError("Gene conversion fraction must be between 0 and 1.")
+        if (gene_conversion_length is not None) and (gene_conversion_length < 1):
+            # values of even 0 are OK for SLiM but not msprime
+            raise ValueError("Gene conversion length must be greater than 1.")
         recomb_map = msprime.RateMap.uniform(length, recombination_rate)
         return Contig(
             recombination_map=recomb_map,
             mutation_rate=mutation_rate,
-            gene_conversion_rate=gene_conversion_rate,
+            bacterial_recombination=bacterial_recombination,
+            gene_conversion_fraction=gene_conversion_fraction,
             gene_conversion_length=gene_conversion_length,
         )
 
@@ -284,8 +361,6 @@ def species_contig(
         length=None,
         mutation_rate=None,
         use_species_gene_conversion=False,
-        gene_conversion_rate=None,
-        gene_conversion_length=None,
         inclusion_mask=None,
         exclusion_mask=None,
         left=None,
@@ -304,24 +379,15 @@ 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 (
+            species.genome.bacterial_recombination is True
+            and use_species_gene_conversion is True
+        ):
+            raise ValueError(
+                "Cannot use species gene conversion with bacterial recombination."
+            )
+        gene_conversion_fraction = None
+        gene_conversion_length = None
         if chromosome is None:
             if left is not None or right is not None:
                 raise ValueError(
@@ -338,34 +404,45 @@ def species_contig(
             L_tot = 0
             r_tot = 0
             u_tot = 0
-            gc_tot = 0
+            gcf_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 chrom_data.gene_conversion_fraction is not None:
+                        gcf_tot += (
+                            chrom_data.length
+                            * chrom_data.recombination_rate
+                            * chrom_data.gene_conversion_fraction
+                        )
+                    if chrom_data.gene_conversion_length is not None:
+                        gcl_tot += (
+                            chrom_data.length
+                            * chrom_data.recombination_rate
+                            * chrom_data.gene_conversion_length
+                        )
             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
+                if gcf_tot > 0:
+                    gene_conversion_fraction = gcf_tot / r_tot
+                    gene_conversion_length = gcl_tot / r_tot
+                    r /= 1 - gene_conversion_fraction
+            if species.genome.bacterial_recombination is True:
+                gene_conversion_length = gcl_tot / r_tot
             contig = Contig.basic_contig(
                 length=length,
                 mutation_rate=mutation_rate,
                 recombination_rate=r,
-                gene_conversion_rate=gene_conversion_rate,
+                bacterial_recombination=species.genome.bacterial_recombination,
+                gene_conversion_fraction=gene_conversion_fraction,
                 gene_conversion_length=gene_conversion_length,
             )
         else:
+            # Named contig:
             if length is not None:
                 raise ValueError("Cannot specify sequence length for named contig")
             if inclusion_mask is not None and exclusion_mask is not None:
@@ -447,15 +524,27 @@ def species_contig(
             if mutation_rate is None:
                 mutation_rate = chrom.mutation_rate
 
+            if species.genome.bacterial_recombination is True:
+                gene_conversion_length = chrom.gene_conversion_length
+
             if use_species_gene_conversion is True:
-                gene_conversion_rate = chrom.gene_conversion_rate
+                gene_conversion_fraction = chrom.gene_conversion_fraction
                 gene_conversion_length = chrom.gene_conversion_length
+                if (
+                    gene_conversion_fraction is not None
+                    and gene_conversion_fraction > 0
+                ):
+                    recomb_map = msprime.RateMap(
+                        position=recomb_map.position,
+                        rate=recomb_map.rate / (1 - gene_conversion_fraction),
+                    )
 
             contig = stdpopsim.Contig(
                 recombination_map=recomb_map,
                 mutation_rate=mutation_rate,
                 genetic_map=gm,
-                gene_conversion_rate=gene_conversion_rate,
+                bacterial_recombination=species.genome.bacterial_recombination,
+                gene_conversion_fraction=gene_conversion_fraction,
                 gene_conversion_length=gene_conversion_length,
                 inclusion_mask=inclusion_intervals,
                 exclusion_mask=exclusion_intervals,
@@ -694,14 +783,16 @@ 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}, gene_conversion_rate={}, "
-            "gene_conversion_length{}, genetic_map={}, origin={}, "
-            "dfe_list={}, interval_list={})"
+            "mutation_rate={:.2G}, bacterial_recombination={}, "
+            "gene_conversion_fraction={}, gene_conversion_length={}, "
+            "genetic_map={}, origin={}, dfe_list={}, "
+            "interval_list={})"
         ).format(
             self.length,
             self.recombination_map.mean_rate,
             self.mutation_rate,
-            self.gene_conversion_rate,
+            self.bacterial_recombination,
+            self.gene_conversion_fraction,
             self.gene_conversion_length,
             gmap,
             self.origin,
diff --git a/stdpopsim/species.py b/stdpopsim/species.py
index a1ce682da..1e12695e4 100644
--- a/stdpopsim/species.py
+++ b/stdpopsim/species.py
@@ -171,8 +171,6 @@ def get_contig(
         length=None,
         mutation_rate=None,
         use_species_gene_conversion=False,
-        gene_conversion_rate=None,
-        gene_conversion_length=None,
         inclusion_mask=None,
         exclusion_mask=None,
         left=None,
@@ -202,17 +200,8 @@ def get_contig(
             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.
+            contigs the gene conversion fraction and length are given by the mean
+            values across all chromosomes.
         :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
@@ -242,8 +231,6 @@ def get_contig(
             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,
             left=left,
diff --git a/validation.py b/validation.py
index d1bf7a03b..49f376664 100755
--- a/validation.py
+++ b/validation.py
@@ -42,7 +42,8 @@ def irradiate(contig, x=20):
     """
     return stdpopsim.Contig(
         recombination_map=contig.recombination_map,
-        gene_conversion_rate=contig.gene_conversion_rate,
+        bacterial_recombination=contig.bacterial_recombination,
+        gene_conversion_fraction=contig.gene_conversion_fraction,
         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_AnaPla.py b/tests/test_AnaPla.py
index 2dc92bfab..415491b10 100644
--- a/tests/test_AnaPla.py
+++ b/tests/test_AnaPla.py
@@ -163,3 +163,6 @@ def test_recombination_rate(self, name, rate):
     )
     def test_mutation_rate(self, name, rate):
         assert rate == pytest.approx(self.genome.get_chromosome(name).mutation_rate)
+
+    def test_bacterial_recombination(self):
+        assert self.genome.bacterial_recombination is False
diff --git a/tests/test_DroMel.py b/tests/test_DroMel.py
index a6098dec5..94fda08f0 100644
--- a/tests/test_DroMel.py
+++ b/tests/test_DroMel.py
@@ -29,8 +29,10 @@ def test_basic_attributes(self):
     def test_non_recombining_chrs(self, chr_id):
         chrom = self.genome.get_chromosome(chr_id)
         assert chrom.recombination_rate == 0
+        assert chrom.gene_conversion_fraction == 0
 
     @pytest.mark.parametrize("chr_id", ["2L", "2R", "3L", "3R", "X"])
     def test_recombining_chrs(self, chr_id):
         chrom = self.genome.get_chromosome(chr_id)
         assert chrom.recombination_rate > 0
+        assert chrom.gene_conversion_fraction > 0
diff --git a/tests/test_EscCol.py b/tests/test_EscCol.py
index 3b8526826..547c0aaa4 100644
--- a/tests/test_EscCol.py
+++ b/tests/test_EscCol.py
@@ -1,5 +1,5 @@
 """
-Tests for the e. coli data definitions.
+Tests for the EscCol data definitions.
 """
 import stdpopsim
 from tests import test_species
@@ -19,7 +19,7 @@ def test_basic_attributes(self):
 
 class TestGenome(test_species.GenomeTestBase):
     """
-    Tests for the e_coli genome.
+    Tests for the EscCol genome.
     """
 
     genome = stdpopsim.get_species("EscCol").genome
@@ -30,8 +30,11 @@ def test_basic_attributes(self):
     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_recombination_rate(self):
+        assert self.genome.get_chromosome("Chromosome").recombination_rate == 8.9e-11
 
     def test_gene_conversion_length(self):
-        assert self.genome.get_chromosome("Chromosome").gene_conversion_length == 345
+        assert self.genome.get_chromosome("Chromosome").gene_conversion_length == 542
+
+    def test_bacterial_recombination(self):
+        assert self.genome.bacterial_recombination is True
diff --git a/tests/test_HomSap.py b/tests/test_HomSap.py
index 85c9e369d..89acbeea4 100644
--- a/tests/test_HomSap.py
+++ b/tests/test_HomSap.py
@@ -37,3 +37,6 @@ def test_recombination_rates(self, chr_id):
         assert chrom.recombination_rate == pytest.approx(
             contig.recombination_map.mean_rate, rel=1e-6
         )
+
+    def test_bacterial_recombination(self):
+        assert self.genome.bacterial_recombination is False
diff --git a/tests/test_StrAga.py b/tests/test_StrAga.py
index eb2949149..a545a8972 100644
--- a/tests/test_StrAga.py
+++ b/tests/test_StrAga.py
@@ -31,7 +31,7 @@ class TestGenomeData(test_species.GenomeTestBase):
     @pytest.mark.parametrize(
         ["name", "rate"],
         {
-            "1": 0,
+            "1": 1.53e-10,
         }.items(),
     )
     def test_recombination_rate(self, name, rate):
diff --git a/tests/test_cli.py b/tests/test_cli.py
index 812474f2d..9a6e5b20a 100644
--- a/tests/test_cli.py
+++ b/tests/test_cli.py
@@ -1022,11 +1022,18 @@ def test_noQC_citations_not_written_verbose(self, caplog):
         )
 
 
-@pytest.mark.parametrize(
-    "species_id", [species.id for species in stdpopsim.all_species()]
-)
+@pytest.mark.parametrize("species_id", [s.id for s in stdpopsim.all_species()])
 def test_species_simulation(species_id):
-    cmd = f"-q {species_id} -L 1 --seed 1234 10"
+    # L must be longer than any mean gene conversion tract length;
+    # since gene conversion length is quite long for a few species
+    # we can't run with L=1000 for all species, so (confusingly) adjust
+    # for those long-GC species here
+    L = 1000
+    species = stdpopsim.get_species(species_id)
+    for chrom in species.genome.chromosomes:
+        if chrom.gene_conversion_length is not None:
+            L = max(L, chrom.gene_conversion_length + 100)
+    cmd = f"-q {species_id} -L {L} --seed 1234 10"
     # Just check to see if the simulation runs
     with mock.patch("sys.stdout", autospec=True) as stdout:
         stdout.buffer = open(os.devnull, "wb")
diff --git a/tests/test_engines.py b/tests/test_engines.py
index 4b2ac9d69..55be1cd3b 100644
--- a/tests/test_engines.py
+++ b/tests/test_engines.py
@@ -143,3 +143,13 @@ def test_non_neutral_contig(self):
         contig.add_dfe(intervals=np.array([[0, 50]]), DFE=dfes[0])
         engine = stdpopsim.get_engine("msprime")
         engine.simulate(model, contig, samples, seed=1)
+
+    def test_gene_conversion(self):
+        species = stdpopsim.get_species("DroMel")
+        model = species.get_demographic_model("African3Epoch_1S16")
+        samples = model.get_samples(10)
+        contig = species.get_contig(length=1000, use_species_gene_conversion=True)
+        assert contig.gene_conversion_fraction > 0
+        assert contig.gene_conversion_length > 0
+        engine = stdpopsim.get_engine("msprime")
+        engine.simulate(model, contig, samples, seed=1)
diff --git a/tests/test_genomes.py b/tests/test_genomes.py
index 7d17b0d8b..e060ac6a7 100644
--- a/tests/test_genomes.py
+++ b/tests/test_genomes.py
@@ -3,7 +3,7 @@
 """
 import numpy as np
 import pytest
-import copy
+import msprime
 import stdpopsim
 
 
@@ -28,8 +28,53 @@ def verify_dfe_breakpoints(self, contig):
 
     def test_default_gc(self):
         contig = stdpopsim.Contig.basic_contig(length=42)
-        assert contig.gene_conversion_rate is None
+        assert contig.gene_conversion_fraction is None
         assert contig.gene_conversion_length is None
+        assert contig.bacterial_recombination is False
+
+    def test_gc_errors(self):
+        # should not error
+        contig = stdpopsim.Contig.basic_contig(
+            length=100, bacterial_recombination=True, gene_conversion_length=14
+        )
+        assert contig.bacterial_recombination is True
+        contig = stdpopsim.Contig.basic_contig(
+            length=100, gene_conversion_fraction=1.0, gene_conversion_length=14
+        )
+        assert contig.gene_conversion_fraction == 1.0
+        assert contig.gene_conversion_length == 14
+        # can't set gc fraction for bacterial recomb
+        with pytest.raises(ValueError, match="Cannot set.*bacterial recomb"):
+            _ = stdpopsim.Contig.basic_contig(
+                length=100,
+                bacterial_recombination=True,
+                gene_conversion_fraction=0.4,
+                gene_conversion_length=20,
+            )
+        # must set gc length for bacterial recomb
+        with pytest.raises(ValueError, match="Must set.*bacterial recomb"):
+            _ = stdpopsim.Contig.basic_contig(length=100, bacterial_recombination=True)
+        # must set gc length and fraction together otherwise
+        with pytest.raises(ValueError, match="without setting"):
+            _ = stdpopsim.Contig.basic_contig(length=100, gene_conversion_fraction=0.4)
+        with pytest.raises(ValueError, match="without setting"):
+            _ = stdpopsim.Contig.basic_contig(length=100, gene_conversion_length=14)
+        # gc_frac must be between 0 and 1 (inclusive)
+        for bad_frac in [-0.1, 1.1]:
+            with pytest.raises(ValueError, match="must be between 0 and 1"):
+                _ = stdpopsim.Contig.basic_contig(
+                    length=100,
+                    gene_conversion_fraction=bad_frac,
+                    gene_conversion_length=14,
+                )
+        # gc_length must be at least 1
+        for bad_len in [-0.1, -10]:
+            with pytest.raises(ValueError, match="must be greater than 1"):
+                _ = stdpopsim.Contig.basic_contig(
+                    length=100,
+                    gene_conversion_fraction=0.8,
+                    gene_conversion_length=bad_len,
+                )
 
     def test_default_dfe(self):
         contig = stdpopsim.Contig.basic_contig(length=100)
@@ -323,38 +368,75 @@ def test_chromosome_segment_exceeds_length(self):
 
 
 class TestGeneConversion(object):
-    species = stdpopsim.get_species("EscCol")
+    def test_mean_gene_conversion(self):
+        dro_mel = stdpopsim.get_species("DroMel")
+        contig = dro_mel.get_contig(length=1000, use_species_gene_conversion=False)
+        assert contig.gene_conversion_fraction is None
+        assert contig.gene_conversion_length is None
+        contig = dro_mel.get_contig(length=1000, use_species_gene_conversion=True)
+        mean_gc = dro_mel.genome.mean_gene_conversion_fraction
+        assert np.isclose(mean_gc, contig.gene_conversion_fraction)
 
-    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_no_mean_gene_conversion(self):
+        # this will need to be changed if we add GC rates to AnaPla
+        ana_pla = stdpopsim.get_species("AnaPla")
+        contig = ana_pla.get_contig(length=1000, use_species_gene_conversion=True)
+        mean_gc = ana_pla.genome.mean_gene_conversion_fraction
+        assert mean_gc == 0.0
+        assert contig.gene_conversion_fraction is None
 
-    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_modifiying_gene_conversion(self):
+        esc_col = stdpopsim.get_species("EscCol")
+        assert esc_col.genome.mean_gene_conversion_fraction == 0
+        contig = esc_col.get_contig(chromosome="Chromosome")
+        assert contig.bacterial_recombination is True
+        assert contig.gene_conversion_fraction is None
+        assert contig.gene_conversion_length > 0
 
-    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
+        rm = msprime.RateMap(position=[0.0, contig.length], rate=[0.1])
+        contig.recombination_map = rm
+        contig.gene_conversion_length = 1
+        assert contig.recombination_map == rm
+        assert contig.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_errors(self):
+        esc_col = stdpopsim.get_species("EscCol")
+        with pytest.raises(ValueError, match="with bacterial recombination"):
+            _ = esc_col.get_contig(
+                chromosome="Chromosome", use_species_gene_conversion=True
+            )
 
     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
+        species = stdpopsim.get_species("AnaPla")
+        contig = species.get_contig("1")
+        assert contig.gene_conversion_fraction is None
+        assert contig.gene_conversion_length is None
+        assert contig.bacterial_recombination is False
+        contig = species.get_contig(length=42, use_species_gene_conversion=True)
+        assert contig.gene_conversion_fraction is None
+        assert contig.gene_conversion_length is None
+        assert contig.bacterial_recombination is False
+
+    def test_get_species_contig_rates(self):
+        species = stdpopsim.get_species("AnaPla")
+        chrom = species.genome.chromosomes[0]
+        gc_frac = 0.7
+        gc_len = 123
+        # this is *not* the recommended workflow; just doing it for testing
+        # see test_modifiying_gene_conversion for the recommended way to modify it
+        chrom.gene_conversion_fraction = gc_frac
+        chrom.gene_conversion_length = gc_len
+        assert species.genome.chromosomes[0].gene_conversion_fraction == gc_frac
+        assert species.genome.chromosomes[0].gene_conversion_length == gc_len
+        # without use_species_gene_conversion
+        contig = species.get_contig(chrom.id)
+        assert contig.gene_conversion_fraction is None
+        assert contig.gene_conversion_length is None
+        assert contig.recombination_map.mean_rate == chrom.recombination_rate
+        # now, with use_species_gene_conversion
+        contig = species.get_contig(chrom.id, use_species_gene_conversion=True)
+        assert contig.gene_conversion_fraction == gc_frac
+        assert contig.gene_conversion_length == gc_len
+        assert np.isclose(
+            contig.recombination_map.mean_rate, chrom.recombination_rate / (1 - gc_frac)
+        )
diff --git a/tests/test_species.py b/tests/test_species.py
index 51d857b65..4ce72eebb 100644
--- a/tests/test_species.py
+++ b/tests/test_species.py
@@ -94,10 +94,9 @@ def test_get_unknown_annotation(self):
             with pytest.raises(ValueError):
                 species.get_annotations(name)
 
-    @pytest.mark.xfail  # HomSap annotation not currently available
     def test_add_duplicate_annotation(self):
         species = stdpopsim.get_species("HomSap")
-        an = species.get_annotations("Ensembl_GRCh38_gff3")
+        an = species.annotations[0]
         with pytest.raises(ValueError):
             species.add_annotations(an)
 
@@ -154,11 +153,27 @@ def test_population_size_defined(self):
     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
+            assert contig.gene_conversion_fraction is None
+            assert (
+                contig.bacterial_recombination
+                == self.species.genome.bacterial_recombination
+            )
+            if contig.bacterial_recombination is True:
+                assert chrom.gene_conversion_length == contig.gene_conversion_length
+            else:
+                assert contig.gene_conversion_length is None
+            if self.species.genome.bacterial_recombination is False:
+                contig = self.species.get_contig(
+                    chrom.id, use_species_gene_conversion=True
+                )
+                assert chrom.gene_conversion_fraction == contig.gene_conversion_fraction
+                assert chrom.gene_conversion_length == contig.gene_conversion_length
+                assert contig.bacterial_recombination is False
+            else:
+                with pytest.raises(ValueError, match="with bacterial recombination"):
+                    contig = self.species.get_contig(
+                        chrom.id, use_species_gene_conversion=True
+                    )
 
 
 class GenomeTestBase:
@@ -290,30 +305,16 @@ def test_contig_options(self):
             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)
-        # TODO: invalid interval
 
     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
+        if self.species.genome.get_chromosome("chr22").gene_conversion_fraction is None:
+            assert contig.gene_conversion_fraction is None
             assert contig.gene_conversion_length is None
         else:
             assert (
-                contig.gene_conversion_rate
-                == self.species.genome.get_chromosome("chr22").gene_conversion_rate
+                contig.gene_conversion_fraction
+                == self.species.genome.get_chromosome("chr22").gene_conversion_fraction
             )
             assert (
                 contig.gene_conversion_length
@@ -322,38 +323,67 @@ def test_use_species_gene_conversion(self):
 
     def test_generic_contig(self):
         L = 1e6
-        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")
-        Ls = [c.length for c in self.species.genome.chromosomes if c.id in chrom_ids]
-        rs = [
-            c.recombination_rate
-            for c in self.species.genome.chromosomes
-            if c.id in chrom_ids
-        ]
-        us = [
-            c.mutation_rate
-            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)
+        for usgc in [True, False]:
+            contig = self.species.get_contig(length=L, use_species_gene_conversion=usgc)
+            assert contig.recombination_map.sequence_length == L
+
+            chrom_ids = np.arange(1, 23).astype("str")
+            Ls = np.array(
+                [c.length for c in self.species.genome.chromosomes if c.id in chrom_ids]
+            )
+            rs = np.array(
+                [
+                    c.recombination_rate
+                    for c in self.species.genome.chromosomes
+                    if c.id in chrom_ids
+                ]
+            )
+            us = np.array(
+                [
+                    c.mutation_rate
+                    for c in self.species.genome.chromosomes
+                    if c.id in chrom_ids
+                ]
+            )
+            gcs = np.array(
+                [
+                    c.gene_conversion_fraction
+                    if c.gene_conversion_fraction is not None
+                    else 0
+                    for c in self.species.genome.chromosomes
+                    if c.id in chrom_ids
+                ]
+            )
+            gcls = np.array(
+                [
+                    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)
+            if usgc:
+                assert contig.recombination_map.mean_rate == np.average(
+                    rs, weights=Ls
+                ) / (1 - np.average(gcs, weights=Ls * rs))
+            else:
+                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
+                or not usgc
+            ):
+                assert contig.gene_conversion_fraction is None
+                assert contig.gene_conversion_length is None
+            else:
+                assert contig.gene_conversion_fraction is not None
+                assert np.isclose(
+                    contig.gene_conversion_fraction, np.average(gcs, weights=Ls * rs)
+                )
+                assert contig.gene_conversion_length is not None
+                assert np.isclose(
+                    contig.gene_conversion_length, np.average(gcls, weights=Ls * rs)
+                )
         # TODO: origin