tewas: unit-level multiple testing for EWAS, and a frame layer for behavioral designs
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
library(tewas)
sim <- tewas_simulate(m = 900, n = 120, seed = 7)
str(sim, max.level = 1)
#> List of 7
#> $ beta : num [1:900, 1:120] 0.45 0.452 0.355 0.486 0.507 ...
#> ..- attr(*, "dimnames")=List of 2
#> $ y : num [1:120] 1.302 0.979 1.786 0.375 0.914 ...
#> $ X : num [1:120, 1:2] 0.0903 -0.4065 -0.444 0.7752 0.4707 ...
#> ..- attr(*, "dimnames")=List of 2
#> $ mod : num [1:120] 2.328 -0.853 0.401 -0.353 -1.946 ...
#> $ chr : chr [1:900] "chr1" "chr1" "chr1" "chr1" ...
#> $ pos : int [1:900] 1202 1495 15631 15752 15910 28443 28829 29135 29368 44545 ...
#> $ planted:List of 3
lengths(sim$planted)
#> linear interaction curvature
#> 4 4 4beta 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.
res <- tewas_v16(sim$beta, sim$y, sim$X, sim$chr, sim$pos, B = 250L, iqr_min = 0, seed = 1L)
#> tewas_v16: no anchor-panel probe is present on this array; running without the knowledge layer (no prior, no lens, no guard).
res$counts
#> samples probes_qc probes_in_scan tested readmitted
#> 120 900 900 900 0
#> ruler_alive units rej_v16 rej_noprior bh
#> 0 426 4 4 4
res$probes[rej_v16 == TRUE, .(cpg, chr, pos, z, p)]
#> cpg chr pos z p
#> <char> <char> <int> <num> <num>
#> 1: cg00000006 chr1 28443 8.224030 1.967765e-16
#> 2: cg00000007 chr1 28829 8.231762 1.844793e-16
#> 3: cg00000008 chr1 29135 9.741523 2.005269e-22
#> 4: cg00000009 chr1 29368 6.666153 2.625951e-11
sim$planted$linear
#> [1] 6 7 8 9The 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.
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)
#>
#> run2000 run500
#> 213 213
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")])
#> k stratum p_unit ihw_w credit
#> 3 4 run500:4-7 1e-300 0.4657659 simes
#> 216 4 run2000:4-7 1e-300 0.3920500 simes
identical(which(fit$rej), which(res$probes$rej_v16))
#> [1] TRUEA 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.
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
#> tewas_frame: predictor frame, 120 samples, 900 probes in scan, 426 units; arms: linear, moderation, quadratic, threshold, mod_quadratic
#> linear DMRs 1 (units 2) | probes rejected 4 | ruler rate NA sigma0 NA | guard - | gain applied NA
#> moderation DMRs 1 (units 2) | probes rejected 4 | ruler rate NA sigma0 NA | guard - | gain applied NA
#> quadratic DMRs 1 (units 2) | probes rejected 4 | ruler rate NA sigma0 NA | guard - | gain applied NA
#> threshold DMRs 0 (units 0) | probes rejected 0 | ruler rate NA sigma0 NA | guard - | gain applied NA
#> mod_quadratic DMRs 0 (units 0) | probes rejected 0 | ruler rate NA sigma0 NA | guard - | gain applied NA
#> shape omnibus: 2 units selected, 8 probesEach arm has its own unit table, its DMRs (selected units merged when they overlap) and its probe table.
names(fr$arms)
#> [1] "linear" "moderation" "quadratic" "threshold"
#> [5] "mod_quadratic"
fr$arms$linear$dmrs[, .(rank, coord, n_probes, p_unit, n_rejected)]
#> rank coord n_probes p_unit n_rejected
#> <int> <char> <int> <num> <int>
#> 1: 1 chr1:28,443-29,368 4 1e-300 4
fr$arms$moderation$dmrs[, .(rank, coord, n_probes, p_unit, n_rejected)]
#> rank coord n_probes p_unit n_rejected
#> <int> <char> <int> <num> <int>
#> 1: 1 chr2:514,223-514,853 4 1e-300 4
fr$arms$quadratic$dmrs[, .(rank, coord, n_probes, p_unit, n_rejected)]
#> rank coord n_probes p_unit n_rejected
#> <int> <char> <int> <num> <int>
#> 1: 1 chr3:983,528-984,588 4 1e-300 4The 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.
list(linear = sim$planted$linear, interaction = sim$planted$interaction, curvature = sim$planted$curvature)
#> $linear
#> [1] 6 7 8 9
#>
#> $interaction
#> [1] 458 459 460 461
#>
#> $curvature
#> [1] 896 897 898 899The run declares what it did, including the standing of each arm as measured by the registered null report NR2:
cat(fr$notes, sep = "\n\n")
#> The continuous phenotype was rank-inverse-normal transformed before the scan (registered rule; heavy tails break the residual permutation).
#>
#> The moderator enters mean-centered on its own scale (not transformed); simple slopes are read at its observed 10th, 50th and 90th percentiles.
#>
#> Run without the knowledge layer: no prior, no lens, no anchors; the selection is the v1.4 (no-prior) selection.
#>
#> Only 0 alive-ruler probes present: no calibration report (rate, sigma0, guard) for this scan.
#>
#> Shape omnibus: departure arms combined per unit by the Cauchy combination and selected by BH over units (unweighted); a curvilinear component is not evidence of a U shape - read the Sasabuchi test and the Fieller interval of the rejected probes.
#>
#> Standing (NR2, registration V18, run 2026-09-14): the moderation (continuous or binary moderator), quadratic, threshold, moderated curvature arms held their null bands on the contest-3b substrate and are calibrated for the regressors tested; the exponential arm did not (10 of 20 permuted-outcome scans with a selection (band <= 0.10; one scan 155 units); unit-p ratio at 1e-3 of 2.24 (band 0.5-2)) and is opt-in.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.
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)
#> [1] TRUEtewas_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.
tewas_granges(fr, "dmrs", "moderation")
#> GRanges object with 1 range and 16 metadata columns:
#> seqnames ranges strand | dmr n_probes
#> <Rle> <IRanges> <Rle> | <integer> <integer>
#> chr2:514,223-514,853 chr2 514223-514853 * | 1 4
#> n_units units unit_types p_unit
#> <integer> <character> <character> <numeric>
#> chr2:514,223-514,853 2 109;322 run2000;run500 1e-300
#> q_unit credit direction mean_z
#> <numeric> <character> <character> <numeric>
#> chr2:514,223-514,853 6.66944e-298 simes + 8.80211
#> family family_frac n_rejected
#> <character> <numeric> <integer>
#> chr2:514,223-514,853 none 0 4
#> rejected_cpgs coord rank
#> <character> <character> <integer>
#> chr2:514,223-514,853 cg00000458;cg0000045.. chr2:514,223-514,853 1
#> -------
#> seqinfo: 1 sequence from an unspecified genome; no seqlengths
gp <- tewas_granges(fr, "probes", "linear"); gp[gp$rejected]
#> GRanges object with 4 ranges and 16 metadata columns:
#> seqnames ranges strand | cpg pos unit
#> <Rle> <IRanges> <Rle> | <character> <integer> <integer>
#> cg00000006 chr1 28443 * | cg00000006 28443 216
#> cg00000007 chr1 28829 * | cg00000007 28829 216
#> cg00000008 chr1 29135 * | cg00000008 29135 216
#> cg00000009 chr1 29368 * | cg00000009 29368 216
#> z p t beta_m se_m beta_e
#> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
#> cg00000006 8.22403 1.96776e-16 9.60798 0.807152 0.0840085 0.549024
#> cg00000007 8.23176 1.84479e-16 9.61992 0.906700 0.0942524 0.489422
#> cg00000008 9.74152 2.00527e-22 12.15481 0.842254 0.0692939 0.665087
#> cg00000009 6.66615 2.62595e-11 7.37709 0.608460 0.0824797 0.524824
#> bh rejected rejected_noprior family tested
#> <logical> <logical> <logical> <character> <logical>
#> cg00000006 TRUE TRUE TRUE none TRUE
#> cg00000007 TRUE TRUE TRUE none TRUE
#> cg00000008 TRUE TRUE TRUE none TRUE
#> cg00000009 TRUE TRUE TRUE none TRUE
#> readmitted ruler
#> <logical> <logical>
#> cg00000006 FALSE FALSE
#> cg00000007 FALSE FALSE
#> cg00000008 FALSE FALSE
#> cg00000009 FALSE FALSE
#> -------
#> seqinfo: 3 sequences from an unspecified genome; no seqlengthsThe 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.
out <- file.path(tempdir(), "tewas_vignette")
files <- tewas_report(fr, out, top = 2, genes = "none")
basename(files)
#> [1] "unit_table_linear.csv" "units_selected_linear.csv"
#> [3] "dmrs_linear.csv" "probe_table_linear.csv"
#> [5] "probes_rejected_linear.csv" "fig_genome_linear.png"
#> [7] "fig_loci_linear.png" "unit_table_moderation.csv"
#> [9] "units_selected_moderation.csv" "dmrs_moderation.csv"
#> [11] "probe_table_moderation.csv" "probes_rejected_moderation.csv"
#> [13] "fig_genome_moderation.png" "fig_loci_moderation.png"
#> [15] "unit_table_quadratic.csv" "units_selected_quadratic.csv"
#> [17] "dmrs_quadratic.csv" "probe_table_quadratic.csv"
#> [19] "probes_rejected_quadratic.csv" "fig_genome_quadratic.png"
#> [21] "fig_loci_quadratic.png" "unit_table_threshold.csv"
#> [23] "units_selected_threshold.csv" "probe_table_threshold.csv"
#> [25] "probes_rejected_threshold.csv" "fig_genome_threshold.png"
#> [27] "unit_table_mod_quadratic.csv" "units_selected_mod_quadratic.csv"
#> [29] "probe_table_mod_quadratic.csv" "probes_rejected_mod_quadratic.csv"
#> [31] "fig_genome_mod_quadratic.png" "calibration.csv"
#> [33] "fig_link_linear.png" "fig_link_quadratic.png"
#> [35] "fig_moderation.png" "shape_table.csv"
#> [37] "shape_probes.csv" "fig_shape.png"
#> [39] "design_notes.txt" "summary.md"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.
sm <- readLines(file.path(out, "summary.md"))
cat(head(sm, 14), sep = "\n")
#> # tewas run summary - 2026-10-15 11:04
#>
#> ## Design
#>
#> - Frame: **predictor** (methylation predicts the phenotype: a positive direction means higher methylation goes with a higher phenotype (beta_e is the phenotype change per M-value unit)).
#> - Phenotype: `phenotype` (continuous, rank-inverse-normal); n = 120 samples after listwise deletion.
#> - Covariates: `age`, `sex`.
#> - Arms: linear; moderation (y x mod); quadratic departure; threshold departure; moderated curvature ((y - c)^2 x mod).
#> - Moderator: `moderator`.
#> - Engine: tewas 0.99.2, v1.6 (frame layer v1.8; NR2 standing per arm in notes); B = 250, alpha = 0.05, alpha_local = 0.05, seed = 1; knowledge layer off.
#> - Probes: 900 after QC, 900 in the scan, 900 tested, 0 re-admitted anchors, 0 alive-ruler probes; 426 units (426 multi-probe).
#>
#> ## Calibration of this scan (read this first)
The link-trajectory panel of the top linear DMR: the covariate-adjusted unit score against the phenotype, with the fitted trajectory.
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.
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.
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)]
#> coord genes
#> <char> <char>
#> 1: chr1:28,443-29,368 GENE1;GENE2;GENE0The 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
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] tewas_0.99.2 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] Matrix_1.7-6 jsonlite_2.0.0
#> [3] compiler_4.6.1 BiocManager_1.30.27
#> [5] SummarizedExperiment_1.43.0 Biobase_2.73.2
#> [7] GenomicRanges_1.65.4 jquerylib_0.1.4
#> [9] png_0.1-9 IRanges_2.47.5
#> [11] Seqinfo_1.3.2 yaml_2.3.12
#> [13] fastmap_1.2.0 lattice_0.23-1
#> [15] XVector_0.53.0 R6_2.6.1
#> [17] S4Arrays_1.13.2 generics_0.1.4
#> [19] knitr_1.52 BiocGenerics_0.59.12
#> [21] DelayedArray_0.39.8 MatrixGenerics_1.25.0
#> [23] maketools_1.3.2 bslib_0.12.0
#> [25] rlang_1.3.0 cachem_1.1.0
#> [27] xfun_0.61 sass_0.4.10
#> [29] sys_3.4.3 viridisLite_0.4.3
#> [31] SparseArray_1.13.4 cli_3.6.6
#> [33] digest_0.6.39 grid_4.6.1
#> [35] lifecycle_1.0.5 S4Vectors_0.51.10
#> [37] evaluate_1.0.5 data.table_1.18.6.1
#> [39] buildtools_1.0.0 abind_1.4-8
#> [41] stats4_4.6.1 rmarkdown_2.32
#> [43] matrixStats_1.5.0 tools_4.6.1
#> [45] htmltools_0.5.9- What tewas is for
- A synthetic array with planted loci
- The linear scan:
tewas_v16() - The pieces, one at a time
- A declared design:
tewas_frame() - Bioconductor
classes:
SummarizedExperimentin,GRangesout - The report:
tewas_report() - Target genes:
tewas_genes() - The standing of each arm
- Session information