---
title: "GDSAnnotator: Whole-genome variant annotation using GDS files"
author:
- name: Xiuwen Zheng
  affiliation: Genomics Research Center, AbbVie, North Chicago, US
date: "`r BiocStyle::doc_date()`"
package: "`r BiocStyle::pkg_ver('GDSAnnotator')`"
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
    toc_depth: 3
    number_sections: true
vignette: >
  %\VignetteIndexEntry{GDSAnnotator: Whole-genome variant annotation using GDS files}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse=TRUE, comment="##")
# the VCF examples need the suggested package VariantAnnotation
has_va <- requireNamespace("VariantAnnotation", quietly=TRUE)
```

# Introduction

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](https://www.bioconductor.org/packages/SeqArray)) 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](https://useast.ensembl.org/info/docs/tools/vep/index.html)),
[SnpEff](https://pcingola.github.io/SnpEff/),
[gnomAD](https://gnomad.broadinstitute.org/) and
[FAVOR](https://favor.genohub.org) (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.

## From aGDS to annotation GDS files

Storing variant annotations in a GDS file was introduced by
[FAVOR](https://doi.org/10.1093/nar/gkac966) 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:

*   the source of the annotations: not only FAVOR, but also the allele
    frequencies of gnomAD, the consequences predicted by Ensembl VEP or SnpEff,
    and any other resource distributed as a VCF file;
*   the content of the file: an annotation GDS file does not need genotypes, so
    a whole-genome annotation database is kept as one GDS file per chromosome
    and queried by position and alleles, without loading it into memory;
*   the use of the file: the annotations of the variants of a study are written
    into its genotype GDS file block by block (with the FAVOR annotations, this
    produces an aGDS file), and the annotations are summarized and reported
    directly from the GDS file, without the annotator.

## Relation to other Bioconductor packages

`GDSAnnotator` builds on `r Biocpkg("gdsfmt")` and `r Biocpkg("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 `r Biocpkg("GenomicRanges")` `GRanges` objects, `r Biocpkg("S4Vectors")`
`DataFrame` objects, VCF files or `VCF` objects of
`r Biocpkg("VariantAnnotation")` (a suggested package, loaded on demand), and
the annotations are returned as a `DataFrame`.

## Key functions

| 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 |

# Installation

`GDSAnnotator` can be installed from Bioconductor (once released):

```{r install-bioc, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("GDSAnnotator")
```

Load the package:

```{r load}
library(GDSAnnotator)
```

# Inspecting an annotation GDS file

`GDSAnnotator` is built on top of the
[SeqArray](https://www.bioconductor.org/packages/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()`:

```{r annot-list}
# FAVOR essential dataset, chromosome 22 (subset)
annot_gds <- system.file("extdata", "favor_chr22_sub.gds",
    package = "GDSAnnotator")

ann <- seqAnnotList(annot_gds)
ann
```

The returned `DataFrame` has four columns: the annotation `name`, the R data
`type`, the GDS storage `trait`, and a textual `description` (when
available).

# Annotating variants

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 `r Biocpkg("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).

## From a character vector of variant IDs

```{r annot-char}
## 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)
```

A subset of annotations can be selected by passing only the names of
interest in `varnm`:

```{r annot-subset}
seqAnnotate(snps, annot_gds,
    varnm = c("apc_conservation", "cadd_phred", "linsight"))
```

## From a `data.frame`

```{r annot-df}
df <- data.frame(
    chr = c("22", "22", "22"),
    pos = c(10510007L, 10510038L, 10510282L),
    ref = c("T", "T", "G"),
    alt = c("C", "C", "C")
)
df

seqAnnotate(df, annot_gds, varnm = c("apc_conservation", "cadd_phred"))
```

## From a `GRanges` object

`GRanges` queries return all variants that fall within the specified ranges:

```{r annot-granges}
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"))
```

## From a VCF file or a `VCF` object

A `VCF` object of the `r Biocpkg("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:

```{r annot-vcf, eval = has_va}
vcf_fn <- system.file("extdata", "gnomad.genomes.v4.chr22_sub.vcf.gz",
    package = "GDSAnnotator")
vcf <- VariantAnnotation::readVcf(vcf_fn)
dim(vcf)

res <- seqAnnotate(vcf, annot_gds,
    varnm = c("rsid", "cadd_phred", "linsight"), verbose = FALSE)
table(found = !is.na(res$cadd_phred))
head(res[!is.na(res$cadd_phred), ])

# the same, directly from the file name
res2 <- seqAnnotateVCF(vcf_fn, annot_gds,
    varnm = c("rsid", "cadd_phred", "linsight"), verbose = FALSE)
identical(res, res2)
```

## Whole-genome annotation with one file per chromosome

For 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:

```{r annot-multi, eval = FALSE}
annot_gds <- sprintf("favor_chr%s.gds", c(1:22, "X"))
seqAnnotate(snps, annot_gds, varnm = c("cadd_phred", "linsight"))
```

## Writing annotations into a SeqArray GDS file

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.

```{r annot-gds, eval = FALSE}
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")
```

# Converting annotation sources to 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()`.

## gnomAD

```{r conv-gnomad}
vcf_fn <- system.file("extdata", "gnomad.genomes.v4.chr22_sub.vcf.gz",
    package = "GDSAnnotator")

out_fn <- tempfile(fileext = ".gds")
seqToGDS_gnomAD(vcf_fn, out_fn, verbose = FALSE)

unlink(out_fn, force = TRUE)
```

## Ensembl VEP

The VEP importer parses the compound `CSQ` INFO field emitted by Ensembl VEP
into individual annotation fields in the GDS file.

```{r conv-vep}
vcf_fn <- system.file("extdata", "example_wgs_sites_chr22_vep.vcf.gz",
    package = "GDSAnnotator")

out_fn <- tempfile(fileext = ".gds")
seqToGDS_VEP(vcf_fn, out_fn, verbose = FALSE)

unlink(out_fn, force = TRUE)
```

## SnpEff

The SnpEff importer parses the `ANN`, `LOF` and `NMD` INFO fields by
default:

```{r conv-snpeff}
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.

# Summarising annotations

## Counting unique values

`seqValueCounts()` is a convenience wrapper that tabulates unique values of
one or more annotation fields without loading the entire annotation into R:

```{r value-counts}
gds_fn <- system.file("extdata", "gnomad.genomes.v4.chr22_sub.gds",
    package = "GDSAnnotator")

f <- SeqArray::seqOpen(gds_fn)

seqValueCounts(f, "variant_type")
seqValueCounts(f, c("variant_type", "allele_type"))

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.

## Summary report

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:

```{r annot-stat}
gds_fn <- system.file("extdata", "example_wgs_sites_chr22_snpeff.gds",
    package = "GDSAnnotator")

st <- seqAnnotStat(gds_fn, verbose = FALSE)
st
```

Each 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:

```{r annot-stat-counts}
st$counts$cons_severe

st$counts$impact
```

`seqAnnotGeneTable()` returns the per-gene counts, the equivalent of the
`genes.txt` file written by SnpEff:

```{r annot-gene-table}
head(seqAnnotGeneTable(st))
```

`seqAnnotReport()` 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.

```{r annot-report-md}
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:

```{r annot-report-embed, results = "asis"}
txt <- readLines(md_fn)
txt <- txt[seq_len(grep("^## ", txt)[4L] - 1L)]   # the title and 3 sections
cat(sub("^#", "###", txt), sep = "\n")
```

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:

```{r annot-report-rmd, eval = FALSE}
seqAnnotReport(st, "annot_summary.Rmd",
    title = "SnpEff annotation of chromosome 22")

rmarkdown::render("annot_summary.Rmd")
```

```{r annot-report-clean, include = FALSE}
unlink(md_fn, force = TRUE)
```

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:

```{r annot-stat-parallel, eval = FALSE}
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:

```{r annot-stat-merge, eval = FALSE}
gds_fn <- sprintf("annot_chr%d.gds", 1:22)

st <- seqAnnotStat(gds_fn)
seqAnnotReport(st, "wgs_summary.html", title = "Whole-genome annotation")
```

# Session information

```{r session-info}
sessionInfo()
```