Guarded and auditable discovery of prognostic RNA-seq signatures.
rnaSentry orchestrates the statistical pipeline that goes from a raw bulk
RNA-seq count matrix to a validated prognostic gene signature. Every
automated decision is reported with the statistic that justified it, and the
pipeline gates rather than silently degrades when its guardrails are not
met. The package is designed to make a biomarker discovery workflow
transparent, reproducible, and reviewable end to end.
rnaSentry is a development release (0.99.x) under active development. The case-study numbers in this README and the package vignettes are illustrative walkthroughs of the pipeline on public datasets. They are not peer-reviewed and should be treated as demonstrations, not clinical claims.
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("rnaSentry")
library(rnaSentry)
library(SummarizedExperiment)
# counts: a gene x sample integer matrix, gene symbols as rownames.
# coldata: a data.frame of sample metadata (one row per sample).
# load_counts() validates the intake and records an audit flag ledger.
se <- load_counts(counts, coldata)
# Run the discovery pipeline (time/event are survival columns in colData).
run <- run_rnaSentry(se, time_col = "time", event_col = "event",
design_vars = c("batch", "age"))
# Optional: parallelize the repeated cross-validation folds.
# run <- run_rnaSentry(se, time_col = "time", event_col = "event",
# design_vars = c("batch", "age"),
# BPPARAM = BiocParallel::SnowParam(2))
# Validate the signature against an independent cohort.
val <- validate_external(run$stages$build_signature, external_se,
cutpoint = run$stages$km_curve$cutpoint)
See vignette("getting-started") to bridge your own CSV files into the
pipeline, vignette("rnaSentry") for a worked walkthrough, and
?run_rnaSentry, ?build_signature, ?validate_external for details.
Most users start with two plain files: a counts matrix (genes in rows,
samples in columns) and a clinical table (one row per sample). rnaSentry
accepts a SummarizedExperiment, so the only bridge you need is the small
block below. It reads both files, keeps exactly the samples present in both
(mismatched sample IDs are the #1 pitfall), and validates the result with
load_counts(), which records any intake issues in a flag ledger on the
object.
library(rnaSentry)
library(SummarizedExperiment)
# 1. Read the two CSV files.
counts <- as.matrix(read.csv("counts.csv", row.names = 1))
clinical <- read.csv("clinical.csv", row.names = 1)
# 2. Keep samples present in BOTH files (the #1 pitfall: mismatched IDs).
common <- intersect(colnames(counts), rownames(clinical))
counts <- counts[, common, drop = FALSE]
clinical <- clinical[common, , drop = FALSE]
# check: dim(counts)[2] == nrow(clinical)
# 3. Build + validate the SummarizedExperiment (issues -> flag ledger).
se <- load_counts(counts, clinical)
S4Vectors::metadata(se)$flags
# 4. Audit candidate confounders before any modeling.
da <- design_audit(se, design_vars = c("batch", "stage"))
da$formula_text
# 5. Run the full pipeline (needs 'time' + 'event' columns in clinical.csv).
# run <- run_rnaSentry(se, time_col = "time", event_col = "event",
# design_vars = c("batch", "stage"))
Notes on your input files:
load_counts() expects a
matrix, so wrap the read.csv() result in as.matrix(). Column names must
be sample IDs; row names should be gene symbols (if your file uses Ensembl
accessions, map them to symbols first -- the pipeline warns otherwise).intersect() keeps only
samples found in both; samples missing from either file are silently
dropped. Check the dimensions after intersecting -- if you lose hundreds of
samples you may have a versioning or delimiter mismatch.load_counts() stores your
matrix under the assay name counts, and the pipeline computes
log2(counts + 1) when it runs. If your file is already normalized
(e.g. CPM, TPM, or VST), name that assay logcounts instead -- either pass
load_counts(..., assay_name = "logcounts") or rename after construction.time column and a 0/1 event indicator in the clinical table.
design_audit() and the other auditing stages do not need them.run_rnaSentry() is single-use per session: after it runs, the signature is
locked and a re-run is refused until the lock is released with
lock_signature(sig, lock = FALSE). This prevents silent re-selection of
signature genes after survival analysis, which would otherwise inflate
reported concordance. Individual stages (for example validate_external())
remain re-runnable on the locked signature.
The concordance reported by build_signature() is a screening-internal
metric: genes are selected on the full cohort before the CV split, so the
CV fold is not fully untouched by selection. On a large candidate panel the
value is optimistic: on GSE20685, identical-protocol builds on
survival-permuted data give CV C = 0.744–0.785 across 6 draws (vs 0.785 for
the real signature — inside, not above, the noise band), versus 0.582 for a
no-selection random-gene control (face-validity run record) and 0.490 for the
same no-selection protocol on simulated null data. Always interpret it
relative to a matched null (a permutation-null or random-gene control),
and treat validate_external() on an independent cohort as the only fully
out-of-sample estimate.
logcounts/vst (RNA-seq) or exprs (microarray log2
intensities).load_counts() to validate the count matrix and metadata before the
pipeline runs; it records any intake issues (duplicate samples,
non-integer counts, all-zero gene rows, non-syntactic gene symbols) in the
object's flag ledger.counts. The pipeline
detects this and warns (possibly_log_scaled), but you should store such
data under the name logcounts (RNA-seq) or exprs (microarray) instead.build_signature() warns
(events_per_parameter). Consider fewer genes or a larger cohort, or a
regularized approach.build_signature() warns when events per signature gene fall
below min_events_per_parameter (default 5, mid-range of the 5-9
events-per-variable guidance; see ?build_signature for the sources).
Treat concordance and hazard ratios from such cohorts as indicative, not
confirmatory.BiocParallelParam (e.g. BiocParallel::SnowParam(2))
to run_rnaSentry()/build_signature() via BPPARAM. Parallel fold
evaluation is bit-identical to serial (same seeded folds), and memory is
usually the binding constraint on large cohorts regardless.Every case study is written as a three-tier narrative: a naive analysis
(what a standard, unguarded pipeline would report), rnaSentry standing
guard (the same data through the pipeline, with the flag ledger), and the
counterfactual (the error prevented and the true signal recovered).
vignette("case-study-impact") is the landing page with the master impact
table. Three studies run on synthetic cohorts so the pipeline's behaviour can
be inspected offline, with a synthetic fallback behind every live-data step;
the BRCA study additionally pulls live GEO data, and the fourth (LUAD) is
reproduced from live GEO and documented here.
| Failure mode | Naive headline | Guarded headline |
|---|---|---|
| Wrong direction (GSE20685) | inverted-convention read (C ~ 0.2) looks like a failure | held-out C = 0.783 vs 0.582 random-gene control (CV alone does not exceed selection-noise controls; see case-study Limitations) |
| Confounded design (synthetic cohort) | six batch-driven genes are "prognostic"; batch + region are "independent" | batch/region redundant (Cramer's V = 0.82); design reduced to ~ batch
|
| Underpowered discovery (30 samples, 23 events) | "mean CV C = 0.74 - decent" | fold range 0.50-1.00, sd 0.14, 4.6 events/parameter |
| False transfer (GSE31210 -> GSE50081) | "CV C = 0.86, validated" | external C = 0.540, p = 0.266: transfer not demonstrated |
A 3000-gene x 327-sample subset of GSE20685 (Li et al., 2010; Affymetrix
GPL570, 83 deaths) is fetched via GEOquery::getGEO("GSE20685") — see
inst/scripts/repro_gse20685.R to regenerate a local
SummarizedExperiment (not committed to inst/extdata; vignettes use a
synthetic offline fallback with the same structure) — and is walked through
in vignette("case-study-brca"). On the live-GEO cohort (3000 genes
x 327 samples, 83 deaths; GEO annotation of September 2026) the
pipeline finds a 20-gene signature with screening-internal CV C = 0.78
(sd 0.03) and log-rank p = 3.3e-17 at the discovery cutpoint; the design
audit flags subtype (eta-squared = 0.79) and recommends adjusting for it.
The narrative contrasts this with the naive default-convention concordance
read (C ~ 0.2, i.e. "below chance"), which the reverse = TRUE convention
and the C + C_rev = 1 parity test prevent -- the exact mechanism behind the
historical concordance-inversion bug this convention was designed against.
The signature is provisional:
at 83 events for 20 genes the events_per_parameter guardrail fires
(4.2 < 5), so this demonstrates the pipeline, not a validated biomarker.
There is no independent breast-cohort validation here; external validation is
shown on an independent cohort in the LUAD narrative (case study 4).
vignette("case-study-confounder-audit") plants a batch effect and a
redundant, near-collinear design variable (Cramer's V = 0.82), plus a
survival signal that runs through the batch. The naive tier shows the six
batch-driven genes ranking as the most "prognostic" genes in an unadjusted
screen; the guarded tier shows design_audit() flagging both variables as
associated with the expression surrogate and recommending the minimal design
~ batch, after which the batch-artifact genes drop out of screening.
vignette("case-study-small-cohort") shows why the events_per_parameter
guardrail exists: on 30 samples / 23 events the naive reading quotes the mean
CV C = 0.74, while the fold-level table exposes a range of 0.50-1.00
(sd 0.14). The guardrail fires (4.6 < 5) rather than quietly reporting
unstable numbers.
A full live-GEO reproduction of a discovery-to-external cross-cohort attempt:
discovery in GSE31210 (LUAD, 35 deaths), external validation in GSE50081
(n = 128, 52 deaths). Mechanics pass end to end, but external C = 0.540 with
log-rank p = 0.266 — transfer was not demonstrated, reported as an honest
negative rather than a tool failure. The companion finding is the
screening-internal CV optimism: on this ~21k-gene panel even
survival-permuted data yields a null CV C ~ 0.8 (not 0.5), so CV concordance
must be read relative to a matched null and only validate_external() is
fully out-of-sample (see "Reading the CV concordance"). A power analysis
in the simulation study (validation/simulate_study.R) quantifies why a
35-event external cohort is structurally underpowered for weak real
signatures: transfer power
< 0.30 at ~35 events, and >= 0.80 only near ~300 events for a C ~ 0.65
signature.
The following case studies are in preparation and will ship as new vignettes. Each demonstrates a failure mode that rnaSentry's guardrails detect automatically at audit time. Results will be added once the analyses are finalized.
| Planned case study | Data (public GEO) | What it will demonstrate |
|---|---|---|
| Clinical-covariate trap (BRCA) | live GSE20685 via fetch_gse20685() (shipped as vignette("case-study-brca"); synthetic offline fallback) |
a naive 10-gene signature looks significant univariately but loses significance after adjusting for age + subtype; design_audit() flags the covariates before modeling |
| Batch catastrophe (identical-platform merge) | GSE31210 (n = 246, GPL570) + GSE30219 (n = 293, GPL570) | merging two same-platform studies creates a batch-only signal; ~50 "prognostic" genes track study/batch; design_audit() exposes the association before modeling |
| Cross-histology false transfer (LUAD → LUSC) | GSE30219 (single cohort, both histologies + survival) | a signature discovered on LUAD fails to transfer across the histology boundary, reported as an honest negative |
rnaSentry's niche is signature discovery with an enforced audit trail and honest failure modes: quality control, sample-identity checks, confounder detection, and explicit gating are first-class product features rather than side effects. Two Bioconductor neighbours are closest:
10.18129/B9.bioc.asuri) — de novo survival-marker
discovery from a SummarizedExperiment via subsampling glmnet + univariate
Cox. The closest methodological analogue; rnaSentry differs by design in its
sparse, directly inspectable linear score (univariate screen + joint Cox +
repeated stratified CV), its concordance/external-validation focus, and its
flag-ledger audit layer.10.18129/B9.bioc.signifinder) — applies 60+
published cancer signatures as single-sample scores (GSVA; single-cell and
spatial support). Complementary rather than competing: it answers "what do
known signatures say about this sample?" where rnaSentry answers "can a new
survival signature be learned from this cohort?". It can be used downstream
to benchmark a discovered signature.Two prior-art anchors for the single-sample scoring approach itself: SurvMarker (Gammune & Gu, 2025, DOI 10.64898/2025.12.31.697184) uses PCA-based weighted scoring, and the mRNAsi stemness index (Malta et al., Cell 2018;173:338-354.e15) uses OCLR-based scoring. Full profiles and a comparison table are in the validation dossier.
Contributions are welcome — see CONTRIBUTING.md for the reporting process,
the developer workflow, and the gate suite that every change must pass.
browseVignettes("rnaSentry") — vignette("getting-started")
bridges two CSV files into a SummarizedExperiment;
vignette("rnaSentry") is the main walkthrough;
vignette("case-study-impact") is the naive-vs-guarded landing page;
vignette("case-study-brca"),
vignette("case-study-confounder-audit"), and
vignette("case-study-small-cohort") are the three synthetic-cohort case
studies (BRCA additionally uses live GEO data).help(package = "rnaSentry")
validation/ directory of the
GitHub repository.Please file bugs and feature requests at https://github.com/Ghrieb/rnaSentry/issues.