--- title: "4. Case study: breast cancer (GSE20685)" author: "Abdelkarim Hani Ghrieb" date: "`r Sys.Date()`" output: BiocStyle::html_document: toc: true toc_float: true vignette: > %\VignetteIndexEntry{4. Case study: breast cancer (GSE20685)} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r brca-setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") ``` # 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. ```{r load-data} library(rnaSentry) library(SummarizedExperiment) library(survival) se_real <- fetch_gse20685() 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 table(se$subtype) sum(se$event) ``` The cohort analyzed below is `r if (using_real) "the live GSE20685 series matrix from GEO (downloaded via GEOquery, which caches the series matrix locally)" else "a synthetic fallback cohort with the same structure (3000 genes x 327 samples, same exprs/time/event/age/subtype layout), used because the GEO download was unavailable"`. It records `r sum(se$event)` 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. ```{r naive-screen, warning = FALSE} 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) ``` `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. ```{r guarded-run, warning = TRUE} 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) sig <- run$stages$build_signature head(sig$genes) sig$cv_summary sig$flags ``` 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: ```{r unlock-signature} 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: `r sum(se$event)` deaths for a 20-gene signature is `r round(sum(se$event) / 20, 1)` 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 = `r round(sig$cv_summary[["mean"]], 2)`) 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 ```{r risk-groups} km <- run$stages$km_curve km$log_rank_p ``` The discovered signature separates the cohort into risk groups whose survival curves differ dramatically at the discovery cutpoint: ```{r brca-km, fig.width = 7, fig.height = 5, fig.align = "center", fig.cap = "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."} plot(run$stages$km_curve) ``` ## 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: ```{r brca-forest, fig.width = 7, fig.height = 5, fig.align = "center", fig.cap = "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)."} plot(run$stages$cox_model) ``` ## 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: ```{r brca-parametric, fig.width = 7, fig.height = 5, fig.align = "center", fig.cap = "Kaplan-Meier curves with parametric survival-model overlay for the low- and high-risk groups."} plot(run$stages$survival_parametric) ``` ## Design audit ```{r design-audit} da <- run$stages$design_audit da$confounder_table da$formula_text ``` `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: ```{r brca-pca-subtype, fig.width = 7, fig.height = 5, fig.align = "center", fig.cap = "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`."} plot_pca_audit(run$stages$pca_audit, color_by = "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 | `r round(sum(se$event) / 20, 1)` < 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: (Affymetrix GPL570, HG-U133 Plus 2.0; downloaded here via `GEOquery::getGEO()`). # Session info ```{r session-info} sessionInfo() ```