Minor things

.github/workflows/tests.yml

@@ -84,7 +84,7 @@ jobs:
         uses: actions/cache@2c8a9bd7457de244a408f35966fab2fb45fda9c8  # v6.0.0
         with:
           path: D:\a\pyslim\pyslim\SLiM
-          key: ${{runner.os}}-${{matrix.sys}}-${{matrix.env}}-v1-key
+          key: ${{runner.os}}-${{matrix.sys}}-${{matrix.env}}-v3-key
 
       - name: Build SLiM (Windows)
         if: matrix.os == 'windows-latest' && steps.cache-slim.outputs.cache-hit != 'true'

CHANGELOG.rst

@@ -47,10 +47,14 @@ https://tskit.dev/pyslim/docs/latest/previous_versions.html
   (`pyslim.INDIVIDUAL_FLAG_MIGRATED`) is now recorded, as `pyslim.INDIVIDUAL_MIGRATED`,
   in `individual.flags` (rather than `individual.metadata['flags']`).
 
+- The default sex ratio for populations is now 0.5 instead of 0.0.
+  (:issue:`339`, :user:`petrelharp`)
+
 **Bug fixes:**
 
 - `pyslim.annotate` now has a `num_chromosomes` argument. Previously it could not be
   easily used to annotate multichromosome simulations with more than 8 chromosomes.
+  (It also now has a `num_traits` argument.)
 
 - In some previous versions, converting files produced by a yet-older version of SLiM
   to the previously-current file version dropped some information from metadata:

docs/previous_versions.md

@@ -72,7 +72,7 @@ replaces the previous ``type`` argument to {class}`msprime.SLiMMutationModel`.
 to pull this information out of top-level metadata using the SLiM ID as a key.
 In brief, if `mut` is a mutation, then you should replace
 `mut.metadata["mutation_list"][j]`
-with `mut_metadata[mut.metadata["derived_states"][j]]`,
+with `mut_metadata[mut.metadata["slim_ids"][j]]`,
 where `mut_metadata` is the output of {func}`.mutation_metadata`.
 For instance, where before you might have done:
 ```python
@@ -85,7 +85,7 @@ Now, you would do:
 ```{code-cell}
 mut_metadata = pyslim.mutation_metadata(ts)
 mut = ts.mutation(0)
-for k in mut.metadata["derived_states"]:
+for k in mut.metadata["slim_ids"]:
     md = mut_metadata[k]
     print(f"SLiM ID: {k}")
     print(f"Metadata: {md}")

docs/traits.md

@@ -192,16 +192,16 @@ mut
 util.pp(mut)
 ```
 We can see which SLiM mutation(s) this mutation represents
-by looking up those SLiM mutation IDs in the mutation's "derived state",
-`mut.metadata["derived_states"]`.
+by looking up those SLiM mutation IDs in the mutation's metadata:
+`mut.metadata["slim_ids"]`.
 (As noted [previously](sec_tutorial_mutation_metadata),
 the same information is stored in `mut.derived_state`,
 but we recommend pulling this information out of metadata.)
 Then, we can find information about those SLiM mutations
 in the top-level mutation metadata (here, `mut_metadata`,
 obtained using {func}`.mutation_metadata`).
 ```{code-cell}
-mut_metadata[mut.metadata["derived_states"][0]]
+mut_metadata[mut.metadata["slim_ids"][0]]
 ```
 So, we can use the `effect_size` and `dominance`
 to calculate genetic effects (we won't need `hemizygous dominance`,

docs/tutorial.md

@@ -437,9 +437,9 @@ to the SLiM mutations with {func}`generate_nucleotides`,
 and (2) move those nucleotides over into the "ancestral state"
 and "derived state" slots of the tree sequence with {func}`convert_alleles`.
 If all your mutations in SLiM were nucleotide mutations, you only need to do (2).
-And, beware that (2) is an irreversible step: if you write the tree sequence
-produced by {func}`convert_alleles` to a file, you can't load that file into SLiM
-without putting those back somehow. So, to do this we'll do:
+Since the information about SLiM mutation IDs is also stored in mutation metadata,
+this entails no loss of information, and the tree sequence can still be read in by SLiM.
+So, to do this we'll do:
 
 ```{code-cell}
 nts = pyslim.generate_nucleotides(ts)
@@ -475,38 +475,46 @@ with open("example_sim2.vcf", "w") as vcffile:
 ## Mutation metadata
 
 Because of mutation stacking (see the SLiM manual),
-each "tskit mutation" can represent a superposition of more than one
-"SLiM mutation".
-This is recorded by setting the derived state of the tskit mutation
-to a comma-separated string of SLiM mutation IDs
-(or the empty string, to denote "no mutations"),
-and the `"derived_states"` entry of the tskit mutation's metadata to this same list.
-(This is redundant, but there are good reasons for it.)
+each "tskit mutation" can represent a superposition of more than one "SLiM mutation".
+This is recorded by recording the SLiM IDs of these mutations
+in the `"slim_ids"` entry of the tskit mutation's metadata.
+(Also, SLiM writes these to the derived state, as a comma-separated string.)
 So, each SLiM mutation can appear in more than one tskit mutation,
 and the [metadata](sec_metadata) about these SLiM mutations is stored in top-level metadata,
 rather than along with the tskit mutations.
 [](sec_tutorial_selected_mutations) has a more in-depth example, but here is a quick overview.
 To print out the information about each SLiM mutation "carried"
 by a given tskit mutation, whose SLiM IDs are stored in the mutation's metadata
-under `"derived_states"`, we'd do:
+under `"slim_ids"`, we'd do:
 ```{code-cell}
 :tags: ["remove-output"]
 mut_metadata = pyslim.mutation_metadata(ts)
 
 mut = ts.mutation(0)
-for x in mut.metadata["derived_states"]:
+for x in mut.metadata["slim_ids"]:
     print(f"SLiM mutation {x}:")
     print(mut_metadata[x])
 ```
 ```{code-cell}
 :tags: ["remove-input"]
-for x in mut.metadata["derived_states"]:
+for x in mut.metadata["slim_ids"]:
     print(f"SLiM mutation {x}:")
     util.pp(mut_metadata[x])
 ```
 See [](sec_tutorial_selected_mutations) for an example where a tskit mutation
 carries more than one (stacked) SLiM mutation.
 
+You may have noticed that the same information goes both into the derived state entry
+of a tskit mutation (e.g., `mut.derived_state`) and into the metadata of that
+same mutation (e.g., `mut.metadata["slim_ids"]`). The reason for this
+redundancy is that tskit uses the `derived_state` entries (and the
+`ancestral_state` entries of sites) to do VCF output and determine allelic
+identity, and so it's helpful to be able to rewrite these
+(as does {func}`.convert_alleles`; see [](sec_output)).
+In fact, SLiM does not use the `derived_state` entry of tskit mutations
+when it reads in a tree sequence (nor the `ancestral_state` entry of sites);
+the definitive source is in metadata.
+
 
 (sec_extracting_individuals)=
 
@@ -1093,7 +1101,7 @@ some number of SLiM mutations, whose SLiM IDs are stored in the mutation's `meta
 For instance, here's which SLiM mutation(s) the first mutation
 in the tree sequence represents:
 ```{code-cell}
-ds = ts.mutation(0).metadata["derived_states"]
+ds = ts.mutation(0).metadata["slim_ids"]
 print(f"SLiM IDs: {ds}")
 ```
 To see the information about these, we pull their information out
@@ -1173,7 +1181,7 @@ Now, mutations have a ``nucleotide`` property in metadata that is not ``-1``:
 :tags: ["remove-output"]
 mut_metadata = pyslim.mutation_metadata(ts)
 m = ts.mutation(0)
-md = [mut_metadata[k] for k in m.metadata["derived_states"]]
+md = [mut_metadata[k] for k in m.metadata["slim_ids"]]
 print(m)
 for x in md:
     print(x)
@@ -1193,7 +1201,7 @@ by indexing the {data}`.NUCLEOTIDES` object:
 for k in range(3):
     m = ts.mutation(k)
     print(f"Mutation {k}: position {ts.site(m.site).position}, time {m.time}")
-    for sid in m.metadata["derived_states"]:
+    for sid in m.metadata["slim_ids"]:
         md = mut_metadata[sid]
         print(f"  nucleotide: {pyslim.NUCLEOTIDES[md['nucleotide']]}")
 ```
@@ -1245,7 +1253,7 @@ Here's the first mutation:
 :tags: ["remove-output"]
 mut_metadata = pyslim.mutation_metadata(ts)
 m = ts.mutation(0)
-md = [mut_metadata[k] for k in m.metadata["derived_states"]]
+md = [mut_metadata[k] for k in m.metadata["slim_ids"]]
 print(m)
 for x in md:
     print(x)
@@ -1273,7 +1281,7 @@ and we can pull up information about that with the `ts.site( )` method:
 s = ts.site(m.site)
 md = [
     mut_metadata[k] for m in s.mutations
-                         for k in m.metadata["derived_states"]
+                    for k in m.metadata["slim_ids"]
 ]
 print(s)
 for x in md:
@@ -1289,7 +1297,7 @@ for x in md:
 
 This mutation occurred at the position along the genome shown in `site.position`,
 which previously had no mutations (since `site.ancestral_state` is the empty string, `''`)
-and was given the SLiM mutation ID shown in `m.metadata["derived_states"]`.
+and was given the SLiM mutation ID shown in `m.metadata["slim_ids"]`.
 The metadata (with `x` the mutation ID, `mut_metadata[x]`, a dict) tells us
 the mutation's selection coefficient and which population and at what SLiM time it occurred.
 This is not a nucleotide model, so the nucleotide entry is `-1`.
@@ -1313,8 +1321,8 @@ for m in ts.mutations():
      break
 
 pm = ts.mutation(m.parent)
-md = [mut_metadata[k] for k in m.metadata["derived_states"]]
-pmd = [mut_metadata[k] for k in pm.metadata["derived_states"]]
+md = [mut_metadata[k] for k in m.metadata["slim_ids"]]
+pmd = [mut_metadata[k] for k in pm.metadata["slim_ids"]]
 
 print(m)
 for x in md:
@@ -1407,7 +1415,7 @@ mut_type = np.zeros(ts.num_sites)
 for j, s in enumerate(ts.sites()):
   mt = []
   for m in s.mutations:
-     for sid in m.metadata["derived_states"]:
+     for sid in m.metadata["slim_ids"]:
         md = mut_metadata[sid]
         mt.append(md["mutation_type"])
   if len(set(mt)) > 1:
@@ -1444,7 +1452,7 @@ Finally, let's pull out information on the allele with the largest selection coe
 :tags: ["remove-output"]
 sel_coeffs = np.array([
         sum(mut_metadata[k]["per_trait"][0]["effect_size"]
-            for k in m.metadata["derived_states"])
+            for k in m.metadata["slim_ids"])
         for m in ts.mutations()
 ])
 which_max = np.argmax(sel_coeffs)

docs/vignette_coalescent_diversity.md

@@ -275,7 +275,7 @@ This runs quickly, since it's only 100 generations.
 First, let's look at what mutations are present.
 ```{code-cell}
 ts = tskit.load("vignette_annotated.trees")
-num_stacked = np.array([len(m.metadata["derived_states"]) for m in ts.mutations()])
+num_stacked = np.array([len(m.metadata["slim_ids"]) for m in ts.mutations()])
 init_time = ts.metadata['SLiM']['tick']
 old_mut = np.array([m.time > init_time - 1 - 1e-12 for m in ts.mutations()])
 assert sum(old_mut) == ots.num_mutations
@@ -311,7 +311,7 @@ p = ts.sample_count_stat(nodes_by_time, lambda x: x/num_nodes, 2, windows='sites
         strict=False, span_normalise=False, polarised=True)
 mut_metadata = pyslim.mutation_metadata(ts)
 s = np.array([sum([sum([mut_metadata[k]["per_trait"][0]["effect_size"]
-                        for k in m.metadata["derived_states"]])
+                        for k in m.metadata["slim_ids"]])
                   for m in site.mutations]) for site in ts.sites()])
 ```
 
@@ -422,7 +422,7 @@ and print a picture of it, with mutations labeled by their type:
 ```{code-cell}
 mut_metadata = pyslim.mutation_metadata(mts)
 for t in mts.trees():
- mt = [max([mut_metadata[k]['mutation_type'] for k in m.metadata["derived_states"]]) for m in t.mutations()]
+ mt = [max([mut_metadata[k]['mutation_type'] for k in m.metadata["slim_ids"]]) for m in t.mutations()]
  if t.num_mutations > 12:
    break
 
@@ -478,26 +478,26 @@ There are indeed:
 :tags: ['remove-output']
 for site in mts.sites():
   if len(site.mutations) > 1:
-     types = [set([mut_metadata[k]["mutation_type"] for k in mut.metadata["derived_states"]])
+     types = [set([mut_metadata[k]["mutation_type"] for k in mut.metadata["slim_ids"]])
               for mut in site.mutations]
      if max(map(len, types)) > 1:
         print(site)
         for mut in site.mutations:
             print(mut)
-            for k in mut.metadata["derived_states"]:
+            for k in mut.metadata["slim_ids"]:
                 print(mut_metadata[k])
 ```
 ```{code-cell}
 :tags: ['remove-input']
 for site in mts.sites():
   if len(site.mutations) > 1:
-     types = [set([mut_metadata[k]["mutation_type"] for k in mut.metadata["derived_states"]])
+     types = [set([mut_metadata[k]["mutation_type"] for k in mut.metadata["slim_ids"]])
               for mut in site.mutations]
      if max(map(len, types)) > 1:
         util.pp(site)
         for mut in site.mutations:
             util.pp(mut)
-            for k in mut.metadata["derived_states"]:
+            for k in mut.metadata["slim_ids"]:
                 util.pp(mut_metadata[k])
 ```
 

pyslim/methods.py

@@ -431,18 +431,18 @@ def add_mutation_metadata(ts, mutation_type=0, remove_unused=False):
     :class:`msprime.SLiMv6MutationModel`. Any information about SLiM mutations already
     in top-level metadata will remain unchanged.
 
-    To do this, this method looks for all SLiM IDs that are found in the derived
-    state of some mutation but are not represented in the top-level metadata
-    (see :func:`.mutation_metadata`). This function then adds entries to that top-level
-    metadata with default values (see :func:`.default_slim_metadata`),
+    To do this, this method looks for all SLiM IDs that are found in the `"slim_ids"`
+    entry of some tskit mutation's metadata but are not represented in the top-level
+    metadata (see :func:`.mutation_metadata`). This function then adds entries to that
+    top-level metadata with default values (see :func:`.default_slim_metadata`),
     except that (a) the ``mutation_type`` can be specified;
     and (b) the ``slim_time`` is set using the ``tick`` value in top-level metadata
     and the ``time`` of the oldest tskit mutation in which the SLiM mutation occurs.
 
     :param tskit.TreeSequence ts: The tree sequence to transform.
     :param int mutation_type: The numeric ID of the mutation type in SLiM.
     :param bool remove_unused: Whether to also remove from metadata information about any
-        mutations not seen in the derived states of the tree sequence.
+        mutations not referenced by mutations in the tree sequence.
     :return tskit.TreeSequence: A copy of the tree sequence with mutation information in
         metadata.
     """
@@ -460,7 +460,7 @@ def add_mutation_metadata_tables(tables, mutation_type=0, remove_unused=False):
     :param tskit.TableCollection tables: The table collection to be modified.
     :param int mutation_type: The numeric ID of the mutation type in SLiM.
     :param bool remove_unused: Whether to also remove from metadata information about any
-        mutations not seen in the derived states of the tree sequence.
+        mutations not referenced by mutations in the tree sequence.
     """
     ts_metadata = tables.metadata
     if (
@@ -475,9 +475,7 @@ def add_mutation_metadata_tables(tables, mutation_type=0, remove_unused=False):
     num_traits = len(ts_metadata["SLiM"]["traits"])
     existing_muts = {x["mutation_id"] for x in ts_metadata["SLiM_mutation_list"]}
     mut_ids = [
-        (int(j), mut.time)
-        for mut in tables.mutations
-        for j in mut.metadata["derived_states"]
+        (int(j), mut.time) for mut in tables.mutations for j in mut.metadata["slim_ids"]
     ]
     if len(mut_ids) > 0:
         mut_ids.sort()
@@ -516,13 +514,13 @@ def convert_alleles(ts):
     their corresponding nucleotides. For sites, SLiM-produced tree sequences
     have "" (the empty string) for the ancestral state at each site; this method
     will replace this with the corresponding nucleotide from the reference sequence.
-    For mutations, SLiM records the 'derived state' as a SLiM mutation ID; this
-    method will replace the derived state with the nucleotide from the mutation's
+    For mutations, SLiM records the 'derived state' as a list of SLiM mutation IDs;
+    this method will replace the derived state with the nucleotide from the mutation's
     metadata.
 
     In SLiM's output the list of mutation IDs is recorded both in each mutation's
-    derived state and metadata, but SLiM uses the derived state for loading files,
-    so the resulting tree sequence will not be loadable by SLiM.
+    derived state and metadata, but SLiM only uses the metadata for loading files,
+    so the resulting tree sequence will still be loadable by SLiM.
 
     The main purpose of this method is for output: for instance, this code will produce
     a VCF file with nucleotide alleles:
@@ -553,7 +551,7 @@ def convert_alleles(ts):
     alleles = np.array([x["nucleotide"] for x in mut_metadata.values()], dtype="int")
     # First, do this for the unstacked mutations quickly
     # mut_index will map from tskit-mutations to slim-mutations
-    mut_index = np.array([mut.metadata["derived_states"][-1] for mut in ts.mutations()])
+    mut_index = np.array([mut.metadata["slim_ids"][-1] for mut in ts.mutations()])
     assert np.all(mut_index >= 0), "This should not occur: please file a bug report."
     nucs = alleles[np.searchsorted(mut_ids, mut_index)]
     # Now, update those where necessary
@@ -605,7 +603,7 @@ def generate_nucleotides(ts, reference_sequence=None, keep=True, seed=None):
 
     Technical note: in the case of stacked mutations, the SLiM mutation that
     determines the nucleotide state of the (tskit) mutation is the last one in the list
-    of "derived states" in the tskit mutation metadata.  This method tries to
+    of "slim_ids" in the tskit mutation metadata.  This method tries to
     assign nucleotides so that each mutation differs from the previous state,
     but this is not always possible if some mutations already have nucleotides and
     others do not.
@@ -659,8 +657,8 @@ def generate_nucleotides(ts, reference_sequence=None, keep=True, seed=None):
                 pds = []
             else:
                 pa = states[mut.parent]
-                pds = ts.mutation(mut.parent).metadata["derived_states"]
-            for i in mut.metadata["derived_states"]:
+                pds = ts.mutation(mut.parent).metadata["slim_ids"]
+            for i in mut.metadata["slim_ids"]:
                 md = mut_info[i]
                 da = md["nucleotide"]
                 if da == -1 or not keep:
@@ -1232,7 +1230,7 @@ def next_slim_mutation_id(ts):
     `next_id` in your :class:`msprime.SLiMv6MutationModel` to be larger than any
     existing mutation IDs. Setting `next_id` equal to the output of this
     function will allow the mutated tree sequence to be read in by SLiM.
-    To do this, recall that the `derived_states` attribute of each mutation's metadata
+    To do this, recall that the `slim_ids` attribute of each mutation's metadata
     is a list of SLiM mutation IDs; this function just parses all these metadata entries
     and returns one larger than the largest integer found.
     """
@@ -1241,7 +1239,7 @@ def next_slim_mutation_id(ts):
         try:
             max_id = functools.reduce(
                 max,
-                (x for mut in ts.mutations() for x in mut.metadata["derived_states"]),
+                (x for mut in ts.mutations() for x in mut.metadata["slim_ids"]),
                 -1,
             )
         except TypeError: