GDSAnnotator 0.99.2
Variant annotation is a critical step in genomic data analysis, facilitating the interpretation of genetic variations and their potential biological significance. As large-scale whole-genome sequencing studies continue to expand, efficient storage, retrieval, and manipulation of variant annotation data have become increasingly important. The Genomic Data Structure (GDS) file format offers a scalable and storage-efficient solution for large genomic datasets, enabling rapid access and computational efficiency.
GDSAnnotator is an R/Bioconductor package designed for fast and
memory-efficient annotation of variants stored in GDS files. It integrates
with Bioconductor tools and utilizes external data resources, such as the
Ensembl Variant Effect Predictor
(VEP),
SnpEff,
gnomAD and
FAVOR (Functional Annotation of Variants Online
Resource), to provide comprehensive functional annotations, including
predicted variant effects and functional characteristics of non-coding
variants. Compared to a plain text file, the GDS format is 30–40 times more
space-efficient when storing whole-genome annotations of FAVOR.
Storing variant annotations in a GDS file was introduced by
FAVOR with the annotated GDS (aGDS)
format: the FAVOR functional annotations of the variants of a study are added
to the SeqArray GDS file that holds the genotypes, so that genotypes and
annotations are read from a single file. GDSAnnotator generalizes the aGDS
format in three directions:
GDSAnnotator builds on gdsfmt and SeqArray:
the annotation files it creates are SeqArray GDS files (with the annotations
stored as INFO fields under annotation/info/, and without genotypes), so
every SeqArray function (seqOpen(), seqGetData(), seqSetFilter(),
seqBlockApply(), …) can be used on them directly. Variants can be supplied
as GenomicRanges GRanges objects, S4Vectors
DataFrame objects, VCF files or VCF objects of
VariantAnnotation (a suggested package, loaded on demand), and
the annotations are returned as a DataFrame.
| Function | Description |
|---|---|
seqFAVOR2GDS |
Convert FAVOR CSV files to GDS |
seqGnomAD2GDS |
Convert gnomAD VCF files to GDS |
seqVEP2GDS |
Convert Ensembl VEP VCF output to GDS |
seqSnpEff2GDS |
Convert SnpEff VCF output to GDS |
seqAnnotate |
Annotate variants using annotations stored in GDS |
seqAnnotateVCF |
Annotate the variants of a VCF file |
seqAnnotateGDS |
Annotate the variants of a SeqArray GDS file |
seqAnnotList |
List the annotations stored in a GDS file |
seqValueCounts |
Count the unique values of annotations |
seqAnnotStat |
Summary statistics of variants and annotations |
seqAnnotReport |
Summary report in HTML, Markdown or R Markdown |
seqAnnotGeneTable |
Per-gene annotation counts, broken down by impact |
GDSAnnotator can be installed from Bioconductor (once released):
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("GDSAnnotator")Load the package:
library(GDSAnnotator)
## Loading required package: gdsfmt
## Loading required package: SeqArrayGDSAnnotator is built on top of the
SeqArray GDS format. Each
annotation GDS file stores per-variant annotations under
annotation/info/<name>. To discover which annotations are available in a
file, use seqAnnotList():
# FAVOR essential dataset, chromosome 22 (subset)
annot_gds <- system.file("extdata", "favor_chr22_sub.gds",
package = "GDSAnnotator")
ann <- seqAnnotList(annot_gds)
ann
## DataFrame with 30 rows and 4 columns
## name type trait description
## <character> <factor> <character> <character>
## 1 apc_conservation Real Float32 NA
## 2 apc_conservation_v2 Real Float32 NA
## 3 apc_epigenetics Real Float32 NA
## 4 apc_epigenetics_active Real Float32 NA
## 5 apc_epigenetics_repr.. Real Float32 NA
## ... ... ... ... ...
## 26 genecode_comprehensi.. Logical Int32 NA
## 27 genehancer Logical Int32 NA
## 28 linsight Real Float32 NA
## 29 cadd_phred Real Float32 NA
## 30 rdhs Logical Int32 NAThe returned DataFrame has four columns: the annotation name, the R data
type, the GDS storage trait, and a textual description (when
available).
The seqAnnotate() generic accepts variants in several forms: a character
vector of chr-pos-ref-alt strings, a data.frame/DataFrame with
chromosome, position, reference and alternative allele columns, a GRanges or
GRangesList, a VCF object of VariantAnnotation, or an open
SeqVarGDSClass object; seqAnnotateVCF() and seqAnnotateGDS() accept a
VCF or GDS file name. The variants are matched by
chromosome, position, reference and alternative allele, and the annotations
are returned as a DataFrame with one row per input variant (NA for the
variants absent from the annotation file).
## Variant IDs are split on `-`, `_` or `:` (the default `split` pattern),
## so all of these formats are accepted:
## "22-10510007-T-C" "22:10510007:T:C" "22_10510007_T_C"
snps <- c(
"22-10510007-T-C",
"22-10510038-T-C",
"22-10510282-G-C",
"22-10510303-G-C" # a fake SNP, returns NA
)
# Use all annotations available in the file
varnm <- ann$name
seqAnnotate(snps, annot_gds, varnm)
## Open 'favor_chr22_sub.gds' [1,000 variants]
## [favor_chr22_sub.gds] # of variants found: 3
## [30] ..............................
## DataFrame with 4 rows and 30 columns
## apc_conservation apc_conservation_v2 apc_epigenetics apc_epigenetics_active
## <numeric> <numeric> <numeric> <numeric>
## 1 2.583288 2.232373 0.0874126 0.226559
## 2 4.766045 4.498693 0.0750365 0.226559
## 3 0.393781 0.323713 0.1819110 0.226559
## 4 NA NA NA NA
## apc_epigenetics_repressed apc_epigenetics_transcription
## <numeric> <numeric>
## 1 0.310309 0.323364
## 2 0.310309 0.323364
## 3 0.310309 0.323364
## 4 NA NA
## apc_local_nucleotide_diversity apc_local_nucleotide_diversity_v2
## <numeric> <numeric>
## 1 4.27664 15.0339
## 2 4.27664 15.0339
## 3 4.27664 15.0271
## 4 NA NA
## apc_local_nucleotide_diversity_v3 apc_mappability apc_micro_rna
## <numeric> <numeric> <numeric>
## 1 15.1012 0.193860 99.4512
## 2 15.1012 0.240450 99.4512
## 3 15.0967 0.184841 99.4512
## 4 NA NA NA
## apc_mutation_density apc_protein_function_v3 apc_proximity_to_coding
## <numeric> <numeric> <numeric>
## 1 14.9277 2.96949 0.136519
## 2 14.9277 2.96949 0.136519
## 3 14.9188 2.96949 0.136519
## 4 NA NA NA
## apc_proximity_to_coding_v2 apc_proximity_to_tsstes apc_transcription_factor
## <numeric> <numeric> <numeric>
## 1 15.1786 0.294951 3.14279
## 2 15.1786 0.295003 3.14279
## 3 15.1786 0.295432 3.14279
## 4 NA NA NA
## cage_promoter cage_tc metasvm_pred rsid fathmm_xf
## <logical> <logical> <logical> <character> <numeric>
## 1 NA NA NA 0.405274
## 2 NA NA NA 0.405274
## 3 NA NA NA 0.405274
## 4 NA NA NA NA NA
## genecode_comprehensive_category genecode_comprehensive_info
## <character> <character>
## 1 intergenic NONE(dist=NONE),RF00..
## 2 intergenic NONE(dist=NONE),RF00..
## 3 intergenic NONE(dist=NONE),RF00..
## 4 NA NA
## genecode_comprehensive_exonic_category genecode_comprehensive_exonic_info
## <logical> <logical>
## 1 NA NA
## 2 NA NA
## 3 NA NA
## 4 NA NA
## genehancer linsight cadd_phred rdhs
## <logical> <numeric> <numeric> <logical>
## 1 NA 0.214926 7.712 NA
## 2 NA 0.214926 10.460 NA
## 3 NA 0.214926 3.944 NA
## 4 NA NA NA NAA subset of annotations can be selected by passing only the names of
interest in varnm:
seqAnnotate(snps, annot_gds,
varnm = c("apc_conservation", "cadd_phred", "linsight"))
## Open 'favor_chr22_sub.gds' [1,000 variants]
## [favor_chr22_sub.gds] # of variants found: 3
## [3] ...
## DataFrame with 4 rows and 3 columns
## apc_conservation cadd_phred linsight
## <numeric> <numeric> <numeric>
## 1 2.583288 7.712 0.214926
## 2 4.766045 10.460 0.214926
## 3 0.393781 3.944 0.214926
## 4 NA NA NAdata.framedf <- data.frame(
chr = c("22", "22", "22"),
pos = c(10510007L, 10510038L, 10510282L),
ref = c("T", "T", "G"),
alt = c("C", "C", "C")
)
df
## chr pos ref alt
## 1 22 10510007 T C
## 2 22 10510038 T C
## 3 22 10510282 G C
seqAnnotate(df, annot_gds, varnm = c("apc_conservation", "cadd_phred"))
## Open 'favor_chr22_sub.gds' [1,000 variants]
## [favor_chr22_sub.gds] # of variants found: 3
## [2] ..
## DataFrame with 3 rows and 2 columns
## apc_conservation cadd_phred
## <numeric> <numeric>
## 1 2.583288 7.712
## 2 4.766045 10.460
## 3 0.393781 3.944GRanges objectGRanges queries return all variants that fall within the specified ranges:
gr <- GenomicRanges::GRanges(
seqnames = c("chr22", "chr22"),
ranges = IRanges::IRanges(
start = c(10510100, 10510225),
end = c(10510125, 10510250)))
seqAnnotate(gr, annot_gds, varnm = c("apc_conservation", "cadd_phred"))
## Open 'favor_chr22_sub.gds' [1,000 variants]
## [2] ..
## DataFrame with 156 rows and 2 columns
## apc_conservation cadd_phred
## <numeric> <numeric>
## 1 1.68708 5.807
## 2 1.68708 5.447
## 3 1.68708 5.203
## 4 6.78621 8.977
## 5 6.78621 8.606
## ... ... ...
## 152 2.45621 6.684
## 153 2.45621 6.992
## 154 2.55586 6.566
## 155 2.55586 6.772
## 156 2.55586 7.076VCF objectA VCF object of the VariantAnnotation package
(CollapsedVCF or ExpandedVCF) is annotated with one row per alternative
allele, and seqAnnotateVCF() does the same for a VCF file, which it reads
with VariantAnnotation::readVcf(). VariantAnnotation is a suggested
package, loaded only when a VCF file or object is annotated. Here the
variants of the gnomAD example VCF are looked up in the FAVOR subset:
vcf_fn <- system.file("extdata", "gnomad.genomes.v4.chr22_sub.vcf.gz",
package = "GDSAnnotator")
vcf <- VariantAnnotation::readVcf(vcf_fn)
dim(vcf)
## [1] 150 0
res <- seqAnnotate(vcf, annot_gds,
varnm = c("rsid", "cadd_phred", "linsight"), verbose = FALSE)
table(found = !is.na(res$cadd_phred))
## found
## FALSE TRUE
## 109 41
head(res[!is.na(res$cadd_phred), ])
## DataFrame with 6 rows and 3 columns
## rsid cadd_phred linsight
## <character> <numeric> <numeric>
## 1 7.847 0.214926
## 2 7.898 0.214926
## 3 rs1290354662 4.273 0.214926
## 4 8.553 0.214926
## 5 rs1325858619 4.889 0.214926
## 6 rs1435837414 4.419 0.214926
# the same, directly from the file name
res2 <- seqAnnotateVCF(vcf_fn, annot_gds,
varnm = c("rsid", "cadd_phred", "linsight"), verbose = FALSE)
identical(res, res2)
## [1] TRUEFor whole-genome resources the annotation is expected to be split into one
GDS file per chromosome (this is how FAVOR and gnomAD are distributed), and
annot_gds accepts a character vector of file names. Each file is opened
once, the query variants are dispatched to the file of their chromosome, and
the results are returned in the order of the input:
annot_gds <- sprintf("favor_chr%s.gds", c(seq_len(22L), "X"))
seqAnnotate(snps, annot_gds, varnm = c("cadd_phred", "linsight"))When the object being annotated is a SeqVarGDSClass genotype GDS file,
annotations can be written directly into the same file (under
annotation/info/) by setting add_to_gds = TRUE, or into a new GDS file by
passing a file name. This avoids materialising a large annotation
DataFrame in R memory for whole-genome studies.
library(SeqArray)
# An existing genotype GDS file
geno_gds <- "my_wgs_geno.gds"
seqAnnotate(seqOpen(geno_gds), annot_gds, varnm,
add_to_gds = "my_wgs_geno_annotated.gds")GDSAnnotator provides dedicated importers for common annotation sources.
All of them output a SeqArray-compatible GDS file, which can subsequently be
queried with seqAnnotate().
vcf_fn <- system.file("extdata", "gnomad.genomes.v4.chr22_sub.vcf.gz",
package = "GDSAnnotator")
out_fn <- tempfile(fileext = ".gds")
seqGnomAD2GDS(vcf_fn, out_fn, verbose = FALSE)
unlink(out_fn, force = TRUE)The VEP importer parses the compound CSQ INFO field emitted by Ensembl VEP
into individual annotation fields in the GDS file.
vcf_fn <- system.file("extdata", "example_wgs_sites_chr22_vep.vcf.gz",
package = "GDSAnnotator")
out_fn <- tempfile(fileext = ".gds")
seqVEP2GDS(vcf_fn, out_fn, verbose = FALSE)
unlink(out_fn, force = TRUE)The SnpEff importer parses the ANN, LOF and NMD INFO fields by
default:
vcf_fn <- system.file("extdata", "example_wgs_sites_chr22_snpeff.vcf.gz",
package = "GDSAnnotator")
out_fn <- tempfile(fileext = ".gds")
seqSnpEff2GDS(vcf_fn, out_fn, verbose = FALSE)
unlink(out_fn, force = TRUE)The resulting GDS file can be inspected with seqAnnotList() and queried with
seqAnnotate() just like the FAVOR example above.
The two files converted above are also shipped with the package, as
example_wgs_sites_chr22_vep.gds and example_wgs_sites_chr22_snpeff.gds,
and are used in the examples that follow.
seqValueCounts() is a convenience wrapper that tabulates unique values of
one or more annotation fields without loading the entire annotation into R:
gds_fn <- system.file("extdata", "gnomad.genomes.v4.chr22_sub.gds",
package = "GDSAnnotator")
f <- SeqArray::seqOpen(gds_fn)
seqValueCounts(f, "variant_type")
##
[..................................................] 0%, ETC: ---
[==================================================] 100%, used 0s
## variant_type
## indel mixed multi-snv snv
## 4 23 54 69
seqValueCounts(f, c("variant_type", "allele_type"))
##
[..................................................] 0%, ETC: ---
[==================================================] 100%, used 0s
## allele_type
## variant_type del ins snv
## indel 3 1 0
## mixed 10 3 10
## multi-snv 0 0 54
## snv 0 0 69
SeqArray::seqClose(f)Set per_variant = TRUE to collapse multi-valued annotations (e.g. multiple
VEP consequences per variant) to a single count per variant, and
parallel = TRUE to enable multi-core processing for large files.
Ensembl VEP and SnpEff each write an HTML summary next to their VCF output,
but that report is lost once the annotations have been converted to GDS,
and it is not available at all for a GDS file assembled from several sources.
seqAnnotStat() and seqAnnotReport() regenerate an equivalent report
directly from the GDS file.
seqAnnotStat() collects every statistic in a single pass over the file,
so the cost does not grow with the number of statistics reported:
gds_fn <- system.file("extdata", "example_wgs_sites_chr22_snpeff.gds",
package = "GDSAnnotator")
st <- seqAnnotStat(gds_fn, verbose = FALSE)
st
## SeqAnnotStat [SnpEff]
## file: example_wgs_sites_chr22_snpeff.gds
## # of variants: 200
## # of annotation records: 1,801
## statistics: chrom, pos_min, pos_max, n_allele, var_class, base_change, indel_len, known, qual_na, n_annot, cons_all, cons_severe, impact, impact_severe, region, func_class, biotype, feature_type, feature, errors, gene_impactEach statistic in st$counts is a named vector of counts. Consequences are
counted twice: once for every annotation record (one record per transcript,
cons_all), and once per variant keeping only the most severe consequence
(cons_severe), following the Sequence Ontology ranking used by VEP:
st$counts$cons_severe
## frameshift_variant splice_acceptor_variant splice_donor_variant
## 93 20 18
## start_lost stop_gained
## 17 52
st$counts$impact
## HIGH MODERATE LOW MODIFIER
## 1140 79 41 541seqAnnotGeneTable() returns the per-gene counts, the equivalent of the
genes.txt file written by SnpEff:
head(seqAnnotGeneTable(st))
## gene total HIGH MODERATE LOW MODIFIER
## 1 CLTCL1 831 660 44 21 106
## 2 BCL2L13 109 28 0 0 81
## 3 GAB4 96 39 7 0 50
## 4 ADA2 93 55 16 0 22
## 5 CDC45 75 43 0 0 32
## 6 SLC25A18 70 55 1 0 14seqAnnotReport() renders the statistics as a report. The output format is
taken from the file extension by default: .md gives GitHub-flavored
Markdown, .Rmd gives R Markdown source, and anything else gives a
self-contained HTML file that needs no JavaScript and no additional package,
so it can be viewed offline. None of the three writes a companion image
directory, so the report is always a single file.
md_fn <- tempfile(fileext = ".md")
seqAnnotReport(st, md_fn, title = "SnpEff annotation of chromosome 22",
verbose = FALSE)Since the report is written in Markdown, it can be embedded in an R Markdown
document such as this vignette: print its lines in a chunk with the option
results = "asis", and they are rendered as part of the document. The first
three sections of the report are shown below, with the headings demoted by
two levels so that they fit under this subsection:
txt <- readLines(md_fn)
txt <- txt[seq_len(grep("^## ", txt)[4L] - 1L)] # the title and 3 sections
cat(sub("^#", "###", txt), sep = "\n")| Item | Value |
|---|---|
| Report generated | 2026-10-23 18:49:51 |
| GDS file | example_wgs_sites_chr22_snpeff.gds |
| File size | 0.1 MB |
| Annotation source | SnpEff |
| SnpEffVersion | "5.2f (build 2025-02-07 08:36), by Pablo Cingolani" |
| SnpEffCmd | "SnpEff -lof -s ./snpEff_summary_chr22.html GRCh38.p14 c22_b0_v1.sites.vcf.gz" |
| source | Minimac4.v1.0.2 |
| # of variants | 200 |
| # of annotation records | 1,801 |
| # of genes with annotation | 56 |
| # of features (transcripts) | 218 |
| Genomic span covered | 9.2 Mb |
| Variant rate | one variant every 45,975 bases |
| Chromosome | Count | Percentage | |
|---|---|---|---|
| 22 | 200 | 100.00% | ████████████████████████ |
| Category | Count | Percentage | |
|---|---|---|---|
| SNV | 100 | 50.00% | ████████████████████████ |
| deletion | 65 | 32.50% | ████████████████ |
| insertion | 35 | 17.50% | ████████ |
The whole report is embedded in the same way by printing all of its lines.
The Markdown charts are drawn with block characters. If the charts should be
real plots instead, e.g., to restyle them or to fold the report into a larger
document, the rmd format writes an R Markdown source in which the counts are
embedded as R code and the charts are produced by barplot() at knit time:
seqAnnotReport(st, "annot_summary.Rmd",
title = "SnpEff annotation of chromosome 22")
rmarkdown::render("annot_summary.Rmd")The counters are additive, and that has two consequences. First, the variants
can be split between processes and the per-process counters merged
afterwards, so parallel gives a summary identical to the serial one:
st <- seqAnnotStat(gds_fn, parallel = 4L)Second, a whole-genome summary can be built from per-chromosome GDS files without ever holding more than one of them open. The memory usage depends on the number of distinct categories (e.g. the number of genes), not on the number of variants:
gds_fn <- sprintf("annot_chr%d.gds", seq_len(22L))
st <- seqAnnotStat(gds_fn)
seqAnnotReport(st, "wgs_summary.html", title = "Whole-genome annotation")sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] GDSAnnotator_0.99.2 SeqArray_1.53.5 gdsfmt_1.49.9
## [4] BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] KEGGREST_1.53.6 SummarizedExperiment_1.43.0
## [3] rjson_0.2.23 xfun_0.62
## [5] bslib_0.12.0 Biobase_2.73.2
## [7] lattice_0.23-1 vctrs_0.7.3
## [9] tools_4.6.1 bitops_1.1-0
## [11] generics_0.1.4 stats4_4.6.1
## [13] curl_8.1.0 parallel_4.6.1
## [15] AnnotationDbi_1.75.2 RSQLite_3.53.3
## [17] blob_1.3.0 BiocBaseUtils_1.15.1
## [19] Matrix_1.7-6 BSgenome_1.81.1
## [21] S4Vectors_0.51.10 cigarillo_1.3.1
## [23] lifecycle_1.0.5 compiler_4.6.1
## [25] Rsamtools_2.29.0 Biostrings_2.81.9
## [27] Seqinfo_1.3.2 codetools_0.2-20
## [29] htmltools_0.5.9 sass_0.4.10
## [31] RCurl_1.98-1.20 yaml_2.3.12
## [33] crayon_1.5.3 jquerylib_0.1.4
## [35] BiocParallel_1.47.0 DelayedArray_0.39.8
## [37] cachem_1.1.0 abind_1.4-8
## [39] digest_0.6.39 restfulr_0.0.17
## [41] bookdown_0.48 VariantAnnotation_1.59.4
## [43] fastmap_1.2.0 grid_4.6.1
## [45] cli_3.6.6 SparseArray_1.13.4
## [47] S4Arrays_1.13.2 GenomicFeatures_1.65.0
## [49] XML_3.99-0.25 bit64_4.8.6
## [51] rmarkdown_2.32 XVector_0.53.0
## [53] httr_1.4.9 matrixStats_1.5.0
## [55] bit_4.6.0 otel_0.2.0
## [57] png_0.1-9 memoise_2.0.1
## [59] evaluate_1.0.5 knitr_1.52
## [61] GenomicRanges_1.65.4 IRanges_2.47.5
## [63] BiocIO_1.23.3 rtracklayer_1.73.0
## [65] rlang_1.3.0 DBI_1.3.0
## [67] BiocManager_1.30.27 BiocGenerics_0.59.12
## [69] jsonlite_2.0.0 R6_2.6.1
## [71] MatrixGenerics_1.25.0 GenomicAlignments_1.49.2