1. Introduction to rnaSentry

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).

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 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.
  • 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?” 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.

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.

Installation

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("rnaSentry")

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.

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_counts

Each 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.

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: S1

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 MISMATCH

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_parametric

The 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_signature

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)
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.

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).

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 on validate_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 name logcounts.
  • One discovery per session. run_rnaSentry() locks the discovered signature and refuses to re-run until lock_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.

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.

Session info

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [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: Etc/UTC
#> 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        R.oo_1.27.1        
#> [13] pkgconfig_2.0.3     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] statmod_1.5.2       htmltools_0.5.9     sys_3.4.3          
#> [25] buildtools_1.0.0    sass_0.4.10         yaml_2.3.12        
#> [28] pillar_1.11.1       jquerylib_0.1.4     tidyr_1.3.2        
#> [31] R.utils_2.13.0      DelayedArray_0.39.8 cachem_1.1.0       
#> [34] limma_3.99.0        abind_1.4-8         tidyselect_1.2.1   
#> [37] digest_0.6.39       dplyr_1.2.1         purrr_1.2.2        
#> [40] labeling_0.4.3      maketools_1.3.2     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     XML_3.99-0.25       readr_2.2.0        
#> [52] withr_3.0.3         scales_1.4.0        rmarkdown_2.32     
#> [55] XVector_0.53.0      otel_0.2.0          R.methodsS3_1.8.2  
#> [58] hms_1.1.4           evaluate_1.0.5      knitr_1.52         
#> [61] rlang_1.3.0         glue_1.8.1          BiocManager_1.30.27
#> [64] xml2_1.6.0          jsonlite_2.0.0      R6_2.6.1