Status note

This case study is an illustrative walkthrough of rnaSentry on the public GSE20685 breast-cancer dataset, provided to show how the pipeline behaves on real data. It is not peer-reviewed and the signature reported here is provisional. Treat every number below as a demonstration, not a clinical claim.

The narrative follows the three-tier structure used across the case studies: a naive analysis first, then rnaSentry standing guard, then the counterfactual – what a researcher would have concluded without the guardrails, and what the guarded analysis actually recovers.

Data

We use GSE20685 (Li et al., 2010), an Affymetrix GPL570 (U133 Plus 2.0) microarray series of 327 primary breast tumors with overall-survival follow-up. When GEOquery is available, the full series matrix is downloaded from GEO, processed to gene-level expression (probe-to-gene collapse by largest mean expression), and subset to the top 3000 most variable genes. The result is a SummarizedExperiment with an exprs assay (microarray log2 intensities, used as-is) plus time (years), event (1 = death), age, and subtype in the sample metadata. GEOquery caches the series matrix itself, so repeat builds reuse the local copy. Only if GEOquery is unavailable or the download fails, an internal synthetic fallback of identical structure (.brca_fallback_cohort(), fixed seed for reproducibility) is used instead so the vignette still builds offline.

library(rnaSentry)
library(SummarizedExperiment)
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#> 
#> Attaching package: 'MatrixGenerics'
#> The following objects are masked from 'package:matrixStats':
#> 
#>     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#>     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#>     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#>     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#>     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#>     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#>     colWeightedMeans, colWeightedMedians, colWeightedSds,
#>     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#>     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#>     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#>     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#>     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#>     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#>     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#>     rowWeightedSds, rowWeightedVars
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: generics
#> 
#> Attaching package: 'generics'
#> The following objects are masked from 'package:base':
#> 
#>     as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#>     setequal, union
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following object is masked from 'package:utils':
#> 
#>     data
#> The following objects are masked from 'package:base':
#> 
#>     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#>     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#>     get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
#>     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
#>     rbind, Reduce, rownames, sapply, saveRDS, scale, sequence, table,
#>     tapply, transform, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#> 
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#> 
#>     findMatches
#> The following objects are masked from 'package:base':
#> 
#>     expand.grid, I, unname
#> Loading required package: IRanges
#> Loading required package: Seqinfo
#> Loading required package: Biobase
#> Welcome to Bioconductor
#> 
#>     Vignettes contain introductory material; view with
#>     'browseVignettes()'. To cite Bioconductor, see
#>     'citation("Biobase")', and for packages 'citation("pkgname")'.
#> 
#> Attaching package: 'Biobase'
#> The following object is masked from 'package:MatrixGenerics':
#> 
#>     rowMedians
#> The following objects are masked from 'package:matrixStats':
#> 
#>     anyMissing, rowMedians
library(survival)

se_real <- fetch_gse20685()
#> Setting options('download.file.method.GEOquery'='auto')
#> Setting options('GEOquery.inmemory.gpl'=FALSE)
#> Downloading GSE20685 from GEO ...
#> Found 1 file(s)
#> GSE20685_series_matrix.txt.gz
#> getGEO() now returns SummarizedExperiment objects by default. Pass returnType = 'ExpressionSet' for the previous behavior.
#> This message is displayed once per session.
using_real <- !is.null(se_real)
if (using_real) {
  se <- se_real
} else {
  message("GEOquery unavailable or download failed; using synthetic fallback.")
  se <- rnaSentry:::.brca_fallback_cohort()
}
se
#> class: SummarizedExperiment 
#> dim: 3000 327 
#> metadata(0):
#> assays(1): exprs
#> rownames(3000): CRISP3 S100A7 ... HYDIN MAGEB3
#> rowData names(0):
#> colnames(327): GSM519117 GSM519118 ... GSM519442 GSM519443
#> colData names(4): time event age subtype
table(se$subtype)
#> 
#>   type I  type II type III  type IV   type V  type VI 
#>       37       34       41       81       41       93
sum(se$event)
#> [1] 83

The cohort analyzed below is the live GSE20685 series matrix from GEO (downloaded via GEOquery, which caches the series matrix locally). It records 83 deaths, so every number printed in this vignette reproduces whichever source was used.

Tier 1: the naive analysis

A standard, unguarded analysis of this cohort would proceed along these lines: pick a few highly variable genes, screen them against survival with univariate Cox regressions, average the survivors into a risk score, and quote the concordance that R prints by default. Let us run exactly that.

lg <- assays(se)$exprs
vars <- apply(lg, 1, stats::var, na.rm = TRUE)
cand <- names(sort(vars, decreasing = TRUE))[seq_len(200)]
stat <- do.call(rbind, lapply(cand, function(g) {
    s <- summary(coxph(Surv(se$time, se$event) ~ lg[g, ]))
    c(beta = s$coefficients[1, 1], p = s$coefficients[1, 5])
}))
rownames(stat) <- cand
top <- stat[order(stat[, "p"]), , drop = FALSE][seq_len(5), , drop = FALSE]
naive_score <- colMeans(sign(top[, "beta"]) * lg[rownames(top), , drop = FALSE])

c_default <- as.numeric(concordance(Surv(se$time, se$event) ~
                                    naive_score)$concordance[1])
c_reverse <- as.numeric(concordance(Surv(se$time, se$event) ~ naive_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.323754  0.676246  1.000000

survival::concordance() defaults to the convention larger score means longer survival (reverse = FALSE). A Cox risk score means the opposite – larger score means shorter survival – so the default read is inverted: C_default = 1 - Harrell's C. The naive analyst sees C_default < 0.5 and concludes “below chance, poor generalization, the signature is a failure.”

This is not hypothetical. The first real-data face-validity run of this very cohort reported CV concordance 0.217 and a naive reading would rationalize it as extreme selection-induced winner’s curse; the correct value is 1 - 0.217 = 0.783, a directionally correct held-out result (see Limitations on what CV alone can and cannot show). The full diagnosis, including the C + C_rev = 1 invariant you can see holding above, is recorded in the repository’s validation dossier. The parity test suite now asserts that invariant on every fold so the convention cannot silently flip again.

Even read in the correct direction, a naive analyst stops at the cross- validated mean and treats it as an out-of-sample estimate. As the LUAD honest-negative case study documents, a screening-internal CV on a large candidate panel is optimistic: even survival-permuted data yields a null CV around 0.8 on ~21k genes, not 0.5. Any concordance below must be read relative to a matched null, and the only fully out-of-sample number is validate_external() on an independent cohort.

Tier 2: rnaSentry standing guard

Now the same cohort through the guarded pipeline: confounder scan, PCA audit, signature construction with repeated cross-validation, and Kaplan-Meier risk groups. Because the cohort metadata records age and subtype, both are supplied as design variables so the audit can check them against the expression surrogate.

run <- run_rnaSentry(se, time_col = "time", event_col = "event",
                     outcome_col = "overall_survival",
                     design_vars = c("age", "subtype"),
                     top_n = 20, repeats = 3, folds = 3, seed = 42,
                     render_report = FALSE)
#> Warning: Only 83 event(s) for a 20-gene signature (4.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.
sig <- run$stages$build_signature
head(sig$genes)
#> [1] "C20orf96" "SPRR3"    "EPHX2"    "CD1E"     "CASP14"   "ARNT2"
sig$cv_summary
#>       mean         sd 
#> 0.78469786 0.03000978
sig$flags
#>                  check severity
#> 1 events_per_parameter  warning
#> 2     signature_locked     info
#>                                                                                                                                                                                                                                detail
#> 1 Only 83 event(s) for a 20-gene signature (4.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.
#> 2                                                                                                                                                                                       Signature locked against downstream mutation.
#>             stage
#> 1 build_signature
#> 2  lock_signature

The run locks the discovered signature, so a second discovery in the same session is refused until the lock is released (see the main vignette). Because all vignettes build in one session, we release it here so the remaining vignettes can run their own pipeline calls:

run$stages$build_signature <-
    lock_signature(run$stages$build_signature, lock = FALSE)

The guardrails that the naive analysis missed are now explicit in the flag ledger:

Direction. Every concordance in the pipeline is computed with reverse = TRUE; the sign convention cannot silently flip (a parity test asserts C + C_rev = 1 on every fold).
Events per parameter. The events_per_parameter warning is expected and matches the documented limitation: 83 deaths for a 20-gene signature is 4.2 events per parameter, below the guardrail of 5. The signature should be read as hypothesis-generating, and the cross-validated concordance and hazard ratios as unstable at this event count.
Screening-internal CV. The cross-validated concordance (mean C = 0.78) is a screening-internal estimate: gene selection happens on the full cohort before the CV split, so it is optimistic. Read it relative to a matched null, as discussed in the main vignette and in the README.
Design. The design audit below decides which confounders the screening adjusts for, instead of an ad-hoc pick.

Risk groups

km <- run$stages$km_curve
km$log_rank_p
#> [1] 3.261912e-17

The discovered signature separates the cohort into risk groups whose survival curves differ dramatically at the discovery cutpoint:

plot(run$stages$km_curve)

Kaplan-Meier survival curves for the low- and high-risk groups at the discovery cutpoint, with the log-rank p-value annotated. The risk groups diverge sharply.

Adjusted hazard ratios

The adjusted Cox model re-estimates each signature gene’s association with survival jointly with the other signature genes, so the forest plot below shows which genes carry their effect even after adjustment for the rest of the signature:

plot(run$stages$cox_model)

Adjusted Cox forest plot for the signature genes: log hazard ratio with 95% confidence interval, ordered by p-value. The vertical dashed line marks no effect (log HR = 0).

Parametric survival model

The parametric stage fits a parametric survival model (exponential and Weibull) to the same risk groups and overlays it on the Kaplan-Meier curves, providing a smooth estimate of the survival function when a parametric shape is a reasonable description:

plot(run$stages$survival_parametric)

Kaplan-Meier curves with parametric survival-model overlay for the low- and high-risk groups.

Design audit

da <- run$stages$design_audit
da$confounder_table
#>   variable        type   n test   statistic     df       p_value effect_size
#> 1      age     numeric 327   lm   0.7932527 1, 325  3.737772e-01 0.002434835
#> 2  subtype categorical 327  aov 245.7217644 5, 321 2.008554e-107 0.792850947
#>   effect_size_type flagged         adj_p
#> 1        r_squared   FALSE  3.737772e-01
#> 2      eta_squared    TRUE 4.017108e-107
da$formula_text
#> [1] "~ age + subtype"

subtype is flagged as associated with the expression surrogate (very large eta-squared) and the recommended design formula therefore keeps it, so that univariate gene screening is adjusted for subtype. age is not flagged. This is exactly the kind of design decision rnaSentry is meant to surface before any differential-expression or signature work.

Coloring the PCA audit scores by subtype shows the same structure the surrogate scan detects: the breast-cancer subtypes occupy distinct regions of expression space, which is why omitting subtype from the design would let subtype-driven expression masquerade as survival association:

plot_pca_audit(run$stages$pca_audit, color_by = "subtype")

PC1-PC2 scatter of the GSE20685 subset colored by breast-cancer subtype. The subtypes separate along the leading PCs, mirroring the design audit’s surrogate-association flag for subtype.

Tier 3: the counterfactual

	Naive analysis (Tier 1)	rnaSentry (Tier 2)
Score direction	default convention reads the score as protective; `C < 0.5` looks like failure	`reverse = TRUE`; `C_rev = 1 - C_default`, direction correct
Headline concordance	`C_default` quoted as gospel	`C_rev` quoted with its convention, and flagged screening-internal
Out-of-sample claim	none / implicit	only `validate_external()` counts; matched-null calibration required
Events per parameter	not reported	4.2 < 5, warned in the ledger
Design	ad-hoc	`subtype` flagged and adjusted for

What a guarded analysis can claim: on the full-cohort face-validity run of this dataset rnaSentry reports a directionally correct held-out concordance of 0.783 against a 0.582 random-gene control (and a 0.490 simulated null) – while the naive pipeline either reports a “failed” signature (inverted convention) or an uncalibrated 0.80 (screening-internal optimism). Caveat (see Limitations): the rerun CV of 0.785 sits inside, not above, a matched permuted-selection noise band (0.744–0.785), so the CV alone is not evidence of discrimination – the signal evidence is the hazard directions and the risk split, and transfer remains unvalidated here. The events-per-parameter limitation is stated, not hidden.

Limitations

No independent breast cohort is validated here. External validation is demonstrated on an independent cohort in the LUAD honest-negative narrative in the README and in the main vignette; this walkthrough deliberately stops at in-cohort analysis so the numbers stay interpretable.
Screening-internal CV. The concordance above is optimistic by design; see the main vignette’s “Assumptions and limitations”.
Internal CV does not clearly beat selection-driven noise. Identical-protocol builds on survival-permuted data (6 draws) give CV C 0.744–0.785 vs this signature’s 0.785 — the real estimate sits inside, not above, the noise band. The CV number alone therefore cannot clearly distinguish genuine signal from screening optimism in this cohort; evidence for real signal rests on MKI67/ESR1 hazard directions matching known biology and risk-group separation (log-rank p = 3.3e-17), and only external validation can confirm the multi-gene signature generalizes.
Low events per parameter. The guardrail warning above is a real signal; this signature is a demonstration, not a published result.
GSE20685 is an Affymetrix series; its exprs are already log2 intensities, used as-is (the microarray assay name, not log-counts).
GPL annotation drift. Gene-level counts may shift slightly across sessions as NCBI updates the GPL570 annotation (23,520 post-collapse genes in August 2026 vs 22,189 in September 2026); core findings — MKI67/ESR1 directions, risk-group separation (median low not reached / high 11.4 y), CV ≈ 0.78 — reproduced across both annotations.

References

Li Y, Zou L, Li Q, et al. A five-gene molecular grade index and HOXB13:IL17BR are complementary prognostic factors in early stage breast cancer. Clinical Cancer Research (2010). GEO accession: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE20685 (Affymetrix GPL570, HG-U133 Plus 2.0; downloaded here via GEOquery::getGEO()).

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

4. Case study: breast cancer (GSE20685)