# MIT License
# Copyright (c) 2018-2021 Tskit Developers
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
Module responsible for computing various statistics on tree sequences.
import numpy as np
Class for calculating `linkage disequilibrium
between pairs of sites in a :class:`TreeSequence`.
.. note:: This interface is deprecated and a replacement is planned.
Please see https://github.com/tskit-dev/tskit/issues/1900 for
more information. Note also that the current implementation is
quite limited (see warning below).
.. warning:: This class does not currently support sites that have more than one
mutation. Using it on such a tree sequence will raise a LibraryError with
an "Only infinite sites mutations supported" message.
Silent mutations are also not supported and will result in a LibraryError.
:param TreeSequence tree_sequence: The tree sequence of interest.
def __init__(self, tree_sequence):
self._tree_sequence = tree_sequence
self._ll_ld_calculator = _tskit.LdCalculator(
# To protect low-level C code, only one method may execute on the
# low-level objects at one time.
self._instance_lock = threading.Lock()
def get_r2(self, a, b):
# Deprecated alias for r2(a, b)
return self.r2(a, b)
[docs] def r2(self, a, b):
Returns the value of the :math:`r^2` statistic between the pair of
sites at the specified indexes. This method is *not* an efficient
method for computing large numbers of pairwise LD values; please use either
:meth:`.r2_array` or :meth:`.r2_matrix` for this purpose.
:param int a: The index of the first site.
:param int b: The index of the second site.
:return: The value of :math:`r^2` between the sites at indexes
``a`` and ``b``.
return self._ll_ld_calculator.get_r2(a, b)
def get_r2_array(self, a, direction=1, max_mutations=None, max_distance=None):
# Deprecated alias for r2_array
[docs] def r2_array(
self, a, direction=1, max_mutations=None, max_distance=None, max_sites=None
Returns the value of the :math:`r^2` statistic between the focal
site at index :math:`a` and a set of other sites. The method
operates by starting at the focal site and iterating over adjacent
sites (in either the forward or backwards direction) until either a
maximum number of other sites have been considered (using the
``max_sites`` parameter), a maximum distance in sequence
coordinates has been reached (using the ``max_distance`` parameter) or
the start/end of the sequence has been reached. For every site
:math:`b` considered, we then insert the value of :math:`r^2` between
:math:`a` and :math:`b` at the corresponding index in an array, and
return the entire array. If the returned array is :math:`x` and
``direction`` is :data:`tskit.FORWARD` then :math:`x` is the
value of the statistic for :math:`a` and :math:`a + 1`, :math:`x`
the value for :math:`a` and :math:`a + 2`, etc. Similarly, if
``direction`` is :data:`tskit.REVERSE` then :math:`x` is the
value of the statistic for :math:`a` and :math:`a - 1`, :math:`x`
the value for :math:`a` and :math:`a - 2`, etc.
:param int a: The index of the focal sites.
:param int direction: The direction in which to travel when
examining other sites. Must be either
:data:`tskit.FORWARD` or :data:`tskit.REVERSE`. Defaults
:param int max_sites: The maximum number of sites to return
:math:`r^2` values for. Defaults to as many sites as
:param int max_mutations: Deprecated synonym for max_sites.
:param float max_distance: The maximum absolute distance between
the focal sites and those for which :math:`r^2` values
:return: An array of double precision floating point values
representing the :math:`r^2` values for sites in the
if max_mutations is not None and max_sites is not None:
raise ValueError("max_mutations is a deprecated synonym for max_sites")
if max_mutations is not None:
max_sites = max_mutations
max_sites = -1 if max_sites is None else max_sites
if max_distance is None:
max_distance = sys.float_info.max
# Deprecated alias for r2_matrix
[docs] def r2_matrix(self):
Returns the complete :math:`m \\times m` matrix of pairwise
:math:`r^2` values in a tree sequence with :math:`m` sites.
:return: An 2 dimensional square array of double precision
floating point values representing the :math:`r^2` values for
all pairs of sites.
m = self._tree_sequence.num_sites
A = np.ones((m, m), dtype=float)
for j in range(m - 1):
a = self.get_r2_array(j)
A[j, j + 1 :] = a
A[j + 1 :, j] = a