Skip to content

Support multiple alleles and mutations per site #1200

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jeromekelleher merged 3 commits into from
Oct 16, 2026
Merged

Support multiple alleles and mutations per site #1200

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
493ca23
Support multiple alleles and mutations per site in HMM matcher
jeromekelleher Oct 16, 2026
b877dca
Pass real allele codes of panel mutations to the matcher
jeromekelleher Oct 16, 2026
664d4ee
Require 0/1 genotypes for ancestor generation
jeromekelleher Oct 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
53 changes: 48 additions & 5 deletions _tsinfermodule.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -426,13 +426,21 @@ MatcherIndexes_init(MatcherIndexes *self, PyObject *args, PyObject *kwds)
int err;
LightweightTableCollection *tables;
PyObject *num_alleles_obj = NULL;
PyObject *ancestral_state_obj = NULL;
PyObject *derived_state_obj = NULL;
PyArrayObject *num_alleles_array = NULL;
PyArrayObject *ancestral_state_array = NULL;
PyArrayObject *derived_state_array = NULL;
tsk_size_t *num_alleles_data = NULL;
static char *kwlist[] = { "tables", "num_alleles", NULL };
allele_t *ancestral_state_data = NULL;
allele_t *derived_state_data = NULL;
static char *kwlist[]
= { "tables", "num_alleles", "ancestral_state", "derived_state", NULL };

self->matcher_indexes = NULL;
if (!PyArg_ParseTupleAndKeywords(args, kwds, "O!|O", kwlist,
&LightweightTableCollectionType, &tables, &num_alleles_obj)) {
if (!PyArg_ParseTupleAndKeywords(args, kwds, "O!|OOO", kwlist,
&LightweightTableCollectionType, &tables, &num_alleles_obj,
&ancestral_state_obj, &derived_state_obj)) {
goto out;
}
if (LightweightTableCollection_check_state(tables) != 0) {
Expand All @@ -444,23 +452,58 @@ MatcherIndexes_init(MatcherIndexes *self, PyObject *args, PyObject *kwds)
if (num_alleles_array == NULL) {
goto out;
}
if (PyArray_DIM(num_alleles_array, 0)
!= (npy_intp) tables->tables->sites.num_rows) {
PyErr_SetString(PyExc_ValueError, "num_alleles must have length num_sites");
goto out;
}
num_alleles_data = PyArray_DATA(num_alleles_array);
}
if (ancestral_state_obj != NULL && ancestral_state_obj != Py_None) {
ancestral_state_array = (PyArrayObject *) PyArray_FROMANY(
ancestral_state_obj, NPY_INT8, 1, 1, NPY_ARRAY_IN_ARRAY);
if (ancestral_state_array == NULL) {
goto out;
}
if (PyArray_DIM(ancestral_state_array, 0)
!= (npy_intp) tables->tables->sites.num_rows) {
PyErr_SetString(
PyExc_ValueError, "ancestral_state must have length num_sites");
goto out;
}
ancestral_state_data = PyArray_DATA(ancestral_state_array);
}
if (derived_state_obj != NULL && derived_state_obj != Py_None) {
derived_state_array = (PyArrayObject *) PyArray_FROMANY(
derived_state_obj, NPY_INT8, 1, 1, NPY_ARRAY_IN_ARRAY);
if (derived_state_array == NULL) {
goto out;
}
if (PyArray_DIM(derived_state_array, 0)
!= (npy_intp) tables->tables->mutations.num_rows) {
PyErr_SetString(
PyExc_ValueError, "derived_state must have length num_mutations");
goto out;
}
derived_state_data = PyArray_DATA(derived_state_array);
}

self->matcher_indexes = PyMem_Calloc(1, sizeof(*self->matcher_indexes));
if (self->matcher_indexes == NULL) {
PyErr_NoMemory();
goto out;
}
err = matcher_indexes_alloc(
self->matcher_indexes, tables->tables, num_alleles_data, 0);
err = matcher_indexes_alloc(self->matcher_indexes, tables->tables, num_alleles_data,
ancestral_state_data, derived_state_data, 0);
if (err != 0) {
handle_library_error(err);
goto out;
}
ret = 0;
out:
Py_XDECREF(num_alleles_array);
Py_XDECREF(ancestral_state_array);
Py_XDECREF(derived_state_array);
return ret;
}

Expand Down
43 changes: 29 additions & 14 deletions lib/ancestor_matcher.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -170,13 +170,21 @@ matcher_indexes_add_mutation(
return ret;
}

/* Copy the site and mutation data into the matcher indexes. The
* ancestral_state (per site) and derived_state (per mutation) arrays
* give the integer-encoded alleles; if NULL, the ancestral state is 0
* and every mutation has derived state 1. Mutations must be in tskit
* order, so that when two mutations at a site are on the same node the
* later one (the child) takes precedence. */
static int
matcher_indexes_copy_mutation_data(matcher_indexes_t *self,
const tsk_table_collection_t *tables, const tsk_size_t *num_alleles)
const tsk_table_collection_t *tables, const tsk_size_t *num_alleles,
const allele_t *ancestral_state, const allele_t *derived_state)
{
int ret = 0;
tsk_size_t j;
tsk_id_t site, last_site;
tsk_id_t site;
allele_t state;
const double *restrict sites_position = tables->sites.position;
const tsk_id_t *restrict mutations_site = tables->mutations.site;
const tsk_id_t *restrict mutations_node = tables->mutations.node;
Expand All @@ -188,23 +196,26 @@ matcher_indexes_copy_mutation_data(matcher_indexes_t *self,
converted_position[j] = (coordinate_t) sites_position[j];
self->sites.mutations[j] = NULL;
self->sites.num_alleles[j] = num_alleles != NULL ? num_alleles[j] : 2;
state = ancestral_state != NULL ? ancestral_state[j] : 0;
if (state < 0 || (tsk_size_t) state >= self->sites.num_alleles[j]) {
ret = TSI_ERR_BAD_ANCESTRAL_STATE;
goto out;
}
self->sites.ancestral_state[j] = state;
}
converted_position[j] = (coordinate_t) tables->sequence_length;

last_site = -1;
for (j = 0; j < total_mutations; j++) {
site = mutations_site[j];
if (site == last_site) {
ret = TSI_ERR_MULTIPLE_MUTATIONS_AT_SITE;
state = derived_state != NULL ? derived_state[j] : 1;
if (state < 0 || (tsk_size_t) state >= self->sites.num_alleles[site]) {
ret = TSI_ERR_BAD_MUTATION_DERIVED_STATE;
goto out;
}

self->sites.mutations[site] = NULL;
ret = matcher_indexes_add_mutation(self, site, mutations_node[j], 1);
ret = matcher_indexes_add_mutation(self, site, mutations_node[j], state);
if (ret != 0) {
goto out;
}
last_site = site;
}

out:
Expand Down Expand Up @@ -252,7 +263,8 @@ matcher_indexes_validate_tables(

int
matcher_indexes_alloc(matcher_indexes_t *self, const tsk_table_collection_t *tables,
const tsk_size_t *num_alleles, tsk_flags_t flags)
const tsk_size_t *num_alleles, const allele_t *ancestral_state,
const allele_t *derived_state, tsk_flags_t flags)
{
int ret = 0;

Expand All @@ -273,11 +285,13 @@ matcher_indexes_alloc(matcher_indexes_t *self, const tsk_table_collection_t *tab
self->right_index_edges = malloc(self->num_edges * sizeof(*self->right_index_edges));
self->sites.mutations = malloc(self->num_sites * sizeof(*self->sites.mutations));
self->sites.num_alleles = malloc(self->num_sites * sizeof(*self->sites.num_alleles));
self->sites.ancestral_state
= malloc(self->num_sites * sizeof(*self->sites.ancestral_state));
self->sites.position
= malloc((self->num_sites + 1) * sizeof(*self->sites.mutations));
if (self->left_index_edges == NULL || self->right_index_edges == NULL
|| self->sites.mutations == NULL || self->sites.position == NULL
|| self->sites.num_alleles == NULL) {
|| self->sites.num_alleles == NULL || self->sites.ancestral_state == NULL) {
ret = TSI_ERR_NO_MEMORY;
goto out;
}
Expand All @@ -290,7 +304,8 @@ matcher_indexes_alloc(matcher_indexes_t *self, const tsk_table_collection_t *tab
if (ret != 0) {
goto out;
}
ret = matcher_indexes_copy_mutation_data(self, tables, num_alleles);
ret = matcher_indexes_copy_mutation_data(
self, tables, num_alleles, ancestral_state, derived_state);
if (ret != 0) {
goto out;
}
Expand All @@ -306,6 +321,7 @@ matcher_indexes_free(matcher_indexes_t *self)
tsk_safe_free(self->sites.mutations);
tsk_safe_free(self->sites.position);
tsk_safe_free(self->sites.num_alleles);
tsk_safe_free(self->sites.ancestral_state);
tsi_blkalloc_free(&self->allocator);
return 0;
}
Expand Down Expand Up @@ -550,8 +566,7 @@ ancestor_matcher_set_allelic_state(
{
mutation_list_node_t *mutation;

/* FIXME assuming that 0 is always the ancestral state */
allelic_state[0] = 0;
allelic_state[0] = self->matcher_indexes->sites.ancestral_state[site];

for (mutation = self->matcher_indexes->sites.mutations[site]; mutation != NULL;
mutation = mutation->next) {
Expand Down
4 changes: 0 additions & 4 deletions lib/err.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -112,10 +112,6 @@ tsi_strerror(int err)
case TSI_ERR_BAD_ANCESTRAL_STATE:
ret = "Ancestral state must be >= 0 and < num_alleles";
break;
case TSI_ERR_MULTIPLE_MUTATIONS_AT_SITE:
ret = "Sites with multiple mutations are not supported "
"in matcher_indexes";
break;
case TSI_ERR_BAD_EDGE_NODE:
ret = "Bad edge: parent or child node out of bounds";
break;
Expand Down
1 change: 0 additions & 1 deletion lib/err.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,6 @@
#define TSI_ERR_ONE_BIT_NON_BINARY -24
#define TSI_ERR_IO -25
#define TSI_ERR_BAD_ANCESTRAL_STATE -26
#define TSI_ERR_MULTIPLE_MUTATIONS_AT_SITE -27
#define TSI_ERR_BAD_EDGE_NODE -28
#define TSI_ERR_NODE_0_NOT_ROOT -29
// clang-format on
Expand Down
Loading
Loading