## ----echo = FALSE, message = FALSE--------------------------------------------
library(knitr)
knitr::opts_chunk$set(
    error = FALSE,
    tidy  = FALSE,
    message = FALSE,
    warning = FALSE,
    fig.align = "center")

## -----------------------------------------------------------------------------
library(rGREAT)

set.seed(123)
gr = randomRegions(nr = 1000, genome = "hg19")

## -----------------------------------------------------------------------------
res = great(gr, "MSigDB:H", "txdb:hg19")
res

## ----eval = FALSE-------------------------------------------------------------
#  great(gr, "GO:BP", "hg19")
#  great(gr, "GO:BP", "TxDb.Hsapiens.UCSC.hg19.knownGene")
#  great(gr, "GO:BP", "RefSeq:hg19")
#  great(gr, "GO:BP", "GREAT:hg19")
#  great(gr, "GO:BP", "Gencode_v19")

## -----------------------------------------------------------------------------
gs = read_gmt(url("http://dsigdb.tanlab.org/Downloads/D2_LINCS.gmt"), 
    from = "SYMBOL", to = "ENTREZ", orgdb = "org.Hs.eg.db")
gs[1:2]
great(gr, gs, "hg19")

## -----------------------------------------------------------------------------
df = read.table(url("https://jokergoo.github.io/rGREAT_suppl/data/GREATv4.genes.hg19.tsv"))
# note there must be a 'gene_id' column
tss = GRanges(seqnames = df[, 2], ranges = IRanges(df[, 3], df[, 3]), 
    strand = df[, 4], gene_id = df[, 5])
head(tss)

## -----------------------------------------------------------------------------
et = extendTSS(tss, genome = "hg19", gene_id_type = "SYMBOL")
great(gr, "msigdb:h", extended_tss = et)

## ----eval = FALSE-------------------------------------------------------------
#  great(gr, gs, tss)  # tss has a `gene_id` meta column

## -----------------------------------------------------------------------------
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
gene = genes(TxDb.Hsapiens.UCSC.hg19.knownGene)
gene = gene[seqnames(gene) %in% paste0("chr", c(1:22, "X", "Y"))]
head(gene)
gl = seqlengths(gene)[paste0("chr", c(1:22, "X", "Y"))]  # restrict to normal chromosomes
head(gl)

## -----------------------------------------------------------------------------
et = extendTSS(gene, gl)
head(et)

## -----------------------------------------------------------------------------
great(gr, gene_sets = gs, extended_tss = et)

## ----eval = FALSE-------------------------------------------------------------
#  gap = getGapFromUCSC("hg19", paste0("chr", c(1:22, "X", "Y")))
#  great(gr, "MSigDB:H", "hg19", exclude = gap)

## ----eval = FALSE-------------------------------------------------------------
#  great(gr, "GO:BP", background = paste0("chr", 1:22))
#  great(gr, "GO:BP", exclude = c("chrX", "chrY"))

## -----------------------------------------------------------------------------
tb = getEnrichmentTable(res)
head(tb)

## ----fig.width = 6, fig.height = 6--------------------------------------------
plotVolcano(res)

## ----fig.width = 10, fig.height = 10/3, fig.align = 'center'------------------
plotRegionGeneAssociations(res)
getRegionGeneAssociations(res)

## ----fig.width = 10, fig.height = 10/3, fig.align = 'center'------------------
plotRegionGeneAssociations(res, term_id = "HALLMARK_APOPTOSIS")
getRegionGeneAssociations(res, term_id = "HALLMARK_APOPTOSIS")

## ----eval = FALSE-------------------------------------------------------------
#  shinyReport(res)

## -----------------------------------------------------------------------------
sessionInfo()