---
title: "tewas: unit-level multiple testing for EWAS, and a frame layer for behavioral designs"
author: "Tsachi Ein-Dor"
date: "`r Sys.Date()`"
package: tewas
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
vignette: >
  %\VignetteIndexEntry{tewas: unit-level multiple testing for EWAS, and a frame layer for behavioral designs}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5)
```

# What tewas is for

An epigenome-wide association study tests hundreds of thousands of CpG probes one at a time and corrects them one at a time, although association signal on methylation arrays is spatially structured and the probes are not exchangeable. tewas tests *units* instead: runs of adjacent probes no more than 500 or 2,000 bp apart, each read by four complementary statistics (a dense sum, a coherence-gated flip, the Simes combination, a kernel matched filter) whose null distributions come from one shared residual-permutation stream. The arm p-values are combined by the Cauchy combination, units are selected by weighted Benjamini–Hochberg with cross-fitted structural weights, and probes inside a selected unit are tested at a local level. On a real 450k or EPIC array the scan also carries an anchor panel: probes that measure whether an external knowledge prior is licensed on this sample (the lens) and whether the scan's null is calibrated (the inflation guard, read on an alive null ruler). The prior is tempered toward zero when it is not licensed, and nothing is ever corrected: every calibration quantity is reported beside the result.

Above the engine sits a *frame layer* for the designs a behavioral epigeneticist asks for: moderation by adversity, curvature, a threshold, acceleration, methylation as the outcome of an exposure rather than its predictor. Every design is one focal regressor after a nuisance design, so every arm runs through the same engine and the same permutation stream. Selected units that overlap are merged into regions with genomic coordinates, the probes and their target genes are listed within, and a summary declares the run: calibration first, then the regions, then the standing of each arm as measured by the registered null report NR2.

This vignette runs everything on a small synthetic array with three planted loci, so it completes in seconds and the truth is known. A synthetic array carries no anchor-panel probe, so the knowledge layer is off throughout; on a real array it is on by default.

# A synthetic array with planted loci

```{r simulate}
library(tewas)
sim <- tewas_simulate(m = 900, n = 120, seed = 7)
str(sim, max.level = 1)
lengths(sim$planted)
```

`beta` is probes × samples with CpG rownames, `y` a continuous phenotype, `X` two covariates, `mod` a continuous moderator, and `planted` the rows of a linear locus, a pure `y × mod` interaction and a pure curvature.

# The linear scan: `tewas_v16()`

`tewas_v16()` is the registered pipeline of the methods paper. It applies the QC of the registered kits, transforms a continuous outcome by rank inverse-normal, computes the per-probe statistics, builds the units, and runs the engine.

```{r v16}
res <- tewas_v16(sim$beta, sim$y, sim$X, sim$chr, sim$pos, B = 250L, iqr_min = 0, seed = 1L)
res$counts
res$probes[rej_v16 == TRUE, .(cpg, chr, pos, z, p)]
sim$planted$linear
```

The rejected probes are the planted linear locus. On a real array `res$ruler` carries the alive ruler's raw rejection rate, empirical-null standard deviation and $\hat\pi_0$; `res$lens` the prior's gain and whether the guard tripped; `res$diagnostic$branch` whether the scan's inflation is low-rank structure a ruler could remove (report-only).

# The pieces, one at a time

The pipeline is three calls the user can also make directly.

```{r pieces}
M <- log2(sim$beta / (1 - sim$beta))
S <- tewas_scan_stats(M, tewas_inverse_normal(sim$y), sim$X)   # per-probe OLS statistics
U <- tewas_run_units(sim$chr, sim$pos)                            # runs at 500 and 2,000 bp
table(U$unit_type)
fit <- tewas_engine2(S$z, S$p, S$yr, S$Xs, S$dfr, U$units, U$weights, sim$chr, sim$pos,
                     unit_type = U$unit_type, B = 250L, lens = FALSE, seed = 1L)
head(fit$unit[fit$unit$selected, c("k", "stratum", "p_unit", "ihw_w", "credit")])
identical(which(fit$rej), which(res$probes$rej_v16))
```

# A declared design: `tewas_frame()`

The frame layer declares the design. Here the phenotype is the outcome of methylation, a moderator is declared, and the departure arms are on.

```{r frame}
fr <- tewas_frame(sim$beta, sim$y, sim$X, sim$chr, sim$pos,
                  moderation = TRUE, mod = sim$mod, type = "non-linear",
                  knowledge = FALSE, ruler = NULL, B = 250L, iqr_min = 0, seed = 1L,
                  diagnostic = FALSE, labels = list(y = "phenotype", mod = "moderator"))
fr
```

Each arm has its own unit table, its DMRs (selected units merged when they overlap) and its probe table.

```{r arms}
names(fr$arms)
fr$arms$linear$dmrs[, .(rank, coord, n_probes, p_unit, n_rejected)]
fr$arms$moderation$dmrs[, .(rank, coord, n_probes, p_unit, n_rejected)]
fr$arms$quadratic$dmrs[, .(rank, coord, n_probes, p_unit, n_rejected)]
```

The linear arm finds the planted linear locus, the moderation arm the planted interaction and the quadratic arm the planted curvature; each arm is silent at the loci planted for the others, which is the specificity of the decomposition.

```{r planted}
list(linear = sim$planted$linear, interaction = sim$planted$interaction, curvature = sim$planted$curvature)
```

The run declares what it did, including the standing of each arm as measured by the registered null report NR2:

```{r notes}
cat(fr$notes, sep = "\n\n")
```

# Bioconductor classes: `SummarizedExperiment` in, `GRanges` out

Both entry points accept a `SummarizedExperiment` in place of the beta matrix, a minfi `GenomicRatioSet` included: the beta values come from the `"Beta"` assay (or `"M"`, converted; or the assay named by `assay =`), the coordinates from `rowRanges()`, and `y`, `X`, `mod` may name columns of `colData()`. The object is unpacked into the inputs of the matrix path, so the run is the same run.

```{r se}
se <- SummarizedExperiment::SummarizedExperiment(
  assays = list(Beta = sim$beta),
  rowRanges = GenomicRanges::GRanges(sim$chr, IRanges::IRanges(sim$pos, sim$pos)),
  colData = data.frame(y = sim$y, sim$X, mod = sim$mod))
fr_se <- tewas_frame(se, "y", c("age", "sex"), moderation = TRUE, mod = "mod", type = "non-linear",
                     knowledge = FALSE, ruler = NULL, B = 250L, iqr_min = 0, seed = 1L, diagnostic = FALSE)
identical(fr_se$arms$linear$units$selected, fr$arms$linear$units$selected)
```

`tewas_granges()` returns the regions, the units or the probes of any arm as a `GRanges` with the table's columns as metadata, for overlap with annotation or export.

```{r granges}
tewas_granges(fr, "dmrs", "moderation")
gp <- tewas_granges(fr, "probes", "linear"); gp[gp$rejected]
```

# The report: `tewas_report()`

`tewas_report()` writes, for every arm, the unit table with coordinates, the DMRs and the probe tables (csv, html, md), the genome-wide unit scan, locus panels, the calibration figure, the link-trajectory figures (the adjusted unit score against the phenotype with the fitted trajectory of every arm that ran), simple slopes and the Johnson–Neyman band for moderation, and `summary.md`.

```{r report}
out <- file.path(tempdir(), "tewas_vignette")
files <- tewas_report(fr, out, top = 2, genes = "none")
basename(files)
```

`summary.md` reads from the top: design, calibration, the DMR tables, how to read them, and a Methods and a Results paragraph for the actual run.

```{r summary}
sm <- readLines(file.path(out, "summary.md"))
cat(head(sm, 14), sep = "\n")
```

```{r fig-link, fig.cap = "The link-trajectory panel of the top linear DMR: the covariate-adjusted unit score against the phenotype, with the fitted trajectory."}
knitr::include_graphics(file.path(out, "fig_link_linear.png"))
```

```{r fig-mod, fig.cap = "The moderation panel of the top moderated DMR: simple slopes at the observed 10th, 50th and 90th percentiles of the moderator, and the Johnson-Neyman band."}
knitr::include_graphics(file.path(out, "fig_moderation.png"))
```

# Target genes: `tewas_genes()`

With `cpgdirection` installed, `tewas_report(genes = "auto")` lists, for every probe, unit and region, the target genes through `cpgdirection::cpg_gene_pairs()`, with the mapping provenance separate from the pair's methylation–expression direction. Any function with the same contract can stand in; the synthetic probe identifiers map to no gene, so a stand-in is used here.

```{r genes}
stub <- function(cpg_ids, tissue = "blood", ...) {
  P <- data.frame(cpg_id = cpg_ids, target_gene = paste0("GENE", seq_along(cpg_ids) %% 3),
                  mapping_primary = "EPICv2_manifest", mapping_sources = "EPICv2_manifest",
                  best_direction = 1, best_evidence = "catalogue", direction_tier = "B",
                  usable = TRUE, probe_masked = FALSE, stringsAsFactors = FALSE)
  attr(P, "qc_excluded_cpgs") <- character(0); P }
fg <- tewas_genes(fr, annotator = stub)
fg$arms$linear$dmrs[, .(coord, genes)]
```

# The standing of each arm

The linear arm inherits the calibration record of the registered contests. The other arms were calibrated by the registered null report NR2 on a real substrate: the moderation (continuous or binary moderator), quadratic, threshold and moderated-curvature arms held their null bands; the exponential arm (`type = "exponential"`) did not and remains opt-in, its reading printed beside any result; a moderator with excess kurtosis above 1 is reported with a caution. Every frame run prints this standing in `fr$notes` and in the Methods paragraph of `summary.md`, so a reader of the output knows what was measured and what was not.

# Session information

```{r session}
sessionInfo()
```