Annotate Variants with FAVOR Database Scores

Source: R/favor_annotator.R

Annotates variants with functional annotation scores from the FAVOR (Functional Annotation of Variants Online Resource) database, read either from the FAVOR v1 CSV chunk files or from the FAVOR 2.0 per-chromosome Parquet files. Every input row receives a matching outcome (match_tier), when matched the key of the database row that supplied its annotation (favor_key), and the rsID evidence for that row (rsid_check). Used in PI estimation to retrieve pathogenicity-relevant features for trait-associated and control variants, and to build the cohort aGDS that GLOW and STAAR scan.

annotate_favor(
  variants,
  favor_db_path,
  favor_split_file = NULL,
  features = .default_favor_features(),
  output_csv = NULL,
  output_agds = NULL,
  update_gds = FALSE,
  variant_filter = NULL,
  match_method = "exact",
  na_allele_method = "average",
  na_handling = "keep",
  use_xsv = TRUE,
  verbose = 1,
  favor_db_format = c("auto", "csv", "parquet"),
  favor_release = NULL,
  rsid_policy = c("require", "record")
)

Arguments

variants

One of the following:

  • A data.frame with VarInfo column (CHR-POS-REF-ALT format), and optionally an rsID column (the input's rsID per row), which feeds the rsID check described in Details

  • An S3 object of class glow_pi_case_data or glow_pi_control_data from prepare_PI_case_data or prepare_PI_control_data

  • A path (character) to a GDS file containing variant data (requires SeqArray package). Variants are extracted automatically, and the GDS annotation/id (the rsID on a chip cohort) is read as the rsID column.

favor_db_path

Character. The database directory: either the FAVOR v1 CSV chunk files (chr1_1.csv, chr1_2.csv, ...) or the FAVOR 2.0 Parquet files (chromosome_1.parquet, ...).

favor_split_file

Character or NULL. CSV backend only. Path to FAVORdatabase_chrsplit.csv, which maps positions to chunk files. If NULL (default), searches in this order: (1) favor_db_path, (2) the package's bundled file.

features

Character vector. Annotation fields to extract, named as the FAVOR v1 columns. The default is the complete annotation content of the FAVOR Essential Database (30 fields): 20 numeric scores (17 annotation principal components, CADD PHRED, LINSIGHT, FATHMM-XF) plus 10 categorical fields (GENCODE category/info and exonic category/info, MetaSVM prediction, GeneHancer, CAGE, rDHS, rsID). The categorical fields are what the built-in coding masks read, so an aGDS built with the default is directly scannable by every built-in variant category. Fields the database does not carry are skipped and reported (FAVOR 2.0 lacks four older APC versions, so the default yields 26 fields there).

output_csv

Character or NULL. If provided, saves annotated results to CSV file at this path

output_agds

Character or NULL. If provided, saves annotated results to aGDS file at this path (requires SeqArray package)

update_gds

Logical. If TRUE and input is a GDS file, write annotations back to the input GDS file, converting it to an aGDS format. Default is FALSE. Ignored if variants is not a GDS file path.

variant_filter

List or NULL. Optional filter to select a subset of variants. Works with all input types (GDS, data.frame, S3 objects). Supported filters:

  • list(chr = "21"): Filter by chromosome (works with all inputs)

  • list(chr = "21", start = 1e6, end = 5e6): Filter by position range (works with all inputs)

  • list(variant_ids = c(1, 5, 10)): Filter by SeqArray variant IDs (GDS input only)

Default is NULL (use all variants). For non-GDS input, chromosome and position are parsed from the VarInfo column.

match_method

Character. How allele-bearing variants are matched:

  • "exact" (default): the key as given only (STAAR-compatible)

  • "flexible": for SNVs, the lookup key is also rewritten against FAVOR's reference base by swapping REF and ALT or complementing both alleles, for cohorts whose alleles may differ from the reference in order or strand (genotyping-chip data). Indels always match by the exact key. See Details.

na_allele_method

Character. How position-only keys (CHR-POS-NA-NA) are annotated from the SNV rows at the position (indel rows are never included; a position with no SNV row is "uncovered"):

  • "average" (default): numeric features averaged over the SNV rows' non-missing values, string features take the first non-empty value

  • "first": every feature from the first SNV row at the position

na_handling

Character. How to handle missing annotation values after matching (never changes match_tier):

  • "keep" (default): Keep NA values as-is

  • "zero": Replace NA with 0 in numeric features

  • "drop": Remove variants with any NA annotations from the returned table and the CSV. The aGDS writers keep every variant of the GDS (the file stays a copy of the cohort) and leave the dropped variants' annotation and tier entries empty.

use_xsv

Logical. CSV backend only. If TRUE (default) and the xsv CLI tool is available, xsv fetches the database rows at the input positions (a streaming join on position); otherwise each chunk is read with data.table::fread() and filtered. Both give the same rows.

verbose

Integer. Verbosity level: 0=silent, 1=messages (default), 2=detailed messages

favor_db_format

Character. "auto" (default), "csv" or "parquet". Under "auto" the format is detected from the files in favor_db_path; a directory holding both kinds, or neither, is an error. (Placed after verbose so that positional calls written for GLOWr 0.1.1 keep their meaning.)

favor_release

Character or NULL. The FAVOR release the database files come from, recorded in the aGDS provenance. NULL records "unknown".

rsid_policy

Character. What the rsID check does with a match:

  • "require" (default): a swapped, flipped or flipped_swapped match whose check is differs or favor_none is withheld as the non-match rsid_conflict (its annotation stays missing; favor_key names the withheld row). Exact matches and matches checked chip_none are kept, so exact matching is unaffected and an input without rsIDs loses nothing.

  • "record": every match is kept and the check is only recorded in rsid_check, for diagnostics.

Value

data.frame with the input columns (for GDS input also variant_id and rsID), one column per served feature, match_tier (the matching outcome of each row), favor_key (the variant_vcf of the database row that supplied the annotation, or of the row withheld under rsid_conflict; NA otherwise) and rsid_check (the rsID evidence of every matched row; NA for non-matches and position-only keys). Order of variants is preserved from input.

Details

Database formats:

  • FAVOR v1 CSV: chr{N}_{K}.csv chunk files keyed by variant_vcf ("CHR-POS-REF-ALT"), with FAVORdatabase_chrsplit.csv mapping positions to chunks (a copy is bundled with GLOWr).

  • FAVOR 2.0 Parquet: chromosome_{N}.parquet, one file per chromosome (1-22, X, Y), sorted by position, with the annotation in nested struct columns, read leaf by leaf one row group at a time. Requires the arrow package. The map from our field names to the Parquet sources is internal (.favor2_field_map()); GeneHancer, GENCODE exonic information and GENCODE information are serialized into the v1 string layouts. FAVOR 2.0 stores scores as float32, which R reads as the exact double value of the float32.

The users obtain the database from the FAVOR team; GLOWr does not redistribute it.

Matching (the lookup key is rewritten, never the data): Each input row is classified as an SNV (REF and ALT each one of A, C, G, T), an indel (A/C/G/T strings, not both of length one), a position-only key (both alleles "NA"), or unsupported (anything else, including symbolic alleles and multiallelic records). For an SNV at position p with FAVOR reference base r:

  1. "exact": the key CHR-POS-REF-ALT exists.

  2. "swapped" (flexible): r equals ALT and CHR-POS-ALT-REF exists.

  3. "flipped" / "flipped_swapped" (flexible): r equals the complement of REF (or of ALT) and the complemented key exists.

Indels match by the exact key only, because a VCF indel's REF is the genome sequence and its reciprocal key names a different variant. Non-matches carry a reason: "unmatched_alt" (the position is covered and the reference base is found, either by REF as given or, under flexible matching, after one of the enabled transformations, but no row carries the resulting REF-ALT pair), "unmatched_ref" (the position is covered, REF as given does not agree with r, and no enabled transformation finds r), "uncovered" (no database row at the position), "no_database" (the database has no file or chunk for the chromosome), "unsupported_allele" and "rsid_conflict" (below). Position-only keys get "position_only", with the features aggregated over the SNV rows at the position only (na_allele_method); a position that carries no SNV row is "uncovered". A tier records the lookup transformation, not a verified strand: a non-palindromic SNV pair at a covered position always resolves through the first four tiers, including a pair of two non-reference alleles. The cohort's own alleles, genotypes and keys are never modified.

The rsID as evidence, never a key. FAVOR carries dbSNP's rsID on the row of the substitution dbSNP lists, and the other rows at the position carry none, so agreement between the input's rsID and the matched row's confirms the site and the allele pair. For every matched row rsid_check records "same", "differs", "chip_none" (the input carries no rsID of the form rs<digits>) or "favor_none" (the database row carries none). The input's rsID is the GDS annotation/id or the data.frame's rsID column; without either, every check is "chip_none". Under rsid_policy = "require", the default, a transformed match checked "differs" or "favor_none" becomes "rsid_conflict"; "record" keeps it. The check reads the database's rsID column whether or not rsid is a requested feature.

Outputs: the aGDS writers align rows to the GDS through variant.id, write one native-typed sub-node per feature under annotation/info/FunctionalAnnotation, the matching outcome as annotation/info/favor_match_tier, the rsID check as annotation/info/favor_rsid_check, and provenance (database format, files, sizes and SHA-256 when a SHA256SUMS file is present, release, GLOWr version, matching settings, the rsID policy and source, the rows the position-only aggregation uses, outcome and check counts, date) as attributes of the FunctionalAnnotation folder.

Computational Complexity: the matching is vectorized and linear in the number of input rows plus the number of database rows at their positions. Reading dominates: the CSV backend streams each chunk that covers the input positions; the Parquet backend reads the position column of every row group and the requested leaf columns of the groups that hold input positions, one group at a time.

References

Zhou, H., Arapoglou, T., Li, X., et al. (2023). FAVOR: functional annotation of variants online resource and annotator for variation across the human genome. Nucleic Acids Research, 51(D1), D1300-D1311. doi:10.1093/nar/gkac966

Zhou, H., Verma, V., Li, X., et al. (2026). FAVOR 2.0: A reengineered functional annotation of variants online resource for interpreting genomic variation. Nucleic Acids Research, 54(D1), D1405-D1414. doi:10.1093/nar/gkaf1217

See also

prepare_PI_case_data for case variant preparation prepare_PI_control_data for control variant preparation get_PI for PI estimation using annotated variants

Examples

if (FALSE) { # \dontrun{
# Annotate a chip cohort's GDS from FAVOR 2.0, rewriting the lookup key for
# allele order and strand, and withholding transformed matches that the
# chip's rsID does not confirm
annotate_favor(
  variants = "data/chr22.gds",
  favor_db_path = "data/favor2-db/FAVOR2.0",
  match_method = "flexible",
  rsid_policy = "require",
  output_agds = "results/chr22_favor2.gds"
)

# Annotate case variants given by position only, from FAVOR v1 CSV
annotated_cases <- annotate_favor(
  variants = cases,
  favor_db_path = "data/FAVOR",
  na_allele_method = "average"
)
table(annotated_cases$match_tier)

# GDS input with a position range filter
annotated <- annotate_favor(
  variants = "data/genotypes.gds",
  favor_db_path = "data/FAVOR",
  variant_filter = list(chr = "21", start = 1e6, end = 5e6)
)
} # }