1. Introduction to rnaSentry
2026-10-02
1 Introduction
rnaSentry automates the statistical pipeline used to go from a raw bulk
RNA-seq count matrix to a validated prognostic gene signature: quality
control and sample-identity auditing, confounder detection prior to
differential expression, signature construction with cross-validation,
survival modeling, and external validation. Every automated decision is
reported with the statistic that justified it. The pipeline is built on
core Bioconductor infrastructure — SummarizedExperiment for data
representation, BiocParallel for optional parallel cross-validation, and
GEOquery for public-cohort reuse — which motivates its
inclusion in Bioconductor over CRAN.
The five-stage workflow exposes eight exported functions:
load_counts() (intake audit), qc_explore() / sex_check() (quality
checks), design_audit() / pca_audit() (confounder detection), and
build_signature() / km_curve() / cox_model() /
survival_parametric() / validate_external() (linear Cox
survival::coxph + survival::survreg risk scores, repeated stratified
repeats/folds cross-validation with Harrell’s concordance
reverse = TRUE, and honest external validation). A flag ledger records
every guardrail decision (e.g. events_per_parameter, possibly_log_scaled,
redundant_variable) for the HTML audit report. Key points: the score is
sparse and inspectable (univariate screen + joint Cox), the cross-validated
concordance is screening-internal and must be calibrated against a null, and
the pipeline never corrects or removes batch effects — it detects and reports
them.
This vignette demonstrates the pipeline’s first stage: quality control and biological-sex verification.
Three companion case studies show the pipeline on real and adversarial data and
are meant to be read alongside this walkthrough:
vignette("case-study-brca") (the full pipeline on a breast-cancer cohort),
vignette("case-study-confounder-audit") (a planted batch confound), and
vignette("case-study-small-cohort") (why small cohorts fail silently).
2 Positioning and prior art
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:
- asuri (Bioc 3.23,
10.18129/B9.bioc.asuri) — de novo survival-marker discovery from aSummarizedExperimentvia subsampling glmnet + univariate Cox. The closest methodological analogue;rnaSentrydiffers 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. - signifinder (Bioc 3.16,
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?” wherernaSentryanswers “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.
| Tool | Scope | Scoring | Audit |
|---|---|---|---|
| rnaSentry | De novo prognostic signature, bulk RNA-seq | Sparse linear: univariate Cox + joint coxph, repeated CV concordance + external validate_external |
Full flag ledger, gated report |
| asuri (Bioc 3.23) | De novo marker via subsampling glmnet + Cox | Penalised glmnet + Cox | Subsampling stability, no ledger |
| signifinder (Bioc 3.16) | Apply 60+ published signatures as single-sample scores (GSVA) | GSVA | N/A (applies known signatures) |
| SurvMarker (2025) | Single-sample prognostic scoring | PCA-weighted | No audit trail |
| mRNAsi (2018) | Stemness single-sample score | OCLR | No survival model |
Full profiles are in the validation dossier.
3 Installation
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("rnaSentry")4 Design decisions
Two deliberate choices in the pipeline are worth stating explicitly, because they are easy to mistake for oversights.
First, build_signature() screens genes adjusted for the design terms
supplied via design_terms (by default), so that genes whose survival
association is driven entirely by a confounder such as batch are not
selected. The final signature weights, however, are estimated from a
gene-only model. This asymmetry is intentional: the risk score must remain
portable to external cohorts that do not record the discovery cohort’s
design covariates. Adjusted inference on the chosen signature is provided
separately by cox_model().
Second, pca_audit() reports both raw and Benjamini-Hochberg adjusted
p-values for its per-PC batch tests but flags on the raw p-value, whereas
design_audit() flags candidate confounders on BH-adjusted p-values. This,
too, is intentional: pca_audit() tests at most a handful of PCs against a
single batch variable, where a multiple-testing correction would be a near
no-op, while design_audit() scans many candidate confounders at once and
needs the correction to control spurious flags.
5 Data intake
load_counts() validates a raw count matrix and its sample metadata before
the pipeline runs. It returns a SummarizedExperiment and records any intake
issues in a flag ledger in the object’s metadata, so problems such as
duplicate sample identifiers or non-integer counts are surfaced up front and
carried into the rendered report’s audit trail.
library(rnaSentry)
library(SummarizedExperiment)
cnt <- matrix(rpois(600, lambda = 200), nrow = 60, ncol = 10,
dimnames = list(paste0("gene", 1:60), paste0("S", 1:10)))
cnt <- cnt * 1.0
cnt[1, ] <- 0L # gene1 has no counts anywhere (empty row)
cnt[2, 1] <- cnt[2, 1] + 0.5 # a non-integer value
colnames(cnt)[10] <- "S1" # a duplicate sample identifier
cd <- S4Vectors::DataFrame(condition = rep(c("tumor", "normal"), each = 5),
row.names = paste0("S", 1:10))
se_in <- load_counts(cnt, cd)
#> Warning: Duplicate sample IDs found: S1. Removing duplicates.
#> Warning: Count matrix contains non-integer values. Downstream stages expect
#> integer counts.
#> Warning: 1 gene(s) have all-zero or all-NA counts and will be retained but may
#> be filtered later.
S4Vectors::metadata(se_in)$flags
#> check severity detail
#> 1 duplicate_samples critical Duplicate sample IDs: S1.
#> 2 non_integer_counts warning Count matrix contains non-integer values.
#> 3 empty_rows warning 1 gene(s) have all-zero or all-NA counts.
#> stage
#> 1 load_counts
#> 2 load_counts
#> 3 load_countsEach row of the ledger is a (check, severity, detail, stage) record; the
duplicate sample is dropped (critical), while the non-integer and empty-row
issues are flagged as warnings and carried forward.
6 Quality control
qc_explore() takes a SummarizedExperiment of raw counts and audits it
for duplicate sample identifiers, missing metadata, non-integer count
values, and library-size outliers.
set.seed(1)
counts <- matrix(rpois(600, lambda = 200), nrow = 60, ncol = 10,
dimnames = list(paste0("gene", 1:60), paste0("S", 1:10)))
counts[, 1] <- counts[, 1] * 50L
coldata <- S4Vectors::DataFrame(condition = rep(c("tumor", "normal"), each = 5),
row.names = paste0("S", 1:10))
se <- SummarizedExperiment(assays = list(counts = counts), colData = coldata)
qc <- qc_explore(se)
qc
#> rnaSentry QC audit: 10 samples x 60 genes
#> 1 issue(s) flagged:
#> [warning] library_size_outlier: Samples with library size > 3 MADs from the median: S17 Biological sex verification
sex_check() compares reported sex in sample metadata against sex inferred
from XIST and Y-chromosome marker gene expression, to catch sample swaps
or metadata errors before they propagate downstream.
genes <- c("XIST", "RPS4Y1", "DDX3Y", "KDM5D")
expr <- matrix(c(950, 5, 5, 5, # sample looks female
5, 400, 380, 410), # sample looks male
nrow = 4, dimnames = list(genes, c("S1", "S2")))
coldata2 <- S4Vectors::DataFrame(sex = c("F", "F"), row.names = c("S1", "S2")) # S2 mislabeled
se2 <- SummarizedExperiment(assays = list(counts = expr), colData = coldata2)
sex_check(se2)
#> sample_id reported_sex inferred_sex status
#> 1 S1 F F OK
#> 2 S2 F M MISMATCH8 The survival pipeline in one run
run_rnaSentry() orchestrates the full pipeline on a single cohort:
confounder scan, PCA audit, signature construction with cross-validation,
Kaplan-Meier risk groups, an adjusted Cox model, and parametric survival
models. Here we build a synthetic cohort with an engineered survival
signal, split it into discovery and validation halves, and run the pipeline
on the discovery half.
set.seed(1)
n_genes <- 120
n_samples <- 120
counts <- matrix(stats::rpois(n_genes * n_samples, lambda = 500),
nrow = n_genes, ncol = n_samples,
dimnames = list(paste0("gene", 1:n_genes),
paste0("S", 1:n_samples)))
sig_expr <- colMeans(counts[1:5, , drop = FALSE])
risk <- scale(sig_expr)[, 1] * 0.4
event_time <- stats::rexp(n_samples, rate = 0.03 * exp(0.8 * risk))
censor_time <- stats::rexp(n_samples, rate = 0.02)
coldata <- S4Vectors::DataFrame(time = pmin(event_time, censor_time),
event = as.integer(event_time < censor_time),
condition = factor(rep(c("A", "B"), length.out = n_samples)),
batch = factor(rep(c("B1", "B2", "B3"), length.out = n_samples)),
row.names = colnames(counts))
se_all <- SummarizedExperiment(assays = list(counts = counts), colData = coldata)
discovery <- se_all[, 1:72]
external <- se_all[, 73:120]
run <- run_rnaSentry(discovery, time_col = "time", event_col = "event",
top_n = 5, repeats = 2, folds = 3, seed = 7,
design_vars = c("condition", "batch"),
design_terms = "condition",
render_report = FALSE)
run
#> rnaSentry pipeline run.
#> - design_audit
#> - pca_audit
#> - build_signature
#> - km_curve
#> - cox_model
#> - survival_parametricThe run locked the discovered signature: run_rnaSentry() is single-use per
session, and a second call in the same session is refused until the lock is
released. This is the reproducibility guardrail against silently
re-selecting signature genes after survival analysis.
# Release the lock so the pipeline can be run again in this session.
run$stages$build_signature <-
lock_signature(run$stages$build_signature, lock = FALSE)The result’s stages list holds each stage’s full result, including the
flags data.frame that records every issue raised:
run$stages$build_signature$flags
#> check severity
#> 1 adjusted_screening info
#> 2 signature_locked info
#> 3 signature_unlocked info
#> detail stage
#> 1 Univariate screening adjusted for design term(s): condition. build_signature
#> 2 Signature locked against downstream mutation. lock_signature
#> 3 Signature unlocked; downstream stages may be rerun. lock_signature8.1 External validation
The discovery cutpoint is recorded by the km_curve stage. Re-scoring the
held-out samples against that same cutpoint, without re-estimating anything
from the external data, gives an honest portability check:
val <- validate_external(run$stages$build_signature, external,
cutpoint = run$stages$km_curve$cutpoint)
val
#> rnaSentry external validation: 5/5 signature gene(s) used, 48 samples.
#> Risk groups: low n = 23 (12 events), high n = 25 (11 events).
#> Cutpoint 55.5 (custom, from discovery).
#> Log-rank p = 0.8303. Concordance = 0.485.
#> Signature not locked.The portability check is also visual: scoring the held-out samples against the locked discovery cutpoint, without re-estimating anything from the external data, splits the external cohort into risk groups whose curves can be compared on truly held-out data:
plot(val)
Figure 1: External-cohort Kaplan-Meier curves for the low- and high-risk groups, scored against the locked discovery cutpoint without re-estimation
The log-rank p-value tests the split on held-out data.
On this synthetic external half the split is not statistically separated (log-rank p and concordance printed above): with 48 held-out samples the toy cohort is far too small to demonstrate transfer. The honest-negative pattern is the point – an out-of-sample split that does not separate is reported as such, and a real transfer claim requires an independent cohort with enough events to detect the effect.
A real discovery-to-external transfer on independent GEO cohorts — where transfer was not demonstrated and the reasons (low event count, CV-optimism) are quantified — is documented in the README’s case study 4 (LUAD, GSE31210 -> GSE50081).
9 Assumptions and limitations
rnaSentry is deliberately narrow. Read this before trusting any output.
- Bulk, right-censored survival only. The pipeline expects bulk RNA-seq and standard right-censored (time, event) end points. Single-cell data, continuous or binary outcomes, competing risks, time-varying covariates, and left truncation are out of scope.
- Linear Cox risk scores. A signature is a weighted linear combination of gene expression. Higher scores mean higher hazard (shorter survival); all concordance estimates use this convention.
- No batch correction, no differential expression. Confounders are detected and reported, and univariate screening can be adjusted for design terms, but the pipeline never corrects or removes effects and makes no DE calls.
- The small-cohort trap. With fewer than ~5 events per signature gene,
cross-validated concordance and hazard ratios are unstable.
build_signature()warns (events_per_parameter); a flagged signature should be treated as hypothesis-generating, not confirmatory. - The CV concordance is screening-internal. Gene selection runs on the
full cohort before the CV split, so
build_signature()’s concordance is optimistic: on a large panel even survival-permuted data yields a high null (≈0.8 on ~21k genes) rather than 0.5. Read it relative to a matched null (permutation or random-gene control), and rely onvalidate_external()on an independent cohort as the fully out-of-sample estimate. - Right input scale. Raw counts are log2-transformed internally. If you
store already-logged data under a name like
counts, it will be double-logged; the pipeline warns (possibly_log_scaled). Store pre-scaled data under the namelogcounts. - One discovery per session.
run_rnaSentry()locks the discovered signature and refuses to re-run untillock_signature(sig, lock = FALSE)releases the lock. This prevents silent gene-set re-selection after survival analysis. To analyze a second cohort in the same session, unlock the first signature first.
9.1 Avoiding common mistakes
Two misuses are easy to hit and easy to avoid.
First, wrong input scale. If your expression matrix is already
log-transformed but stored under a counts-named assay, every downstream
stage sees a double-logged matrix. The pipeline catches this and warns:
prelogged <- log2(as.matrix(assay(se_all, "counts")) + 1)
se_wrong <- se_all
assay(se_wrong, "counts") <- prelogged
build_signature(se_wrong, "time", "event", top_n = 5, repeats = 1,
folds = 2, seed = 7, min_events_per_parameter = 1)
#> Warning: The first assay looks already log-transformed (non-integer values with
#> a maximum below 40). rnaSentry will log2-transform it again, double-logging
#> expression. Rename the assay to "logcounts" (RNA-seq) or "exprs" (microarray),
#> or supply raw integer counts.
#> rnaSentry signature: 5 gene(s), outcome 'overall_survival'.
#> Selection: top_n (5 gene(s)). CV concordance: mean 0.626, sd 0.014.
#> Not locked.
#> 1 issue(s) flagged:
#> [warning] possibly_log_scaled: The first assay looks already log-transformed (non-integer values with a maximum below 40). rnaSentry will log2-transform it again, double-logging expression. Rename the assay to "logcounts" (RNA-seq) or "exprs" (microarray), or supply raw integer counts.The fix is to name such an assay logcounts, which the pipeline uses as-is
without a second transform.
Second, underpowered discovery. On a small cohort, build_signature()
warns that the event count per signature gene is low:
small <- se_all[, 1:40]
build_signature(small, "time", "event", top_n = 10, repeats = 1,
folds = 2, seed = 7)
#> Warning: Only 22 event(s) for a 10-gene signature (2.2 events per parameter),
#> below 'min_events_per_parameter' = 5. Cross-validated concordance and hazard
#> ratios are unstable at this event count; consider fewer genes or a larger
#> cohort.
#> rnaSentry signature: 10 gene(s), outcome 'overall_survival'.
#> Selection: top_n (10 gene(s)). CV concordance: mean 0.647, sd 0.073.
#> Not locked.
#> 1 issue(s) flagged:
#> [warning] events_per_parameter: Only 22 event(s) for a 10-gene signature (2.2 events per parameter), below 'min_events_per_parameter' = 5. Cross-validated concordance and hazard ratios are unstable at this event count; consider fewer genes or a larger cohort.Both warnings are also recorded in each stage’s flags data.frame, so they
propagate into the rendered report’s audit trail.
10 Session info
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] survival_3.8-12 SummarizedExperiment_1.43.0
#> [3] Biobase_2.73.2 GenomicRanges_1.65.4
#> [5] Seqinfo_1.3.2 IRanges_2.47.5
#> [7] S4Vectors_0.51.10 BiocGenerics_0.59.12
#> [9] generics_0.1.4 MatrixGenerics_1.25.0
#> [11] matrixStats_1.5.0 rnaSentry_0.99.19
#> [13] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 ggplot2_4.0.3 xfun_0.61
#> [4] bslib_0.12.0 httr2_1.3.0 lattice_0.23-1
#> [7] tzdb_0.5.0 vctrs_0.7.3 tools_4.6.1
#> [10] curl_8.0.0 tibble_3.3.1 pkgconfig_2.0.3
#> [13] R.oo_1.27.1 Matrix_1.7-6 data.table_1.18.6.1
#> [16] RColorBrewer_1.1-3 rentrez_1.2.4 S7_0.2.2
#> [19] lifecycle_1.0.5 farver_2.1.2 compiler_4.6.1
#> [22] tinytex_0.61 statmod_1.5.2 htmltools_0.5.9
#> [25] sass_0.4.10 yaml_2.3.12 pillar_1.11.1
#> [28] jquerylib_0.1.4 tidyr_1.3.2 R.utils_2.13.0
#> [31] DelayedArray_0.39.8 cachem_1.1.0 limma_3.99.0
#> [34] magick_2.9.1 abind_1.4-8 tidyselect_1.2.1
#> [37] digest_0.6.39 dplyr_1.2.1 purrr_1.2.2
#> [40] bookdown_0.48 labeling_0.4.3 splines_4.6.1
#> [43] fastmap_1.2.0 grid_4.6.1 GEOquery_2.81.28
#> [46] cli_3.6.6 SparseArray_1.13.4 magrittr_2.0.5
#> [49] S4Arrays_1.13.2 dichromat_2.0-1 XML_3.99-0.25
#> [52] readr_2.2.0 withr_3.0.3 scales_1.4.0
#> [55] rmarkdown_2.32 XVector_0.53.0 otel_0.2.0
#> [58] R.methodsS3_1.8.2 hms_1.1.4 evaluate_1.0.5
#> [61] knitr_1.52 rlang_1.3.0 Rcpp_1.1.2
#> [64] glue_1.8.1 BiocManager_1.30.27 xml2_1.6.0
#> [67] jsonlite_2.0.0 R6_2.6.1