Parsing the SFS

The Parser reads variants and counts derived-allele frequencies into a site-frequency spectrum. For a sample of n haplotypes it returns counts in bins 0, 1, ..., n, projecting down when more haplotypes are present. Polarisation relies on the ancestral allele: by default the AA info field is consulted, and sites where it is undefined are skipped (see skip_non_polarized). When no ancestral state is available the spectrum can instead be folded, at the cost of conflating minor and major allele counts.

A single spectrum is rarely enough on its own. More often we want to separate sites into classes that behave differently and read off one spectrum per class. This is expressed by passing Stratification objects to the Parser. The example below contrasts 0-fold and 4-fold degenerate sites using a VCF for Betula spp.

library(sfsutils)

su <- load_sfsutils()

# instantiate parser
p <- su$Parser(
  n = 8,
  source = "resources/genome/betula/biallelic.polarized.subset.10000.vcf.gz",
  fasta = "resources/genome/betula/genome.subset.20.fasta",
  gff = "resources/genome/betula/genome.gff.gz",
  annotations = list(
    su$DegeneracyAnnotation()
  ),
  stratifications = list(su$DegeneracyStratification())
)

# parse SFS
spectra <- p$parse()

INFO:Parser: Using stratification: [neutral, selected].
INFO:Parser: Loading VCF file
INFO:Parser: Loading GFF file
INFO:sfsutils: Unzipping resources/genome/betula/genome.gff.gz to /var/folders/n4/m5q2jgw91zv0tp1c4j9bh48m0000gn/T/tmpn1bg9eo8.gff
INFO:Parser: Loading FASTA file
Parser>Counting sites: 10000it [00:00, 394958.76it/s]
Parser>Processing sites: 100%|██████████| 10000/10000 [00:03<00:00, 2771.90it/s]
INFO:PolyAllelicFiltration: Filtered out 0 sites.
INFO:DegeneracyStratification: Number of sites with valid type: 2538
INFO:DegeneracyAnnotation: Annotated 3494 sites.
INFO:Parser: Skipped 401 sites without ancestral allele information.
INFO:Parser: Included 2538 out of 10000 sites in total from the input.

# visualize SFS
p <- spectra$plot()

sfsutils relies here on VCF info tags to determine the degeneracy of a site but this behavior can be customized (cf. DegeneracyStratification).

Stacked stratification

Stratifications split sites into categories, so the parser produces a separate spectrum for each, such as neutral versus selected sites. Several can be combined by passing a list: here we stratify the SFS by degeneracy as well as ancestral base. See the Site stratification reference for the complete list of available stratifications.

# instantiate parser
p <- su$Parser(
  n = 10,
  source = "resources/genome/betula/biallelic.polarized.subset.10000.vcf.gz",
  fasta = "resources/genome/betula/genome.subset.20.fasta",
  gff = "resources/genome/betula/genome.gff.gz",
  annotations = list(
    su$DegeneracyAnnotation()
  ),
  stratifications = list(
    su$DegeneracyStratification(),
    su$AncestralBaseStratification()
  )
)

# parse SFS
spectra <- p$parse()

INFO:Parser: Using stratification: [neutral, selected].[A, C, G, T].
INFO:Parser: Loading VCF file
INFO:Parser: Loading GFF file
INFO:sfsutils: Unzipping resources/genome/betula/genome.gff.gz to /var/folders/n4/m5q2jgw91zv0tp1c4j9bh48m0000gn/T/tmph87iy77_.gff
INFO:Parser: Loading FASTA file
Parser>Counting sites: 10000it [00:00, 390622.03it/s]
Parser>Processing sites: 100%|██████████| 10000/10000 [00:03<00:00, 2719.35it/s]
INFO:PolyAllelicFiltration: Filtered out 0 sites.
INFO:DegeneracyStratification: Number of sites with valid type: 2538
INFO:AncestralBaseStratification: Number of sites with valid type: 2538
INFO:DegeneracyAnnotation: Annotated 3494 sites.
INFO:Parser: Skipped 401 sites without ancestral allele information.
INFO:Parser: Included 2538 out of 10000 sites in total from the input.

# visualize SFS
p <- spectra$plot()