---
title: "1. Introduction to rnaSentry"
author: "Abdelkarim Hani Ghrieb"
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
vignette: >
  %\VignetteIndexEntry{1. Introduction to rnaSentry}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r rnaSentry-setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

# 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

```{r install, eval = FALSE}
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.

```{r data-intake}
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)
S4Vectors::metadata(se_in)$flags
```

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.

```{r qc-explore}
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
```

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

```{r sex-check}
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)
```

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

```{r run-pipeline}
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
```

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.

```{r unlock-signature}
# 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:

```{r stage-flags}
run$stages$build_signature$flags
```

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

```{r ext-validate}
val <- validate_external(run$stages$build_signature, external,
                         cutpoint = run$stages$km_curve$cutpoint)
val
```

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:

```{r intro-ext-km, fig.width = 7, fig.height = 5, fig.align = "center", fig.cap = "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."}
plot(val)
```

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:

```{r wrong-scale, warning = TRUE}
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)
```

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:

```{r underpowered, warning = TRUE}
small <- se_all[, 1:40]
build_signature(small, "time", "event", top_n = 10, repeats = 1,
                folds = 2, seed = 7)
```

Both warnings are also recorded in each stage's `flags` data.frame, so they
propagate into the rendered report's audit trail.

# Session info

```{r session-info}
sessionInfo()
```