Documentation for IBD segments

Closes #2013
Closes #1716
Closes #1682
Closes #1657

Author: jeromekelleher

docs/_toc.yml

@@ -13,6 +13,7 @@ parts:
   chapters:
   - file: stats
   - file: topological-analysis
+  - file: ibd
 - caption: Interfaces
   chapters:
   - file: python-api

docs/ibd.md

@@ -0,0 +1,201 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.12
+    jupytext_version: 1.9.1
+kernelspec:
+  display_name: Python 3
+  language: python
+  name: python3
+---
+
+```{currentmodule} tskit
+```
+
+
+(sec_identity)=
+
+# Identity by descent
+
+The {meth}`.TreeSequence.ibd_segments` allows us to compute
+segments of identity by descent.
+
+:::{note}
+This documentation page is preliminary
+:::
+
+## Examples
+
+Let's take a simple tree sequence to illustrate the {meth}`.TreeSequence.ibd_segments`
+method and associated {ref}`sec_python_api_reference_identity`:
+
+```{code-cell}
+:tags: [hide-input]
+
+import tskit
+import io
+from IPython.display import SVG
+
+nodes = io.StringIO(
+    """\
+    id      is_sample   time
+    0       1           0
+    1       1           0
+    2       1           0
+    3       0           1
+    4       0           2
+    5       0           3
+    """
+)
+edges = io.StringIO(
+    """\
+    left    right   parent  child
+    2     10     3       0
+    2     10     3       2
+    0     10     4       1
+    0     2      4       2
+    2     10     4       3
+    0     2      5       0
+    0     2      5       4
+    """
+)
+ts = tskit.load_text(nodes=nodes, edges=edges, strict=False)
+
+SVG(ts.draw_svg())
+```
+
+### Definition
+
+A pair of nodes ``(u, v)`` has an IBD segment with a left and right
+coordinate ``[left, right)`` and ancestral node ``a`` iff the most
+recent common ancestor of the segment ``[left, right)`` in nodes ``u``
+and ``v`` is ``a``, and the segment has been inherited along the same
+genealogical path (ie. it has not been broken by recombination). The
+segments returned are the longest possible ones.
+
+Consider the IBD segments that we get from our example tree sequence:
+
+```{code-cell}
+segs = ts.ibd_segments(store_segments=True)
+for pair, segment_list in segs.items():
+    print(pair, list(segment_list))
+```
+
+Each of the sample pairs (0, 1), (0, 2) and (1, 2) is associated with
+two IBD segments, representing the different paths from these sample
+pairs to their common ancestor. Note in particular that (1, 2) has
+**two** IBD segments rather than one: even though the MRCA is
+4 in both cases, the paths from the samples to the MRCA are different
+in the left and right trees.
+
+
+### Data structures
+
+The result of calling {meth}`.TreeSequence.ibd_segments` is an
+{class}`.IdentitySegments` class:
+
+```{code-cell}
+segs = ts.ibd_segments()
+print(segs)
+```
+
+By default this class only stores the high-level summaries of the
+IBD segments discovered. As we can see in this example, we have a
+total of six segments and
+the total span (i.e., the sum lengths of the genomic intervals spanned
+by IBD segments) is 30.
+
+If required, we can get more detailed information about particular
+segment pairs and the actual segments using the ``store_pairs``
+and ``store_segments`` arguments.
+
+:::{warning}
+Only use the ``store_pairs`` and ``store_segments`` arguments if you
+really need this information! The number of IBD segments can be
+very large and storing them all requires a lot of memory. It is
+also much faster to just compute the overall summaries, without
+needing to store the actual lists.
+:::
+
+
+```{code-cell}
+segs = ts.ibd_segments(store_pairs=True)
+for pair, value in segs.items():
+    print(pair, "::", value)
+```
+
+Now we can see the more detailed breakdown of how the identity segments
+are distributed among the sample pairs. The {class}`.IdentitySegments`
+class behaves like a dictionary, such that ``segs[(a, b)]`` will return
+the {class}`.IdentitySegmentList` instance for that pair of samples:
+
+```{code-cell}
+seglist = segs[(0, 1)]
+print(seglist)
+```
+
+If we want to access the detailed information about the actual
+identity segments, we must use the ``store_segments`` argument:
+
+```{code-cell}
+segs = ts.ibd_segments(store_pairs=True, store_segments=True)
+segs[(0, 1)]
+```
+
+The {class}`.IdentitySegmentList` behaves like a Python list,
+where each element is an instance of {class}`.IdentitySegment`.
+
+:::{warning}
+The order of segments in an {class}`.IdentitySegmentList`
+is arbitrary, and may change in future versions.
+:::
+
+
+```{eval-rst}
+.. todo:: More examples using the other bits of the IdentitySegments
+    API here
+```
+
+### Controlling the sample sets
+
+By default we get the IBD segments between all pairs of
+:ref:`sample<sec_data_model_definitions_samples>` nodes.
+
+#### IBD within a sample set
+We can reduce this to pairs within a specific set using the
+``within`` argument::
+
+
+```{eval-rst}
+.. todo:: More detail and better examples here.
+```
+
+```{code-cell}
+segs = ts.ibd_segments(within=[0, 2], store_pairs=True)
+print(list(segs.keys()))
+```
+
+#### IBD between sample sets
+
+We can also compute IBD **between** sample sets:
+
+```{code-cell}
+segs = ts.ibd_segments(between=[[0,1], [2]], store_pairs=True)
+print(list(segs.keys()))
+```
+
+:::{seealso}
+See the {meth}`.TreeSequence.ibd_segments` documentation for
+more details.
+:::
+
+### Constraints on the segments
+
+The ``max_time`` and ``min_length`` arguments allow us to constrain the
+segments that we consider.
+
+```{eval-rst}
+.. todo:: Add examples for these arguments.
+```

docs/python-api.md

@@ -220,6 +220,20 @@ which perform the same actions but modify the {class}`TableCollection` in place.
   TreeSequence.trim
 ```
 
+(sec_python_api_tree_sequences_ibd)=
+
+#### Identity by descent
+
+The {meth}`.TreeSequence.ibd_segments` method allows us to compute
+identity relationships between pairs of samples. See the
+{ref}`sec_identity` section for more details and examples
+and the {ref}`sec_python_api_reference_identity` section for
+API documentation on the associated classes.
+
+```{eval-rst}
+.. autosummary::
+  TreeSequence.ibd_segments
+```
 
 (sec_python_api_tree_sequences_tables)=
 
@@ -471,7 +485,7 @@ high performance interface which can be used in conjunction with the equivalent
 
 Moving around within a tree usually involves visiting the tree nodes in some sort of
 order. Often, given a particular order, it is convenient to iterate over each node
-using the :meth:`Tree.nodes` method. However, for high performance algorithms, it
+using the {meth}`Tree.nodes` method. However, for high performance algorithms, it
 may be more convenient to access the node indices for a particular order as
 an array, and use this, for example, to index into one of the node arrays (see
 {ref}`sec_topological_analysis_traversal`).
@@ -646,7 +660,7 @@ Other properties
 These methods act in-place to transform the contents of a {class}`TableCollection`,
 either by modifying the underlying tables (removing, editing, or adding to them) or
 by adjusting the table collection so that it meets the
-{ref}`sec_valid_tree_sequence_requirements.
+{ref}`sec_valid_tree_sequence_requirements`.
 
 
 (sec_tables_api_modification)=
@@ -707,14 +721,12 @@ Indexing
       TableCollection.drop_index
 ```
 
-
 #### Miscellaneous methods
 
 ```{eval-rst}
 .. autosummary::
   TableCollection.copy
   TableCollection.equals
-  TableCollection.ibd_segments
   TableCollection.link_ancestors
 ```
 
@@ -1413,6 +1425,34 @@ basic class, where each attribute matches an identically named attribute in the
     :inherited-members:
 ```
 
+(sec_python_api_reference_identity)=
+
+### Identity classes
+
+The classes documented in this section are associated with summarising
+identity relationships between pairs of samples. See the {ref}`sec_identity`
+section for more details and examples.
+
+#### The {class}`IdentitySegments` class
+
+```{eval-rst}
+.. autoclass:: IdentitySegments()
+    :members:
+```
+
+#### The {class}`IdentitySegmentList` class
+
+```{eval-rst}
+.. autoclass:: IdentitySegmentList()
+    :members:
+```
+
+#### The {class}`IdentitySegment` class
+
+```{eval-rst}
+.. autoclass:: IdentitySegment()
+    :members:
+```
 
 ### Miscellaneous classes
 

python/CHANGELOG.rst

@@ -55,6 +55,10 @@
 
 **Features**
 
+- Add the ``ibd_segments`` method and associated classes to compute, summarise
+  and store segments of identity by descent from a tree sequence
+  (:user:`gtsambos`, :user:`jeromekelleher`).
+
 - Allow skipping of site and mutation tables in ``TableCollection.sort``
   (:user:`benjeffery`, :issue:`1475`, :pr:`1826`).
 
@@ -114,6 +118,7 @@
 
 - tskit now supports python 3.10 (:user:`benjeffery`, :issue:`1895`, :pr:`1949`)
 
+
 **Fixes**
 
 - `dump_tables` omitted individual parents. (:user:`benjeffery`, :issue:`1828`, :pr:`1884`)