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/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 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