Skip to contents

Annotate Variants with Nearest Gene and High-Resolution Genomic Context

Source: R/04_annotation.R
av_annotate_genes.Rd

Maps the exact 1-bp variant position (parsed from obj$data$variant_id, format chr:pos:REF>ALT) to a reference transcript database (TxDb). Bioconductor lookups run on the deduplicated set of variant positions, so AlphaGenome inputs with many tracks per variant only pay one lookup cost per unique variant. The four annotation columns are joined back onto every row of $data:

  • nearest_gene: The Entrez ID of the closest gene.

  • gene_symbol: The HGNC symbol (resolved via org.Hs.eg.db). Falls back to Entrez ID if the annotation package is unavailable.

  • dist_to_tss: The distance in base pairs to the nearest TSS, calculated from the exact 1-bp TSS coordinate respecting strand orientation (+ strand: gene start; - strand: gene end).

  • region_type: Hierarchical classification using winner-takes-all logic: "cds" > "5utr" > "3utr" > "intron" > "promoter" > "intergenic".

The four whole-genome TxDb extractions (genes, CDS, 5' UTR, 3' UTR) depend on the annotation database only, so they are cached per TxDb for the duration of the R session. The first call pays the full lookup cost, later calls do not.

Usage

av_annotate_genes(
  obj,
  txdb_package = "TxDb.Hsapiens.UCSC.hg38.knownGene",
  promoter_upstream = 2000,
  promoter_downstream = 200
)

Arguments

obj

An object of class AlphaVarSet. Must carry a variant_id column in the AlphaGenome format chr:pos:REF>ALT.

txdb_package

Character string. Name of the TxDb package to use (default: "TxDb.Hsapiens.UCSC.hg38.knownGene"). Can also be a TxDb object directly.

promoter_upstream

Numeric. Distance upstream of the TSS to define the promoter region (default: 2000).

promoter_downstream

Numeric. Distance downstream of the TSS to define the promoter region (default: 200).

Value

An AlphaVarSet object with columns nearest_gene, gene_symbol, dist_to_tss, and region_type joined onto $data by variant_id.

Details

Assigns biological context to variants by identifying the nearest gene (TSS), resolving gene symbols, calculating the strand-aware distance to the Transcription Start Site (TSS), and classifying the genomic region using a fine-grained hierarchical scheme.

Examples

example_file <- system.file(
  "extdata", "bcl11a_rs1427407_example.csv",
  package = "AlphaVaR"
)
obj <- av_import_csv(example_file)
#> Reading bcl11a_rs1427407_example.csv ...

# TxDb.Hsapiens.UCSC.hg38.knownGene is a ~1 GB Suggests package.
if (requireNamespace("TxDb.Hsapiens.UCSC.hg38.knownGene", quietly = TRUE)) {
  obj <- av_annotate_genes(
    obj,
    txdb_package = "TxDb.Hsapiens.UCSC.hg38.knownGene"
  )

  print(table(obj$data$region_type))

  # One annotation row per variant, repeated across all its tracks.
  print(unique(
    obj$data[, c("variant_id", "gene_symbol", "dist_to_tss", "region_type")]
  ))
}
#> Extracting gene information from annotation database...
#>   2169 genes were dropped because they have exons located on both strands of
#>   the same reference sequence or on more than one reference sequence, so cannot
#>   be represented by a single genomic range.
#>   Use 'single.strand.genes.only=FALSE' to get all the genes in a GRangesList
#>   object, or use suppressMessages() to suppress this message.
#> Calculating strand-aware distances to nearest TSS...
#> Resolving gene symbols...
#> 
#> 'select()' returned 1:1 mapping between keys and columns
#> Classifying genomic regions (hierarchical)...
#> Annotation complete.
#> 
#> intron 
#>   2676 
#> # A tibble: 2 × 4
#>   variant_id        gene_symbol  dist_to_tss region_type
#>   <chr>             <chr>              <dbl> <chr>      
#> 1 chr2:60490908:T>C LOC102724142       51066 intron     
#> 2 chr2:60490908:T>G LOC102724142       51066 intron