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

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

## ----load---------------------------------------------------------------------
library(GDSAnnotator)

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

ann <- seqAnnotList(annot_gds)
ann

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

## ----annot-subset-------------------------------------------------------------
seqAnnotate(snps, annot_gds,
    varnm = c("apc_conservation", "cadd_phred", "linsight"))

## ----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"))

## ----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"))

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

## ----annot-multi, eval = FALSE------------------------------------------------
# annot_gds <- sprintf("favor_chr%s.gds", c(seq_len(22L), "X"))
# seqAnnotate(snps, annot_gds, varnm = c("cadd_phred", "linsight"))

## ----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")

## ----conv-gnomad--------------------------------------------------------------
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)

## ----conv-vep-----------------------------------------------------------------
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)

## ----conv-snpeff--------------------------------------------------------------
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)

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

## ----annot-stat---------------------------------------------------------------
gds_fn <- system.file("extdata", "example_wgs_sites_chr22_snpeff.gds",
    package = "GDSAnnotator")

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

## ----annot-stat-counts--------------------------------------------------------
st$counts$cons_severe

st$counts$impact

## ----annot-gene-table---------------------------------------------------------
head(seqAnnotGeneTable(st))

## ----annot-report-md----------------------------------------------------------
md_fn <- tempfile(fileext = ".md")
seqAnnotReport(st, md_fn, title = "SnpEff annotation of chromosome 22",
    verbose = FALSE)

## ----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")

## ----annot-report-rmd, eval = FALSE-------------------------------------------
# seqAnnotReport(st, "annot_summary.Rmd",
#     title = "SnpEff annotation of chromosome 22")
# 
# rmarkdown::render("annot_summary.Rmd")

## ----annot-report-clean, include = FALSE--------------------------------------
unlink(md_fn, force = TRUE)

## ----annot-stat-parallel, eval = FALSE----------------------------------------
# st <- seqAnnotStat(gds_fn, parallel = 4L)

## ----annot-stat-merge, eval = FALSE-------------------------------------------
# gds_fn <- sprintf("annot_chr%d.gds", seq_len(22L))
# 
# st <- seqAnnotStat(gds_fn)
# seqAnnotReport(st, "wgs_summary.html", title = "Whole-genome annotation")

## ----session-info-------------------------------------------------------------
sessionInfo()