## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

## ----vignetteSetup, echo=FALSE, message=FALSE, warning = FALSE----------------
## For links
library("BiocStyle")

## Track time spent on making the vignette
startTime <- Sys.time()

## Bib setup
library("RefManageR")

## Write bibliography information
bib <- c(
    R = citation(),
    BiocStyle = citation("BiocStyle"),
    GenomeInfoDb = RefManageR::BibEntry(
        bibtype = "manual",
        key = "GenomeInfoDb",
        author = "Sonali Arora and Martin Morgan and Marc Carlson and H. Pagès",
        title = "GenomeInfoDb: Utilities for manipulating chromosome and other 'seqname' identifiers",
        year = 2017, doi = "10.18129/B9.bioc.GenomeInfoDb"
    ),
    knitr = citation("knitr")[3],
    GenomicState = citation("GenomicState")[1],
    RefManageR = citation("RefManageR")[1],
    rmarkdown = citation("rmarkdown")[1],
    rtracklayer = citation("rtracklayer")[1],
    sessioninfo = citation("sessioninfo"),
    testthat = citation("testthat"),
    GenomicFeatures = citation("GenomicFeatures"),
    bumphunter = citation("bumphunter")[1],
    derfinder = citation("derfinder")[1],
    AnnotationDbi = citation("AnnotationDbi"),
    IRanges = citation("IRanges"),
    org.Hs.eg.db = citation("org.Hs.eg.db"),
    glue = citation("glue"),
    AnnotationHub = citation("AnnotationHub"),
    AnnotationHubData = citation("AnnotationHubData"),
    GenomicRanges = citation("GenomicRanges")
)

## ----'install', eval = FALSE--------------------------------------------------
#  if (!requireNamespace("BiocManager", quietly = TRUE)) {
#        install.packages("BiocManager")
#    }
#  
#  BiocManager::install("GenomicState")
#  
#  ## Check that you have a valid Bioconductor installation
#  BiocManager::valid()

## ----'citation'---------------------------------------------------------------
## Citation info
citation("GenomicState")

## ----setup, message = FALSE, warning = FALSE----------------------------------
library("GenomicState")

## ----'annotation_hub'---------------------------------------------------------
## Query AnnotationHub for the GenomicState object for Gencode v31 on
## hg19 coordinates
hub_query_gs_gencode_v31_hg19 <- GenomicStateHub(
    version = "31",
    genome = "hg19",
    filetype = "GenomicState"
)
hub_query_gs_gencode_v31_hg19


## Check the metadata
mcols(hub_query_gs_gencode_v31_hg19)

## Access the file through AnnotationHub
if (length(hub_query_gs_gencode_v31_hg19) == 1) {
    hub_gs_gencode_v31_hg19 <- hub_query_gs_gencode_v31_hg19[[1]]

    hub_gs_gencode_v31_hg19
}

## ----'build_ex_objects'-------------------------------------------------------
## Load the example TxDb object
## or start from scratch with:
## txdb_v31_hg19_chr21 <- gencode_txdb(version = '31', genome = 'hg19',
##     chrs = 'chr21')
txdb_v31_hg19_chr21 <- AnnotationDbi::loadDb(
    system.file("extdata", "txdb_v31_hg19_chr21.sqlite",
        package = "GenomicState"
    )
)

## Build the GenomicState and annotated genes
genes_v31_hg19_chr21 <- gencode_annotated_genes(txdb_v31_hg19_chr21)
gs_v31_hg19_chr21 <- gencode_genomic_state(txdb_v31_hg19_chr21)

## ----'obtain_annotation_hub'--------------------------------------------------
## Create the AnnotationHub object once and re-use it to speed up things
ah <- AnnotationHub::AnnotationHub()

## Find the TxDb object for hg19 Gencode version 31
hub_query_txdb_gencode_v31_hg19 <- GenomicStateHub(
    version = "31",
    genome = "hg19",
    filetype = "TxDb", ah = ah
)
hub_query_txdb_gencode_v31_hg19

## Now the Annotated Genes for hg19 Gencode v31
hub_query_genes_gencode_v31_hg19 <- GenomicStateHub(
    version = "31",
    genome = "hg19",
    filetype = "AnnotatedGenes", ah = ah
)
hub_query_genes_gencode_v31_hg19

## And finally the GenomicState for hg19 Gencode v31
hub_query_gs_gencode_v31_hg19 <- GenomicStateHub(
    version = "31",
    genome = "hg19",
    filetype = "GenomicState", ah = ah
)
hub_query_gs_gencode_v31_hg19

## If you want to access the files use the double bracket AnnotationHub syntax
## to retrieve the R objects from the web.
if (FALSE) {
    hub_txdb_gencode_v31_hg19 <- hub_query_txdb_gencode_v31_hg19[[1]]
    hub_genes_gencode_v31_hg19 <- hub_query_genes_gencode_v31_hg19[[1]]
    hub_gs_gencode_v31_hg19 <- hub_query_gs_gencode_v31_hg19[[1]]
}

## ----'load_derfinder', messages = FALSE, warnings = FALSE---------------------
## Load external packages
library("derfinder")
library("derfinderPlot")
library("bumphunter")
library("GenomicRanges")

## ----'example_plot_prep'------------------------------------------------------
## Some example regions from derfinder (set the chromosome lengths)
regions <- genomeRegions$regions[1:2]
seqlengths(regions) <- seqlengths(txdb_v31_hg19_chr21)[
    names(seqlengths(regions))
]

## Annotate them
nearestAnnotation <- matchGenes(x = regions, subject = genes_v31_hg19_chr21)
annotatedRegions <- annotateRegions(
    regions = regions,
    genomicState = gs_v31_hg19_chr21$fullGenome, minoverlap = 1
)

## Obtain fullCov object
fullCov <- list("chr21" = genomeDataRaw$coverage)
regionCov <- getRegionCoverage(fullCov = fullCov, regions = regions)

## ----'example_plot'-----------------------------------------------------------
## now make the plot
plotRegionCoverage(
    regions = regions, regionCoverage = regionCov,
    groupInfo = genomeInfo$pop, nearestAnnotation = nearestAnnotation,
    annotatedRegions = annotatedRegions, whichRegions = 1:2,
    txdb = txdb_v31_hg19_chr21, verbose = FALSE
)

## ----'local_metadata'---------------------------------------------------------
## Get the local metadata
meta <- local_metadata()

## Subset to the data of interest, lets say hg19 TxDb for v31
interest <- subset(meta, RDataClass == "TxDb" & Tags == "Gencode:v31:hg19")

## Next you can load the data
if (file.exists(interest$RDataPath)) {
    ## This only works at JHPCE
    eval(parse(text = interest$loadCode))

    ## Explore the loaded object (would be gencode_v31_hg19_txdb in this case)
    gencode_v31_hg19_txdb
}

## ----'build_local', eval = FALSE----------------------------------------------
#  outdir <- "gencode"
#  dir.create(outdir, showWarnings = FALSE)
#  
#  ## Build and save the TxDb object
#  gencode_v23_hg19_txdb <- gencode_txdb("23", "hg19")
#  saveDb(gencode_v23_hg19_txdb,
#      file = file.path(outdir, "gencode_v23_hg19_txdb.sqlite")
#  )
#  
#  ## Build and save the annotateTranscripts output
#  gencode_v23_hg19_annotated_genes <- gencode_annotated_genes(
#      gencode_v23_hg19_txdb
#  )
#  save(gencode_v23_hg19_annotated_genes,
#      file = file.path(outdir, "gencode_v23_hg19_annotated_genes.rda")
#  )
#  
#  ## Build and save the GenomicState
#  gencode_v23_hg19_GenomicState <- gencode_genomic_state(
#      gencode_v23_hg19_txdb
#  )
#  save(gencode_v23_hg19_GenomicState,
#      file = file.path(outdir, "gencode_v23_hg19_GenomicState.rda")
#  )

## ----'build_source'-----------------------------------------------------------
## R commands for building the files:
system.file("scripts", "make-data_gencode_human.R",
    package = "GenomicState"
)
## The above file was created by this one:
system.file("scripts", "generate_make_data_gencode_human.R",
    package = "GenomicState"
)

## ----createVignette, eval=FALSE-----------------------------------------------
#  ## Create the vignette
#  library("rmarkdown")
#  system.time(render("GenomicState.Rmd"))
#  
#  ## Extract the R code
#  library("knitr")
#  knit("GenomicState.Rmd", tangle = TRUE)

## ----reproduce1, echo=FALSE---------------------------------------------------
## Date the vignette was generated
Sys.time()

## ----reproduce2, echo=FALSE---------------------------------------------------
## Processing time in seconds
totalTime <- diff(c(startTime, Sys.time()))
round(totalTime, digits = 3)

## ----reproduce3, echo=FALSE-------------------------------------------------------------------------------------------
## Session info
library("sessioninfo")
options(width = 120)
session_info()

## ----vignetteBiblio, results = 'asis', echo = FALSE, warning = FALSE, message = FALSE---------------------------------
## Print bibliography
PrintBibliography(bib, .opts = list(hyperlink = "to.doc", style = "html"))