Site annotation

Site annotation#

An Annotation adds site-level information such as degeneracy, synonymy or the ancestral allele. Annotations can be applied on the fly by the Parser while it builds a spectrum, or run through the Annotator to write the annotated variants to a file.

`Annotator`	Annotate the input with the given annotations.
`Annotation`	Base class for annotations.
`DegeneracyAnnotation`	Degeneracy annotation.
`SynonymyAnnotation`	Synonymy annotation.
`AncestralAlleleAnnotation`	Base class for ancestral allele annotation.
`MaximumParsimonyAncestralAnnotation`	Annotation of ancestral alleles using maximum parsimony.
`SubstitutionModel`	Base class for substitution models.
`JCSubstitutionModel`	Jukes-Cantor substitution model.
`K2SubstitutionModel`	Kimura 2-parameter substitution model.
`SiteConfig`	Ancestral allele site configuration for a single subsample.
`SiteInfo`	Ancestral allele information on a single site.
`BaseType`	The base type, either major or minor.
`PolarizationPrior`	Base class for priors used with `MaximumLikelihoodAncestralAnnotation`.
`KingmanPolarizationPrior`	Prior based on the standard Kingman coalescent.
`AdaptivePolarizationPrior`	Adaptive prior.
`MaximumLikelihoodAncestralAnnotation`	Annotation of ancestral alleles following the probabilistic model of EST-SFS (https://doi.org/10.1534/genetics.118.301120).
`AdHocAncestralAnnotation`	Ad-hoc ancestral allele annotation using simple rules.

Annotator#

class Annotator(source: str | os.PathLike | 'tskit.TreeSequence' | VariantReader | Iterable[Site] | None = None, output: str = None, annotations: List[Annotation] = None, gff: str | None = None, fasta: str | None = None, 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

Annotate the input with the given annotations.

Example usage:

import sfsutils as su

ann = su.Annotator(
    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",
    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",
    output='sapiens.chr21.degeneracy.vcf.gz',
    annotations=[su.DegeneracyAnnotation()],
    aliases=dict(chr21=['21'])
)

ann.annotate()

Methods:

`__init__`	Create a new annotator instance.
`annotate`	Annotate the input.
`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.

__init__(source: str | os.PathLike | 'tskit.TreeSequence' | VariantReader | Iterable[Site] | None = None, output: str = None, annotations: List[Annotation] = None, gff: str | None = None, fasta: str | None = None, 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 annotator 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 path to the output file
annotations (List[Annotation]) – The annotations to apply.
gff (str | None) – The path to the GFF file, can be gzipped, urls are also supported. Required for annotations that require a GFF file.
fasta (str | None) – The path to the FASTA file, can be gzipped, urls are also supported. Required for annotations that require a FASTA file.
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. 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.
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.

output: str#: The path to the output file.

annotations: List[Annotation]#: The annotations to apply.

annotate()[source]#: Annotate the input.

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

Annotation#

class Annotation[source]#

Bases: ABC

Base class for annotations.

Methods:

`__init__`	Create a new annotation instance.
`annotate_site`	Annotate a single site.
`count_target_sites`	Count the number of target sites in a GFF file.

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

n_annotated: int#: The number of annotated sites.

abstractmethod annotate_site(variant: Site)[source]#

Annotate a single site.

Parameters:: variant (Site) – The variant to annotate.
Returns:: The annotated variant.

static count_target_sites(file: str, remove_overlaps: bool = False, contigs: List[str] = None)[source]#

Count the number of target sites in a GFF file.

Parameters:

file (str) – The path to The GFF file path, possibly gzipped or a URL
remove_overlaps (bool) – Whether to remove overlapping target sites.
contigs (List[str]) – The contigs to count the target sites for.

Return type:

Dict[str, int]

Returns:

The number of target sites per chromosome/contig.

DegeneracyAnnotation#

class DegeneracyAnnotation[source]#

Bases: Annotation

Degeneracy annotation. We annotate the degeneracy by looking at each codon for coding variants. This also annotates mono-allelic sites.

This annotation adds the info fields Degeneracy and Degeneracy_Info, which hold the degeneracy of a site (0, 2, 4) and extra information about the degeneracy, respectively. To be used with DegeneracyStratification.

For this annotation to work, we require a FASTA and GFF file (passed to Parser or Annotator).

Example usage:

import sfsutils as su

ann = su.Annotator(
    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",
    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",
    output='sapiens.chr21.degeneracy.vcf.gz',
    annotations=[su.DegeneracyAnnotation()],
    aliases=dict(chr21=['21'])
)

ann.annotate()

Methods:

`__init__`	Create a new annotation instance.
`annotate_site`	Annotate a single site.
`count_target_sites`	Count the number of target sites in a GFF file.

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

mismatches: List['cyvcf2.Variant']#: The variants that could not be annotated correctly.

n_skipped: int#: The variant that were skipped because they were not in coding regions.

errors: List['cyvcf2.Variant']#: The variants for which the codon could not be determined.

annotate_site(v: Site)[source]#

Annotate a single site.

Parameters:: v (Site) – The variant to annotate.

static count_target_sites(file: str, remove_overlaps: bool = False, contigs: List[str] = None)#

Count the number of target sites in a GFF file.

Parameters:

file (str) – The path to The GFF file path, possibly gzipped or a URL
remove_overlaps (bool) – Whether to remove overlapping target sites.
contigs (List[str]) – The contigs to count the target sites for.

Return type:

Dict[str, int]

Returns:

The number of target sites per chromosome/contig.

n_annotated: int#: The number of annotated sites.

SynonymyAnnotation#

class SynonymyAnnotation[source]#

Bases: DegeneracyAnnotation

Synonymy annotation. This class annotates a variant with the synonymous/non-synonymous status.

This annotation adds the info fields Synonymy and Synonymy_Info, which hold the synonymy (synonymous [1] or non-synonymous [0]) and the codon information, respectively. To be used with SynonymyStratification.

For this annotation to work, we require a FASTA and GFF file (passed to Parser or Annotator).

This class was tested against VEP and SnpEff and provides the same annotations in almost all cases.

Warning

Not recommended for use with Parser as we also need to annotate mono-allelic sites. Consider using DegeneracyAnnotation and DegeneracyStratification instead.

Methods:

`__init__`	Create a new annotation instance.
`mutate`	Mutate the codon at the given position to the given alternative allele.
`is_synonymous`	Check if two codons are synonymous.
`annotate_site`	Annotate a single site.
`count_target_sites`	Count the number of target sites in a GFF file.

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

vep_mismatches: List['cyvcf2.Variant']#: The number of sites that did not match with VEP.

snpeff_mismatches: List['cyvcf2.Variant']#: The number of sites that did not math with the annotation provided by SnpEff

n_vep_comparisons: int#: The number of sites that were concordant with VEP.

n_snpeff_comparisons: int#: The number of sites that were concordant with SnpEff.

static mutate(codon: str, alt: str, pos: int)[source]#

Mutate the codon at the given position to the given alternative allele.

Parameters:

codon (str) – The codon to mutate.
alt (str) – The alternative allele.
pos (int) – The position to mutate.

Return type:

str

Returns:

Mutated codon.

static is_synonymous(codon1: str, codon2: str)[source]#

Check if two codons are synonymous.

Parameters:

codon1 (str) – The first codon.
codon2 (str) – The second codon.

Return type:

bool

Returns:

True if the codons are synonymous, False otherwise.

annotate_site(v: Site)[source]#

Annotate a single site.

Parameters:: v (Site) – The variant to annotate.
Returns:: The annotated variant.

static count_target_sites(file: str, remove_overlaps: bool = False, contigs: List[str] = None)#

Count the number of target sites in a GFF file.

Parameters:

file (str) – The path to The GFF file path, possibly gzipped or a URL
remove_overlaps (bool) – Whether to remove overlapping target sites.
contigs (List[str]) – The contigs to count the target sites for.

Return type:

Dict[str, int]

Returns:

The number of target sites per chromosome/contig.

mismatches: List['cyvcf2.Variant']#: The variants that could not be annotated correctly.

n_skipped: int#: The variant that were skipped because they were not in coding regions.

errors: List['cyvcf2.Variant']#: The variants for which the codon could not be determined.

n_annotated: int#: The number of annotated sites.

AncestralAlleleAnnotation#

class AncestralAlleleAnnotation[source]#

Bases: Annotation, ABC

Base class for ancestral allele annotation.

Methods:

`__init__`	Create a new annotation instance.
`annotate_site`	Annotate a single site.
`count_target_sites`	Count the number of target sites in a GFF file.

__init__()#: Create a new annotation instance.

abstractmethod annotate_site(variant: Site)#

Annotate a single site.

Parameters:: variant (Site) – The variant to annotate.
Returns:: The annotated variant.

static count_target_sites(file: str, remove_overlaps: bool = False, contigs: List[str] = None)#

Count the number of target sites in a GFF file.

Parameters:

file (str) – The path to The GFF file path, possibly gzipped or a URL
remove_overlaps (bool) – Whether to remove overlapping target sites.
contigs (List[str]) – The contigs to count the target sites for.

Return type:

Dict[str, int]

Returns:

The number of target sites per chromosome/contig.

n_annotated: int#: The number of annotated sites.

MaximumParsimonyAncestralAnnotation#

class MaximumParsimonyAncestralAnnotation(samples: List[str] = None)[source]#

Bases: AncestralAlleleAnnotation

Annotation of ancestral alleles using maximum parsimony. To be used with AncestralBaseStratification and Annotator or Parser.

Note that maximum parsimony is not a reliable way to determine ancestral alleles, so it is recommended to use this annotation together with the ancestral misidentification parameter eps or to fold spectra altogether. Alternatively, you can use DeviantOutgroupFiltration to filter out sites where the major allele among outgroups does not coincide with the major allele among ingroups. This annotation has the advantage of requiring no outgroup data. If outgroup data is available consider using :class`MLEAncestralAlleleAnnotation` instead.

Methods:

`__init__`	Create a new ancestral allele annotation instance.
`annotate_site`	Annotate a single site.
`count_target_sites`	Count the number of target sites in a GFF file.

__init__(samples: List[str] = None)[source]#

Create a new ancestral allele annotation instance.

Parameters:: samples (List[str]) – The samples to consider when determining the ancestral allele. If None, all samples are considered.

samples: List[str] | None#: The samples to consider when determining the ancestral allele.

samples_mask: ndarray | None#

annotate_site(variant: Site)[source]#

Annotate a single site.

Parameters:: variant (Site) – The variant to annotate.
Returns:: The annotated variant.

static count_target_sites(file: str, remove_overlaps: bool = False, contigs: List[str] = None)#

Count the number of target sites in a GFF file.

Parameters:

file (str) – The path to The GFF file path, possibly gzipped or a URL
remove_overlaps (bool) – Whether to remove overlapping target sites.
contigs (List[str]) – The contigs to count the target sites for.

Return type:

Dict[str, int]

Returns:

The number of target sites per chromosome/contig.

n_annotated: int#: The number of annotated sites.

SubstitutionModel#

class SubstitutionModel(bounds: Dict[str, Tuple[float, float]] = {}, pool_branch_rates: bool = False, fixed_params: Dict[str, float] = {})[source]#

Bases: ABC

Base class for substitution models.

Methods:

`__init__`	Create a new substitution model instance.
`cache`	Cache the probabilities for the given parameters.
`get_x0`	Get the initial values for the parameters.
`get_bounds`	Get the bounds for the parameters.
`validate_bounds`	Make sure the lower bounds are positive and the upper bounds are larger than the lower bounds.

__init__(bounds: Dict[str, Tuple[float, float]] = {}, pool_branch_rates: bool = False, fixed_params: Dict[str, float] = {})[source]#

Create a new substitution model instance.

Parameters:

bounds (Dict[str, Tuple[float, float]]) – The bounds for the parameters.
pool_branch_rates (bool) – Whether to pool the branch rates. By default, each branch has its own rate which is optimized using MLE. If True, the branch rates are pooled and a single rate is optimized. This is useful if the number of sites used is small.
fixed_params (Dict[str, float]) – The fixed parameters. Parameters that are not fixed are optimized using MLE.

pool_branch_rates: bool#: Whether to pool the branch rates.

fixed_params: Dict[str, float]#: The fixed parameters.

bounds: Dict[str, Tuple[float, float]]#: Parameter bounds.

cache(params: Dict[str, float], n_branches: int)[source]#

Cache the probabilities for the given parameters.

Parameters:

params (Dict[str, float]) – The parameters.
n_branches (int) – The number of branches.

static get_x0(bounds: Dict[str, Tuple[float, float]], rng: Generator)[source]#

Get the initial values for the parameters.

Parameters:

bounds (Dict[str, Tuple[float, float]]) – The bounds for the parameters.
rng (Generator) – The random number generator.

Return type:

Dict[str, float]

Returns:

The initial values.

get_bounds(n_outgroups: int)[source]#

Get the bounds for the parameters.

Parameters:: n_outgroups (int) – The number of outgroups.
Return type:: Dict[str, Tuple[float, float]]
Returns:: The bounds.

static validate_bounds(bounds: Dict[str, Tuple[float, float]])[source]#

Make sure the lower bounds are positive and the upper bounds are larger than the lower bounds.

Parameters:: bounds (Dict[str, Tuple[float, float]]) – The bounds to validate
Raises:: ValueError – If the bounds are invalid

JCSubstitutionModel#

class JCSubstitutionModel(bounds: Dict[str, Tuple[float, float]] = {'K': (1e-05, 10)}, pool_branch_rates: bool = False, fixed_params: Dict[str, float] = {})[source]#

Bases: SubstitutionModel

Jukes-Cantor substitution model.

Methods:

`__init__`	Create a new substitution model instance.
`get_bound`	Get the bounds for a parameter.
`get_bounds`	Get the bounds for the parameters.
`cache`	Cache the probabilities for the given parameters.
`get_x0`	Get the initial values for the parameters.
`validate_bounds`	Make sure the lower bounds are positive and the upper bounds are larger than the lower bounds.

__init__(bounds: Dict[str, Tuple[float, float]] = {'K': (1e-05, 10)}, pool_branch_rates: bool = False, fixed_params: Dict[str, float] = {})[source]#

Create a new substitution model instance.

Parameters:

bounds (Dict[str, Tuple[float, float]]) – The bounds for the parameters. K is the branch rate.
pool_branch_rates (bool) – Whether to pool the branch rates. By default, each branch has its own rate which is optimized using MLE. If True, the branch rates are pooled and a single rate is optimized. This is useful if the number of sites used is small. If False, each branch has its own rate denoted by “K{i}”, where i is the branch index. If True, the branch rate is denoted by “K”.
fixed_params (Dict[str, float]) – The fixed parameters. Parameters that are not fixed are optimized using MLE.

get_bound(param: str)[source]#

Get the bounds for a parameter.

Parameters:: param (str) – The parameter.
Return type:: Tuple[float, float]
Returns:: The lower and upper bounds.

get_bounds(n_outgroups: int)[source]#

Get the bounds for the parameters.

Parameters:: n_outgroups (int) – The number of outgroups.
Return type:: Dict[str, Tuple[float, float]]
Returns:: The lower and upper bounds.

cache(params: Dict[str, float], n_branches: int)#

Cache the probabilities for the given parameters.

Parameters:

params (Dict[str, float]) – The parameters.
n_branches (int) – The number of branches.

static get_x0(bounds: Dict[str, Tuple[float, float]], rng: Generator)#

Get the initial values for the parameters.

Parameters:

bounds (Dict[str, Tuple[float, float]]) – The bounds for the parameters.
rng (Generator) – The random number generator.

Return type:

Dict[str, float]

Returns:

The initial values.

static validate_bounds(bounds: Dict[str, Tuple[float, float]])#

Make sure the lower bounds are positive and the upper bounds are larger than the lower bounds.

Parameters:: bounds (Dict[str, Tuple[float, float]]) – The bounds to validate
Raises:: ValueError – If the bounds are invalid

pool_branch_rates: bool#: Whether to pool the branch rates.

fixed_params: Dict[str, float]#: The fixed parameters.

bounds: Dict[str, Tuple[float, float]]#: Parameter bounds.

K2SubstitutionModel#

class K2SubstitutionModel(bounds: Dict[str, Tuple[float, float]] = {'K': (1e-05, 10), 'k': (0.1, 10)}, pool_branch_rates: bool = False, fixed_params: Dict[str, float] = {}, fix_transition_transversion_ratio: bool = False)[source]#

Bases: JCSubstitutionModel

Kimura 2-parameter substitution model.

Methods:

`__init__`	Create a new substitution model instance.
`get_bounds`	Get the bounds for the parameters.
`cache`	Cache the probabilities for the given parameters.
`get_bound`	Get the bounds for a parameter.
`get_x0`	Get the initial values for the parameters.
`validate_bounds`	Make sure the lower bounds are positive and the upper bounds are larger than the lower bounds.

__init__(bounds: Dict[str, Tuple[float, float]] = {'K': (1e-05, 10), 'k': (0.1, 10)}, pool_branch_rates: bool = False, fixed_params: Dict[str, float] = {}, fix_transition_transversion_ratio: bool = False)[source]#

Create a new substitution model instance.

Parameters:

bounds (Dict[str, Tuple[float, float]]) – The bounds for the parameters. K{i} are the branch rates. k is the transition/transversion ratio.
pool_branch_rates (bool) – Whether to pool the branch rates. By default, each branch has its own rate which is optimized using MLE. If True, the branch rates are pooled and a single rate is optimized. This is useful if the number of sites used is small.
fixed_params (Dict[str, float]) – The fixed parameters. Parameters that are not fixed are optimized using MLE.
fix_transition_transversion_ratio (bool) – Whether to fix the transition/transversion ratio to the ratio observed in the data.

fix_transition_transversion_ratio: bool#: Whether to fix the transition/transversion ratio to the ratio observed in the data.

get_bounds(n_outgroups: int)[source]#

Get the bounds for the parameters.

Parameters:: n_outgroups (int) – The number of outgroups.
Return type:: Dict[str, Tuple[float, float]]
Returns:: The lower and upper bounds.

cache(params: Dict[str, float], n_branches: int)#

Cache the probabilities for the given parameters.

Parameters:

params (Dict[str, float]) – The parameters.
n_branches (int) – The number of branches.

get_bound(param: str)#

Get the bounds for a parameter.

Parameters:: param (str) – The parameter.
Return type:: Tuple[float, float]
Returns:: The lower and upper bounds.

static get_x0(bounds: Dict[str, Tuple[float, float]], rng: Generator)#

Get the initial values for the parameters.

Parameters:

bounds (Dict[str, Tuple[float, float]]) – The bounds for the parameters.
rng (Generator) – The random number generator.

Return type:

Dict[str, float]

Returns:

The initial values.

static validate_bounds(bounds: Dict[str, Tuple[float, float]])#

Make sure the lower bounds are positive and the upper bounds are larger than the lower bounds.

Parameters:: bounds (Dict[str, Tuple[float, float]]) – The bounds to validate
Raises:: ValueError – If the bounds are invalid

pool_branch_rates: bool#: Whether to pool the branch rates.

fixed_params: Dict[str, float]#: The fixed parameters.

bounds: Dict[str, Tuple[float, float]]#: Parameter bounds.

SiteConfig#

class SiteConfig(n_major: int, major_base: int, minor_base: int, outgroup_bases: ndarray, multiplicity: float = 1.0, sites: ndarray = None, p_minor: float = nan, p_major: float = nan)[source]#

Bases: object

Ancestral allele site configuration for a single subsample.

Methods:

__init__

Create a new site configuration instance.

__init__(n_major: int, major_base: int, minor_base: int, outgroup_bases: ndarray, multiplicity: float = 1.0, sites: ndarray = None, p_minor: float = nan, p_major: float = nan)[source]#: Create a new site configuration instance.

n_major: int#: The number of major alleles.

major_base: int#: The major allele base index.

minor_base: int#: The minor base index.

outgroup_bases: ndarray#: The outgroup base indices.

multiplicity: float#: The multiplicity of the site.

sites: ndarray#: The site indices.

p_minor: float#: The probability of the minor allele.

p_major: float#: The probability of the major allele.

SiteInfo#

class SiteInfo(n_major: Dict[int, float], major_base: str, minor_base: str, outgroup_bases: List[str], p_minor: float = nan, p_major: float = nan, p_major_ancestral: float = nan, major_ancestral: str = '.', p_bases_first_node: Dict[str, float] = None, p_first_node_ancestral: float = nan, first_node_ancestral: str = '.', rate_params: Dict[str, float] = None)[source]#

Bases: object

Ancestral allele information on a single site.

Methods:

`__init__`	Initialize the site information.
`plot_tree`	Plot the tree for a site.

__init__(n_major: Dict[int, float], major_base: str, minor_base: str, outgroup_bases: List[str], p_minor: float = nan, p_major: float = nan, p_major_ancestral: float = nan, major_ancestral: str = '.', p_bases_first_node: Dict[str, float] = None, p_first_node_ancestral: float = nan, first_node_ancestral: str = '.', rate_params: Dict[str, float] = None)[source]#: Initialize the site information.

n_major: Dict[int, float]#: Dictionary mapping number of major alleles to its probability of observation.

major_base: str#: The major allele base.

minor_base: str#: The minor base index.

outgroup_bases: List[str]#: The outgroup base indices.

p_minor: float#: The probability of the minor allele being the ancestral allele (without prior).

p_major: float#: The probability of the major allele being the ancestral allele (without prior).

p_major_ancestral: float#: The probability of the major allele being the ancestral allele rather than the minor allele (possibly with prior if specified).

major_ancestral: str#: The predicted ancestral base based on comparing major and minor allele.

p_bases_first_node: Dict[str, float]#: The probability of each base being the ancestral base for the first internal node.

p_first_node_ancestral: float#: The probability that the mostly likely base for the first internal node is the ancestral base.

first_node_ancestral: str#: The ancestral base index for the first internal node.

rate_params: Dict[str, float]#: The branch rates.

plot_tree(ax: plt.Axes = None, show: bool = True)[source]#

Plot the tree for a site. Only Python visualization is supported.

Parameters:

self – The site information.
ax (plt.Axes) – Axes to plot on.
show (bool) – Whether to show the plot.

BaseType#

class BaseType(value)[source]#

Bases: Enum

The base type, either major or minor.

MINOR: int = 0#

MAJOR: int = 1#

PolarizationPrior#

class PolarizationPrior(allow_divergence: bool = False)[source]#

Bases: ABC

Base class for priors used with MaximumLikelihoodAncestralAnnotation. These priors incorporate information about the general probability of the major allele being ancestral across all sites with the same minor allele count. Prior thus take ingroup allele frequencies into account, when making predictions about the ancestral state of a site. This is useful because it enhances ancestral allele probability estimates, especially when outgroup information is unavailable for a particular site. Knowing the likelihood of the major allele being ancestral in general allows for more informed estimations.

Methods:

`__init__`	Create a new instance.
`plot`	Visualize the polarization probabilities using a scatter plot.

__init__(allow_divergence: bool = False)[source]#

Create a new instance.

Parameters:

allow_divergence (bool) –

Whether to allow divergence. If True, the probability of the minor allele being ancestral, which is not contained in the ingroup subsample but rather in all specified ingroups or among the outgroup, is taken to be the same as if it was present in the ingroup subsample with frequency 1. This is a hack, but allows us to consider alleles that are not present in the ingroup subsample.

Warning

Setting this to True greatly increases the probability of high-frequency derived alleles which introduces a strong bias in the distribution of frequency counts, e.g., the SFS. Only use this if you’re interested in handling divergence counts, i.e., sites where the ingroup is mono-allelic.

allow_divergence: bool#: Whether to allow divergence.

probabilities: ndarray | None#: The polarization probabilities.

plot(file: str = None, show: bool = True, title: str = 'polarization probabilities', scale: Literal['lin', 'log'] = 'lin', ax: plt.Axes = None, ylabel: str = 'p')[source]#

Visualize the polarization probabilities using a scatter plot.

Parameters:

scale (Literal['lin', 'log']) – y-scale of the plot.
title (str) – Plot title.
file (str) – File to save plot to.
show (bool) – Whether to show plot.
ax (plt.Axes) – Axes to plot on. Only for Python visualization backend.
ylabel (str) – y-axis label.

Return type:

plt.Axes

Returns:

Axes object

KingmanPolarizationPrior#

class KingmanPolarizationPrior(allow_divergence: bool = False)[source]#

Bases: PolarizationPrior

Prior based on the standard Kingman coalescent. To be used with MaximumLikelihoodAncestralAnnotation.

Methods:

`__init__`	Create a new instance.
`plot`	Visualize the polarization probabilities using a scatter plot.

__init__(allow_divergence: bool = False)#

Create a new instance.

Parameters:

allow_divergence (bool) –

Warning

plot(file: str = None, show: bool = True, title: str = 'polarization probabilities', scale: Literal['lin', 'log'] = 'lin', ax: plt.Axes = None, ylabel: str = 'p')#

Visualize the polarization probabilities using a scatter plot.

Parameters:

scale (Literal['lin', 'log']) – y-scale of the plot.
title (str) – Plot title.
file (str) – File to save plot to.
show (bool) – Whether to show plot.
ax (plt.Axes) – Axes to plot on. Only for Python visualization backend.
ylabel (str) – y-axis label.

Return type:

plt.Axes

Returns:

Axes object

allow_divergence: bool#: Whether to allow divergence.

probabilities: ndarray | None#: The polarization probabilities.

AdaptivePolarizationPrior#

class AdaptivePolarizationPrior(n_runs: int = 1, parallelize: bool = True, allow_divergence: bool = False, seed: int | None = 0)[source]#

Bases: PolarizationPrior

Adaptive prior. To be used with MaximumLikelihoodAncestralAnnotation. This is the same prior as used in the EST-SFS paper. This prior is adaptive in the sense that the most likely polarization probabilities given the site configurations are found. This is the most accurate prior, but requires a lot of sites in order to work properly. You can check that the polarization probabilities are smooth enough across frequency counts by calling plot(). If they are not smooth enough, you can increase the number of sites, decrease the number of ingroups, or use KingmanPolarizationPrior instead.

Note

In practice, this prior provides very similar results to KingmanPolarizationPrior in most cases.

Methods:

`__init__`	Create a new adaptive prior instance.
`plot`	Visualize the polarization probabilities using a scatter plot.

__init__(n_runs: int = 1, parallelize: bool = True, allow_divergence: bool = False, seed: int | None = 0)[source]#

Create a new adaptive prior instance.

Parameters:

n_runs (int) – The number of runs to perform when determining the polarization parameters. One run should be sufficient as only one parameter is optimized.
parallelize (bool) – Whether to parallelize the optimization.
allow_divergence (bool) – Whether to allow divergence. See PolarizationPrior for details.
seed (int | None) – The seed for the random number generator.

n_runs: int#: The number of runs to use for the adaptive prior.

parallelize: bool#: Whether to parallelize the optimization.

seed: int | None#: The seed for the random number generator.

rng: Generator#: The random number generator.

plot(file: str = None, show: bool = True, title: str = 'polarization probabilities', scale: Literal['lin', 'log'] = 'lin', ax: plt.Axes = None, ylabel: str = 'p')#

Visualize the polarization probabilities using a scatter plot.

Parameters:

scale (Literal['lin', 'log']) – y-scale of the plot.
title (str) – Plot title.
file (str) – File to save plot to.
show (bool) – Whether to show plot.
ax (plt.Axes) – Axes to plot on. Only for Python visualization backend.
ylabel (str) – y-axis label.

Return type:

plt.Axes

Returns:

Axes object

allow_divergence: bool#: Whether to allow divergence.

probabilities: ndarray | None#: The polarization probabilities.

MaximumLikelihoodAncestralAnnotation#

class MaximumLikelihoodAncestralAnnotation(outgroups: List[str], n_ingroups: int = 11, ingroups: List[str] | None = None, exclude: List[str] | None = None, n_runs: int = 10, model: SubstitutionModel = None, parallelize: bool = True, prior: PolarizationPrior | None = '', max_sites: int = 10000, seed: int | None = 0, confidence_threshold: float = 0, n_target_sites: int | None = None, n_samples_target_sites: int | None = 100000, adjust_target_sites: bool = True, subsample_mode: Literal['random', 'probabilistic'] = 'probabilistic')[source]#

Bases: _OutgroupAncestralAlleleAnnotation

Annotation of ancestral alleles following the probabilistic model of EST-SFS (https://doi.org/10.1534/genetics.118.301120). Note that only bi-allelic SNPs are supported by this model. By default, the info field AA (see Annotator.info_ancestral) is added to the VCF file, which holds the ancestral allele. To be used with Annotator or Parser. The info field AA_prob holds the probability that the predicted ancestral allele (AA tag) is correct, as opposed to the other allele. This probability can be also be used by Parser to polarize the SFS according to the ancestral allele probability. In addition to the ancestral allele, the info field AA_info is added, which contains additional information about the ancestral allele (see SiteInfo for an overview of the available information). This class can also be used independently, see the from_dataframe(), from_data() and from_est_sfs() methods.

Initially, the branch rates are determined using MLE. Similar to Parser, we can also specify the number of mutational target sites (see the n_target_sites argument) in case our input does not contain the full set of monomorphic sites. This is necessary to obtain realistic branch rate estimates. You can also choose a prior for the polarization probabilities (see PolarizationPrior). Eventually, for every site, the probability that the major allele is ancestral is calculated.

When annotating the variants of an input, we check the most likely ancestral allele against a naive ad-hoc ancestral allele annotation, and record the sites for which we have disagreement. You might want to sanity-check the mismatches to make sure the model has been properly specified (see mismatches).

Note

The polarization prior corresponds to the Kingman coalescent probability by default. Using an adaptive prior, as in the EST-SFS paper, is also possible, but this is only recommended if the number of sites used for the inference is large (see prior).
The model can only handle sites that have at most 2 alleles across the in- and outgroups, so sites with more than 2 alleles are ignored. Only variants that are at most bi-allelic in the provided in- and outgroups are annotated.
The model determines the probability of the major allele being ancestral opposed to the minor allele. This can be problematic if the actual ancestral allele is not contained in the ingroup (possibly due to subsampling). To avoid this issue, we also keep track of potential minor alleles at frequency 0. If we were to ignore this, it would be impossible to infer divergence, i.e. fixed derived allele that are no longer observed in the ingroups (see PolarizationPrior.allow_divergence). That said, if only the SFS is of interest, allow_divergence should not be set to True.
The model assumes a single coalescent topology for all sites, in which all outgroups coalesce first with the ingroup and not with each other. It is important to specify the outgroups in order of increasing divergence and not to select outgroups that are not much more closely related to each other than to the ingroup (as this would give rise to a different coalescent topology than the one assumed). You can call get_outgroup_divergence() after the inference to check the estimated branch rates for each outgroup. The assumption of a single fixed topology should be good enough provided that in- and outgroups are sufficiently diverged.

Example usage:

import sfsutils as su

ann = su.Annotator(
    source="path/to/variants.with_outgroups.vcf.gz",
    annotations=[su.MaximumLikelihoodAncestralAnnotation(
        outgroups=["ERR2103730"],
        n_ingroups=15
    )],
    output="genome.polarized.vcf.gz"
)

ann.annotate()

Methods:

`__init__`	Create a new ancestral allele annotation instance.
`from_est_sfs`	Create instance from EST-SFS input file.
`to_est_sfs`	Write the object's state to an EST-SFS formatted file.
`to_file`	Save object to file.
`from_file`	Load object from file.
`to_json`	Serialize object.
`from_json`	Load object from file.
`from_data`	Create an instance by passing the data directly.
`from_dataframe`	Create an instance from a dataframe.
`infer`	Infer the ancestral allele probabilities for the data provided.
`set_mle_params`	Set the MLE parameters and update the cache and site configurations.
`is_monotonic`	Whether the outgroups are monotonically increasing in divergence.
`get_p_tree`	Get the probability of a tree.
`get_p_config`	Get the probability for a site configuration.
`get_p_configs`	Get the probabilities for each site configuration.
`evaluate_likelihood`	Evaluate the likelihood function for the rate parameters.
`get_inferred_site_info`	Get the site information for the sites included in the parsing process.
`get_site_info`	Get information on the specified sites using the inferred parameters.
`annotate_site`	Annotate a single site.
`plot_likelihoods`	Visualize the likelihoods of the rate optimization runs using a scatter plot.
`get_folded_spectra`	Get the folded spectra for the parsed sites (used to estimate the parameters).
`get_outgroup_divergence`	Get the inferred branch rates between the ingroup and outgroups by combining the inferred branch rates.
`get_observed_transition_transversion_ratio`	Get the observed transition/transversion ratio.
`count_target_sites`	Count the number of target sites in a GFF file.
`get_base_index`	Get base index/indices from base string(s).
`get_base_string`	Get base string(s) from base index/indices.

max_sampling_attempts_factor: int = 100#: Number of random draws per requested monomorphic sample, bounding the sampling loop so an all-masked (non-ACGT) window cannot hang. Tune up if a reference is heavily masked.

__init__(outgroups: List[str], n_ingroups: int = 11, ingroups: List[str] | None = None, exclude: List[str] | None = None, n_runs: int = 10, model: SubstitutionModel = None, parallelize: bool = True, prior: PolarizationPrior | None = '', max_sites: int = 10000, seed: int | None = 0, confidence_threshold: float = 0, n_target_sites: int | None = None, n_samples_target_sites: int | None = 100000, adjust_target_sites: bool = True, subsample_mode: Literal['random', 'probabilistic'] = 'probabilistic')[source]#

Create a new ancestral allele annotation instance.

Parameters:

outgroups (List[str]) – The outgroup samples to consider when determining the ancestral allele in the order of increasing divergence. A list of sample names as they appear in the input. The order of the outgroups is important as it determines the order of the branches in the tree, whose rates are optimized, and whose topology is predetermined. The first outgroup is the closest outgroup to the ingroups, and the last outgroup is the most distant outgroup. More outgroups lead to a more accurate inference of the ancestral allele, but also increase the computational cost. Using more than 1 outgroup is recommended, but more than 3 is likely not necessary. Sites where these outgroups are not present are not included when optimizing the rate parameters. Due to assumptions on the tree topology connecting the in- and outgroups, it is important that the outgroups are not much more closely related to each other than to the ingroups. Ideally, the optimized branch rates are show markedly different values, and in any case, they should be monotonically increasing with the outgroups (see get_outgroup_divergence()).
n_ingroups (int) – The minimum number of ingroups that must be present at a site for it to be considered for ancestral allele inference. The ingroup subsampling is necessary since our model requires an equal number of ingroups for all sites. Note that a larger number of ingroups does not necessarily improve the accuracy of the ancestral allele inference (see prior). A larger number of ingroups can lead to a large variance in the polarization probabilities, across different frequency counts. n_ingroups should thus only be large if the number of sites used for the inference is also large. A sensible value for a reasonably large number of sites (a few thousand) is 10 or perhaps 20 for a larger numbers of sites. Very small values can lead to the ingroup subsamples not being representative of the actual allele frequencies at a site, especially when not using probabilistic subsampling (see subsample_mode). This value also influences the number of frequency bins used for the polarization probabilities, and should thus not be too small. Note that if ingroups is an even number, the major allele is chosen arbitrarily if the number of major alleles is equal to the number of minor alleles. To avoid this, you can use an odd number of ingroups.
ingroups (Optional[List[str]]) – The ingroup samples to consider when determining the ancestral allele. If None, all (non-outgroup) samples are considered. A list of sample names as they appear in the input. Has to be at least as large as n_ingroups.
exclude (Optional[List[str]]) – Samples to exclude from the ingroup. A list of sample names as they appear in the input.
n_runs (int) – The number of optimization runs to perform when determining the branch rates. You can check that the likelihoods of the different runs are similar by calling plot_likelihoods().
model (SubstitutionModel) – The substitution model to use. By default, K2SubstitutionModel is used.
parallelize (bool) – Whether to parallelize the computation across multiple cores.
prior (PolarizationPrior | None) – The prior to use for the polarization probabilities. See KingmanPolarizationPrior and AdaptivePolarizationPrior for more information. By default, KingmanPolarizationPrior is used. Use None for no prior.
max_sites (int) – The maximum number of sites to consider. This is useful if the number of sites is very large. Choosing a reasonably large subset of sites (on the order of a few thousand bi-allelic sites) can speed up the computation considerably as parsing can be slow. This subset is then used to calibrate the rate parameters, and possibly the polarization priors.
seed (int | None) – The seed for the random number generator. If None, a random seed is chosen. By default, the seed is set to 0.
confidence_threshold (float) –
The confidence threshold for the ancestral allele annotation. Only if the probability of the major allele being ancestral as opposed to the minor allele is not within ((1 - confidence_threshold) / 2, 1 - (1 - confidence_threshold) / 2), the ancestral allele is annotated. This is useful to avoid annotating sites where the ancestral allele state is not clear. Use values close to 0 to annotate as many sites as possible, and values close to 1 to annotate only sites where the ancestral allele state is very clear.

Warning

This threshold introduces a bias by excluding more sites with high-frequency derived alleles and should thus be kept at 0 if the distribution of frequency counts is important, e.g., if the SFS is to be determined.
n_target_sites (int | None) – The total number of target sites if this class is used in conjunction with Parser or Annotator. This is useful if the provided set of sites only consists of bi-allelic sites. Specify here the total number of sites underlying the given dataset, i.e., both mono- and bi-allelic sites. Ignoring mono-allelic sites will lead to overestimation of the rate parameters. For this to work, a FASTA file must be provided from which the mono-allelic sites can be sampled. Sampling takes place between the variants of the last and first site on every contig considered in the input. Use None to disable this feature. Note that the number of target sites is automatically carried over if not specified and this class is used together with Parser. In order to use this feature, you also need to specify a FASTA file to Parser or Annotator. Also note that by default we extrapolate the number of mono-allelic sites to be sampled from the FASTA file based on the ratio of sites with called outgroup bases parsed from the input (adjust_target_sites).
n_samples_target_sites (int | None) – The number of sites to sample from the FASTA file when determining the number of target sites (n_target_sites). From this the total number of target sites is extrapolated.
adjust_target_sites (bool) – Whether to adjust the number of target sites based on the parsed sites relative to the total number of sites in the input. Defaults to True.
subsample_mode (Literal['random', 'probabilistic']) – The subsampling mode. For random, we draw once without replacement from the set of all available ingroup genotypes per site. For probabilistic, we integrate over the hypergeometric distribution when parsing and computing the ancestral probabilities. Probabilistic subsampling requires a bit more time, but produces much more stable results, while requiring far fewer sites, so it is highly recommended.

parallelize: bool#: Whether to parallelize the computation.

max_sites: int#: Maximum number of sites to consider

confidence_threshold: float#: The confidence threshold for the ancestral allele annotation.

prior: PolarizationPrior | None#: The prior to use for the polarization probabilities.

n_runs: int#: Number of random ML starts when determining the rate parameters

model: SubstitutionModel#: The substitution model.

configs: pd.DataFrame | None#: The data frame holding all site configurations.

p_bins: Dict[str, np.ndarray | None]#: The probability of all sites per frequency bin.

n_sites: int | None#: The total number of valid sites parsed (including sites not considered for ancestral allele inference).

param_names: List[str]#: The parameter names in the order they are passed to the optimizer.

likelihoods: np.ndarray | None#: The log likelihoods for the different runs when optimizing the rate parameters.

likelihood: float | None#: The best log likelihood when optimizing the rate parameters.

result: 'scipy.optimize.OptimizeResult' | None#: Optimization result of the best run.

params_mle: Dict[str, float] | None#: The MLE parameters.

params_mle_runs: pd.DataFrame | None#: The MLE parameters for all runs.

mismatches: List[SiteInfo]#

n_target_sites: int | None#: The total number of target sites.

n_samples_target_sites: int#: The number of sites to sample from the FASTA file when determining the number of target sites.

adjust_target_sites: bool#: Whether to adjust the number of target sites based on the number of sites parsed from the input.

classmethod from_est_sfs(file: str, prior: PolarizationPrior | None = '', n_runs: int = 10, model: SubstitutionModel = None, parallelize: bool = True, seed: int = 0, chunk_size: int = 100000)[source]#

Create instance from EST-SFS input file.

Parameters:

file (str) – File containing EST-SFS-formatted input data.
prior (PolarizationPrior | None) – The prior to use for the polarization probabilities (see __init__()).
n_runs (int) – Number of runs for rate estimation (see __init__()).
model (SubstitutionModel) – The substitution model (see __init__()).
parallelize (bool) – Whether to parallelize the runs (see __init__()).
seed (int) – The seed to use for the random number generator.
chunk_size (int) – The chunk size for reading the file.

Return type:

MaximumLikelihoodAncestralAnnotation

Returns:

The instance.

to_est_sfs(file: str)[source]#

Write the object’s state to an EST-SFS formatted file.

Parameters:: file (str) – The output file name.

to_file(file: str)[source]#

Save object to file.

Note

References to the handler and the reader are discarded.

Parameters:: file (str) – File path.

classmethod from_file(file: str)[source]#

Load object from file.

Note

The handler and the reader are not restored, so serialization is mostly useful for analyzing the results in detail, not for further processing.

Parameters:: file (str) – File path.
Return type:: MaximumLikelihoodAncestralAnnotation
Returns:: The object.

to_json()[source]#

Serialize object.

Note

References to the handler and the reader are discarded.

Return type:: str
Returns:: JSON string

classmethod from_json(json: str)[source]#

Load object from file.

Note

The handler and the reader are not restored, so serialization is mostly useful for analyzing the results in detail, not for further processing.

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

classmethod from_data(n_major: Iterable[int], major_base: Iterable[str | int], minor_base: Iterable[str | int], outgroup_bases: Iterable[Iterable[str | int]], n_ingroups: int, n_runs: int = 10, model: SubstitutionModel = None, parallelize: bool = True, prior: PolarizationPrior | None = '', seed: int = 0, pass_indices: bool = False, confidence_threshold: float = 0)[source]#

Create an instance by passing the data directly.

Parameters:

n_major (Iterable[int]) – The number of major alleles per site. Note that this number has to be lower than n_ingroups, as we consider the number of major alleles of subsamples of size n_ingroups.
major_base (Iterable[str | int]) – The major allele per site. A string representation of the base or the base index according to ['A', 'C', 'G', 'T'] if pass_indices is True. Use None if the base is not defined when pass_indices is False and -1 when pass_indices is True.
minor_base (Iterable[str | int]) – The minor allele per site. A string representation of the base or the base index according to ['A', 'C', 'G', 'T'] if pass_indices is True. Use None if the base is not defined when pass_indices is False and -1 when pass_indices is True.
outgroup_bases (Iterable[Iterable[str | int]]) – The outgroup alleles per site. A string representation of the base or the base index if pass_indices is True. This should be a list of lists, where the outer list corresponds to the sites and the inner list to the outgroups per site. All sites are required to have the same number of outgroups. Use None if the base is not defined when pass_indices is False and -1 when pass_indices is True.
n_ingroups (int) – The number of ingroup samples (see __init__()).
n_runs (int) – The number of runs for rate estimation (see __init__()).
model (SubstitutionModel) – The substitution model (see __init__()).
parallelize (bool) – Whether to parallelize the runs.
prior (PolarizationPrior | None) – The prior to use for the polarization probabilities (see __init__()).
seed (int) – The seed for the random number generator.
pass_indices (bool) – Whether to pass the base indices instead of the bases.
confidence_threshold (float) – The confidence threshold for the ancestral allele annotation (see __init__()).

Return type:

MaximumLikelihoodAncestralAnnotation

Returns:

The instance.

classmethod from_dataframe(data: DataFrame, n_ingroups: int, n_runs: int = 10, model: SubstitutionModel = None, parallelize: bool = True, prior: PolarizationPrior | None = '', seed: int = 0, grouped: bool = False, confidence_threshold: float = 0)[source]#

Create an instance from a dataframe.

Parameters:

data (DataFrame) – Dataframe with the columns: major_base, minor_base, outgroup_bases, n_major of type int, int, list and int, respectively. The outgroup bases should have the same length for every site.
n_ingroups (int) – The number of ingroups (see __init__()).
n_runs (int) – Number of runs for rate estimation (see __init__()).
model (SubstitutionModel) – The substitution model (see __init__()).
parallelize (bool) – Whether to parallelize computations.
prior (PolarizationPrior | None) – The prior to use for the polarization probabilities (see __init__()).
seed (int) – The seed for the random number generator. If None, a random seed is chosen.
grouped (bool) – Whether the dataframe is already grouped by all columns (used for internal purposes).
confidence_threshold (float) – The confidence threshold for the ancestral allele annotation (see __init__()).

Return type:

MaximumLikelihoodAncestralAnnotation

Returns:

The instance.

infer()[source]#: Infer the ancestral allele probabilities for the data provided. This method is only supposed to be called manually if the data is provided directly, e.g. using from_data(), from_dataframe() or from_est_sfs(). If the data is provided using an input source, this method is called automatically.

set_mle_params(params: Dict[str, float])[source]#

Set the MLE parameters and update the cache and site configurations. Use this method if you want to use different parameters for the annotation.

Parameters:: params (Dict[str, float]) – The new parameters.

is_monotonic()[source]#

Whether the outgroups are monotonically increasing in divergence.

Return type:: bool
Returns:: Whether the outgroups are monotonically increasing in divergence.

property p_polarization: ndarray | None#: Get the polarization probabilities or None if prior is no.

static get_p_tree(base: int, n_outgroups: int, internal_nodes: List[int] | ndarray, outgroup_bases: List[int] | ndarray, params: Dict[str, float], model: SubstitutionModel)[source]#

Get the probability of a tree.

Parameters:

base (int) – An observed ingroup base index.
n_outgroups (int) – The number of outgroups.
internal_nodes (Union[List[int], ndarray]) – The internal nodes of the tree. We have n_outgroups - 1 internal nodes.
outgroup_bases (Union[List[int], ndarray]) – The observed base indices for the outgroups.
params (Dict[str, float]) – The parameters of the model.
model (SubstitutionModel) – The model to use. Either ‘K2’ or ‘JC’.

Return type:

float

classmethod get_p_config(config: ~sfsutils.annotation.SiteConfig, base_type: ~sfsutils.annotation.BaseType, params: ~typing.Dict[str, float], model: ~sfsutils.annotation.SubstitutionModel = <sfsutils.annotation.K2SubstitutionModel object>, internal: ~numpy.ndarray | None = None)[source]#

Get the probability for a site configuration.

Parameters:

config (SiteConfig) – The site configuration.
base_type (BaseType) – The base type.
params (Dict[str, float]) – The parameters for the substitution model.
model (SubstitutionModel) – The substitution model to use.
internal (ndarray | None) – Base indices of internal nodes of the tree if fixed. If None, the internal nodes are considered as free parameters. -1 also indicates a free parameter. The number of internal nodes is the number of outgroups minus one.

Return type:

float

Returns:

The probability for a site.

classmethod get_p_configs(configs: DataFrame, model: SubstitutionModel, base_type: BaseType, params: Dict[str, float])[source]#

Get the probabilities for each site configuration.

Parameters:

configs (DataFrame) – The site configurations.
model (SubstitutionModel) – The substitution model.
base_type (BaseType) – The base type.
params (Dict[str, float]) – A dictionary of the rate parameters.

Return type:

ndarray

Returns:

The probability for each site.

evaluate_likelihood(params: Dict[str, float])[source]#

Evaluate the likelihood function for the rate parameters.

Parameters:: params (Dict[str, float]) – A dictionary of parameters.
Return type:: float
Returns:: The log likelihood.

get_inferred_site_info()[source]#

Get the site information for the sites included in the parsing process. The sites are in the same order as parsed. You can use get_site_info() to get the site information for a specific site.

Return type:: Generator[SiteInfo, None, None]
Returns:: A generator yielding a dictionary with the site information (see get_site_info()).
Raises:: RuntimeError – If the subsample mode is probabilistic.

get_site_info(n_major: int, major_base: int | str, minor_base: int | str, outgroup_bases: List[int | str] | ndarray, pass_indices: bool = False)[source]#

Get information on the specified sites using the inferred parameters.

Parameters:

n_major (int) – The number of copies of the major allele.
major_base (int | str) – The major bases indices or strings.
minor_base (int | str) – The minor bases indices or strings.
outgroup_bases (Union[List[int | str], ndarray]) – The outgroup base indices or strings.
pass_indices (bool) – Whether to pass the indices as strings or convert them to integers.

Return type:

SiteInfo

Returns:

The site information.

annotate_site(variant: Site)[source]#

Annotate a single site.

Parameters:: variant (Site) – The variant to annotate.
Returns:: The annotated variant.

plot_likelihoods(file: str = None, show: bool = True, title: str = 'rate likelihoods', scale: Literal['lin', 'log'] = 'lin', ax: plt.Axes = None, ylabel: str = 'lnl')[source]#

Visualize the likelihoods of the rate optimization runs using a scatter plot.

Parameters:

scale (Literal['lin', 'log']) – y-scale of the plot.
title (str) – Plot title.
file (str) – File to save plot to.
show (bool) – Whether to show plot.
ax (plt.Axes) – Axes to plot on. Only for Python visualization backend.
ylabel (str) – Label for y-axis.

Return type:

plt.Axes

Returns:

Axes object

get_folded_spectra(groups: List[Literal['major_base', 'minor_base', 'outgroup_bases']] = ['major_base'])[source]#

Get the folded spectra for the parsed sites (used to estimate the parameters).

Parameters:: groups (List[Literal['major_base', 'minor_base', 'outgroup_bases']]) – The groups to group the spectra by.
Return type:: Spectra
Returns:: Spectra object

get_outgroup_divergence()[source]#

Get the inferred branch rates between the ingroup and outgroups by combining the inferred branch rates.

Return type:: ndarray
Returns:: One rate for each outgroup.

get_observed_transition_transversion_ratio()[source]#

Get the observed transition/transversion ratio. Note that this may differ from the estimated ratio for K2SubstitutionModel.

Return type:: float
Returns:: The observed transition/transversion ratio.

static count_target_sites(file: str, remove_overlaps: bool = False, contigs: List[str] = None)#

Count the number of target sites in a GFF file.

Parameters:

file (str) – The path to The GFF file path, possibly gzipped or a URL
remove_overlaps (bool) – Whether to remove overlapping target sites.
contigs (List[str]) – The contigs to count the target sites for.

Return type:

Dict[str, int]

Returns:

The number of target sites per chromosome/contig.

classmethod get_base_index(base_string: str | ndarray)#

Get base index/indices from base string(s).

Parameters:: base_string (str | ndarray) – The base string(s).
Return type:: int | ndarray
Returns:: Base index/indices.

static get_base_string(indices: int | ndarray)#

Get base string(s) from base index/indices.

Parameters:: indices (int | ndarray) – The base index/indices.
Return type:: str | ndarray
Returns:: Base string(s).

ingroups: List[str] | None#: The ingroup samples to consider when determining the ancestral allele.

exclude: List[str]#: The samples excluded from the ingroup.

outgroups: List[str]#: The outgroup samples to consider when determining the ancestral allele.

n_ingroups: int#: The number of ingroups.

n_outgroups: int#: The number of outgroups.

seed: int | None#: The seed for the random number generator.

subsample_mode: Literal['random', 'probabilistic']#: The subsampling mode.

rng: np.random.Generator#: The random number generator.

n_annotated: int#: The number of annotated sites.

AdHocAncestralAnnotation#

class AdHocAncestralAnnotation(outgroups: List[str], n_ingroups: int, ingroups: List[str] | None = None, exclude: List[str] = [], seed: int | None = 0, subsample_mode: Literal['random', 'probabilistic'] = 'random')[source]#

Bases: _OutgroupAncestralAlleleAnnotation

Ad-hoc ancestral allele annotation using simple rules. Used for testing and sanity checking.

Methods:

`annotate_site`	Annotate a single site.
`__init__`	Create a new ancestral allele annotation instance.
`count_target_sites`	Count the number of target sites in a GFF file.
`get_base_index`	Get base index/indices from base string(s).
`get_base_string`	Get base string(s) from base index/indices.

annotate_site(variant: Site)[source]#

Annotate a single site. Mono-allelic sites are assigned the major allele as ancestral. Sites with more than two alleles are ignored.

Parameters:: variant (Site) – The variant to annotate.
Returns:: The annotated variant.

__init__(outgroups: List[str], n_ingroups: int, ingroups: List[str] | None = None, exclude: List[str] = [], seed: int | None = 0, subsample_mode: Literal['random', 'probabilistic'] = 'random')#

Create a new ancestral allele annotation instance.

Parameters:

outgroups (List[str]) – The outgroup samples to consider when determining the ancestral allele. A list of sample names as they appear in the input.
n_ingroups (int) – The minimum number of ingroups that must be present at a site for it to be considered for ancestral allele inference.
ingroups (Optional[List[str]]) – The ingroup samples to consider when determining the ancestral allele. A list of sample names as they appear in the input. If None, all samples except the outgroups are considered.
exclude (List[str]) – Samples to exclude from the ingroup. A list of sample names as they appear in the input.
seed (int | None) – The seed for the random number generator.
subsample_mode (Literal['random', 'probabilistic']) – The subsampling mode. Either ‘random’ or ‘probabilistic’.

static count_target_sites(file: str, remove_overlaps: bool = False, contigs: List[str] = None)#

Count the number of target sites in a GFF file.

Parameters:

file (str) – The path to The GFF file path, possibly gzipped or a URL
remove_overlaps (bool) – Whether to remove overlapping target sites.
contigs (List[str]) – The contigs to count the target sites for.

Return type:

Dict[str, int]

Returns:

The number of target sites per chromosome/contig.

classmethod get_base_index(base_string: str | ndarray)#

Get base index/indices from base string(s).

Parameters:: base_string (str | ndarray) – The base string(s).
Return type:: int | ndarray
Returns:: Base index/indices.

static get_base_string(indices: int | ndarray)#

Get base string(s) from base index/indices.

Parameters:: indices (int | ndarray) – The base index/indices.
Return type:: str | ndarray
Returns:: Base string(s).

ingroups: List[str] | None#: The ingroup samples to consider when determining the ancestral allele.

exclude: List[str]#: The samples excluded from the ingroup.

outgroups: List[str]#: The outgroup samples to consider when determining the ancestral allele.

n_ingroups: int#: The number of ingroups.

n_outgroups: int#: The number of outgroups.

seed: int | None#: The seed for the random number generator.

subsample_mode: Literal['random', 'probabilistic']#: The subsampling mode.

rng: Generator#: The random number generator.

n_annotated: int#: The number of annotated sites.

Site annotation

Contents

Site annotation#

Annotator#

Annotation#

DegeneracyAnnotation#

SynonymyAnnotation#

AncestralAlleleAnnotation#

MaximumParsimonyAncestralAnnotation#

SubstitutionModel#

JCSubstitutionModel#

K2SubstitutionModel#

SiteConfig#

SiteInfo#

BaseType#

PolarizationPrior#

KingmanPolarizationPrior#

AdaptivePolarizationPrior#

MaximumLikelihoodAncestralAnnotation#

AdHocAncestralAnnotation#