## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)
# Build this vignette against a throwaway cache. Without it, knitting
# would overwrite whatever CTD data the reader has already imported
# with the ten-chemical sample used for the examples below.
options(ctdR.cache = tempfile())

## ----import-real, eval=FALSE--------------------------------------------------
# library(ctdR)
# import_CTD("~/Downloads/CTD_chem_gene_ixns.csv")

## ----import-sample------------------------------------------------------------
library(ctdR)
sample_file <- system.file(
    "extdata", "CTD_chem_gene_ixns_sample.csv",
    package = "ctdR"
)
import_CTD(sample_file)

## ----provenance---------------------------------------------------------------
# With no argument: what is in the cache right now.
ctd_provenance()

## ----prepare-genes------------------------------------------------------------
genes <- data.frame(
    EntrezID = c(
        "7124", "3569", "7157", "672", "1956",
        "4609", "3845", "207", "5290", "3553"
    ),
    pvalue = c(
        0.001, 0.003, 0.005, 0.008, 0.01,
        0.02, 0.03, 0.04, 0.05, 0.06
    )
)

## ----ora----------------------------------------------------------------------
ora_results <- enrichment_CTD(genes, method = "ORA")
head(ora_results)

## ----gsea---------------------------------------------------------------------
# Quick-start: no signed statistic available in this
# toy gene list, so the pvalue fallback is used.
gsea_results <- enrichment_CTD(genes, method = "GSEA")
head(gsea_results)

## ----load-real----------------------------------------------------------------
library(SummarizedExperiment)

se <- readRDS(system.file(
    "extdata", "GSE311566_subset.rds", package = "ctdR"
))
se
table(se$group)

## ----de-real------------------------------------------------------------------
expr <- assay(se)
grp  <- se$group

de <- t(apply(expr, 1, function(y) {
    tt <- stats::t.test(y ~ grp)
    c(log2FC = unname(diff(tt$estimate)),
      pvalue = tt$p.value)
}))
de <- as.data.frame(de)
de$padj <- stats::p.adjust(de$pvalue, method = "BH")
de$EntrezID <- rownames(de)
de <- de[, c("EntrezID", "log2FC", "pvalue", "padj")]
head(de[order(de$padj), ])
c(padj_lt_05 = sum(de$padj < 0.05),
  pvalue_lt_05 = sum(de$pvalue < 0.05))

## ----gsea-real----------------------------------------------------------------
gsea_real <- enrichment_CTD(
    data.frame(EntrezID = de$EntrezID,
               pvalue   = de$pvalue,
               stat     = de$log2FC),
    method = "GSEA"
)
head(gsea_real[, c("ChemicalID", "ChemicalName",
                   "PValue", "PValueAdjusted",
                   "NormalizedEnrichmentScore")])

## ----ora-is-subset------------------------------------------------------------
## Set is_subset <- FALSE to run ORA on the full CTD download.
## Requires CTD_chem_gene_ixns.csv.gz from ctdbase.org, decompressed locally.
is_subset <- TRUE

## ----ora-real-----------------------------------------------------------------
sig <- de[de$pvalue < 0.05, c("EntrezID", "pvalue")]

if (is_subset) {
    ora_real <- enrichment_CTD(sig, method = "ORA")
    if (nrow(ora_real)) {
        head(ora_real[, c("ChemicalID", "ChemicalName",
                          "PValue", "PValueAdjusted", "Count")])
    } else {
        message(
            "ORA returned 0 results on the toy CTD sample (expected).\n",
            "Set is_subset <- FALSE and provide CTD_chem_gene_ixns.csv",
            " for a powered analysis."
        )
    }
} else {
    ## Full CTD path -----------------------------------------------------------
    ## 1. Download https://ctdbase.org/reports/CTD_chem_gene_ixns.csv.gz
    ## 2. gunzip CTD_chem_gene_ixns.csv.gz
    ## 3. Set ctd_path to the decompressed file location
    ctd_path <- "~/Downloads/CTD_chem_gene_ixns.csv"  # adjust to your path
    import_CTD(ctd_path)
    ora_full <- enrichment_CTD(sig, method = "ORA")
    if (nrow(ora_full)) {
        head(ora_full[order(ora_full$PValue),
                      c("ChemicalID", "ChemicalName",
                        "PValue", "PValueAdjusted",
                        "FoldEnrichment", "Count")])
    } else {
        message("No chemicals enriched — check CTD file path and gene list.")
    }
}

## ----camera-------------------------------------------------------------------
design <- model.matrix(~ se$group)

camera_results <- enrichment_CTD(
    se,
    method   = "CAMERA",
    design   = design,
    contrast = 2  # column 2 of design = Dex vs DMSO
)
head(camera_results)

## ----gsva---------------------------------------------------------------------
gsva_scores <- enrichment_CTD(se, method = "GSVA")
dim(gsva_scores)  # chemicals x samples
head(assay(gsva_scores), 3)

# the treatment groups are still there, next to the scores
colData(gsva_scores)

## ----gsva-tune----------------------------------------------------------------
gsva_strict <- enrichment_CTD(
    se, method = "GSVA",
    minSize = 5, maxSize = 500
)
nrow(gsva_strict)

## ----plot-bar, fig.width=7, fig.height=5--------------------------------------
# ORA bar plot of top enriched chemicals
plot_CTD(ora_results, type = "bar")

## ----plot-dot, fig.width=7, fig.height=5--------------------------------------
# ORA dot plot (size = gene count, color = adjusted p-value)
plot_CTD(ora_results, type = "dot", n = 10)

## ----plot-camera, fig.width=7, fig.height=5-----------------------------------
# CAMERA bar plot: x-axis = -log10(padj), fill = Direction
plot_CTD(camera_results, type = "bar")

## ----plot-gsva, fig.width=7, fig.height=5-------------------------------------
# GSVA heatmap: top-variance chemicals across samples
plot_CTD(gsva_scores)

## ----padj-methods-------------------------------------------------------------
# Bonferroni correction (more conservative)
ora_bonf <- enrichment_CTD(genes, method = "ORA",
    pAdjustMethod = "bonferroni")

# No adjustment
ora_raw <- enrichment_CTD(genes, method = "ORA",
    pAdjustMethod = "none")

# Compare the adjusted p-value of the top hit
data.frame(
    method = c("BH (default)", "bonferroni", "none"),
    top_padj = c(min(ora_results$padj),
        min(ora_bonf$padj), min(ora_raw$padj))
)

## ----size-universe------------------------------------------------------------
# ORA: restrict the background to expressed genes. Both size thresholds
# are left at their defaults; maxGSSize is shown here only to make the
# point that capping it is a deliberate choice, not a starting position.
ora_expressed <- enrichment_CTD(
    genes, method = "ORA",
    universe   = c(genes$EntrezID, "7422", "836")  # expressed genes
)

# GSEA: filter gene sets by size
gsea_filtered <- enrichment_CTD(
    genes, method = "GSEA",
    minSize = 3,
    maxSize = 300
)

## ----ctdfile, eval=FALSE------------------------------------------------------
# library(BiocIO)
# ctd <- import(CTDFile("~/Downloads/CTD_chem_gene_ixns.csv"))
# head(ctd)

## ----asgenesets---------------------------------------------------------------
gene_sets <- as_genesets_CTD("entrez")
length(gene_sets)
gene_sets[[1]]

## ----sbea, eval=FALSE---------------------------------------------------------
# # Illustrative: route the CTD gene sets through EnrichmentBrowser's
# # set-based engine. sbea() covers all four of ctdR's methods, GSVA
# # included, plus eight more (see EnrichmentBrowser::sbeaMethods()).
# library(EnrichmentBrowser)
# res <- sbea(method = "camera", se = my_se, gs = as_genesets_CTD("entrez"))
# gsRanking(res)

## ----reimport, eval=FALSE-----------------------------------------------------
# import_CTD("~/Downloads/CTD_chem_gene_ixns.csv")

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