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

Issue:
Add a way to retrieve a (possibly random) interval from a `Contig`
As discussed in #670 and #401, and started in #402.
Add get_chunk(), to get a subset of the genome.
Should eventually close #401.

I\m opening this PR early to verify that the basic approach is sound.
update ubuntu in tests
Github warns us about our tests now:
```
The ubuntu-18.04 environment is deprecated, consider switching to ubuntu-20.04(ubuntu-latest), or ubuntu-22.04 instead. For more details see https://github.com/actions/virtual-environments/issues/6002
```
... I figure we should just change to [20.04 here](https://github.com/popsim-consortium/stdpopsim/blob/main/.github/workflows/tests.yml#L20).


Discussion / hints:

# [Codecov](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/402?src=pr&el=h1) Report
> Merging [#402](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/402?src=pr&el=desc) into [master](https://codecov.io/gh/popsim-consortium/stdpopsim/commit/b3f117dbf91d578fe95023087d371a6cfee17f46?src=pr&el=desc) will **increase** coverage by `0.02%`.
> The diff coverage is `100%`.

[![Impacted file tree graph](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/402/graphs/tree.svg?width=650&token=AOkcaVfvhX&height=150&src=pr)](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/402?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master     #402      +/-   ##
==========================================
+ Coverage   99.49%   99.51%   +0.02%     
==========================================
  Files          15       15              
  Lines        1179     1243      +64     
  Branches      121      136      +15     
==========================================
+ Hits         1173     1237      +64     
  Misses          2        2              
  Partials        4        4
```


| [Impacted Files](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/402?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [stdpopsim/cli.py](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/402/diff?src=pr&el=tree#diff-c3RkcG9wc2ltL2NsaS5weQ==) | `97.93% <100%> (+0.17%)` | :arrow_up: |
| [stdpopsim/species.py](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/402/diff?src=pr&el=tree#diff-c3RkcG9wc2ltL3NwZWNpZXMucHk=) | `100% <100%> (ø)` | :arrow_up: |
| [stdpopsim/genomes.py](https://codecov.io/gh/popsim-consortium/stdpopsim/pull/402/diff?src=pr&el=tree#diff-c3RkcG9wc2ltL2dlbm9tZXMucHk=) | `100% <100%> (ø)` | :arrow_up: |

------

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

Added CLI parameters and docstrings. I decided against using `chr:from-to` syntax on the CLI.
I converted this to a draft to make it easy see what PRs are live for the next release @grahamgower - feel free to change back if you disagree.
Is it simplest to close this @grahamgower or would you like to keep it open and update to use the msprime 1.0 apis?
I think I'd prefer to keep this open for now, so I don't forget about it. As you say, it should be possible now to use the msprime 1.0 api.
hey @grahamgower -- it would be great to get this feature in to stdpopsim at this point. want to revisit it? if you're busy i'm sure someone else could pick it up
It does seem like all the pieces are now in place for this get finished up. I don't have the time/motivation to work through this right now, but I'd be happy to advise if someone else wanted to pick it up. The pr needs to be split into two:
- The ability to specify a genomic interval, so that the recombination map and annotations correspond to that interval.
- Randomly choosing a genomic interval of a given length. This is a little trickier because we now support masking of recombination maps, and we definitely want to respect those masks.
I can take this over, if that's OK with you @grahamgower.
That would be great, thanks @nspope!


Fix patch:
diff --git a/docs/tutorial.rst b/docs/tutorial.rst
index f5c0d782f..51b1705b9 100644
--- a/docs/tutorial.rst
+++ b/docs/tutorial.rst
@@ -337,8 +337,30 @@ To simulate, however, only a small chunk of contig with selection––as in a s
 
 The DFE defined by Gamma_K17 will be incorporated only in positions between 1000 and 100000.
 
+Finally, for the sake of efficiency it might be desirable to simulate a region
+of a chromosome rather than the entire chromosome. Regions are specified via
+``--left`` and ``--right``:
+
+.. code-block:: console
+
+    $ stdpopsim -e slim --slim-scaling-factor 10 HomSap \
+        -c chr22 --dfe Gamma_K17 -o foo.ts -d OutOfAfrica_2T12 \
+        --dfe-interval 200000,300000 --left 150000 --right 350000 2 4
+
+This will simulate a 200 Kb contig that uses the genetic map from coordinates
+150000 to 350000 along `chr22`; with a DFE applied between coordinates 200000
+and 300000. Note that the coordinate system used for ``--dfe-interval`` is that
+of the original chromosome: relative to the **start of the contig**, the DFE
+spans from coordinate 50000 to 150000.
+
 See also the python API to incorporate `selection <https://popsim-consortium.github.io/stdpopsim-docs/latest/tutorial.html#incorporating-selection>`__.
 
+.. warning::
+
+    Simulating a region under selection is **not** the same as simulating a
+    chromosome under selection and clipping to the region. This is because
+    selected mutations outside of the region can influence ancestry within
+    the region, due to linkage.
 
 Debugging output from SLiM
 ==========================
@@ -1177,6 +1199,26 @@ but nonetheless there is lower diversity in exons than outside of them:
     # and outside of exons it is 0.211 differences per Kb.
 
 
+For the sake of efficiency, we might want to simulate a particular region
+rather than the entire annotated chromosome. This can be done by supplying
+`left` and `right` to :meth:`Species.get_contig`. For example, to simulate the
+region of chromosome 20 spanning from 30 Mb to 40 Mb (containing roughly 300
+identified genes):
+
+.. code-block:: python
+
+    contig_region = species.get_contig("chr20", left=30e6, right=40e6)
+    contig_region.add_dfe(intervals=exon_intervals, DFE=dfe)
+
+The annotation (as well as any supplied masks) will be automatically clipped to
+the region of interest.
+
+.. warning::
+
+    Simulating a region under selection is **not** the same as simulating a
+    chromosome under selection and clipping to the region. This is because
+    selected mutations outside of the region can influence ancestry within
+    the region, due to linkage.
 
 .. _sec_tute_selective_sweep:
 
diff --git a/stdpopsim/cli.py b/stdpopsim/cli.py
index 165fda711..89706c2ac 100644
--- a/stdpopsim/cli.py
+++ b/stdpopsim/cli.py
@@ -556,9 +556,23 @@ def add_simulate_species_parser(parser, species):
         "--length-multiplier",
         default=1,
         type=float,
-        help="Simulate a sequence of length l times the named chromosome's length, "
+        help="Simulate a sequence of length l times the named chromosome's length "
         "using the named chromosome's mutation and recombination rates.",
     )
+    species_parser.add_argument(
+        "--left",
+        default=None,
+        type=float,
+        help="The leftmost coordinate (inclusive) of the contig on the named "
+        "chromosome. Defaults to 0 (the beginning of the chromosome).",
+    )
+    species_parser.add_argument(
+        "--right",
+        default=None,
+        type=float,
+        help="The rightmost coordinate (exclusive) of the contig on the named "
+        "chromosome. Defaults to the length of the chromosome.",
+    )
     species_parser.add_argument(
         "-s",
         "--seed",
@@ -602,11 +616,11 @@ def add_simulate_species_parser(parser, species):
     )
     interval_help = (
         "Specify the interval where selection (given a DFE) is simulated. "
-        "Anything outside of the interval is simulated as neutral "
+        "Anything outside of the interval is simulated as neutral. "
         "If no interval is specified, "
-        "selection is simulated across the entire region (contig). "
-        "Interval (for now) is writen as 'left,right' (separated by a comma, "
-        "with no space, for instance: --dfe-interval 1000,2000"
+        "selection is simulated across the entire contig. "
+        "The interval is written as 'start,end' (separated by a comma, "
+        "with no space, for instance: --dfe-interval 1000,2000. "
         # say something here about multiple dfes?
     )
     species_parser.add_argument(
@@ -618,7 +632,7 @@ def add_simulate_species_parser(parser, species):
     )
     bed_help = (
         "A bed file specifing the intervals where selection (given a DFE) is simulated. "
-        "Non-overlaping intervals belonging to the same chromosome are required. "
+        "Non-overlapping intervals belonging to the same chromosome are required. "
         "If no interval is specified, "
         "selection is simulated across the entire contig. "
         "See also --dfe-interval and --dfe-annotation."
@@ -709,6 +723,8 @@ def run_simulation(args):
             inclusion_mask=args.inclusion_mask,
             exclusion_mask=args.exclusion_mask,
             mutation_rate=model.mutation_rate,
+            left=args.left,
+            right=args.right,
         )
         engine = stdpopsim.get_engine(args.engine)
         logger.info(
@@ -765,9 +781,8 @@ 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
-                intervals = np.array(
-                    [[0, int(contig.recombination_map.sequence_length)]]
-                )
+                _, left, right = contig.original_coordinates
+                intervals = np.array([[left, right]])
                 intervals_summary_str = f"[{intervals[0][0]}, {intervals[0][1]})"
 
             dfe = species.get_dfe(args.dfe)
@@ -861,7 +876,9 @@ def write_simulation_summary(
             f"using mutation rate from species contig ({contig.mutation_rate})"
         )
     contig_len = contig.recombination_map.sequence_length
+    contig_origin = contig.origin
     dry_run_text += "Contig Description:\n"
+    dry_run_text += f"{indent}Contig origin: {contig_origin}\n"
     dry_run_text += f"{indent}Contig length: {contig_len}\n"
     dry_run_text += f"{indent}Mean recombination rate: {mean_recomb_rate}\n"
     dry_run_text += f"{indent}Mean mutation rate: {mut_rate}\n"
diff --git a/stdpopsim/genetic_maps.py b/stdpopsim/genetic_maps.py
index 2bb40f367..102549839 100644
--- a/stdpopsim/genetic_maps.py
+++ b/stdpopsim/genetic_maps.py
@@ -122,8 +122,10 @@ def get_chromosome_map(self, id):
             rates = np.append(recomb_map.rate, 0)
             recomb_map = msprime.RateMap(position=positions, rate=rates)
         elif map_length > chrom.length:
+            # TODO: consider making this an error, and deprecating genetic maps
+            # that do not match the assembly.
             warnings.warn(
                 f"Recombination map has length {map_length}, which is longer than"
-                f" chromosome length {chrom.length}. The former will be used."
+                f" chromosome length {chrom.length}. The latter will be used."
             )
         return recomb_map
diff --git a/stdpopsim/genomes.py b/stdpopsim/genomes.py
index fc14a1708..4780d2089 100644
--- a/stdpopsim/genomes.py
+++ b/stdpopsim/genomes.py
@@ -189,6 +189,9 @@ 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 gene_conversion_rate: The rate of gene conversion initiation per base per
@@ -207,12 +210,17 @@ class Contig:
     :ivar exclude: If True, ``mask_intervals`` specify regions to exclude. If False,
         ``mask_intervals`` specify regions in keep.
     :vartype exclude: bool
-    :ivar dfe_list: a list of :class:`.DFE` objects.
+    :ivar dfe_list: A list of :class:`.DFE` objects.
         By default, the only DFE is completely neutral.
     :vartype dfe_list: list
     :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,
+        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
 
     .. note::
         To run stdpopsim simulations with alternative, user-specified mutation
@@ -238,10 +246,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)
 
     def __attrs_post_init__(self):
+        if self.original_coordinates is None:
+            self.original_coordinates = (None, 0, int(self.length))
+        _, left, right = self.original_coordinates
         self.add_dfe(
-            np.array([[0, int(self.length)]]),
+            [[left, right]],
             stdpopsim.dfe.neutral_dfe(),
         )
 
@@ -276,6 +288,8 @@ def species_contig(
         gene_conversion_length=None,
         inclusion_mask=None,
         exclusion_mask=None,
+        left=None,
+        right=None,
     ):
         """
         Build a Contig for a species.
@@ -309,6 +323,10 @@ def species_contig(
                     "gene conversion rate"
                 )
         if chromosome is None:
+            if left is not None or right is not None:
+                raise ValueError(
+                    "Cannot use left or right coordinates with generic contig"
+                )
             if genetic_map is not None:
                 raise ValueError("Cannot use genetic map with generic contig")
             if length_multiplier != 1:
@@ -352,12 +370,44 @@ def species_contig(
                 raise ValueError("Cannot specify sequence length for named contig")
             if inclusion_mask is not None and exclusion_mask is not None:
                 raise ValueError("Cannot specify both inclusion and exclusion masks")
+            if length_multiplier != 1 and (left is not None or right is not None):
+                raise ValueError(
+                    "Cannot use length multiplier when specifying left or "
+                    "right coordinates of the contig on a named chromosome."
+                )
+
             chrom = species.genome.get_chromosome(chromosome)
+            if left is None:
+                left = 0
+            else:
+                left = round(left)
+                if not 0 <= left < chrom.length:
+                    raise ValueError(
+                        f"Left coordinate {left} falls outside chromosome "
+                        f"{chrom.id} with length {chrom.length}."
+                    )
+            if right is None:
+                right = round(chrom.length * length_multiplier)
+            else:
+                right = round(right)
+                if not left < right <= chrom.length:
+                    raise ValueError(
+                        f"Right coordinate {right} falls outside chromosome "
+                        f"{chrom.id} with length {chrom.length}."
+                    )
+
             if genetic_map is None:
-                logger.debug(f"Making flat chromosome {length_multiplier} * {chrom.id}")
                 gm = None
+                if length_multiplier != 1:
+                    logger.debug(
+                        f"Making flat chromosome {length_multiplier} * {chrom.id}"
+                    )
+                else:
+                    logger.debug(
+                        f"Making flat contig of length {right - left} from {chrom.id}"
+                    )
                 recomb_map = msprime.RateMap.uniform(
-                    round(chrom.length * length_multiplier), chrom.recombination_rate
+                    right - left, chrom.recombination_rate
                 )
             else:
                 if length_multiplier != 1:
@@ -365,6 +415,7 @@ def species_contig(
                 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)
 
             inclusion_intervals = None
             exclusion_intervals = None
@@ -377,6 +428,9 @@ def species_contig(
                     )
                 else:
                     inclusion_intervals = inclusion_mask
+                inclusion_intervals = stdpopsim.utils.clip_and_shift_intervals(
+                    inclusion_intervals, left, right
+                )
             if exclusion_mask is not None:
                 if length_multiplier != 1:
                     raise ValueError("Cannot use length multiplier with mask")
@@ -386,6 +440,9 @@ def species_contig(
                     )
                 else:
                     exclusion_intervals = exclusion_mask
+                exclusion_intervals = stdpopsim.utils.clip_and_shift_intervals(
+                    exclusion_intervals, left, right
+                )
 
             if mutation_rate is None:
                 mutation_rate = chrom.mutation_rate
@@ -402,6 +459,7 @@ def species_contig(
                 gene_conversion_length=gene_conversion_length,
                 inclusion_mask=inclusion_intervals,
                 exclusion_mask=exclusion_intervals,
+                original_coordinates=(chromosome, left, right),
             )
 
         return contig
@@ -410,7 +468,19 @@ def species_contig(
     def length(self):
         return self.recombination_map.sequence_length
 
-    def dfe_breakpoints(self):
+    @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
+        if chromosome is None:
+            return None
+        else:
+            return f"{chromosome}:{left}-{right}"
+
+    def dfe_breakpoints(self, relative_coordinates=True):
         """
         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,
@@ -433,6 +503,10 @@ def dfe_breakpoints(self):
         ``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).
         """
         breaks = np.unique(
@@ -441,6 +515,9 @@ def dfe_breakpoints(self):
         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):
@@ -468,13 +545,19 @@ def add_dfe(self, intervals, DFE):
         contig.add_dfe(a1, dfe1)
         contig.add_dfe(a2, dfe2)
         ```
-        then ``dfe1`` applies to the region from 0 to 50 and ``dfe2`` applies to the
-        region from 50 to 120.
+        then ``dfe1`` applies to the region from 0 to 50 and ``dfe2`` applies
+        to the region from 50 to 120.
+
+        Any of the ``intervals`` that fall outside of the contig will be
+        clipped to the contig boundaries. If no ``intervals`` overlap the
+        contig, an "empty" DFE will be added with a warning.
 
         :param array intervals: A valid set of intervals.
         :param DFE dfe: A DFE object.
         """
-        stdpopsim.utils._check_intervals_validity(intervals)
+        _, 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)
         for j, ints in enumerate(self.interval_list):
             self.interval_list[j] = stdpopsim.utils.mask_intervals(ints, intervals)
         self.dfe_list.append(DFE)
@@ -612,8 +695,8 @@ def __str__(self):
         s = (
             "Contig(length={:.2G}, recombination_rate={:.2G}, "
             "mutation_rate={:.2G}, gene_conversion_rate={}, "
-            "gene_conversion_length{}, genetic_map={}, dfe_list={}, "
-            "interval_list={})"
+            "gene_conversion_length{}, genetic_map={}, origin={}, "
+            "dfe_list={}, interval_list={})"
         ).format(
             self.length,
             self.recombination_map.mean_rate,
@@ -621,6 +704,7 @@ def __str__(self):
             self.gene_conversion_rate,
             self.gene_conversion_length,
             gmap,
+            self.origin,
             self.dfe_list,
             self.interval_list,
         )
diff --git a/stdpopsim/species.py b/stdpopsim/species.py
index 716b0ac1d..a1ce682da 100644
--- a/stdpopsim/species.py
+++ b/stdpopsim/species.py
@@ -175,6 +175,8 @@ def get_contig(
         gene_conversion_length=None,
         inclusion_mask=None,
         exclusion_mask=None,
+        left=None,
+        right=None,
     ):
         """
         Returns a :class:`.Contig` instance describing a section of genome that
@@ -223,6 +225,12 @@ def get_contig(
             length of genome sequence for this contig. For a generic contig, mutation
             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.
+        :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.
         :rtype: :class:`.Contig`
         :return: A :class:`.Contig` describing the section of the genome.
         """
@@ -238,6 +246,8 @@ def get_contig(
             gene_conversion_length=gene_conversion_length,
             inclusion_mask=inclusion_mask,
             exclusion_mask=exclusion_mask,
+            left=left,
+            right=right,
         )
 
     def get_demographic_model(self, id):
diff --git a/stdpopsim/utils.py b/stdpopsim/utils.py
index f4cb2ea67..f85ad8471 100644
--- a/stdpopsim/utils.py
+++ b/stdpopsim/utils.py
@@ -9,6 +9,7 @@
 import tarfile
 import contextlib
 import numpy as np
+import warnings
 
 
 def is_valid_demographic_model_id(model_id):
@@ -249,3 +250,21 @@ def mask_intervals(intervals, mask):
             left = max(next_left, next_right)
     out = np.array(out, dtype=intervals.dtype).reshape((len(out), 2))
     return out
+
+
+def clip_and_shift_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.
+    """
+    assert 0 <= left < right
+    intervals = np.array(intervals).astype(int)
+    _check_intervals_validity(intervals)
+    out_of_bounds = np.logical_or(intervals[:, 1] <= left, intervals[:, 0] >= 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
diff --git a/stdpopsim/warning_categories.py b/stdpopsim/warning_categories.py
index 785c45b74..a4d7365a7 100644
--- a/stdpopsim/warning_categories.py
+++ b/stdpopsim/warning_categories.py
@@ -22,3 +22,7 @@ class SLiMOddSampleWarning(UserWarning):
 
 class UnspecifiedSLiMWarning(UserWarning):
     pass
+
+
+class EmptyIntervalsWarning(UserWarning):
+    pass


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 98aeb688d..2fdd1d6a6 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -17,7 +17,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        os: [ubuntu-18.04, macos-10.15, windows-latest]
+        os: [ubuntu-20.04, macos-10.15, windows-latest]
         python: [3.7, "3.10"]
     env:
       CONDA_ENV_NAME: stdpopsim
diff --git a/tests/test_genomes.py b/tests/test_genomes.py
index 862f8aae3..7d17b0d8b 100644
--- a/tests/test_genomes.py
+++ b/tests/test_genomes.py
@@ -110,7 +110,6 @@ def test_add_dfe(self):
                 assert np.all(a1 == a2)
             self.verify_dfe_breakpoints(contig)
 
-    @pytest.mark.skip(reason="TODO allow more flexible inputs")
     def test_add_dfe_interval_formats(self):
         L = 50818
         for intervals in (
@@ -121,6 +120,8 @@ def test_add_dfe_interval_formats(self):
             contig = stdpopsim.Contig.basic_contig(length=L)
             contig.add_dfe(intervals=intervals, DFE=self.example_dfe)
             np.testing.assert_array_equal(intervals, contig.interval_list[1])
+        with pytest.raises(ValueError, match="must be a numpy array"):
+            stdpopsim.utils._check_intervals_validity([[0, 1]])
 
     def test_is_neutral(self):
         for neutral in (True, False):
@@ -156,6 +157,170 @@ def test_is_neutral(self):
                 # exponential with mean zero doesn't count as neutral!
                 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}"
+
+    def test_chromosome_segment_with_genetic_map(self):
+        chr_id = "chr2"
+        left = 1000001
+        right = 3000000
+        species = stdpopsim.get_species("HomSap")
+        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(
+            left
+        )
+        assert contig.recombination_map.get_rate(
+            contig.length - 1
+        ) == chrom.recombination_map.get_rate(right - 1)
+
+    def test_chromosome_segment_with_mask(self):
+        chr_id = "chr2"
+        left = 1000001
+        right = 3000000
+        length = right - left
+        offset = int(length * 0.1)
+        mask = [
+            (left - 2 * offset, left - offset),
+            (left - offset, left + offset),
+            (left + offset, left + 2 * offset),
+            (right - offset, right + offset),
+            (right + offset, right + 2 * offset),
+        ]
+        clipped_mask = [(0, offset), (offset, 2 * offset), (length - offset, length)]
+        species = stdpopsim.get_species("HomSap")
+        contig_inclusion = species.get_contig(
+            chromosome=chr_id,
+            left=left,
+            right=right,
+            inclusion_mask=mask,
+        )
+        contig_exclusion = species.get_contig(
+            chromosome=chr_id,
+            left=left,
+            right=right,
+            exclusion_mask=mask,
+        )
+        np.testing.assert_array_equal(
+            contig_inclusion.inclusion_mask,
+            clipped_mask,
+        )
+        np.testing.assert_array_equal(
+            contig_exclusion.exclusion_mask,
+            clipped_mask,
+        )
+
+    def test_generic_contig_origin(self):
+        length = 42
+        contig = stdpopsim.Contig.basic_contig(length=length)
+        chromosome, left, right = contig.original_coordinates
+        assert chromosome is None
+        assert left == 0
+        assert right == length
+        assert contig.origin is None
+
+    def test_add_single_site_coordinate_system(self):
+        chrom = "chr2"
+        species = stdpopsim.get_species("HomSap")
+        interval = [100000, 200000]
+        original_sweep_coord = 100100
+        shifted_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,
+            )
+
+    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]]
+        contig = species.get_contig(
+            chrom, left=contig_interval[0], right=contig_interval[1]
+        )
+        # Correctly specified coordinate system
+        contig.add_dfe(
+            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,
+                DFE=self.example_dfe,
+            )
+
+    def test_dfe_breakpoints_coordinate_system(self):
+        chrom = "chr2"
+        species = stdpopsim.get_species("HomSap")
+        contig_interval = [200000, 300000]
+        dfe_interval = [[190000, 210000], [290000, 310000]]
+        contig = species.get_contig(
+            chrom, left=contig_interval[0], right=contig_interval[1]
+        )
+        contig.add_dfe(
+            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
+
+    def test_chromosome_segment_fails_with_length_multiplier(self):
+        chrom = "chr2"
+        species = stdpopsim.get_species("HomSap")
+        with pytest.raises(ValueError, match="specifying left or right"):
+            species.get_contig(
+                chrom,
+                left=200000,
+                right=300000,
+                length_multiplier=8.2,
+            )
+
+    def test_no_chromosome_segment_with_generic_contig(self):
+        species = stdpopsim.get_species("HomSap")
+        interval = [200000, 300000]
+        with pytest.raises(ValueError, match="coordinates with generic contig"):
+            species.get_contig(
+                left=interval[0],
+                right=interval[1],
+            )
+
+    def test_chromosome_segment_exceeds_length(self):
+        chr_id = "chr22"
+        species = stdpopsim.get_species("HomSap")
+        for interval in [[50e6, 80e6], [80e6, 100e6]]:
+            with pytest.raises(ValueError, match="falls outside chromosome"):
+                species.get_contig(
+                    chromosome=chr_id,
+                    left=interval[0],
+                    right=interval[1],
+                )
+
 
 class TestGeneConversion(object):
     species = stdpopsim.get_species("EscCol")
diff --git a/tests/test_masking.py b/tests/test_masking.py
index 9476b1991..3453ff03d 100644
--- a/tests/test_masking.py
+++ b/tests/test_masking.py
@@ -1,7 +1,6 @@
 """
 Tests for the genetic maps management.
 """
-import os
 import sys
 
 import numpy as np
@@ -14,18 +13,19 @@
 
 
 class TestMasking:
-    def test_load_intervals(self):
+    @pytest.mark.usefixtures("tmp_path")
+    def test_load_intervals(self, tmp_path):
         intervals_in = {
             "chr1": [(10, 10000)],
             "chr22": [(100, 1000), (2000, 5000), (6000, 9000)],
         }
-        with open("temp_file.bed", "w+") as fout:
+        bedfile = str(tmp_path / "temp_file.bed")
+        with open(bedfile, "w+") as fout:
             for c, i in intervals_in.items():
                 for (l, r) in i:
                     fout.write(f"{c}\t{l}\t{r}\n")
-        intervals_chr1 = stdpopsim.utils.read_bed("temp_file.bed", "chr1")
-        intervals_chr22 = stdpopsim.utils.read_bed("temp_file.bed", "chr22")
-        os.remove("temp_file.bed")
+        intervals_chr1 = stdpopsim.utils.read_bed(bedfile, "chr1")
+        intervals_chr22 = stdpopsim.utils.read_bed(bedfile, "chr22")
         for interval in intervals_chr1:
             assert tuple(interval) in intervals_in["chr1"]
         for interval in intervals_chr22:
@@ -35,13 +35,13 @@ def test_load_intervals(self):
 
     def test_length_interval_invalid(self):
         species = stdpopsim.get_species("HomSap")
-        with pytest.raises(ValueError):
-            contig = species.get_contig(
+        with pytest.raises(ValueError, match="Cannot use length multiplier"):
+            species.get_contig(
                 "chr22", length_multiplier=0.1, inclusion_mask="test.bed"
             )
-            print(contig.recombination_map.get_sequence_length())
 
-    def test_read_bed(self):
+    @pytest.mark.usefixtures("tmp_path")
+    def test_read_bed(self, tmp_path):
         def intervals_to_keep_test_func(mask_fpath, chrom):
             intervals = []
             with open(mask_fpath, "r") as fin:
@@ -52,13 +52,13 @@ def intervals_to_keep_test_func(mask_fpath, chrom):
             return intervals
 
         intervals_in = {"chr1": [(0, 10), (100, 1000), (2000, 5000), (6000, 10000)]}
-        with open("temp_file.bed", "w+") as fout:
+        bedfile = str(tmp_path / "temp_file.bed")
+        with open(bedfile, "w+") as fout:
             for c, i in intervals_in.items():
                 for (l, r) in i:
                     fout.write(f"{c}\t{l}\t{r}\n")
-        intervals_test = intervals_to_keep_test_func("temp_file.bed", "chr1")
-        intervals_utils = stdpopsim.utils.read_bed("temp_file.bed", "chr1")
-        os.remove("temp_file.bed")
+        intervals_test = intervals_to_keep_test_func(bedfile, "chr1")
+        intervals_utils = stdpopsim.utils.read_bed(bedfile, "chr1")
         for i1, i2 in zip(intervals_utils, intervals_test):
             assert i1[0] == i2[0]
             assert i1[1] == i2[1]
@@ -81,20 +81,21 @@ def test_multiple_masks(self):
                 "chr22", exclusion_mask=[(0, 16.5e6)], inclusion_mask=[(16.5e6, 50e6)]
             )
 
-    def test_read_masks_from_bed(self):
+    @pytest.mark.usefixtures("tmp_path")
+    def test_read_masks_from_bed(self, tmp_path):
         intervals_in = {"chr1": [(0, 10), (100, 1000), (2000, 5000), (6000, 10000)]}
-        with open("temp_file.bed", "w+") as fout:
+        bedfile = str(tmp_path / "temp_file.bed")
+        with open(bedfile, "w+") as fout:
             for c, i in intervals_in.items():
                 for (l, r) in i:
                     fout.write(f"{c}\t{l}\t{r}\n")
         species = stdpopsim.get_species("HomSap")
-        contig = species.get_contig("chr1", inclusion_mask="temp_file.bed")
+        contig = species.get_contig("chr1", inclusion_mask=bedfile)
         assert contig.exclusion_mask is None
         assert len(contig.inclusion_mask) == 4
-        contig = species.get_contig("chr1", exclusion_mask="temp_file.bed")
+        contig = species.get_contig("chr1", exclusion_mask=bedfile)
         assert contig.inclusion_mask is None
         assert contig.exclusion_mask[1][0] == 100
-        os.remove("temp_file.bed")
 
     def test_mask_tree_sequence(self):
         intervals = np.array([[0, 10], [100, 200], [500, 1000]])
diff --git a/tests/test_slim_engine.py b/tests/test_slim_engine.py
index c879bb19b..6b9e73ec1 100644
--- a/tests/test_slim_engine.py
+++ b/tests/test_slim_engine.py
@@ -371,6 +371,76 @@ def test_dfe_bed_file(self, tmp_path):
         ts = tskit.load(fname)
         self.verify_slim_sim(ts, num_samples=10)
 
+    @pytest.mark.filterwarnings("ignore::stdpopsim.SLiMScalingFactorWarning")
+    def test_chromosomal_segment(self):
+        left = 101024
+        right = 200111
+        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"--dfe Gamma_K17 -s 24 10"
+            ).split()
+            capture_output(stdpopsim.cli.stdpopsim_main, cmd)
+            ts = tskit.load(f.name)
+        assert ts.sequence_length == right - left
+        self.verify_slim_sim(ts, num_samples=10)
+
+    @pytest.mark.filterwarnings("ignore::stdpopsim.SLiMScalingFactorWarning")
+    def test_chromosomal_segment_with_dfe_interval(self):
+        left = 101024
+        right = 200111
+        dfe_left = left - 5000
+        dfe_right = left + 90000
+        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"--dfe Gamma_K17 -s 984 --dfe-interval {dfe_left},{dfe_right} 10"
+            ).split()
+            capture_output(stdpopsim.cli.stdpopsim_main, cmd)
+            ts = tskit.load(f.name)
+        assert ts.sequence_length == right - left
+        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
+        lines = [
+            "\t".join(["chr22", "100000", "145000"]),
+            "\t".join(["chr22", "150000", "302425"]),
+        ]
+        bedfile = open(tmp_path / "ex.bed", "w")
+        for lin in lines:
+            bedfile.write(lin + "\n")
+        bedfile.close()
+        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"--dfe Gamma_K17 -s 183 --dfe-bed-file {tmp_path / 'ex.bed'} 10"
+        ).split()
+        capture_output(stdpopsim.cli.stdpopsim_main, cmd)
+        ts = tskit.load(fname)
+        assert ts.sequence_length == right - left
+        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
+        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"--dfe Gamma_K17 -s 913 --dfe-annotation ensembl_havana_104_CDS 10"
+            ).split()
+            capture_output(stdpopsim.cli.stdpopsim_main, cmd)
+            ts = tskit.load(f.name)
+        assert ts.sequence_length == right - left
+        self.verify_slim_sim(ts, num_samples=10)
+
     # tmp_path is a pytest fixture
     @pytest.mark.filterwarnings("ignore::stdpopsim.SLiMScalingFactorWarning")
     def test_errors(self, tmp_path):
@@ -780,9 +850,13 @@ def test_engine_available(self):
         assert IS_WINDOWS == ("slim" not in all_engines)
 
 
-def get_test_contig(spp="HomSap", chrom="chr22", length_multiplier=0.001):
+def get_test_contig(
+    spp="HomSap", chrom="chr22", length_multiplier=0.001, left=None, right=None
+):
     species = stdpopsim.get_species(spp)
-    contig = species.get_contig(chrom, length_multiplier=length_multiplier)
+    contig = species.get_contig(
+        chrom, length_multiplier=length_multiplier, left=left, right=right
+    )
     return contig
 
 
@@ -848,11 +922,10 @@ def verify_recombination_map(self, contig, ts):
         assert np.all(rm.rate[i - 1] > 0)
 
     @pytest.mark.parametrize("Q", [1, 12])
-    @pytest.mark.filterwarnings("ignore:Recombination map has length:UserWarning")
     def test_chr1(self, Q):
         engine = stdpopsim.get_engine("slim")
         species = stdpopsim.get_species("HomSap")
-        contig = species.get_contig("chr1", genetic_map="HapMapII_GRCh37")
+        contig = species.get_contig("chr1", genetic_map="HapMapII_GRCh38")
         model = stdpopsim.PiecewiseConstantSize(100)
         samples = model.get_samples(10)
         ts = engine.simulate(
@@ -1329,6 +1402,48 @@ def test_neutral_dfe_slim_proportions(self):
             [contig.mutation_rate * 0.1 * Q],
         )
 
+    def test_warn_when_dfe_intervals_outside_contig(self):
+        contig = get_test_contig()
+        dfe = stdpopsim.DFE(
+            id="test",
+            description="non-neutral",
+            long_description="",
+            mutation_types=[
+                stdpopsim.MutationType(),
+                stdpopsim.MutationType(
+                    distribution_type="lp", distribution_args=[0.01, 0.2]
+                ),
+            ],
+            proportions=[0.5, 0.5],
+        )
+        with pytest.warns(UserWarning, match="No intervals remain"):
+            contig.add_dfe(np.array([[1e6, 1.1e6]], dtype="int"), dfe)
+        assert contig.interval_list[1].shape[0] == 0
+
+    def test_chromosomal_segment(self):
+        left = 100101
+        right = 201024
+        contig = get_test_contig(length_multiplier=1, left=left, right=right)
+        L = contig.length
+        example_dfes = self.get_example_dfes()
+        for j, dfe in enumerate(example_dfes):
+            contig.add_dfe(
+                left
+                + np.array([[10 * (j + 1), 100 * (10 - j)], [L / 2, L]], dtype="int"),
+                dfe,
+            )
+        assert len(contig.dfe_list) == 1 + len(example_dfes)
+        engine = stdpopsim.get_engine("slim")
+        ts = engine.simulate(
+            demographic_model=self.model,
+            contig=contig,
+            samples=self.samples,
+            verbosity=3,  # to get metadata output
+        )
+        self.verify_genomic_elements(contig, ts)
+        self.verify_mutation_rates(contig, ts)
+        assert int(ts.sequence_length) == right - left
+
 
 @pytest.mark.skipif(IS_WINDOWS, reason="SLiM not available on windows")
 class TestLogfile(PiecewiseConstantSizeMixin):
diff --git a/tests/test_species.py b/tests/test_species.py
index 202457cc9..51d857b65 100644
--- a/tests/test_species.py
+++ b/tests/test_species.py
@@ -303,6 +303,7 @@ def test_contig_options(self):
         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)
@@ -355,3 +356,4 @@ def test_generic_contig(self):
         else:
             assert contig.gene_conversion_rate == np.average(gcs, weights=Ls)
             assert contig.gene_conversion_length == np.average(gcls, weights=Ls)
+        # TODO: origin