Methylation reference data: EWAS Catalog and cross-reactive probes

data-raw/build-cross-reactive-probes.R

@@ -0,0 +1,52 @@
+## Build CROSS_REACTIVE_PROBES: Illumina methylation probes that map to more
+## than one genome location, so their signal is not the CpG it is attributed
+## to. Excluded from an EWAS by the "Cross-reactive" mask in the Methylome
+## Profiler; see mp_xreactive_probes() in components/app_methylome.
+##
+## Why extdata and not data/: same reason as the EWAS catalog next door - a
+## LazyData object stays in the namespace for the life of the session once
+## touched, and this is read once per fit. As an .rds the reader owns it.
+##
+## Why here and not in the app repo: it is published reference data, not app
+## code. It lived at components/app_methylome/inst/masking/ as a 2.7 MB CSV
+## behind a .gitignore negation, and the loader tried a list of candidate
+## paths because the app's working directory varies - which resolved to
+## nothing once, silently, leaving 53,498 probes tested while the UI reported
+## a mask had been applied. playdata::get_file() is system.file(mustWork =
+## TRUE), so that failure mode is gone.
+##
+## SOURCE   Three published blacklists, carried in the `source` column:
+##            chen2013_450k       29,233 rows  (450K)
+##            mccartney2016_epic  44,210 rows  (EPIC)
+##            zhou2017_mask       30,681 rows  (non-SNP general mask)
+##          104,124 rows, 53,498 unique probes - a probe flagged by two lists
+##          appears twice, and the column is kept so provenance survives.
+##
+## PROVENANCE GAP: the download URLs were never recorded. The CSV arrived
+##          pre-assembled in omicsplayground 2825f3bc0 with only these three
+##          labels, and no build script. The labels identify the papers, but
+##          the exact files and accession dates are not known, so this script
+##          converts what we have rather than rebuilding from source. Before
+##          the next refresh, resolve each label to its published supplement
+##          and rewrite this as a real download.
+## RE-RUN   only after that gap is closed.
+
+csv <- "../omicsplayground-methylome-app/components/app_methylome/inst/masking/cross_reactive_probes.csv"
+stopifnot(file.exists(csv))
+
+d <- utils::read.csv(csv, stringsAsFactors = FALSE)
+stopifnot(identical(colnames(d), c("probe", "source")))
+
+CROSS_REACTIVE_PROBES <- data.frame(
+  probe  = factor(d$probe),
+  source = factor(d$source),
+  stringsAsFactors = FALSE
+)
+attr(CROSS_REACTIVE_PROBES, "sources") <- levels(CROSS_REACTIVE_PROBES$source)
+
+message(nrow(CROSS_REACTIVE_PROBES), " rows, ",
+        nlevels(CROSS_REACTIVE_PROBES$probe), " unique probes, from ",
+        nlevels(CROSS_REACTIVE_PROBES$source), " lists")
+
+saveRDS(CROSS_REACTIVE_PROBES, "inst/extdata/cross-reactive-probes.rds",
+        compress = "xz")

data-raw/build-ewas-catalog.R

@@ -0,0 +1,70 @@
+## Build EWAS_CATALOG: which traits each CpG has already been associated with,
+## and in how many study analyses - the "is this hit known or new?" annotation
+## on an EWAS hit list.
+##
+## Why bundled: the app used to ask ewascatalog.org/api/ at runtime, one HTTP
+## GET per CpG. That needs outbound network from the container, is slow enough
+## to need a button, and pins us to whatever is live. The catalog publishes the
+## whole thing as a flat file, so we ship a version we control.
+##
+## Why extdata and not data/: a LazyData object stays in the package namespace
+## for the life of the R session once touched - 129 MB resident in every Shiny
+## worker, for a table read once per hit list. As an .rds the reader owns it and
+## the GC takes it back; see playbase.epigenetics::ewas_catalog_lookup().
+##
+## SOURCE   https://www.ewascatalog.org/download/
+## VERSION  the results/studies pair Last-Modified 2026-04-09 (the catalog
+##          publishes no version string; the mtime is the release marker)
+## RE-RUN   when the catalog publishes a new release (roughly yearly)
+##
+## Reduction: results is 174 MB gzipped / 8.0M rows of per-association summary
+## statistics. Nothing but the CpG and the study it came from is ever rendered,
+## so beta/se/p/position/gene are dropped, the study id is resolved to its trait
+## through the studies metadata file, and the pairs are counted. 8.0M rows of 13
+## columns collapse to 6.3M rows of 3.
+
+base <- "https://www.ewascatalog.org/static//docs/"
+options(timeout = 3600) # the results file is 174 MB; the 60s default never finishes
+
+fetch <- function(f) {
+  p <- file.path(tempdir(), f)
+  if (!file.exists(p)) {
+    message("downloading ", f, " ...")
+    utils::download.file(paste0(base, f), p, quiet = TRUE, mode = "wb")
+  }
+  p
+}
+
+## study_id -> trait. study_id is one analysis within a paper ("22232023_smoking
+## _current_vs_never_smoking"), which is the unit the old API counted, so
+## counting distinct study_ids reproduces the numbers the app used to show.
+studies <- data.table::fread(fetch("ewascatalog-studies.txt.gz"),
+                             select = c("study_id", "trait"), sep = "\t",
+                             quote = "", colClasses = "character")
+
+## Column 1 is cpg and column 6 is study_id; the file repeats the cpg column at
+## position 7, so select by index rather than by name.
+res <- data.table::fread(fetch("ewascatalog-results.txt.gz"),
+                         select = c(1L, 6L), col.names = c("cpg", "study_id"),
+                         sep = "\t", quote = "", colClasses = "character")
+
+res <- unique(res) # one CpG can appear twice in an analysis
+res[studies, trait := i.trait, on = "study_id"]
+res <- res[!is.na(trait) & nzchar(trait)]
+
+EWAS_CATALOG <- res[, .(n = data.table::uniqueN(study_id)), by = .(cpg, trait)]
+data.table::setorder(EWAS_CATALOG, cpg, -n, trait)
+EWAS_CATALOG <- data.frame(
+  cpg = factor(EWAS_CATALOG$cpg),
+  trait = factor(EWAS_CATALOG$trait),
+  n = as.integer(EWAS_CATALOG$n),
+  stringsAsFactors = FALSE
+)
+attr(EWAS_CATALOG, "source") <- "https://www.ewascatalog.org/download/"
+attr(EWAS_CATALOG, "version") <- "2026-04-09"
+
+message(nrow(EWAS_CATALOG), " CpG-trait pairs, ",
+        nlevels(EWAS_CATALOG$cpg), " CpGs, ",
+        nlevels(EWAS_CATALOG$trait), " traits")
+
+saveRDS(EWAS_CATALOG, "inst/extdata/ewas-catalog.rds", compress = "xz")