4. Case study: breast cancer (GSE20685)
2026-10-02
1 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.
2 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':
#>
#> Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append,
#> as.data.frame, basename, cbind, colnames, dirname, do.call,
#> duplicated, eval, evalq, get, grep, grepl, is.unsorted, lapply,
#> mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#> rank, rbind, 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':
#>
#> I, expand.grid, 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] 83The 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.
3 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.000000survival::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.
4 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_signatureThe 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 assertsC + C_rev = 1on every fold). - Events per parameter. The
events_per_parameterwarning 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.
4.1 Risk groups
km <- run$stages$km_curve
km$log_rank_p
#> [1] 3.261912e-17The discovered signature separates the cohort into risk groups whose survival curves differ dramatically at the discovery cutpoint:
plot(run$stages$km_curve)
Figure 1: 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.
4.2 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)
Figure 2: 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).
4.3 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)
Figure 3: Kaplan-Meier curves with parametric survival-model overlay for the low- and high-risk groups
4.4 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")
Figure 4: 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.
5 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.
6 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
exprsare 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.
7 References
Appendix
- 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()).
A 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