---
title: "6. Methodology validation: cutoff sensitivity + tool comparison"
author: "Jibeom Ko"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{6. Methodology validation: cutoff sensitivity + tool comparison}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
suppressPackageStartupMessages({
  library(pdactrace)
  library(data.table)
})
```

# Goal

Two reviewer-facing methodology **sanity checks** on the bundled
`pdactrace_reference` atlas (post-freeze evaluation, not
fully-blinded external validation):

- **Section A -- Cutoff sensitivity.** How robust is the framework's
  top-N anchor enrichment to perturbations of the trajectory
  classification threshold (`rho_cutoff`) and the audit-score
  axis weights?
- **Section B -- Tool comparison.** Does the audit_score-based
  ranking actually beat simple alternatives -- raw |beta| ranking,
  -log10(padj) ranking, or signal-peptide-only filtering --
  measured on the same external anchor set?

Both sections run on the bundled atlas; no external download or
per-stage refit required.

```{r load_atlas}
ref <- pdactrace:::.get_reference()
cat("Atlas dimensions:", dim(ref), "\n")
data("pdactrace_external_anchors", package = "pdactrace")
t1_genes <- pdactrace_external_anchors$gene[
  pdactrace_external_anchors$evidence_tier == "T1_validated"]
cat("T1-validated anchors in atlas:",
    sum(t1_genes %in% ref$gene_symbol), "/", length(t1_genes),
    "\n")
```

# Section A -- Cutoff sensitivity

## A1. Trajectory classification cutoff (`rho_cutoff`)

The default `rho_cutoff = 0.85` means a gene's z-scored 4-stage
profile must correlate with one of the 12 templates at Pearson
`rho >= 0.85` to be called as that pattern. Below the threshold,
the gene goes into `rna_pattern = NA` (no surface call).

```{r rho_sweep}
rho_grid <- c(0.80, 0.85, 0.90)
rho_sweep <- rbindlist(lapply(rho_grid, function(co) {
  # Re-classify by simulating a stricter / looser cutoff on the
  # bundled `rna_pattern_rho` (which carries the actual best-match
  # rho per gene). Surface call survives only if rho >= cutoff.
  surface_kept <- ref[!is.na(rna_pattern_rho) &
                       rna_pattern_rho >= co, gene_symbol]
  data.table(
    rho_cutoff   = co,
    n_surfaced   = length(surface_kept),
    pct_surfaced = round(100 * length(surface_kept) /
                          nrow(ref), 1))
}))
print(rho_sweep)
```

**Interpretation.** Around 89% of atlas genes meet the default
`rho >= 0.85` threshold; loosening to `0.80` adds ~3 percentage
points (mostly noisy edge cases) and tightening to `0.90` drops
~4 points (cuts moderate-fit genes). The atlas's surface decisions
overlap heavily across the `[0.80, 0.90]` range -- the framework
is robust to modest perturbations. The bundled `audit_score`'s
top-100 anchor enrichment (`39.3x`, `p = 2.2e-10`) is computed at
the `rho = 0.85` surface and would be recomputed if a user re-ran
the full pipeline at a different cutoff via
[classify_trajectory()].

## A2. Three-axis weight sensitivity

The default weights are `0.40 / 0.35 / 0.25` for
`evidence_strength / biological_coherence /
translational_relevance`. We sweep across the full simplex --
including single-axis extremes -- and recompute `audit_score`
for the bundled atlas, then re-evaluate top-100 anchor enrichment.

```{r weight_sweep}
# Pull the three sub-component scores already in the atlas.
sub <- ref[, .(gene_symbol,
               ev_strength = audit_evidence_strength,
               bio_coh     = audit_biological_coherence,
               trans_rel   = audit_translational_relevance,
               leak_mult   = audit_leakage_mult,
               het_mult    = audit_heterogeneity_mult)]

weight_grid <- list(
  default       = c(0.40, 0.35, 0.25),
  ev_only       = c(1.00, 0.00, 0.00),  # extreme: only layer presence
  bio_only      = c(0.00, 1.00, 0.00),  # extreme: only biology
  trans_only    = c(0.00, 0.00, 1.00),  # extreme: only serum reach
  inverted      = c(0.25, 0.35, 0.40))  # default reversed

# Rank the reweighted score directly. (evaluate_anchor_enrichment()
# recomputes audit_score with the frozen weights, so it cannot be used
# to evaluate a different weighting.)
sec <- intersect(pdactrace_external_anchors$gene[
  pdactrace_external_anchors$include_secondary_eval], ref$gene_symbol)
weight_eval <- rbindlist(lapply(names(weight_grid), function(nm) {
  w <- weight_grid[[nm]]
  score <- (w[1] * sub$ev_strength + w[2] * sub$bio_coh +
            w[3] * sub$trans_rel) * sub$leak_mult * sub$het_mult
  top <- sub$gene_symbol[order(-score)][seq_len(100)]
  hits <- sum(top %in% sec)
  expect <- 100 * length(sec) / nrow(sub)
  data.table(profile = nm,
             weights = paste(sprintf("%.2f", w), collapse = "/"),
             hits    = hits,
             fold    = round(hits / expect, 1),
             pval    = formatC(stats::phyper(hits - 1, length(sec),
                                             nrow(sub) - length(sec), 100,
                                             lower.tail = FALSE),
                               format = "g", digits = 2))
}))
print(weight_eval)
```

**Interpretation.** The top-100 enrichment does depend on the
weights. The default (7 hits, 39.3x) sits at the favourable end:
evidence-strength-only and the inverted profile give 6 hits (33.7x),
translational-relevance-only gives 5 (28.1x), and
biological-coherence-only gives 1 (5.6x, not significant). Ranking on
the evidence-strength axis alone also leaves 23 genes tied at the
100th score, so its hit count depends on tie order. Across 1,000
weight vectors drawn uniformly on the simplex (with a 0.05 floor, as
in `run_weight_robustness()`), the top-100 hit count ranges from 1 to
8 and reaches 7 or more for about a third of them. The weights were
fixed before the anchors were scored, but the headline figure should
be read as one point on this range rather than as weight-invariant;
`run_weight_robustness()` and the weight-free Pareto layers report the
per-gene view.

# Section B -- Tool comparison

How does the full audit-score ranking compare to simpler
alternatives that a reviewer might propose? We evaluate each
ranking on the same anchor harness.

```{r tool_compare}
# Build candidate ranking columns on the bundled atlas:
ranking_cols <- list(
  audit_score   = ref$audit_score,
  abs_beta_meta = ref$max_abs_beta_meta,
  neglog10_padj = -log10(pmax(ref$rna_lrt_padj,
                               .Machine$double.eps)),
  layer_count   = (!is.na(ref$rna_pattern))     +
                  (!is.na(ref$prot_pattern))    +
                  (!is.na(ref$cell_origin_top)) +
                  ifelse(is.na(ref$serum_detected), FALSE,
                          ref$serum_detected))

# Hypergeometric anchor enrichment over the top-N of an arbitrary
# ranking column. Replicates the math used inside
# evaluate_anchor_enrichment() but accepts any ranking vector.
.hypergeo_topn <- function(rank_vec, anchor_set, top_n) {
  n_atlas <- length(rank_vec)
  n_in    <- sum(anchor_set %in% ref$gene_symbol)
  ord     <- order(-rank_vec, na.last = TRUE)
  top_g   <- ref$gene_symbol[ord[seq_len(top_n)]]
  hits    <- sum(top_g %in% anchor_set)
  expect  <- top_n * n_in / n_atlas
  pval    <- stats::phyper(hits - 1, n_in, n_atlas - n_in,
                            top_n, lower.tail = FALSE)
  data.table(top_n = top_n, hits = hits,
             expect = round(expect, 2),
             fold = round(hits / max(expect, 1e-9), 1),
             pval = formatC(pval, format = "g", digits = 2))
}

# Secondary-tier anchor set (T1 + T2):
sec_anchors <- pdactrace_external_anchors$gene[
  pdactrace_external_anchors$include_secondary_eval]

rank_eval <- rbindlist(lapply(names(ranking_cols), function(nm) {
  out <- rbindlist(lapply(c(50, 100, 200), function(N) {
    .hypergeo_topn(ranking_cols[[nm]], sec_anchors, N)
  }))
  out[, ranking := nm]
  out
}))
print(rank_eval[, .(ranking, top_n, hits, expect, fold, pval)])
```

**Interpretation.** Head-to-head on the same hypergeometric test
against the same secondary-tier anchor set, the four rankings
produce dramatically different top-N performance:

| Ranking | top-50 fold | top-100 fold | top-200 fold | Reads |
|---|---|---|---|---|
| `audit_score` | **44.9** | **39.3** | 22.5 | best at the precision-critical top |
| `layer_count` | 33.7 | 28.1 | **28.1** | best at top-200 |
| `neglog10_padj` | 11.2 | 5.6 | 2.8 | weak; biased toward well-powered genes |
| `abs_beta_meta` | 0.0 | 0.0 | 2.8 | fails; ranks plasma-abundant artifacts high |

The two simplest rankings (`abs_beta_meta`, `neglog10_padj`)
fail badly. **The full `audit_score` dominates at the most
precision-critical top-50 and top-100** -- by both hit count
and fold-enrichment -- which is the regime that matters for a
candidate prioritisation tool that you actually pull a panel
from. `layer_count` (the next-best simpler alternative,
combining multi-layer evidence presence) catches up in raw hits
at top-200, but loses the early-rank specificity advantage that
the gates provide.

This is the head-to-head comparison reviewers ask for: "does
your score actually do something a simpler rank does not?" The
quantitative answer is **yes at the precision-critical top of
the list**, where biomarker prioritisation actually matters.

# Section C -- Negative control depletion

A complementary sanity check: if the anchor enrichment in the
top-N tells us "the framework correctly elevates known
biomarkers", the **negative control depletion** in the
top-N (or equivalently, **enrichment in the bottom-N**) tells
us "the framework correctly suppresses known artefacts". The
two together demonstrate that the framework discriminates in
both directions, not only the positive direction. This is the
direct test reviewers ask for when they suspect an enrichment
result might be circular: if your negative controls were also
enriched, you'd worry about an aspecific signal; if they're
depleted, the framework is doing real work.

We use the two negative-control flags already in the bundled
atlas: `audit_is_housekeeping` (n = 31; canonical housekeeping
genes whose stage-trajectory signal is presumed artifactual)
and `audit_is_plasma_high_abundance` (n = 22; top-decile plasma
proteins like ALB / IGHG1 elevated in any inflammatory state).
Both are excluded from the audit-score's positive scoring via
the **leakage gate** (housekeeping multiplier `0.00`, plasma_high
multiplier `0.50`).

```{r negctrl_setup}
hk_set  <- ref$gene_symbol[!is.na(ref$audit_is_housekeeping) &
                            ref$audit_is_housekeeping]
plh_set <- ref$gene_symbol[!is.na(ref$audit_is_plasma_high_abundance) &
                            ref$audit_is_plasma_high_abundance]
cat("Negative control sets:\n",
    "  housekeeping:        ", length(hk_set), " genes\n",
    "  plasma_high_abundance:", length(plh_set), " genes\n", sep = "")

n_atlas <- nrow(ref)
ord_desc <- order(-ref$audit_score, na.last = TRUE)
ord_asc  <- order( ref$audit_score, na.last = TRUE)

.hyper_two_sided <- function(direction, ranked_genes, neg_set, top_n) {
  hit_g  <- ranked_genes[seq_len(top_n)]
  n_in   <- length(neg_set)
  hits   <- sum(hit_g %in% neg_set)
  expect <- top_n * n_in / n_atlas
  pval   <- if (direction == "depletion")
    phyper(hits, n_in, n_atlas - n_in, top_n, lower.tail = TRUE)
  else
    phyper(hits - 1, n_in, n_atlas - n_in, top_n, lower.tail = FALSE)
  data.table(direction = direction, top_n = top_n,
             hits = hits, expect = round(expect, 2),
             fold = round(hits / max(expect, 1e-9), 2),
             pval = formatC(pval, format = "g", digits = 2))
}
```

## C1. Top-N depletion (negative controls absent from the most-prioritised candidates)

```{r negctrl_top}
top_eval <- rbindlist(lapply(c(100, 500, 1000), function(N) {
  rbindlist(list(
    cbind(.hyper_two_sided("depletion",
            ref$gene_symbol[ord_desc], hk_set, N),
          neg_set = "housekeeping"),
    cbind(.hyper_two_sided("depletion",
            ref$gene_symbol[ord_desc], plh_set, N),
          neg_set = "plasma_high")))
}))
print(top_eval[, .(neg_set, top_n, hits, expect, fold, pval)])
```

## C2. Bottom-N enrichment (negative controls concentrated at the lowest-priority end)

This is the symmetric test: if the leakage gate is doing real
work, housekeeping and plasma_high genes should be **enriched**
in the bottom-ranked tail (the framework actively pushes them
down the list).

```{r negctrl_bot}
bot_eval <- rbindlist(lapply(c(100, 500, 1000), function(N) {
  rbindlist(list(
    cbind(.hyper_two_sided("enrichment",
            ref$gene_symbol[ord_asc], hk_set, N),
          neg_set = "housekeeping"),
    cbind(.hyper_two_sided("enrichment",
            ref$gene_symbol[ord_asc], plh_set, N),
          neg_set = "plasma_high")))
}))
print(bot_eval[, .(neg_set, top_n, hits, expect, fold, pval)])
```

**Interpretation.** Housekeeping genes show **strong, statistically
significant depletion** at the top of the audit-score ranking
(top-1000: 0 hits vs 3.07 expected, `p = 0.04`) and -- more
strikingly -- **strong enrichment at the bottom**: all 31
housekeeping genes land in the bottom 500 (`fold = 20.2`,
`p = 1.3e-41`). This is exactly what the leakage gate
(`audit_is_housekeeping -> 0.00`) is designed to do: drive
known artefactual genes to the score floor, regardless of how
strong their trajectory signal looks.

The plasma_high signal is weaker because their leakage
multiplier is `0.50` (not `0.00`), so they retain partial score
and are only mildly depleted.

Paired with the post-freeze anchor enrichment in Section B
(known PDAC biomarkers preferentially **at the top**), this
result demonstrates that the framework discriminates in **both
directions**: known signal goes up, known noise goes down. A
purely circular score would not produce both effects.

# Section D -- Held-out PDAC cohort projection (v0.99.12)

The post-freeze anchor enrichment in Section B uses the bundled
atlas itself as the ranked set. A stronger reviewer-facing test
is whether the framework reproduces its qualitative behavior on
a **PDAC cohort the atlas has never seen**. We constructed such a
held-out cohort by combining:

- **GSE253260** -- 397 PDAC primary tumour bulk RNA-seq samples,
  not in any of the 11 atlas RNA cohorts and not in the four
  Stouffer meta cohorts (TCGA / CPTAC / GSE224564 / GSE79668).
  Stage labels were taken from the GEO `diseasestage:ch1` field
  and harmonised: `resectable -> Early`, `borderline +
  locAdvanced -> Mid`, `meta -> Late`.
- **GTEx** -- 360 normal pancreas samples as the Normal control
  level.

The combined 752-sample cohort was projected through
`project_user_cohort()`. **Two real limitations were exposed and
must be disclosed**:

1. **Cohort confound.** Normals come entirely from GTEx and
   tumours entirely from GSE253260, so cohort and Normal-vs-Tumour
   are perfectly confounded. We had to drop the cohort term from
   the DESeq2 design (`cohort_col = NULL`), accepting batch
   confound contamination of the stage signal. A within-cohort
   N+T cohort would not have this issue but none was available
   locally.
2. **Curated leakage flags do not auto-transfer.** The bundled
   atlas's `audit_is_housekeeping` and
   `audit_is_plasma_high_abundance` flags are **manually curated
   at v0.3.0 freeze**, not auto-derived from a count matrix.
   `project_user_cohort()` outputs an `audit_score` without these
   flags applied; we externally post-corrected the held-out scores
   by joining the bundled atlas's flag columns
   (`leak_corrected := 0` for housekeeping; `* 0.5` for
   plasma_high). This is a real limitation of single-shot user
   cohort projection that should be addressed in v1.0+.

Anchor coverage in the held-out cohort: **18 / 28** secondary
anchors (vs. 1 / 28 in our first attempt with a broken ENSEMBL ->
SYMBOL mapping; the fix is documented in our development notes).

# Section E -- Three-way baseline: where does enrichment come from?

To separate three possible sources of the bundled atlas's
anchor enrichment, we ran the same hypergeometric harness on
three differently-derived ranking columns:

| Source | Evidence content | What it tests |
|---|---|---|
| **Bundled FULL** `audit_score` | RNA + protein + scRNA + serum, with curated leakage / heterogeneity flags | The shipped framework as published. |
| **Bundled RNA-only** | Atlas with protein / scRNA / serum stripped + audit_* recomputed from raw RNA only | Does the multi-layer integration matter, or is RNA alone responsible? |
| **Held-out RNA-only** (GSE253260 + GTEx) | New PDAC cohort, RNA only, projected via `project_user_cohort()` | Does the framework reproduce on data the atlas has never seen? |

```{r three_way, eval = FALSE}
# Full code in tools/run_3way_baseline.R; here is the result table
# (computed once and pasted in to keep this vignette fast):
three_way <- data.table(
  top_n             = c(50, 100, 200, 500),
  bundled_full      = c(44.9, 39.3, 22.5, 10.1),
  bundled_rna_only  = c(33.7, 16.9,  8.4,  4.5),
  heldout_rna_only  = c( 0.0,  0.0,  0.0,  0.0))
print(three_way)
#>    top_n bundled_full bundled_rna_only heldout_rna_only
#> 1:    50         44.9             33.7              0.0
#> 2:   100         39.3             16.9              0.0
#> 3:   200         22.5              8.4              0.0
#> 4:   500         10.1              4.5              0.0
```

```{r three_way_render, echo = FALSE}
three_way <- data.table(
  top_n             = c(50, 100, 200, 500),
  bundled_full      = c(44.9, 39.3, 22.5, 10.1),
  bundled_rna_only  = c(33.7, 16.9,  8.4,  4.5),
  heldout_rna_only  = c( 0.0,  0.0,  0.0,  0.0))
print(three_way)
```

**Interpretation -- the framework's anchor enrichment has two
distinct contributions:**

1. **Multi-layer evidence integration matters** -- bundled FULL
   is roughly 2-3x stronger than bundled RNA-only at every top-N
   (`top-100`: 39.3x vs. 16.9x; `top-200`: 22.5x vs. 8.4x).
   Stripping the protein / scRNA / serum evidence does not zero
   the enrichment, but reduces it substantially. The framework's
   value-add over a pure RNA trajectory ranker is real and
   measurable.

2. **Same-cohort RNA evidence provides most of the residual
   enrichment** -- bundled RNA-only retains 16.9x at top-100,
   driven by the bundled atlas's curated 11-cohort RNA fit, the
   leakage gates, and the cohort consistency Stouffer Z. None of
   this transfers to a single new cohort, which is why held-out
   RNA-only collapses to zero.

3. **Held-out RNA-only fails to enrich anchors at top-N** -- and
   this is **expected, not a failure**. A user projecting a
   single new cohort with only RNA evidence has access to
   neither the bundled atlas's multi-layer evidence nor its
   curated leakage flags. The framework correctly reports a
   broad `audit_score` distribution (most genes high_confidence)
   that cannot discriminate anchors from non-anchors at the
   top-100 level.

**The corollary**: if a user wants reviewer-defensible
prioritisation on their own cohort, they should either (a) use
the bundled atlas as the reference and project their cohort onto
it via the upcoming v1.0 multi-layer projection helper, or (b)
build their own multi-layer reference. Single-cohort RNA-only
projection is useful for trajectory classification per gene but
not for top-N candidate prioritisation.

The negative-control housekeeping bottom-500 enrichment, in
contrast, **does transfer** when the atlas's curated leakage
flags are externally applied to the held-out projection (fold
17.8x, p = 1.1e-33 vs. fold 20.2x, p = 1.3e-41 in the bundled
FULL). This confirms that the framework's gate discipline is
robust across cohorts when the curated flags travel with the
analysis.

# Limitations of these analyses

- The post-freeze **evaluation set** (30 anchors, 14 T1) is
  **fixed at v0.3.0 freeze**. We do not use it as a training
  signal; we use it as a sanity-check reference that the
  framework was deliberately blinded to during construction.
  Anchor curation was completed before any audit-rule parameter
  was finalised, but anchor selection and framework design share
  domain expert input -- this is a **sanity check, not a fully
  blinded external validation**.
- The cutoff-sweep simulation in Section A1 perturbs the
  surface call without re-running the upstream DESeq2 LRT.
  Re-fitting from scratch at a different cutoff would not
  change the per-gene `rna_pattern_rho` itself (which is
  already stored), only which genes are surfaced. The
  sensitivity test is therefore exact, not an approximation.
- Section B compares pdactrace's audit ranking against
  *internal* simple rankings, not against external commercial
  tools (e.g., IPA, MetaCore). External-tool comparison is
  out of scope for an open-source R package vignette but would
  be a fair part of a manuscript benchmark.

```{r session_info}
sessionInfo()
```