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.

library(sfsutils)

su <- load_sfsutils()

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

cyvcf2 <- reticulate::import("cyvcf2")

vcf <- "resources/genome/betula/biallelic.polarized.subset.10000.vcf.gz"
samples <- cyvcf2$VCF(vcf)$samples

pops <- list()
for (loc in c("ASP", "PIT")) {
  pops[[loc]] <- samples[startsWith(samples, loc)]
}

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, 388923.26it/s]
Parser>Processing sites: 100%|██████████| 10000/10000 [00:01<00:00, 5966.52it/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.

p <- sfs["all"]$plot()

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

p <- 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.

p <- su$Spectrum(sfs["all"]$marginalize(list(0))$data)$plot()  # the SFS of population ASP