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 |
|---|---|
seqToGDS_FAVOR |
Convert FAVOR CSV files to GDS |
seqToGDS_gnomAD |
Convert gnomAD VCF files to GDS |
seqToGDS_VEP |
Convert Ensembl VEP VCF output to GDS |
seqToGDS_SnpEff |
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:
GDSAnnotator 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:
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.
GDSAnnotator provides dedicated importers for common
annotation sources. All of them output a SeqArray-compatible GDS file,
which can subsequently be queried with seqAnnotate().
The VEP importer parses the compound CSQ INFO field
emitted by Ensembl VEP into individual annotation fields in the GDS
file.
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")
seqToGDS_SnpEff(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-09-27 08:29:12 |
| 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:
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:
sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 26.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [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: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] GDSAnnotator_0.99.1 SeqArray_1.53.3 gdsfmt_1.49.8
## [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.61
## [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.0.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] cigarillo_1.3.1 S4Vectors_0.51.10
## [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 sys_3.4.3
## [31] buildtools_1.0.0 sass_0.4.10
## [33] RCurl_1.98-1.20 yaml_2.3.12
## [35] crayon_1.5.3 jquerylib_0.1.4
## [37] BiocParallel_1.47.0 DelayedArray_0.39.7
## [39] cachem_1.1.0 abind_1.4-8
## [41] digest_0.6.39 restfulr_0.0.17
## [43] VariantAnnotation_1.59.4 maketools_1.3.2
## [45] fastmap_1.2.0 grid_4.6.1
## [47] cli_3.6.6 SparseArray_1.13.3
## [49] S4Arrays_1.13.1 GenomicFeatures_1.65.0
## [51] XML_3.99-0.24 bit64_4.8.6
## [53] rmarkdown_2.32 XVector_0.53.0
## [55] httr_1.4.9 matrixStats_1.5.0
## [57] bit_4.6.0 otel_0.2.0
## [59] png_0.1-9 memoise_2.0.1
## [61] evaluate_1.0.5 knitr_1.52
## [63] GenomicRanges_1.65.4 IRanges_2.47.5
## [65] BiocIO_1.23.3 rtracklayer_1.73.0
## [67] rlang_1.3.0 DBI_1.3.0
## [69] BiocManager_1.30.27 BiocGenerics_0.59.12
## [71] jsonlite_2.0.0 R6_2.6.1
## [73] GenomicAlignments_1.49.2 MatrixGenerics_1.25.0