Prepare BREAD 0.99.0 for Bioconductor submission - #1
Open
BacZemin wants to merge 17 commits into
Open
Conversation
Bioconductor requires new submissions to be versioned x.99.z, and until now nothing verified the package on a clean machine -- the only workflow was pkgdown. - R-CMD-check.yaml: ubuntu/release only for now; more platforms once green. - bioc-check.yaml: runs R CMD check + BiocCheck(new-package = TRUE) in the bioconductor/bioconductor_docker:devel container. The devel R/Bioc versions are printed by the job rather than hardcoded anywhere. - DESCRIPTION: 0.99.0, Depends: R (>= 4.4.0), drop the unused GenomeInfoDb import (R CMD check NOTE). - .Rbuildignore: data-raw, build_preview and CLAUDE.md were all being packaged into the tarball. - .gitignore: keep CLAUDE.md and tools/teaching out of a public repo.
BiocCheck ERROR: at least 80% of man pages documenting exported objects must have runnable examples. Only bread_kycg had one, and it was a dontrun block over undefined objects. - Every exported object now has an example built on the packaged EPICv2 data. A full 500-region fit is ~0.7s, so each example is self-contained rather than sharing a cached fit. Examples that call the lower-level helpers pass assay_name = "betas" / input_scale = "Beta" explicitly, because only fit_bread() auto-detects them. Plot and draw examples always name a region: the NULL default would facet 500 panels / return 2M rows. - bread_kycg keeps a dontrun block for the enrichment call itself, which downloads KnowYourCG reference databases, but the surrounding fit now runs. - posterior_summary() accepts a BreadFit directly, so callers no longer reach into the model slot. - Add a BreadResults() constructor and show method. The class had no way to construct one anywhere in the package. - Remove report_feature_set(), an unimplemented stop() stub; plot_feature_set() covers the same ground. - posterior_draws() uses withr::with_seed() instead of set.seed(). BiocCheck flags set.seed in package code, and the bare call also leaked its reseed into the caller RNG stream. (with_seed, not local_seed: the latter does not restore the previous state in withr 3.0.2.) - Add missing \value sections flagged by BiocCheck.
…xtdata Six R files had no test file at all, including posterior.R and classify.R -- the statistical core. - test-classify.R drives the decision rule with hand-built posterior tables so the cutoff boundary is exact: >= is inclusive, a stricter cutoff can only move regions to inconclusive, NA rows never yield NA classifications. - test-posterior.R pins the 13-column contract, checks that a wider ci widens intervals without moving point estimates, that p_pos + p_neg == 1, and that the new BreadFit fast-path is identical to passing the model list. - test-methods.R covers accessor/name alignment, seed reproducibility, and that posterior_draws() restores the caller RNG state. - test-classes.R pins the S4 slot layout, which the plotting helpers and accessors depend on by name. - test-utils.R covers the M/beta round trip, boundary clamping, and the documented assay-name priority order. - test-kycg.R covers argument validation only -- no network calls. - test-extdata-contract.R pins the packaged data shape. Every new example is written against those two objects, so a rename or a dropped colData column would otherwise surface as an opaque R CMD check failure. 284 tests pass.
The README hard-coded its example output as a hand-typed comment, and an earlier draft had already shipped the wrong classification labels. It is now knitted from README.Rmd, so the 76/39/385 counts are generated rather than asserted. - inst/CITATION, picking the version up from DESCRIPTION. No DOI yet; add one via bibentry(doi = ) once there is a preprint. - .github/CONTRIBUTING.md, SUPPORT.md, and CODE_OF_CONDUCT.md (adopting the Bioconductor Code of Conduct), plus issue templates that route usage questions to the Bioconductor support site. - NEWS.md: heading is now a bare version, which the Bioconductor markdown NEWS parser requires. - _pkgdown.yml: the reference section claimed the backend fit helpers were exposed, but they are internal and unexported; the news block still advertised 0.0.0.9000.
Dropping GenomeInfoDb from Imports broke vignette building on both CI
workflows: on Bioconductor devel the Seqinfo class has been split into its own
package but is still registered as belonging to GenomeInfoDb, so
findOverlaps() -> merge(seqinfo) -> methods::is() calls
.requirePackage("GenomeInfoDb") and fails if it was never installed.
R CMD check's "Namespace in Imports field not imported from" NOTE was about
the namespace, not the dependency -- the dependency is real and indirect. It
is now a genuine import: map_probes_to_features() already told users to check
seqlevels() when nothing overlapped, so the error now reports them instead of
just naming the function.
The HPC preflight did not catch this because GenomeInfoDb is installed there
as a transitive dependency of the current GenomicRanges.
BiocCheck on Bioc devel reported "Unable to read CITATION file". The file
fell back to packageVersion("BREAD") when meta was NULL, but BiocCheck
validates the CITATION against an uninstalled package, so that call errored.
The fallback is now a literal. It passed on the HPC only because a stale
BREAD 0.0.0.9000 happened to be installed there.
Also:
- DESCRIPTION Description expanded to three sentences (BiocCheck NOTE),
covering the two backends and the knowYourCG hand-off.
- bioc-check now fails on any BiocCheck ERROR or WARNING beyond the two that
cannot be resolved in code (support-site registration; the case-insensitive
name collision with CRAN's unrelated 'bread'). Previously the workflow
reported findings but always went green.
BiocCheck on devel warns when a CITATION file has no doi argument, and
Bioconductor's own guidance is to include one only once a preprint or
publication exists. There is none yet, so the file goes; citation("BREAD")
still works, generated from DESCRIPTION.
Add it back with bibentry(doi = ...) once there is a preprint or a Zenodo DOI.
DESCRIPTION gains a Date field so the generated citation carries a year
instead of (????); the README renders citation() live and was showing the
placeholder.
R-CMD-check failed on test-fit-brms.R with "Eigen not found". The existing
skip_if_not_installed("brms") / ("rstan") guards all pass on GitHub runners
-- both packages install fine -- but Stan then cannot compile a model because
the RcppEigen headers are absent, so the test ran and errored.
Installing a working Stan toolchain on every CI run costs about ten minutes
and is flaky, so the test is skipped on CI instead. It still runs on the HPC,
and on build systems that do not set CI, which is where the brms backend is
actually exercised.
p_gt_delta and p_lt_neg_delta read as p-values, which is exactly the inference they are not. Renamed to prob_hyper / prob_hypo, matching the classification factor levels, and p_pos / p_neg to prob_pos / prob_neg. These are posterior probabilities of the parameter given the data -- prob_hyper = P(effect > +delta), prob_hypo = P(effect < -delta) -- not tail probabilities of a statistic under a null. Clean break, no deprecation shim: 0.99.0 is unreleased.
Several GRanges may share one name -- the only way to pin a region to an exact probe set, since a bounding interval sweeps in neighbours. split() handled this correctly all along but the counters did not: n_features_in counted ranges (length(features)) show() printed length(object@features) as n_regions which produced "n_regions: 788 (of 790 input)" for what was really 30 regions built from 790 probes. Results were always right; only the counters lied. Adds .make_toy_features_dup() (3 ranges, 2 regions) and pins the invariant n_features_out == nrow(results(fit)).
The headline is a fourth classification level. `inconclusive` was absorbing two different situations: a region whose posterior sits tightly inside the ROPE (strong evidence of no change) and one whose posterior says nothing at all. Collapsing them throws away the single claim a p-value structurally cannot make. The mousearray_609G CArG result is the case in point -- REGULON_promoter carries P(|effect| < delta) = 0.958 and was labelled `inconclusive`, indistinguishable from noise. classification levels are now hypermethylated / hypomethylated / unchanged / inconclusive. `unchanged` is inserted third so levels()[1:2] stay the directional pair. Membership of the two directional classes is bit-identical for any prob_cutoff > 0.5 -- the new level strictly partitions the old `inconclusive` bucket -- so directional filters are unaffected. posterior_summary() gains six columns: prob_rope plus a beta-scale translation (ref_beta, mean_dbeta, dbeta_lo, dbeta_hi, delta_beta). The translation uses one linearisation multiplier per region, anchored at the region's own mean methylation, so the beta comparison can never contradict the M-scale classification printed beside it. delta stays M-scale only: a beta-defined threshold is necessarily per-region and would silently become a 20-fold wider ROPE at the methylation extremes. bread_delta_beta() and bread_delta_m() convert explicitly instead. Also: * rope_cutoff, separate from prob_cutoff but defaulting to it. Equivalence needs the whole posterior inside +/-delta, a far stricter demand than a directional call; at n = 6+6 a shared knob ships the feature unusable. * ci is now a fit_bread() argument. It was hardcoded at 0.95, distinct from prob_cutoff but unreachable. * Matrix input. openSesame() returns a matrix and sesameData examples are list(betas=, sampleInfo=); BREAD rejected both, so the sesame workflow died at step one. New bread_se() / .as_bread_se() accept them with colData plus either rowRanges or platform. The platform is never inferred from probe IDs -- cg-numbers are shared across HM450/EPIC/MM285 and a wrong guess yields wrong coordinates with no error. EPICv2 suffix mismatches get a specific diagnosis. * refit_bread(). Permutation calibration needs the region matrix computed once and only the fit repeated; without a public entry point the only route was BREAD:::fit_bread_summary(). Re-thresholding short-circuits the model fit entirely. * Rank-deficient designs now warn. fit_summary adds Lambda0 = 0.01*I before the Cholesky, so collinear designs never errored -- they silently returned prior-driven estimates. Every caller so far hand-rolled its own qr() guard. * bread_kycg(): the default pattern required a literal . after TFBS, so it could not match KYCG.MM285.TFBSconsensus.20220116 and mouse users got a silently empty data.frame. Now a per-platform table, a warning that lists what is available, and mtc_by_group/mtc_method passthrough filtered against the installed signature. Also fixes an rbind error reachable today whenever one classification level yields a stub and another a real table. * The toy fixture had sex perfectly collinear with group, so every test using ~ group + sex was silently fitting a rank-deficient design. 412 tests pass, up from 284.
Documents the four-level classification, the six new results() columns, the n_features_in region-vs-range fix, the fit_bread() first-argument rename, and the palette change, all as breaking changes -- 0.99.0 is unreleased so they land without a deprecation period.
README gains an explicit section on the unchanged class and on why delta is M-scale only, plus the matrix-input and refit_bread workflows. The vitc vignette n=2 paragraph is rewritten rather than patched: at n=2 almost nothing reaches prob_rope >= 0.95 either, so those regions are unresolved rather than shown to be flat -- which is a sharper statement than the old undifferentiated inconclusive bucket, and the best advertisement for the feature.
…nchanged At n=4 nothing reaches prob_rope >= 0.95, which is the honest result and the concrete argument for rope_cutoff being a separate knob. Better shown than hidden -- it also demonstrates the unresolved-vs-flat distinction the class exists to draw.
…d README - README described `n` as the probe count per region; it is the number of samples contributing to the fit after NA-dropping. Probe counts are in `mapping$n_probes`. Say so in posterior_summary()'s @return too, which had listed `n` without defining it. - Document how far the first-order beta translation drifts from an exact back-transform: median 0.0005 and p99 0.021 in beta units on the packaged vitamin C example, ~2% relative below |dM| 0.25 rising to ~10% above 0.5. No classification is affected, but mean_dbeta should not be quoted to three decimals for a strongly changing region. - README.md was stale: it still documented the three-level API and none of prob_rope / unchanged / the beta columns, because 5f60905 and 4c3c3c0 edited README.Rmd without rerunning build_readme(). Rebuilt. - ROPE interval was written \[-delta, +delta\]; pandoc read the escaped brackets as LaTeX and emitted $$..$$ display math into the GFM output. Reworded to avoid the brackets.
The conjugate update a_n = a0 + n/2 makes nu = 2*a_n a function of the
sample size alone and never of the number of coefficients p. That is the
literal NIG result and correct given the stated prior, but the default
prior (Lambda0 = 0.01 I) was only ever meant to be uninformative, and the
reference-prior answer -- the one lm() gives -- is n - p. The posterior
scale is too small by exactly sqrt(n / (n - p)), so the error grows with
p/n and bites hardest on interaction designs at small n.
Measured on the packaged vitamin C data with ~ condition + passage
(n = 8, p = 3): BREAD reported df = 8.002 and scale 0.04715 where lm gives
df = 5 and se 0.05559. Switching to df_mode = "residual" moves 24 of 93
hypermethylated and 5 of 50 hypomethylated regions to inconclusive -- 20%
of the directional calls were an artefact of the optimistic df.
Two changes:
* fit_bread(df_mode = c("conjugate", "residual")). "conjugate" stays the
default so existing results reproduce; "residual" uses
a_n = a0 + (n - p)/2 and reproduces the classical t_{n-p} marginal,
matching lm() exactly as Lambda0 -> 0 (tested at lambda0 = 1e-8).
Recorded in fit@params$df_mode and fit@model$df_mode. Ignored by the
brms backend, which samples sigma^2 directly.
* Regions with n <= p are dropped with
error = "n <= number of coefficients" rather than fitted. With no
residual information the residuals are identically zero, b_n collapses
to b0, and the scale collapses with it: on pure noise at n = p = 4 the
old path returned a scale ~7x tighter than at n = 8 and a credible
interval excluding zero. Only n < 2 was guarded before.
Also warns once per fit -- not once per region -- when any fitted region
has fewer than 3 residual degrees of freedom.
test-classes.R pins the fit@model field names as a deliberate contract;
updated for the added df_mode field.
man/figures/logo.png is picked up automatically by pkgdown for the navbar brand, the home-page header and the og:image card, so _pkgdown.yml needs no change. The README uses the anchor-wrapped form; pkgdown strips both the bare and anchored variants from the homepage body (tweak_homepage_html matches .//h1/img and .//h1/a/img) before inserting its own sidebar logo, so the logo is not duplicated there. The vignettes embed the logo as a base64 data URI via knitr::image_uri() rather than linking ../man/figures/logo.png. A plain relative link would work in the tarball -- html_vignette is self-contained, so pandoc inlines it -- but it breaks on the pkgdown site: copy_article_images() only copies resources whose path is a child of the vignette directory, and a parent- escaping path is dropped. Resolving the file at knit time sidesteps this; pkgdown's check_missing_images() skips any src with a non-empty URI scheme and handles data:image explicitly in its alt-text check. Verified: R CMD build ships BREAD/man/figures/logo.png and nothing from data-raw/ or pkgdown/; an isolated html_vignette render of the logo chunk emits one data URI, zero leftover relative paths, and the alt attribute. README.md was hand-mirrored rather than regenerated -- the edit is a single line above every code chunk and github_document passes raw HTML through verbatim, so build_readme() would only have re-run the live model fits and churned unrelated numbers into a logo commit.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Prepares BREAD for Bioconductor submission as 0.99.0. 8 commits, 65 files, +1839/−104.
Not merged, not tagged, not released, not submitted — see "Deliberately not done".
CI status
R-CMD-check(ubuntu / R release)bioc-check(Bioc devel container)Neither workflow existed before this branch.
R CMD checkon the Bioconductordevel container reports
Status: OK— 0 errors, 0 warnings, 0 notes.Devel toolchain, read from the container rather than assumed:
R 4.6.1 (2026-06-24) · Bioconductor 3.24 · BiocCheck 1.49.30 · x86_64-pc-linux-gnu
Test suite: 284 passing, 0 failures. Coverage: 74.73%
(
classify.Randclasses.Rat 100%;fit_brms.Rat 0% because the Stan testis skipped where there is no toolchain).
BiocCheck: baseline → now
Baseline on unmodified
mainwas 3 ERRORS | 4 WARNINGS | 11 NOTES.Now 1 ERROR | 1 WARNING | 10 NOTES, and both remaining items are outside
code control.
x.99.z(0.0.0.9000)0.99.0set.seed()in package code (R/methods.R:60)withr::with_seed()\valuesections (6 man pages)Depends: R (>= 4.4.0)ClassificationsuppressWarnings, function length, indentation,dontrun)R CMD checkseparately lost two NOTEs: non-standard top-level files and theunused
GenomeInfoDbimport. Only the environment-onlyqpdfwarning and"unable to verify current time" remain on the HPC, and neither reproduces on CI.
⛔ Decision needed: the package name collides with CRAN
BiocCheck's
checkIsPackageNameAlreadyInUse()comparestolower(pkgName) %in% tolower(CRAN packages)— case-insensitive by design.CRAN has
breadv0.4.1,"Analyze Big Files Without Loading Them in Memory", entirely unrelated.
This is not a BiocCheck quirk. On case-insensitive filesystems (macOS, Windows)
BREADandbreadcannot coexist in one library, so the rule is real.Bioconductor's stated remedy — "must be removed from CRAN before the next
release" — assumes the CRAN package is yours. It is not. The practical
outcome is that Bioconductor will require a rename.
A rename touches: the GitHub repo name, the pkgdown URL,
library()calls inboth vignettes, the thesis text, and the "bread & brie" pairing with the
sibling
BRIEpackage. Worth knowing: no CRAN or Bioconductor package isnamed
brie— onlyBREADcollides.I did not act on this. It is the single decision that should be made before
anything else here is merged, because a rename would touch most of this diff.
⛔ Other things only you can do
DESCRIPTION(jaemin.park@vai.org) — https://support.bioconductor.org/accounts/signup/.This is the last remaining BiocCheck ERROR.
BiocCheck cannot verify this without admin credentials (NOTE).
comment = c(ORCID = "...")toAuthors@R. I did not invent one — afabricated identifier in a public package is worse than a missing one.
What changed
CI and packaging
R-CMD-check.yaml— ubuntu/release only for now; add macOS and Windows onceit has been green for a while.
bioc-check.yaml—R CMD build→R CMD check→BiocCheck(new-package = TRUE)in
bioconductor/bioconductor_docker:devel. The devel R/Bioc versions areprinted by the job, never hardcoded. It fails on any BiocCheck ERROR or
WARNING beyond the two known-blocked ones, so it is a real regression gate
rather than a green badge.
pkgdown.yamluntouched — it triggers only onpush: main, so the live sitewas never rebuilt from this branch.
0.99.0,Depends: R (>= 4.4.0),Date, biocViewClassification,three-sentence
Description.Examples — the main event. All 17 exported objects now have runnable
examples on the packaged EPICv2 data. A full 500-region fit is ~0.7 s, so each
example is self-contained rather than sharing a cached fit. Three traps the
examples had to route around, all of which would have shipped as broken
documentation:
validate_bread_input()andsummarize_features()default toassay_name = "M", but the packaged assay is"betas"— onlyfit_bread()auto-detects, so those examples pass both arguments explicitly.
plot_region_posterior()andposterior_draws()default toregion_id = NULL,which means 500 facets and 2,000,000 rows respectively. Every example names a
region.
defensively with a fallback.
API
posterior_summary()accepts aBreadFitdirectly instead of requiringfit@model.BreadResults()constructor andshowmethod. The class waspreviously unconstructable — nothing in the package ever called
new("BreadResults", ...), and it had noshowmethod. See the openquestion below.
report_feature_set(), astop("not yet implemented")stub that wasshipping with a live man page.
plot_feature_set()covers the same ground.posterior_draws()useswithr::with_seed()instead ofset.seed(). Thisfixes a real defect as well as the BiocCheck warning: the bare
set.seed()leaked its reseed into the caller's RNG stream. Note
withr::local_seed()does not restore state in withr 3.0.2 — verified, hence
with_seed().Tests — 7 new files covering the six R files that had none, including
posterior.Randclassify.R, the statistical core.test-extdata-contract.Rpins the packaged data shape so a rename or droppedcolDatacolumn surfaces as a named failure instead of an opaqueR CMD checkerror in 12 examples at once.Docs —
README.mdis now generated fromREADME.Rmd. It previouslyhard-coded
76 / 39 / 385as a hand-typed comment, and an earlier draft hadalready shipped the wrong classification labels; those counts are now
generated. Added
.github/CONTRIBUTING.md,SUPPORT.md,CODE_OF_CONDUCT.md(adopting the Bioconductor CoC), and issue templates thatroute usage questions to the Bioconductor support site. Fixed
_pkgdown.yml,which claimed the backend fit helpers were exposed (they are internal) and
still advertised
0.0.0.9000.Two things worth a second look
CLAUDE.mdwas being packaged. The baselineR CMD checklisted it undernon-standard top-level files — it was neither
.gitignored nor.Rbuildignored, so it was going into every built tarball. Both are fixed, andI audited the branch history to confirm nothing leaked into a commit.
GenomeInfoDbis a real dependency.R CMD checkreported"Namespace in Imports field not imported from", so I dropped it — and both CI
workflows immediately failed to build vignettes. On Bioc devel the
Seqinfoclass has moved to its own package but is still registered as belonging to
GenomeInfoDb, sofindOverlaps()→merge(seqinfo)→methods::is()calls.requirePackage("GenomeInfoDb"). The NOTE was about the namespace, not thedependency. It is restored, and now genuinely imported:
map_probes_to_features()already told users to checkseqlevels()whennothing overlapped, so it reports them. The HPC preflight missed this because
GenomeInfoDbis installed there as a transitive dependency.Deferred, with reasons
inst/CITATIONremoved. BiocCheck warns when a CITATION has nodoi,and Bioconductor's guidance is to ship one only once a preprint or
publication exists.
citation("BREAD")still works, generated fromDESCRIPTION. Add it back with
bibentry(doi = ...)when there is a DOI.Depends: R (>= 4.4.0)not bumped to 4.6.0. BiocCheck raises this onlyas a NOTE. Bumping it would make the package uninstallable on the HPC
(R 4.4.3) and break the preflight loop. Worth revisiting at submission.
fit_breadis 125). NOTE only; refactoringthe main entry point unattended was not a good trade.
stylerrun would burythis diff.
suppressWarnings/suppressMessages. Each is deliberate; the one inmap_probes_to_features()is now explained in a comment.\dontrun{}onbread_kycg(). BiocCheck suggestsdonttestinstead,but
donttestblocks are run byR CMD check, and that call downloadsKnowYourCG reference databases. The surrounding fit does run.
RcppEigen headers Stan needs, so the test cleared every
skip_if_not_installed()and then died with "Eigen not found". A full Stantoolchain per run costs ~10 minutes and is flaky. It still runs on the HPC.
Open question
BreadResultsis a thin wrapper with no consumer anywhere in the package —its only plausible caller was the
report_feature_set()stub this PR deletes.I gave it a constructor and a
showmethod, which is additive and reversible.The alternative is deleting the class before 0.99.0 freezes the API, which
is arguably the better long-term call. It is announced in
NEWS.mdand_pkgdown.yml, so it is your decision, not one to make unattended.Deliberately not done
No tag. No GitHub release. No Zenodo. No issue on
Bioconductor/Contributions.Not merged.
origin/mainis untouched atcc229eb.Repo settings changed outside this branch (the only such changes):
homepage→https://baczemin.github.io/BREAD/(it pointed at the repoitself), 7 topics added, and the empty wiki and projects tabs disabled.