Joint SFS

Given several populations, the Parser counts derived alleles jointly rather than pooling them, yielding the joint SFS: an array whose entry \((i, j, \dots)\) is the number of sites with \(i\) derived alleles in the first population, \(j\) in the second, and so on. Populations are supplied as a mapping from name to sample names via pops, and n sets the sample size each population is projected to. Parsing then returns a JointSpectra, holding one JointSFS per type.

The example VCF for Betula spp. carries samples from several localities, encoded in the sample names. We take two of them as populations.

import sfsutils as su
from cyvcf2 import VCF

vcf = "resources/genome/betula/biallelic.polarized.subset.10000.vcf.gz"
samples = VCF(vcf).samples

pops = {loc: [s for s in samples if s.startswith(loc)] for loc in ["ASP", "PIT"]}

sfs = su.Parser(source=vcf, pops=pops, n=10).parse()

INFO:Parser: Using stratification: [all].
INFO:Parser: Loading VCF file
Parser>Counting sites: 10000it [00:00, 398327.03it/s]
Parser>Processing sites: 100%|██████████| 10000/10000 [00:01<00:00, 6036.74it/s]
INFO:PolyAllelicFiltration: Filtered out 0 sites.
INFO:Parser: Skipped 345 sites without ancestral allele information.
INFO:Parser: Included 9496 out of 10000 sites in total from the input.

For two populations the joint SFS is a matrix, plotted as a heatmap of derived-allele counts.

sfs["all"].plot();

When the ancestral state is unknown, the joint SFS is folded so that each allele count is combined with its complement.

sfs["all"].fold().plot();

Marginalizing over one population recovers the single-population SFS of the other, which lets the joint spectrum be reconciled with a one-dimensional parse.

su.Spectrum(sfs["all"].marginalize([0]).data).plot();  # the SFS of population ASP