## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
ext <- system.file("extdata", "yang2024", package = "ctdR")
read_tsv <- function(f) utils::read.delim(file.path(ext, f),
    stringsAsFactors = FALSE, check.names = FALSE)
genes <- read_tsv("yang2024_genes.tsv")
rep <- read_tsv("yang2024_replication.tsv")
top <- read_tsv("yang2024_discovery_top30.tsv")
pos <- read_tsv("yang2024_discovery_yang_chemicals.tsv")
summ <- read_tsv("yang2024_discovery_summary.tsv")
yang_vs_ctdR <- read_tsv("yang2024_comparison.tsv")
fmt <- function(p) ifelse(is.na(p), "", formatC(p, format = "g", digits = 3))

## ----download, eval = FALSE---------------------------------------------------
# library(ctdR)
# ctd_file <- "CTD_chem_gene_ixns_2023-05-23.csv.gz"
# utils::download.file(paste0("https://web.archive.org/web/20230607020439id_/",
#     "https://ctdbase.org/reports/CTD_chem_gene_ixns.csv.gz"), ctd_file,
#     mode = "wb")

## ----import, eval = FALSE-----------------------------------------------------
# import_CTD(ctd_file, organism = "all")
# ctd_provenance()
# #> CTD provenance
# #>   Report created: Tue May 23 14:45:41 EDT 2023
# #>   Retained:       14,572 chemicals, 1,523,407 chemical-gene pairs
# #>   Organism:       all species

## ----mapping------------------------------------------------------------------
table(rep$Mapping, useNA = "ifany")

## ----same-sets----------------------------------------------------------------
m <- rep[!is.na(rep$ChemicalID) & nzchar(rep$ChemicalID), ]
c(chemicals = nrow(m),
  N_identical = sum(m$N_ctdR == m$N_paper, na.rm = TRUE),
  N_within_2pct = sum(abs(m$N_ctdR / m$N_paper - 1) <= 0.02, na.rm = TRUE),
  k_identical = sum(m$k_ctdR == m$k_paper, na.rm = TRUE))

## ----same-sets-plot, fig.width = 7, fig.height = 3.6--------------------------
op <- par(mfrow = c(1, 2), mar = c(4.2, 4.2, 2, 1))
plot(m$N_paper, m$N_ctdR, log = "xy", pch = 19, col = "#1f77b4",
    xlab = "N, Yang et al.", ylab = "N, ctdR", main = "Set size")
abline(0, 1, lty = 2)
plot(m$k_paper, m$k_ctdR, pch = 19, col = "#1f77b4",
    xlab = "k, Yang et al.", ylab = "k, ctdR", main = "Dyslexia genes in set")
abline(0, 1, lty = 2)
par(op)

## ----same-test----------------------------------------------------------------
ok <- !is.na(m$PValue_paper) & m$k_paper > 0 &
    abs(m$N_ctdR / m$N_paper - 1) <= 0.02 & m$k_ctdR == m$k_paper
cmp <- m[ok, ]
c(chemicals = nrow(cmp),
  median_abs_log10_ratio = round(median(abs(log10(
      cmp$PValue_sameTest_ctdR / cmp$PValue_paper))), 3),
  same_call_at_0.05 = sum((cmp$PValue_sameTest_ctdR < 0.05) ==
      (cmp$PValue_paper < 0.05)))

## ----same-test-table----------------------------------------------------------
show <- cmp[order(cmp$PValue_paper), c("Class", "Chemical", "N_paper",
    "N_ctdR", "k_paper", "PValue_paper", "PValueFDR_paper",
    "PValue_sameTest_ctdR")]
show[, 6:8] <- lapply(show[, 6:8], fmt)
knitr::kable(show, row.names = FALSE, col.names = c("Class", "Chemical",
    "N paper", "N ctdR", "k", "p paper", "FDR paper", "p ctdR sets"))

## ----discovery-run, eval = FALSE----------------------------------------------
# entrez <- AnnotationDbi::mapIds(org.Hs.eg.db::org.Hs.eg.db,
#     genes$GeneSymbol, "ENTREZID", "SYMBOL")
# coding <- unique(AnnotationDbi::select(org.Hs.eg.db::org.Hs.eg.db,
#     keys = "protein-coding", keytype = "GENETYPE", columns = "SYMBOL")$SYMBOL)
# discovery <- enrichment_CTD(
#     data.frame(EntrezID = unname(entrez[!is.na(entrez)]), pvalue = 0.01),
#     method = "ORA", universe = coding, minGSSize = 2)

## ----discovery-summary--------------------------------------------------------
knitr::kable(summ, col.names = c("Analysis", "Chemicals tested",
    "Adjusted p < 0.05", "Unadjusted p < 0.05"))

## ----discovery-top------------------------------------------------------------
t30 <- top
t30$PValue <- fmt(t30$PValue)
t30$PValueAdjusted <- fmt(t30$PValueAdjusted)
knitr::kable(t30[, c("Rank", "ChemicalName", "Count", "GeneSetSize",
    "FoldEnrichment", "PValue", "PValueAdjusted", "InYangList")],
    col.names = c("Rank", "Chemical", "k", "N", "Fold", "p",
        "Adjusted p", "In Yang's list"))

## ----discovery-yang-----------------------------------------------------------
p35 <- pos[pos$Analysis == "main", c("Class", "Chemical", "Rank", "Count",
    "GeneSetSize", "PValue", "PValueAdjusted")]
p35$PValue <- fmt(p35$PValue)
p35$PValueAdjusted <- fmt(p35$PValueAdjusted)
knitr::kable(p35, row.names = FALSE, col.names = c("Class", "Chemical",
    "Rank", "k", "N", "p", "Adjusted p"))

## ----discovery-yang-count-----------------------------------------------------
sapply(split(pos, pos$Analysis), function(d) c(
    tested = sum(!is.na(d$Rank)),
    adjusted_below_0.05 = sum(d$PValueAdjusted < 0.05, na.rm = TRUE),
    unadjusted_below_0.05 = sum(d$PValue < 0.05, na.rm = TRUE)))

## ----comparison-count---------------------------------------------------------
outcome_levels <- c("significant in both", "significant only in Yang et al.",
    "significant only in ctdR", "significant in neither", "not tested by ctdR")
cbind(
    unadjusted = table(factor(yang_vs_ctdR$Outcome_unadjusted, outcome_levels)),
    adjusted = table(factor(yang_vs_ctdR$Outcome_adjusted, outcome_levels)))

## ----comparison-table---------------------------------------------------------
agree <- c("significant in both", "significant in neither")
disagree <- yang_vs_ctdR[!(yang_vs_ctdR$Outcome_unadjusted %in% agree &
    yang_vs_ctdR$Outcome_adjusted %in% agree), ]
disagree <- disagree[, c("Chemical", "k_paper", "N_paper", "PValue_paper",
    "PValue_discovery", "PValueFDR_paper", "PValueAdjusted_discovery",
    "Outcome_unadjusted", "Outcome_adjusted")]
disagree[, 4:7] <- lapply(disagree[, 4:7], fmt)
knitr::kable(disagree, row.names = FALSE, col.names = c("Chemical", "k", "N",
    "p paper", "p ctdR", "FDR paper", "Adjusted p ctdR", "Unadjusted",
    "Adjusted"))

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