Package 'dmsa'

Description

The four 2026 Project Alpha builds, transcribed from the student pack's covariate_sets.csv including its never-list. Supply focal for the analysis at hand; everything else is fixed by the build.

Usage

alpha_design(build, focal, drop = character())

Arguments

Details

Build 2 and 3 are longitudinal: their contract asks for the time x predictor interaction, so pass e.g. focal = "time:BSI_Total_c".

Value

A dmsa_design.

Examples

d <- alpha_design(1, focal = "IC_T1_c")
d
# builds 2 and 3 are longitudinal: the contract asks for the time interaction
alpha_design(2, focal = "time:BSI_Total_c")
# a deviation from the contract is recorded on the object, never silent
alpha_design(1, focal = "IC_T1_c", drop = "Epi_T1")$dropped

Description

Read the Project Alpha gene codebook shipped with the package: 549 genes assigned to 30 biological systems.

Usage

alpha_gene_systems()

Value

data.frame with system_id, system, gene, n_kept, usability

Examples

gs <- alpha_gene_systems()
dim(gs)
length(unique(gs$system_id))
head(gs[, c("system_id", "system", "gene", "n_kept", "usability")], 4)

Description

Curated w_g per gene: +1 gene product increases the system's activation tone, -1 decreases it (e.g. NR3C1 in the HPA axis: more GR = stronger negative feedback = LOWER axis tone), 0 = not on the activation axis (readouts, specificity controls, ambiguous).

Usage

alpha_polarity()

Details

Returns the audited 2026c table: 1,234 genes across all 30 systems. A legacy 115-gene draft covering six systems also ships, and is used only as a fallback if the 2026c file is missing from the installation.

The 2026c table is why a system-level alignment works at all: an alpha_polarity() that covered only six systems made dmsa_align(level = "system") fail on most of the panel, since a gene with no w_g entry stops the alignment by design.

Value

data.frame with system_id, system, gene, w_g, role, confidence and, for the 2026c table, module, w_g_source, anchor, evidence and citation. The "status" attribute is "audited" or "draft".

Examples

pol <- alpha_polarity()
nrow(pol)
head(pol[, c("system", "gene", "w_g", "role")], 4)
# CRH drives the HPA axis (w_g = +1); NR3C1 is its brake, so more GR means
# stronger negative feedback and LOWER axis tone (w_g = -1)
pol[pol$gene %in% c("CRH", "NR3C1") & pol$system_id == 2,
    c("gene", "w_g", "role")]

Description

Reads the audited 2026c gene codebook that ships with the package and returns it as a dmsa_reference: one row per gene within its system and module, with the curated gene-to-system polarity attached.

Usage

alpha_reference(modules = NULL)

Arguments

Details

The reference is BIOLOGICAL - it defines system > module > gene and stops there. Which CpGs exist is a property of the methylation data a user supplies, not of the biology, so a gene with no probe in any particular cohort still belongs to its system here. Genes with no curated polarity are left UNRESOLVED; they are never silently treated as +1.

Value

A dmsa_reference.

Examples

ref <- alpha_reference()
ref
# the polarity dmsa_align() consumes for one system (2 = HPA axis)
head(dmsa_polarity_for(ref, "2"), 4)

Description

Derives the system, module and gene counts from the shipped file and checks the reference is well formed: every gene in exactly one system and module, no probe-derived filtering, and polarity either curated or explicitly unresolved. Counts are DERIVED, never asserted against literals.

Usage

alpha_reference_check(verbose = TRUE)

Arguments

Value

Invisibly, a list with the derived counts and the check results.

Examples

alpha_reference_check(verbose = FALSE)$counts

Description

Build the aligned sign s_j for a set of probes

Usage

dmsa_align(
  direction,
  genes = NULL,
  level = c("gene", "system"),
  polarity = "alpha",
  system_id = NULL,
  missing_polarity = c("error", "zero", "drop"),
  tissue = "blood"
)

Arguments

Value

data.frame: probe, gene, d, p_plus, w_g, s (= d * w_g at system level, d at gene level), p_s_plus = P(s = +1) (chained), usable (logical), and reason for every non-usable probe.

Examples

# The short way: probe IDs alone. Directions and gene mapping come from
# the bundled map, so nothing has to be obtained first.
al <- dmsa_align(c("cg00052046", "cg00176879", "cg00308631"))
al[, c("probe", "gene", "d", "s", "usable")]

# Which map produced this alignment travels with it.
attr(al, "direction_map")

# Pin each probe to one gene when the panel says which gene it belongs to.
dmsa_align(c("cg00052046", "cg00176879"), genes = c("AVP", "AVP"))$s

# Direction is tissue-specific, so which layer to read is part of the
# analysis, not a default to leave alone. An adult saliva sample is about
# 82% immune and takes blood calls; a neonatal buccal swab is near-pure
# epithelium and takes epithelial ones.
bu <- dmsa_align("cg00176879", tissue = "epithelium")
attr(bu, "direction_map")

# d is the CpG -> expression sign, w_g the gene -> system-tone sign.
# A promoter CpG silences NR3C1 (d = -1) and GR is the HPA brake (w_g = -1),
# so the probe votes for HIGHER axis tone: s = d * w_g = +1.
g <- c("CRH", "NR3C1")
dcall <- data.frame(cpg = c("cg_crh", "cg_nr3c1"), d = c(1, -1),
                    p_plus = c(0.9, 0.1))
pol <- data.frame(gene = g, w_g = c(1, -1))
dmsa_align(dcall, genes = g, level = "gene")
dmsa_align(dcall, genes = g, level = "system", polarity = pol)

Description

The honest-reporting companion: how many probes each system can actually align, on which side of the activation axis. A system whose brake genes have no usable probes silently degrades to an activation-gene test - this table is how that is caught and reported.

Usage

dmsa_balance(alignment)

Arguments

Value

one-row data.frame per call

Examples

dcall <- data.frame(cpg = paste0("p", 1:4), d = c(1, -1, 1, NA),
                    p_plus = c(0.9, 0.1, 0.8, NA))
genes <- c("CRH", "POMC", "AVPR2", "NR3C1")
pol <- data.frame(gene = genes, w_g = c(1, 1, 0, -1))
al <- dmsa_align(dcall, genes, level = "system", polarity = pol)
# the only brake gene abstained, so this "system" test is really an
# activation-gene test - which is what the balance table makes visible
dmsa_balance(al)

Description

Largest set whose MEAN local fdr does not exceed q.

Usage

dmsa_bfdr_select(lfdr, q = 0.05)

Arguments

Value

integer indices of selected units

Examples

lfdr <- c(rep(0.001, 10), rep(0.5, 10), rep(0.99, 10))
sel <- dmsa_bfdr_select(lfdr, q = 0.05)
sel
mean(lfdr[sel])
# nothing is selected when no unit is credible enough to fit the budget
dmsa_bfdr_select(rep(0.9, 20), q = 0.05)

Description

Gated hierarchical DMSA with level-specific inference

Usage

dmsa_cascade(
  b,
  se,
  alignment,
  tree,
  nulls,
  engine = NULL,
  gate = c("both", "dense", "sparse"),
  q = 0.05,
  method = c("expected", "fixed"),
  borrow = 4
)

Arguments

Value

An object of class dmsa_cascade.

Examples

set.seed(1)
tree <- expand.grid(probe = 1:3, gene = 1:2, system = 1:4)
tree <- tree[, c("system", "gene")]      # outermost level first
P <- nrow(tree); n <- 60; x <- rnorm(n)
M <- matrix(rnorm(n * P), n, P) + outer(x, 0.9 * (tree$system == 1))
al <- dmsa_align(data.frame(cpg = paste0("p", 1:P), d = 1, p_plus = 0.95),
                 genes = paste0("g", tree$gene), level = "gene")
fit <- lm(M ~ x); b <- coef(fit)["x", ]
se <- sqrt(colSums(residuals(fit)^2) / (n - 2) / sum((x - mean(x))^2))
nulls <- dmsa_cascade_null(b, se, al, tree, exposure = x, M = M, B = 49)
cs <- dmsa_cascade(b, se, al, tree, nulls)
cs$tables$system

Description

Simulates the pooled aligned statistic under the global null for each unit shape present in the tree, so that dmsa_cascade() can convert observed statistics into calibrated p-values. Uses the supplied residual covariance structure by resampling the exposure within blocks.

Usage

dmsa_cascade_null(
  b,
  se,
  alignment,
  tree,
  exposure,
  M,
  block = NULL,
  B = 999L,
  method = c("expected", "fixed")
)

Arguments

Value

named list suitable for the nulls argument.

Examples

set.seed(1)
tree <- expand.grid(probe = 1:3, gene = 1:2, system = 1:4)
tree <- tree[, c("system", "gene")]      # outermost level first
P <- nrow(tree); n <- 60; x <- rnorm(n)
M <- matrix(rnorm(n * P), n, P)
al <- dmsa_align(data.frame(cpg = paste0("p", 1:P), d = 1, p_plus = 0.95),
                 genes = paste0("g", tree$gene), level = "gene")
fit <- lm(M ~ x); b <- coef(fit)["x", ]
se <- sqrt(colSums(residuals(fit)^2) / (n - 2) / sum((x - mean(x))^2))
# B = 49 keeps the example quick; a real calibration wants B >= 999
nulls <- dmsa_cascade_null(b, se, al, tree, exposure = x, M = M, B = 49)
vapply(nulls, length, integer(1))

Description

Tests Delta-S x exposure (x mod2) per unit through the composite machinery (dmsa_model, Freedman-Lane, family blocks), with maxT across the units of the level-local family on one shared permutation stream.

Usage

dmsa_change(
  M_before,
  M_after,
  data,
  outcome = NULL,
  exposure,
  mod2 = NULL,
  units,
  alignment,
  covariates = NULL,
  baseline = TRUE,
  role = c("predictor", "outcome"),
  block = NULL,
  B = 1999,
  correction = c("maxT", "minP"),
  weighting = c("combined", "reliability", "flat"),
  w_floor = 1.5,
  seed = 1
)

Arguments

Value

data.frame of class "dmsa_change": unit, n, n_probes, b, t, p, p_adj; attributes term, role, formula.

Examples

set.seed(5)
n <- 80; g <- rep(c("g1", "g2"), each = 4); P <- length(g)
d <- rep(c(1, -1), length.out = P); E <- rnorm(n); y <- rnorm(n)
M0 <- matrix(rnorm(n * P), n, P)
M1 <- 0.6 * M0 + matrix(rnorm(n * P), n, P)
# plant a change in g1 that tracks y x E through the alignment
M1[, g == "g1"] <- M1[, g == "g1"] + outer(0.5 * y * E, d[g == "g1"])
al <- dmsa_align(data.frame(cpg = paste0("cg", 1:P), d = d,
                            p_plus = ifelse(d > 0, 0.9, 0.1)), genes = g)
dat <- data.frame(y = y, E = E, dEpi = rnorm(n))
dmsa_change(M0, M1, dat, outcome = "y", exposure = "E", units = g,
            alignment = al, covariates = "dEpi", B = 99, seed = 1)

Description

Checks the things that silently ruin a methylation set analysis: a covariate that is constant in this subsample, a covariate that proxies the focal term, a grouping factor with too few levels to estimate, a grouping factor aliased with the focal term, and declared independence in data that plainly repeats.

Usage

dmsa_check_design(design, data, strict = TRUE)

Arguments

Value

Invisibly, a list with n, problems, notes and the per-term diagnostics.

Examples

set.seed(1)
d <- data.frame(exposure = rnorm(80), sex_c = rep(c(-0.5, 0.5), 40),
                age = rnorm(80, 40, 5), Fib = 0,
                chip = rep(1:10, each = 8), cID = rep(1:40, each = 2))
des <- dmsa_design("exposure", c("sex_c", "age"),
                   random = c("chip", "cID"), exchangeable = "cID")
dmsa_check_design(des, d)$n
# Fib is constant in this subsample: absorbed silently by lm, caught here
bad <- dmsa_design("exposure", c("sex_c", "Fib"), random = "chip",
                   exchangeable = "cID")
dmsa_check_design(bad, d, strict = FALSE)$problems

Description

Correct p-values across the contrasts of one categorical exposure

Usage

dmsa_contrast_adjust(x, p, direction = NULL)

Arguments

Value

A data frame with raw and adjusted p-values.

Examples

set.seed(1)
d <- data.frame(style = factor(sample(c("secure", "anxious", "avoidant"), 90,
                                      TRUE)),
                age = rnorm(90), cID = rep(1:45, each = 2))
des <- dmsa_design(focal = "style", fixed = "age", exchangeable = "cID")
cc <- dmsa_contrasts(d, des, reference = "secure")
# one p-value per directional contrast, corrected within the family
dmsa_contrast_adjust(cc, c(anxious = 0.01, avoidant = 0.04),
                     direction = c("+", "-"))

Description

DMSA pools signed per-probe effects, so a focal factor with k > 2 levels cannot be tested as one directional term. This builds the k - 1 single-contrast designs, each of which is directional, and states the multiplicity correction that applies across them.

Usage

dmsa_contrasts(
  data,
  design,
  reference = NULL,
  correct = c("holm", "BH", "bonferroni", "none")
)

Arguments

Value

An object of class dmsa_contrasts: a list of dmsa_design objects in $designs, the data frame with the contrast columns added in $data, and the correction in $correct. Feed each design to dmsa_fit() and pass the resulting p-values to dmsa_contrast_adjust().

Examples

set.seed(1)
d <- data.frame(style = factor(sample(c("secure", "anxious", "avoidant"), 40,
                                      TRUE)),
                age = rnorm(40), cID = rep(1:20, each = 2))
des <- dmsa_design(focal = "style", fixed = "age", exchangeable = "cID")
dmsa_contrasts(d, des)

Description

What the biological reference asked for, what the data at hand can test, and why every discovered CpG x gene pair was or was not used. Only available for frames built on the pair path (direction_source = "cpgdirection"); a frame on the bundled snapshot has no pair ledger and gets an error saying so.

Usage

dmsa_coverage(frame, level = c("system", "gene", "pairs"))

Arguments

Value

A data.frame (invisible its print). For "system" the columns are system_id, system, n_genes_reference, n_genes_testable, n_genes_untestable, n_measurements, n_pairs_used, n_pairs_discovered, family_size_gene (= n_genes_testable, the size of the maxT/minP family at the gene level).

Examples

## needs cpgdirection; see the vignette for a full run

Description

Returns the CpG-site column names of x that survive an include/exclude filter on the NAMES. Patterns are fixed substrings (not regular expressions), matched case-sensitively; a column survives include if it contains ANY of the include patterns, and is then removed by exclude if it contains ANY of the exclude patterns. Only CpG-site columns (names starting cg...) are considered - outcomes and covariates are never touched. The same filter is available directly inside dmsa_frame() as cpgs_include / cpgs_exclude.

Usage

dmsa_cpg_columns(x, include = NULL, exclude = NULL)

Arguments

Value

Character vector of surviving CpG column names.

Examples

cols <- c("cg00247334_BC21_T4_OXTR",
          "cg00247334_BC21_maternal_T1_OXTR",
          "cg00247334_BC21_paternal_T1_OXTR",
          "EPDS_Total_T4")
## the child's probes: "not the parents'"
dmsa_cpg_columns(cols, exclude = c("_maternal", "_paternal"))
## the mother's probes
dmsa_cpg_columns(cols, include = "_maternal")

Description

Declare a DMSA design: covariates, dependence, and permutation blocks

Usage

dmsa_design(
  focal,
  fixed,
  random = NULL,
  exchangeable = NULL,
  forbid = character(),
  label = NULL,
  id = NULL,
  format = c(NA, "wide", "long")
)

Arguments

Value

An object of class dmsa_design.

Examples

dmsa_design(focal = "IC_c",
            fixed = c("sex_c", "age", "Fib", "ctrlSV3", "ctrlSV5"),
            random = c("chip", "cID"), exchangeable = "cID",
            forbid = c("ctrlSV4", "plate"))

Description

Returns per-probe methylation-to-expression direction calls from the maps bundled with dmsa, so an alignment can be built without obtaining a direction resource separately.

Usage

dmsa_directions(probes = NULL, genes = NULL, tissue = "blood")

Arguments

Details

Two layers of cpgdirection are bundled, each restricted to the probe-gene pairs that carry an actual call. About seven in ten rows of the full table abstain, and an abstention cannot be aligned to, so those rows are not shipped - which is why the reduction costs nothing DMSA could have used.

"blood"

Whole blood: 488,204 pairs over 237,761 probes and 5,376 genes. The externally validated layer, and the right calls for an immune-dominated sample - adult saliva by passive drool runs about 82 immune.

"epithelium"

Nasal epithelium: 688,843 pairs over 277,830 probes and 4,447 genes. The right calls for an epithelium-dominated sample - neonatal and infant buccal swabs run .89 to 1.00 epithelial, where blood calls would be the wrong tissue.

Choose the layer by what the sample is made of, not by what it is called. The failure mode this guards against is quiet: an epithelial sample scored against immune-derived calls reports full coverage while every direction is drawn from the wrong tissue.

The brain bridge, the SMR causal layer, the remaining tissues and the full 3.1M-pair table live in cpgdirection (doi:10.5281/zenodo.22024185), which remains the complete resource.

A probe may map to more than one gene - a median of two, up to seven - and every mapped pair is returned, because which gene a probe is read against is a modelling decision rather than a lookup.

Value

A data.frame with probe, gene, d (-1 or +1), p_plus (the probability the direction is +1) and tier (evidence grade). The map's provenance and version are attached as the "dmsa_direction_map" attribute.

Examples

# What the bundled blood map knows about one gene's probes.
avp <- dmsa_directions(genes = "AVP")
nrow(avp)
head(avp, 3)

# Coverage of a probe set you care about.
probes <- c("cg00000029", "cg00000165", "not_a_real_probe")
got <- dmsa_directions(probes)
got

# The layers are not nested: each calls pairs the other declines to call,
# so coverage of the same probe set differs by tissue.
dmsa_directions(probes, tissue = "epithelium")

# Provenance travels with the map.
attr(dmsa_directions(), "dmsa_direction_map")[c("tissue", "version")]

Description

Every module label in the bundled cascade was checked against the literature. This returns what that check concluded, so a module-level finding can be reported with the strength of its own definition attached: whether the label was retained or revised, whether the locked membership is homogeneous, the evidence tier, and the citation keys behind it.

Usage

dmsa_evidence(x, which = c("all", "moderate", "flagged"))

Arguments

Value

data.frame of class dmsa_evidence, one row per module.

Examples

dmsa_evidence(dmsa_select(systems = "immune"), which = "flagged")

Description

Fit, pool and permute a DMSA set test under a declared design

Usage

dmsa_fit(
  data,
  probes,
  alignment,
  design,
  method = c("expected", "fixed"),
  B = 999L,
  engine = c("auto", "lm", "lmer"),
  beta_input = TRUE,
  check = TRUE,
  seed = NULL
)

Arguments

Value

An object of class dmsa_fit.

Examples

set.seed(1)
n <- 80; K <- 6; x <- rnorm(n); d <- sample(c(-1, 1), K, TRUE)
M <- matrix(rnorm(n * K), n, K) + outer(x, 0.5 * d)
dat <- data.frame(cID = rep(1:40, each = 2), chip = rep(1:10, each = 8),
                  x = x, sex_c = rep(c(-0.5, 0.5), 40))
probes <- paste0("cg", 1:K); dat[probes] <- 1 / (1 + 2^(-M))
al <- dmsa_align(data.frame(cpg = probes, d = d,
                            p_plus = ifelse(d > 0, 0.95, 0.05)),
                 genes = rep(c("g1", "g2"), each = 3))
des <- dmsa_design("x", "sex_c", random = c("chip", "cID"),
                   exchangeable = "cID")
dmsa_fit(dat, probes, al, des, B = 99, engine = "lm", seed = 1)

Description

Call 1 of the two-call interface. Declares the data, the levels (system / module / gene / probe), the outcomes, the covariates and blocks, optional moderation, and the direction map - validates all of it, runs the test drive (with autofix = TRUE correcting what is safely correctable, LOUDLY), and returns a dmsa_frame for dmsa_report.

Usage

dmsa_frame(
  data,
  methylation = NULL,
  map = "alpha",
  reference = "alpha",
  direction_source = c("auto", "cpgdirection", "bundled"),
  tissue = "blood",
  replicate_policy = c("collapse_site", "keep_correlated"),
  cpgs_include = NULL,
  cpgs_exclude = NULL,
  sample_id = NULL,
  sample_align = c("id", "positional"),
  missing_methylation = c("error", "drop_probes", "common_complete_rows"),
  methylation_scale = c("auto", "beta", "M"),
  system = TRUE,
  module = FALSE,
  gene = TRUE,
  probe = TRUE,
  systems = NULL,
  sets = "alpha",
  ...,
  outcomes = NULL,
  predictors = NULL,
  covariates = "contract",
  blocks = NULL,
  random_effects = NULL,
  chip = TRUE,
  chip_effect = c("random", "fixed", "none"),
  outcome_levels = NULL,
  predictor_levels = NULL,
  moderation = FALSE,
  mod = "",
  mod2 = "",
  outcome_label = NULL,
  predictor_labels = NULL,
  mod_label = NULL,
  mod2_label = NULL,
  covariate_labels = NULL,
  mod_levels = NULL,
  type = c("linear", "non-linear", "exponential"),
  outcome_type = c("gaussian", "logistic", "multinomial"),
  frame_role = c("predictor", "outcome"),
  cpg_map = c("full", "confidence"),
  B = 1999,
  alpha = 0.05,
  seed = 1,
  correction = c("maxT", "minP"),
  weighting = c("combined", "reliability", "flat"),
  w_floor = 1.5,
  palette = "viridis",
  plots = TRUE,
  tables = TRUE,
  summary = TRUE,
  gene_models = "auto",
  progress = interactive(),
  beep = interactive(),
  plot_type = c("png", "pdf"),
  table_type = c("html", "docx", "rtf"),
  outdir = "dmsa_output",
  autofix = TRUE
)

Arguments

Value

object of class dmsa_frame.

Examples

set.seed(42)
map <- data.frame(gene = rep(paste0("G", 1:3), each = 3), system_id = 1L,
                  system = "Sim system")
map$probe <- sprintf("cg%07d", seq_len(nrow(map)))
map$column <- paste0(map$probe, "_", map$gene)
map$best_direction <- rep(c(-1, 1), length.out = nrow(map))
map$p_plus <- ifelse(map$best_direction > 0, 0.9, 0.1)
d <- data.frame(out1 = rnorm(60), cov1 = rnorm(60),
                cID = rep(1:30, each = 2))
sig <- 0.5 * outer(d$out1, map$best_direction * (map$gene == "G1"))
d[map$column] <- plogis(matrix(rnorm(60 * nrow(map)), 60) + sig)
dmsa_frame(d, map = map, outcomes = "out1", covariates = "cov1",
           blocks = "cID", B = 99, outdir = tempfile("dmsa_ex"))

Description

Gated hierarchical DMSA: one score, one model, level-local correction

Usage

dmsa_gate(
  data,
  M,
  map,
  alignment,
  formula,
  term,
  roots,
  block = NULL,
  alpha = 0.05,
  B = 1999,
  winsor = 3,
  gate = c("both", "sparse", "dense"),
  sparse_reach = "children",
  seed = NULL
)

Arguments

Value

Object of class dmsa_gate.

Examples

set.seed(1)
n <- 60; probes <- paste0("cg", 1:6)
map <- data.frame(gene = rep(c("OXTR", "AVP"), each = 3), probe = probes)
al <- dmsa_align(data.frame(cpg = probes, d = rep(c(1, -1), 3),
                            p_plus = rep(c(.9, .1), 3)),
                 genes = map$gene, level = "gene")
M <- matrix(rnorm(n * 6), n, 6, dimnames = list(NULL, probes))
d <- data.frame(age = rnorm(n), cID = rep(seq_len(n / 2), each = 2))
d$y <- as.numeric(scale(M %*% (2 * al$p_s_plus - 1))) + rnorm(n)
## the two genes are the pre-registered family; a probe is corrected
## inside its own gene only, never against all six
dmsa_gate(d, M, map, al, y ~ S + age, term = "S", roots = c("OXTR", "AVP"),
          block = d$cID, B = 99, seed = 1)

Description

Returns one row per drawable feature - coding exon slice, 5' UTR, 3' UTR, or an exon whose coding status the annotation did not resolve - for every transcript of a gene, with strand and genomic start/end. This is what dmsa_plot_locus draws the gene from, and it is deliberately a plain data frame so it can be built, inspected, corrected and cached like any other reference layer in this package.

Usage

dmsa_gene_model(
  gene,
  source = c("auto", "ensembl", "gff", "txdb", "ensdb", "table"),
  file = NULL,
  db = NULL,
  orgdb = NULL,
  table = NULL,
  species = "homo_sapiens",
  genome = "hg38",
  cache = NULL,
  host = "https://rest.ensembl.org",
  quiet = FALSE
)

Arguments

Value

A data frame with the columns in .GM_COLS, class dmsa_gene_model. Zero rows if the model could not be obtained - never a fabricated one.

See Also

dmsa_gene_model_check, dmsa_plot_locus

Examples

## Not run: 
## Needs network: one lookup plus one region query per gene.
gm <- dmsa_gene_model("NR3C1")
subset(gm, canonical)[, c("feature", "start", "end", "exon_rank")]

## Or from an annotation you already have, with no network at all:
gm <- dmsa_gene_model("NR3C1", source = "gff", file = "Homo_sapiens.gtf.gz")

## End(Not run)

Description

Check a gene model before drawing from it

Usage

dmsa_gene_model_check(gm, verbose = TRUE)

Arguments

Value

Invisibly, a list with ok and one entry per check.

Examples

## AVP sits on the minus strand, so its 5' UTR is the HIGH coordinate
gm <- data.frame(
  gene = "AVP", gene_id = "ENSG00000101200", chr = "chr20",
  gene_start = 3082556, gene_end = 3084724, strand = "-",
  transcript = "ENST00000380293", transcript_biotype = "protein_coding",
  canonical = TRUE, feature = c("utr5", "cds", "cds", "utr3"),
  start = c(3084665, 3084555, 3082700, 3082556),
  end   = c(3084724, 3084664, 3082802, 3082699),
  exon_rank = c(1, 1, 3, 3), source = "gff", genome = "hg38")
dmsa_gene_model_check(gm)$ok
## a feature outside the gene span means a neighbour's exons survived
bad <- gm; bad$start[1] <- 3000000; bad$end[1] <- 3000100
dmsa_gene_model_check(bad, verbose = FALSE)$ok

Description

Takes the end object of any mainstream methylation preprocessing pipeline (minfi, sesame, meffil, ewastools, ChAMP, RnBeads, methylprep - see Details) and returns a dmsa_import object holding a samples-by-probes beta matrix and an aligned phenotype data.frame. Pass the result straight to dmsa_frame as its first argument.

Usage

dmsa_import(x, ...)

## S3 method for class 'matrix'
dmsa_import(
  x,
  pheno = NULL,
  align_by = NULL,
  values = c("auto", "beta", "m"),
  orientation = c("auto", "probes_rows", "probes_cols"),
  id = NULL,
  format = NULL,
  source = "matrix",
  ...
)

## S3 method for class 'data.frame'
dmsa_import(
  x,
  pheno = NULL,
  align_by = NULL,
  values = c("auto", "beta", "m"),
  orientation = c("auto", "probes_rows", "probes_cols"),
  id = NULL,
  format = NULL,
  source = "data.frame",
  ...
)

## S3 method for class 'character'
dmsa_import(x, ...)

## S3 method for class 'list'
dmsa_import(x, pheno = NULL, align_by = NULL, id = NULL, format = NULL, ...)

## Default S3 method:
dmsa_import(
  x,
  pheno = NULL,
  align_by = NULL,
  values = c("auto", "beta", "m"),
  id = NULL,
  format = NULL,
  source = NULL,
  ...
)

Arguments

Value

An object of class dmsa_import: a list with data (the aligned phenotype data.frame), methylation (samples x probes beta matrix), source, values (scale detected), align (how samples were matched), id, format and notes. dmsa_frame() accepts it directly as its first argument.

Examples

## a probes-x-samples beta matrix, the shape most pipelines produce
set.seed(1)
B <- matrix(runif(60), nrow = 6,
            dimnames = list(sprintf("cg%08d", 1:6), paste0("S", 1:10)))
ph <- data.frame(Sample_Name = paste0("S", 1:10), age = rnorm(10))
imp <- dmsa_import(B, pheno = ph)
imp

Description

Gated hierarchical DMSA on the pooled aligned statistic

Usage

dmsa_levels(
  M,
  data,
  rhs,
  term,
  map,
  alignment,
  roots,
  block = NULL,
  alpha = 0.05,
  B = 1999,
  winsor = 3,
  method = c("expected", "fixed"),
  gate = c("both", "sparse", "dense"),
  family_correction = c("maxT", "minP"),
  sparse_reach = "children",
  seed = NULL
)

Arguments

Value

Object of class dmsa_gate.

Examples

set.seed(1)
n <- 60; probes <- paste0("cg", 1:6)
map <- data.frame(gene = rep(c("OXTR", "AVP"), each = 3), probe = probes)
al <- dmsa_align(data.frame(cpg = probes, d = rep(c(1, -1), 3),
                            p_plus = rep(c(.9, .1), 3)),
                 genes = map$gene, level = "gene")
d <- data.frame(anx = rnorm(n), age = rnorm(n),
                cID = rep(seq_len(n / 2), each = 2))
## the probes are the RESPONSES: their aligned coefficients are pooled,
## so cross-probe concordance is what carries the evidence
M <- matrix(rnorm(n * 6), n, 6) + outer(0.3 * d$anx, 2 * al$p_s_plus - 1)
dmsa_levels(M, d, rhs = c("anx", "age"), term = "anx", map = map,
            alignment = al, roots = c("OXTR", "AVP"), block = d$cID,
            B = 99, seed = 1)

Description

Efron-style two-group model on z-values: the null centre and scale are fitted on an iteratively trimmed core (robust to a minority signal slab), pi0 is the observed mass inside the window that holds 50 the fitted null, and the alternative is a wider normal. Returns one local fdr per unit.

Usage

dmsa_lfdr(z, prior_odds_multiplier = 1, min_n = 50L)

Arguments

Details

Assumptions, stated plainly (calibrated by simulation, 2026-08-29): signal must be a MINORITY of units - power is good to roughly a third signal and collapses conservatively (toward selecting nothing) beyond that; weak signal (|z| ~ 2) is largely invisible to the EB arm (also conservative); and the null must be approximately GAUSSIAN. A null with both tails significantly heavier than the fitted Gaussian - kurtosis, which the two-group model would misread as signal - is detected by a two-radius symmetric-tail-excess guard and the function returns NULL with a warning, the callers' documented fall-back to the frequentist rule. The guard keys on BOTH tails being in excess, so one-sided signal does not trip it (~3 symmetric two-sided signal DOES trip it and is handled by the frequentist arm instead - a conservative, never a silent, outcome. Residual risk: with only a few hundred units a heavy-tailed null is detected in about two runs of three; below that scale prefer the frequentist engine.

Value

numeric vector of local fdr values, or NULL if not estimable.

Examples

set.seed(1)
## 300 null genes and 30 carrying a real aligned effect
z <- c(rnorm(300), rnorm(30, mean = 4))
lfdr <- dmsa_lfdr(z)
round(quantile(lfdr, c(0, .5, 1)), 3)
mean(lfdr[301:330]) < mean(lfdr[1:300])
## too few units to fit the two-group model: NULL asks the caller to
## fall back to a frequentist rule rather than guess
is.null(dmsa_lfdr(rnorm(10)))

Description

Tests a moderated-mediation claim in which the predictors, moderator and mediator are latent SYSTEM TONES, each measured by its module tones as reflective indicators, and the outcome is an observed variable. The structural model runs on measurement-error-corrected moments (structural-after-measurement / Croon family), so the interaction paths are not attenuated by tone unreliability; inference is permutation-native (Freedman-Lane per path, joint-significance for products, Holm over the declared index family).

Usage

dmsa_mediate(
  data,
  methylation = NULL,
  map = "alpha",
  x,
  m,
  w = NULL,
  outcome,
  covariates = character(),
  measurement = c("tau", "congeneric", "single"),
  min_module_probes = 3L,
  single_rel = NULL,
  module_weights = c("equal", "sqrt_probes", "probes"),
  measure_covariates = NULL,
  chip = NULL,
  block = NULL,
  perm_blocks = c("none", "chip"),
  model = c("7", "8"),
  map_level = c("full", "confidence"),
  B = 1999L,
  boot = 0L,
  alpha_ssc = 0,
  conditional_at = c(-1, 0, 1),
  beta_input = TRUE,
  seed = NULL
)

Arguments

Value

An object of class dmsa_mediate: per-system measurement tables, path coefficient tables with permutation p per path, indices of moderated mediation with joint-significance + Holm p, conditional indirect effects, bootstrap CIs when requested, and the corrected moment machinery for inspection.

Examples

## synthetic end-to-end run at analysis scale (see tests for the
## simulation-gate harness)

set.seed(1)
n <- 120
eta <- matrix(rnorm(n * 3), n, 3)                 # hpa, oxytocin, kyn
mk <- function(f, p) sapply(seq_len(p), function(j)
  0.7 * f + rnorm(n, sd = sqrt(1 - 0.49)))
tones <- list(hpa = mk(eta[, 1], 6), oxytocin = mk(eta[, 2], 8),
              kyn = mk(0.4 * eta[, 1] * eta[, 2] + 0.3 * eta[, 1], 5))
dat <- data.frame(y = 0.5 * (0.4 * eta[, 1] * eta[, 2] + 0.3 * eta[, 1]) +
                    rnorm(n), age = rnorm(n))
fit <- dmsa:::.med_core(tones, dat, x = "hpa", m = "kyn", w = "oxytocin",
                        outcome = "y", covariates = "age",
                        measurement = "tau", B = 199, blocks = NULL,
                        model8 = FALSE, conditional_at = c(-1, 0, 1),
                        alpha_ssc = 0)
fit$indices

Description

The upper level of DMSA 2.0. Fits formula once, then permutes the reduced-model residuals within permutation blocks of equal size (Freedman-Lane) to obtain an exact p-value for one named coefficient - including an interaction. dmsa_moderate is the special case of this for score * moderator with a declared dmsa_design(); use that when you want the design contract enforced and simple slopes reported, and this when the question is an arbitrary term.

Usage

dmsa_model(
  formula,
  data,
  term,
  block = NULL,
  B = 1999,
  seed = NULL,
  nulls = FALSE
)

Arguments

Value

An object of class dmsa_model.

Examples

set.seed(1)
n <- 80
d <- data.frame(S = rnorm(n), ACE = rnorm(n),
                cID = rep(seq_len(n / 2), each = 2))   # couples
d$y <- 0.5 * d$S * d$ACE + rnorm(n)
## rows are relabelled only within equal-size families, so the couple
## dependence is preserved under the null
dmsa_model(y ~ S * ACE, d, term = "S:ACE", block = d$cID, B = 99, seed = 1)

Description

Refits the same model once per projection, substituting each into the column named score_col, and combines the p-values with the aggregated Cauchy test. ACAT is valid under arbitrary dependence, so hedging across projections costs no multiplicity - unlike a Bonferroni over arms, which is what the old min(1, 2*min(p)) gate amounted to.

Usage

dmsa_model_omnibus(
  scores,
  formula,
  data,
  term,
  block = NULL,
  score_col = "S",
  B = 1999,
  seed = NULL
)

Arguments

Value

list: p_by_flavour, the omnibus p, and the fits.

Examples

set.seed(1)
n <- 80
M <- matrix(rnorm(n * 6), n, 6)
al <- dmsa_align(data.frame(cpg = paste0("cg", 1:6), d = c(1, 1, -1, 1, -1, 1),
                            p_plus = rep(.9, 6)),
                 genes = rep("OXTR", 6), level = "gene")
sc <- dmsa_scores(M, al)
d <- data.frame(age = rnorm(n), cID = rep(seq_len(n / 2), each = 2))
d$y <- 0.6 * sc$aligned + rnorm(n)
## one p per projection, combined by ACAT - hedging costs no multiplicity
o <- dmsa_model_omnibus(sc, y ~ S + age, d, term = "S", block = d$cID,
                        B = 99, seed = 1)
round(c(o$p_by_flavour, omnibus = o$p), 3)

Description

Moderated DMSA (mDMSA)

Usage

dmsa_moderate(
  data,
  probes,
  alignment,
  outcome,
  moderator,
  design,
  method = c("expected", "fixed"),
  center = TRUE,
  scale_outcome = TRUE,
  B = 999L,
  engine = c("auto", "lm", "lmer"),
  beta_input = TRUE,
  check = TRUE,
  seed = NULL
)

Arguments

Value

An object of class dmsa_moderate.

Examples

set.seed(1)
n <- 100; probes <- paste0("cg", 1:6)
dat <- data.frame(sex_c = rep(c(-.5, .5), length.out = n),
                  chip = rep(1:10, each = 10),
                  cID = rep(seq_len(n / 2), each = 2), mod = rnorm(n))
dat[probes] <- plogis(matrix(rnorm(n * 6), n, 6))
al <- dmsa_align(data.frame(cpg = probes, d = rep(c(1, -1), 3),
                            p_plus = rep(c(.9, .1), 3)),
                 genes = rep("OXTR", 6), level = "gene")
dat$out <- 0.6 * dmsa_score(dat, probes, al) * dat$mod + rnorm(n)
des <- dmsa_design("placeholder", "sex_c", random = "chip", exchangeable = "cID")
dmsa_moderate(dat, probes, al, outcome = "out", moderator = "mod",
              design = des, B = 99, engine = "lm", seed = 3)

Description

M-value transform with clamping

Usage

dmsa_mvalues(beta, eps = 1e-04)

Arguments

Value

A numeric matrix or vector of M-values, the same shape as beta, on the logit scale.

Examples

beta <- matrix(runif(24, 0.05, 0.95), nrow = 6,
               dimnames = list(NULL, paste0("cg", 1:4)))
round(head(dmsa_mvalues(beta), 3), 2)
## clamping is what keeps a beta of exactly 0 or 1 from becoming +/-Inf
dmsa_mvalues(c(0, 0.5, 1))

Description

Pools the per-probe F statistic for the k - 1 contrast columns and calibrates it by the same block permutation the directional tests use. This answers "does this factor move the set at all" and gives up the direction, which is the property that separates DMSA from the methods it is compared against. Use it as a screen, or when the levels have no ordering and no reference level is defensible - and report dmsa_contrasts() alongside it.

Usage

dmsa_omnibus(data, Y, design, beta_input = TRUE, B = 999, seed = 1L)

Arguments

Value

A list with the pooled statistic, the permutation p-value, and an explicit directional = FALSE.

Examples

set.seed(1)
n <- 90
d <- data.frame(style = factor(sample(c("secure", "anxious", "avoidant"),
                                      n, TRUE)),
                age = rnorm(n), cID = rep(seq_len(n / 2), each = 2))
Y <- matrix(plogis(rnorm(n * 8)), n, 8)     # beta values for 8 probes
des <- dmsa_design(focal = "style", fixed = "age", exchangeable = "cID")
## three attachment styles with no defensible reference level: the pooled
## F asks only whether the factor moves the set, and gives up the direction
o <- dmsa_omnibus(d, Y, des, B = 99, seed = 5)
round(o$p_perm, 3)
o$directional    # always FALSE: a pooled F on 2 df has no sign

Description

The subject-level arm of DMSA for outcomes that are not continuous: binary (logistic), ordered-categorical (proportional odds), unordered-categorical (multinomial), or counts (Poisson / negative binomial). Optionally moderated.

Usage

dmsa_outcome(
  data,
  outcome,
  score,
  family = c("gaussian", "binomial", "ordinal", "multinomial", "poisson"),
  moderator = NULL,
  covariates = character(),
  block = NULL,
  center = TRUE,
  B = 999L,
  seed = NULL
)

Arguments

Details

Direction. A binary, count or ordinal outcome yields ONE coefficient for the tone score, so DMSA's directional claim survives intact - the sign of the log-odds (or log-rate, or latent-scale slope) is the expression-aligned direction. A multinomial outcome yields K-1 coefficients and there is no single sign: the test becomes an omnibus "tone predicts category" and the per-category log-odds are reported separately. If the categories are in any sense ordered, family = "ordinal" keeps the direction and is strictly preferable.

Value

An object of class dmsa_outcome.

Examples

set.seed(1)
n <- 120
d <- data.frame(S = rnorm(n), cv = rnorm(n),
                cid = rep(seq_len(n / 2), each = 2))
d$y <- rbinom(n, 1, plogis(0.9 * d$S))
## a binary outcome yields ONE coefficient, so the sign of the log-odds
## is still the expression-aligned direction
dmsa_outcome(d, outcome = "y", score = "S", family = "binomial",
             covariates = "cv", block = "cid", B = 99, seed = 2)

Description

Perturbs the analysis set around a fixed finding and records what happens to each engine's p-value. The focal unit's own probes are never touched, so any movement is attributable to the analysis set rather than to the evidence.

Usage

dmsa_pdance(
  engines,
  set,
  focal,
  pool = NULL,
  form = c("dropout", "addition"),
  grid = NULL,
  R = 100L,
  alpha = 0.05,
  seed = NULL
)

Arguments

Details

Two forms, both reported in the same way:

"dropout"

Removes a random fraction of the non-focal probes. The omitted probes carry no focal evidence, so a finding that depends on them was never about its own unit.

"addition"

Adds probes drawn from pool, which sits outside the analysis set. This is the direction a competitive engine is most sensitive to, because added probes enter its background.

The two forms move competitive and whole-set engines in opposite directions, which is why running only one of them understates the problem.

Value

A data.frame with one row per engine, grid point and replicate: engine, form, size (the fraction dropped or the count added), rep, p, and flip - whether that replicate's verdict at alpha differs from the unperturbed verdict. The unperturbed run is included as size = 0, rep = 0. The baseline p-value per engine is attached as the "baseline" attribute.

Examples

set.seed(1)
n <- 120L; K <- 40L
f <- rnorm(n)
M <- sapply(seq_len(K), function(j) if (j <= 6) 0.7 * f + rnorm(n) else rnorm(n))
colnames(M) <- sprintf("cg%04d", seq_len(K))
y <- 0.6 * f + rnorm(n)
focal <- colnames(M)[1:6]

# Two engines on identical evidence. The competitive one scores the focal
# unit against whatever else happens to be loaded; the unit-only one uses
# the focal probes and nothing else.
competitive <- function(u) {
  rest <- setdiff(u, focal)
  if (length(rest) < 2) return(NA_real_)
  t.test(abs(cor(M[, focal], y)), abs(cor(M[, rest], y)))$p.value
}
unit_only <- function(u) {
  summary(lm(rowMeans(M[, focal]) ~ y))$coefficients["y", 4]
}

pd <- dmsa_pdance(list(competitive = competitive, unit_only = unit_only),
                  set = colnames(M), focal = focal,
                  grid = c(0.1, 0.3, 0.5), R = 20, seed = 1)

# Both start from the same evidence, and it is never touched.
round(attr(pd, "baseline"), 4)

# How far each engine's p travels when only the OTHER probes change.
round(tapply(pd$p, pd$engine, function(x) diff(range(x, na.rm = TRUE))), 4)

# And how often the verdict at .05 reverses. The unit-only engine cannot
# move, because nothing it reads has changed.
tapply(pd$flip, pd$engine, sum)

Description

Permutation p-value for a DMSA statistic

Usage

dmsa_perm_pvalue(observed, null_stats)

Arguments

Value

two-sided permutation p with the +1 correction (resolution floor 1/(B+1)).

Examples

set.seed(1)
null_z <- rnorm(999)          # z from B family-wise permuted datasets
dmsa_perm_pvalue(3.1, null_z)
## the +1 correction puts a floor of 1/(B+1) on any permutation p
dmsa_perm_pvalue(50, null_z)

Description

Panel 2: gene level inside one system

Usage

dmsa_plot_genes(
  genes,
  null_max = NULL,
  title = "",
  alpha = 0.05,
  file = NULL,
  width = 1500,
  height = 1500,
  res = 190
)

Arguments

Value

Called for its side effect of drawing a plot on the active graphics device. Returns NULL invisibly.

Examples

set.seed(1)
genes <- data.frame(gene = paste0("G", 1:8), n_probes = sample(3:9, 8, TRUE),
                    z = c(rnorm(7), 3.6))
genes$p_adj <- 2 * pnorm(-abs(genes$z))
## the permutation null of max|z| over this same family of genes
null_max <- apply(matrix(rnorm(200 * 8), 200, 8), 1, function(r) max(abs(r)))
dmsa_plot_genes(genes, null_max = null_max, title = "HPA axis",
                file = tempfile(fileext = ".png"))

Description

Genomic coordinates are optional. Supply pos and the CpGs are placed to scale and the chromosome strip is drawn; supply gene_region instead and they are ordered along the gene (TSS1500 to 3'UTR) and spaced evenly; supply neither and they keep the order you gave. The axis states which of the three you got. Nothing else in the figure changes.

Usage

dmsa_plot_locus(
  probes,
  gene = "",
  chrom = NA,
  gene_start = NA,
  gene_end = NA,
  chrom_length = NA,
  signal_p = 0.05,
  invert = c("+1", "-1", "none"),
  order_by = c("auto", "pos", "region", "given"),
  gene_model = NULL,
  transcripts = c("canonical", "all"),
  context = "",
  file = NULL,
  width = 1800,
  height = NULL,
  res = 190
)

Arguments

Value

Invisibly, the plotted data frame, with the columns the figure actually used (x, sig, shown, flipped).

See Also

dmsa_gene_model to build the model, dmsa_plot_locus_gviz for the Gviz engine, dmsa_probe_coords for probe positions.

Examples

probes <- data.frame(probe = paste0("cg", 1:5), chr = "chr20",
                     pos = c(3082600, 3082750, 3083050, 3084600, 3084700),
                     b = c(-.04, .02, .03, -.05, -.06), se = rep(.015, 5),
                     d = c(-1, 1, 1, -1, -1), p = c(.01, .3, .04, .002, .001))
## invert = "+1" (the default) reflects the d = +1 probes, so a positive
## effect always means "moved methylation the way that LOWERS expression"
out <- dmsa_plot_locus(probes, gene = "AVP", context = "attachment anxiety",
                       file = tempfile(fileext = ".png"))
out[, c("probe", "d", "shown", "flipped")]

Description

The same probes and the same gene model as dmsa_plot_locus, handed to Gviz instead of drawn here. Use this when you want what Gviz gives and this package does not: a real cytoband ideogram, transcript collapsing, and the ability to stack any other GdObject you already build - CpG islands, ATAC peaks, ChIP tracks, conservation - in the same coordinate frame.

Usage

dmsa_plot_locus_gviz(
  probes,
  gene = "",
  gene_model = NULL,
  chrom = NA,
  genome = "hg38",
  extra = list(),
  ideogram = FALSE,
  file = NULL,
  width = 1800,
  height = 1100,
  res = 190
)

Arguments

Details

Gviz, GenomicRanges and IRanges are Suggests, not dependencies. This function errors with the install line if they are absent; everything else in the package keeps working without them.

Value

Invisibly, the list of tracks plotted.

See Also

dmsa_plot_locus for the dependency-free panel.

Examples

## Not run: 
gm <- dmsa_gene_model("NR3C1")
dmsa_plot_locus_gviz(pr, gene = "NR3C1", gene_model = gm)

## End(Not run)

Description

The panel that declares a direction. Each probe's slope is reflected according to its cpgdirection call so that every probe is read on one axis; reflected probes are drawn as open symbols and the convention is printed under the title.

Usage

dmsa_plot_probes(
  probes,
  invert = c("+1", "-1", "none"),
  gene_summary = NULL,
  title = "",
  axis_label = NULL,
  alpha = 0.05,
  file = NULL,
  width = 1700,
  height = 1400,
  res = 190
)

Arguments

Value

Called for its side effect of drawing a plot on the active graphics device. Returns NULL invisibly.

Examples

## Every probe is read on one axis. A CpG whose methylation predicts HIGHER
## expression (d = +1) is reflected, so positive always means the exposure
## moved methylation the way that LOWERS this gene's expression.
pr <- data.frame(
  gene  = "NR3C1",
  probe = paste0("cg", 10001:10006),
  b     = c(-0.021, -0.014, 0.018, -0.009, 0.026, -0.017),
  se    = c(0.008, 0.007, 0.009, 0.006, 0.010, 0.008),
  d     = c(-1, -1, 1, -1, 1, -1))
dmsa_plot_probes(pr, gene_summary = data.frame(gene = "NR3C1", z = -2.9),
                 title = "HPA axis", file = tempfile(fileext = ".png"))

Description

Distance from the diagonal is how concentrated a system's response is: a system null on the pooled arm and strong on the max-gene arm carries its signal in one gene.

Usage

dmsa_plot_systems(
  systems,
  alpha = 0.05,
  file = NULL,
  width = 1500,
  height = 1400,
  res = 190
)

Arguments

Value

Invisibly, the plotted data frame.

Examples

## Distance from the diagonal says how concentrated a system's response is:
## strong on the sparse arm only means the signal sits in a single gene.
sys <- data.frame(
  system   = c("HPA axis", "Oxytocin", "Serotonin", "Inflammation"),
  outcome  = "PTSD severity",
  p_dense  = c(0.0004, 0.08, 0.41, 0.62),
  p_sparse = c(0.0012, 0.0003, 0.33, 0.55),
  primary  = c(TRUE, TRUE, FALSE, FALSE),
  top_gene = c("NR3C1", "OXTR", "TPH2", "IL6"))
dmsa_plot_systems(sys, file = tempfile(fileext = ".png"))

Description

Load a gene-to-system polarity table

Usage

dmsa_polarity(x = "alpha", sets = NULL)

Arguments

Details

Signs are read against the system's declared activation tone, not against a gene's immediate partner. A gene that inhibits a brake therefore carries w_g = +1: FKBP5 restrains GR, GR mediates HPA negative feedback, so more FKBP5 means more axis drive. The role vocabulary names this explicitly (brake-of-brake, feedback-enabler) because it is the error that inverts a module.

Value

An object of class dmsa_polarity, or NULL.

See Also

dmsa_polarity_review() for the rows needing adjudication, dmsa_polarity_fetch() to draft signs for a panel of your own, dmsa_polarity_check() to validate a candidate table.

Examples

pol <- dmsa_polarity()
pol
dmsa_polarity_review(pol)

Description

Validate a candidate polarity table

Usage

dmsa_polarity_check(x = "alpha", sets = NULL, verbose = TRUE)

Arguments

Value

Invisibly, a list with ok and one entry per check.

Examples

## An anchor gene is what defines "more activation" for its system, so an
## anchor that brakes its own system is a contradiction, not a small error.
cand <- data.frame(system_id = "HPA", gene = c("CRH", "POMC", "NR3C1"),
                   w_g = c(-1, 1, -1), anchor = c(TRUE, FALSE, FALSE),
                   w_g_source = "curated")
chk <- dmsa_polarity_check(cand)
chk$ok
dmsa_polarity_check("alpha", verbose = FALSE)$ok

Description

Queries public signed-relationship resources for the genes of a cascade and returns a DRAFT polarity table for a curator to check. It does not decide anything on its own, and it is deliberately unwilling to guess.

Usage

dmsa_polarity_fetch(
  sets,
  anchors = NULL,
  sources = c("go", "omnipath"),
  max_genes = 200L,
  quiet = FALSE
)

Arguments

Details

What each resource can and cannot do is worth knowing before you rely on it. go is the only one that gives a gene-to-process sign directly, through positive/negative regulation of X terms, so it needs no path chaining - but its coverage is uneven gene to gene. signor carries explicit gene-to-phenotype rows, with a cell-biological vocabulary that fits metabolic and immune systems better than behavioural ones. omnipath has the broadest coverage and the weakest signs: CollecTRI and DoRothEA, two of its inputs, both document assigning an activating mode BY DEFAULT when they do not know, so a consensus_stimulation at curation_effort = 1 should be read as sign-unknown. trrust is the most honest transcription-factor resource because it emits Unknown rather than guessing.

Nothing returned here is a polarity. It is a proposal, and the needs_review column is TRUE for every row.

Value

A data.frame of drafted rows with w_g, w_g_source, evidence, citation and needs_review, plus an attr(, "unreachable") listing resources that did not respond.

Examples

## Not run: 
## Needs network access: each gene is one HTTP call to the chosen resources.
sel <- dmsa_select(systems = "hpa")
draft <- dmsa_polarity_fetch(sel, anchors = list("2" = c("CRH", "POMC")))
utils::write.csv(draft, file.path(tempdir(), "polarity_draft.csv"),
                 row.names = FALSE)

## End(Not run)

Description

Polarity table for one system, in the form dmsa_align() wants

Usage

dmsa_polarity_for(reference, system_id)

Arguments

Value

data.frame with gene and w_g, or NULL if the bundle carries no polarity.

Examples

f <- tempfile(fileext = ".csv")
dmsa_reference_template(f)
ref <- dmsa_reference_csv(f, quiet = TRUE)
## w_g is the gene's sign toward its system's activation tone: +1 raises it,
## -1 opposes it. dmsa_align() multiplies w_g by each CpG's
## methylation-to-expression direction to get the probe's aligned sign.
dmsa_polarity_for(ref, "HPA axis")

Description

Three kinds of row should never be used without someone looking at them: a sign the curator and a database disagree about, a sign no evidence could settle, and a non-zero sign held at low confidence. This returns them.

Usage

dmsa_polarity_review(
  x = "alpha",
  which = c("all", "disagreement", "unresolved", "low_confidence")
)