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

Issue:
When using `left`/`right`, output coordinate system is relative to start of contig
The options `left`, `right` in `Species.get_contig` allow a portion of a chromosome to be simulated (in contrast to simulating an entire chromosome and "masking" everything outside of `[left, right]`). When adding DFEs/inclusion masks/sweeps/etc to the contig, the coordinate system used is that of the original chromosome.

However, the coordinate system of the output tree sequence will be relative to the start of the contig, not of the original chromosome. So for example, if I simulate chr22 from 30 Mb to 50 Mb then dump a VCF, the coordinates of variants won't reflect their positions on the original chromosome.

Instead, should the output tree sequence span the entire chromosome and have missing data outside of `[left, right]`?
When using `left`/`right`, output coordinate system is relative to start of contig
The options `left`, `right` in `Species.get_contig` allow a portion of a chromosome to be simulated (in contrast to simulating an entire chromosome and "masking" everything outside of `[left, right]`). When adding DFEs/inclusion masks/sweeps/etc to the contig, the coordinate system used is that of the original chromosome.

However, the coordinate system of the output tree sequence will be relative to the start of the contig, not of the original chromosome. So for example, if I simulate chr22 from 30 Mb to 50 Mb then dump a VCF, the coordinates of variants won't reflect their positions on the original chromosome.

Instead, should the output tree sequence span the entire chromosome and have missing data outside of `[left, right]`?
Restarting simulations in SLiM engine (community.allScriptBlocks.active property)
Hi all, first, I want to congratulate you on such a fantastic project you're working on. 

I've been using Stdpopsim to automatize sweep simulations given the time (_t_) at which beneficial mutation arises, the final frequency (_f_), and the selection coefficient (_s_) range that can produce such a sweep. In my approach, I solved the Otto and Whitelock 2013 Eq. 5 to approximate the time of fixation or the desired final frequency to obtain proper _t_ and _s_ combinations producing the expected sweep. Because SLiM + Stdpopsim restarts the simulation as much as needed, such an approximation works well (I've tested constant population size and OutOfAfrica_2T12 demography models).

I was using SLiM 4.0.1, and everything worked fine. Today, I upgraded SLiM to version 4.2.1, and the simulations don't work. I dumped the SLiM recipe and tested it in different environments in addition to trying to debug it in SLiMgui. Nonetheless, it does not work on SLiM 4.2.1 but still works on 4.0.1.

The error I obtained is the following one:

```bash 
ERROR (EidosPropertySignature::CheckAssignedValue): value cannot be type logical for read-write property active.  
Error on script line 25, character 33 (inside runtime script block):  
    sb[sb.type == "late"].active = F;  
                                 ^
```

I checked the SLiM manual, and now the `community.allScriptBlocks.active` property is defined as an integer. The problem should be related to assigning a boolean value instead of an integer. 

I changed the line ` sb[sb.type == "late"].active = F;` to  sb[sb.type == "late"].active = 0;` and tested in SLiM and SLiMgui, and it works now as expected. There is no error, and I tracked the mutation frequency over time using the SLiMgui, and the beneficial mutation reached the proper frequency.

**How to reproduce the problem:**
I am attaching a piece of code to simulate a sweep with constant population size and uniform recombination rate in a 1.2MB region. Timing, selection coefficient, and initial and final frequency are defined to produce an ongoing sweep that started ~1300 generations ago and reaches a final frequency in the range [0.8,0.9].

```python
import stdpopsim
import msprime
import numpy as np
from contextlib import redirect_stdout

slim_script = "/home/jmurgamoreno/test.slim"
# selection coefficient, time to beneficial, end time to partial sweep if needed, burning, iteration number
(_s, _t, _t_end, _f_i, _f_t, _burn_in) = (
    0.0119,
    1384,
    0,
    0.0001,
    0.8,
    0.15000000000000002,
)


# Retrieve species information
species = stdpopsim.get_species("HomSap")
model = stdpopsim.PiecewiseConstantSize(species.population_size)
model.generation_time = species.generation_time
# Force to sample to pop0 is PiecewiseConstantSize selected
sample = {"pop_0": 50}

contig = species.get_contig("chr2", mutation_rate=model.mutation_rate, right=1.2e6)

# Change manually to test og Flex-sweep
contig.mutation_rate = np.random.uniform(2 * 1e-9, 5.2 * 1e-8)
contig.recombination_map = msprime.RateMap(
    position=contig.recombination_map.position,
    rate=np.array([np.random.exponential(1e-8)]),
)

# adaptive_position = np.random.randint(1.2e6*0.25,1.2e6*0.75)
adaptive_position = sum(contig.original_coordinates[1:]) / 2

contig.add_single_site(id="hard", coordinate=adaptive_position)


extended_events = stdpopsim.ext.selective_sweep(
    single_site_id="hard",
    population=list(sample.keys())[0],
    selection_coeff=_s,
    min_freq_at_end=_f_t,
    mutation_generation_ago=_t,
    end_generation_ago=_t_end,
)

if _f_t != 1:
    end_condition = stdpopsim.ext.ConditionOnAlleleFrequency(
        start_time=_t_end,
        end_time=_t_end,
        single_site_id="hard",
        population=list(sample.keys())[0],
        op="<",
        allele_frequency=_f_t + 0.1,
    )

    extended_events.append(end_condition)

engine = stdpopsim.get_engine("slim")


with open(slim_script, "w") as f:
    with redirect_stdout(f):
        ts_sweep = engine.simulate(
            model,
            contig,
            sample,
            extended_events=extended_events,
            slim_scaling_factor=1,
            slim_burn_in=_burn_in,
            slim_script=True,
            verbosity=0,
        )

ts_sweep = engine.simulate(
    model,
    contig,
    sample,
    extended_events=extended_events,
    slim_scaling_factor=1,
    slim_burn_in=_burn_in,
)
```
**Suggested fixes:**
The `community.allScriptBlocks.active` property seems unable to deal with boolean values now. Using 0 values as `False` solves the problem. So, changing line 219 on `slim_engine.py` is a potential solution in case I am not doing anything wrong.

Best,

Jesús Murga-Moreno





Discussion / hints:
> Instead, should the output tree sequence span the entire chromosome and have missing data outside of [left, right]?

Probably yes I think, it's the less bad option. 
Should be able to do this by removing `trim=True` from [here](https://github.com/popsim-consortium/stdpopsim/blob/b1383cce4eac2a0fea92b660fdca88b8c4ad5928/stdpopsim/genomes.py#L507), and removing the coordinate "shifting" for DFEs and exclusion masks, I think.
We discussed this a bit in the hackathon at the meeting before prob gen and my argument was that for many uses people will be interested in generating simulated data that "look like" some region and then plugging that into some pipeline that would want coordinates starting with 1. (E.g. calculating summary stats on a bunch of sims before running ABC.)

So one option is to leave the current behavior as is, while warning the user that the rest of the chrom is not simulated, and maybe add in another flag for having the tree seq span the whole chrom and doing the missing data thing as proposed above.
If we're outputting the original coordinate system with "flanks" outside of the simulated window set to missing, then trimming to the window (thus resetting the coordinate system) is as simple as ``ts_new_coord = ts.trim()``. So, I'm not sure we'd need to provide an option for this in the python API, we could just mention ``trim`` method in the docstring.

We could add a ``--trim-missing`` flag to the CLI, that'd remove missing/masked flanks (this would also trim off parts of chromosomes where recombination maps are NaN, like the first 15Mb of chr22).
I have just encountered this problem and it confused me, because I wanted to add a reference sequence to the simulated genome. I agree with @nspope and @jeromekelleher that we should not "trim by default", and we should simply have empty flanks (empty trees with `tree.num_edges=0`) if we specify `left>0` and `right < seq_len`. This will also preserve the coordinate system of the simulated chromosome, in case the positions are used e.g. for annotations etc.
In doing this, we should rename the `original_coordinates` attribute of Contigs (and maybe keep `.original_coordinates` as a deprecated alias?)
> Instead, should the output tree sequence span the entire chromosome and have missing data outside of [left, right]?

Probably yes I think, it's the less bad option. 
Should be able to do this by removing `trim=True` from [here](https://github.com/popsim-consortium/stdpopsim/blob/b1383cce4eac2a0fea92b660fdca88b8c4ad5928/stdpopsim/genomes.py#L507), and removing the coordinate "shifting" for DFEs and exclusion masks, I think.
We discussed this a bit in the hackathon at the meeting before prob gen and my argument was that for many uses people will be interested in generating simulated data that "look like" some region and then plugging that into some pipeline that would want coordinates starting with 1. (E.g. calculating summary stats on a bunch of sims before running ABC.)

So one option is to leave the current behavior as is, while warning the user that the rest of the chrom is not simulated, and maybe add in another flag for having the tree seq span the whole chrom and doing the missing data thing as proposed above.
If we're outputting the original coordinate system with "flanks" outside of the simulated window set to missing, then trimming to the window (thus resetting the coordinate system) is as simple as ``ts_new_coord = ts.trim()``. So, I'm not sure we'd need to provide an option for this in the python API, we could just mention ``trim`` method in the docstring.

We could add a ``--trim-missing`` flag to the CLI, that'd remove missing/masked flanks (this would also trim off parts of chromosomes where recombination maps are NaN, like the first 15Mb of chr22).
I have just encountered this problem and it confused me, because I wanted to add a reference sequence to the simulated genome. I agree with @nspope and @jeromekelleher that we should not "trim by default", and we should simply have empty flanks (empty trees with `tree.num_edges=0`) if we specify `left>0` and `right < seq_len`. This will also preserve the coordinate system of the simulated chromosome, in case the positions are used e.g. for annotations etc.
In doing this, we should rename the `original_coordinates` attribute of Contigs (and maybe keep `.original_coordinates` as a deprecated alias?)
Thanks for the very thorough bug report @jmurga!

@petrelharp do you see any potential issues with the solution above? e.g. do we need to handle back-compatibility with earlier slim versions?
Looking at the help (in SLiMgui), the documentation says:
```
If this evaluates to logical F (i.e., is equal to 0), the script block is inactive and will not be called.  The value of active for all registered script blocks is reset to -1 at the beginning of each tick, prior to script events being called, thus activating all blocks (except callbacks associated with a species that is not active in that tick, which are deactivated as part of the deactivation of the species).  Any integer value other than -1 may be used instead of -1 to represent that a block is active; for example, active may be used as a counter to make a block execute a fixed number of times in each tick.  This value is not cached by SLiM; if it is changed, the new value takes effect immediately.  For example, a callback might be activated and inactivated repeatedly during a single tick.
```
... so, it looks like the type got changed to allow more behavior (in particular, using `-1` for something). This might actually be a bit of unintended backwards incompatibility? FYI, @bhaller.

But - I don't see any problem with this solution, no. Thanks very much for the report!
> Looking at the help (in SLiMgui), the documentation says:
> 
> ```
> If this evaluates to logical F (i.e., is equal to 0), the script block is inactive and will not be called.  The value of active for all registered script blocks is reset to -1 at the beginning of each tick, prior to script events being called, thus activating all blocks (except callbacks associated with a species that is not active in that tick, which are deactivated as part of the deactivation of the species).  Any integer value other than -1 may be used instead of -1 to represent that a block is active; for example, active may be used as a counter to make a block execute a fixed number of times in each tick.  This value is not cached by SLiM; if it is changed, the new value takes effect immediately.  For example, a callback might be activated and inactivated repeatedly during a single tick.
> ```
> 
> ... so, it looks like the type got changed to allow more behavior (in particular, using `-1` for something). This might actually be a bit of unintended backwards incompatibility? FYI, @bhaller.
> 
> But - I don't see any problem with this solution, no. Thanks very much for the report!

From the 4.2 release notes:

> Assignment into a property of an object must now match the type of the property; type promotion will no longer occur. Existing code that now errors should use asInteger(), asFloat(), etc. as needed.

It broke a few patterns like this; but `active` has actually always been type `integer`, not type `logical`.  Tightening up implicit type coercion in Eidos a bit.  Sorry!

Fix patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 8687235f5..9e13be4ec 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -2,8 +2,28 @@
 [0.2.1] - 2023-XX-XX
 --------------------
 
+**Bug fixes**:
+
+- Updates to SLiM support: updated the `active` flag in the SLiM code to be integer.`
+
 **Breaking changes**:
 
+- The `time_units` attribute of tree sequence metadata is now set to "generations"
+    for SLiM output, avoiding the "time units mismatch warning".
+    (:user:`nspope`, :pr:`1567`)
+
+- Previously, Contigs specified with `left` and `right` arguments produced
+  simulations having coordinates shifted so they were relative to `left`. This
+  is no longer the case: coordinates in resulting simulations retain the
+  original chromosome's coordinates. As a result, regions outside `[left,
+  right)` will now have missing data. To produce simulation output with the
+  previous behavior, use the `.trim()` method of the resulting TreeSequence.
+  So, the Contig's `.original_coordinates` attribute is now called
+  `.coordinates`, and several methods (e.g., `Contig.basic_contig()`) now take
+  `left` and `right` arguments as well as `length`. Finally, arguments such as
+  `relative_coordinates=True/False` to `Contig.dfe_breakpoints()` are no longer
+  necessary and are deprecated. (:user:`petrelharp`, :pr:`1570`)
+
 - To add the possibility for a relationship between dominance and selection
     coefficient, now each stdpopsim MutationType might have more than one
     underlying SLiM mutation type; so, where this is recorded in top-level
diff --git a/stdpopsim/cli.py b/stdpopsim/cli.py
index 89747a544..c27268745 100644
--- a/stdpopsim/cli.py
+++ b/stdpopsim/cli.py
@@ -813,9 +813,9 @@ def run_simulation(args):
                 intervals_summary_str = f"[{left}, {right})"
             if intervals_summary_str is None:
                 # case where no intervals specified but we have a DFE
-                _, left, right = contig.original_coordinates
+                _, left, right = contig.coordinates
                 intervals = np.array([[left, right]])
-                intervals_summary_str = f"[{intervals[0][0]}, {intervals[0][1]})"
+                intervals_summary_str = f"[{left}, {right})"
 
             dfe = species.get_dfe(args.dfe)
             contig.add_dfe(
diff --git a/stdpopsim/genomes.py b/stdpopsim/genomes.py
index b73d8911d..5e9257227 100644
--- a/stdpopsim/genomes.py
+++ b/stdpopsim/genomes.py
@@ -14,6 +14,25 @@
 logger = logging.getLogger(__name__)
 
 
+def _uniform_rate_map(left, right, length, rate):
+    """
+    Makes a RateMap that is
+    """
+    pos = [
+        0.0,
+    ]
+    x = []
+    if left > 0:
+        pos.append(left)
+        x.append(np.nan)
+    pos.append(right)
+    x.append(rate)
+    if right < length:
+        pos.append(length)
+        x.append(np.nan)
+    return msprime.RateMap(position=pos, rate=x)
+
+
 @attr.s
 class Genome:
     """
@@ -225,9 +244,6 @@ class Contig:
     and ``interval_list``. These must be of the same length, and the k-th DFE
     applies to the k-th interval; see :meth:`.add_dfe` for more information.
 
-    The contig may be a segment of a named chromosome. If so, the original
-    coordinate system is used for :meth:`.add_dfe` and :meth:`add_single_site`.
-
     :ivar mutation_rate: The rate of mutation per base per generation.
     :vartype mutation_rate: float
     :ivar ploidy: The ploidy of the contig. Defaults to 2 (diploid).
@@ -267,11 +283,11 @@ class Contig:
     :ivar interval_list: A list of :class:`np.array` objects containing integers.
         By default, the inital interval list spans the whole chromosome with the
         neutral DFE.
-    :ivar original_coordinates: The location of the contig on a named chromosome,
+    :ivar coordinates: The location of the contig on a named chromosome,
         as a tuple of the form `(chromosome, left, right)`. If `None`, the contig
         is assumed to be generic (i.e. it does not inherit a coordinate system
         from a larger chromosome).
-    :vartype original_coordinates: tuple
+    :vartype coordinates: tuple
 
     .. note::
         To run stdpopsim simulations with alternative, user-specified mutation,
@@ -300,14 +316,14 @@ class Contig:
     exclusion_mask = attr.ib(default=None)
     dfe_list = attr.ib(factory=list)
     interval_list = attr.ib(factory=list)
-    original_coordinates = attr.ib(default=None, type=tuple)
+    coordinates = attr.ib(default=None, type=tuple)
 
     def __attrs_post_init__(self):
-        if self.original_coordinates is None:
-            self.original_coordinates = (None, 0, int(self.length))
-        _, left, right = self.original_coordinates
+        if self.coordinates is None:
+            self.coordinates = (None, 0, int(self.length))
+        _, left, right = self.coordinates
         self.add_dfe(
-            [[left, right]],
+            np.array([[left, right]]),
             stdpopsim.dfe.neutral_dfe(),
         )
 
@@ -491,20 +507,24 @@ def species_contig(
                     logger.debug(
                         f"Making flat chromosome {length_multiplier} * {chrom.id}"
                     )
+                    recomb_map = msprime.RateMap.uniform(
+                        round(length_multiplier * chrom.length),
+                        chrom.recombination_rate,
+                    )
                 else:
                     logger.debug(
-                        f"Making flat contig of length {right - left} from {chrom.id}"
+                        f"Making flat contig from {left} to {right} for {chrom.id}"
+                    )
+                    recomb_map = _uniform_rate_map(
+                        left, right, chrom.length, chrom.recombination_rate
                     )
-                recomb_map = msprime.RateMap.uniform(
-                    right - left, chrom.recombination_rate
-                )
             else:
                 if length_multiplier != 1:
                     raise ValueError("Cannot use length multiplier with empirical maps")
                 logger.debug(f"Getting map for {chrom.id} from {genetic_map}")
                 gm = species.get_genetic_map(genetic_map)
                 recomb_map = gm.get_chromosome_map(chrom.id)
-                recomb_map = recomb_map.slice(left=left, right=right, trim=True)
+                recomb_map = recomb_map.slice(left=left, right=right)
 
             inclusion_intervals = None
             exclusion_intervals = None
@@ -517,7 +537,7 @@ def species_contig(
                     )
                 else:
                     inclusion_intervals = inclusion_mask
-                inclusion_intervals = stdpopsim.utils.clip_and_shift_intervals(
+                inclusion_intervals = stdpopsim.utils.clip_intervals(
                     inclusion_intervals, left, right
                 )
             if exclusion_mask is not None:
@@ -529,7 +549,7 @@ def species_contig(
                     )
                 else:
                     exclusion_intervals = exclusion_mask
-                exclusion_intervals = stdpopsim.utils.clip_and_shift_intervals(
+                exclusion_intervals = stdpopsim.utils.clip_intervals(
                     exclusion_intervals, left, right
                 )
 
@@ -562,7 +582,7 @@ def species_contig(
                 gene_conversion_length=gene_conversion_length,
                 inclusion_mask=inclusion_intervals,
                 exclusion_mask=exclusion_intervals,
-                original_coordinates=(chromosome, left, right),
+                coordinates=(chromosome, left, right),
                 ploidy=ploidy,
             )
 
@@ -572,19 +592,30 @@ def species_contig(
     def length(self):
         return self.recombination_map.sequence_length
 
+    @property
+    def original_coordinates(self):
+        warnings.warn(
+            stdpopsim.DeprecatedFeatureWarning(
+                "The original_coordinates property is deprecated, "
+                "since contigs are no longer shifted to be relative to their "
+                "left endpoint: use Contig.coordinates instead."
+            )
+        )
+        return self.coordinates
+
     @property
     def origin(self):
         """
         The location of the contig on a named chromosome as a string with format,
         "chromosome:left-right"; or None if a generic contig.
         """
-        chromosome, left, right = self.original_coordinates
+        chromosome, left, right = self.coordinates
         if chromosome is None:
             return None
         else:
             return f"{chromosome}:{left}-{right}"
 
-    def dfe_breakpoints(self, relative_coordinates=True):
+    def dfe_breakpoints(self, *, relative_coordinates=None):
         """
         Returns two things: the sorted vector of endpoints of all intervals across
         all DFEs in the contig, and a vector of integer labels for these intervals,
@@ -607,21 +638,23 @@ def dfe_breakpoints(self, relative_coordinates=True):
         ``breaks[k]``.  Some intervals may not be covered by a DFE, in which
         case they will have the label ``-1`` (beware of python indexing!).
 
-        :param bool relative_coordinates: If True, the returned breakpoints
-            will be relative to the start of the contig, rather than to the
-            start of the chromosome to which to contig belongs.
-
         :return: A tuple (breaks, dfe_labels).
         """
+        if relative_coordinates is not None:
+            warnings.warn(
+                stdpopsim.DeprecatedFeatureWarning(
+                    "The relative_coordinates argument is deprecated "
+                    "(and no longer needed),"
+                    "because contigs are no longer shifted to be relative to their "
+                    "left endpoint."
+                )
+            )
         breaks = np.unique(
             np.vstack(self.interval_list + [[[0, int(self.length)]]])  # also sorted
         )
         dfe_labels = np.full(len(breaks) - 1, -1, dtype="int")
         for j, intervals in enumerate(self.interval_list):
             dfe_labels[np.isin(breaks[:-1], intervals[:, 0], assume_unique=True)] = j
-        if not relative_coordinates:
-            _, left, _ = self.original_coordinates
-            breaks += left
         return breaks, dfe_labels
 
     def clear_dfes(self):
@@ -659,9 +692,9 @@ def add_dfe(self, intervals, DFE):
         :param array intervals: A valid set of intervals.
         :param DFE dfe: A DFE object.
         """
-        _, left, right = self.original_coordinates
-        intervals = stdpopsim.utils.clip_and_shift_intervals(intervals, left, right)
-        stdpopsim.utils._check_intervals_validity(intervals, start=0, end=self.length)
+        _, left, right = self.coordinates
+        intervals = stdpopsim.utils.clip_intervals(intervals, left, right)
+        stdpopsim.utils._check_intervals_validity(intervals, start=left, end=right)
         for j, ints in enumerate(self.interval_list):
             self.interval_list[j] = stdpopsim.utils.mask_intervals(ints, intervals)
         self.dfe_list.append(DFE)
diff --git a/stdpopsim/slim_engine.py b/stdpopsim/slim_engine.py
index 1fa867d2d..5c1a114ce 100644
--- a/stdpopsim/slim_engine.py
+++ b/stdpopsim/slim_engine.py
@@ -216,7 +216,7 @@
      * start of the next tick.
      */
     sb = community.allScriptBlocks;
-    sb[sb.type == "late"].active = F;
+    sb[sb.type == "late"].active = 0;
 
     /*
      * 2. The late events below were run in the save() tick,
diff --git a/stdpopsim/species.py b/stdpopsim/species.py
index e477ccf62..1fef23878 100644
--- a/stdpopsim/species.py
+++ b/stdpopsim/species.py
@@ -184,6 +184,17 @@ def get_contig(
         is to be simulated based on empirical information for a given species
         and chromosome.
 
+        *Coordinates:* If a particular chunk of a given chromosome is obtained
+        (by specifying a `chromosome` ID and `left` and `right` coordinates),
+        the resulting genomes will have coordinates matching the original
+        (full) genome, and portions of the genome outside of `[left, right)`
+        will be masked (and so contain missing data). So, coordinates of
+        `genetic_map` and masks should be in coordinates of the original genome
+        (so, not shifted to be relative to `left`). If it is desired to have
+        output coordinates relative to `left`, use the
+        :meth:`tskit.TreeSequence.trim` method on the result (for instance,
+        `ts_shifted = ts.trim()`).
+
         :param str chromosome: The ID of the chromosome to simulate.
             A complete list of chromosome IDs for each species can be found in the
             "Genome" subsection for the species in the :ref:`sec_catalog`.
@@ -218,11 +229,11 @@ def get_contig(
             and recombination rates are equal to the genome-wide average across all
             autosomal chromosomes.
         :param float left: The left coordinate (inclusive) of the region to
-            keep on the chromosome. Defaults to 0.  Genetic maps and masks are
-            clipped to this boundary.
+            keep on the chromosome. Defaults to 0. Remaining regions will have
+            missing data when simulated.
         :param float right: The right coordinate (exclusive) of the region to
             keep on the chromosome. Defaults to the length of the chromosome.
-            Genetic maps and masks are clipped to this boundary.
+            Remaining regions will have missing data when simulated.
         :rtype: :class:`.Contig`
         :return: A :class:`.Contig` describing the section of the genome.
         """
diff --git a/stdpopsim/utils.py b/stdpopsim/utils.py
index 714e18bf0..4b3ea9794 100644
--- a/stdpopsim/utils.py
+++ b/stdpopsim/utils.py
@@ -254,12 +254,11 @@ def mask_intervals(intervals, mask):
     return out
 
 
-def clip_and_shift_intervals(intervals, left, right):
+def clip_intervals(intervals, left, right):
     """
     Intersect each interval in ``intervals`` with ``[left, right]``.
     ``intervals`` should be a numpy array with two columns that are the
-    beginning and end coordinates, respectively. The coordinates are then
-    shifted such that ``left`` becomes 0.
+    beginning and end coordinates, respectively.
     """
     assert 0 <= left < right
     intervals = np.array(intervals).astype(int)
@@ -268,7 +267,6 @@ def clip_and_shift_intervals(intervals, left, right):
     intervals = np.clip(intervals[~out_of_bounds], left, right)
     if intervals.shape[0] == 0:
         warnings.warn(f"No intervals remain after clipping to [{left}, {right}]")
-    intervals -= left
     return intervals
 
 


Test patch:
diff --git a/tests/test_genomes.py b/tests/test_genomes.py
index 61beb6a1b..48500329b 100644
--- a/tests/test_genomes.py
+++ b/tests/test_genomes.py
@@ -203,16 +203,17 @@ def test_is_neutral(self):
                 assert contig.is_neutral is (neutral and dist == "f")
 
     def test_chromosome_segment(self):
-        chrom = "chr2"
-        interval = [100, 130]
-        species = stdpopsim.get_species("HomSap")
-        contig = species.get_contig(chrom, left=interval[0], right=interval[1])
-        chromosome, left, right = contig.original_coordinates
-        assert chromosome == chrom
-        assert left == interval[0]
-        assert right == interval[1]
-        assert contig.length == interval[1] - interval[0]
-        assert contig.origin == f"{chromosome}:{left}-{right}"
+        species = stdpopsim.get_species("AnaPla")
+        chrom = species.genome.chromosomes[1]
+        length = chrom.length
+        for interval in [(0, 1341), (5020, 12850), (4249, length), (0, length)]:
+            contig = species.get_contig(chrom.id, left=interval[0], right=interval[1])
+            chromosome, left, right = contig.coordinates
+            assert chromosome == chrom.id
+            assert left == interval[0]
+            assert right == interval[1]
+            assert contig.origin == f"{chromosome}:{left}-{right}"
+            assert contig.length == length
 
     def test_chromosome_segment_with_genetic_map(self):
         chr_id = "chr2"
@@ -222,11 +223,11 @@ def test_chromosome_segment_with_genetic_map(self):
         genmap = "HapMapII_GRCh38"
         contig = species.get_contig(chr_id, left=left, right=right, genetic_map=genmap)
         chrom = species.get_contig(chr_id, genetic_map=genmap)
-        assert contig.recombination_map.get_rate(0) == chrom.recombination_map.get_rate(
+        assert contig.recombination_map.get_rate(
             left
-        )
+        ) == chrom.recombination_map.get_rate(left)
         assert contig.recombination_map.get_rate(
-            contig.length - 1
+            right - 1
         ) == chrom.recombination_map.get_rate(right - 1)
 
     def test_chromosome_segment_with_mask(self):
@@ -242,7 +243,11 @@ def test_chromosome_segment_with_mask(self):
             (right - offset, right + offset),
             (right + offset, right + 2 * offset),
         ]
-        clipped_mask = [(0, offset), (offset, 2 * offset), (length - offset, length)]
+        clipped_mask = [
+            (left, left + offset),
+            (left + offset, left + 2 * offset),
+            (right - offset, right),
+        ]
         species = stdpopsim.get_species("HomSap")
         contig_inclusion = species.get_contig(
             chromosome=chr_id,
@@ -268,7 +273,7 @@ def test_chromosome_segment_with_mask(self):
     def test_generic_contig_origin(self):
         length = 42
         contig = stdpopsim.Contig.basic_contig(length=length)
-        chromosome, left, right = contig.original_coordinates
+        chromosome, left, right = contig.coordinates
         assert chromosome is None
         assert left == 0
         assert right == length
@@ -279,28 +284,39 @@ def test_add_single_site_coordinate_system(self):
         species = stdpopsim.get_species("HomSap")
         interval = [100000, 200000]
         original_sweep_coord = 100100
-        shifted_sweep_coord = original_sweep_coord - interval[0]
+        bad_sweep_coord = original_sweep_coord - interval[0]
         contig = species.get_contig(chrom, left=interval[0], right=interval[1])
         # Correctly specified coordinate system
         contig.add_single_site(
             id="sweep",
             coordinate=original_sweep_coord,
         )
-        assert contig.interval_list[1][0, 0] == shifted_sweep_coord
         # Coordinate system mismatch. Shifted coordinate falls outside of
         # `interval` so the added DFE has empty intervals: throw a warning
         with pytest.warns(UserWarning, match="No intervals remain"):
             contig.add_single_site(
                 id="sweep",
-                coordinate=shifted_sweep_coord,
+                coordinate=bad_sweep_coord,
             )
 
+    def test_original_coordinates(self):
+        species = stdpopsim.get_species("AnaPla")
+        contig = species.get_contig("2", left=10000, right=392342)
+        c, x, y = contig.coordinates
+        with pytest.warns(
+            stdpopsim.DeprecatedFeatureWarning, match="no longer shifted"
+        ):
+            oc, ox, oy = contig.original_coordinates
+        assert c == oc
+        assert x == ox
+        assert y == oy
+
     def test_add_dfe_coordinate_system(self):
         chrom = "chr2"
         species = stdpopsim.get_species("HomSap")
         contig_interval = [200000, 300000]
-        original_dfe_interval = [[190000, 210000], [290000, 310000]]
-        shifted_dfe_interval = [[0, 10000], [90000, 100000]]
+        original_dfe_interval = np.array([[190000, 210000], [290000, 310000]])
+        bad_dfe_interval = np.array([[0, 10000], [90000, 100000]])
         contig = species.get_contig(
             chrom, left=contig_interval[0], right=contig_interval[1]
         )
@@ -309,12 +325,11 @@ def test_add_dfe_coordinate_system(self):
             intervals=original_dfe_interval,
             DFE=self.example_dfe,
         )
-        np.testing.assert_array_equal(contig.interval_list[1], shifted_dfe_interval)
         # Coordinate system mismatch. Shifted DFE intervals all fall outside of
         # `contig.origin` so the added DFE is empty: throw a warning
         with pytest.warns(UserWarning, match="No intervals remain"):
             contig.add_dfe(
-                intervals=shifted_dfe_interval,
+                intervals=bad_dfe_interval,
                 DFE=self.example_dfe,
             )
 
@@ -330,10 +345,16 @@ def test_dfe_breakpoints_coordinate_system(self):
             intervals=dfe_interval,
             DFE=self.example_dfe,
         )
-        shifted_breakpoints, _ = contig.dfe_breakpoints(relative_coordinates=True)
-        original_breakpoints, _ = contig.dfe_breakpoints(relative_coordinates=False)
-        for x, y in zip(shifted_breakpoints, original_breakpoints):
-            assert x + contig_interval[0] == y
+        breakpoints, _ = contig.dfe_breakpoints()
+        for x, y in zip(breakpoints, [0, 200000, 210000, 290000, 300000]):
+            assert x == y
+
+        with pytest.warns(
+            stdpopsim.DeprecatedFeatureWarning, match="relative_coordinates argument"
+        ):
+            ob, _ = contig.dfe_breakpoints(relative_coordinates=True)
+        for x, y in zip(breakpoints, ob):
+            assert x == y
 
     def test_chromosome_segment_fails_with_length_multiplier(self):
         chrom = "chr2"
diff --git a/tests/test_slim_engine.py b/tests/test_slim_engine.py
index b31ab9422..16a24e09f 100644
--- a/tests/test_slim_engine.py
+++ b/tests/test_slim_engine.py
@@ -176,7 +176,7 @@ def test_simulate(self):
             slim_burn_in=0,
         )
         assert ts.num_samples == 10
-        assert all(tree.num_roots == 1 for tree in ts.trees())
+        assert all(tree.num_roots <= 1 for tree in ts.trees(root_threshold=2))
 
     @pytest.mark.filterwarnings("ignore::stdpopsim.SLiMScalingFactorWarning")
     def test_simulate_verbosity(self):
@@ -195,7 +195,7 @@ def test_simulate_verbosity(self):
                 verbosity=v,
             )
             assert ts.num_samples == 10
-            assert all(tree.num_roots == 1 for tree in ts.trees())
+            assert all(tree.num_roots <= 1 for tree in ts.trees(root_threshold=2))
 
     @pytest.mark.filterwarnings("ignore::stdpopsim.SLiMScalingFactorWarning")
     @pytest.mark.filterwarnings("ignore::stdpopsim.DeprecatedFeatureWarning")
@@ -464,7 +464,8 @@ def test_simulate(self):
 
     def verify_slim_sim(self, ts, num_samples):
         assert ts.num_samples == num_samples
-        assert all(tree.num_roots == 1 for tree in ts.trees())
+        # regions of missing data will have 0 roots
+        assert all(tree.num_roots <= 1 for tree in ts.trees(root_threshold=2))
         n_mut_types = count_mut_types(ts)
         assert n_mut_types[0] > 0
         assert n_mut_types[1] > 0
@@ -481,6 +482,18 @@ def test_dfe_no_demography(self):
             ts = tskit.load(f.name)
         self.verify_slim_sim(ts, num_samples=10)
 
+    @pytest.mark.filterwarnings("ignore::stdpopsim.SLiMScalingFactorWarning")
+    def test_dfe_no_interval(self):
+        with tempfile.NamedTemporaryFile(mode="w") as f:
+            cmd = (
+                f"-q -e slim --slim-scaling-factor 20 --slim-path {slim_path} "
+                f"HomSap -c chr22 -l 0.01 -o {f.name} --dfe Gamma_K17 -s 984 "
+                f"pop_0:5"
+            ).split()
+            capture_output(stdpopsim.cli.stdpopsim_main, cmd)
+            ts = tskit.load(f.name)
+        self.verify_slim_sim(ts, num_samples=10)
+
     @pytest.mark.filterwarnings("ignore::stdpopsim.SLiMScalingFactorWarning")
     def test_dfe_interval(self):
         with tempfile.NamedTemporaryFile(mode="w") as f:
@@ -544,15 +557,19 @@ def test_dfe_bed_file(self, tmp_path):
     def test_chromosomal_segment(self):
         left = 101024
         right = 200111
+        sp = "HomSap"
+        chrom = "chr22"
         with tempfile.NamedTemporaryFile(mode="w") as f:
             cmd = (
                 f"-q -e slim --slim-scaling-factor 20 --slim-path {slim_path} "
-                f"HomSap -c chr22 --left {left} --right {right} -o {f.name} "
+                f"{sp} -c {chrom} --left {left} --right {right} -o {f.name} "
                 f"--dfe Gamma_K17 -s 24 pop_0:5"
             ).split()
             capture_output(stdpopsim.cli.stdpopsim_main, cmd)
             ts = tskit.load(f.name)
-        assert ts.sequence_length == right - left
+        species = stdpopsim.get_species(sp)
+        contig = species.get_contig(chrom)
+        assert ts.sequence_length == contig.length
         self.verify_slim_sim(ts, num_samples=10)
 
     @pytest.mark.filterwarnings("ignore::stdpopsim.SLiMScalingFactorWarning")
@@ -561,22 +578,28 @@ def test_chromosomal_segment_with_dfe_interval(self):
         right = 200111
         dfe_left = left - 5000
         dfe_right = left + 90000
+        sp = "HomSap"
+        chrom = "chr21"
         with tempfile.NamedTemporaryFile(mode="w") as f:
             cmd = (
                 f"-q -e slim --slim-scaling-factor 20 --slim-path {slim_path} "
-                f"HomSap -c chr22 --left {left} --right {right} -o {f.name} "
+                f"{sp} -c {chrom} --left {left} --right {right} -o {f.name} "
                 f"--dfe Gamma_K17 -s 984 --dfe-interval {dfe_left},{dfe_right} "
                 f"pop_0:5"
             ).split()
             capture_output(stdpopsim.cli.stdpopsim_main, cmd)
             ts = tskit.load(f.name)
-        assert ts.sequence_length == right - left
+        species = stdpopsim.get_species(sp)
+        contig = species.get_contig(chrom)
+        assert ts.sequence_length == contig.length
         self.verify_slim_sim(ts, num_samples=10)
 
     @pytest.mark.filterwarnings("ignore::stdpopsim.SLiMScalingFactorWarning")
     def test_chromosomal_segment_with_dfe_bed_file(self, tmp_path):
         left = 101024
         right = 300111
+        sp = "HomSap"
+        chrom = "chr19"
         lines = [
             "\t".join(["chr22", "100000", "145000"]),
             "\t".join(["chr22", "150000", "302425"]),
@@ -588,29 +611,35 @@ def test_chromosomal_segment_with_dfe_bed_file(self, tmp_path):
         fname = tmp_path / "sim.trees"
         cmd = (
             f"-q -e slim --slim-scaling-factor 20 --slim-path {slim_path} "
-            f"HomSap -c chr22 -s 1234 --left {left} --right {right} -o {fname} "
+            f"{sp} -c {chrom} -s 1234 --left {left} --right {right} -o {fname} "
             f"--dfe Gamma_K17 -s 183 --dfe-bed-file {tmp_path / 'ex.bed'} "
             f"pop_0:5"
         ).split()
         capture_output(stdpopsim.cli.stdpopsim_main, cmd)
         ts = tskit.load(fname)
-        assert ts.sequence_length == right - left
+        species = stdpopsim.get_species(sp)
+        contig = species.get_contig(chrom)
+        assert ts.sequence_length == contig.length
         self.verify_slim_sim(ts, num_samples=10)
 
     @pytest.mark.filterwarnings("ignore::stdpopsim.SLiMScalingFactorWarning")
     def test_chromosomal_segment_with_dfe_annotation(self):
         left = 37500000
-        right = 38000000
+        right = 48000000
+        sp = "HomSap"
+        chrom = "chr18"
         with tempfile.NamedTemporaryFile(mode="w") as f:
             cmd = (
                 f"-q -e slim --slim-scaling-factor 20 --slim-path {slim_path} "
-                f"HomSap -c chr22 --left {left} --right {right} -o {f.name} "
+                f"{sp} -c {chrom} --left {left} --right {right} -o {f.name} "
                 f"--dfe Gamma_K17 -s 913 --dfe-annotation ensembl_havana_104_CDS "
                 f"pop_0:5"
             ).split()
             capture_output(stdpopsim.cli.stdpopsim_main, cmd)
             ts = tskit.load(f.name)
-        assert ts.sequence_length == right - left
+        species = stdpopsim.get_species(sp)
+        contig = species.get_contig(chrom)
+        assert ts.sequence_length == contig.length
         self.verify_slim_sim(ts, num_samples=10)
 
     # tmp_path is a pytest fixture
@@ -1224,6 +1253,7 @@ def test_slim_simulates_dfes(self):
         breaks = np.linspace(0, contig.length, n + 1)
         for j, dfe in enumerate(dfes):
             interval = [breaks[j], breaks[j + 1]]
+            print(j, interval)
             contig.add_dfe(intervals=np.array([interval], dtype="int"), DFE=dfe)
         engine = stdpopsim.get_engine("slim")
         ts = engine.simulate(
@@ -1694,7 +1724,7 @@ def test_chromosomal_segment(self):
         )
         self.verify_genomic_elements(contig, ts)
         self.verify_mutation_rates(contig, ts)
-        assert int(ts.sequence_length) == right - left
+        assert int(ts.sequence_length) == contig.length
 
     def test_dominance_coeff_list(self):
         # test that discretized h-s relationships work