Manipulating the SFS

Parsing yields Spectrum and Spectra objects. A Spectrum holds a single site-frequency spectrum; a Spectra holds a named collection, one entry per type, such as the neutral and selected spectra produced by a stratified parse. Both expose the operations needed to inspect and reshape spectra, including indexing, grouping, folding, resampling, serialisation, and plotting. The examples below build spectra directly to illustrate these operations, but in practice they are usually returned by the Parser.

import sfsutils as su
import matplotlib.pyplot as plt

# create spectra with two subtypes and two types
spectra = su.Spectra.from_spectra({
    "subtype1.type1": su.Spectrum.standard_kingman(10) * 1,
    "subtype1.type2": su.Spectrum.standard_kingman(10) * 2,
    "subtype2.type1": su.Spectrum.standard_kingman(10) * 3,
})

# plot spectra
spectra.plot();

We access types by their index from which we obtain a Spectrum object.

sfs: su.Spectrum = spectra["subtype1.type1"]

sfs.plot();

We can also use wildcards to access multiple types at once.

spectra["subtype1.*"].plot();

Grouping

To get rid of the subtypes, we can merge the spectra over the specified number of groups.

spectra.merge_groups(1).plot();

All subtypes for each type are merged into a single spectrum by adding them up.

Serialization

We can also save the spectra to a file and restore them again.

spectra.to_file("out/spectra.csv")

spectra2 = su.Spectra.from_file("out/spectra.csv")

Prefixing

Here we prefix the spectra with a string to distinguish them and then combine them into a single spectra object.

spectra.prefix('original').combine(spectra2.prefix('restored')).plot();

For a complete reference of the available methods and properties, see Spectra and Spectrum.

Folded spectra

Spectrum and Spectra objects can also be folded by collapsing the bins corresponding to the derived allele counts onto the bins corresponding to the ancestral allele counts. Folding discards information, which is particularly noticeable when beneficial mutations are present. However, folded spectra are easier to obtain, and are robust to misspecification of the ancestral state, which is often unknown. A spectrum reports whether it has been folded through its is_folded() method.

# create subplots
fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(7, 3))

# fold spectra object
spectra.fold().plot(ax=ax1)

# fold spectrum object
sfs.fold().plot(ax=ax2);