diff --git a/data-raw/build-cross-reactive-probes.R b/data-raw/build-cross-reactive-probes.R index 1424ab6..b170014 100644 --- a/data-raw/build-cross-reactive-probes.R +++ b/data-raw/build-cross-reactive-probes.R @@ -29,7 +29,9 @@ ## 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. +## RE-RUN only after that gap is closed, and then re-run +## build-epicv2-reference.R, which appends the fourth list +## (peters2024_epicv2, EPIC v2) to the file this one writes. csv <- "../omicsplayground-methylome-app/components/app_methylome/inst/masking/cross_reactive_probes.csv" stopifnot(file.exists(csv)) diff --git a/data-raw/build-epicv2-reference.R b/data-raw/build-epicv2-reference.R new file mode 100644 index 0000000..0343f09 --- /dev/null +++ b/data-raw/build-epicv2-reference.R @@ -0,0 +1,194 @@ +## Build the two EPIC v2 reference tables playbase.epigenetics reads: +## +## inst/extdata/epicv2-replicates.rds which replicate probe to keep per CpG +## inst/extdata/cross-reactive-probes.rds + the EPIC v2 cross-hybridising +## list, source "peters2024_epicv2" +## +## Run from the package root, after build-cross-reactive-probes.R when that one +## is re-run (this script replaces only its own source's rows, so either order +## of re-runs converges). +## +## ---- epicv2-replicates.rds ------------------------------------------------- +## +## WHY EPIC v2 names every probe _ (cg00000029_TC21) and +## measures ~5,200 ids with 2-10 replicate probes. Clocks, cell +## references, the EWAS Catalog and the masks are all keyed by the bare +## id, so playbase.epigenetics::collapse_epicv2_replicates() keeps one +## replicate per id. Never the mean: replicates differ in design type +## and in how well they track true methylation, so their average is a +## probe nobody characterised. +## +## SHAPE one row per probe whose bare id is replicated in the manifest +## (11,616 probes over 5,222 ids, ordered by id then rank): +## probe_id full IlmnID (factor) +## cpg bare id (factor) +## recommended keep this one (logical, exactly one per cpg) +## evidence Peters 2024 verdict (factor, NA when none) +## attr "manifest" = "20a1": the Illumina manifest release, as the +## IlluminaHumanMethylationEPICv2anno.20a1.hg38 package names it. +## playbase.epigenetics refuses a table whose manifest differs from its +## annotation package, so a manifest bump needs a paired release here. +## attr "source" names both inputs. +## +## WHICH Illumina's manifest does NOT designate a preferred replicate: its +## only replicate field is Rep_Num, a synthesis counter. The one +## published per-probe evaluation is Peters et al. 2024 (below), which +## compared each replicate against matched EPIC v1 and WGBS and labels +## every replicated probe in "Rep_results_by_NAME". Ranked, first wins: +## 1 "Superior probe" / "Superior by WGBS" best on both counts +## 2 "Best sensitivity" tracks methylation change +## 3 "Best precision" +## 4 no verdict, "Insufficient evidence", or a verdict that only the +## probe-set MEAN is best ("... group mean ...") - averaging is +## ruled out above, so these fall through to the tie-break +## 5 "Inferior probe" / "Inferior by WGBS" +## Ties go to a probe the annotation package can place (two sets mix +## a chr0 replicate with a mapped one), then the lowest replicate +## number, then the probe id. This is DMRcate's filter.strategy = +## "sensitivity" (same author), made deterministic: DMRcate breaks +## ties with sample(). +## +## ---- cross-reactive-probes.rds, source peters2024_epicv2 -------------------- +## +## Peters 2024 BLAT-aligned every EPIC v2 probe; "CH_BLAT" = Y marks the 30,627 +## in-silico cross-hybridising probes (47+ matching bases off target), and the +## paper adds the 6,889 probes Illumina could not place ("CHR" = chr0, 97% +## multimappers) to reach its "total of 37,346". Both are taken, matching the +## in-silico definition of the three lists already in the file. Ids are bare, +## and a replicated id is flagged by the replicate that survives the collapse +## above - the probe whose betas a dataset actually holds. SNP-affected probes +## are not added here: the Methylome app masks those from the manifest's own +## SNP columns per array. +## +## SOURCE Peters TJ et al. (2024) Characterisation and reproducibility of the +## HumanMethylationEPIC v2.0 BeadChip for DNA methylation profiling. +## BMC Genomics 25:251. doi:10.1186/s12864-024-10027-5 +## Additional file 4 (augmented EPIC v2 manifest, 1.2 GB CSV). +## Bioconductor's EPICv2manifest package carries the same table, but +## only from Bioconductor 3.19 and via AnnotationHub; the published +## supplement is the citable artefact, pinned by checksum below. +## Probe universe: the installed IlluminaHumanMethylationEPICv2manifest; +## placement: IlluminaHumanMethylationEPICv2anno.20a1.hg38 (both +## Bioconductor 3.19, both install on R >= 4.2 / minfi 1.48). +## RE-RUN when Illumina ships a new EPIC v2 manifest (new annotation package +## name) or a revised replicate evaluation is published. + +suppressPackageStartupMessages({ + stopifnot(requireNamespace("data.table"), requireNamespace("minfi")) +}) +ANNO <- "IlluminaHumanMethylationEPICv2anno.20a1.hg38" +stopifnot(requireNamespace(ANNO)) +MANIFEST <- strsplit(ANNO, ".", fixed = TRUE)[[1]][2] # "20a1" +SUFFIX <- "_[TB][CO][12][0-9]+$" + +PETERS_URL <- paste0( + "https://static-content.springer.com/esm/art%3A10.1186%2Fs12864-024-10027-5/", + "MediaObjects/12864_2024_10027_MOESM4_ESM.csv" +) +PETERS_MD5 <- "259eb024cf0ab064eb348c376c40dbd2" + +options(timeout = 3600) # 1.2 GB +f <- file.path(tempdir(), "peters2024_epicv2_manifest.csv") +if (!file.exists(f)) { + message("downloading Peters 2024 Additional file 4 ...") + utils::download.file(PETERS_URL, f, quiet = TRUE, mode = "wb") +} +stopifnot(unname(tools::md5sum(f)) == PETERS_MD5) +peters <- data.table::fread( + f, select = c("IlmnID", "CHR", "CH_BLAT", "Rep_results_by_NAME") +) +peters <- as.data.frame(peters) +rownames(peters) <- peters$IlmnID + +## ---- replicate table -------------------------------------------------------- + +## The universe is every probe minfi can return a beta for: the manifest +## package's type I and II probes (not the rs genotyping probes). +MANIFEST_PKG <- "IlluminaHumanMethylationEPICv2manifest" +stopifnot(requireNamespace(MANIFEST_PKG)) +mf <- getExportedValue(MANIFEST_PKG, MANIFEST_PKG) +ids <- c(minfi::getProbeInfo(mf, type = "I")$Name, minfi::getProbeInfo(mf, type = "II")$Name) +stopifnot(all(grepl(SUFFIX, ids)), !anyDuplicated(ids)) +bare <- sub(SUFFIX, "", ids) +is_rep <- bare %in% bare[duplicated(bare)] + +## Probes the annotation package can place. It leaves out Illumina's chr0 +## (unmapped) and chrM probes, so their betas never reach an annotated matrix. +library(ANNO, character.only = TRUE) # getAnnotation() reads the search path +placed <- ids %in% rownames(minfi::getAnnotation(get(ANNO))) + +## Every replicated probe is one Peters evaluated. +stopifnot(all(ids[is_rep] %in% peters$IlmnID)) + +evidence <- peters[ids[is_rep], "Rep_results_by_NAME"] +evidence[!nzchar(evidence)] <- NA_character_ +tier <- ifelse( + grepl("^Superior (probe|by WGBS)$", evidence), 1L, + ifelse(evidence %in% "Best sensitivity", 2L, + ifelse(evidence %in% "Best precision", 3L, + ifelse(grepl("^Inferior", evidence), 5L, 4L) + ) + ) +) + +reps <- data.frame( + probe_id = ids[is_rep], cpg = bare[is_rep], evidence = evidence, + tier = tier, unplaced = !placed[is_rep], + rep_num = as.integer(sub("^.*_[TB][CO][12]", "", ids[is_rep])), + stringsAsFactors = FALSE +) +reps <- reps[order(reps$cpg, reps$tier, reps$unplaced, reps$rep_num, reps$probe_id), ] +reps$recommended <- !duplicated(reps$cpg) + +stopifnot( + all(table(reps$cpg[reps$recommended]) == 1), + setequal(reps$cpg[reps$recommended], unique(reps$cpg)) +) + +EPICV2_REPLICATES <- data.frame( + probe_id = factor(reps$probe_id), + cpg = factor(reps$cpg), + recommended = reps$recommended, + evidence = factor(reps$evidence), + row.names = NULL +) +attr(EPICV2_REPLICATES, "manifest") <- MANIFEST +attr(EPICV2_REPLICATES, "source") <- paste0( + MANIFEST_PKG, " ", utils::packageVersion(MANIFEST_PKG), " and ", + ANNO, " ", utils::packageVersion(ANNO), " (Illumina EPIC v2.0 manifest ", + MANIFEST, "); replicate choice from Peters et al. 2024, BMC Genomics 25:251, ", + "Additional file 4 (md5 ", PETERS_MD5, "), Rep_results_by_NAME; built ", + Sys.Date() +) + +message( + nrow(EPICV2_REPLICATES), " replicate probes over ", + nlevels(EPICV2_REPLICATES$cpg), " ids; recommended by verdict:" +) +print(table(reps$evidence[reps$recommended], useNA = "ifany")) + +saveRDS(EPICV2_REPLICATES, "inst/extdata/epicv2-replicates.rds", compress = "xz") + +## ---- cross-reactive list ---------------------------------------------------- + +survives <- !is_rep +survives[match(reps$probe_id, ids)] <- reps$recommended +flagged <- peters$IlmnID[peters$CH_BLAT %in% "Y" | peters$CHR %in% "chr0"] +xr <- sort(unique(bare[survives & ids %in% flagged])) +message(length(flagged), " flagged probes -> ", length(xr), " bare ids after the collapse") + +old <- readRDS("inst/extdata/cross-reactive-probes.rds") +old <- old[old$source != "peters2024_epicv2", , drop = FALSE] +CROSS_REACTIVE_PROBES <- data.frame( + probe = factor(c(as.character(old$probe), xr)), + source = factor(c(as.character(old$source), rep("peters2024_epicv2", length(xr)))), + 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") diff --git a/inst/extdata/cross-reactive-probes.rds b/inst/extdata/cross-reactive-probes.rds index 08f9c19..0b0286f 100644 Binary files a/inst/extdata/cross-reactive-probes.rds and b/inst/extdata/cross-reactive-probes.rds differ diff --git a/inst/extdata/epicv2-replicates.rds b/inst/extdata/epicv2-replicates.rds new file mode 100644 index 0000000..7bba6ae Binary files /dev/null and b/inst/extdata/epicv2-replicates.rds differ