Site annotation

sfsutils can add site-level annotations such as the ancestral allele, site degeneracy, and synonymy. The Parser applies them on the fly while it builds a spectrum, while the Annotator applies them and writes the annotated variants to a file. Polarising against an ancestral state is only needed for an unfolded spectrum: the Parser reads it from the AA field by default (this can be customised), and otherwise the spectrum can be folded.

Degeneracy Annotation

DegeneracyAnnotation annotates the SFS by the degeneracy of the site. This annotation requires information from a FASTA and GFF file and is useful for stratifying the SFS by 0-fold and 4-fold degenerate sites, a common way of contrasting putatively neutral and selected sites (see DegeneracyStratification).

library(sfsutils)

su <- load_sfsutils()

# example for degeneracy annotation
ann <- su$Annotator(
  source = "resources/genome/betula/biallelic.subset.10000.vcf.gz",
  fasta = "resources/genome/betula/genome.subset.20.fasta",
  gff = "resources/genome/betula/genome.gff.gz",
  annotations = list(su$DegeneracyAnnotation()),
  output = "genome.deg.vcf.gz"
)

ann$annotate()

INFO:Annotator: Start annotating
INFO:Annotator: Loading GFF file
INFO:sfsutils: Unzipping resources/genome/betula/genome.gff.gz to /var/folders/n4/m5q2jgw91zv0tp1c4j9bh48m0000gn/T/tmpkqfmbk5g.gff
INFO:Annotator: Loading FASTA file
INFO:Annotator: Loading VCF file
Annotator>Counting sites: 10000it [00:00, 480887.87it/s]
Annotator>Processing sites: 100%|██████████| 10000/10000 [00:00<00:00, 18445.35it/s]
INFO:DegeneracyAnnotation: Annotated 3566 sites.

Ancestral Allele Annotation

Currently, two ancestral allele annotations are available: MaximumParsimonyAncestralAnnotation and MaximumLikelihoodAncestralAnnotation. The former is error-prone and not recommended. Alternatively, if outgroups are missing, the spectra can be folded, though this discards information and yields a less informative spectrum. Ideally, we would like to use MaximumLikelihoodAncestralAnnotation, which is more sophisticated and requires one or several outgroup to be specified. Its underlying model is based on EST-SFS. The maximum-likelihood model estimates branch rates from monomorphic sites, so ideally these are present in the input. When there are few or none, supply n_target_sites, the total number of sites (segregating and monomorphic) underlying the data, along with a FASTA reference from which the monomorphic sites are sampled.

# example for ancestral allele annotation with outgroups
ann <- su$Annotator(
  source = "resources/genome/betula/all.with_outgroups.subset.10000.vcf.gz",
  fasta = "resources/genome/betula/genome.subset.20.fasta",
  annotations = list(su$MaximumLikelihoodAncestralAnnotation(
    outgroups = list("ERR2103730"),
    n_ingroups = 10,
    n_target_sites = 200000
  )),
  output = "genome.aa.vcf.gz"
)

ann$annotate()

INFO:Annotator: Start annotating
INFO:Annotator: Loading VCF file
INFO:Annotator: Loading VCF file
INFO:MaximumLikelihoodAncestralAnnotation: Subsampling 10 ingroup haplotypes probabilistically from 378 individuals in total.
INFO:MaximumLikelihoodAncestralAnnotation: Using 1 outgroup samples (ERR2103730).
Annotator>Counting sites: 10000it [00:00, 599880.43it/s]
MaximumLikelihoodAncestralAnnotation>Parsing sites: 100%|██████████| 10000/10000 [00:06<00:00, 1591.03it/s]
INFO:MaximumLikelihoodAncestralAnnotation: Sampling mono-allelic sites.
MaximumLikelihoodAncestralAnnotation>Sampling mono-allelic sites:   0%|          | 0/100000 [00:00<?, ?it/s]INFO:Annotator: Loading FASTA file
MaximumLikelihoodAncestralAnnotation>Sampling mono-allelic sites: 100%|██████████| 100000/100000 [00:00<00:00, 589771.84it/s]
INFO:MaximumLikelihoodAncestralAnnotation: Extrapolating to 127090 mutational target sites based on the number of sites parsed.
INFO:MaximumLikelihoodAncestralAnnotation: Included 133779 sites for the inference (1027 polymorphic, 132752 monomorphic).
MaximumLikelihoodAncestralAnnotation>Optimizing rates: 100%|██████████| 10/10 [00:00<00:00, 12.14it/s]
Annotator>Processing sites: 100%|██████████| 10000/10000 [00:02<00:00, 4686.13it/s]
INFO:MaximumLikelihoodAncestralAnnotation: Annotated 10000 sites.
INFO:MaximumLikelihoodAncestralAnnotation: There were 195 mismatches between the most likely ancestral allele and the ad-hoc ancestral allele annotation.