+
+
+[](https://github.com/BacZemin/BREAD/actions/workflows/R-CMD-check.yaml)
+[](https://github.com/BacZemin/BREAD/actions/workflows/bioc-check.yaml)
+[](https://github.com/BacZemin/BREAD/actions/workflows/pkgdown.yaml)
+[](https://lifecycle.r-lib.org/articles/stages.html#experimental)
+[](https://opensource.org/licenses/MIT)
+[](https://baczemin.github.io/BREAD/)
+
+
+**BREAD** (Bayesian Region-specific DNA methylation inference) provides
+targeted Bayesian inference for **predefined** DNA methylation regions from
+array data — a `SummarizedExperiment`, or a plain probe-by-sample matrix as
+`sesame::openSesame()` returns. For each region you supply, BREAD fits a
+Bayesian model, computes posterior probabilities of methylation change, and
+classifies the region as **hypermethylated**, **hypomethylated**,
+**unchanged**, or **inconclusive** at user-configurable effect-size and
+probability thresholds.
+
+Unlike genome-wide DMR callers that scan for regions, BREAD answers a
+different question: *given regions I already care about (PRC2 targets, CGIs,
+LADs, a chromHMM state, a custom BED), what is the posterior evidence for
+methylation change in each one, and how confident am I?* Output is a
+per-region posterior — effect size, credible interval, and probabilities —
+not just a p-value.
+
+### The `unchanged` class
+
+Most tools give you two states: significant, and not-significant. That
+conflates *flat* with *underpowered*, which at the sample sizes typical of
+experimental epigenetics is most of your genome.
+
+BREAD reports `prob_rope`, the posterior mass inside the region of practical
+equivalence spanning −`delta` to +`delta`, and calls a region **`unchanged`** when
+that mass clears `rope_cutoff`. `inconclusive` then means only what its name
+says. Being able to state that a region *demonstrably did not move* — with a
+credible ceiling on how much it could have — is the one claim a p-value is
+structurally unable to make, and it is often the claim a pathway or cascade
+argument actually needs.
+
+## Installation
+
+BREAD is in development. Install the latest version from GitHub:
+
+``` r
+# install.packages("pak")
+pak::pak("BacZemin/BREAD")
+
+# or with remotes:
+# install.packages("remotes")
+remotes::install_github("BacZemin/BREAD")
+```
+
+BREAD depends on Bioconductor packages (`SummarizedExperiment`,
+`GenomicRanges`, `S4Vectors`). The optional `brms` backend additionally needs
+`brms` plus a working Stan toolchain.
+
+## Example
+
+BREAD ships a small packaged dataset so you can run the whole pipeline out of
+the box: 8 EPICv2 arrays from a fibroblast passage-aging × vitamin C
+experiment, plus 500 predefined regions spanning five feature classes
+(PMD, PRC-CGI, bivalent, ...).
+
+```{r example, message = FALSE}
+library(BREAD)
+suppressPackageStartupMessages(library(SummarizedExperiment))
+
+se <- readRDS(system.file("extdata", "vitc_ag06561.rds", package = "BREAD"))
+reg <- readRDS(system.file("extdata", "vitc_regions.rds", package = "BREAD"))
+
+# Which regions change methylation with passage in the control fibroblasts?
+se_ctrl <- se[, se$condition == "ctrl"]
+
+fit <- fit_bread(
+ se_ctrl,
+ features = reg,
+ design = ~ passage,
+ feature_class_col = "feature_class"
+)
+
+fit # summary: n regions, backend, classification counts
+```
+
+```{r example-results}
+res <- results(fit) # one row per region
+table(res$classification)
+```
+
+That is the whole pattern: a `SummarizedExperiment` of array data, a
+`GRanges` of regions to test, and a model formula referencing columns of
+`colData(se)`. `fit_bread()` auto-detects the assay and whether values are
+on the beta or M scale, so the three arguments above are usually all you
+need. `feature_class_col` is optional — supply it when your regions carry a
+grouping column you want summarized.
+
+### Reading the results
+
+`results(fit)` returns one row per region:
+
+| column | meaning |
+|---|---|
+| `region_id` | region identifier |
+| `n` | samples contributing to the fit (non-NA); probe counts per region live in `fit@mapping$n_probes` |
+| `mean_effect` | posterior mean methylation change (M-scale) |
+| `ci_lo`, `ci_hi` | 95% credible interval |
+| `prob_hyper` | P(effect > +delta) — evidence for hypermethylation |
+| `prob_hypo` | P(effect < -delta) — evidence for hypomethylation |
+| `prob_rope` | P(\|effect\| <= delta) — evidence of *no* change |
+| `mean_dbeta`, `dbeta_lo`, `dbeta_hi` | the effect and interval on the beta scale |
+| `ref_beta`, `delta_beta` | the anchor used for that translation, and `delta` in beta units |
+| `classification` | `hypermethylated` / `hypomethylated` / `unchanged` / `inconclusive` |
+
+```{r results-head}
+head(res[, c("region_id", "n", "mean_effect", "ci_lo", "ci_hi",
+ "prob_rope", "classification")])
+```
+
+The rule: **hypermethylated** if `prob_hyper >= prob_cutoff`,
+**hypomethylated** if `prob_hypo >= prob_cutoff`, **unchanged** if
+`prob_rope >= rope_cutoff`, otherwise **inconclusive**. The three
+probabilities partition the posterior, so at any sensible cutoff exactly one
+can apply. Defaults are `delta = 0.10` (M-scale), `prob_cutoff = 0.95`, and
+`rope_cutoff = prob_cutoff`.
+
+`rope_cutoff` is separate for a reason, and the example above demonstrates it:
+**zero** regions are called `unchanged`. Concluding equivalence requires the
+whole posterior to fit inside ±`delta`, which is far stricter than a
+directional call, and this packaged dataset has too few samples to certify
+anything as flat:
+
+```{r rope-distribution}
+summary(res$prob_rope)
+```
+
+Those 385 `inconclusive` regions are genuinely *unresolved* — BREAD can
+neither detect a `delta = 0.10` effect nor exclude one. That is a more
+specific statement than a non-significant q-value, which looks the same
+whether the region is flat or the study is underpowered. A shared cutoff
+would make the equivalence class unreachable exactly where it is most
+wanted, so relax `rope_cutoff` on its own when you want to see the
+gradient — without loosening the bar for discovery.
+
+`classifications(fit)` returns just the per-region calls, and
+`posterior_draws(fit)` gives posterior samples for downstream summaries.
+
+### `delta` is on the M scale
+
+`delta = 0.10` in M-units is a beta change of about **0.017** at
+mid-methylation, and *less* toward the extremes (about 0.006 at β = 0.1) —
+the local slope dβ/dM = β(1−β)·ln2 is not constant. So a single beta-scale
+threshold does not exist; BREAD reports a per-region `delta_beta` instead,
+and `bread_delta_beta()` / `bread_delta_m()` convert explicitly when you want
+to choose `delta` from a target Δβ.
+
+```{r delta-scale}
+bread_delta_beta(0.10) # at beta = 0.5
+bread_delta_beta(0.10, ref_beta = c(0.5, 0.2, 0.1))
+bread_delta_m(0.02) # delta for a 2-point change
+```
+
+### Re-thresholding and permutation nulls
+
+`refit_bread()` reuses a fit's region-by-sample matrix. Changing only
+thresholds skips the model fit entirely; supplying new `colData` re-fits
+without recomputing the probe-to-region mapping, which is what makes
+label-permutation calibration practical.
+
+```{r refit, eval = FALSE}
+# Sweep the equivalence bar for free
+table(results(refit_bread(fit, rope_cutoff = 0.80))$classification)
+
+# One draw from a label-permutation null
+cd <- as.data.frame(colData(se_ctrl))
+cd$passage <- sample(cd$passage)
+results(refit_bread(fit, colData = cd))$prob_hyper
+```
+
+### Matrix input
+
+If your pipeline hands you a matrix rather than a `SummarizedExperiment` —
+`sesame::openSesame()` does — pass it directly with `colData` and either
+`rowRanges` or a `platform` to look the manifest up through `sesameData`:
+
+```{r matrix-input, eval = FALSE}
+fit_bread(betas, features = reg, design = ~ condition,
+ colData = pheno, platform = "EPICv2")
+```
+
+The platform is never inferred from probe IDs: `cg`-numbers are shared across
+HM450, EPIC and MM285, so a wrong guess would give wrong coordinates silently.
+
+### Backends
+
+- `backend = "conjugate"` (default) — analytic Normal-Inverse-Gamma
+ posterior, no MCMC. Hundreds of regions fit in well under a second.
+- `backend = "brms"` — full MCMC via Stan; compiles once, then reuses the
+ compiled model across regions. Use when you need the flexibility of a
+ full Bayesian fit.
+
+### KnowYourCG enrichment
+
+`bread_kycg()` takes the probes in your hyper- or hypo-classified regions
+and runs `knowYourCG::testEnrichment()` against curated CpG databases, so
+you can ask what genomic features your called regions are enriched for.
+
+## Vignettes
+
+Two worked examples on real data (rendered on the
+[documentation site](https://baczemin.github.io/BREAD/)):
+
+- **Getting started** (`bread-intro`) — TCGA HM450 matched normal/tumour
+ pairs over chromHMM chromatin-state regions, end-to-end through
+ `bread_kycg()`.
+- **Vitamin C EPICv2** (`bread-vitc`) — the packaged fibroblast
+ passage-aging × vitamin C experiment, recovering classic PMD-hypo /
+ PRC-CGI-hyper region signatures.
+
+## Status
+
+Milestone 1 (MVP) is complete: the full `fit_bread()` pipeline, conjugate
+and brms backends, S4 classes with accessors, plotting helpers, KYCG
+integration, two real-data vignettes, and a live pkgdown site. Partial
+pooling across regions, feature-class priors, random-effects designs, and
+contrast vectors are on the roadmap. See `NEWS.md` for the changelog.
+
+## Citation
+
+BREAD does not have an associated publication yet. For now, cite the
+package itself:
+
+```{r citation, comment = ""}
+citation("BREAD")
+```
+
+## License
+
+MIT © Jaemin Park. See `LICENSE`.
diff --git a/README.md b/README.md
index 583b95b..28b278e 100644
--- a/README.md
+++ b/README.md
@@ -1,27 +1,51 @@
-# BREAD
+
+
+
+# BREAD
-[](https://lifecycle.r-lib.org/articles/stages.html#experimental)
-[](https://opensource.org/licenses/MIT)
+
+[](https://github.com/BacZemin/BREAD/actions/workflows/R-CMD-check.yaml)
+[](https://github.com/BacZemin/BREAD/actions/workflows/bioc-check.yaml)
[](https://github.com/BacZemin/BREAD/actions/workflows/pkgdown.yaml)
+[](https://lifecycle.r-lib.org/articles/stages.html#experimental)
+[](https://opensource.org/licenses/MIT)
[](https://baczemin.github.io/BREAD/)
**BREAD** (Bayesian Region-specific DNA methylation inference) provides
-targeted Bayesian inference for **predefined** DNA methylation regions from
-array data stored in `SummarizedExperiment` objects. For each region you
-supply, BREAD fits a Bayesian model, computes the posterior probability of a
-directional methylation change, and classifies the region as
-**hypermethylated**, **hypomethylated**, or **inconclusive** at
+targeted Bayesian inference for **predefined** DNA methylation regions
+from array data — a `SummarizedExperiment`, or a plain probe-by-sample
+matrix as `sesame::openSesame()` returns. For each region you supply,
+BREAD fits a Bayesian model, computes posterior probabilities of
+methylation change, and classifies the region as **hypermethylated**,
+**hypomethylated**, **unchanged**, or **inconclusive** at
user-configurable effect-size and probability thresholds.
Unlike genome-wide DMR callers that scan for regions, BREAD answers a
-different question: *given regions I already care about (PRC2 targets, CGIs,
-LADs, a chromHMM state, a custom BED), what is the posterior evidence for
-methylation change in each one, and how confident am I?* Output is a
-per-region posterior — effect size, credible interval, and directional
+different question: *given regions I already care about (PRC2 targets,
+CGIs, LADs, a chromHMM state, a custom BED), what is the posterior
+evidence for methylation change in each one, and how confident am I?*
+Output is a per-region posterior — effect size, credible interval, and
probabilities — not just a p-value.
+### The `unchanged` class
+
+Most tools give you two states: significant, and not-significant. That
+conflates *flat* with *underpowered*, which at the sample sizes typical
+of experimental epigenetics is most of your genome.
+
+BREAD reports `prob_rope`, the posterior mass inside the region of
+practical equivalence spanning −`delta` to +`delta`, and calls a region
+**`unchanged`** when that mass clears `rope_cutoff`. `inconclusive` then
+means only what its name says. Being able to state that a region
+*demonstrably did not move* — with a credible ceiling on how much it
+could have — is the one claim a p-value is structurally unable to make,
+and it is often the claim a pathway or cascade argument actually needs.
+
## Installation
BREAD is in development. Install the latest version from GitHub:
@@ -36,18 +60,19 @@ remotes::install_github("BacZemin/BREAD")
```
BREAD depends on Bioconductor packages (`SummarizedExperiment`,
-`GenomicRanges`, `S4Vectors`, `GenomeInfoDb`). The optional `brms` backend
-additionally needs `brms` + a working Stan toolchain.
+`GenomicRanges`, `S4Vectors`). The optional `brms` backend additionally
+needs `brms` plus a working Stan toolchain.
## Example
-BREAD ships a small packaged dataset so you can run the whole pipeline out of
-the box: 8 EPICv2 arrays from a fibroblast passage-aging × vitamin C
-experiment, plus 500 predefined regions spanning five feature classes
-(PMD, PRC-CGI, bivalent, ...).
+BREAD ships a small packaged dataset so you can run the whole pipeline
+out of the box: 8 EPICv2 arrays from a fibroblast passage-aging ×
+vitamin C experiment, plus 500 predefined regions spanning five feature
+classes (PMD, PRC-CGI, bivalent, …).
``` r
library(BREAD)
+suppressPackageStartupMessages(library(SummarizedExperiment))
se <- readRDS(system.file("extdata", "vitc_ag06561.rds", package = "BREAD"))
reg <- readRDS(system.file("extdata", "vitc_regions.rds", package = "BREAD"))
@@ -61,61 +86,180 @@ fit <- fit_bread(
design = ~ passage,
feature_class_col = "feature_class"
)
+#> Warning: 496 of 496 fitted region(s) have fewer than 3 residual degrees of
+#> freedom (n - p < 3). Intervals there are optimistic; consider df_mode =
+#> "residual".
fit # summary: n regions, backend, classification counts
+#> }}\preformatted{nulls <- vapply(permutations, function(g) \{
+ cd <- coldata; cd$genotype <- g
+ results(refit_bread(fit, colData = cd))$prob_hyper[i]
+\}, numeric(1))
+}\if{html}{\out{
}}
+}
+
+\section{Re-thresholding is free}{
+
+When only \code{delta}, \code{prob_cutoff}, \code{rope_cutoff}, \code{ci} or \code{ref_beta} change,
+the model is not re-fitted at all — the stored posterior is re-summarized
+and re-classified. So sweeping a delta x cutoff grid costs essentially
+nothing.
+}
+
+\examples{
+suppressPackageStartupMessages(library(SummarizedExperiment))
+
+se <- readRDS(system.file("extdata", "vitc_ag06561.rds", package = "BREAD"))
+reg <- readRDS(system.file("extdata", "vitc_regions.rds", package = "BREAD"))
+se_ctrl <- se[, se$condition == "ctrl"]
+
+fit <- fit_bread(se_ctrl, reg, ~ passage)
+
+# Re-threshold without re-fitting anything
+table(results(refit_bread(fit, delta = 0.25))$classification)
+
+# Relax only the equivalence bar
+table(results(refit_bread(fit, rope_cutoff = 0.80))$classification)
+
+# Re-fit against shuffled labels (one draw from a permutation null)
+cd <- as.data.frame(colData(se_ctrl))
+cd$passage <- sample(cd$passage)
+head(results(refit_bread(fit, colData = cd))$prob_hyper)
+}
+\seealso{
+\code{\link[=fit_bread]{fit_bread()}}, \code{\link[=posterior_summary]{posterior_summary()}}, \code{\link[=classify_regions]{classify_regions()}}
+}
diff --git a/man/report_feature_set.Rd b/man/report_feature_set.Rd
deleted file mode 100644
index a93b356..0000000
--- a/man/report_feature_set.Rd
+++ /dev/null
@@ -1,21 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/report.R
-\name{report_feature_set}
-\alias{report_feature_set}
-\title{Feature-set level summary report}
-\usage{
-report_feature_set(fit, feature_class_col = NULL)
-}
-\arguments{
-\item{fit}{A \link{BreadFit}.}
-
-\item{feature_class_col}{Column in \code{mcols(features)} defining feature class.}
-}
-\value{
-A data frame with one row per feature class.
-}
-\description{
-Aggregates region-level classifications into a feature-set / feature-class
-summary (counts and proportions of hyper / hypo / inconclusive).
-}
-\keyword{internal}
diff --git a/man/results.Rd b/man/results.Rd
index 7b1d8b3..f5d1f62 100644
--- a/man/results.Rd
+++ b/man/results.Rd
@@ -20,3 +20,16 @@ A data frame with one row per region.
\description{
Extract the region-level results table from a \link{BreadFit}
}
+\examples{
+suppressPackageStartupMessages(library(SummarizedExperiment))
+
+se <- readRDS(system.file("extdata", "vitc_ag06561.rds", package = "BREAD"))
+reg <- readRDS(system.file("extdata", "vitc_regions.rds", package = "BREAD"))
+se_ctrl <- se[, se$condition == "ctrl"]
+
+fit <- fit_bread(se_ctrl, reg, ~ passage)
+
+res <- results(fit)
+head(res)
+colnames(res)
+}
diff --git a/man/summarize_features.Rd b/man/summarize_features.Rd
index fb7df65..957c2ee 100644
--- a/man/summarize_features.Rd
+++ b/man/summarize_features.Rd
@@ -50,3 +50,19 @@ M-value interpretation under \code{"pc1"}.
}
}
+\examples{
+suppressPackageStartupMessages(library(SummarizedExperiment))
+
+se <- readRDS(system.file("extdata", "vitc_ag06561.rds", package = "BREAD"))
+reg <- readRDS(system.file("extdata", "vitc_regions.rds", package = "BREAD"))
+se_ctrl <- se[, se$condition == "ctrl"]
+
+mapping <- map_probes_to_features(se_ctrl, reg)
+
+# As above, the packaged assay is "betas" on the beta scale; values
+# are converted to M-values internally before summarizing.
+mat <- summarize_features(se_ctrl, mapping,
+ assay_name = "betas", input_scale = "Beta")
+dim(mat)
+mat[1:3, 1:3]
+}
diff --git a/man/validate_bread_input.Rd b/man/validate_bread_input.Rd
index a38e12b..6c49c7d 100644
--- a/man/validate_bread_input.Rd
+++ b/man/validate_bread_input.Rd
@@ -37,3 +37,16 @@ coordinates, that \code{features} is a non-empty \link[GenomicRanges:GRanges-cla
that \code{assay_name} is present, and that \code{contrast} (when non-NULL) resolves
to a coefficient of the design's model matrix.
}
+\examples{
+suppressPackageStartupMessages(library(SummarizedExperiment))
+
+se <- readRDS(system.file("extdata", "vitc_ag06561.rds", package = "BREAD"))
+reg <- readRDS(system.file("extdata", "vitc_regions.rds", package = "BREAD"))
+se_ctrl <- se[, se$condition == "ctrl"]
+
+# The packaged data keeps beta values in an assay named "betas", so
+# both arguments are given explicitly here. fit_bread() detects them
+# for you; this lower-level helper does not.
+validate_bread_input(se_ctrl, reg, ~ passage,
+ assay_name = "betas", input_scale = "Beta")
+}
diff --git a/tests/testthat/helper-toy.R b/tests/testthat/helper-toy.R
index 6dd70b8..fbd69b6 100644
--- a/tests/testthat/helper-toy.R
+++ b/tests/testthat/helper-toy.R
@@ -15,7 +15,10 @@
cd <- S4Vectors::DataFrame(
group = factor(rep(c("young", "old"), length.out = n_samples),
levels = c("young", "old")),
- sex = factor(rep(c("F", "M"), length.out = n_samples),
+ # Period 4 against group's period 2, so `~ group + sex` is full rank.
+ # (When both alternated, the two were perfectly collinear and any design
+ # using both silently fell back on the prior.)
+ sex = factor(rep(c("F", "F", "M", "M"), length.out = n_samples),
levels = c("F", "M")),
row.names = colnames(m)
)
@@ -41,6 +44,24 @@
gr
}
+# Toy features where ONE region is defined by several disjoint ranges sharing
+# a name -- the only way to pin a region to an exact probe set, since a single
+# bounding interval would sweep in the probes between them.
+# regD = probes 1-3 (1..2500) + probes 8-10 (7001..9500) = 6 probes, 2 ranges.
+# regE = probes 16-20, 1 range. So: 2 distinct regions from 3 ranges.
+.make_toy_features_dup <- function() {
+ gr <- GenomicRanges::GRanges(
+ seqnames = "chr1",
+ ranges = IRanges::IRanges(
+ start = c( 1L, 7001L, 15001L),
+ end = c( 2500L, 9500L, 20000L)
+ ),
+ feature_class = c("PRC", "PRC", "LAD")
+ )
+ names(gr) <- c("regD", "regD", "regE")
+ gr
+}
+
# Toy SE with injected signal: regA becomes hyper, regC becomes hypo,
# in "old" vs "young" under ~ group. Sample size large enough for
# prob_cutoff = 0.95 classification to recover truth.
diff --git a/tests/testthat/test-classes.R b/tests/testthat/test-classes.R
new file mode 100644
index 0000000..10ddad1
--- /dev/null
+++ b/tests/testthat/test-classes.R
@@ -0,0 +1,57 @@
+# The S4 slot layout is a public contract: the plotting helpers and the
+# accessors reach into these slots by name. Pin them.
+
+test_that("BreadFit slots are the documented set, in order", {
+ expect_identical(
+ methods::slotNames("BreadFit"),
+ c("call", "params", "mode", "assay_name", "input_scale",
+ "mapping", "features", "model", "posterior", "results", "diagnostics")
+ )
+})
+
+test_that("BreadResults slots are the documented set", {
+ expect_identical(methods::slotNames("BreadResults"), c("table", "params"))
+})
+
+test_that("a fitted BreadFit populates every slot it promises", {
+ fit <- fit_bread(.make_toy_signal_se(), .make_toy_features(), ~ group)
+ expect_s4_class(fit, "BreadFit")
+ expect_true(is.call(fit@call))
+ expect_identical(fit@mode, "summary")
+ expect_type(fit@params, "list")
+ expect_type(fit@diagnostics, "list")
+ expect_s3_class(fit@mapping, "data.frame")
+ expect_s3_class(fit@results, "data.frame")
+ expect_s3_class(fit@posterior, "data.frame")
+ expect_type(fit@model, "list")
+ expect_identical(fit@diagnostics$backend, "conjugate")
+})
+
+test_that("BreadResults() round-trips the results table", {
+ fit <- fit_bread(.make_toy_signal_se(), .make_toy_features(), ~ group)
+ br <- BreadResults(fit)
+ expect_s4_class(br, "BreadResults")
+ expect_equal(methods::slot(br, "table"), results(fit))
+ expect_equal(methods::slot(br, "params"), fit@params)
+})
+
+test_that("BreadResults() rejects non-BreadFit input", {
+ expect_error(BreadResults(data.frame(x = 1)), "must be a BreadFit")
+ expect_error(BreadResults(NULL), "must be a BreadFit")
+})
+
+test_that("the model list shape is a stable contract", {
+ # refit_bread(), posterior_summary() and plot_region_data() all reach into
+ # these names. Pin them so a backend refactor cannot quietly break them.
+ fit <- suppressMessages(
+ fit_bread(.make_toy_signal_se(), .make_toy_features(), ~ group)
+ )
+ expect_identical(
+ names(fit@model),
+ c("fits", "design_matrix", "coef_names", "contrast", "contrast_idx",
+ "region_ids", "prior", "df_mode", "region_mat", "design", "coldata")
+ )
+ expect_identical(fit@model$df_mode, "conjugate")
+ expect_true(is.matrix(fit@model$region_mat))
+ expect_identical(rownames(fit@model$region_mat), fit@model$region_ids)
+})
diff --git a/tests/testthat/test-classify.R b/tests/testthat/test-classify.R
new file mode 100644
index 0000000..9176389
--- /dev/null
+++ b/tests/testthat/test-classify.R
@@ -0,0 +1,175 @@
+# classify_regions() applies the decision rule; it does not recompute
+# posterior probabilities. These tests drive it with hand-built posterior
+# tables so the boundaries are exact rather than approximate.
+#
+# prob_hyper, prob_hypo and prob_rope partition the posterior, so fixtures
+# must sum to 1 -- an incoherent triple can exercise states the sampler can
+# never reach and would let a real bug hide.
+
+`%|NULL|%` <- function(a, b) if (is.null(a)) b else a
+
+.fake_post <- function(p_gt, p_lt, p_rope = NULL, ids = NULL) {
+ n <- length(p_gt)
+ out <- data.frame(
+ region_id = ids %|NULL|% sprintf("r%02d", seq_len(n)),
+ prob_hyper = p_gt,
+ prob_hypo = p_lt,
+ stringsAsFactors = FALSE
+ )
+ if (!is.null(p_rope)) out$prob_rope <- p_rope
+ out
+}
+
+# Uniform draws from the 2-simplex (Dirichlet(1,1,1)).
+.coherent_post <- function(n, seed = 42L) {
+ withr::with_seed(seed, {
+ e <- matrix(stats::rexp(3L * n), ncol = 3L)
+ p <- e / rowSums(e)
+ })
+ .fake_post(p[, 1L], p[, 2L], p[, 3L])
+}
+
+
+test_that("classification levels are exact and in a fixed order", {
+ out <- classify_regions(
+ .fake_post(c(0.99, 0.01, 0.02, 0.34),
+ c(0.01, 0.99, 0.02, 0.33),
+ c(0.00, 0.00, 0.96, 0.33))
+ )
+ expect_s3_class(out$classification, "factor")
+ expect_identical(
+ levels(out$classification),
+ c("hypermethylated", "hypomethylated", "unchanged", "inconclusive")
+ )
+ expect_identical(
+ as.character(out$classification),
+ c("hypermethylated", "hypomethylated", "unchanged", "inconclusive")
+ )
+})
+
+test_that("the motivating case is called unchanged, not inconclusive", {
+ # REGULON_promoter from the mousearray_609G CArG analysis: the posterior
+ # sits almost entirely inside the ROPE. Before the fourth level existed
+ # this was labelled `inconclusive`, indistinguishable from a region whose
+ # posterior spanned the whole line.
+ out <- classify_regions(.fake_post(0.021, 0.021, 0.958), prob_cutoff = 0.95)
+ expect_identical(as.character(out$classification), "unchanged")
+})
+
+test_that("the cutoff comparison is inclusive (>=), not strict", {
+ at <- classify_regions(.fake_post(0.95, 0.0, 0.05), prob_cutoff = 0.95)
+ below <- classify_regions(.fake_post(0.95 - 1e-12, 0.0, 0.05),
+ prob_cutoff = 0.95)
+ expect_identical(as.character(at$classification), "hypermethylated")
+ expect_identical(as.character(below$classification), "inconclusive")
+
+ rope_at <- classify_regions(.fake_post(0.02, 0.03, 0.95),
+ rope_cutoff = 0.95)
+ rope_below <- classify_regions(.fake_post(0.02, 0.03, 0.95 - 1e-12),
+ rope_cutoff = 0.95)
+ expect_identical(as.character(rope_at$classification), "unchanged")
+ expect_identical(as.character(rope_below$classification), "inconclusive")
+})
+
+test_that("at sane cutoffs the four classes are mutually exclusive", {
+ post <- .coherent_post(10000L)
+ cls <- classify_regions(post, prob_cutoff = 0.95,
+ rope_cutoff = 0.95)$classification
+ n_clearing <- (post$prob_hyper >= 0.95) +
+ (post$prob_hypo >= 0.95) +
+ (post$prob_rope >= 0.95)
+ expect_true(all(n_clearing <= 1L))
+ expect_false(any(is.na(cls)))
+ # Every row that clears something is labelled, and nothing else is
+ expect_identical(cls != "inconclusive", n_clearing == 1L)
+})
+
+test_that("tightening both bars only ever moves regions to inconclusive", {
+ post <- .coherent_post(500L)
+ loose <- classify_regions(post, prob_cutoff = 0.60)$classification
+ strict <- classify_regions(post, prob_cutoff = 0.99)$classification
+ moved <- loose != strict
+ expect_true(all(strict[moved] == "inconclusive"))
+})
+
+test_that("rope_cutoff moves only the equivalence bar", {
+ post <- .coherent_post(500L)
+ base <- classify_regions(post, prob_cutoff = 0.95)$classification
+ looser <- classify_regions(post, prob_cutoff = 0.95,
+ rope_cutoff = 0.70)$classification
+ # Directional calls are untouched ...
+ dir_base <- base %in% c("hypermethylated", "hypomethylated")
+ dir_loose<- looser %in% c("hypermethylated", "hypomethylated")
+ expect_identical(dir_base, dir_loose)
+ # ... and the only movement is inconclusive -> unchanged
+ moved <- base != looser
+ expect_true(all(base[moved] == "inconclusive"))
+ expect_true(all(looser[moved] == "unchanged"))
+ expect_true(any(moved))
+
+ # It defaults to prob_cutoff
+ expect_identical(
+ classify_regions(post, prob_cutoff = 0.8)$classification,
+ classify_regions(post, prob_cutoff = 0.8, rope_cutoff = 0.8)$classification
+ )
+})
+
+test_that("failed fits (NA probabilities) become inconclusive, never NA", {
+ out <- classify_regions(.fake_post(c(NA, 0.99), c(NA, 0.001)))
+ expect_false(any(is.na(out$classification)))
+ expect_identical(as.character(out$classification)[1], "inconclusive")
+})
+
+test_that("an NA posterior is never called unchanged", {
+ # A missing posterior is absence of evidence, not evidence of equivalence.
+ out <- classify_regions(.fake_post(NA_real_, NA_real_, 1.0))
+ expect_identical(as.character(out$classification), "inconclusive")
+})
+
+test_that("when several classes clear a low cutoff the largest wins", {
+ out <- classify_regions(
+ .fake_post(c(0.40, 0.31, 0.31), c(0.31, 0.40, 0.31), c(0.29, 0.29, 0.38)),
+ prob_cutoff = 0.30, rope_cutoff = 0.30
+ )
+ expect_identical(as.character(out$classification),
+ c("hypermethylated", "hypomethylated", "unchanged"))
+
+ # Exact ties resolve hyper > hypo > unchanged, as before the fourth level
+ tie <- classify_regions(.fake_post(1/3, 1/3, 1/3), prob_cutoff = 0.30)
+ expect_identical(as.character(tie$classification), "hypermethylated")
+})
+
+test_that("prob_rope is derived when absent and clamped when incoherent", {
+ # Hand-built two-column tables remain a legitimate input
+ derived <- classify_regions(.fake_post(0.01, 0.01), rope_cutoff = 0.95)
+ expect_identical(as.character(derived$classification), "unchanged")
+
+ # hyper + hypo > 1 cannot happen from a real posterior, but must not
+ # produce a negative ROPE mass or an NA class
+ bad <- classify_regions(.fake_post(0.8, 0.8), prob_cutoff = 0.95)
+ expect_identical(as.character(bad$classification), "inconclusive")
+
+ # A supplied prob_rope is authoritative over the complement
+ supplied <- classify_regions(.fake_post(0.01, 0.01, 0.10),
+ rope_cutoff = 0.95)
+ expect_identical(as.character(supplied$classification), "inconclusive")
+})
+
+test_that("delta and both cutoffs are recorded as attributes", {
+ out <- classify_regions(.fake_post(0.99, 0.0), delta = 0.25,
+ prob_cutoff = 0.9, rope_cutoff = 0.7)
+ expect_equal(attr(out, "delta"), 0.25)
+ expect_equal(attr(out, "prob_cutoff"), 0.9)
+ expect_equal(attr(out, "rope_cutoff"), 0.7)
+})
+
+test_that("bad input is rejected with an informative error", {
+ expect_error(classify_regions("nope"), "data.frame")
+ expect_error(classify_regions(data.frame(x = 1)), "missing required columns")
+ expect_error(classify_regions(.fake_post(0.9, 0.1), prob_cutoff = 1),
+ "`prob_cutoff` must be in \\(0, 1\\)")
+ expect_error(classify_regions(.fake_post(0.9, 0.1), prob_cutoff = 0),
+ "`prob_cutoff` must be in \\(0, 1\\)")
+ expect_error(classify_regions(.fake_post(0.9, 0.1), rope_cutoff = 1),
+ "`rope_cutoff` must be in \\(0, 1\\)")
+})
diff --git a/tests/testthat/test-coerce.R b/tests/testthat/test-coerce.R
new file mode 100644
index 0000000..b8b56cd
--- /dev/null
+++ b/tests/testthat/test-coerce.R
@@ -0,0 +1,187 @@
+# Matrix input has to reach exactly the same answer as SummarizedExperiment
+# input -- otherwise the sesame workflow is a second, subtly different code
+# path rather than a convenience.
+
+.toy_parts <- function() {
+ se <- .make_toy_signal_se()
+ list(
+ se = se,
+ mat = SummarizedExperiment::assay(se, "M"),
+ cd = as.data.frame(SummarizedExperiment::colData(se)),
+ gr = SummarizedExperiment::rowRanges(se)
+ )
+}
+
+test_that("matrix input and SE input give identical results", {
+ p <- .toy_parts()
+ gr_feat <- .make_toy_features()
+
+ fit_se <- suppressMessages(fit_bread(p$se, gr_feat, ~ group))
+ fit_m <- suppressMessages(
+ fit_bread(p$mat, gr_feat, ~ group, colData = p$cd, rowRanges = p$gr)
+ )
+
+ expect_equal(results(fit_m), results(fit_se))
+ expect_equal(fit_m@posterior, fit_se@posterior)
+ expect_equal(fit_m@mapping, fit_se@mapping)
+ expect_identical(fit_m@input_scale, fit_se@input_scale)
+})
+
+test_that("bread_se() round-trips a decomposed SummarizedExperiment", {
+ p <- .toy_parts()
+ se2 <- bread_se(p$mat, colData = p$cd, rowRanges = p$gr)
+ expect_s4_class(se2, "SummarizedExperiment")
+ expect_equal(SummarizedExperiment::assay(se2, "M"), p$mat)
+ expect_equal(nrow(SummarizedExperiment::colData(se2)), ncol(p$mat))
+ expect_identical(names(SummarizedExperiment::rowRanges(se2)), rownames(p$mat))
+})
+
+test_that("a sesameData-style list is unwrapped", {
+ p <- .toy_parts()
+ se2 <- bread_se(list(betas = p$mat, sampleInfo = p$cd), rowRanges = p$gr)
+ expect_s4_class(se2, "SummarizedExperiment")
+ expect_true("group" %in% colnames(SummarizedExperiment::colData(se2)))
+})
+
+test_that("the assay is named from the value range, as fit_bread expects", {
+ p <- .toy_parts()
+ se_m <- bread_se(p$mat, colData = p$cd, rowRanges = p$gr)
+ expect_identical(SummarizedExperiment::assayNames(se_m), "M")
+
+ betas <- 2^p$mat / (2^p$mat + 1)
+ se_b <- bread_se(betas, colData = p$cd, rowRanges = p$gr)
+ expect_identical(SummarizedExperiment::assayNames(se_b), "betas")
+})
+
+
+# ---- rowRanges alignment ---------------------------------------------------
+
+test_that("a named manifest is subset and reordered to the matrix rows", {
+ p <- .toy_parts()
+ shuffled <- p$gr[sample(length(p$gr))]
+ extra <- p$gr[1:3]
+ names(extra) <- paste0("zz", 1:3)
+ manifest <- c(shuffled, extra) # longer, out of order
+
+ se2 <- bread_se(p$mat, colData = p$cd, rowRanges = manifest)
+ expect_identical(names(SummarizedExperiment::rowRanges(se2)), rownames(p$mat))
+ expect_equal(SummarizedExperiment::rowRanges(se2), p$gr)
+})
+
+test_that("probes absent from the manifest are dropped with a message", {
+ p <- .toy_parts()
+ partial <- p$gr[1:15] # 5 probes have no coordinates
+ expect_message(
+ se2 <- bread_se(p$mat, colData = p$cd, rowRanges = partial),
+ "Dropping 5 of 20 probes"
+ )
+ expect_equal(nrow(se2), 15L)
+})
+
+test_that("an entirely mismatched manifest errors rather than dropping all", {
+ p <- .toy_parts()
+ wrong <- p$gr
+ names(wrong) <- paste0("nope", seq_along(wrong))
+ expect_error(bread_se(p$mat, colData = p$cd, rowRanges = wrong),
+ "None of `rownames\\(x\\)`")
+})
+
+test_that("an unnamed rowRanges must match the row count exactly", {
+ p <- .toy_parts()
+ unnamed <- p$gr; names(unnamed) <- NULL
+ se2 <- bread_se(p$mat, colData = p$cd, rowRanges = unnamed)
+ expect_identical(names(SummarizedExperiment::rowRanges(se2)), rownames(p$mat))
+
+ expect_error(bread_se(p$mat, colData = p$cd, rowRanges = unnamed[1:5]),
+ "one range per row")
+})
+
+
+# ---- colData alignment -----------------------------------------------------
+
+test_that("colData is reordered by rownames, not trusted positionally", {
+ p <- .toy_parts()
+ shuffled <- p$cd[sample(nrow(p$cd)), , drop = FALSE]
+
+ ordered <- bread_se(p$mat, colData = p$cd, rowRanges = p$gr)
+ reorder <- bread_se(p$mat, colData = shuffled, rowRanges = p$gr)
+ expect_equal(SummarizedExperiment::colData(reorder),
+ SummarizedExperiment::colData(ordered))
+})
+
+test_that("colData that does not cover every sample errors", {
+ p <- .toy_parts()
+ expect_error(
+ bread_se(p$mat, colData = p$cd[1:4, , drop = FALSE], rowRanges = p$gr),
+ "no row for"
+ )
+})
+
+test_that("colData without rownames is accepted but warns", {
+ p <- .toy_parts()
+ bare <- p$cd; rownames(bare) <- NULL
+ expect_warning(bread_se(p$mat, colData = bare, rowRanges = p$gr),
+ "assuming its rows are in the same order")
+})
+
+
+# ---- refusals --------------------------------------------------------------
+
+test_that("coordinates are never guessed from probe IDs", {
+ p <- .toy_parts()
+ expect_error(bread_se(p$mat, colData = p$cd), "does not guess the platform")
+})
+
+test_that("matrix-only arguments are rejected alongside an SE", {
+ p <- .toy_parts()
+ expect_error(bread_se(p$se, colData = p$cd), "supplied alongside")
+ expect_error(bread_se(p$se, rowRanges = p$gr), "supplied alongside")
+ expect_error(bread_se(p$se, platform = "EPIC"),"supplied alongside")
+ expect_s4_class(bread_se(p$se), "SummarizedExperiment")
+})
+
+test_that("a matrix without dimnames is rejected", {
+ p <- .toy_parts()
+ m <- p$mat; rownames(m) <- NULL
+ expect_error(bread_se(m, colData = p$cd, rowRanges = p$gr), "rownames")
+
+ m2 <- p$mat; colnames(m2) <- NULL
+ expect_error(bread_se(m2, colData = p$cd, rowRanges = p$gr), "colnames")
+})
+
+test_that("unsupported input still names SummarizedExperiment in the error", {
+ # test-smoke.R relies on this string
+ expect_error(bread_se(NULL), "SummarizedExperiment")
+ expect_error(fit_bread(NULL, NULL, ~ 1), "SummarizedExperiment")
+})
+
+test_that("the platform route reaches sesameData", {
+ skip_on_cran()
+ skip_on_ci()
+ skip_if_not_installed("sesameData")
+ # Probe IDs are taken from the manifest itself, so this tests the lookup
+ # rather than whether some other dataset happens to share its ID
+ # convention. (The packaged vitc subset does not: it carries stripped
+ # EPICv2 IDs while the manifest keeps the replicate suffix.)
+ man <- sesameData::sesameData_getManifestGRanges("EPICv2")
+ skip_if(length(man) == 0L, "EPICv2 manifest unavailable offline")
+
+ ids <- names(man)[seq_len(50L)]
+ mat <- matrix(stats::runif(50L * 4L), nrow = 50L,
+ dimnames = list(ids, sprintf("S%d", 1:4)))
+ cd <- data.frame(group = rep(c("a", "b"), 2), row.names = colnames(mat))
+
+ se2 <- suppressMessages(bread_se(mat, colData = cd, platform = "EPICv2"))
+ expect_s4_class(se2, "SummarizedExperiment")
+ expect_equal(nrow(se2), 50L)
+ expect_identical(names(SummarizedExperiment::rowRanges(se2)), ids)
+ expect_identical(SummarizedExperiment::assayNames(se2), "betas")
+})
+
+test_that("stripped EPICv2 suffixes get a specific diagnosis", {
+ p <- .toy_parts()
+ manifest <- p$gr
+ names(manifest) <- paste0(names(manifest), "_BC11")
+ expect_error(bread_se(p$mat, colData = p$cd, rowRanges = manifest),
+ "EPICv2 replicate suffixes")
+})
diff --git a/tests/testthat/test-df-mode.R b/tests/testthat/test-df-mode.R
new file mode 100644
index 0000000..f061b7f
--- /dev/null
+++ b/tests/testthat/test-df-mode.R
@@ -0,0 +1,154 @@
+# Degrees-of-freedom handling: the n <= p guard and the "residual" mode.
+#
+# Motivation: `a_n = a0 + n/2` makes nu = 2*a_n a function of n alone, never of
+# p. Under the weak default prior that overstates precision by exactly
+# sqrt(n / (n - p)) on the posterior scale, and at n == p the residuals are
+# identically zero so the scale collapses to the prior floor.
+
+# --- helpers ---------------------------------------------------------------
+
+# Region x sample matrix with a known design; y is pure noise unless `beta` set.
+sim_mat <- function(n_per_cell, n_regions = 5L, sd = 1, beta = 0, seed = 1L) {
+ set.seed(seed)
+ g <- rep(c("a", "b"), each = 2L * n_per_cell)
+ t <- rep(rep(c("x", "y"), each = n_per_cell), 2L)
+ cd <- data.frame(g = factor(g), t = factor(t))
+ X <- stats::model.matrix(~ g * t, cd)
+ k <- which(colnames(X) == "gb:ty")
+ eta <- as.numeric(X[, k]) * beta
+ m <- matrix(rep(eta, each = n_regions) + stats::rnorm(n_regions * nrow(cd), sd = sd),
+ nrow = n_regions, dimnames = list(paste0("R", seq_len(n_regions)), rownames(cd)))
+ list(mat = m, cd = cd, p = ncol(X), contrast = "gb:ty")
+}
+
+fit_one <- function(s, df_mode = "conjugate", prior = NULL) {
+ BREAD:::fit_bread_summary(s$mat, s$cd, ~ g * t, s$contrast,
+ prior = prior, df_mode = df_mode)
+}
+
+# --- the n <= p guard ------------------------------------------------------
+
+test_that("regions with n == p are dropped rather than fitted", {
+ s <- sim_mat(n_per_cell = 1L) # n = 4, p = 4
+ expect_identical(s$p, 4L)
+ f <- suppressWarnings(fit_one(s))
+ errs <- vapply(f$fits, function(z) z$error, character(1))
+ expect_true(all(errs == "n <= number of coefficients"))
+ expect_true(all(vapply(f$fits, function(z) is.na(z$a_n), logical(1))))
+})
+
+test_that("the n <= p guard applies under both df_mode settings", {
+ s <- sim_mat(n_per_cell = 1L)
+ for (dm in c("conjugate", "residual")) {
+ f <- suppressWarnings(fit_one(s, df_mode = dm))
+ expect_true(all(vapply(f$fits, function(z) z$error, character(1)) ==
+ "n <= number of coefficients"))
+ }
+})
+
+test_that("n < 2 still reports the pre-existing reason", {
+ s <- sim_mat(n_per_cell = 2L)
+ s$mat[1, ] <- NA_real_
+ s$mat[1, 1] <- 0.5 # a single non-NA sample
+ f <- suppressWarnings(fit_one(s))
+ expect_identical(f$fits[[1]]$error, "too few non-NA samples")
+})
+
+test_that("n > p fits normally and carries no error", {
+ s <- sim_mat(n_per_cell = 3L) # n = 12, p = 4
+ f <- fit_one(s)
+ expect_true(all(is.na(vapply(f$fits, function(z) z$error, character(1)))))
+})
+
+# --- low residual df warns once, not per region ---------------------------
+
+test_that("fewer than 3 residual df warns exactly once for the whole fit", {
+ s <- sim_mat(n_per_cell = 2L, n_regions = 10L) # n = 8, p = 4 -> n - p = 4
+ expect_silent(fit_one(s))
+
+ s2 <- sim_mat(n_per_cell = 2L, n_regions = 10L)
+ s2$cd$z <- factor(rep(c("u", "v"), length.out = nrow(s2$cd)))
+ # ~ g * t + z -> p = 5, n = 8, n - p = 3 -> still silent
+ f5 <- BREAD:::fit_bread_summary(s2$mat, s2$cd, ~ g * t + z, "gb:ty")
+ expect_true(is.list(f5$fits))
+
+ s3 <- sim_mat(n_per_cell = 2L, n_regions = 10L)
+ s3$mat[, 1:3] <- NA_real_ # n drops to 5, p = 4 -> n - p = 1
+ w <- capture_warnings(fit_one(s3))
+ expect_length(w, 1L)
+ expect_match(w, "residual degrees of freedom", fixed = FALSE)
+})
+
+test_that("the warning names df_mode = residual only in conjugate mode", {
+ s <- sim_mat(n_per_cell = 2L, n_regions = 4L)
+ s$mat[, 1:3] <- NA_real_
+ expect_match(capture_warnings(fit_one(s, df_mode = "conjugate")),
+ "df_mode")
+ expect_false(any(grepl("df_mode",
+ capture_warnings(fit_one(s, df_mode = "residual")))))
+})
+
+# --- residual mode reproduces the classical answer ------------------------
+
+test_that("df_mode = 'residual' matches lm() df and standard error", {
+ s <- sim_mat(n_per_cell = 4L, n_regions = 3L, seed = 7L) # n = 16, p = 4
+ # weak prior so the posterior should collapse onto OLS
+ pr <- bread_prior(lambda0 = 1e-8, a0 = 1e-8, b0 = 1e-8)
+ f <- fit_one(s, df_mode = "residual", prior = pr)
+ k <- f$contrast_idx
+
+ for (i in seq_len(nrow(s$mat))) {
+ ml <- stats::lm(s$mat[i, ] ~ g * t, data = s$cd)
+ fo <- f$fits[[i]]
+ scale_b <- sqrt((fo$b_n / fo$a_n) * fo$Lambda_n_inv[k, k])
+ expect_equal(2 * fo$a_n, ml$df.residual, tolerance = 1e-5)
+ expect_equal(scale_b, summary(ml)$coefficients[k, 2], tolerance = 1e-4)
+ expect_equal(fo$mu_n[k], unname(coef(ml)[k]), tolerance = 1e-5)
+ }
+})
+
+test_that("conjugate vs residual differ by exactly sqrt(n / (n - p))", {
+ s <- sim_mat(n_per_cell = 4L, n_regions = 3L, seed = 11L) # n = 16, p = 4
+ pr <- bread_prior(lambda0 = 1e-8, a0 = 1e-8, b0 = 1e-8)
+ fc <- fit_one(s, df_mode = "conjugate", prior = pr)
+ fr <- fit_one(s, df_mode = "residual", prior = pr)
+ k <- fc$contrast_idx
+
+ sc <- vapply(fc$fits, function(z) sqrt((z$b_n / z$a_n) * z$Lambda_n_inv[k, k]), 0)
+ sr <- vapply(fr$fits, function(z) sqrt((z$b_n / z$a_n) * z$Lambda_n_inv[k, k]), 0)
+ n <- 16L; p <- 4L
+ expect_equal(unname(sr / sc), rep(sqrt(n / (n - p)), length(sc)), tolerance = 1e-5)
+})
+
+test_that("residual mode widens credible intervals", {
+ s <- sim_mat(n_per_cell = 4L, n_regions = 5L, seed = 3L)
+ pr <- bread_prior(lambda0 = 1e-8, a0 = 1e-8, b0 = 1e-8)
+ pc <- posterior_summary(fit_one(s, "conjugate", pr))
+ prr <- posterior_summary(fit_one(s, "residual", pr))
+ expect_true(all((prr$ci_hi - prr$ci_lo) > (pc$ci_hi - pc$ci_lo)))
+ expect_equal(prr$mean_effect, pc$mean_effect, tolerance = 1e-6)
+})
+
+# --- end-to-end through fit_bread() ---------------------------------------
+
+test_that("fit_bread() accepts df_mode and records it in params", {
+ skip_if_not_installed("SummarizedExperiment")
+ se <- .make_toy_signal_se(); gr <- .make_toy_features()
+ fit_c <- fit_bread(se, gr, ~ group)
+ fit_r <- fit_bread(se, gr, ~ group, df_mode = "residual")
+
+ expect_identical(fit_c@params$df_mode, "conjugate")
+ expect_identical(fit_r@params$df_mode, "residual")
+
+ rc <- results(fit_c); rr <- results(fit_r)
+ ok <- is.na(rc$error) & is.na(rr$error)
+ expect_true(any(ok))
+ expect_true(all(rr$df[ok] < rc$df[ok]))
+ expect_true(all(rr$scale[ok] >= rc$scale[ok]))
+ expect_equal(rr$mean_effect[ok], rc$mean_effect[ok], tolerance = 1e-8)
+})
+
+test_that("df_mode is rejected when misspelled", {
+ se <- .make_toy_signal_se(); gr <- .make_toy_features()
+ expect_error(fit_bread(se, gr, ~ group, df_mode = "residuals"))
+})
diff --git a/tests/testthat/test-extdata-contract.R b/tests/testthat/test-extdata-contract.R
new file mode 100644
index 0000000..f74db98
--- /dev/null
+++ b/tests/testthat/test-extdata-contract.R
@@ -0,0 +1,80 @@
+# Every @examples block in this package is written against the two packaged
+# extdata objects. If their shape drifts -- an assay rename, a dropped
+# colData column, a change in the feature classes -- the examples break at
+# R CMD check time with an opaque error. These tests pin the contract so the
+# failure lands here instead, with a message that says what changed.
+
+.se <- function() {
+ readRDS(system.file("extdata", "vitc_ag06561.rds", package = "BREAD"))
+}
+.reg <- function() {
+ readRDS(system.file("extdata", "vitc_regions.rds", package = "BREAD"))
+}
+
+test_that("both extdata files are installed and loadable", {
+ expect_true(nzchar(system.file("extdata", "vitc_ag06561.rds",
+ package = "BREAD")))
+ expect_true(nzchar(system.file("extdata", "vitc_regions.rds",
+ package = "BREAD")))
+})
+
+test_that("the packaged SE has the assay name and scale the examples assume", {
+ se <- .se()
+ expect_s4_class(se, "RangedSummarizedExperiment")
+ expect_identical(SummarizedExperiment::assayNames(se), "betas")
+
+ x <- SummarizedExperiment::assay(se, "betas")
+ expect_true(all(x >= 0 & x <= 1, na.rm = TRUE))
+ # Examples rely on fit_bread() auto-detecting both of these.
+ expect_identical(BREAD:::.detect_input_scale(x), "Beta")
+ expect_identical(BREAD:::.detect_assay_name(se), "betas")
+})
+
+test_that("the colData columns the examples subset on are present", {
+ cd <- SummarizedExperiment::colData(.se())
+
+ expect_true(all(c("condition", "passage") %in% colnames(cd)))
+ expect_s3_class(cd$condition, "factor")
+ expect_s3_class(cd$passage, "factor")
+ expect_identical(levels(cd$condition), c("ctrl", "aa57"))
+ expect_identical(levels(cd$passage), c("early", "late"))
+
+ # `se[, se$condition == "ctrl"]` must leave both passage levels populated,
+ # otherwise `~ passage` is not estimable.
+ ctrl <- cd[cd$condition == "ctrl", , drop = FALSE]
+ expect_gt(nrow(ctrl), 0L)
+ expect_setequal(as.character(unique(ctrl$passage)), c("early", "late"))
+})
+
+test_that("the packaged regions carry the feature_class column", {
+ reg <- .reg()
+ expect_s4_class(reg, "GRanges")
+ expect_gt(length(reg), 0L)
+ expect_true("feature_class" %in% colnames(S4Vectors::mcols(reg)))
+ expect_false(is.null(names(reg)))
+ expect_identical(anyDuplicated(names(reg)), 0L)
+})
+
+test_that("regions and probes overlap at the example threshold", {
+ # If this fails, every example calling fit_bread() on the packaged data
+ # errors with "no regions retained".
+ mapping <- map_probes_to_features(.se(), .reg(), min_probes = 3L)
+ expect_gt(nrow(mapping), 0L)
+ expect_gt(length(unique(mapping$region_id)), 0L)
+ expect_true("feature_class" %in% colnames(mapping))
+})
+
+test_that("the documented example fit runs end to end", {
+ se <- .se()
+ reg <- .reg()
+ se_ctrl <- se[, se$condition == "ctrl"]
+
+ fit <- fit_bread(se_ctrl, reg, ~ passage,
+ feature_class_col = "feature_class")
+ expect_s4_class(fit, "BreadFit")
+
+ res <- results(fit)
+ expect_gt(nrow(res), 0L)
+ expect_true(all(c("region_id", "classification") %in% colnames(res)))
+ expect_true(any(!is.na(res$prob_hyper)))
+})
diff --git a/tests/testthat/test-fit-bread.R b/tests/testthat/test-fit-bread.R
index e8015bf..c183758 100644
--- a/tests/testthat/test-fit-bread.R
+++ b/tests/testthat/test-fit-bread.R
@@ -48,7 +48,7 @@ test_that("results() and classifications() accessors return expected objects", {
min_probes = 3L))
r <- results(fit)
expect_s3_class(r, "data.frame")
- expect_true(all(c("region_id","classification","p_gt_delta","p_lt_neg_delta")
+ expect_true(all(c("region_id","classification","prob_hyper","prob_hypo")
%in% colnames(r)))
cls <- classifications(fit)
expect_type(cls, "character")
diff --git a/tests/testthat/test-fit-brms.R b/tests/testthat/test-fit-brms.R
index 53a2e7c..3b50804 100644
--- a/tests/testthat/test-fit-brms.R
+++ b/tests/testthat/test-fit-brms.R
@@ -1,7 +1,14 @@
-# Slow: Stan compile + sampling. ~60s on HPC. Skipped on CRAN, on systems
-# without brms, and when $_R_CHECK_FORCE_SUGGESTS_ is FALSE and brms missing.
+# Slow: Stan compile + sampling. ~60s on HPC. Skipped on CRAN, on CI, and
+# on systems without a working brms/rstan Stan toolchain.
test_that("fit_bread(backend = 'brms') recovers injected signal", {
skip_on_cran()
+ # skip_on_ci() is load-bearing and not redundant with the guards below:
+ # GitHub runners install brms and rstan happily but lack the headers Stan
+ # needs to compile a model (RcppEigen), so the test cleared every
+ # skip_if_not_installed() and then died with "Eigen not found".
+ # Installing a full Stan toolchain per CI run costs ~10 min and is flaky.
+ # The test still runs on the HPC and anywhere CI is unset.
+ skip_on_ci()
skip_if_not_installed("brms")
skip_if_not_installed("rstan")
@@ -41,7 +48,7 @@ test_that("fit_bread(backend = 'brms') recovers injected signal", {
res <- results(fit)
expect_true(is.na(res$df[1L])) # df undefined for empirical path
expect_true(all(!is.na(res$mean_effect)))
- expect_true(all(res$p_pos + res$p_neg >= 0.999))
+ expect_true(all(res$prob_pos + res$prob_neg >= 0.999))
# posterior_draws returns actual MCMC draws (subsampled to n)
d <- posterior_draws(fit, region_id = "regA", n = 200L, seed = 1L)
diff --git a/tests/testthat/test-fit-summary.R b/tests/testthat/test-fit-summary.R
index 2f5ef22..3a22f9e 100644
--- a/tests/testthat/test-fit-summary.R
+++ b/tests/testthat/test-fit-summary.R
@@ -44,21 +44,27 @@ test_that("posterior_summary recovers sign and orders true effects", {
expect_true(all(post$mean_effect[5:6] < -0.25))
expect_true(all(abs(post$mean_effect[1:2]) < 0.25))
# Probs in [0,1]
- for (col in c("p_pos","p_neg","p_gt_delta","p_lt_neg_delta"))
+ for (col in c("prob_pos","prob_neg","prob_hyper","prob_hypo"))
expect_true(all(post[[col]] >= 0 & post[[col]] <= 1))
- # p_pos + p_neg == 1 (within tolerance)
- expect_equal(post$p_pos + post$p_neg, rep(1, nrow(post)),
+ # prob_pos + prob_neg == 1 (within tolerance)
+ expect_equal(post$prob_pos + post$prob_neg, rep(1, nrow(post)),
tolerance = 1e-8)
})
test_that("classify_regions recovers truth with strong signal", {
+ # The ROPE call needs the whole posterior inside +/- delta, which is a much
+ # tighter demand than a directional call. At the original sigma = 0.2 the
+ # contrast SE was ~0.052 and a true null landed on prob_rope ~0.948, two
+ # thousandths under the cutoff -- a coin flip. sigma = 0.08 puts the SE at
+ # ~0.021, so even a null that happens to sit 2 SE off zero still carries
+ # >0.99 of its mass inside the ROPE.
sim <- .sim_region_mat(
n_samples = 60L,
true_betas = c(rep(0, 3), # nulls
rep(0.8, 3), # hyper
rep(-0.8, 3), # hypo
rep(0.02, 3)), # weak (below delta)
- seed = 42L, sigma = 0.2
+ seed = 42L, sigma = 0.08
)
fit <- fit_bread_summary(sim$mat, sim$coldata,
design = ~ group, contrast = "groupold")
@@ -66,10 +72,12 @@ test_that("classify_regions recovers truth with strong signal", {
cls <- classify_regions(post, delta = 0.10, prob_cutoff = 0.95)
got <- as.character(cls$classification)
- expect_true(all(got[1:3] == "inconclusive"), info = paste(got[1:3], collapse=","))
+ # The nulls and the sub-delta regions were always *scientifically*
+ # unchanged; calling them `inconclusive` was the defect this class fixes.
+ expect_true(all(got[1:3] == "unchanged"), info = paste(got[1:3], collapse=","))
expect_true(all(got[4:6] == "hypermethylated"), info = paste(got[4:6], collapse=","))
expect_true(all(got[7:9] == "hypomethylated"), info = paste(got[7:9], collapse=","))
- expect_true(all(got[10:12] == "inconclusive"), info = paste(got[10:12],collapse=","))
+ expect_true(all(got[10:12] == "unchanged"), info = paste(got[10:12],collapse=","))
expect_identical(attr(cls, "delta"), 0.10)
expect_identical(attr(cls, "prob_cutoff"), 0.95)
})
diff --git a/tests/testthat/test-kycg.R b/tests/testthat/test-kycg.R
new file mode 100644
index 0000000..d74438d
--- /dev/null
+++ b/tests/testthat/test-kycg.R
@@ -0,0 +1,76 @@
+# bread_kycg() reaches out to KnowYourCG reference databases, so these tests
+# deliberately cover only the validation that happens before any network or
+# annotation-hub access.
+
+.toy_fit <- function() {
+ fit_bread(.make_toy_signal_se(), .make_toy_features(), ~ group)
+}
+
+test_that("a non-BreadFit is rejected before any database access", {
+ expect_error(bread_kycg(data.frame(x = 1)), "BreadFit")
+ expect_error(bread_kycg(NULL), "BreadFit")
+})
+
+test_that("`which` and `platform` are matched against their allowed values", {
+ skip_if_not_installed("knowYourCG")
+ fit <- .toy_fit()
+ expect_error(bread_kycg(fit, platform = "NotAPlatform"))
+ expect_error(bread_kycg(fit, which = "sideways"))
+})
+
+test_that("the fit exposes the mapping columns bread_kycg() reads", {
+ # Not a network test: pins the columns bread_kycg() depends on, so a
+ # refactor of map_probes_to_features() cannot silently break it.
+ fit <- .toy_fit()
+ expect_true(all(c("probe_id", "region_id") %in% colnames(fit@mapping)))
+ expect_type(fit@mapping$probe_id, "character")
+})
+
+# The real group titles as registered by knowYourCG (verified against the
+# installed release). Used as a fixture so default-database selection is
+# testable without touching ExperimentHub.
+.KYCG_TITLES_MM285 <- c(
+ "KYCG.MM285.chromHMM.20210210",
+ "KYCG.MM285.chromosome.mm10.20210630",
+ "KYCG.MM285.designGroup.20210210",
+ "KYCG.MM285.HMconsensus.20220116",
+ "KYCG.MM285.Mask.20220123",
+ "KYCG.MM285.metagene.20220126",
+ "KYCG.MM285.probeType.20210630",
+ "KYCG.MM285.seqContext.20210630",
+ "KYCG.MM285.seqContextN.20210630",
+ "KYCG.MM285.TFBSconsensus.20220116",
+ "KYCG.MM285.tissueSignature.20211211"
+)
+
+test_that("mouse default databases actually match the real MM285 titles", {
+ # The regression: the old pattern required a literal "." after "TFBS", so
+ # it could never match "KYCG.MM285.TFBSconsensus.20220116" and mouse users
+ # silently received an empty data.frame.
+ dbs <- BREAD:::.kycg_default_dbs("MM285", .KYCG_TITLES_MM285)
+
+ expect_true("KYCG.MM285.TFBSconsensus.20220116" %in% dbs)
+ expect_true("KYCG.MM285.chromHMM.20210210" %in% dbs)
+ expect_true("KYCG.MM285.HMconsensus.20220116" %in% dbs)
+ expect_length(dbs, 6L)
+})
+
+test_that("technical annotation groups are deliberately excluded", {
+ dbs <- BREAD:::.kycg_default_dbs("MM285", .KYCG_TITLES_MM285)
+ for (junk in c("Mask", "chromosome", "probeType", "seqContext")) {
+ expect_false(any(grepl(junk, dbs, fixed = TRUE)), info = junk)
+ }
+})
+
+test_that("default selection returns nothing for an unlisted platform", {
+ expect_length(BREAD:::.kycg_default_dbs("EPIC", .KYCG_TITLES_MM285), 0L)
+ expect_length(BREAD:::.kycg_default_dbs("MM285", character(0)), 0L)
+})
+
+test_that("human platforms select the documented families", {
+ titles <- c("KYCG.EPIC.TFBS.20210210", "KYCG.EPIC.chromHMM.20211020",
+ "KYCG.EPIC.CGI.20210713", "KYCG.EPIC.Mask.20220123")
+ dbs <- BREAD:::.kycg_default_dbs("EPIC", titles)
+ expect_length(dbs, 3L)
+ expect_false(any(grepl("Mask", dbs, fixed = TRUE)))
+})
diff --git a/tests/testthat/test-mapping.R b/tests/testthat/test-mapping.R
index 34fc847..26f1a69 100644
--- a/tests/testthat/test-mapping.R
+++ b/tests/testthat/test-mapping.R
@@ -52,3 +52,51 @@ test_that("mapping rejects non-positive min_probes", {
expect_error(map_probes_to_features(se, gr, min_probes = 0L),
"positive integer")
})
+
+test_that("several ranges sharing a region_id collapse into one region", {
+ se <- .make_toy_se(); gr <- .make_toy_features_dup()
+ expect_length(gr, 3L) # three ranges ...
+ expect_length(unique(names(gr)), 2L) # ... but two regions
+
+ m <- suppressMessages(map_probes_to_features(se, gr, min_probes = 3L))
+
+ expect_setequal(unique(m$region_id), c("regD", "regE"))
+ # regD spans probes 1-3 and 8-10 across two disjoint ranges
+ expect_equal(sum(m$region_id == "regD"), 6L)
+ expect_true(all(m$n_probes[m$region_id == "regD"] == 6L))
+
+ # Counts are of distinct region IDs, never of ranges. This is the
+ # regression: n_features_in used to report 3 here.
+ expect_equal(attr(m, "n_features_in"), 2L)
+ expect_equal(attr(m, "n_features_out"), 2L)
+})
+
+test_that("dropped_regions is deduplicated and the message counts regions", {
+ se <- .make_toy_se()
+ gr <- .make_toy_features_dup()
+ # Raise the bar so regD (6 probes across 2 ranges) survives but regE (5) does not
+ expect_message(
+ m <- map_probes_to_features(se, gr, min_probes = 6L),
+ "Dropped 1 of 2 regions"
+ )
+ expect_equal(attr(m, "dropped_regions"), "regE")
+
+ # And when the multi-range region itself is dropped, it appears once
+ m2 <- suppressMessages(map_probes_to_features(se, gr, min_probes = 20L))
+ expect_equal(sort(attr(m2, "dropped_regions")), c("regD", "regE"))
+})
+
+test_that("fit diagnostics count regions, not ranges", {
+ se <- .make_toy_se(); gr <- .make_toy_features_dup()
+ fit <- suppressMessages(fit_bread(se, gr, ~ group, min_probes = 3L))
+
+ # The invariant that would have caught "n_regions: 788 (of 790 input)"
+ expect_equal(fit@diagnostics$n_features_out, nrow(results(fit)))
+ expect_equal(fit@diagnostics$n_features_in, 2L)
+
+ # Documented property: the features slot keeps every range of a surviving
+ # region, so it is longer than the results table when IDs repeat.
+ expect_gt(length(fit@features), nrow(results(fit)))
+
+ expect_output(show(fit), "n_regions : 2 \\(of 2 input\\)")
+})
diff --git a/tests/testthat/test-methods.R b/tests/testthat/test-methods.R
new file mode 100644
index 0000000..29a38b5
--- /dev/null
+++ b/tests/testthat/test-methods.R
@@ -0,0 +1,79 @@
+# Accessors and show methods for BreadFit / BreadResults.
+
+.toy_fit <- function() {
+ fit_bread(.make_toy_signal_se(), .make_toy_features(), ~ group)
+}
+
+test_that("results() returns the region table with a classification column", {
+ res <- results(.toy_fit())
+ expect_s3_class(res, "data.frame")
+ expect_true(all(c("region_id", "classification") %in% colnames(res)))
+ expect_gt(nrow(res), 0L)
+})
+
+test_that("classifications() names line up with results() region ids", {
+ fit <- .toy_fit()
+ cls <- classifications(fit)
+ res <- results(fit)
+ expect_type(cls, "character")
+ expect_identical(names(cls), as.character(res$region_id))
+ expect_identical(unname(cls), as.character(res$classification))
+ expect_true(all(cls %in% c("hypermethylated", "hypomethylated",
+ "unchanged", "inconclusive")))
+})
+
+test_that("posterior_draws() is reproducible for a fixed seed", {
+ fit <- .toy_fit()
+ rid <- results(fit)$region_id[1]
+ a <- posterior_draws(fit, region_id = rid, n = 200L, seed = 7L)
+ b <- posterior_draws(fit, region_id = rid, n = 200L, seed = 7L)
+ expect_equal(a, b)
+ expect_identical(colnames(a), c("region_id", "draw", "value"))
+ expect_identical(nrow(a), 200L)
+ expect_identical(unique(a$region_id), rid)
+})
+
+test_that("posterior_draws() differs across seeds", {
+ fit <- .toy_fit()
+ rid <- results(fit)$region_id[1]
+ a <- posterior_draws(fit, region_id = rid, n = 200L, seed = 1L)
+ b <- posterior_draws(fit, region_id = rid, n = 200L, seed = 2L)
+ expect_false(isTRUE(all.equal(a$value, b$value)))
+})
+
+test_that("posterior_draws() restores the caller's RNG state", {
+ # local_seed() must not leak its reseed into the calling session.
+ fit <- .toy_fit()
+ rid <- results(fit)$region_id[1]
+ set.seed(99)
+ before <- runif(1)
+ set.seed(99)
+ invisible(posterior_draws(fit, region_id = rid, n = 10L, seed = 123L))
+ after <- runif(1)
+ expect_equal(before, after)
+})
+
+test_that("posterior_draws() defaults to every region", {
+ fit <- .toy_fit()
+ n_regions <- nrow(results(fit))
+ d <- posterior_draws(fit, n = 10L, seed = 1L)
+ expect_identical(nrow(d), as.integer(n_regions * 10L))
+})
+
+test_that("an unknown region_id errors and names the offender", {
+ expect_error(posterior_draws(.toy_fit(), region_id = "no_such_region"),
+ "not found")
+})
+
+test_that("show() prints the expected BreadFit header", {
+ fit <- .toy_fit()
+ expect_output(show(fit), "