documentation edits from Ben

CHANGELOG.rst

@@ -40,6 +40,10 @@
   may have affected users using `pyslim.convert(ts)` in a previous version of pyslim
   on a tree sequence `ts` with SLiM file version prior to 0.9.
 
+- Since verison 1.1, the value of `pyslim.INDIVIDUAL_FLAG_MIGRATED` has been 2,
+  when in fact it should have been 1, so code using this flag to detect
+  migrants would have been wrong (and should have found no migrants, ever).
+
 **New features**:
 
 - SLiM now includes in metadata information about the effects of mutations on

docs/.gitignore

@@ -4,3 +4,4 @@ spatial_sim_genotypes.vcf
 example_snps.vcf
 sims.make
 example_sim.vcf
+.ipynb_checkpoints

docs/_config.yml

@@ -3,7 +3,7 @@
 
 title: PySLiM manual
 author: Tskit Developers
-copyright: "2021"
+copyright: "2026"
 only_build_toc_files: true
 # logo: logo.png
 favicon: favicon.ico
@@ -55,3 +55,4 @@ sphinx:
       myst_enable_extensions:
       - colon_fence
       - deflist
+      nb_merge_streams: True

docs/_static/recomb_rates.tsv

@@ -1,5 +1,6 @@
-end_position	rate(cM/Mb)
-15000000	3.2
-50000000	2.5
-85000000	0.25
-99999999	2.8
+position	rate(cM/Mb)
+0	3.2
+15000000	2.5
+50000000	0.25
+85000000	2.8
+100000000	0

docs/installation.md

@@ -38,4 +38,4 @@ If you find a bug in ``pyslim`` or want to suggest an improvement, please
 If you have a question about using tree sequences,
 please ask it at [the tskit discussion page](https://github.com/tskit-dev/tskit/discussions).
 Finally, questions about SLiM should be directed to
-[the SLiM mailing list](https://groups.google.com/forum/#!forum/slim-discuss).
+[the slim-discuss mailing list](https://groups.google.com/forum/#!forum/slim-discuss).

docs/introduction.md

@@ -16,7 +16,8 @@ kernelspec:
 # Introduction
 
 This is the documentation for pyslim, a Python API
-for reading and modifying {ref}`tskit<tskit:sec_introduction>` tree sequence files
+for reading and modifying {ref}`tskit<tskit:sec_introduction>`
+[tree sequence](https://tskit.dev/learn/) files
 produced by [SLiM](https://messerlab.org/slim/), 
 or modifying files produced by other programs (e.g.,
 {ref}`msprime<msprime:sec_intro>`,
@@ -25,13 +26,24 @@ and [tsinfer](https://tsinfer.readthedocs.io/)) for use in SLiM.
 
 SLiM can read and write *tree sequences*, which store genetic genealogies
 for entire populations. These can be used to efficiently encode both the state of the
-population at various points during a simulation *as well as* the complete genomic
+population at various times during a simulation *as well as* the complete genomic
 ancestry. Furthermore, SLiM can "load" a saved tree sequence
 file to recreate the exact state of the population at the time it was saved.
-To do this, SLiM stores some additional information in the basic
-tree sequence file.
+To do this, SLiM stores some additional information (as {ref}`metadata<tskit:sec_metadata>`)
+in the basic tree sequence file.
 
 First, you probably want to read the {ref}`sec_overview`
 for a description of what (and *who*)
 SLiM records in the tree sequence.
 Next, the {ref}`sec_tutorial` walks through the common workflows.
+
+(sec_citation)=
+
+# Citation
+
+If you use pyslim, please cite:
+
+> *Bridging forward-in-time and coalescent simulations using pyslim*, by
+> Shyamalika Gopalan, Murillo F. Rodrigues, Peter L. Ralph, and Benjamin C. Haller
+> [bioRxiv:2025.09.30.679676](https://doi.org/10.1101/2025.09.30.679676). 
+

docs/metadata.md

@@ -66,13 +66,13 @@ and `ts.metadata["SLiM"]` contains information about the simulation:
 - `traits`: a list of information for each of the traits:
     * `index`: the index of the trait in SLiM
     * `name`: the name in SLiM for the trait
-    * `type`: additive, multiplicative, or logistic
+    * `type`: `"additive"`, `"multiplicative"`, or `"logistic"`
     * `baselineOffsetFromUser`: a value added to all traits
     * `baselineOffsetFromSubstitutions`: the total effect of all substitutions on the traits
       at the time the tree sequence was saved
     * `baselineAccumulation`: whether the effect of substitutions accumulate in that value
     * `directFitnessEffect`: whether the trait has a direct effect on fitness
-    * `individualOffsetMean`, `individualOffsetSD`: parameters governing the individual-level offsets
+    * `individualOffsetMean`, `individualOffsetSD`: parameters governing the distribution of individual-level offsets
       (i.e., "environment" effects)
 
 **Populations:**
@@ -83,7 +83,7 @@ in a nonspatial simulation). The metadata may be `None` for populations
 that SLiM did not use. The keys that SLiM uses are:
 
 - `slim_id`: the ID of this population in SLiM
-- `name`: the name of the population (by default, `p0`, `p1`, etcetera)
+- `name`: the name of the population (by default, `"p0"`, `"p1"`, etcetera)
 - `description`: a string describing the population
 - `selfing_fraction`, `female_cloning_fraction`, `male_cloning_fraction`, and `sex_ratio`: only present when applicable (e.g., in WF simulations)
 - `bounds_x0`, `bounds_x1`, `bounds_y0`, `bounds_y1`, `bounds_z0`, and `bounds_z1`: the spatial bounds, when applicable
@@ -95,7 +95,7 @@ Each individual produced by SLiM contains the following metadata:
 - `pedigree_id`: the "pedigree ID", unique within the SLiM simulation
 - `pedigree_p1`, `pedigree_p2`: the pedigree IDs of the individuals' two
   parents (they may be equal in the case of selfing, or `-1` to indicate no
-  parent, in the case of the initial generation or for cloning)
+  parent, in the case of the initial generation)
 - `age`: the `.age` property within SLiM at the time the file was written out
 - `subpopulation`: the subpopulation within SLiM the individual was in at the time the file was written out
 - `sex`: the sex of the individual (either {data}`.INDIVIDUAL_TYPE_FEMALE`, {data}`.INDIVIDUAL_TYPE_MALE`, or {data}`.INDIVIDUAL_TYPE_HERMAPHRODITE`)
@@ -106,21 +106,21 @@ Each individual produced by SLiM contains the following metadata:
   the purpose of `tagLX_set` is to record whether the tag has been set in the simulation
 - `per_trait`: a list of information about the trait values for this individual; these are in the same order
   as the traits listed in top-level metadata;
-    * `phenotype`: the trait value
+    * `phenotype`: the trait value (which is `NAN` if the phenotype has not been calculated, or was invalidated)
     * `offset`: the individual's offset (i.e., the "environmental effect")
 
 **Nodes:**
-Each "node" produced by SLiM (i.e., "genome" within SLiM) has:
+Each "node" produced by SLiM (i.e., "haplosome" within SLiM) has:
 
-- `slim_id`: the unique ID associated with the genome by SLiM
-- `is_vacant`: records the genome is a "vacant" genome (in which case it isn't
+- `slim_id`: the unique ID associated with the haplosome by SLiM
+- `is_vacant`: records that the node is a "vacant" node (in which case it isn't
   really there, so shouldn't have any mutations or relationships in the tree
   sequence!) - see [](sec_overview_vacant_nodes) for more explanation
 
 **Mutations:**
 Prior to SLiM 6.0, mutation metadata was associated with the tskit mutation objects.
 Now, this is stored in top-level metadata, under ``ts.metadata["SLiM_mutation_list"]``.
-Each entry
+Each entry contains:
 
 - `mutation_id`: the numeric ID of mutation in SLiM
 - `mutation_type`: the numeric ID of the `MutationType` within SLiM
@@ -231,8 +231,10 @@ md["SLiM"]["model_type"] = "nonWF"
 tables.metadata = md
 ```
 Modifying the top-level metadata
-could be used to set spatial bounds on an annotated msprime simulation, for instance.
-(This is recorded in the population metadata.)
+could be used to declare that an annotated msprime simulation
+is to be a spatial simulation in SLiM,
+by changing the "dimensionality" property in top-level metadata
+(and perhaps also the bounds in population metadata as well).
 
 
 ### Modifying SLiM metadata in tables
@@ -254,7 +256,7 @@ tables.individuals.clear()
 for ind in ts.individuals():
     md = ind.metadata
     md["age"] = random.choice([1,2,3,4])
-    _ = tables.individuals.append(
+    tables.individuals.append(
         ind.replace(metadata=md)
     )
 
@@ -273,8 +275,11 @@ mod_ts.dump("modified_ts.trees")
 
 ### Metadata entries
 
-SLiM records additional information in the metadata columns of Individual, Node, and Mutation tables,
+SLiM records additional information in the metadata columns of Individual and Node, and Mutation tables,
+and about mutations in top-level metadata,
 in a binary format using the python ``struct`` module.
+However, this is transparently taken care of by tskit's metadata module,
+so the user does not have to work directly with binary information or the struct module itself.
 See {ref}`tskit's metadata documentation <tskit:sec_metadata>`
 for details on how this works.
 Nothing besides this binary information can be stored in the metadata of these tables if the tree sequence is to be used by SLiM,

docs/overview.md

@@ -16,7 +16,7 @@ kernelspec:
 # Overview
 
 A tree sequence is a way of storing both the full genetic history and the genotypes
-of a bunch of genomes.
+of a potentially large collection of genome sequences.
 See {ref}`the tskit documentation<tskit:sec_introduction>`
 for more description of the tree sequence and underlying data structure,
 and definitions of the important terms.
@@ -25,10 +25,10 @@ and the "focal" nodes are called *sample nodes* or simply *samples*.
 Many operations on tree sequences act on the sample nodes by default (see the
 {ref}`tskit data model<tskit:sec_nodes_or_individuals>`
 for more on this topic), and the tree sequence always describes the genealogy of the
-entire genome of all the samples, at at least over the simulated time period.
+entire genome of all the samples, at least over the simulated time period.
 (Other nodes in the tree sequence represent ancestral genomes
 about which we might have only partial information).
-SLiM simulates diploid organisms, so each *individual* usually has two nodes;
+For simulated diploid organisms, each *individual* has two nodes, so
 many operations you might want to do involve first finding the individuals you want,
 and then looking at their nodes.
 
@@ -39,8 +39,10 @@ and then looking at their nodes.
 Suppose we've run a very small simulation with SLiM. The genetic relationships between
 the various diploid individuals who were alive over the course of the simulation might
 look something like the picture on the left below. Note that individuals (circles) are
-diploid, so that each contains two chromosomes or *nodes* (shaded rectangles),
+diploid, so that each contains two copies of the chromosome or *nodes* (shaded rectangles),
 and that relationships are between the *nodes*, not the individuals.
+(The term SLiM uses for each of these copies-of-chromosomes is "haplosome",
+a term we'll continue to use below.)
 
 At the end of the simulation we are typically only interested in the genetic
 relationships between the nodes in those individuals which are still alive;
@@ -57,7 +59,10 @@ name: pedigree01
 ---
 A conceptual diagram of
 **(A)** relationships between chromosomes of diploid individuals in a SLiM simulation, and
-**(B)** which information is returned in the tree sequence.
+**(B)** which information is retained in the tree sequence.
+Only pedigree relationships are shown, not genetic relationships:
+the diagram does not show which portions of genetic material were inherited
+along each relationship.
 ```
 
 (sec_left_in_tree_sequence)=
@@ -90,18 +95,25 @@ which individuals they were in (the tutorial explains ways to
 [retain this information](sec_retaining_individuals).
 
 As well as the historical individuals, many historical *nodes*  have been removed too,
-along with with their genealogical relationships (i.e. the lines, which in
+along with with their genealogical relationships (i.e., the lines, which in
 tree-sequence-speak are known as "edges"). The deleted nodes are simply those that are
 not needed to reconstruct the relationships between the sample nodes. For example, we
-remove nodes leading to a dead end (e.g. in individuals who had no offspring). Similarly,
-as time goes on, recombination events in conjunction with genetic drift can gradually
-reduce the genetic contribution of parts of older genomes to the current generation.
+remove nodes leading to a dead end (e.g., in individuals who had no offspring).
+Similarly, as time goes on, recombination events in conjunction with genetic drift can
+gradually reduce the genetic contribution of parts of older genomes to the current
+generation.
 The generated tree sequence therefore need not contain historical nodes whose genetic
-contribution to the samples has been whittled down to zero. Finally, to
-reconstruct relationships between samples, strictly we only need to keep a node if it
-represents the genetic *most recent common ancestor* (MRCA) of at least two samples. So
-by default, we also remove historical nodes that are only "on the line to" a sample, but
-do not represent a branching point (i.e. coalescent event) on the tree.
+contribution to the samples has been whittled down to zero.
+Finally, to reconstruct relationships between samples, we only need to keep a node if it
+represents the genetic *most recent common ancestor* (MRCA) of at least two samples
+at some point on the genome.
+We will also need to keep nodes in the founding generation
+to describe relationships of the samples to this initial time.
+So by default, SLiM also removes historical nodes that are only "on the line to" a sample,
+but do not represent a either a branching point (i.e. coalescent event) or a root
+of a the marginal tree describing genetic relationships at some point on the genome.
+(See [Kelleher et al 2018](https://doi.org/10.1371/journal.pcbi.1006581)
+for a description of what relationships, exactly, are retained.)
 
 
 ## What else can I find out from the tree sequence?
@@ -115,9 +127,8 @@ Most of this is stored as *metadata*: see [](sec_metadata).
 
 ## Vacant nodes: sex chromosomes and haploidy
 
-Under the hood in SLiM, all individuals are diploid,
-and so for each individual there are two *nodes*,
-representing their two haplosomes
+Under the hood, SLiM reserves two *nodes* for each individual
+to represent their haplosomes
 (i.e., chromosome copies; see the SLiM manual).
 So, any individual that is *not* diploid for whichever chromosome
 the tree sequence represents will have in the tree sequence
@@ -166,9 +177,10 @@ in all the tree sequences in a trees archive are identical.
 This is quite handy if you're going to analyze data from more than one chromosome,
 but also creates some potential pitfalls.
 
-First, there will probably be nodes in the tree sequence
-that are not represented in *any* of the trees.
-This is simply because these are nodes that are needed for other chromosomes:
+First, in one particular tree sequence of a trees archive,
+there will probably be nodes that are not represented in *any* of the trees.
+This is simply because these are nodes that are referenced by
+the tree sequence for other chromosomes:
 perhaps the node in question is not ancestral to any of the samples
 on this chromosome, but it is on other chromosomes.
 The presence of these nodes is harmless, however: it will not affect
@@ -193,7 +205,7 @@ if both were from the same multichromosome simulation,
 and vice-versa.
 
 Finally, this explains why vacant nodes must be included in the tree sequence:
-in the presence of autosomes, every diploid individual needs two nodes,
+in the presence of diploid autosomes, every diploid individual is allocated two nodes,
 but for an individual that has less than two copies of a given chromosome,
 some of those nodes will not represent an actual chromosome copy.
 Those are "vacant", as described above.

docs/phenotypes.slim

@@ -0,0 +1,66 @@
+initialize() {
+	setSeed(123);
+	initializeTreeSeq();
+	defineConstant("START1", -5.0);
+	defineConstant("START2", 5.0);
+	defineConstant("OPT1", 20.0);
+	defineConstant("OPT2", -20.0);
+	initializeTrait("trait1", "additive", baselineOffset=START1, individualOffsetSD=1.0);
+	initializeTrait("trait2", "additive", baselineOffset=START2, individualOffsetSD=1.0);
+
+	initializeMutationType("m1", NAN, "f", 0.0);   // neutral
+	initializeMutationType("m2", NAN, "f", 0.0);   // QTLs
+	m2.color = "red";
+	m2.logMutationData(enable=T, effectSize=T);
+	m2.convertToSubstitution = F; // necessary to disentagle from baseline offset
+
+	// g1 is a neutral region, g2 is a QTL
+	initializeGenomicElementType("g1", m1, 1.0);
+	initializeGenomicElementType("g2", c(m1,m2), c(1.0, 0.1));
+
+	// chromosome of length 100 kb with two QTL regions
+	initializeGenomicElement(g1, 0, 39999);
+	initializeGenomicElement(g2, 40000, 49999);
+	initializeGenomicElement(g1, 50000, 79999);
+	initializeGenomicElement(g2, 80000, 89999);
+	initializeGenomicElement(g1, 90000, 99999);
+	initializeRecombinationRate(1e-8);
+	initializeMutationRate(1e-7);
+
+	// QTL-related constants used below
+	defineConstant("QTL_mu", c(0.0, 0.0));
+	defineConstant("QTL_cov", 0.25);
+	defineConstant("QTL_sigma", matrix(c(1,QTL_cov,QTL_cov,1), nrow=2));
+	}
+
+1 late() {
+	sim.addSubpop("p1", 500);
+}
+
+mutation(m2) {
+	// draw mutational effects for the new m2 mutation
+	effects = rmvnorm(1, QTL_mu, QTL_sigma);
+	mut.setEffectSizeForTrait(NULL, effects);
+	return T;
+}
+
+late() {
+	sim.demandPhenotype(NULL, NULL);
+
+	inds = sim.subpopulations.individuals;
+	effects1 = 1.0 + dnorm(inds.trait1, OPT1, 15.0) * 10.0;
+	effects2 = 1.0 + dnorm(inds.trait2, OPT2, 15.0) * 10.0;
+	inds.fitnessScaling = effects1 * effects2;
+}
+
+1:1000000 late() {
+	if (sim.cycle % 5000 == 0) {
+		inds = p1.individuals;
+		sim.treeSeqRememberIndividuals(inds);
+		// keep running until we get close to both optima
+		if ((abs(mean(inds.trait1) - OPT1) <= 1.0) & (abs(mean(inds.trait2) - OPT2) <= 1.0)) {
+	   	sim.treeSeqOutput("phenotypes.trees", metadata=Dictionary("all_mutations", m2.loggedData("values")));
+			sim.simulationFinished();
+		}
+	}
+}

docs/phenotypes2.slim

@@ -0,0 +1,24 @@
+initialize() {
+	setSeed(123);
+	initializeTreeSeq();
+	initializeTrait("mult", "multiplicative", individualOffsetSD=1.0);
+	initializeTrait("logistic", "logistic", individualOffsetSD=1.0);
+
+	initializeMutationType("m1", NAN, "n", 0.0, 0.01); 
+	m1.convertToSubstitution = F; // necessary to disentagle from baseline offset
+
+	initializeGenomicElementType("g1", m1, 1.0);
+
+	initializeGenomicElement(g1, 0, 9999999);
+	initializeRecombinationRate(1e-8);
+	initializeMutationRate(1e-7);
+}
+
+1 late() {
+	sim.addSubpop("p1", 100);
+}
+
+1000 late() {
+	sim.demandPhenotype(NULL, NULL);
+	sim.treeSeqOutput("phenotypes2.trees");
+}