Parsing

The Parser reads variants into a site-frequency spectrum. When the input contains no monomorphic sites, a TargetSiteCounter recovers the number of mutational target sites, so that the monomorphic bins, and the overall scale of the spectrum, are correct. Sites can be split into categories with a Stratification (see Site stratification).

Parser	Parse site-frequency spectra from VCF files, VCF-Zarr stores, or tskit tree sequences.
TargetSiteCounter	Class for counting the number of target sites when parsing an input that does not contain monomorphic sites.

Parser

class Parser(source: str | os.PathLike | 'tskit.TreeSequence' | VariantReader | None = None, n: int | Dict[str, int] | List[int] = None, pops: Dict[str, List[str]] | None = None, gff: str | None = None, fasta: str | None = None, info_ancestral: str = 'AA', info_ancestral_prob: str = 'AA_prob', skip_non_polarized: bool = True, stratifications: List[Stratification] = [], annotations: List[Annotation] = [], filtrations: List[Filtration] = None, include_samples: List[str] = None, exclude_samples: List[str] = None, max_sites: int = inf, seed: int | None = 0, cache: bool = True, aliases: Dict[str, List[str]] = {}, target_site_counter: TargetSiteCounter = None, subsample_mode: Literal['random', 'probabilistic'] = 'probabilistic', polarize_probabilistically: bool = False, two_sfs: bool = False, d: int = 1000, two_sfs_offset: int = 0, vcf: str | os.PathLike | 'tskit.TreeSequence' | VariantReader | None = None)

Bases: MultiHandler

Parse site-frequency spectra from VCF files, VCF-Zarr stores, or tskit tree sequences.

By default, the parser looks at the AA tag in the VCF file’s info field to retrieve the correct polarization. Polymorphic sites for which this tag is not well-defined are by default ignored (see skip_non_polarized).

This class also offers on-the-fly annotation of the sites such as site degeneracy and ancestral allele state. This is done by providing a list of annotations to the parser which are applied in the order they are provided.

The parser also allows to filter sites based on site properties which is done by passing a list of filtrations. By default, we filter out poly-allelic sites as sites are assumed to be at most bi-allelic.

In addition, the parser allows to stratify the SFS by providing a list of stratifications. This is useful to obtain the SFS separately for different types of sites.

To correctly determine the number of target sites when parsing an input that does not contain monomorphic sites, we can use a TargetSiteCounter. This class is used in conjunction with the parser and samples sites from the given FASTA file that are found in between variants on the same contig that were parsed from the input.

Note that we assume the sites in the input to be sorted by position in ascending order (per contig).

Example usage:

import sfsutils as su

# Parse selected and neutral SFS from human chromosome 1.
p = su.Parser(
    source="https://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/"
        "hgdp_wgs.20190516.full.chr21.vcf.gz",
    fasta="http://ftp.ensembl.org/pub/release-109/fasta/homo_sapiens/"
          "dna/Homo_sapiens.GRCh38.dna.chromosome.21.fa.gz",
    gff="http://ftp.ensembl.org/pub/release-109/gff3/homo_sapiens/"
        "Homo_sapiens.GRCh38.109.chromosome.21.gff3.gz",
    aliases=dict(chr21=['21']),  # mapping for contig names
    n=10,  # SFS sample size
    # we use a target site counter to infer the number of target sites.
    target_site_counter=su.TargetSiteCounter(
        n_samples=1000000,
        # determine number of target sites by looking at total length of coding sequences
        n_target_sites=su.Annotation.count_target_sites(
            "http://ftp.ensembl.org/pub/release-109/gff3/homo_sapiens/"
            "Homo_sapiens.GRCh38.109.chromosome.21.gff3.gz"
        )['21']
    ),
    # add degeneracy annotation for sites
    annotations=[
        su.DegeneracyAnnotation()
    ],
    filtrations=[
        # exclude non-SNPs as we infer monomorphic sites with target site counter
        su.SNPFiltration(),
        # filter out sites not in coding sequences
        su.CodingSequenceFiltration()
    ],
    # stratify by 4-fold/0-fold degeneracy
    stratifications=[su.DegeneracyStratification()],
    info_ancestral='AA_ensembl'
)

sfs = p.parse()

sfs.plot()

Methods:

__init__	Initialize the parser.
count_sites	Count the number of sites in the source.
download_file	Download a file from a URL.
download_if_url	Download the file if it is a URL.
get_aliases	Get all aliases for the given contig alias including the primary alias.
get_contig	Get the contig from the FASTA file.
get_contig_names	Get the names of the contigs in the FASTA file.
get_filename	Return the file extension of a URL.
get_pbar	Return a progress bar for the number of sites.
hash	Return a truncated SHA1 hash of a string.
is_url	Check if the given path is a URL.
load_fasta	Load a FASTA file into a dictionary.
load_vcf	Open a VCF file for streaming.
remove_overlaps	Remove overlapping coding sequences.
unzip_if_zipped	If the given file is gzipped, unzip it and return the path to the unzipped file.
parse	Parse the site-frequency spectrum from the configured source (VCF, VCF-Zarr, or tree sequence).

__init__(source: str | os.PathLike | 'tskit.TreeSequence' | VariantReader | None = None, n: int | Dict[str, int] | List[int] = None, pops: Dict[str, List[str]] | None = None, gff: str | None = None, fasta: str | None = None, info_ancestral: str = 'AA', info_ancestral_prob: str = 'AA_prob', skip_non_polarized: bool = True, stratifications: List[Stratification] = [], annotations: List[Annotation] = [], filtrations: List[Filtration] = None, include_samples: List[str] = None, exclude_samples: List[str] = None, max_sites: int = inf, seed: int | None = 0, cache: bool = True, aliases: Dict[str, List[str]] = {}, target_site_counter: TargetSiteCounter = None, subsample_mode: Literal['random', 'probabilistic'] = 'probabilistic', polarize_probabilistically: bool = False, two_sfs: bool = False, d: int = 1000, two_sfs_offset: int = 0, vcf: str | os.PathLike | 'tskit.TreeSequence' | VariantReader | None = None)

Initialize the parser.

Parameters:

• source (str | os.PathLike | ‘tskit.TreeSequence’ | VariantReader | None) – The variant source: a path to a VCF file (gzipped or a URL), a path to a VCF-Zarr store (a .vcz or .zarr directory), a tskit tree sequence (a .trees file or an in-memory tskit.TreeSequence), or a pre-built VariantReader. A non-path, non-tree-sequence source must be a VariantReader (it exposes the sample names, contig names and site count the parser needs, and supports a fresh iteration pass); a bare iterable of sites or a raw cyvcf2.VCF object is not accepted and raises a TypeError. VCF-Zarr requires the optional zarr package and tree sequences the optional tskit package.

• gff (str | None) – The path to the GFF file, possibly gzipped or a URL. This file is optional and depends on the stratifications, annotations and filtrations that are used.

• fasta (str | None) – The path to the FASTA file, possibly gzipped or a URL. This file is optional and depends on the annotations and filtrations that are used.

• n (Union[int, Dict[str, int], List[int]]) – The size of the resulting SFS. We down-sample to this number by drawing without replacement from the set of all available genotypes per site. Sites with fewer than n genotypes are skipped. For a joint (multi-population) SFS (see pops), this is the per-population sample size, given either as a single int applied to every population, a list aligned with pops (in insertion order), or a dictionary keyed by population name.

• pops (Optional[Dict[str, List[str]]]) – Mapping of population name to the list of sample names making up that population. When given, parse() returns a JointSpectra holding the joint SFS across these populations (one JointSFS per stratification type) instead of the one-dimensional Spectra. The joint SFS is obtained by down-projecting each population independently and accumulating the outer product of the per-population projections, which is the exact hypergeometric down-projection under independent sampling within populations. When None (default), a single-population, one-dimensional SFS is parsed as before. Note that pops supersedes include_samples / exclude_samples. A target_site_counter is supported: monomorphic sites map to the all-ancestral origin of the joint SFS, which is scaled to the target-site count.

• info_ancestral (str) – The tag in the INFO field that contains ancestral allele information. Consider using an ancestral allele annotation if this information is not available yet.

• info_ancestral_prob (str) – The tag in the INFO field that holds the per-site probability that info_ancestral names the true ancestral allele, used when polarize_probabilistically is set.

• skip_non_polarized (bool) – Whether to skip poly-morphic sites that are not polarized, i.e., without a valid info tag providing the ancestral allele. If False, we use the reference allele as ancestral allele (only recommended if working with folded spectra).

• stratifications (List[Stratification]) – List of stratifications to use.

• annotations (List[Annotation]) – List of annotations to use.

• filtrations (List[Filtration]) – List of filtrations to use. By default, we use PolyAllelicFiltration.

• include_samples (List[str]) – List of sample names to consider when determining the SFS. If None, all samples are used. Note that this restriction does not apply to the annotations and filtrations.

• exclude_samples (List[str]) – List of sample names to exclude when determining the SFS. If None, no samples are excluded. Note that this restriction does not apply to the annotations and filtrations.

• max_sites (int) – Maximum number of sites to parse from the input, which must be positive.

• seed (int | None) – Seed for the random number generator. Use None for no seed.

• cache (bool) – Whether to cache files downloaded from URLs.

• aliases (Dict[str, List[str]]) – Dictionary of aliases for the contigs in the input, e.g. {'chr1': ['1']}. This is used to match the contig names in the input with the contig names in the FASTA file and GFF file.

• target_site_counter (TargetSiteCounter) – The target site counter, used to recover the number of monomorphic sites when the input contains only polymorphic sites. If None, no target sites are added. It applies to the one-dimensional SFS (sampling monomorphic sites from the FASTA), the joint SFS (scaling the all-ancestral corner to the target-site count), and the unstratified two-SFS (extrapolating the monomorphic-involving pairs from the target-site count, which makes the resulting cov/corr only approximate).

• subsample_mode (Literal['random', 'probabilistic']) – The subsampling mode. For random, we draw once without replacement from the set of all available genotypes per site. For probabilistic, we add up the hypergeometric distribution for all sites. This will produce a smoother SFS, especially when a small number of sites is considered.

• polarize_probabilistically (bool) – Whether to probabilistically polarize sites. In addition to the AA tag (see info_ancestral), we use the AA_prob tag (see info_ancestral_prob) to polarize sites probabilistically. For example, if the ancestral allele is A with a probability of 0.8 and the derived allele is G, we assign 0.8 probability mass to the ancestral allele and 0.2 to the derived allele. This should enhance accuracy, especially for small datasets. Whenever the ancestral probability tag is not present, we assume a probability of 1 for the ancestral allele.

• two_sfs (bool) – Whether to parse the two-dimensional (two-site) SFS instead of the ordinary SFS. When True, parse() returns a square TwoSFS whose entry (i, j) counts pairs of sites, on the same contig and within the distance window of one another, where one site has i and the other j derived alleles (down-projected to n). As for the one-dimensional SFS, monomorphic sites are retained (contributing to the zero-frequency row and column); add a SNPFiltration to restrict the spectrum to segregating sites. Each site’s contribution is the outer product of the two per-site down-projection vectors, so it is exact when subsample_mode='random' at the full sample size and smoother under 'probabilistic'. The matrix is symmetrized. Not compatible with pops. Stratifications are supported (counting only within-stratum pairs, into a TwoSpectra). The unstratified two-SFS accepts a target_site_counter, which extrapolates the monomorphic-involving pairs from the target-site count; the resulting cov() / corr() are then only approximate, so a real all-sites input (monomorphic sites retained) is preferred for an accurate covariance, while fpmi() needs no monomorphic sites.

• d (int) – The width (in base pairs) of the distance window over which the two sites of a pair are separated when two_sfs=True. Together with two_sfs_offset it defines the window (two_sfs_offset, two_sfs_offset + d]; with the default two_sfs_offset=0 this is simply the maximum separation. Restricting the window restricts the spectrum to (approximately) linked pairs.

• two_sfs_offset (int) – The minimum genomic distance (in base pairs, exclusive) between the two sites of a pair when two_sfs=True; pairs are formed for separations in (two_sfs_offset, two_sfs_offset + d]. Defaults to 0 (pairs at any separation up to d).

• vcf (str | os.PathLike | ‘tskit.TreeSequence’ | VariantReader | None) – Deprecated alias for source, kept for backward compatibility. Provide either source or vcf, not both.

target_site_counter: TargetSiteCounter | None

The target site counter

pops: Dict[str, List[str]] | None

Mapping of population name to sample names for the joint SFS, or None for a single-population SFS

n: int | None

The per-population sample size (single population)

include_samples: List[str] | None

The list of samples to include

exclude_samples: List[str] | None

The list of samples to exclude

skip_non_polarized: bool

Whether to skip sites that are not polarized, i.e., without a valid info tag providing the ancestral allele

stratifications: List[Stratification]

List of stratifications to use

annotations: List[Annotation]

List of annotations to use

filtrations: List[Filtration]

List of filtrations to use

n_skipped: int

The number of sites that were skipped for various reasons

n_aa_prob: int

n_no_ancestral: int

The number of sites that were skipped because they had no valid ancestral allele

sfs: Dict[str, np.ndarray]

Dictionary of SFS indexed by (stratification) type. Each value is a 1-D array of length n + 1 for a single-population SFS, or a joint SFS array of shape _joint_shape when pops is given.

subsample_mode: Literal['random', 'probabilistic']

The subsampling mode

info_ancestral_prob: str

The tag in the INFO field that contains the ancestral allele probability

polarize_probabilistically: bool

Whether to probabilistically polarize sites

two_sfs: bool

Whether to parse the two-dimensional (two-site) SFS

d: int

Maximum genomic distance (bp) between the two sites of a pair for the two-SFS

two_sfs_offset: int

Minimum genomic distance (bp, exclusive) between the two sites of a pair for the two-SFS

property contig_lengths: Dict[str, int] | None

The declared length (bp) of each contig of the source: the ##contig headers of a VCF, the contig_length array of a VCF-Zarr store, or the length of the genome of a tree sequence. Only the contigs whose length is known are present, and a source declaring none gives None. The span of the observed variants understates the region on a sparsely covered contig, so this is what a spectrum extrapolating monomorphic sites should be sized against.

Returns:

The per-contig lengths, or None where the source declares none.

count_sites()

Count the number of sites in the source.

Return type:

int

Returns:

Number of sites

classmethod download_file(url: str, cache: bool = True, desc: str = 'Downloading file')

Download a file from a URL.

Parameters:

• cache (bool) – Whether to cache the file.

• url (str) – The URL to download the file from.

• desc (str) – Description for the progress bar

Return type:

str

Returns:

The path to the downloaded file.

download_if_url(path: str)

Download the file if it is a URL.

Parameters:

path (str) – The path to the file.

Return type:

str

Returns:

The path to the downloaded file or the original path.

get_aliases(contig: str)

Get all aliases for the given contig alias including the primary alias.

Parameters:

contig (str) – The contig.

Return type:

List[str]

Returns:

The aliases.

get_contig(aliases, rewind: bool = True, notify: bool = True)

Get the contig from the FASTA file. The contig is looked up in the header index and read by seeking to it, so the sites may visit the contigs in any order at the same cost.

Note that pyfaidx would be more efficient here, but there were problems when running it in parallel.

Parameters:

• aliases – The contig aliases.

• rewind (bool) – Unused, the lookup does not depend on the position of a cursor.

• notify (bool) – Unused, the lookup does not depend on the position of a cursor.

Return type:

SeqRecord

Returns:

The contig.

Raises:

LookupError – Where the FASTA carries none of the aliases.

get_contig_names()

Get the names of the contigs in the FASTA file.

Return type:

List[str]

Returns:

The contig names.

static get_filename(url: str)

Return the file extension of a URL.

Parameters:

url (str) – The URL to get the file extension from.

Returns:

The file extension.

get_pbar(desc: str = 'Processing sites', total: int | None = 0)

Return a progress bar for the number of sites.

Parameters:

• desc (str) – Description for the progress bar

• total (int | None) – Total number of items

Return type:

tqdm

Returns:

tqdm

static hash(s: str)

Return a truncated SHA1 hash of a string.

Parameters:

s (str) – The string to hash.

Return type:

str

Returns:

The SHA1 hash.

static is_url(path: str)

Check if the given path is a URL.

Parameters:

path (str) – The path to check.

Return type:

bool

Returns:

True if the path is a URL, False otherwise.

load_fasta(file: str)

Load a FASTA file into a dictionary.

Parameters:

file (str) – The path to The FASTA file path, possibly gzipped or a URL

Return type:

FastaIterator

Returns:

Iterator over the sequences.

load_vcf()

Open a VCF file for streaming.

Return type:

cyvcf2.VCF

Returns:

The VCF reader.

property n_sites: int

Get the number of sites in the input.

Returns:

Number of sites

static remove_overlaps(df: DataFrame)

Remove overlapping coding sequences.

Parameters:

df (DataFrame) – The coding sequences.

Return type:

DataFrame

Returns:

The coding sequences without overlaps.

static unzip_if_zipped(file: str)

If the given file is gzipped, unzip it and return the path to the unzipped file. If the file is not gzipped, return the path to the original file.

Parameters:

file (str) – The path to the file.

Returns:

The path to the unzipped file, or the original file if it was not gzipped.

info_ancestral: str

The tag in the INFO field that contains the ancestral allele

max_sites: int

Maximum number of sites to consider

seed: int | None

Seed for the random number generator

fasta: str

The path to the FASTA file.

cache: bool

Whether to cache files that are downloaded from URLs

source

The variant source. Alias of vcf, which is kept for backward compatibility.

vcf

The variant source (a path or a tskit TreeSequence object)

rng

Random generator instance

gff

The GFF file path

aliases

The contig mappings

parse()

Parse the site-frequency spectrum from the configured source (VCF, VCF-Zarr, or tree sequence).

The return type is fixed by the parsing mode, and this mapping is a stable part of the API. Every mode returns a collection keyed by stratification type (a single 'all' key when no stratifications were given), so the return type is consistent across modes:

• one-dimensional SFS (the default): a Spectra of one-dimensional spectra;

• multi-population (pops given): a JointSpectra of joint spectra, keyed by (sorted) stratification type;

• two-SFS (two_sfs set, with or without stratifications): a TwoSpectra of per-type two-SFS matrices. Without stratifications this holds the single 'all' entry, so the two-SFS is reached as parse()['all'].

Return type:

Spectra | JointSpectra | TwoSpectra

Returns:

A Spectra, JointSpectra, or TwoSpectra as described above.

TargetSiteCounter

class TargetSiteCounter(n_target_sites: int, n_samples: int = 100000)

Bases: object

Class for counting the number of target sites when parsing an input that does not contain monomorphic sites. This class is used in conjunction with Parser and samples sites from the given fasta file that are found in between variants on the same contig that were parsed from the input. Ideally, we obtain the SFS by parsing inputs that contain both mono- and polymorphic sites. This is because we need to know about the number of mutational opportunities for synonymous and non-synonymous sites which contain plenty of information on the strength of selection. It is recommended to use a SNPFiltration when using this class to avoid biasing the result by monomorphic sites present in the input.

Warning

This class is not compatible with stratifications based on info tags that are pre-defined in the input, as opposed to those added dynamically using the annotations argument of the parser. We also need to stratify mono-allelic sites which, in this case, won’t be present in the input so that they have no info tags when sampling from the FASTA file, and are thus ignored by the stratifications. However, using the annotations argument of the parser, the info tags the stratifications are based on are added on-the-fly, also for monomorphic sites sampled from the FASTA file.

Note

With the unstratified two-SFS (two_sfs=True) this counter extrapolates the monomorphic-involving pairs from the target-site count. That anchors the marginal only approximately, so the resulting cov() / corr are approximate (a real all-sites input is preferred); fpmi() needs no monomorphic sites. The stratified two-SFS with a counter is not supported.

Methods:

__init__	Initialize counter.
count	Count the number of target sites.

__init__(n_target_sites: int, n_samples: int = 100000)

Initialize counter.

Parameters:

• n_target_sites (int) – The total number of sites (mono- and polymorphic) that would be present in the input if it contained monomorphic sites. This number should be considerably larger than the number of polymorphic sites in the input. this value is not extremely important for downstream inference, the ratio of synonymous to non-synonymous sites being more informative, but the order of magnitude should be correct, in any case.

• n_samples (int) – The number of sites to sample from the fasta file. Many sampled sites will not be valid as they are non-coding. To obtain good estimates, a few thousand sites should be sampled per type of site (depending on the stratifications used).

n_target_sites: int | None

The total number of sites considered when parsing the input

n_samples: int

Number of samples

count()

Count the number of target sites.

Returns:

The number of target sites