Site filtration

Site filtration#

A Filtration drops sites that violate downstream modelling assumptions. Filtrations can be applied on the fly by the Parser while it builds a spectrum, or run through the Filterer to write the retained sites to a file.

`Filterer`	Filter the input using a list of filtrations.
`Filtration`	Base class for filtering sites based on certain criteria.
`MaskedFiltration`	Filter sites based on a samples mask, where `None` selects every sample.
`SNPFiltration`	Only keep SNPs.
`SNVFiltration`	Only keep single site variants (discard indels and MNPs but keep monomorphic sites).
`PolyAllelicFiltration`	Filter out poly-allelic sites.
`AllFiltration`	Filter out all sites.
`NoFiltration`	Do not filter out any sites.
`CodingSequenceFiltration`	Filter out sites that are not in coding sequences.
`DeviantOutgroupFiltration`	Filter out sites where the major allele of the specified outgroup samples differs from the major allele of the ingroup samples.
`ExistingOutgroupFiltration`	Filter out sites for which at least `n_missing` of the specified outgroup samples have no called base.
`BiasedGCConversionFiltration`	Only retain A<->T and G<->C substitutions (which are unaffected by biased gene conversion, see [CITGB]).
`CpGFiltration`	Filter out sites whose reference base is in a CpG dinucleotide context.
`ContigFiltration`	Filter out sites that are not on the specified contigs.

Filterer#

class Filterer(source: str | os.PathLike | 'tskit.TreeSequence' | VariantReader | Iterable[Site] | None = None, output: str = None, gff: str | None = None, fasta: str | None = None, filtrations: List[Filtration] = [], info_ancestral: str = 'AA', max_sites: int = inf, seed: int | None = 0, cache: bool = True, aliases: Dict[str, List[str]] = {}, vcf: str | os.PathLike | 'tskit.TreeSequence' | VariantReader | Iterable[Site] | None = None)[source]#

Bases: MultiHandler

Filter the input using a list of filtrations.

Example usage:

import sfsutils as su

# only keep variants in coding sequences
f = su.Filterer(
    source="http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/"
        "1000_genomes_project/release/20181203_biallelic_SNV/"
        "ALL.chr21.shapeit2_integrated_v1a.GRCh38.20181129.phased.vcf.gz",
    gff="http://ftp.ensembl.org/pub/release-109/gff3/homo_sapiens/"
        "Homo_sapiens.GRCh38.109.chromosome.21.gff3.gz",
    output='sapiens.chr21.coding.vcf.gz',
    filtrations=[su.CodingSequenceFiltration()],
    aliases=dict(chr21=['21'])
)

f.filter()

Methods:

`__init__`	Create a new filter instance.
`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.
`is_filtered`	Whether the given variant is kept.
`filter`	Filter the input.

__init__(source: str | os.PathLike | 'tskit.TreeSequence' | VariantReader | Iterable[Site] | None = None, output: str = None, gff: str | None = None, fasta: str | None = None, filtrations: List[Filtration] = [], info_ancestral: str = 'AA', max_sites: int = inf, seed: int | None = 0, cache: bool = True, aliases: Dict[str, List[str]] = {}, vcf: str | os.PathLike | 'tskit.TreeSequence' | VariantReader | Iterable[Site] | None = None)[source]#

Create a new filter instance.

Parameters:

source (str | os.PathLike | ‘tskit.TreeSequence’ | VariantReader | Iterable[Site] | None) – The variant source: a VCF file (gzipped or a URL), 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 / iterable of sites. Read through the same streamed site interface as all handlers.
output (str) – The output file.
gff (str | None) – The GFF file, possibly gzipped or a URL. This argument is required for some filtrations.
fasta (str | None) – The FASTA reference file, possibly gzipped or a URL. This argument is required for filtrations that depend on the reference sequence (e.g. base context).
filtrations (List[Filtration]) – The filtrations.
info_ancestral (str) – The info field for the ancestral allele.
max_sites (int) – The maximum number of sites to process.
seed (int | None) – The 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']}.
vcf (str | os.PathLike | ‘tskit.TreeSequence’ | VariantReader | Iterable[Site] | None) – Deprecated alias for source, kept for backward compatibility. Provide either source or vcf, not both.

Raises:

ValueError – If max_sites is not positive.

filtrations: List[Filtration]#: The filtrations.

output: str#: The output file.

n_filtered: int#: The number of sites that did not pass the filters.

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

is_filtered(variant: Site)[source]#

Whether the given variant is kept.

Parameters:: variant (Site) – The variant to check.
Return type:: bool
Returns:: True if the variant is kept, False otherwise.

filter()[source]#: Filter the input.

Filtration#

class Filtration[source]#

Bases: ABC

Base class for filtering sites based on certain criteria.

Methods:

`__init__`	Initialize filtration.
`filter_site`	Filter site.

n_filtered: int = 0#: The number of sites that didn’t pass the filter.

__init__()[source]#: Initialize filtration.

abstractmethod filter_site(variant: Site)[source]#

Filter site.

Parameters:: variant (Site) – The variant to filter.
Return type:: bool
Returns:: True if the variant should be kept, False otherwise.

MaskedFiltration#

class MaskedFiltration(use_parser: bool = True, include_samples: List[str] | None = None, exclude_samples: List[str] | None = None)[source]#

Bases: Filtration, ABC

Filter sites based on a samples mask, where None selects every sample.

Where the input carries no samples at all, as a sites-only VCF or store does, the verdict is taken from the declared alleles, with one warning on setup.

Methods:

`__init__`	Create a new filtration instance.
`filter_site`	Filter site.

__init__(use_parser: bool = True, include_samples: List[str] | None = None, exclude_samples: List[str] | None = None)[source]#

Create a new filtration instance.

Parameters:

use_parser (bool) – Whether to use the samples mask from the parser, if used together with parser.
include_samples (Optional[List[str]]) – The samples to include, defaults to all samples.
exclude_samples (Optional[List[str]]) – The samples to exclude, defaults to no samples.

use_parser: bool#: Whether to use the samples mask from the parser, if used together with parser.

include_samples: List[str] | None#: The samples to include.

exclude_samples: List[str] | None#: The samples to exclude.

abstractmethod filter_site(variant: Site)#

Filter site.

Parameters:: variant (Site) – The variant to filter.
Return type:: bool
Returns:: True if the variant should be kept, False otherwise.

n_filtered: int = 0#: The number of sites that didn’t pass the filter.

SNPFiltration#

class SNPFiltration(use_parser: bool = True, include_samples: List[str] | None = None, exclude_samples: List[str] | None = None)[source]#

Bases: MaskedFiltration

Only keep SNPs. Note that this entails discarding mono-morphic sites, monomorphism being judged from the alleles the included samples actually carry rather than from the ALT field.

Methods:

`filter_site`	Filter site.
`__init__`	Create a new filtration instance.

filter_site(variant: Site)[source]#

Filter site.

Parameters:: variant (Site) – The variant to filter.
Return type:: bool
Returns:: True if the variant is an SNP that is polymorphic among the included samples, False otherwise. Poly-allelic SNPs are retained; use PolyAllelicFiltration to drop those.

__init__(use_parser: bool = True, include_samples: List[str] | None = None, exclude_samples: List[str] | None = None)#

Create a new filtration instance.

Parameters:

use_parser (bool) – Whether to use the samples mask from the parser, if used together with parser.
include_samples (Optional[List[str]]) – The samples to include, defaults to all samples.
exclude_samples (Optional[List[str]]) – The samples to exclude, defaults to no samples.

n_filtered: int = 0#: The number of sites that didn’t pass the filter.

use_parser: bool#: Whether to use the samples mask from the parser, if used together with parser.

include_samples: List[str] | None#: The samples to include.

exclude_samples: List[str] | None#: The samples to exclude.

SNVFiltration#

class SNVFiltration[source]#

Bases: Filtration

Only keep single site variants (discard indels and MNPs but keep monomorphic sites).

Methods:

`filter_site`	Filter site.
`__init__`	Initialize filtration.

filter_site(variant: Site)[source]#

Filter site.

Parameters:: variant (Site) – The variant to filter.
Return type:: bool
Returns:: True if the variant is kept, False otherwise.

__init__()#: Initialize filtration.

n_filtered: int = 0#: The number of sites that didn’t pass the filter.

PolyAllelicFiltration#

class PolyAllelicFiltration(use_parser: bool = True, include_samples: List[str] | None = None, exclude_samples: List[str] | None = None)[source]#

Bases: MaskedFiltration

Filter out poly-allelic sites.

Methods:

`filter_site`	Filter site.
`__init__`	Create a new filtration instance.

filter_site(variant: Site)[source]#

Filter site. A site is poly-allelic where three or more distinct alleles are called among the included samples, which are all of them unless include_samples / exclude_samples or a parser’s populations restrict them. An alternate allele that the ALT field declares but no included sample carries therefore does not make a site poly-allelic.

Parameters:: variant (Site) – The variant to filter.
Return type:: bool
Returns:: True if the variant is not poly-allelic, False otherwise.

__init__(use_parser: bool = True, include_samples: List[str] | None = None, exclude_samples: List[str] | None = None)#

Create a new filtration instance.

Parameters:

use_parser (bool) – Whether to use the samples mask from the parser, if used together with parser.
include_samples (Optional[List[str]]) – The samples to include, defaults to all samples.
exclude_samples (Optional[List[str]]) – The samples to exclude, defaults to no samples.

n_filtered: int = 0#: The number of sites that didn’t pass the filter.

use_parser: bool#: Whether to use the samples mask from the parser, if used together with parser.

include_samples: List[str] | None#: The samples to include.

exclude_samples: List[str] | None#: The samples to exclude.

AllFiltration#

class AllFiltration[source]#

Bases: Filtration

Filter out all sites. Only useful for testing purposes.

Methods:

`filter_site`	Filter site.
`__init__`	Initialize filtration.

filter_site(variant: Site)[source]#

Filter site.

Parameters:: variant (Site) – The variant to filter.
Return type:: bool
Returns:: False.

__init__()#: Initialize filtration.

n_filtered: int = 0#: The number of sites that didn’t pass the filter.

NoFiltration#

class NoFiltration[source]#

Bases: Filtration

Do not filter out any sites. Only useful for testing purposes.

Methods:

`filter_site`	Filter site.
`__init__`	Initialize filtration.

filter_site(variant: Site)[source]#

Filter site.

Parameters:: variant (Site) – The variant to filter.
Return type:: bool
Returns:: True.

__init__()#: Initialize filtration.

n_filtered: int = 0#: The number of sites that didn’t pass the filter.

CodingSequenceFiltration#

class CodingSequenceFiltration[source]#

Bases: Filtration

Filter out sites that are not in coding sequences. This filter should find frequent use when parsing spectra for which only sites in coding sequences should be considered. By using it, the annotation and parsing of unnecessary sites can be avoided which increases the speed. Note that we assume here that within contigs, sites in the GFF file are sorted by position in ascending order.

For this filtration to work, we require a GFF file (passed to Parser or Filterer).

Methods:

`__init__`	Create a new filtration instance.
`filter_site`	Filter site by whether it is in a coding sequence.

__init__()[source]#: Create a new filtration instance.

cd: Series | None#: The coding sequence enclosing the current variant or the closest one downstream.

n_processed: int#: The number of processed sites.

filter_site(v: Site)[source]#

Filter site by whether it is in a coding sequence.

Parameters:: v (Site) – The variant to filter.
Return type:: bool
Returns:: True if the variant is in a coding sequence, False otherwise.

n_filtered: int = 0#: The number of sites that didn’t pass the filter.

DeviantOutgroupFiltration#

class DeviantOutgroupFiltration(outgroups: List[str], ingroups: List[str] = None, strict_mode: bool = True, retain_monomorphic: bool = True)[source]#

Bases: Filtration

Filter out sites where the major allele of the specified outgroup samples differs from the major allele of the ingroup samples.

Methods:

`__init__`	Construct DeviantOutgroupFiltration.
`filter_site`	Filter site.

__init__(outgroups: List[str], ingroups: List[str] = None, strict_mode: bool = True, retain_monomorphic: bool = True)[source]#

Construct DeviantOutgroupFiltration.

Parameters:

outgroups (List[str]) – The name of the outgroup samples to consider.
ingroups (List[str]) – The name of the ingroup samples to consider, defaults to all samples but the outgroups.
strict_mode (bool) – Whether to filter out sites where no outgroup sample is present, defaults to True.
retain_monomorphic (bool) – Whether to retain monomorphic sites, defaults to True, which is faster.

ingroups: List[str] | None#: The ingroup samples.

outgroups: List[str]#: The outgroup samples.

strict_mode: bool#: Whether to filter out sites where no outgroup sample is present.

retain_monomorphic: bool#: Whether to retain monomorphic sites.

samples: ndarray | None#: The samples found in the input.

ingroup_mask: ndarray | None#: The ingroup mask.

outgroup_mask: ndarray | None#: The outgroup mask.

filter_site(variant: Site)[source]#

Filter site.

Parameters:: variant (Site) – The variant to filter.
Return type:: bool
Returns:: True if the variant should be kept, False otherwise.

n_filtered: int = 0#: The number of sites that didn’t pass the filter.

ExistingOutgroupFiltration#

class ExistingOutgroupFiltration(outgroups: List[str], n_missing: int = 1)[source]#

Bases: Filtration

Filter out sites for which at least n_missing of the specified outgroup samples have no called base.

Methods:

`__init__`	Construct ExistingOutgroupFiltration.
`filter_site`	Filter site.

__init__(outgroups: List[str], n_missing: int = 1)[source]#

Construct ExistingOutgroupFiltration.

Parameters:

outgroups (List[str]) – The names of the outgroup samples considered.
n_missing (int) – The number of outgroup samples that need to be missing to fail the filter.

outgroups: List[str]#: The outgroup samples.

n_missing: int#: Minimum number of missing outgroups required to filter out a site.

samples: ndarray | None#: The samples found in the input.

outgroup_mask: ndarray | None#: The outgroup mask.

filter_site(variant: Site)[source]#

Filter site.

Parameters:: variant (Site) – The variant to filter.
Return type:: bool
Returns:: True if the variant should be kept, False otherwise.

n_filtered: int = 0#: The number of sites that didn’t pass the filter.

BiasedGCConversionFiltration#

class BiasedGCConversionFiltration[source]#

Bases: Filtration

Only retain A<->T and G<->C substitutions (which are unaffected by biased gene conversion, see [CITGB]).

Mono-allelic sites are always retained, and we assume sites are at most bi-allelic. Note that the number of mutational target sites is reduced by this filtration.

[CITGB] (1,2)

Pouyet et al., ‘Background selection and biased gene conversion affect more than 95% of the human genome and bias demographic inferences.’, Elife, 7:e36317, 2018

Methods:

`filter_site`	Remove bi-allelic sites that are not A<->T or G<->C mutations.
`__init__`	Initialize filtration.

filter_site(variant: Site)[source]#

Remove bi-allelic sites that are not A<->T or G<->C mutations.

Parameters:: variant (Site) – The variant to filter.
Return type:: bool
Returns:: True if the variant should be kept, False otherwise.

__init__()#: Initialize filtration.

n_filtered: int = 0#: The number of sites that didn’t pass the filter.

CpGFiltration#

class CpGFiltration[source]#

Bases: Filtration

Filter out sites whose reference base is in a CpG dinucleotide context. CpG sites are hypermutable (the cytosine is prone to deamination), so they are commonly excluded to avoid mutation-rate heterogeneity. A site is in CpG context iff:

the reference base is C and the next base on the same strand is G, or
the reference base is G and the previous base on the same strand is C.

Like CodingSequenceFiltration, this filtration requires a FASTA reference (passed to Parser or Filterer), which is used for the ±1 base lookup. Sites on a contig the FASTA carries no sequence for, and sites the FASTA sequence does not reach, cannot be typed and are kept, with one warning per contig.

Methods:

`__init__`	Create a new filtration instance.
`filter_site`	Filter site by whether its reference base is in a CpG context.

__init__()[source]#: Create a new filtration instance.

filter_site(variant: Site)[source]#

Filter site by whether its reference base is in a CpG context.

Parameters:: variant (Site) – The variant to filter.
Return type:: bool
Returns:: True if the site should be kept (not CpG), False otherwise.

n_filtered: int = 0#: The number of sites that didn’t pass the filter.

ContigFiltration#

class ContigFiltration(contigs: List[str])[source]#

Bases: Filtration

Filter out sites that are not on the specified contigs.

Methods:

`__init__`	Construct ContigFiltration.
`filter_site`	Filter site.

__init__(contigs: List[str])[source]#

Construct ContigFiltration.

Parameters:: contigs (List[str]) – The contigs to retain.

contigs: List[str]#: The contigs to retain.

filter_site(variant: Site)[source]#

Filter site.

Parameters:: variant (Site) – The variant to filter.
Return type:: bool
Returns:: True if the variant is on one of the specified contigs, False otherwise.

n_filtered: int = 0#: The number of sites that didn’t pass the filter.

Site filtration

Contents

Site filtration#

Filterer#

Filtration#

MaskedFiltration#

SNPFiltration#

SNVFiltration#

PolyAllelicFiltration#

AllFiltration#

NoFiltration#

CodingSequenceFiltration#

DeviantOutgroupFiltration#

ExistingOutgroupFiltration#

BiasedGCConversionFiltration#

CpGFiltration#

ContigFiltration#