The case studies are written as three acts. Tier 1 runs the analysis
the way a standard, unguarded pipeline would, and prints the headline a
researcher would have believed. Tier 2 runs the same data
through rnaSentry and shows which guardrail fired, with the
statistic that justified it. Tier 3 states the counterfactual: the error
the guardrails prevented, and the true signal the guarded analysis still
recovered.
Four failure modes – the ways a prognostic-signature study most often goes wrong – are covered by the case studies:
| # | Failure mode | Where it is demonstrated |
|---|---|---|
| 1 | Wrong direction: an inverted concordance convention turns a strong signature into a “failure” | vignette("case-study-brca") |
| 2 | Confounded design: batch-driven genes masquerade as prognostic; collinear terms are double-counted | vignette("case-study-confounder-audit") |
| 3 | Underpowered discovery: a mean C-index hides coin-flip folds | vignette("case-study-small-cohort") |
| 4 | False transfer: screening-internal CV optimism is mistaken for cross-cohort validity | README / main vignette (LUAD honest negative, below) |
| Failure mode | Naive headline | Guarded headline | Guardrails that fired | What a guarded analysis can claim |
|---|---|---|---|---|
| Wrong direction (GSE20685) | “CV C = 0.217, signature fails to generalize (winner’s curse)” | held-out C = 0.783, direction correct, far above the 0.582 random-gene control | reverse = TRUE convention; C + C_rev = 1
parity tests; direction-sensitive regression tests |
the signature is not a 0.22 collapse but a directionally correct held-out signal (CV alone does not exceed selection-noise controls; see BRCA Limitations) |
| Confounded design (synthetic batch cohort) | “six batch-driven genes are prognostic; batch and region are two independent risk factors” | batch/region flagged redundant (Cramer’s V
= 0.82); batch surrogate eta-squared 0.97; minimal design
~ batch; adjusted screening drops the batch-artifact
genes |
redundant_variable, surrogate-association flags,
adjusted_screening |
the discovered signature is not a batch echo |
| Underpowered discovery (30 samples, 23 events) | “cross-validated C = 0.74 – decent” | fold range 0.50-1.00, sd 0.14, 4.6 events/parameter < 5 | events_per_parameter, fold-level reporting |
the concordance is an average over coin-flip folds; hypothesis-generating only |
| False transfer (GSE31210 -> GSE50081, LUAD) | “CV C = 0.86, the signature is validated” | CV C = 0.862 vs permutation null 0.836 (delta +0.026); external C = 0.540, log-rank p = 0.266 – transfer not demonstrated | permutation-null calibration; honest-negative external validation; power analysis | no false claim of cross-cohort validity; the negative is a power finding (35 events), not a tool failure |
A sign-convention bug is the one error class that flips a paper’s
conclusion without any red flag in the data.
survival::concordance() defaults to larger score means
longer survival; a Cox risk score means the opposite. Calling the
default on a risk score returns 1 - Harrell's C. The
identity below is what the parity test suite now asserts on every
fold:
library(rnaSentry)
library(SummarizedExperiment)
library(survival)
# Synthetic cohort for the sign-convention demo (real GEO via GEOquery
# shown in case-study-brca; synthetic used here for offline builds).
# NOTE: seed deliberately differs from any GEO accession so future
# `git log -S "set.seed(20685)"` audits stay unambiguous.
set.seed(42)
n_genes <- 3000
n_samples <- 327
genes <- paste0("GENE", seq_len(n_genes))
samples <- paste0("Sample", seq_len(n_samples))
mat <- matrix(rnorm(n_genes * n_samples, mean = 6, sd = 1.5),
nrow = n_genes, ncol = n_samples,
dimnames = list(genes, samples))
risk0 <- colMeans(mat[1:5, , drop = FALSE])
event_time <- rexp(n_samples, rate = 0.08 * exp(0.6 * scale(risk0)[, 1]))
censor_time <- rexp(n_samples, rate = 0.04)
time <- pmin(event_time, censor_time)
event <- as.integer(event_time < censor_time)
coldata <- S4Vectors::DataFrame(time = time, event = event,
row.names = samples)
se <- SummarizedExperiment(assays = list(logcounts = mat), colData = coldata)
score <- colMeans(assays(se)$logcounts[1:5, , drop = FALSE])
c_default <- as.numeric(concordance(Surv(se$time, se$event) ~
score)$concordance[1])
c_reverse <- as.numeric(concordance(Surv(se$time, se$event) ~ score,
reverse = TRUE)$concordance[1])
c(c_default = c_default, c_reverse = c_reverse, sum = c_default + c_reverse)
#> c_default c_reverse sum
#> 0.3494893 0.6505107 1.0000000A naive reading of the default concordance convention reports 0.349
and would rationalize it as winner’s curse;
1 - 0.349 = 0.651 is the correct, biologically expected
value all along. The same class of error has appeared in published
survival literature; rnaSentry pins the convention in code, in parity
tests, and in the rendered report.
Running rnaSentry on the synthetic cohort above
demonstrates the guardrail flags in action:
run <- run_rnaSentry(se, time_col = "time", event_col = "event",
top_n = 5, repeats = 2, folds = 3, seed = 42,
render_report = FALSE)
run$stages$build_signature$flags
#> check severity detail
#> 1 signature_locked info Signature locked against downstream mutation.
#> stage
#> 1 lock_signatureEvery flag in the ledger is backed by a statistic. The
events_per_parameter warning, the
reverse = TRUE convention, and the fold-level concordance
distribution are the guardrails that prevent the three failure modes
demonstrated in the case studies above.
The fourth failure mode is demonstrated on independent GEO cohorts in
the README and main vignette; the numbers are reproduced via
validation/repro_live_geo.R in the GitHub repository.
Because the data are downloaded live, the walkthrough is narrative here
rather than a build-time vignette:
validation/simulate_study.R) quantifies why: transfer
power < 0.30 at ~35 discovery events, and >= 0.80 only near ~300
events for a C ~ 0.65 signature.The same pattern – high screening-internal CV but poor external transfer – can be reproduced on a small synthetic cohort:
set.seed(7)
n_genes <- 120
n_samples <- 80
cnt <- matrix(rpois(n_genes * n_samples, lambda = 500),
nrow = n_genes, ncol = n_samples,
dimnames = list(paste0("gene", seq_len(n_genes)),
paste0("S", seq_len(n_samples))))
sig_expr <- colMeans(cnt[1:5, , drop = FALSE])
risk <- scale(sig_expr)[, 1] * 0.4
event_time <- rexp(n_samples, rate = 0.03 * exp(0.8 * risk))
censor_time <- rexp(n_samples, rate = 0.02)
cd <- S4Vectors::DataFrame(time = pmin(event_time, censor_time),
event = as.integer(event_time < censor_time),
row.names = colnames(cnt))
se_null <- SummarizedExperiment(assays = list(counts = cnt), colData = cd)
discovery <- se_null[, seq_len(60)]
external <- se_null[, 61:80]
sig_null <- build_signature(discovery, "time", "event", top_n = 5,
repeats = 2, folds = 3, seed = 7)
cat("Screening-internal CV C:", mean(sig_null$cv_results$c_index), "\n")
#> Screening-internal CV C: 0.7235437
km_null <- km_curve(sig_null, discovery)
val_null <- validate_external(sig_null, external,
cutpoint = km_null$cutpoint)
cat("External C:", val_null$concordance, "\n")
#> External C: 0.4341085
cat("External log-rank p:", val_null$log_rank_p, "\n")
#> External log-rank p: 0.1220304The screening-internal CV is optimistic (high C) but the external validation on the held-out half is poor – exactly the honest-negative pattern the tool is designed to surface rather than hide.
GEOquery (which caches the series matrix locally), with an
evaluated synthetic fallback of identical structure if the download is
unavailable, and the live-GEO LUAD numbers are reproduced by
validation/repro_live_geo.R in the GitHub repository.(check, severity, detail, stage)
record backed by a statistic, and each has a regression test asserting
the condition that fires it.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