From ffc6cd13da66932100b1b90c691e3fd9ad348e31 Mon Sep 17 00:00:00 2001 From: peter Date: Thu, 27 Aug 2026 11:02:29 -0700 Subject: [PATCH] added new slim v6 model --- CHANGELOG.md | 9 ++ docs/api.md | 5 + docs/mutations.md | 170 ++++++++++++++++---- lib/msprime.h | 4 + lib/mutgen.c | 244 ++++++++++++++++++++++------ lib/tests/test_mutations.c | 318 ++++++++++++++++++++++++++++++++++++- msprime/__init__.py | 2 + msprime/_msprimemodule.c | 115 +++++++++++++- msprime/mutations.py | 49 +++++- tests/test_mutations.py | 295 ++++++++++++++++++++++++++++------ tests/test_python_c.py | 33 ++++ 11 files changed, 1105 insertions(+), 139 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e92ec4573..de0d42023 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,15 @@ In development. +**New features**: + +- A new mutation model, ``SLiMv6MutationModel``, is introduced to match changes in + SLiM 6.0: this mutation behaves as ``SLiMMutationModel``, but without the ``type`` + and ``slim_generation`` arguments, and stores different information in metadata + (the list of mutation IDs, in binary). ``SLiMMutationModel`` will remain for + backwards compatibility, but code for use with SLiM versions 6.0 and above should + use ``SLiMv6MutationModel`` ({issue}`2534`, {pr}`2539`, {user}`petrelharp`) + ## [1.4.2] - 2026-05-20 Bugfix release. diff --git a/docs/api.md b/docs/api.md index 94b772ba0..e1ecd58ef 100644 --- a/docs/api.md +++ b/docs/api.md @@ -48,6 +48,7 @@ for discussion and examples of individual features. TPM EL2 InfiniteAlleles + SLiMv6MutationModel SLiMMutationModel ``` @@ -217,6 +218,10 @@ for discussion and examples of individual features. .. autoclass:: msprime.InfiniteAlleles() ``` +```{eval-rst} +.. autoclass:: msprime.SLiMv6MutationModel() +``` + ```{eval-rst} .. autoclass:: msprime.SLiMMutationModel() ``` diff --git a/docs/mutations.md b/docs/mutations.md index 33d4d6638..6775948e2 100644 --- a/docs/mutations.md +++ b/docs/mutations.md @@ -105,8 +105,11 @@ for more information. {class}`.InfiniteAlleles` (Integers) : A generic infinite-alleles mutation model +{class}`.SLiMv6MutationModel` (Integers) +: An infinite-alleles model producing SLiM-style mutations for SLiM v6+ + {class}`.SLiMMutationModel` (Integers) -: An infinite-alleles model producing SLiM-style mutations +: An infinite-alleles model producing SLiM-style mutations for older versions of SLiM {class}`.MatrixMutationModel` (General finite state model) : Superclass of mutation models with a finite set of states @@ -485,7 +488,8 @@ Here are the available models; they are documented in more detail below. - {class}`.TPM`: Two-phase mutation model for microsatellite repeat copy number. DiRienzo et al. ('94) - {class}`.EL2`: Two-phase mutation model, equal rate, linear bias model for microsatellite repeat copy number. Garza et al. ('95) - {class}`.InfiniteAlleles`: A generic infinite-alleles mutation model -- {class}`.SLiMMutationModel`: An infinite-alleles model of mutation producing SLiM-style mutations +- {class}`.SLiMv6MutationModel`: An infinite-alleles model of mutation producing SLiM-style mutations for SLiM v6+ +- {class}`.SLiMMutationModel`: An infinite-alleles model of mutation producing SLiM-style mutations for older versions of SLiM (sec_mutations_matrix_mutations_models)= @@ -855,30 +859,34 @@ to set the starting allele appropriately, and to make sure the results make sens A special class of infinite alleles model is provided for use with [SLiM](), to agree with the underlying mutation model in SLiM. -As with the InfiniteAlleles model, it assigns each new mutation a unique integer, -by keeping track of the `next_id` and incrementing it each time a new mutation appears. +As with the {class}`.InfiniteAlleles` model, it assigns each new mutation +a unique integer, by keeping track of the `next_id` and incrementing it each +time a new mutation appears. For more information, +see {ref}`the pyslim documentation`. -This differs from the {class}`.InfiniteAlleles` because mutations +This differs from the {class}`.InfiniteAlleles` model because mutations in SLiM can "stack": new mutations can add to the existing state, rather than replacing the previous state. So, derived states are comma-separated lists of mutation IDs, and the ancestral state is always the empty string. For instance, if a new mutation with ID 5 occurs at a site, and then later another mutation appears with ID 64, the sequence of alleles moving along this line of descent -would be `""`, then `"5"`, and finally `"5,64"`. Furthermore, the mutation -model adds SLiM metadata to each mutation, which records, among other things, -the SLiM mutation type of each mutation, and the selection coefficient (which -is always 0.0, since adding mutations in this way only makes sense if they are -neutral). For this reason, the model has one required parameter: the `type` -of the mutation, a nonnegative integer. If, for instance, you specify -`type=1`, then the mutations in SLiM will be of type `m1`. For more -information, and for how to modify the metadata (e.g., changing the selection -coefficients), see -{ref}`the pyslim documentation`. -For instance, +would be `""`, then `"5"`, and finally `"5,64"`. The same list of integers is +recorded in metadata, as binary integers rather than ASCII text. +This redundancy between derived state and metadata is useful +because the model behaves as an infinite-alleles model out-of-the-box; but if +ancestral and derived states are modified (for instance, to be VCF-compliant) +then there is no loss of information. + +There are in fact two mutation models for SLiM: +{class}`.SLiMv6MutationModel` and {class}`.SLiMMutationModel`, +because the format that SLiM uses to store mutations changed with version 6.0. +Both store the same thing in derived state; the difference is in metadata. +So, you probably want to use {class}`.SLiMv6MutationModel`; +see {ref}`sec_mutations_mutation_slim_mutations_versions` for more on the differences. ```{code-cell} python -model = msprime.SLiMMutationModel(type=1) +model = msprime.SLiMv6MutationModel() mts = msprime.sim_mutations( ts, rate=1, random_seed=1, model=model) t = mts.first() @@ -888,14 +896,16 @@ SVG(t.draw_svg(mutation_labels=ml, node_labels={}, size=(400, 300))) ``` These resulting alleles show how derived states are built. +(We're looking at derived states for convenience, but could equivalently +look at metadata; see {ref}`sec_mutations_mutation_slim_mutations_metadata`.) The behaviour of this mutation model when used to add mutations to a previously mutated tree sequence can be subtle. Let's look at a simple example. -Here, we first lay down mutations of type 1, starting from ID 0: +Here, we first lay down mutations starting from ID 0 (the default): ```{code-cell} python -model_1 = msprime.SLiMMutationModel(type=1) +model_1 = msprime.SLiMv6MutationModel() mts_1 = msprime.sim_mutations(ts, rate=0.5, random_seed=2, model=model_1) t = mts_1.first() ml = {m.id: m.derived_state for m in mts_1.mutations()} @@ -903,15 +913,16 @@ SVG(t.draw_svg(mutation_labels=ml, node_labels={}, size=(400, 300))) ``` -Next, we lay down mutations of type 2. -These we assign starting from ID 100, -to make it easy to see which are which: -in general just need to make sure that we start at an ID greater than any +Now suppose we wanted to add more mutations. +(This is essentially what happens when we use msprime to add mutations +to a tree sequence produced by SLiM.) +We assign the new mutations starting from ID 100 to make it easy to see which are which: +in general we just need to make sure that we start at an ID greater than any previously assigned. ```{code-cell} python -model_2 = msprime.SLiMMutationModel(type=2, next_id=100) +model_2 = msprime.SLiMv6MutationModel(next_id=100) mts = msprime.sim_mutations( mts_1, rate=0.5, random_seed=3, model=model_2, keep=True) t = mts.first() @@ -928,7 +939,110 @@ with IDs `100` and `102`, between these two mutations. These were added to mutation `0`, obtaining alleles `0,100` and `0,100,102`. But then, moving down the branch, we come upon the mutation with ID `3`. This was already present in the tree sequence, so its derived state is not modified: -`0,3`. We can rationalise this, post-hoc, by saying that the type 1 mutation `3` -has "erased" the type 2 mutations `100` and `102`. -If you want a different arrangement, -you can go back and edit the derived states (and metadata) as you like. +`0,3`. We can rationalise this, post-hoc, by saying that the mutation with ID `3` +from the earlier batch of mutations has "erased" the mutations with IDs `100` and `102` +from the later batch. If you want a different arrangement, +you will need to edit the derived states and metadata directly. + +(sec_mutations_mutation_slim_mutations_metadata)= + +#### Metadata + +The most common use of this model is to add mutations to a tree sequence produced +by SLiM, or by {func}`pyslim.annotate`; if so, the metadata will already be decoded. +(In general, we recommend using {func}`pyslim.annotate` after +{func}`.sim_ancestry` but before {func}`.sim_mutations` if you're using msprime +to create a tree sequence for SLiM.) + +For completeness, here is a self-contained example +that shows how to decode the metadata directly +(but we emphasize that you probably *don't* need this in your code; +it's provided for pedagogical reasons only!). +If we look at the metadata in the example we've produced so far, we simply see binary +(the metadata is not decoded because a metadata schema is not present): + +```{code-cell} python +mut = mts.mutation(10) +mut.metadata +``` + +To see the actual values, we need a {class}`tskit.MetadataSchema`, +which we can add as follows: + +```{code-cell} python +metadata_schema = tskit.MetadataSchema( + { + "codec": "struct", + "type": "object", + "properties": { + "derived_states": { + "items": { + "binaryFormat": "q", + "type": "number", + }, + "noLengthEncodingExhaustBuffer": True, + "type": "array", + } + }, + "required": ["derived_states"], + "additionalProperties": False, + } +) + +tables = mts.dump_tables() +tables.mutations.metadata_schema = metadata_schema +mts = tables.tree_sequence() +``` + +(This is a simplified version of the schema used by SLiM 6, +so it matches the metadata produced by the {class}`.SLiMv6MutationModel`.) +Now, the metadata will be properly decoded, into lists of integer IDs: + +```{code-cell} python +mut = mts.mutation(10) +mut.metadata +``` + +(sec_mutations_mutation_slim_mutations_stacking)= + +#### More on "mutation stacking" + +Both SLiM models simply create the derived state by appending to the inherited state, +and so if you use {func}`.sim_mutations` with ``keep=True`` to add SLiM mutations +to a tree sequence that already has non-SLiM mutations, you may end up with +mutations with odd derived states. For instance, if a SLiM mutations occurs +at the same site as, and inheriting from, a previous mutation +with derived state ``"A"``, the SLiM mutation might have derived state ``"A,345"``. +For a mutation without a parent, this inherited state is the ancestral state +of the site. If the site was not already present in the tree sequence, this model +will create a new site with the empty string as the ancestral state, +but if the site is already present and has a nonempty ancestral state, +then this ancestral state will appear at the beginning of any derived states +produced by these models at that site. + +The same thing is true of metadata: metadata is simply appended (in binary) +to the "inherited" metadata, which is from the parent mutation (if any) or from +the site (otherwise). So, if in your tree sequence your sites have nonempty metadata, +then any SLiM mutations that occur on those sites will have site metadata +prepended to their metadata (and so the metadata schema above will not work). +However, these SLiM mutation models are designed specifically to work with SLiM-format +tree sequences, which don't have any site metadata; +if you're trying to use them in another situation then you probably +don't have the right tool for the job. + +(sec_mutations_mutation_slim_mutations_versions)= + +#### Previous versions + +The technical difference between +{class}`.SLiMv6MutationModel` and {class}`.SLiMMutationModel` is as follows. +{class}`.SLiMv6MutationModel` stores the SLiM IDs in binary metadata +as well as in the derived state, +so that the actual metadata about mutations can be looked up from where +SLiM expects it in top-level metadata. +The older version, {class}`.SLiMMutationModel`, +stores the metadata about SLiM mutations directly in each tskit mutation's metadata +(and since individual SLiM mutations can appear in more than one tskit mutation, +this leads to some redundancy; thus the change +to how mutation metadata is handled in SLiMv6). For more information, see +{ref}`the pyslim documentation`. diff --git a/lib/msprime.h b/lib/msprime.h index 742d17564..3a47b34e6 100644 --- a/lib/msprime.h +++ b/lib/msprime.h @@ -439,6 +439,8 @@ typedef struct { double *transition_matrix; } mutation_matrix_t; +// this one is used by both SLiM mutation models +// (with mutation_type_id and slim_generation unused by v6) typedef struct { int32_t mutation_type_id; // following SLiM's MutationMetadataRec int64_t next_mutation_id; // following SLiM's slim_mutationid_t @@ -584,6 +586,8 @@ int matrix_mutation_model_alloc(mutation_model_t *self, size_t num_alleles, double *transition_matrix); int slim_mutation_model_alloc(mutation_model_t *self, int32_t mutation_type_id, int64_t next_mutation_id, int32_t slim_generation, size_t block_size); +int slim_v6_mutation_model_alloc( + mutation_model_t *self, int64_t next_mutation_id, size_t block_size); int infinite_alleles_mutation_model_alloc( mutation_model_t *self, uint64_t start_allele, tsk_flags_t options); int mutation_model_free(mutation_model_t *self); diff --git a/lib/mutgen.c b/lib/mutgen.c index 946e1dec8..820168972 100644 --- a/lib/mutgen.c +++ b/lib/mutgen.c @@ -255,7 +255,94 @@ mutation_matrix_free(mutation_model_t *self) } /*********************** - * SLiM mutation model */ + * SLiM mutation models */ + +/* Code shared between the SLiMMutationModel and SLiMv6MutationModel */ + +static int +slim_both_mutator_choose_root_state( + mutation_model_t *MSP_UNUSED(self), gsl_rng *MSP_UNUSED(rng), site_t *site) +{ + site->ancestral_state = NULL; + site->ancestral_state_length = 0; + return 0; +} + +static int +slim_both_mutator_transition(mutation_model_t *self, const char *parent_allele, + tsk_size_t parent_allele_length, const char *parent_metadata, + tsk_size_t parent_metadata_length, mutation_t *mutation, + const tsk_size_t new_metadata_length, char **metadata_buffer) +{ + /* This function performs shared work for both SLiM models: + * - copying the inherited state into the derived state, + * and appending the SLiM ID to that derived state + * - copying over the previous metadata + * but NOT appending the new metadata (since that differs between models); + * the *metadata_buffer argument here exists so that the model-specific + * function knows where to put the new bit of metadata + * - it is that client code's responsibility to make sure that no + * more than parent_metadata_length + new_metadata_length + * gets copied into *metadata_buffer + * - it is also the client code's responsibility to increment + * params->next_mutation_id, since that is needed for the new metadata + */ + int ret = 0; + slim_mutator_t *params = &self->params.slim_mutator; + char *buff = NULL; + int len; + /* The maximum number of digits for a signed 64 bit integer (including + * the leading "-") */ + const size_t max_digits = 20; + /* We allow for a possible comma to separate the previous element + * in the list, as well as a NULL byte added by snprintf. We don't bother + * trying to alloc the exact number of bytes needed. */ + const size_t alloc_size = parent_allele_length + max_digits + 2; + const char *sep = parent_allele_length == 0 ? "" : ","; + const tsk_size_t metadata_length = parent_metadata_length + new_metadata_length; + + /* Append to derived_state */ + buff = tsk_blkalloc_get(¶ms->allocator, alloc_size); + if (buff == NULL) { + ret = MSP_ERR_NO_MEMORY; + goto out; + } + len = snprintf(buff, alloc_size, "%.*s%s%" PRId64, (int) parent_allele_length, + parent_allele, sep, params->next_mutation_id); + if (len < 0) { + /* Technically this can happen. Returning TSK_ERR_IO should result + * in Python code checking errno. */ + ret = TSK_ERR_IO; + goto out; + } + tsk_bug_assert(len < (int) alloc_size); + mutation->derived_state = buff; + mutation->derived_state_length = (tsk_size_t) len; + + /* Append to metadata */ + buff = tsk_blkalloc_get(¶ms->allocator, metadata_length); + if (buff == NULL) { + ret = MSP_ERR_NO_MEMORY; + goto out; + } + memcpy(buff, parent_metadata, parent_metadata_length); + // tell client code where to put the new metadata + *metadata_buffer = buff; + mutation->metadata_length = (tsk_size_t) (metadata_length); +out: + return ret; +} + +static int +slim_both_mutator_free(mutation_model_t *self) +{ + slim_mutator_t params = self->params.slim_mutator; + tsk_blkalloc_free(¶ms.allocator); + return 0; +} + +/******************************* + * (Older) SLiM mutation model */ /* Typedefs from MutationMetadataRec in slim_sim.h: * line 125 in v3.4, git hash b2c2b634199f35e53c4e7e513bd26c91c6d99fd9 @@ -337,79 +424,97 @@ slim_mutator_check_validity(slim_mutator_t *self) return ret; } -static int -slim_mutator_choose_root_state( - mutation_model_t *MSP_UNUSED(self), gsl_rng *MSP_UNUSED(rng), site_t *site) -{ - site->ancestral_state = NULL; - site->ancestral_state_length = 0; - return 0; -} - static int slim_mutator_transition(mutation_model_t *self, gsl_rng *MSP_UNUSED(rng), const char *parent_allele, tsk_size_t parent_allele_length, const char *parent_metadata, tsk_size_t parent_metadata_length, mutation_t *mutation) { int ret = 0; + char *buff; slim_mutator_t *params = &self->params.slim_mutator; - char *buff = NULL; - int len; - /* The maximum number of digits for a signed 64 bit integer (including - * the leading "-") */ - const size_t max_digits = 20; - /* We allow for a possible comma to separate the previous element - * in the list, as well as a NULL byte added by snprintf. We don't bother - * trying to alloc the exact number of bytes needed. */ - const size_t alloc_size = parent_allele_length + max_digits + 2; - const char *sep = parent_allele_length == 0 ? "" : ","; + const tsk_size_t new_metadata_length = SLIM_MUTATION_METADATA_SIZE; - /* Append to derived_state */ - buff = tsk_blkalloc_get(¶ms->allocator, alloc_size); - if (buff == NULL) { - ret = MSP_ERR_NO_MEMORY; - goto out; - } - len = snprintf(buff, alloc_size, "%.*s%s%" PRId64, (int) parent_allele_length, - parent_allele, sep, params->next_mutation_id); - if (len < 0) { - /* Technically this can happen. Returning TSK_ERR_IO should result - * in Python code checking errno. */ - ret = TSK_ERR_IO; + ret = slim_both_mutator_transition(self, parent_allele, parent_allele_length, + parent_metadata, parent_metadata_length, mutation, new_metadata_length, &buff); + if (ret < 0) { goto out; } - tsk_bug_assert(len < (int) alloc_size); + + // copy the final bit of metadata in (this is different than the v6 model) + copy_slim_mutation_metadata(params, buff + parent_metadata_length, mutation->time); + mutation->metadata = buff; + + // finally, set up for the next one if (params->next_mutation_id == INT64_MAX) { ret = MSP_ERR_MUTATION_ID_OVERFLOW; goto out; } params->next_mutation_id++; - mutation->derived_state = buff; - mutation->derived_state_length = (tsk_size_t) len; - /* Append to metadata */ - buff = tsk_blkalloc_get( - ¶ms->allocator, parent_metadata_length + SLIM_MUTATION_METADATA_SIZE); - if (buff == NULL) { - ret = MSP_ERR_NO_MEMORY; +out: + return ret; +} + +/************************* + * SLiMv6 mutation model */ + +/* To be fully SLiM-compatible, there also needs to be some information in + * top-level metadata. But, it is not our job to add it in here, and that + * should be added after the fact with the method pyslim.add_mutation_metadata(). + * */ + +static void +slim_v6_mutator_print_state(mutation_model_t *self, FILE *out) +{ + slim_mutator_t params = self->params.slim_mutator; + fprintf(out, "SLiMv6 mutation model :: next mutation ID = %d\n", + (int) params.next_mutation_id); +} + +static int MSP_WARN_UNUSED +slim_v6_mutator_check_validity(slim_mutator_t *self) +{ + int ret = 0; + + if (self->next_mutation_id < 0) { + ret = MSP_ERR_BAD_SLIM_PARAMETERS; goto out; } - memcpy(buff, parent_metadata, parent_metadata_length); - copy_slim_mutation_metadata(params, buff + parent_metadata_length, mutation->time); - mutation->metadata = buff; - mutation->metadata_length - = (tsk_size_t) (parent_metadata_length + SLIM_MUTATION_METADATA_SIZE); out: return ret; } static int -slim_mutator_free(mutation_model_t *self) +slim_v6_mutator_transition(mutation_model_t *self, gsl_rng *MSP_UNUSED(rng), + const char *parent_allele, tsk_size_t parent_allele_length, + const char *parent_metadata, tsk_size_t parent_metadata_length, mutation_t *mutation) { - slim_mutator_t params = self->params.slim_mutator; - tsk_blkalloc_free(¶ms.allocator); - return 0; + int ret = 0; + char *buff; + slim_mutator_t *params = &self->params.slim_mutator; + const tsk_size_t new_metadata_length = sizeof(int64_t); + + ret = slim_both_mutator_transition(self, parent_allele, parent_allele_length, + parent_metadata, parent_metadata_length, mutation, new_metadata_length, &buff); + if (ret < 0) { + goto out; + } + + // copy the final bit of metadata in (this is different than the non-v6 model) + memcpy( + buff + parent_metadata_length, ¶ms->next_mutation_id, new_metadata_length); + mutation->metadata = buff; + + // finally, set up for the next one + if (params->next_mutation_id == INT64_MAX) { + ret = MSP_ERR_MUTATION_ID_OVERFLOW; + goto out; + } + params->next_mutation_id++; + +out: + return ret; } /*********************** @@ -558,10 +663,10 @@ slim_mutation_model_alloc(mutation_model_t *self, int32_t mutation_type_id, memset(self, 0, sizeof(*self)); - self->choose_root_state = &slim_mutator_choose_root_state; + self->choose_root_state = &slim_both_mutator_choose_root_state; self->transition = &slim_mutator_transition; self->print_state = &slim_mutator_print_state; - self->free = &slim_mutator_free; + self->free = &slim_both_mutator_free; if (block_size == 0) { /* 8K is a good default, but we need to have the * block_size argument here because the size of the allocations @@ -588,6 +693,43 @@ slim_mutation_model_alloc(mutation_model_t *self, int32_t mutation_type_id, return ret; } +int MSP_WARN_UNUSED +slim_v6_mutation_model_alloc( + mutation_model_t *self, int64_t next_mutation_id, size_t block_size) +{ + int ret = 0; + slim_mutator_t *params = &self->params.slim_mutator; + + memset(self, 0, sizeof(*self)); + + self->choose_root_state = &slim_both_mutator_choose_root_state; + self->transition = &slim_v6_mutator_transition; + self->print_state = &slim_v6_mutator_print_state; + self->free = &slim_both_mutator_free; + if (block_size == 0) { + /* 8K is a good default, but we need to have the + * block_size argument here because the size of the allocations + * we need to make are in principle unbounded since SLiM + * copies the entire parent state for every mutation as it + * goes down along the tree. + */ + block_size = 8192; + } + ret = tsk_blkalloc_init(¶ms->allocator, block_size); + if (ret != 0) { + ret = msp_set_tsk_error(ret); + goto out; + } + params->next_mutation_id = next_mutation_id; + + ret = slim_v6_mutator_check_validity(params); + if (ret != 0) { + goto out; + } +out: + return ret; +} + int MSP_WARN_UNUSED infinite_alleles_mutation_model_alloc( mutation_model_t *self, uint64_t start_allele, tsk_flags_t MSP_UNUSED(options)) diff --git a/lib/tests/test_mutations.c b/lib/tests/test_mutations.c index ccc8d78f4..48c5e5abf 100644 --- a/lib/tests/test_mutations.c +++ b/lib/tests/test_mutations.c @@ -37,6 +37,9 @@ insert_single_tree(tsk_table_collection_t *tables, int alphabet) "0 1 6 4,5\n"; */ int ret; + const char *mut_metadata = "lmnop"; + const char *site_metadata = "abc"; + tables->sequence_length = 1.0; ret = tsk_node_table_add_row( &tables->nodes, TSK_NODE_IS_SAMPLE, 0.0, 0, TSK_NULL, NULL, 0); @@ -73,18 +76,18 @@ insert_single_tree(tsk_table_collection_t *tables, int alphabet) ret = tsk_population_table_add_row(&tables->populations, NULL, 0); CU_ASSERT_FATAL(ret == 0); - /* Add a site and a mutation */ + /* Add a site and a mutation with metadata, above node 4 */ if (alphabet == ALPHABET_BINARY) { - ret = tsk_site_table_add_row(&tables->sites, 0.1, "0", 1, NULL, 0); + ret = tsk_site_table_add_row(&tables->sites, 0.1, "0", 1, site_metadata, 3); CU_ASSERT_FATAL(ret >= 0); ret = tsk_mutation_table_add_row( - &tables->mutations, 0, 0, -1, 0.0, "1", 1, NULL, 0); + &tables->mutations, 0, 4, -1, 2.5, "1", 1, mut_metadata, 5); CU_ASSERT_FATAL(ret >= 0); } else if (alphabet == ALPHABET_NUCLEOTIDE) { - ret = tsk_site_table_add_row(&tables->sites, 0.1, "A", 1, NULL, 0); + ret = tsk_site_table_add_row(&tables->sites, 0.1, "A", 1, site_metadata, 3); CU_ASSERT_FATAL(ret >= 0); ret = tsk_mutation_table_add_row( - &tables->mutations, 0, 0, -1, 0.0, "C", 1, NULL, 0); + &tables->mutations, 0, 4, -1, 2.5, "C", 1, mut_metadata, 5); CU_ASSERT_FATAL(ret >= 0); } @@ -864,13 +867,19 @@ parse_text_int64(char *ds, tsk_size_t n) static void verify_slim_mutation_ids(int64_t *mut_ids, size_t mut_ids_length, int64_t min_mut_id) { - int j; + int j = 0; CU_ASSERT_FATAL(mut_ids_length > 0); qsort(mut_ids, mut_ids_length, sizeof(*mut_ids), &cmp_int64); - CU_ASSERT_EQUAL_FATAL(mut_ids[0], min_mut_id); + // skip any -1s we might have put in + while (mut_ids[j] < 0) { + j++; + } + CU_ASSERT_EQUAL_FATAL(mut_ids[j], min_mut_id); + j++; if (mut_ids_length > 1) { - for (j = 1; j < mut_ids_length; j++) { + while (j < mut_ids_length) { CU_ASSERT_EQUAL_FATAL(mut_ids[j] - mut_ids[j - 1], 1); + j++; } } } @@ -1056,6 +1065,265 @@ test_mutgen_slim_mutation_large_values(void) gsl_rng_free(rng); } +static void +verify_slim_v6_derived_states(tsk_site_table_t *sites, tsk_mutation_table_t *mutations, + int64_t *all_mut_ids, size_t all_mut_ids_length, int skip_mut) +{ + int j, k, s; + size_t len, parent_len; + char *ds; + int64_t mut_id; + char *parent_allele; + CU_ASSERT_FATAL(mutations->num_rows <= all_mut_ids_length); + for (j = 0; j < mutations->num_rows; j++) { + ds = mutations->derived_state + mutations->derived_state_offset[j]; + len = (mutations->derived_state_offset[j + 1] + - mutations->derived_state_offset[j]); + if (j == skip_mut) { + parent_len = 0; + all_mut_ids[j] = -1; + } else { + k = mutations->parent[j]; + if (k == TSK_NULL) { + s = mutations->site[j]; + parent_len = (sites->ancestral_state_offset[s + 1] + - sites->ancestral_state_offset[s]); + parent_allele + = sites->ancestral_state + sites->ancestral_state_offset[s]; + } else { + parent_len = (mutations->derived_state_offset[k + 1] + - mutations->derived_state_offset[k]); + parent_allele + = mutations->derived_state + mutations->derived_state_offset[k]; + } + CU_ASSERT_FATAL(len > parent_len); + if (parent_len > 0) { + CU_ASSERT_EQUAL_FATAL(memcmp(ds, parent_allele, parent_len), 0); + CU_ASSERT_EQUAL_FATAL((ds + parent_len)[0], 44); // 44 is ',' in ascii + parent_len++; + } + CU_ASSERT_FATAL(len > parent_len); + mut_id = parse_text_int64(ds + parent_len, len - parent_len); + all_mut_ids[j] = mut_id; + } + } +} + +static void +verify_slim_v6_metadata(tsk_site_table_t *sites, tsk_mutation_table_t *mutations, + int64_t *all_mut_ids, size_t all_mut_ids_length, int skip_mut) +{ + int j, k, s; + size_t len, parent_len; + int64_t mut_id; + char *md, *pmd; + CU_ASSERT_FATAL(mutations->num_rows <= all_mut_ids_length); + for (j = 0; j < mutations->num_rows; j++) { + md = mutations->metadata + mutations->metadata_offset[j]; + len = (mutations->metadata_offset[j + 1] - mutations->metadata_offset[j]); + if (j == skip_mut) { + parent_len = 0; + all_mut_ids[j] = -1; + } else { + k = mutations->parent[j]; + if (k == TSK_NULL) { + s = mutations->site[j]; + pmd = sites->metadata + sites->metadata_offset[s]; + parent_len = (sites->metadata_offset[s + 1] - sites->metadata_offset[s]); + } else { + pmd = mutations->metadata + mutations->metadata_offset[k]; + parent_len = (mutations->metadata_offset[k + 1] + - mutations->metadata_offset[k]); + } + CU_ASSERT_FATAL(len > parent_len); + if (parent_len > 0) { + CU_ASSERT_EQUAL_FATAL(memcmp(md, pmd, parent_len), 0); + } + CU_ASSERT_FATAL(len == parent_len + sizeof(int64_t)); + mut_id = (int64_t) md[parent_len]; + all_mut_ids[j] = mut_id; + } + } +} + +static void +test_mutgen_slim_v6_mutations(void) +{ + int ret = 0; + int j; + mutgen_t mutgen; + gsl_rng *rng = gsl_rng_alloc(gsl_rng_default); + tsk_table_collection_t tables; + mutation_model_t mut_model; + int64_t *all_mut_ids; + int64_t next_mutation_id = 23; + + CU_ASSERT_FATAL(rng != NULL); + ret = slim_v6_mutation_model_alloc(&mut_model, next_mutation_id, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + + ret = tsk_table_collection_init(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + insert_single_tree(&tables, -1); + + ret = mutgen_alloc(&mutgen, rng, &tables, &mut_model, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = mutgen_set_rate(&mutgen, 1); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = mutgen_generate(&mutgen, MSP_DISCRETE_SITES); + CU_ASSERT_EQUAL_FATAL(ret, 0); + CU_ASSERT_FATAL(tables.mutations.num_rows > 0); + + // should have empty ancestral states + for (j = 0; j < tables.sites.num_rows; j++) { + CU_ASSERT_EQUAL_FATAL(tables.sites.ancestral_state_offset[j], 0); + } + // check that derived states append unique integers, + // counting up from next_mutation_id + // in both derived state (text) and metadata (binary) + all_mut_ids = malloc(tables.mutations.num_rows * sizeof(int64_t)); + CU_ASSERT_FATAL(all_mut_ids != NULL); + // derived state + verify_slim_v6_derived_states( + &tables.sites, &tables.mutations, all_mut_ids, tables.mutations.num_rows, -1); + verify_slim_mutation_ids(all_mut_ids, tables.mutations.num_rows, next_mutation_id); + // metadata + verify_slim_v6_metadata( + &tables.sites, &tables.mutations, all_mut_ids, tables.mutations.num_rows, -1); + verify_slim_mutation_ids(all_mut_ids, tables.mutations.num_rows, next_mutation_id); + + mutgen_print_state(&mutgen, _devnull); + + mutgen_free(&mutgen); + free(all_mut_ids); + mutation_model_free(&mut_model); + tsk_table_collection_free(&tables); + gsl_rng_free(rng); +} + +static void +test_mutgen_slim_v6_mutations_keeps(void) +{ + int ret = 0; + int j; + mutgen_t mutgen; + gsl_rng *rng = gsl_rng_alloc(gsl_rng_default); + tsk_table_collection_t tables, orig_tables; + mutation_model_t mut_model; + int64_t *all_mut_ids; + int64_t next_mutation_id = 23; + + CU_ASSERT_FATAL(rng != NULL); + ret = slim_v6_mutation_model_alloc(&mut_model, next_mutation_id, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + + ret = tsk_table_collection_init(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + insert_single_tree(&tables, ALPHABET_NUCLEOTIDE); + // this should have a site and a mutation and only one position + // so we know future mutations will hit it + CU_ASSERT_EQUAL_FATAL(tables.sites.num_rows, 1); + // change the position of that site to 0.0 + tables.sites.position[0] = 0.0; + CU_ASSERT_FATAL(tables.mutations.num_rows > 0); + CU_ASSERT_EQUAL_FATAL(tables.sequence_length, 1.0); + + ret = tsk_table_collection_copy(&tables, &orig_tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + CU_ASSERT_EQUAL_FATAL(tables.sequence_length, 1.0); + + ret = mutgen_alloc(&mutgen, rng, &tables, &mut_model, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = mutgen_set_rate(&mutgen, 1); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = mutgen_generate(&mutgen, MSP_DISCRETE_SITES | MSP_KEEP_SITES); + CU_ASSERT_EQUAL_FATAL(ret, 0); + CU_ASSERT_FATAL(tables.mutations.num_rows > orig_tables.mutations.num_rows); + CU_ASSERT_EQUAL_FATAL(tables.sites.num_rows, 1); + + // check site metadata has been preserved + ret = tsk_site_table_equals(&tables.sites, &orig_tables.sites, 0); + CU_ASSERT_FATAL(ret); + + // check for derived states and metadata that each mutation + // is something appended onto its parent's + all_mut_ids = malloc(tables.mutations.num_rows * sizeof(int64_t)); + CU_ASSERT_FATAL(all_mut_ids != NULL); + for (j = 0; j < tables.mutations.num_rows; j++) { + // find the mutation that's kept; it doesn't follow the rules + if (tables.mutations.time[j] == orig_tables.mutations.time[0]) + break; + } + verify_slim_v6_derived_states( + &tables.sites, &tables.mutations, all_mut_ids, tables.mutations.num_rows, j); + verify_slim_mutation_ids(all_mut_ids, tables.mutations.num_rows, next_mutation_id); + verify_slim_v6_metadata( + &tables.sites, &tables.mutations, all_mut_ids, tables.mutations.num_rows, j); + verify_slim_mutation_ids(all_mut_ids, tables.mutations.num_rows, next_mutation_id); + + free(all_mut_ids); + mutgen_free(&mutgen); + mutation_model_free(&mut_model); + tsk_table_collection_free(&tables); + tsk_table_collection_free(&orig_tables); + gsl_rng_free(rng); +} + +static void +test_mutgen_slim_v6_mutation_large_values(void) +{ + int ret = 0; + mutgen_t mutgen; + tsk_mutation_t mut; + gsl_rng *rng = gsl_rng_alloc(gsl_rng_default); + tsk_table_collection_t tables; + mutation_model_t mut_model; + char value[21]; /* longest 64 bit int has 20 digits */ + int64_t next_mutation_id; + + CU_ASSERT_FATAL(rng != NULL); + ret = tsk_table_collection_init(&tables, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + insert_single_tree(&tables, -1); + + /* Trying to generate mutations that overflow raises an error */ + ret = slim_v6_mutation_model_alloc(&mut_model, INT64_MAX, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = mutgen_alloc(&mutgen, rng, &tables, &mut_model, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = mutgen_set_rate(&mutgen, 2); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = mutgen_generate(&mutgen, MSP_DISCRETE_SITES); + CU_ASSERT_EQUAL_FATAL(ret, MSP_ERR_MUTATION_ID_OVERFLOW); + mutgen_free(&mutgen); + mutation_model_free(&mut_model); + + tsk_mutation_table_clear(&tables.mutations); + tsk_site_table_clear(&tables.sites); + /* Try out with a large value that doesn't hit the ceiling */ + next_mutation_id = INT64_MAX - 100; + ret = slim_v6_mutation_model_alloc(&mut_model, next_mutation_id, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = mutgen_alloc(&mutgen, rng, &tables, &mut_model, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = mutgen_set_rate(&mutgen, 2); + CU_ASSERT_EQUAL_FATAL(ret, 0); + ret = mutgen_generate(&mutgen, MSP_DISCRETE_SITES); + CU_ASSERT_EQUAL_FATAL(ret, 0); + CU_ASSERT_TRUE(tables.mutations.num_rows > 10); + mutgen_free(&mutgen); + mutation_model_free(&mut_model); + + ret = tsk_mutation_table_get_row(&tables.mutations, 0, &mut); + CU_ASSERT_EQUAL_FATAL(ret, 0); + CU_ASSERT_EQUAL(mut.derived_state_length, 19); + sprintf(value, "%" PRId64, next_mutation_id); + CU_ASSERT_NSTRING_EQUAL(value, mut.derived_state, mut.derived_state_length); + + tsk_table_collection_free(&tables); + gsl_rng_free(rng); +} + static void test_mutgen_infinite_alleles(void) { @@ -1306,6 +1574,33 @@ test_slim_mutation_model_properties(void) mutation_model_free(&model); } +static void +test_slim_v6_mutation_model_errors(void) +{ + int ret; + mutation_model_t model; + int64_t next_mutation_id = 0; + + next_mutation_id--; + ret = slim_v6_mutation_model_alloc(&model, next_mutation_id, 0); + CU_ASSERT_EQUAL_FATAL(ret, MSP_ERR_BAD_SLIM_PARAMETERS); + mutation_model_free(&model); +} + +static void +test_slim_v6_mutation_model_properties(void) +{ + int ret; + mutation_model_t model; + int64_t next_mutation_id = 2; + + ret = slim_v6_mutation_model_alloc(&model, next_mutation_id, 0); + CU_ASSERT_EQUAL_FATAL(ret, 0); + CU_ASSERT_EQUAL_FATAL(model.params.slim_mutator.next_mutation_id, 2); + + mutation_model_free(&model); +} + int main(int argc, char **argv) { @@ -1331,6 +1626,10 @@ main(int argc, char **argv) { "test_mutgen_slim_mutations", test_mutgen_slim_mutations }, { "test_mutgen_slim_mutation_large_values", test_mutgen_slim_mutation_large_values }, + { "test_mutgen_slim_v6_mutations", test_mutgen_slim_v6_mutations }, + { "test_mutgen_slim_v6_mutations_keeps", test_mutgen_slim_v6_mutations_keeps }, + { "test_mutgen_slim_v6_mutation_large_values", + test_mutgen_slim_v6_mutation_large_values }, { "test_mutgen_infinite_alleles", test_mutgen_infinite_alleles }, { "test_mutgen_infinite_alleles_large_values", test_mutgen_infinite_alleles_large_values }, @@ -1339,6 +1638,9 @@ main(int argc, char **argv) test_matrix_mutation_model_properties }, { "test_slim_mutation_model_errors", test_slim_mutation_model_errors }, { "test_slim_mutation_model_properties", test_slim_mutation_model_properties }, + { "test_slim_v6_mutation_model_errors", test_slim_v6_mutation_model_errors }, + { "test_slim_v6_mutation_model_properties", + test_slim_v6_mutation_model_properties }, CU_TEST_INFO_NULL, }; diff --git a/msprime/__init__.py b/msprime/__init__.py index 11acef934..1878d7661 100644 --- a/msprime/__init__.py +++ b/msprime/__init__.py @@ -81,6 +81,7 @@ NUCLEOTIDES, PAM, SLiMMutationModel, + SLiMv6MutationModel, sim_mutations, ) @@ -151,6 +152,7 @@ "RateMap", "RecombinationMap", "SLiMMutationModel", + "SLiMv6MutationModel", "Sample", "SampleSet", "SimpleBottleneck", diff --git a/msprime/_msprimemodule.c b/msprime/_msprimemodule.c index f3cf82d94..e9b38d673 100644 --- a/msprime/_msprimemodule.c +++ b/msprime/_msprimemodule.c @@ -65,6 +65,11 @@ typedef struct { mutation_model_t *mutation_model; } SLiMMutationModel; +typedef struct { + PyObject_HEAD + mutation_model_t *mutation_model; +} SLiMv6MutationModel; + typedef struct { PyObject_HEAD mutation_model_t *mutation_model; @@ -796,6 +801,98 @@ static PyTypeObject SLiMMutationModelType = { .tp_new = PyType_GenericNew, }; +/*=================================================================== + * SLiM mutation model, v6 + *=================================================================== + */ + +static int +SLiMv6MutationModel_check_state(SLiMv6MutationModel *self) +{ + int ret = 0; + if (self->mutation_model == NULL) { + PyErr_SetString(PyExc_SystemError, "SLiMv6MutationModel not initialised"); + ret = -1; + } + return ret; +} + +static void +SLiMv6MutationModel_dealloc(SLiMv6MutationModel *self) +{ + if (self->mutation_model != NULL) { + mutation_model_free(self->mutation_model); + PyMem_Free(self->mutation_model); + self->mutation_model = NULL; + } + Py_TYPE(self)->tp_free((PyObject *) self); +} + +static int +SLiMv6MutationModel_init(SLiMv6MutationModel *self, PyObject *args, PyObject *kwds) +{ + int ret = -1; + int err; + static char *kwlist[] = { "next_id", "block_size", NULL }; + long long next_id = 0; + Py_ssize_t block_size = 0; + + self->mutation_model = NULL; + if (!PyArg_ParseTupleAndKeywords(args, kwds, "|Ln", kwlist, &next_id, &block_size)) { + goto out; + } + + /* Note: it's important we zero out mutation_model here because + * we can error before we can mutation_model_alloc, leaving the + * object in an uninitialised state */ + self->mutation_model = PyMem_Calloc(1, sizeof(*self->mutation_model)); + if (self->mutation_model == NULL) { + PyErr_NoMemory(); + goto out; + } + err = slim_v6_mutation_model_alloc( + self->mutation_model, (int64_t) next_id, (size_t) block_size); + if (err != 0) { + handle_library_error(err); + goto out; + } + ret = 0; +out: + return ret; +} + +static PyObject * +SLiMv6MutationModel_get_next_id(SLiMv6MutationModel *self, void *closure) +{ + slim_mutator_t *params; + PyObject *ret = NULL; + + if (SLiMv6MutationModel_check_state(self) != 0) { + goto out; + } + params = &self->mutation_model->params.slim_mutator; + ret = Py_BuildValue("L", params->next_mutation_id); +out: + return ret; +} + +static PyGetSetDef SLiMv6MutationModel_getsetters[] = { + { "next_id", (getter) SLiMv6MutationModel_get_next_id, NULL, + "Return the next mutation id" }, + { NULL } /* Sentinel */ +}; + +static PyTypeObject SLiMv6MutationModelType = { + .tp_name = "_msprime.SLiMv6MutationModel", + .tp_basicsize = sizeof(SLiMv6MutationModel), + .tp_dealloc = (destructor) SLiMv6MutationModel_dealloc, + .tp_flags = Py_TPFLAGS_DEFAULT | Py_TPFLAGS_BASETYPE, + .tp_doc = "SLiMv6MutationModel objects", + .tp_getset = SLiMv6MutationModel_getsetters, + .tp_init = (initproc) SLiMv6MutationModel_init, + .tp_new = PyType_GenericNew, +}; + /*=================================================================== * Infinite alleles mutation model *=================================================================== @@ -3014,6 +3111,7 @@ parse_mutation_model(PyObject *py_model) mutation_model_t *model = NULL; MatrixMutationModel *matrix_mutation_model = NULL; SLiMMutationModel *slim_mutation_model = NULL; + SLiMv6MutationModel *slim_v6_mutation_model = NULL; InfiniteAllelesMutationModel *infinite_alleles_model = NULL; if (PyObject_TypeCheck(py_model, &MatrixMutationModelType)) { @@ -3028,6 +3126,12 @@ parse_mutation_model(PyObject *py_model) goto out; } model = slim_mutation_model->mutation_model; + } else if (PyObject_TypeCheck(py_model, &SLiMv6MutationModelType)) { + slim_v6_mutation_model = (SLiMv6MutationModel *) py_model; + if (SLiMv6MutationModel_check_state(slim_v6_mutation_model) != 0) { + goto out; + } + model = slim_v6_mutation_model->mutation_model; } else if (PyObject_TypeCheck(py_model, &InfiniteAllelesMutationModelType)) { infinite_alleles_model = (InfiniteAllelesMutationModel *) py_model; if (InfiniteAllelesMutationModel_check_state(infinite_alleles_model) != 0) { @@ -3037,7 +3141,7 @@ parse_mutation_model(PyObject *py_model) } else { PyErr_SetString(PyExc_TypeError, "model must be an instance of MatrixMutationModel, " - "SLiMMutationModel or InfiniteAllelesMutationModel."); + "SLiMMutationModel, SLiMv6MutationModel or InfiniteAllelesMutationModel."); goto out; } out: @@ -3256,6 +3360,15 @@ PyInit__msprime(void) Py_INCREF(&SLiMMutationModelType); PyModule_AddObject(module, "SLiMMutationModel", (PyObject *) &SLiMMutationModelType); + /* SLiMv6MutationModel type */ + SLiMv6MutationModelType.tp_base = &BaseMutationModelType; + if (PyType_Ready(&SLiMv6MutationModelType) < 0) { + return NULL; + } + Py_INCREF(&SLiMv6MutationModelType); + PyModule_AddObject( + module, "SLiMv6MutationModel", (PyObject *) &SLiMv6MutationModelType); + /* InfiniteAllelesMutationModel type */ InfiniteAllelesMutationModelType.tp_base = &BaseMutationModelType; if (PyType_Ready(&InfiniteAllelesMutationModelType) < 0) { diff --git a/msprime/mutations.py b/msprime/mutations.py index 55c58305d..4d7fbdfa8 100644 --- a/msprime/mutations.py +++ b/msprime/mutations.py @@ -108,7 +108,9 @@ def __str__(self): class SLiMMutationModel(_msprime.SLiMMutationModel, MutationModel): """ - An infinite-alleles model of mutation producing "SLiM-style" mutations. + An infinite-alleles model of mutation producing "SLiM-style" mutations + for versions of SLiM before 6.0. For SLiM version 6.0 and beyond, + use {class}`.SLiMv6MutationModel`. To agree with mutations produced by SLiM, the ancestral state of each new site is set to the empty string, and each derived state is produced by @@ -144,6 +146,48 @@ def __str__(self): ) +class SLiMv6MutationModel(_msprime.SLiMv6MutationModel, MutationModel): + """ + An infinite-alleles model of mutation producing "SLiM-style" mutations + for versions of SLiM 6.0 and beyond. For older versions of SLiM, + use {class}`.SLiMMutationModel`. + + To agree with mutations produced by SLiM, the ancestral state of each new + site is set to the empty string, and each derived state is produced by + appending the "next allele" to the previous state. The result is that each + allele is a comma-separated string of all mutations that have occurred up + to the root. Alleles are numeric IDs, starting with ``next_id``. + + The same list of numeric IDs is stored in metadata, in binary, as int64s. + It is actually these IDs in metadata that SLiM uses when reading in a tree + sequence, not the text-based version in the derived state. + + This model works as following: for each mutation, if the next ID is ``n``, + then it sets the derived state to be the string ``"n"`` appended to its + inherited state (which is the derived state of the parent mutation, if any, + or the ancestral state otherwise), with an intervening ``","`` if the inherited + state is nonempty. Also, it sets the metadata to be the 64-bit integer n + appended to the bytes of the inherited metadata (which similarly comes from + the parent mutation, if any, or the site). So, if this mutation model is + applied to a tree sequence with existing mutations of a different sort (especially + if the existing sites have nonempty metadata) the result may be surprising: + see {ref}`sec_mutations_mutation_slim_mutations_stacking`. + + To update your code to use this model instead of {class}`.SLiMMutationModel`, + remove the ``type`` and ``slim_generation`` arguments, and follow + up {func}`.sim_mutations` with a call to ``pyslim.add_mutation_metadata``. + + :param int next_id: The nonnegative integer to start assigning alleles from. + (default: 0) + :param int block_size: The block size for allocating derived states. + You do not need to change this unless you get an "out of memory" error + due to a very large number of stacked mutations. + """ + + def __str__(self): + return f"Mutation model for SLiM mutations, v6+. Next ID: {self.next_id}\n" + + # NOTE: we use a hacky workaround here for documenting the next_allele instance # variable because of the nasty sphinx bug: # https://github.com/sphinx-doc/sphinx/issues/2549 @@ -1512,7 +1556,8 @@ def mutate( "jc69": JC69, "blosum62": BLOSUM62, "pam": PAM, - # "slim": SLiMMutationModel(), Needs type argument so can't be string init'd + # "slim": SLiMMutationModel(), or SLiMv6MutationModel(), not clear + # and maybe needs the next_id argument # "hky": HKY(), Needs kappa argument # "f84": F84(), Needs kappa argument # "gtr": GTR(), Needs relative_rates argument diff --git a/tests/test_mutations.py b/tests/test_mutations.py index 569b5d9a2..1deae6608 100644 --- a/tests/test_mutations.py +++ b/tests/test_mutations.py @@ -1551,7 +1551,57 @@ class SlimMetadata: nucleotide: int -class TestSLiMMutationModel: +class SLiMModelMixin: + """ + Common tests for both SLiM mutation models + """ + + def test_binary_n_4_low_rate(self): + ts = msprime.sim_ancestry(4, sequence_length=10, random_seed=5) + mts = self.run_mutate(ts, rate=0.1, random_seed=23) + assert mts.num_mutations > 1 + self.validate_slim_mutations(mts) + + def test_binary_n_4_high_rate(self): + ts = msprime.sim_ancestry(4, sequence_length=2, random_seed=5) + mts = self.run_mutate(ts, rate=2.0, random_seed=23) + assert mts.num_mutations > 10 + self.validate_slim_mutations(mts) + + def test_binary_n_8_low_rate(self): + ts = msprime.sim_ancestry(8, sequence_length=10, random_seed=50) + mts = self.run_mutate(ts, rate=0.1, random_seed=342) + assert mts.num_mutations > 1 + self.validate_slim_mutations(mts) + + def test_binary_n_8_high_rate(self): + ts = msprime.sim_ancestry(8, sequence_length=10, random_seed=5) + mts = self.run_mutate(ts, rate=2.0, random_seed=23) + assert mts.num_mutations > 10 + self.validate_slim_mutations(mts) + + def test_binary_incomplete_trees(self): + ts = msprime.sim_ancestry(8, sequence_length=5, random_seed=50, end_time=0.1) + assert ts.first().num_roots > 1 + mts = self.run_mutate(ts, rate=2.0, random_seed=23) + assert mts.num_mutations > 10 + self.validate_slim_mutations(mts) + + def test_binary_many_trees(self): + ts = msprime.sim_ancestry( + 8, + sequence_length=5, + recombination_rate=5, + random_seed=50, + discrete_genome=False, + ) + assert ts.num_trees > 20 + mts = self.run_mutate(ts, rate=2.0, random_seed=23) + assert mts.num_mutations > 10 + self.validate_slim_mutations(mts) + + +class TestSLiMMutationModel(SLiMModelMixin): """ Tests for the SLiM mutation generator. """ @@ -1619,18 +1669,10 @@ def run_mutate( t1 = mts1.dump_tables() t2 = mts2.dump_tables() - assert t1.sites == t2.sites # Drop the mutation metadata - we're validating that elsewhere and # it's not worth complicating the Python generator with it. - t1.mutations.set_columns( - site=t1.mutations.site, - node=t1.mutations.node, - parent=t1.mutations.parent, - time=t1.mutations.time, - derived_state=t1.mutations.derived_state, - derived_state_offset=t1.mutations.derived_state_offset, - ) - assert t1.mutations == t2.mutations + t1.sites.assert_equals(t2.sites) + t1.mutations.assert_equals(t2.mutations, ignore_metadata=True) return mts1 def test_slim_mutation_type(self): @@ -1651,49 +1693,194 @@ def test_slim_generation(self): assert mts.num_mutations > 10 self.validate_slim_mutations(mts, slim_generation=slim_generation) - def test_binary_n_4_low_rate(self): - ts = msprime.sim_ancestry(4, sequence_length=10, random_seed=5) - mts = self.run_mutate(ts, rate=0.1, random_seed=23) - assert mts.num_mutations > 1 - self.validate_slim_mutations(mts) - def test_binary_n_4_high_rate(self): - ts = msprime.sim_ancestry(4, sequence_length=2, random_seed=5) - mts = self.run_mutate(ts, rate=2.0, random_seed=23) - assert mts.num_mutations > 10 - self.validate_slim_mutations(mts) +class TestSLiMv6MutationModel(SLiMModelMixin): + """ + Tests for the SLiMv6 mutation generator. + """ - def test_binary_n_8_low_rate(self): - ts = msprime.sim_ancestry(8, sequence_length=10, random_seed=50) - mts = self.run_mutate(ts, rate=0.1, random_seed=342) - assert mts.num_mutations > 1 - self.validate_slim_mutations(mts) + metadata_schema = tskit.MetadataSchema( + { + "codec": "struct", + "type": "object", + "properties": { + "derived_states": { + "items": { + "binaryFormat": "q", + "type": "number", + }, + "noLengthEncodingExhaustBuffer": True, + "type": "array", + } + }, + "required": ["derived_states"], + "additionalProperties": False, + } + ) - def test_binary_n_8_high_rate(self): - ts = msprime.sim_ancestry(8, sequence_length=10, random_seed=5) - mts = self.run_mutate(ts, rate=2.0, random_seed=23) - assert mts.num_mutations > 10 - self.validate_slim_mutations(mts) + def validate_slim_mutations(self, ts): + # slim alleles should be lists of integers + # and ancestral states the empty string + # while metadata should be lists of integers + t = ts.dump_tables() + t.mutations.metadata_schema = self.metadata_schema + ts = t.tree_sequence() + for site in ts.sites(): + assert site.ancestral_state == "" + alleles = {} + for mutation in site.mutations: + a = list(map(int, mutation.derived_state.split(","))) + md = mutation.metadata["derived_states"] + assert len(a) == len(md) + for x, y in zip(a, md): + assert x == y + alleles[mutation.id] = a + if mutation.parent == tskit.NULL: + assert len(a) == 1 + else: + parent_allele = alleles[mutation.parent] + assert a[:-1] == parent_allele - def test_binary_incomplete_trees(self): - ts = msprime.sim_ancestry(8, sequence_length=5, random_seed=50, end_time=0.1) - assert ts.first().num_roots > 1 - mts = self.run_mutate(ts, rate=2.0, random_seed=23) - assert mts.num_mutations > 10 - self.validate_slim_mutations(mts) + def run_mutate( + self, + ts, + rate=1, + random_seed=42, + mutation_id=0, + keep=True, + ): + model = msprime.SLiMv6MutationModel(next_id=mutation_id) + mts1 = msprime.sim_mutations( + ts, + rate=rate, + random_seed=random_seed, + model=model, + discrete_genome=True, + keep=keep, + ) + if ts.num_mutations == 0 or not keep: + assert mts1.num_mutations == model.next_id - def test_binary_many_trees(self): - ts = msprime.sim_ancestry( - 8, - sequence_length=5, - recombination_rate=5, - random_seed=50, - discrete_genome=False, + model = PythonSLiMv6MutationModel(next_id=mutation_id) + mts2 = py_sim_mutations( + ts, + rate=rate, + random_seed=random_seed, + model=model, + discrete_genome=True, + keep=keep, ) - assert ts.num_trees > 20 + + t1 = mts1.dump_tables() + t2 = mts2.dump_tables() + # Drop the mutation metadata - we're validating that elsewhere and + # it's not worth complicating the Python generator with it. + t1.sites.assert_equals(t2.sites) + t1.mutations.assert_equals(t2.mutations, ignore_metadata=True) + return mts1 + + metadata_schema_with_extra = tskit.MetadataSchema( + { + "codec": "struct", + "type": "object", + "properties": { + "tag": { + "index": 0, + "type": "string", + "binaryFormat": "24s", + "nullTerminated": True, + }, + "derived_states": { + "index": 1, + "items": { + "binaryFormat": "q", + "type": "number", + }, + "noLengthEncodingExhaustBuffer": True, + "type": "array", + }, + }, + "required": ["tag", "derived_states"], + "additionalProperties": False, + } + ) + + def test_keep(self): + ts = msprime.sim_ancestry(8, sequence_length=10, random_seed=5) mts = self.run_mutate(ts, rate=2.0, random_seed=23) - assert mts.num_mutations > 10 - self.validate_slim_mutations(mts) + mt = mts.dump_tables() + # there can be weird interactions with site metadata, so add some + mt.sites.clear() + mt.sites.metadata_schema = self.metadata_schema_with_extra + for s in mts.sites(): + md = { + "tag": f"lmnop{s.id}", + "derived_states": [], + } + mt.sites.append(s.replace(metadata=md)) + mt.mutations.clear() + mt.mutations.metadata_schema = self.metadata_schema_with_extra + for m in mts.mutations(): + md = { + "tag": f"lmnop{m.site}", + "derived_states": [int(x) for x in m.derived_state.split(",")], + } + mt.mutations.append(m.replace(metadata=md)) + mts = mt.tree_sequence() + num_muts = {s.position: len(s.mutations) for s in mts.sites()} + mut_id_start = mts.num_mutations + 1 + model = msprime.SLiMv6MutationModel(next_id=mut_id_start) + mmts = msprime.sim_mutations( + mts, + rate=2.0, + random_seed=1234, + model=model, + discrete_genome=True, + keep=True, + ) + # check we've stacked on top of some previous ones + more_hits = False + for s in mmts.sites(): + if s.position in num_muts: + if len(s.mutations) > num_muts[s.position]: + more_hits = True + break + assert more_hits + for mut in mmts.mutations(): + site = mmts.site(mut.site) + mp = mut.parent + # mutations we put down the first time won't have inherited the state + # from mutations we put down the second time, so traverse up until + # we find the correct one we should have inherited from + mut_id = mut.metadata["derived_states"][-1] + if mut_id < mut_id_start: + while ( + mp >= 0 + and mmts.mutation(mp).metadata["derived_states"][-1] >= mut_id_start + ): + mp = mmts.mutation(mp).parent + if mp == -1: + parent_allele = site.ancestral_state + parent_metadata = site.metadata + else: + parent_allele = mmts.mutation(mp).derived_state + parent_metadata = mmts.mutation(mp).metadata + n = len(parent_allele) + assert n < len(mut.derived_state) + assert parent_allele == mut.derived_state[:n] + if n > 0: + assert mut.derived_state[n] == "," + n += 1 + assert mut_id == int(mut.derived_state[n:]) + assert parent_metadata["tag"] == mut.metadata["tag"] + n = len(parent_metadata["derived_states"]) + assert n + 1 == len(mut.metadata["derived_states"]) + if n > 0: + assert ( + parent_metadata["derived_states"] + == mut.metadata["derived_states"][:n] + ) + assert mut.metadata["derived_states"][-1] == mut_id class TestInfiniteAllelesMutationModel: @@ -1984,8 +2171,7 @@ def transition_allele(self, rng, current_allele): @dataclasses.dataclass -class PythonSLiMMutationModel(PythonMutationModel): - mutation_type: int = 0 +class PythonSLiMv6MutationModel(PythonMutationModel): next_id: int = 0 def root_allele(self, rng): @@ -2000,6 +2186,11 @@ def transition_allele(self, rng, current_allele): return out +@dataclasses.dataclass +class PythonSLiMMutationModel(PythonSLiMv6MutationModel): + mutation_type: int = 0 + + @dataclasses.dataclass class PythonInfiniteAllelesMutationModel(PythonMutationModel): start_allele: int = 0 @@ -2302,6 +2493,7 @@ def test_bad_models(self): def test_returns_mutation_model_instances_without_copying(self): models = [ msprime.SLiMMutationModel(0, 0), + msprime.SLiMv6MutationModel(0), msprime.InfiniteAlleles(), msprime.BinaryMutationModel(), msprime.JC69(), @@ -2324,6 +2516,11 @@ def test_slim(self): assert m.slim_generation == 9 assert str(m) == "Mutation model for SLiM mutations of type m1\n next ID: 2\n" + def test_slim_v6(self): + m = msprime.SLiMv6MutationModel(next_id=2) + assert m.next_id == 2 + assert str(m) == "Mutation model for SLiM mutations, v6+. Next ID: 2\n" + def test_infinite_alleles(self): m = msprime.InfiniteAlleles(start_allele=1) assert ( diff --git a/tests/test_python_c.py b/tests/test_python_c.py index 5f12e5e7e..105b35c09 100644 --- a/tests/test_python_c.py +++ b/tests/test_python_c.py @@ -2808,6 +2808,30 @@ def test_generation(self): assert model.slim_generation == 1 +class TestSLiMv6MutationModel: + """ + Tests for the slim v6 mutation model class. + """ + + def test_constructor_errors(self): + for bad_id in ["sdr", 0.222, None]: + with pytest.raises(TypeError): + _msprime.SLiMv6MutationModel(next_id=bad_id) + + with pytest.raises(_msprime.LibraryError): + _msprime.SLiMv6MutationModel(next_id=-1) + + def test_uninitialised(self): + model = _msprime.SLiMv6MutationModel.__new__(_msprime.SLiMv6MutationModel) + with pytest.raises(SystemError): + _ = model.next_id + + def test_next_id(self): + for next_id in [0, 10, 2**63 - 1]: + model = _msprime.SLiMv6MutationModel(next_id=next_id) + assert model.next_id == next_id + + class TestInfiniteAllelesMutationModel: """ Tests for the infinite alleles mutation model class. @@ -2984,6 +3008,7 @@ def test_model(self): _msprime.sim_mutations(tables, rng, rate_map=imap, model=bad_type) model_classes = [ _msprime.SLiMMutationModel, + _msprime.SLiMv6MutationModel, _msprime.InfiniteAllelesMutationModel, _msprime.MatrixMutationModel, ] @@ -3002,6 +3027,14 @@ def test_slim_model(self): _msprime.sim_mutations(tables, rng, imap, model) assert model.next_id == 5678 + def test_slim_v6_model(self): + rng = _msprime.RandomGenerator(1) + imap = uniform_rate_map(1) + tables = _msprime.LightweightTableCollection(1.0) + model = _msprime.SLiMv6MutationModel(5678) + _msprime.sim_mutations(tables, rng, imap, model) + assert model.next_id == 5678 + def test_infinite_alleles_model(self): rng = _msprime.RandomGenerator(1) imap = uniform_rate_map(1, 1)