Add 'mutation' mode; switch genetic_relatedness_matrix and

divergence_matrix to this mode; document divergence_matrix

Author: petrelharp

c/tests/test_stats.c

@@ -1315,6 +1315,12 @@ test_general_stat_input_errors(void)
         &ts, 1, &W, 0, general_stat_sum, NULL, 0, NULL, 0, &result);
     CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_BAD_RESULT_DIMS);
 
+    /* Unsupported mode */
+    /* TODO: change when STAT_MUTATION is supported */
+    ret = tsk_treeseq_general_stat(
+        &ts, 1, &W, 1, general_stat_sum, NULL, 0, NULL, TSK_STAT_MUTATION, &result);
+    CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_UNSUPPORTED_STAT_MODE);
+
     /* Multiple stats*/
     ret = tsk_treeseq_general_stat(&ts, 1, &W, 1, general_stat_sum, NULL, 0, NULL,
         TSK_STAT_SITE | TSK_STAT_BRANCH, &result);
@@ -1464,18 +1470,26 @@ test_single_tree_divergence_matrix(void)
     tsk_treeseq_from_text(&ts, 1, single_tree_ex_nodes, single_tree_ex_edges, NULL,
         single_tree_ex_sites, single_tree_ex_mutations, NULL, NULL, 0);
 
+    ret = tsk_treeseq_divergence_matrix(
+        &ts, 0, NULL, NULL, 0, NULL, TSK_STAT_SITE, result);
+    CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_UNSUPPORTED_STAT_MODE);
+
+    ret = tsk_treeseq_divergence_matrix(
+        &ts, 0, NULL, NULL, 0, NULL, TSK_STAT_NODE, result);
+    CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_UNSUPPORTED_STAT_MODE);
+
     ret = tsk_treeseq_divergence_matrix(
         &ts, 0, NULL, NULL, 0, NULL, TSK_STAT_BRANCH, result);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
     assert_arrays_almost_equal(16, result, D_branch);
 
     ret = tsk_treeseq_divergence_matrix(
-        &ts, 0, NULL, NULL, 0, NULL, TSK_STAT_SITE, result);
+        &ts, 0, NULL, NULL, 0, NULL, TSK_STAT_MUTATION, result);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
     assert_arrays_almost_equal(16, result, D_site);
 
     ret = tsk_treeseq_divergence_matrix(
-        &ts, 2, sample_set_sizes, NULL, 0, NULL, TSK_STAT_SITE, result);
+        &ts, 2, sample_set_sizes, NULL, 0, NULL, TSK_STAT_MUTATION, result);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
 
     ret = tsk_treeseq_divergence_matrix(
@@ -1488,15 +1502,15 @@ test_single_tree_divergence_matrix(void)
         &ts, 2, sample_set_sizes, NULL, 0, NULL, TSK_STAT_BRANCH, result);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
     ret = tsk_treeseq_divergence_matrix(
-        &ts, 2, sample_set_sizes, NULL, 0, NULL, TSK_STAT_SITE, result);
+        &ts, 2, sample_set_sizes, NULL, 0, NULL, TSK_STAT_MUTATION, result);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
 
     /* assert_arrays_almost_equal(4, result, D_site); */
 
     verify_divergence_matrix(&ts, TSK_STAT_BRANCH);
     verify_divergence_matrix(&ts, TSK_STAT_BRANCH | TSK_STAT_SPAN_NORMALISE);
-    verify_divergence_matrix(&ts, TSK_STAT_SITE);
-    verify_divergence_matrix(&ts, TSK_STAT_SITE | TSK_STAT_SPAN_NORMALISE);
+    verify_divergence_matrix(&ts, TSK_STAT_MUTATION);
+    verify_divergence_matrix(&ts, TSK_STAT_MUTATION | TSK_STAT_SPAN_NORMALISE);
 
     tsk_treeseq_free(&ts);
 }
@@ -1542,7 +1556,7 @@ test_single_tree_divergence_matrix_internal_samples(void)
     CU_ASSERT_EQUAL_FATAL(ret, 0);
     assert_arrays_almost_equal(16, result, D);
     ret = tsk_treeseq_divergence_matrix(
-        &ts, 0, NULL, NULL, 0, NULL, TSK_STAT_SITE, result);
+        &ts, 0, NULL, NULL, 0, NULL, TSK_STAT_MUTATION, result);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
     assert_arrays_almost_equal(16, result, D);
 
@@ -1551,7 +1565,7 @@ test_single_tree_divergence_matrix_internal_samples(void)
     CU_ASSERT_EQUAL_FATAL(ret, 0);
     assert_arrays_almost_equal(16, result, D);
     ret = tsk_treeseq_divergence_matrix(
-        &ts, 4, sizes, samples, 0, NULL, TSK_STAT_SITE, result);
+        &ts, 4, sizes, samples, 0, NULL, TSK_STAT_MUTATION, result);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
     assert_arrays_almost_equal(16, result, D);
 
@@ -1560,14 +1574,14 @@ test_single_tree_divergence_matrix_internal_samples(void)
     CU_ASSERT_EQUAL_FATAL(ret, 0);
     assert_arrays_almost_equal(16, result, D);
     ret = tsk_treeseq_divergence_matrix(
-        &ts, 4, NULL, samples, 0, NULL, TSK_STAT_SITE, result);
+        &ts, 4, NULL, samples, 0, NULL, TSK_STAT_MUTATION, result);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
     assert_arrays_almost_equal(16, result, D);
 
     verify_divergence_matrix(&ts, TSK_STAT_BRANCH);
     verify_divergence_matrix(&ts, TSK_STAT_BRANCH | TSK_STAT_SPAN_NORMALISE);
-    verify_divergence_matrix(&ts, TSK_STAT_SITE);
-    verify_divergence_matrix(&ts, TSK_STAT_SITE | TSK_STAT_SPAN_NORMALISE);
+    verify_divergence_matrix(&ts, TSK_STAT_MUTATION);
+    verify_divergence_matrix(&ts, TSK_STAT_MUTATION | TSK_STAT_SPAN_NORMALISE);
 
     tsk_treeseq_free(&ts);
     free(result);
@@ -1616,8 +1630,8 @@ test_single_tree_divergence_matrix_multi_root(void)
 
     verify_divergence_matrix(&ts, TSK_STAT_BRANCH);
     verify_divergence_matrix(&ts, TSK_STAT_BRANCH | TSK_STAT_SPAN_NORMALISE);
-    verify_divergence_matrix(&ts, TSK_STAT_SITE);
-    verify_divergence_matrix(&ts, TSK_STAT_SITE | TSK_STAT_SPAN_NORMALISE);
+    verify_divergence_matrix(&ts, TSK_STAT_MUTATION);
+    verify_divergence_matrix(&ts, TSK_STAT_MUTATION | TSK_STAT_SPAN_NORMALISE);
 
     tsk_treeseq_free(&ts);
 }
@@ -2557,8 +2571,8 @@ test_paper_ex_divergence_matrix(void)
 
     verify_divergence_matrix(&ts, TSK_STAT_BRANCH);
     verify_divergence_matrix(&ts, TSK_STAT_BRANCH | TSK_STAT_SPAN_NORMALISE);
-    verify_divergence_matrix(&ts, TSK_STAT_SITE);
-    verify_divergence_matrix(&ts, TSK_STAT_SITE | TSK_STAT_SPAN_NORMALISE);
+    verify_divergence_matrix(&ts, TSK_STAT_MUTATION);
+    verify_divergence_matrix(&ts, TSK_STAT_MUTATION | TSK_STAT_SPAN_NORMALISE);
 
     tsk_treeseq_free(&ts);
 }
@@ -3856,7 +3870,7 @@ test_simplest_divergence_matrix(void)
     assert_arrays_almost_equal(4, D_site, result);
 
     ret = tsk_treeseq_divergence_matrix(
-        &ts, 2, NULL, sample_ids, 0, NULL, TSK_STAT_SITE, result);
+        &ts, 2, NULL, sample_ids, 0, NULL, TSK_STAT_MUTATION, result);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
     assert_arrays_almost_equal(4, D_site, result);
 
@@ -3866,7 +3880,7 @@ test_simplest_divergence_matrix(void)
     assert_arrays_almost_equal(4, D_branch, result);
 
     ret = tsk_treeseq_divergence_matrix(
-        &ts, 0, NULL, NULL, 0, NULL, TSK_STAT_SITE, result);
+        &ts, 0, NULL, NULL, 0, NULL, TSK_STAT_MUTATION, result);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
     assert_arrays_almost_equal(4, D_site, result);
 
@@ -3879,7 +3893,7 @@ test_simplest_divergence_matrix(void)
     CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_STAT_POLARISED_UNSUPPORTED);
 
     ret = tsk_treeseq_divergence_matrix(
-        &ts, 0, NULL, NULL, 0, NULL, TSK_STAT_SITE | TSK_STAT_BRANCH, result);
+        &ts, 0, NULL, NULL, 0, NULL, TSK_STAT_MUTATION | TSK_STAT_BRANCH, result);
     CU_ASSERT_EQUAL_FATAL(ret, TSK_ERR_MULTIPLE_STAT_MODES);
 
     sample_ids[0] = -1;
@@ -3935,7 +3949,7 @@ test_simplest_divergence_matrix_windows(void)
 
     /* Windows for the second half */
     ret = tsk_treeseq_divergence_matrix(
-        &ts, 2, NULL, sample_ids, 1, windows + 1, TSK_STAT_SITE, result);
+        &ts, 2, NULL, sample_ids, 1, windows + 1, TSK_STAT_MUTATION, result);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
     assert_arrays_almost_equal(4, D_site, result);
     ret = tsk_treeseq_divergence_matrix(
@@ -3985,7 +3999,7 @@ test_simplest_divergence_matrix_internal_sample(void)
     assert_arrays_almost_equal(9, D_branch, result);
 
     ret = tsk_treeseq_divergence_matrix(
-        &ts, 3, NULL, sample_ids, 0, NULL, TSK_STAT_SITE, result);
+        &ts, 3, NULL, sample_ids, 0, NULL, TSK_STAT_MUTATION, result);
     CU_ASSERT_EQUAL_FATAL(ret, 0);
     assert_arrays_almost_equal(9, D_site, result);
 
@@ -4002,8 +4016,8 @@ test_multiroot_divergence_matrix(void)
 
     verify_divergence_matrix(&ts, TSK_STAT_BRANCH);
     verify_divergence_matrix(&ts, TSK_STAT_BRANCH | TSK_STAT_SPAN_NORMALISE);
-    verify_divergence_matrix(&ts, TSK_STAT_SITE);
-    verify_divergence_matrix(&ts, TSK_STAT_SITE | TSK_STAT_SPAN_NORMALISE);
+    verify_divergence_matrix(&ts, TSK_STAT_MUTATION);
+    verify_divergence_matrix(&ts, TSK_STAT_MUTATION | TSK_STAT_SPAN_NORMALISE);
 
     tsk_treeseq_free(&ts);
 }

c/tskit/trees.c

@@ -2041,13 +2041,18 @@ tsk_treeseq_general_stat(const tsk_treeseq_t *self, tsk_size_t state_dim,
     bool stat_site = !!(options & TSK_STAT_SITE);
     bool stat_branch = !!(options & TSK_STAT_BRANCH);
     bool stat_node = !!(options & TSK_STAT_NODE);
+    bool stat_mutation = !!(options & TSK_STAT_MUTATION);
     double default_windows[] = { 0, self->tables->sequence_length };
     tsk_size_t row_size;
 
     /* If no mode is specified, we default to site mode */
-    if (!(stat_site || stat_branch || stat_node)) {
+    if (!(stat_site || stat_branch || stat_node || stat_mutation)) {
         stat_site = true;
     }
+    if (stat_mutation) {
+        ret = tsk_trace_error(TSK_ERR_UNSUPPORTED_STAT_MODE);
+        goto out;
+    }
     /* It's an error to specify more than one mode */
     if (stat_site + stat_branch + stat_node > 1) {
         ret = tsk_trace_error(TSK_ERR_MULTIPLE_STAT_MODES);
@@ -8751,11 +8756,11 @@ remap_to_sample_sets(const tsk_size_t num_samples, const tsk_id_t *restrict samp
 }
 
 static int
-tsk_treeseq_divergence_matrix_site(const tsk_treeseq_t *self, tsk_size_t num_sample_sets,
-    const tsk_id_t *restrict sample_set_index_map, const tsk_size_t num_samples,
-    const tsk_id_t *restrict samples, tsk_size_t num_windows,
-    const double *restrict windows, tsk_flags_t TSK_UNUSED(options),
-    double *restrict result)
+tsk_treeseq_divergence_matrix_mutation(const tsk_treeseq_t *self,
+    tsk_size_t num_sample_sets, const tsk_id_t *restrict sample_set_index_map,
+    const tsk_size_t num_samples, const tsk_id_t *restrict samples,
+    tsk_size_t num_windows, const double *restrict windows,
+    tsk_flags_t TSK_UNUSED(options), double *restrict result)
 {
     int ret = 0;
     tsk_size_t i;
@@ -8913,20 +8918,21 @@ tsk_treeseq_divergence_matrix(const tsk_treeseq_t *self, tsk_size_t num_sample_s
     bool stat_site = !!(options & TSK_STAT_SITE);
     bool stat_branch = !!(options & TSK_STAT_BRANCH);
     bool stat_node = !!(options & TSK_STAT_NODE);
+    bool stat_mutation = !!(options & TSK_STAT_MUTATION);
     tsk_id_t *sample_set_index_map
         = tsk_malloc(num_nodes * sizeof(*sample_set_index_map));
     tsk_size_t j;
 
-    if (stat_node) {
+    if (stat_node || stat_site) {
         ret = tsk_trace_error(TSK_ERR_UNSUPPORTED_STAT_MODE);
         goto out;
     }
-    /* If no mode is specified, we default to site mode */
-    if (!(stat_site || stat_branch)) {
-        stat_site = true;
+    /* If no mode is specified, we default to mutation mode */
+    if (!(stat_mutation || stat_branch)) {
+        stat_mutation = true;
     }
     /* It's an error to specify more than one mode */
-    if (stat_site + stat_branch > 1) {
+    if (stat_mutation + stat_branch > 1) {
         ret = tsk_trace_error(TSK_ERR_MULTIPLE_STAT_MODES);
         goto out;
     }
@@ -8982,8 +8988,8 @@ tsk_treeseq_divergence_matrix(const tsk_treeseq_t *self, tsk_size_t num_sample_s
         ret = tsk_treeseq_divergence_matrix_branch(self, N, sample_set_sizes,
             sample_sets, num_windows, windows, options, result);
     } else {
-        tsk_bug_assert(stat_site);
-        ret = tsk_treeseq_divergence_matrix_site(self, N, sample_set_index_map,
+        tsk_bug_assert(stat_mutation);
+        ret = tsk_treeseq_divergence_matrix_mutation(self, N, sample_set_index_map,
             total_samples, sample_sets, num_windows, windows, options, result);
     }
     if (ret != 0) {

c/tskit/trees.h

@@ -47,6 +47,7 @@ extern "C" {
 #define TSK_STAT_SITE               (1 << 0)
 #define TSK_STAT_BRANCH             (1 << 1)
 #define TSK_STAT_NODE               (1 << 2)
+#define TSK_STAT_MUTATION           (1 << 3)
 
 /* Leave room for other stat types */
 #define TSK_STAT_POLARISED               (1 << 10)

docs/python-api.md

@@ -323,6 +323,7 @@ Single site
       TreeSequence.genetic_relatedness_weighted
       TreeSequence.genetic_relatedness_vector
       TreeSequence.genetic_relatedness_matrix
+      TreeSequence.divergence_matrix
       TreeSequence.general_stat
       TreeSequence.segregating_sites
       TreeSequence.sample_count_stat

docs/stats.md

@@ -204,29 +204,43 @@ statistics API in use.
 
 #### Mode
 
-There are three **modes** of statistic: `site`, `branch`, and `node`,
+There are four **modes** of statistic: `site`, `branch`, `mutation`, and `node`,
 that each summarize aspects of the tree sequence in different but related ways.
 Roughly speaking, these answer the following sorts of question:
 
 site
-: How many mutations differentiate these two genomes?
+: At how many sites do these two genomes differ?
 
 branch
 : How long since these genomes' common ancestor?
 
+branch
+: How many mutations have occurred since these genomes' common ancestor?
+
 node
 : On how much of the genome is each node an ancestor of only one of these genomes, but not both?
 
-These three examples can all be answered in the same way with the tree sequence:
+These examples can all be answered in the same way with the tree sequence:
 first, draw all the paths from one genome to the other through the tree sequence
 (back up to their common ancestor and back down in each marginal tree).
 Then,
-(`site`) count the number of mutations falling on the paths,
+(`site`) check if the allelic state at the endpoints is the same, or
 (`branch`) measure the length of the paths, or
+(`mutation`) count the number of mutations falling on the paths, or
 (`node`) count how often the path goes through each node.
 There is more discussion of this correspondence in the paper describing these statistics,
 and precise definitions are given in each statistic.
 
+Furthermore, since mutations are rare, if two genomes differ in state at a site,
+then it's probably because there was a single mutation.
+This means that `site` and `mutation` are *almost* the same:
+and, they will be identical if there is at most one mutation per site
+and no mutations are silent.
+So, we tend to think of `site` in the same way, as counting numbers of mutations;
+with a small correction due to multiple mutations.
+Furthermore, `mutation` mode is not implemented for many statistics at this point,
+so use `site` to answer the same questions.
+
 Here's an example of using the {meth}`~TreeSequence.diversity` statistic to return the
 average branch length between all pairs of samples:
 
@@ -235,7 +249,7 @@ ts.diversity(mode="branch")
 ```
 
 One important thing to know is that `node` statistics have somewhat different output.
-While `site` and `branch` statistics naturally return one number
+While `site`, `branch`, and `mutation` statistics naturally return one number
 for each portion of the genome (and thus incorporates information about many nodes: see below),
 the `node` statistics return one number **for each node** in the tree sequence (and for each window).
 There can be a lot of nodes in the tree sequence, so beware.
@@ -795,8 +809,9 @@ pairwise) and Patterson's f statistics (which are four-way).
 ### Derived statistics
 
 Most statistics have the property that `mode="branch"` and
-`mode="site"` are "dual" in the sense that they are equal, on average, under
-a high neutral mutation rate. {meth}`~TreeSequence.Fst` and {meth}`~TreeSequence.Tajimas_D`
+`mode="mutation"` are "dual" in the sense that they are equal, on average, 
+and so is `mode="site"` under the infinite sites model of mutation.
+{meth}`~TreeSequence.Fst` and {meth}`~TreeSequence.Tajimas_D`
 do not have this property (since both are ratios of statistics that do have this property).
 
 

python/CHANGELOG.rst

@@ -14,6 +14,10 @@
   is not the case.
   (:user:`benjeffery`, :issue:`2729`, :issue:`2732`, :pr:`3212`).
 
+- The previous behavior of ``genetic_relatedness_matrix`` with ``mode="site"`` is
+  now provided by ``mode="branch"``, and ``mode="site"`` is not supported.
+  (:user:`petrelharp`, :pr:`3314`)
+
 - Drop Python 3.9 support, require Python >= 3.10 (:pr:`3267`, :user:`benjeffery`)
 
 
@@ -59,6 +63,12 @@
   allowing greater flexibility in "disjoint union" situations.
   (:user:`hyanwong`, :user:`petrelharp`, :issue:`3181`)
 
+- Add a new statistics mode, ``mode="mutation"``, that counts numbers of mutations
+  rather than branch length (``mode="branch"``) or allelic state (``mode="site"``)
+  (:user:`petrelharp`, :pr:`3314`).
+
+- Add ``TreeSequence.divergence_matrix``, which was previously undocumented.
+
 **Bugfixes**
 
 - In some tables with mutations out-of-order ``TableCollection.sort`` did not re-order
@@ -102,6 +112,7 @@
 - ``ltrim``, ``rtrim``, ``trim`` and ``shift`` raise an error if used on a tree sequence
   containing a reference sequence (:user:`hyanwong`, :pr:`3210`, :issue:`2091`)
 
+
 --------------------
 [0.6.4] - 2025-05-21
 --------------------

python/_tskitmodule.c

@@ -6442,6 +6442,8 @@ parse_stats_mode(char *mode, tsk_flags_t *ret)
         value = TSK_STAT_BRANCH;
     } else if (strcmp(mode, "node") == 0) {
         value = TSK_STAT_NODE;
+    } else if (strcmp(mode, "mutation") == 0) {
+        value = TSK_STAT_MUTATION;
     } else {
         PyErr_SetString(PyExc_ValueError, "Unrecognised stats mode");
         return -1;