Spectrum classes

Spectrum classes#

The spectrum classes and their named-collection types.

`Spectrum`	Class for holding and manipulating a site-frequency spectrum.
`Spectra`	Class for holding and manipulating site-frequency spectra of multiple types.
`JointSFS`	A joint (multi-population) site-frequency spectrum.
`JointSpectra`	A collection of joint (multi-population) site-frequency spectra keyed by type, the multi-population analogue of `Spectra`.
`TwoSFS`	A 2-dimensional site-frequency spectrum, i.e. a square matrix whose entry `(i, j)` relates the number of derived alleles at a pair of frequency classes `i` and `j` of a single population (for example the second moment of the SFS).
`TwoLocusSFS`	The two-locus site-frequency spectrum under recombination: a square matrix whose entry `(i, j)` is the expected product of the branch length subtending `i` samples at locus 0 and `j` samples at locus 1, for two loci separated by a given recombination rate `r`.
`TwoSpectra`	A collection of two-site (two-locus) site-frequency spectra keyed by type, the two-dimensional analogue of `Spectra`.
`AbstractSpectrum`	Abstract base class for site-frequency spectrum containers.
`AbstractSpectra`	Abstract base class for a collection of site-frequency spectra keyed by type, for example the SFS stratified into neutral and selected sites.

Spectrum#

class Spectrum(data: Sequence[float])[source]#

Bases: AbstractSpectrum

Class for holding and manipulating a site-frequency spectrum.

Methods:

`__init__`	Initialize spectrum.
`to_list`	Convert to list.
`to_spectra`	Convert to Spectra object.
`to_file`	Save object to file.
`from_file`	Load object from file.
`to_numpy`	Convert to array.
`fold`	Fold the site-frequency spectrum.
`misidentify`	Introduce ancestral misidentification at rate epsilon.
`subsample`	Subsample spectrum to a given sample size.
`resample`	Resample SFS assuming independent Poisson counts.
`is_folded`	Check if the site-frequency spectrum is folded.
`normalize`	Normalize SFS so that all non-monomorphic counts add up to 1.
`copy`	Copy the spectrum.
`from_polymorphic`	Create Spectrum from polymorphic counts only.
`from_list`	Create Spectrum from list.
`from_polydfe`	Create Spectrum from polyDFE specification which treats the number of mutational target sites and the divergence counts separately.
`plot`	Plot spectrum.
`kingman`	The standard (Kingman) neutral site-frequency spectrum for a sample of size `n`.
`standard_kingman`	Alias of `kingman()`.
`get_neutral`	Obtain a standard neutral SFS for a given theta and number of sites.
`scale_theta`	Scale the spectrum to a different theta value by
`from_json`	Load from a JSON string.
`to_json`	Convert to a JSON string.

__init__(data: Sequence[float])[source]#

Initialize spectrum.

Parameters:: data (Sequence[float]) – SFS counts

data: ndarray#: The underlying array holding the spectrum.

property n: int#

The sample size.

Returns:: Sample size

property n_sites: float#

The total number of sites.

Returns:: Total number of sites

property n_div: float#

Number of divergence counts.

Returns:: Number of divergence counts

property has_div: bool#

Whether n_div was specified.

Returns:: Whether n_div was specified

property n_monomorphic: float#

Number of monomorphic sites.

Returns:: Number of monomorphic sites

property polymorphic: ndarray#

Get the polymorphic counts.

Returns:: Polymorphic counts

property n_polymorphic: float#

The total number of polymorphic sites.

Returns:: Total number of polymorphic sites

to_list()[source]#

Convert to list.

Return type:: list
Returns:: SFS counts

to_spectra()[source]#

Convert to Spectra object.

Return type:: Spectra
Returns:: Spectra object

to_file(file: str)[source]#

Save object to file.

Parameters:: file (str) – File name

static from_file(file: str)[source]#

Load object from file.

Parameters:: file (str) – File name
Return type:: Spectrum
Returns:: Spectrum object

to_numpy()[source]#

Convert to array.

Return type:: ndarray
Returns:: SFS counts

property theta: float#: Calculate site-wise population mutation rate using Watterson’s estimator. Note that theta is given per site, i.e. Watterson’s estimator is divided by the total number of sites (n_sites).

property Theta: float#: Calculate genome-wide population mutation rate using Watterson’s estimator.

Note

Property Theta is not normalized by the total number of sites, unlike theta.

fold()[source]#

Fold the site-frequency spectrum.

Return type:: Spectrum
Returns:: Folded spectrum

misidentify(epsilon: float)[source]#

Introduce ancestral misidentification at rate epsilon. Note that monomorphic counts won’t be affected.

Parameters:: epsilon (float) – Misidentification rate (0 <= epsilon <= 1)
Return type:: Spectrum
Returns:: Spectrum with misidentification applied
Raises:: ValueError – If epsilon is not between 0 and 1

subsample(n: int, mode: Literal['random', 'probabilistic'] = 'probabilistic', seed: int | Generator = None)[source]#

Subsample spectrum to a given sample size.

Warning

If using the ‘random’ mode, The SFS counts are cast to integers before subsampling so this will only provide sensible results if the SFS counts are integers or if they are large enough to be approximated well by integers. The ‘probabilistic’ mode does not have this limitation.

Parameters:

n (int) – Sample size
mode (Literal['random', 'probabilistic']) – Subsampling mode. Either ‘random’ or ‘probabilistic’.
seed (int | Generator) – Random state or seed. Only for ‘random’ mode.

Return type:

Spectrum

Returns:

Subsampled spectrum

resample(seed: int | Generator = None)[source]#

Resample SFS assuming independent Poisson counts.

Parameters:: seed (int | Generator) – Random state or seed
Return type:: Spectrum
Returns:: Resampled spectrum.

is_folded()[source]#

Check if the site-frequency spectrum is folded.

Return type:: bool
Returns:: True if folded, False otherwise

normalize()[source]#

Normalize SFS so that all non-monomorphic counts add up to 1.

Return type:: Spectrum
Returns:: Normalized spectrum

copy()[source]#

Copy the spectrum.

Return type:: Spectrum
Returns:: Copy of the spectrum

static from_polymorphic(data: Sequence)[source]#

Create Spectrum from polymorphic counts only.

Parameters:: data (Sequence) – Polymorphic counts
Return type:: Spectrum
Returns:: Spectrum

static from_list(data: Sequence)[source]#

Create Spectrum from list.

Parameters:: data (Sequence) – SFS counts
Return type:: Spectrum
Returns:: Spectrum

static from_polydfe(polymorphic: Sequence, n_sites: float, n_div: float)[source]#

Create Spectrum from polyDFE specification which treats the number of mutational target sites and the divergence counts separately.

Parameters:

polymorphic (Sequence) – Polymorphic counts
n_sites (float) – Total number of sites
n_div (float) – Number of divergence counts

Return type:

Spectrum

Returns:

Spectrum

plot(show: bool = True, file: str = None, title: str = None, log_scale: bool = False, show_monomorphic: bool = False, kwargs_legend: dict = {'prop': {'size': 8}}, ax: plt.Axes = None)[source]#

Plot spectrum.

Parameters:

show (bool) – Whether to show plot.
file (str) – File to save plot to.
title (str) – Title of plot.
log_scale (bool) – Whether to use log scale on y-axis.
show_monomorphic (bool) – Whether to show monomorphic counts.
kwargs_legend (dict) – Keyword arguments passed to plt.legend(). Only for Python visualization backend.
ax (plt.Axes) – Axes to plot on. Only for Python visualization backend.

Return type:

plt.Axes

Returns:

Axes

static kingman(n: int, n_monomorphic: int = 0)[source]#

The standard (Kingman) neutral site-frequency spectrum for a sample of size n.

Parameters:

n (int) – sample size
n_monomorphic (int) – number of monomorphic sites placed in the zero-frequency bin

Return type:

Spectrum

Returns:

the Kingman SFS

static standard_kingman(n: int, n_monomorphic: int = 0)[source]#

Alias of kingman().

Deprecated since version 1.0.0: Use kingman() instead.

Parameters:

n (int) – sample size
n_monomorphic (int) – Number of monomorphic sites.

Return type:

Spectrum

Returns:

Standard Kingman SFS

static get_neutral(theta: float, n_sites: float, n: int, r: Sequence[float] = None)[source]#

Obtain a standard neutral SFS for a given theta and number of sites.

Parameters:

theta (float) – Population mutation rate
n_sites (float) – Number of total sites
n (int) – Number of frequency classes
r (Sequence[float]) – Nuisance parameters that account for demography. An array of length n-1 whose elements are multiplied element-wise with the polymorphic counts of the Kingman SFS. By default, no demography effects are considered which is equivalent to r = [1] * (n-1). Note that non-default values of r will also affect estimates of the population mutation rate.

Return type:

Spectrum

Returns:

Neutral SFS

scale_theta(theta: float)[source]#

Scale the spectrum to a different theta value by

Parameters:: theta (float) – New theta value
Return type:: Spectrum
Returns:: Scaled spectrum

classmethod from_json(json: str)#

Load from a JSON string.

Parameters:: json (str) – JSON string.
Return type:: Self
Returns:: Spectrum
Raises:: ValueError – If the payload does not decode to this class.

property shape: Tuple[int, ...]#: The shape of the underlying array.

to_json()#

Convert to a JSON string.

Return type:: str
Returns:: JSON string

Spectra#

class Spectra(data: Dict[str, Iterable])[source]#

Bases: AbstractSpectra

Class for holding and manipulating site-frequency spectra of multiple types.

Methods:

`__init__`	Initialize spectra.
`from_list`	Create from array of spectra.
`normalize`	Normalize spectra by sum of all entries.
`to_file`	Save object to file.
`to_spectra`	Convert to dictionary of spectrum objects.
`to_dataframe`	Get representation as dataframe.
`to_numpy`	Convert to numpy array.
`to_list`	Convert to nested list.
`to_dict`	Convert to dictionary.
`select`	Select types.
`copy`	Copy object.
`get_empty`	Get a Spectra object with zero counts but having the same shape and types as self.
`merge_groups`	Group over given levels and sum up spectra so the spectra are summed over the levels that were not specified.
`has_dots`	Check whether column names contain dots.
`replace_dots`	Replace dots in column names with a given string.
`combine`	Merge types of two Spectra objects.
`from_dict`	Load from nested dictionary first indexed by types and then by samples.
`from_dataframe`	Load Spectra object from dataframe.
`from_file`	Save object to file.
`from_spectra`	Create from dict of spectrum objects indexed by type.
`from_spectrum`	Create from single spectrum object.
`to_spectrum`	Convert to Spectrum object by summing over all types.
`plot`	Visualize spectra.
`drop_empty`	Remove types whose spectra have no counts.
`drop_zero_entries`	Remove types whose spectra have some zero entries.
`drop_sparse`	Remove types whose spectra have fewer than equal `n_polymorphic` polymorphic sites.
`rename`	Rename types.
`prefix`	Prefix types, i.e. 'type' -> 'prefix.type' for all types.
`reorder_levels`	Reorder levels.
`print`	Print spectra.
`fold`	Fold spectra.
`subsample`	Subsample spectra to a given sample size.
`resample`	Resample SFS assuming independent Poisson counts.
`is_folded`	Check whether spectra are folded.
`sort_types`	Sort types alphabetically.

__init__(data: Dict[str, Iterable])[source]#

Initialize spectra.

Parameters:: data (Dict[str, Iterable]) – Dictionary of SFS counts keyed by type

property n: int#

The sample size.

Returns:: Sample size

property k: int#

The number of types.

Returns:: Number of types

property n_monomorphic: Series#

The number of monomorphic sites.

Returns:: Number of monomorphic sites

property polymorphic: ndarray#

The polymorphic counts.

Returns:: Polymorphic counts

property n_polymorphic: Series#

The total number of polymorphic sites per type.

Returns:: Total number of polymorphic sites for each type

static from_list(data: Sequence, types: List)[source]#

Create from array of spectra. Note that data.ndim needs to be 2.

Parameters:

data (Sequence) – Array of spectra
types (List) – Types

Return type:

Spectra

Returns:

Spectra

property types: List[str]#

The types.

Returns:: Types

property n_sites: Series#

The number of mutational target sites which is the sum of all SFS entries.

Returns:: Number of mutational target sites for each type

property n_div: Series#

The number of divergence counts.

Returns:: Number of divergence counts for each type

property has_div: Series#

Whether n_div was specified.

Returns:: Whether n_div was specified for each type

property theta: Series#: Calculate site-wise population mutation rate using Watterson’s estimator. Note that theta is given per site, i.e. Watterson’s estimator is divided by the total number of sites (n_sites).

property Theta: Series#: Calculate genome-wide population mutation rate using Watterson’s estimator.

Note

Property Theta is not normalized by the total number of sites, unlike theta.

normalize()[source]#

Normalize spectra by sum of all entries.

Return type:: Spectra
Returns:: Normalized spectra

to_file(file: str)[source]#

Save object to file.

Parameters:: file (str) – File name

to_spectra()[source]#

Convert to dictionary of spectrum objects.

Return type:: Dict[str, Spectrum]
Returns:: Dictionary of spectrum objects

to_dataframe()[source]#

Get representation as dataframe.

Return type:: DataFrame
Returns:: Dataframe

to_numpy()[source]#

Convert to numpy array.

Return type:: ndarray
Returns:: Numpy array

to_list()[source]#

Convert to nested list.

Return type:: list
Returns:: Nested list

to_dict()[source]#

Convert to dictionary.

Return type:: dict
Returns:: Dictionary of lists

select(keys: str | List[str] | ndarray | tuple, use_regex: bool = True)[source]#

Select types. Alias for __getitem__.

Parameters:

keys (Union[str, List[str], ndarray, tuple]) – String or list of strings, possibly regex to match type names
use_regex (bool) – Whether to use regex to match type names

Return type:

Spectra

Returns:

Spectrum or Spectra depending on the number of matches

copy()[source]#

Copy object.

Return type:: Spectra
Returns:: Copy of object

get_empty()[source]#

Get a Spectra object with zero counts but having the same shape and types as self.

Return type:: Spectra
Returns:: Spectra object with zero counts

merge_groups(level: List[int] | int = 0)[source]#

Group over given levels and sum up spectra so the spectra are summed over the levels that were not specified.

Parameters:: level (Union[List[int], int]) – Level(s) to group over
Return type:: Spectra
Returns:: Spectra object with merged groups
Raises:: ValueError – If some type has no name at one of the requested levels

has_dots()[source]#

Check whether column names contain dots.

Return type:: bool
Returns:: True if column names contain dots, False otherwise

replace_dots(replacement: str = '_')[source]#

Replace dots in column names with a given string.

Parameters:: replacement (str) – Replacement string
Return type:: Spectra
Returns:: Spectra object with replaced dots

property all: Spectrum#

The ‘all’ type equals the sum of all spectra.

Returns:: Spectrum object

combine(s: Spectra)[source]#

Merge types of two Spectra objects.

Parameters:: s (Spectra) – Other Spectra object
Return type:: Spectra
Returns:: Merged Spectra object

static from_dict(data: dict)[source]#

Load from nested dictionary first indexed by types and then by samples.

Parameters:: data (dict) – Dictionary of lists indexed by types
Return type:: Spectra
Returns:: Spectra object

static from_dataframe(data: DataFrame)[source]#

Load Spectra object from dataframe.

Parameters:: data (DataFrame) – Dataframe
Return type:: Spectra
Returns:: Spectra object

classmethod from_file(file: str)[source]#

Save object to file.

Parameters:: file (str) – Path to file, possibly URL
Return type:: Spectra
Returns:: Spectra object

static from_spectra(spectra: Dict[str, Spectrum])[source]#

Create from dict of spectrum objects indexed by type.

Parameters:: spectra (Dict[str, Spectrum]) – Dictionary of spectrum objects indexed by type
Return type:: Spectra
Returns:: Spectra object

static from_spectrum(sfs: Spectrum)[source]#

Create from single spectrum object. The type of the spectrum is set to ‘all’.

Parameters:: sfs (Spectrum) – Spectrum
Return type:: Spectra
Returns:: Spectra object

to_spectrum()[source]#

Convert to Spectrum object by summing over all types.

Return type:: Spectrum
Returns:: Spectrum object

plot(show: bool = True, file: str = None, title: str = None, log_scale: bool = False, use_subplots: bool = False, show_monomorphic: bool = False, kwargs_legend: dict = {'prop': {'size': 8}}, ax: plt.Axes = None)[source]#

Visualize spectra.

Parameters:

show (bool) – Whether to show the plot.
file (str) – File name to save the plot to.
title (str) – Plot title.
log_scale (bool) – Whether to use log scale on y-axis.
use_subplots (bool) – Whether to use subplots. Only for Python visualization backend.
show_monomorphic (bool) – Whether to show monomorphic sites.
kwargs_legend (dict) – Keyword arguments passed to plt.legend(). Only for Python visualization backend.
ax (plt.Axes) – Axes to plot on. Only for Python visualization backend and if use_subplots is False.

Return type:

plt.Axes

Returns:

Axes

drop_empty()[source]#

Remove types whose spectra have no counts.

Return type:: Spectra
Returns:: Spectra with non-empty types

drop_zero_entries()[source]#

Remove types whose spectra have some zero entries. Note that we ignore zero counts in the last entry i.e. fixed derived alleles.

Return type:: Spectra
Returns:: Spectra with non-zero entries

drop_sparse(n_polymorphic: int)[source]#

Remove types whose spectra have fewer than equal n_polymorphic polymorphic sites.

Parameters:: n_polymorphic (int) – Minimum number of polymorphic sites a type must exceed to be kept.
Return type:: Spectra
Returns:: Spectra

rename(names: List[str])[source]#

Rename types.

Parameters:: names (List[str]) – New names
Return type:: Spectra
Returns:: Spectra with renamed types

prefix(prefix: str)[source]#

Prefix types, i.e. ‘type’ -> ‘prefix.type’ for all types.

Parameters:: prefix (str) – Prefix
Return type:: Spectra
Returns:: Spectra with prefixed types

reorder_levels(levels: List[int])[source]#

Reorder levels.

Parameters:: levels (List[int]) – New order of levels
Return type:: Spectra
Returns:: Spectra with reordered levels

print()[source]#: Print spectra.

fold()[source]#

Fold spectra.

Returns:: Folded spectra

subsample(n: int, mode: Literal['random', 'probabilistic'] = 'probabilistic', seed: int | Generator = None)[source]#

Subsample spectra to a given sample size.

Warning

If using the ‘random’ mode, The SFS counts are cast to integers before subsampling so this will only provide sensible results if the SFS counts are integers or if they are large enough to be approximated by integers. The ‘probabilistic’ mode does not have this limitation.

Parameters:

n (int) – Sample size
mode (Literal['random', 'probabilistic']) – Subsampling mode. Either ‘random’ or ‘probabilistic’.
seed (int | Generator) – Random state or seed. Only for ‘random’ mode.

Return type:

Spectra

Returns:

Subsampled spectra

resample(seed: int | Generator = None)[source]#

Resample SFS assuming independent Poisson counts.

Parameters:: seed (int | Generator) – Random state or seed
Return type:: Spectra
Returns:: Resampled spectra.

is_folded()[source]#

Check whether spectra are folded.

Return type:: Dict[str, bool]
Returns:: Dictionary of types and whether they are folded

sort_types()[source]#

Sort types alphabetically.

Return type:: Spectra
Returns:: Sorted spectra object

property is_empty: bool#

Whether the collection holds no parsed mass, either because there are no types or because every spectrum is all zero. Counting types alone does not answer this for the two-site and joint collections, which always carry an all type even when no site made it into the spectra.

Returns:: True if nothing was counted.

JointSFS#

class JointSFS(data: ndarray | list, pop_names: List[str] = None)[source]#

Bases: AbstractSpectrum

A joint (multi-population) site-frequency spectrum.

The data is a P-dimensional array of shape (n_0 + 1, ..., n_{P-1} + 1) where P is the number of populations and entry (k_0, ..., k_{P-1}) counts sites (or branch length) with k_p derived alleles in population p. For two populations this is a 2-dimensional array (analogous to but generally rectangular, unlike the square TwoSFS); for three populations it is a 3-dimensional array, and so on.

Methods:

`__init__`	Construct from a data array.
`marginalize`	Marginalize the joint SFS onto a subset of populations by summing over the other populations.
`fold`	Fold the joint SFS by the total derived count across populations, the dadi/moments convention: ancestral and derived is one site-level property, so every axis is reflected at once.
`is_folded`	Check if the joint SFS is folded.
`plot`	Plot the joint SFS as a 2-dimensional heatmap.
`plot_surface`	Plot the joint SFS of two populations as a surface, with the two allele-count axes on the horizontal plane and the number of sites as height.
`copy`	Create a deep copy.
`from_file`	Load from file.
`from_json`	Load from a JSON string.
`to_file`	Save to file (in JSON format).
`to_json`	Convert to a JSON string.

__init__(data: ndarray | list, pop_names: List[str] = None)[source]#

Construct from a data array.

Parameters:

data (ndarray | list) – A P-dimensional array.
pop_names (List[str]) – Optional names of the P populations (one per axis), used for plot axis labels. Defaults to pop_0, ..., pop_{P-1} when not given.

Raises:

ValueError – If the data is not at least 1-dimensional or if the number of population names does not match the number of axes.

data: ndarray#: The joint SFS array.

pop_names: List[str] = None#

property n_pops: int#: Number of populations (dimensions of the joint SFS).

marginalize(pops: Sequence[int])[source]#

Marginalize the joint SFS onto a subset of populations by summing over the other populations. This is useful for example to obtain a 2-dimensional view of a higher-dimensional joint SFS.

Parameters:: pops (Sequence[int]) – The population indices to keep, in the desired axis order.
Return type:: JointSFS
Returns:: A joint SFS over the specified populations.
Raises:: ValueError – If any population index is out of range.

fold()[source]#

Fold the joint SFS by the total derived count across populations, the dadi/moments convention: ancestral and derived is one site-level property, so every axis is reflected at once. Entries above half the total sample size are emptied into their reflection, and entries on the boundary are shared between the two. This only makes sense for counts or frequencies, and folding a folded joint SFS is a no-op.

Return type:: JointSFS
Returns:: Folded joint SFS.

is_folded()[source]#

Check if the joint SFS is folded.

Return type:: bool
Returns:: Whether the joint SFS is folded.

plot(pops: Tuple[int, int] = (0, 1), ax: plt.Axes = None, title: str = None, log_scale: bool = True, mask_monomorphic: bool = True, cbar_kws: Dict = None, show: bool = True)[source]#

Plot the joint SFS as a 2-dimensional heatmap. For more than two populations, the joint SFS is first marginalized onto the two requested populations (summing over the others). The colour scale is logarithmic by default, since the joint SFS is heavily skewed toward the low-frequency corner.

Parameters:

pops (Tuple[int, int]) – The two population indices to plot (y-axis, x-axis).
ax (plt.Axes) – Axes to plot on.
title (str) – Title of the plot.
log_scale (bool) – Whether to use a logarithmic color scale (default True).
mask_monomorphic (bool) – Whether to mask the monomorphic corners (all-zero and all-derived).
cbar_kws (Dict) – Keyword arguments for the color bar.
show (bool) – Whether to show the plot.

Return type:

plt.Axes

Returns:

Axes.

Raises:

ValueError – If not exactly two populations are requested or the marginalized data is not 2-dimensional.

plot_surface(pops: Tuple[int, int] = (0, 1), ax: plt.Axes = None, title: str = None, log_scale: bool = False, mask_monomorphic: bool = True, cmap: str = 'viridis', show: bool = True)[source]#

Plot the joint SFS of two populations as a surface, with the two allele-count axes on the horizontal plane and the number of sites as height. For more than two populations, the joint SFS is first marginalized onto the two requested populations (summing over the others).

Parameters:

pops (Tuple[int, int]) – The two population indices to plot (y-axis, x-axis).
ax (plt.Axes) – Axes to plot on.
title (str) – Title of the plot.
log_scale (bool) – Whether to use a logarithmic color scale.
mask_monomorphic (bool) – Whether to mask the monomorphic corners (all-zero and all-derived).
cmap (str) – The colormap.
show (bool) – Whether to show the plot.

Return type:

plt.Axes

Returns:

Axes.

Raises:

ValueError – If not exactly two populations are requested or the marginalized data is not 2-dimensional.

copy()#

Create a deep copy.

Return type:: Self
Returns:: Deep copy.

classmethod from_file(file: str)#

Load from file.

Parameters:: file (str) – File path.
Return type:: Self
Returns:: Spectrum

classmethod from_json(json: str)#

Load from a JSON string.

Parameters:: json (str) – JSON string.
Return type:: Self
Returns:: Spectrum
Raises:: ValueError – If the payload does not decode to this class.

property n_sites: float#: The total number of sites, i.e. the sum of all entries.

property shape: Tuple[int, ...]#: The shape of the underlying array.

to_file(file: str)#

Save to file (in JSON format).

Parameters:: file (str) – File path.
Return type:: None

to_json()#

Convert to a JSON string.

Return type:: str
Returns:: JSON string

JointSpectra#

class JointSpectra(data: Dict[str, JointSFS | ndarray], pop_names: List[str] = None)[source]#

Bases: _DictSpectraSerialization, AbstractSpectra

A collection of joint (multi-population) site-frequency spectra keyed by type, the multi-population analogue of Spectra. This is the return type of parse() when the parser is given populations, with one JointSFS per stratification type (or a single all type when no stratifications are used).

Methods:

`__init__`	Construct from a dictionary of joint spectra.
`marginalize`	Marginalize every joint SFS onto a subset of populations (see `JointSFS.marginalize()`).
`to_dict`	Get the joint spectra as a dictionary keyed by type.
`from_file`	Load from file.
`from_json`	Load from a JSON string.
`to_file`	Save to file (in JSON format).
`to_json`	Convert to a JSON string.

__init__(data: Dict[str, JointSFS | ndarray], pop_names: List[str] = None)[source]#

Construct from a dictionary of joint spectra.

Parameters:

data (Dict[str, Union[JointSFS, ndarray]]) – Dictionary of joint spectra (or plain arrays) keyed by type.
pop_names (List[str]) – Optional names of the populations (one per axis). Only used for entries given as plain arrays; entries already given as JointSFS keep their own names.

data: Dict[str, JointSFS]#: The joint spectra keyed by type.

property types: List[str]#: The types.

property pop_names: List[str]#

Names of the populations (one per axis), taken from the first type.

Raises:: ValueError – If the collection is empty.

property n_pops: int#

Number of populations (dimensions of each joint SFS).

Raises:: ValueError – If the collection is empty.

property shape: Tuple[int, ...]#

Shape of each joint SFS.

Raises:: ValueError – If the collection is empty.

property all: JointSFS#

The all type, equal to the sum of all joint spectra.

Raises:: ValueError – If the collection is empty.

marginalize(pops: Sequence[int])[source]#

Marginalize every joint SFS onto a subset of populations (see JointSFS.marginalize()).

Parameters:: pops (Sequence[int]) – The population indices to keep, in the desired axis order.
Return type:: JointSpectra
Returns:: Marginalized joint spectra.

to_dict()[source]#

Get the joint spectra as a dictionary keyed by type.

Return type:: Dict[str, JointSFS]
Returns:: Dictionary of joint spectra keyed by type.

classmethod from_file(file: str)#

Load from file.

Parameters:: file (str) – File path.
Return type:: Self
Returns:: Spectra.

classmethod from_json(json: str)#

Load from a JSON string.

Parameters:: json (str) – JSON string.
Return type:: Self
Returns:: Spectra.
Raises:: ValueError – If the payload does not decode to this class.

property is_empty: bool#

Returns:: True if nothing was counted.

to_file(file: str)#

Save to file (in JSON format).

Parameters:: file (str) – File path.
Return type:: None

to_json()#

Convert to a JSON string.

Return type:: str
Returns:: JSON string.

TwoSFS#

class TwoSFS(data: ndarray | list)[source]#

Bases: AbstractSpectrum

A 2-dimensional site-frequency spectrum, i.e. a square matrix whose entry (i, j) relates the number of derived alleles at a pair of frequency classes i and j of a single population (for example the second moment of the SFS). For the joint spectrum across populations, which is generally rectangular, see JointSFS.

Methods:

`__init__`	Construct from a data matrix.
`is_folded`	Check if the 2-SFS is folded.
`fold`	Fold the 2-SFS by adding up `i` and `n - i` for both axes.
`symmetrize`	Symmetrize the 2-SFS so that `i, j` and `j, i` are the same.
`interior`	The interior (segregating) block of the 2-SFS: pairs for which both sites carry between `1` and `n - 1` derived alleles, i.e. both are polymorphic.
`cov`	The class-resolved branch-length covariance `Cov(L_i, L_j)` between two linked sites, returned over the segregating interior of a full-size `TwoSFS` (monomorphic bins zeroed).
`corr`	The class-resolved branch-length correlation corresponding to `cov()`: `R[i, j] = Cov(L_i, L_j) / sqrt(Var(L_i) Var(L_j))` over the segregating interior, standardized by the branch-length variances (the interior diagonal of the full-spectrum covariance).
`fpmi`	The frequency pointwise mutual information (fPMI) of two linked sites, over the segregating interior.
`fill_monomorphic`	Fill the monomorphic entries (first and last row and column) of the 2-SFS.
`plot`	Plot the 2-SFS as a heatmap.
`plot_surface`	Plot the 2-SFS as a surface.
`mask_diagonal`	Mask both the primary and secondary diagonal entries of the 2-SFS matrix.
`get_max_abs`	Get the maximum absolute entry of the 2-SFS matrix.
`mask_upper`	Mask the upper triangular entries of the 2-SFS matrix.
`copy`	Create a deep copy.
`from_file`	Load from file.
`from_json`	Load from a JSON string.
`to_file`	Save to file (in JSON format).
`to_json`	Convert to a JSON string.

__init__(data: ndarray | list)[source]#

Construct from a data matrix.

Parameters:: data (ndarray | list) – A square 2-dimensional array.
Raises:: ValueError – If the data is not a square 2-dimensional array.

n: int = None#: The matrix dimension along one axis (the sample size plus one, i.e. the number of derived-count bins).

w: int = None#: The width of one folded half.

data: ndarray#: The 2-SFS matrix.

is_folded()[source]#

Check if the 2-SFS is folded.

Return type:: bool
Returns:: Whether the 2-SFS is folded.

fold()[source]#

Fold the 2-SFS by adding up i and n - i for both axes. Note that this only makes sense for counts or frequencies.

Return type:: TwoSFS
Returns:: Folded 2-SFS.

symmetrize()[source]#

Symmetrize the 2-SFS so that i, j and j, i are the same.

Return type:: TwoSFS
Returns:: Symmetric 2-SFS.

interior(normalize: bool = False)[source]#

The interior (segregating) block of the 2-SFS: pairs for which both sites carry between 1 and n - 1 derived alleles, i.e. both are polymorphic. The two monomorphic bins (all-ancestral and all-derived) are excluded, matching the polymorphic block of the standard site-frequency spectrum.

Parameters:: normalize (bool) – If True, return the conditional joint distribution P(i, j | both polymorphic) instead of the raw pair counts.
Return type:: ndarray
Returns:: The (n - 1) x (n - 1) interior block, or its normalization.
Raises:: ValueError – If normalize is requested but the interior is empty.

cov()[source]#

The class-resolved branch-length covariance Cov(L_i, L_j) between two linked sites, returned over the segregating interior of a full-size TwoSFS (monomorphic bins zeroed). Entry (i, j) is the covariance of the branch lengths subtending i and j derived alleles at the two loci; a positive entry means the two classes co-vary positively across linked sites. It is the deviation of the full joint class distribution from independence, P(i, j) - P(i) P(j), normalized over all pairs so that P(i) is the true site-frequency spectrum. It therefore matches PhaseGen’s sfs2.mean - outer(sfs.mean) (up to the constant mutational scale) and reproduces the multiple-merger signal, unlike a polymorphic-only normalization.

This requires the monomorphic-site counts: parse an all-sites input with the monomorphic (invariant) sites specified. With only a TargetSiteCounter (extrapolated, not real, monomorphic counts) the result is approximate and can be unreliable, as the interior residual is hypersensitive to the target-site count; use a real all-sites input for an accurate covariance, or fpmi() for a monomorphic-free statistic on polymorphic-only data. For a valid two-SFS the interior diagonal is a nonnegative branch-length variance; an arbitrary (non-coalescent) input can give a non-positive-semidefinite result.

Return type:: TwoSFS
Returns:: The branch-length covariance over the interior, embedded in a full-size TwoSFS.
Raises:: ValueError – If the 2-SFS is empty, non-finite, or carries no monomorphic-involving pairs.

corr()[source]#

The class-resolved branch-length correlation corresponding to cov(): R[i, j] = Cov(L_i, L_j) / sqrt(Var(L_i) Var(L_j)) over the segregating interior, standardized by the branch-length variances (the interior diagonal of the full-spectrum covariance). For a valid (coalescent) two-SFS this is a proper correlation, matching PhaseGen’s branch-length correlation, with entries in [-1, 1] and a unit diagonal. Like cov(), it requires the monomorphic-site counts and is only approximate under a TargetSiteCounter (see cov(); fpmi() needs no monomorphic sites). Classes with negligible branch-length variance are returned as zero, and because the underlying quantity is a cross-covariance between the two loci (not a within-locus covariance matrix), the entries are clipped to [-1, 1] as a safeguard for arbitrary input.

Return type:: TwoSFS
Returns:: The branch-length correlation over the interior, embedded in a full-size TwoSFS.
Raises:: ValueError – If the 2-SFS is empty, non-finite, or carries no monomorphic-involving pairs.

fpmi()[source]#

The frequency pointwise mutual information (fPMI) of two linked sites, over the segregating interior. For the normalized joint distribution p(i, j) of the derived-allele classes of paired polymorphic sites, entry (i, j) is log[p(i, j) / (p(i) p(j))], the log-ratio of the observed joint to the product of its marginals: positive where classes i and j co-occur more often than under independence, negative where they co-occur less, and zero for independent loci.

Unlike cov() / corr(), fPMI is a ratio computed purely from the polymorphic interior, so it needs no monomorphic sites and is exactly invariant to them: a polymorphic-only (SNP) spectrum gives the same result as the all-sites spectrum. It is the statistic of Fenton, Rice, Novembre and Desai (2025, Genetics 229(4):iyaf023, https://doi.org/10.1093/genetics/iyaf023) for detecting departures from Kingman coalescence; its low-frequency associations shift toward positive under multiple-merger genealogies.

Return type:: TwoSFS
Returns:: The fPMI over the interior, embedded in a full-size TwoSFS; classes absent from the spectrum are returned as zero.
Raises:: ValueError – If the 2-SFS is empty or non-finite, or has no polymorphic pairs.

fill_monomorphic(fill_value=nan)[source]#

Fill the monomorphic entries (first and last row and column) of the 2-SFS.

Parameters:: fill_value – Value to fill the monomorphic entries with.
Return type:: TwoSFS
Returns:: 2-SFS

plot(ax: plt.Axes = None, title: str = None, max_abs: float = None, log_scale: bool = False, cbar_kws: Dict = None, show: bool = True)[source]#

Plot the 2-SFS as a heatmap.

Parameters:

ax (plt.Axes) – Axes to plot on.
title (str) – Title of the plot.
max_abs (float) – Maximum absolute value to plot.
log_scale (bool) – Whether to use a logarithmic scale.
cbar_kws (Dict) – Keyword arguments for the color bar.
show (bool) – Whether to show the plot.

Return type:

plt.Axes

Returns:

Axes.

plot_surface(ax: plt.Axes = None, title: str = None, max_abs: float = None, vmin: float = None, vmax: float = None, show: bool = True)[source]#

Plot the 2-SFS as a surface.

Parameters:

ax (plt.Axes) – Axes to plot on.
title (str) – Title of the plot.
max_abs (float) – Maximum absolute value to plot.
vmin (float) – Minimum value to plot.
vmax (float) – Maximum value to plot.
show (bool) – Whether to show the plot.

Return type:

plt.Axes

Returns:

Axes.

mask_diagonal(fill_value=nan)[source]#

Mask both the primary and secondary diagonal entries of the 2-SFS matrix.

The primary diagonal runs from the top-left to the bottom-right, and the secondary diagonal runs from the top-right to the bottom-left.

Parameters:: fill_value – The value to fill the diagonal entries with.
Return type:: TwoSFS
Returns:: A new 2-SFS with both diagonals masked.

get_max_abs()[source]#

Get the maximum absolute entry of the 2-SFS matrix.

Return type:: float
Returns:: The maximum absolute entry.

mask_upper(fill_value=nan)[source]#

Mask the upper triangular entries of the 2-SFS matrix.

Parameters:: fill_value – The value to fill the upper triangular entries with.
Return type:: TwoSFS
Returns:: A new 2-SFS with upper triangular entries masked.

copy()#

Create a deep copy.

Return type:: Self
Returns:: Deep copy.

classmethod from_file(file: str)#

Load from file.

Parameters:: file (str) – File path.
Return type:: Self
Returns:: Spectrum

classmethod from_json(json: str)#

Load from a JSON string.

Parameters:: json (str) – JSON string.
Return type:: Self
Returns:: Spectrum
Raises:: ValueError – If the payload does not decode to this class.

property n_sites: float#: The total number of sites, i.e. the sum of all entries.

property shape: Tuple[int, ...]#: The shape of the underlying array.

to_file(file: str)#

Save to file (in JSON format).

Parameters:: file (str) – File path.
Return type:: None

to_json()#

Convert to a JSON string.

Return type:: str
Returns:: JSON string

TwoLocusSFS#

class TwoLocusSFS(data: ndarray | list)[source]#

Bases: TwoSFS

The two-locus site-frequency spectrum under recombination: a square matrix whose entry (i, j) is the expected product of the branch length subtending i samples at locus 0 and j samples at locus 1, for two loci separated by a given recombination rate r. It interpolates between the within-tree cross-moment of the SFS at r = 0 (fully linked) and the outer product of the marginal SFS as r tends to infinity (independent loci).

Methods:

`__init__`	Construct from a data matrix.
`copy`	Create a deep copy.
`corr`	The class-resolved branch-length correlation corresponding to `cov()`: `R[i, j] = Cov(L_i, L_j) / sqrt(Var(L_i) Var(L_j))` over the segregating interior, standardized by the branch-length variances (the interior diagonal of the full-spectrum covariance).
`cov`	The class-resolved branch-length covariance `Cov(L_i, L_j)` between two linked sites, returned over the segregating interior of a full-size `TwoSFS` (monomorphic bins zeroed).
`fill_monomorphic`	Fill the monomorphic entries (first and last row and column) of the 2-SFS.
`fold`	Fold the 2-SFS by adding up `i` and `n - i` for both axes.
`fpmi`	The frequency pointwise mutual information (fPMI) of two linked sites, over the segregating interior.
`from_file`	Load from file.
`from_json`	Load from a JSON string.
`get_max_abs`	Get the maximum absolute entry of the 2-SFS matrix.
`interior`	The interior (segregating) block of the 2-SFS: pairs for which both sites carry between `1` and `n - 1` derived alleles, i.e. both are polymorphic.
`is_folded`	Check if the 2-SFS is folded.
`mask_diagonal`	Mask both the primary and secondary diagonal entries of the 2-SFS matrix.
`mask_upper`	Mask the upper triangular entries of the 2-SFS matrix.
`plot`	Plot the 2-SFS as a heatmap.
`plot_surface`	Plot the 2-SFS as a surface.
`symmetrize`	Symmetrize the 2-SFS so that `i, j` and `j, i` are the same.
`to_file`	Save to file (in JSON format).
`to_json`	Convert to a JSON string.

__init__(data: ndarray | list)#

Construct from a data matrix.

Parameters:: data (ndarray | list) – A square 2-dimensional array.
Raises:: ValueError – If the data is not a square 2-dimensional array.

copy()#

Create a deep copy.

Return type:: Self
Returns:: Deep copy.

corr()#

Return type:: TwoSFS
Returns:: The branch-length correlation over the interior, embedded in a full-size TwoSFS.
Raises:: ValueError – If the 2-SFS is empty, non-finite, or carries no monomorphic-involving pairs.

cov()#

Return type:: TwoSFS
Returns:: The branch-length covariance over the interior, embedded in a full-size TwoSFS.
Raises:: ValueError – If the 2-SFS is empty, non-finite, or carries no monomorphic-involving pairs.

fill_monomorphic(fill_value=nan)#

Fill the monomorphic entries (first and last row and column) of the 2-SFS.

Parameters:: fill_value – Value to fill the monomorphic entries with.
Return type:: TwoSFS
Returns:: 2-SFS

fold()#

Fold the 2-SFS by adding up i and n - i for both axes. Note that this only makes sense for counts or frequencies.

Return type:: TwoSFS
Returns:: Folded 2-SFS.

fpmi()#

Return type:: TwoSFS
Returns:: The fPMI over the interior, embedded in a full-size TwoSFS; classes absent from the spectrum are returned as zero.
Raises:: ValueError – If the 2-SFS is empty or non-finite, or has no polymorphic pairs.

classmethod from_file(file: str)#

Load from file.

Parameters:: file (str) – File path.
Return type:: Self
Returns:: Spectrum

classmethod from_json(json: str)#

Load from a JSON string.

Parameters:: json (str) – JSON string.
Return type:: Self
Returns:: Spectrum
Raises:: ValueError – If the payload does not decode to this class.

get_max_abs()#

Get the maximum absolute entry of the 2-SFS matrix.

Return type:: float
Returns:: The maximum absolute entry.

interior(normalize: bool = False)#

Parameters:: normalize (bool) – If True, return the conditional joint distribution P(i, j | both polymorphic) instead of the raw pair counts.
Return type:: ndarray
Returns:: The (n - 1) x (n - 1) interior block, or its normalization.
Raises:: ValueError – If normalize is requested but the interior is empty.

is_folded()#

Check if the 2-SFS is folded.

Return type:: bool
Returns:: Whether the 2-SFS is folded.

mask_diagonal(fill_value=nan)#

Mask both the primary and secondary diagonal entries of the 2-SFS matrix.

The primary diagonal runs from the top-left to the bottom-right, and the secondary diagonal runs from the top-right to the bottom-left.

Parameters:: fill_value – The value to fill the diagonal entries with.
Return type:: TwoSFS
Returns:: A new 2-SFS with both diagonals masked.

mask_upper(fill_value=nan)#

Mask the upper triangular entries of the 2-SFS matrix.

Parameters:: fill_value – The value to fill the upper triangular entries with.
Return type:: TwoSFS
Returns:: A new 2-SFS with upper triangular entries masked.

n: int = None#: The matrix dimension along one axis (the sample size plus one, i.e. the number of derived-count bins).

property n_sites: float#: The total number of sites, i.e. the sum of all entries.

plot(ax: plt.Axes = None, title: str = None, max_abs: float = None, log_scale: bool = False, cbar_kws: Dict = None, show: bool = True)#

Plot the 2-SFS as a heatmap.

Parameters:

ax (plt.Axes) – Axes to plot on.
title (str) – Title of the plot.
max_abs (float) – Maximum absolute value to plot.
log_scale (bool) – Whether to use a logarithmic scale.
cbar_kws (Dict) – Keyword arguments for the color bar.
show (bool) – Whether to show the plot.

Return type:

plt.Axes

Returns:

Axes.

plot_surface(ax: plt.Axes = None, title: str = None, max_abs: float = None, vmin: float = None, vmax: float = None, show: bool = True)#

Plot the 2-SFS as a surface.

Parameters:

ax (plt.Axes) – Axes to plot on.
title (str) – Title of the plot.
max_abs (float) – Maximum absolute value to plot.
vmin (float) – Minimum value to plot.
vmax (float) – Maximum value to plot.
show (bool) – Whether to show the plot.

Return type:

plt.Axes

Returns:

Axes.

property shape: Tuple[int, ...]#: The shape of the underlying array.

symmetrize()#

Symmetrize the 2-SFS so that i, j and j, i are the same.

Return type:: TwoSFS
Returns:: Symmetric 2-SFS.

to_file(file: str)#

Save to file (in JSON format).

Parameters:: file (str) – File path.
Return type:: None

to_json()#

Convert to a JSON string.

Return type:: str
Returns:: JSON string

w: int = None#: The width of one folded half.

data: ndarray#: The 2-SFS matrix.

TwoSpectra#

class TwoSpectra(data: Dict[str, TwoSFS | ndarray])[source]#

Bases: _DictSpectraSerialization, AbstractSpectra

A collection of two-site (two-locus) site-frequency spectra keyed by type, the two-dimensional analogue of Spectra. This is the return type of parse() when the two-SFS is parsed with stratifications, holding one TwoSFS per stratification type. Because the two-SFS pairs sites, stratified parsing counts only within-stratum pairs, so summing the per-type spectra does not in general recover the unstratified two-SFS (cross-stratum pairs are not counted).

Methods:

`__init__`	Construct from a dictionary of two-site spectra.
`from_file`	Load from file.
`from_json`	Load from a JSON string.
`to_file`	Save to file (in JSON format).
`to_json`	Convert to a JSON string.
`to_dict`	Get the two-site spectra as a dictionary keyed by type.

__init__(data: Dict[str, TwoSFS | ndarray])[source]#

Construct from a dictionary of two-site spectra.

Parameters:: data (Dict[str, Union[TwoSFS, ndarray]]) – Dictionary of two-site spectra (or plain square arrays) keyed by type.

data: Dict[str, TwoSFS]#: The two-site spectra keyed by type.

property types: List[str]#: The types.

property shape: Tuple[int, ...]#

Shape of each two-site SFS.

Raises:: ValueError – If the collection is empty.

property all: TwoSFS#

The all type, equal to the sum of the per-type two-site spectra (the within-stratum pairs pooled over strata; cross-stratum pairs remain uncounted).

Raises:: ValueError – If the collection is empty.

classmethod from_file(file: str)#

Load from file.

Parameters:: file (str) – File path.
Return type:: Self
Returns:: Spectra.

classmethod from_json(json: str)#

Load from a JSON string.

Parameters:: json (str) – JSON string.
Return type:: Self
Returns:: Spectra.
Raises:: ValueError – If the payload does not decode to this class.

property is_empty: bool#

Returns:: True if nothing was counted.

to_file(file: str)#

Save to file (in JSON format).

Parameters:: file (str) – File path.
Return type:: None

to_json()#

Convert to a JSON string.

Return type:: str
Returns:: JSON string.

to_dict()[source]#

Get the two-site spectra as a dictionary keyed by type.

Return type:: Dict[str, TwoSFS]
Returns:: Dictionary of two-site spectra keyed by type.

AbstractSpectrum#

class AbstractSpectrum[source]#

Bases: ABC

Abstract base class for site-frequency spectrum containers.

A concrete spectrum wraps a numpy array in data: the one-dimensional Spectrum, the two-dimensional TwoSFS (with its TwoLocusSFS specialization), and the multi-population JointSFS. This base supplies the shared array interface and JSON serialization; subclasses add dimension-specific behaviour such as folding and plotting.

Methods:

`copy`	Create a deep copy.
`to_file`	Save to file (in JSON format).
`to_json`	Convert to a JSON string.
`from_file`	Load from file.
`from_json`	Load from a JSON string.
`plot`	Plot the spectrum.

data: ndarray#: The underlying array holding the spectrum.

property shape: Tuple[int, ...]#: The shape of the underlying array.

property n_sites: float#: The total number of sites, i.e. the sum of all entries.

copy()[source]#

Create a deep copy.

Return type:: Self
Returns:: Deep copy.

to_file(file: str)[source]#

Save to file (in JSON format).

Parameters:: file (str) – File path.
Return type:: None

to_json()[source]#

Convert to a JSON string.

Return type:: str
Returns:: JSON string

classmethod from_file(file: str)[source]#

Load from file.

Parameters:: file (str) – File path.
Return type:: Self
Returns:: Spectrum

classmethod from_json(json: str)[source]#

Load from a JSON string.

Parameters:: json (str) – JSON string.
Return type:: Self
Returns:: Spectrum
Raises:: ValueError – If the payload does not decode to this class.

abstractmethod plot(*args, **kwargs)[source]#

Plot the spectrum.

Return type:: plt.Axes
Returns:: Axes.

AbstractSpectra#

class AbstractSpectra[source]#

Bases: ABC

Abstract base class for a collection of site-frequency spectra keyed by type, for example the SFS stratified into neutral and selected sites. The concrete backings share this interface but not their storage: Spectra is DataFrame-backed and one-dimensional, while JointSpectra is dict-backed and holds multi-population JointSFS objects. Code needing only the common collection operations can be written against this base regardless of dimensionality.

Methods:

`to_dict`	Convert to a dictionary keyed by type.
`to_file`	Save to file.

abstract property types: List[str]#: The types.

abstract property all: AbstractSpectrum#: The all type, equal to the sum of all spectra.

property is_empty: bool#

Returns:: True if nothing was counted.

abstractmethod to_dict()[source]#

Convert to a dictionary keyed by type.

Return type:: dict
Returns:: Dictionary keyed by type

abstractmethod to_file(file: str)[source]#

Save to file.

Parameters:: file (str) – File path.
Return type:: None

Spectrum classes

Contents

Spectrum classes#

Spectrum#

Spectra#

JointSFS#

JointSpectra#

TwoSFS#

TwoLocusSFS#

TwoSpectra#

AbstractSpectrum#

AbstractSpectra#