Application scenarios: atlas, HRD and treatment forensics
Source:vignettes/scenarios.Rmd
scenarios.Rmdmsuiter’s paper pillar 3 is the application-scenario layer:
end-to-end workflows that adjudicate against planted truth or an
independent oracle and persist their evidence tables. The three
executable drivers live under bench/scenarios/
(run_atlas.R, run_hrd.R,
run_treatment.R; committed evidence under
bench/scenarios/out/, each driver ends in its own OVERALL
PASS or OVERALL FAIL verdict line). This page runs one reduced, runnable
snippet of each against the same committed fixture generators the
drivers use, and cross-checks the live output against the committed
evidence table whenever that table is on disk. The browsable census of
the committed evidence is the scenario
gallery; this page’s own HRD
dual-oracle leg pins those tables byte for byte.
Atlas mode: the full cohort adjudication
The atlas scenario plants a synthetic cohort — five variant classes,
known channel plan, known arm differential — then drives the public
faces and adjudicates against the a-priori generator truth. The snippet
reproduces the ingestion-and-tally core; the committed driver adds the
class roster, extraction, hungarian fidelity and association legs with
their verdict lines. The cohort generator below is the very one the
driver and the governance tests share
(tests/testthat/helper-atlas.R), so the numbers are the
same bytes, not a lookalike.
if (is.na(PKG)) {
cat("rendered outside the package root: the fixture generators are not reachable here\n")
} else {
source(file.path(PKG, "tests", "testthat", "helper-tally.R"))
source(file.path(PKG, "tests", "testthat", "helper-atlas.R"))
fx <- .atlas_cohort(dir = file.path(tempdir(), "tour-atlas"))
truth <- fx$truth
mv <- ms_variants(fx$variants, genome = "synth", caller = "atlas-synth",
matched_normal = "none")
cat1 <- ms_tally(mv, genome = fx$twobit, mode = "SBS96")
cat("samples:", ncol(cat1@counts), " channels:", nrow(cat1@counts),
" planted:", paste(.ATLAS_PLANTED, collapse = ","), "\n")
cat("SBS96 plan exactness:",
identical(rowSums(cat1@counts)[truth$sbs_labels],
colSums(truth$sbs_plan)), "\n")
ev <- file.path(PKG, "bench", "scenarios", "out", "counts_sbs.csv")
if (file.exists(ev)) {
ref <- read.csv(ev, row.names = 1L, check.names = FALSE)
cat("live run vs committed evidence", ev, ": maxabs delta =",
max(abs(as.matrix(cat1@counts) - as.matrix(as.matrix(ref)))), "\n")
} else {
cat("committed evidence table not on disk (check tarball): the driver output lives in the repository under bench/scenarios/out/\n")
}
}
#> samples: 24 channels: 96 planted: SBS1,SBS5,SBS22a,SBS2,SBS13
#> SBS96 plan exactness: TRUE
#> live run vs committed evidence /home/runner/work/msuiter/msuiter/bench/scenarios/out/counts_sbs.csv : maxabs delta = 0The committed tables for this leg:
bench/scenarios/out/counts_sbs.csv,
counts_dbs.csv, qc_report.csv,
compare_sweep.csv, pbr.csv,
prevalence.csv, association.csv. The headline
committed table counts_sbs.csv is the 96-channel x
24-sample SBS grid — rows and columns re-derived from the witness bytes
themselves by the gallery census face, each cell adjudicated against
planted synthetic truth.
Byte-exact HRD: adjudication against the dual oracles
The HRD leg reproduces Cortes-Castro’s published values (PMID
31151873)
pLOH-by-segment-integration and the aneuploidy count, and every
adjudicated row must be byte-identical to the independent
oracle implementation carried in
tests/testthat/helper-hrd.R. The planted geometry encodes
the composite claim: only the double-positive sample (pLOH-high
and aneuploid) scores HRD.
if (is.na(PKG)) {
cat("rendered outside the package root: the oracle fixtures are not reachable here\n")
} else {
source(file.path(PKG, "tests", "testthat", "helper-hrd.R"))
chr_size <- .HRD_TEST_CHR_SIZE
arms <- do.call(rbind, lapply(c("chr1", "chr2", "chr13", "chr17", "chr22"),
function(ch) data.frame(chrom = c(ch, ch), location = c("p", "q"),
arm_start = c(1L, chr_size %/% 2L + 1L),
arm_end = c(chr_size %/% 2L, chr_size), stringsAsFactors = FALSE)))
whole <- function(sample, chrom, state, minor)
.make_hrd_segments(sample, chrom, start = 1L, end = chr_size,
state = state, minor_cn = minor)
segs <- do.call(rbind, c(
lapply(c("chr1", "chr2", "chr13", "chr17", "chr22"),
function(ch) whole("H1", ch, 5L, 0L)),
lapply(c("chr1", "chr2", "chr13", "chr17", "chr22"),
function(ch) whole("H3", ch, 3L, 1L)),
list(.make_hrd_segments("H2", "chr1", start = 1L, end = 450000L,
state = 2L, minor_cn = 0L),
whole("H4", "chr1", 2L, 1L))))
report <- ms_hrd_report(segs, genome = "GRCh37", caller = "ascat",
assay = "wgs", chr_size = chr_size,
arms = arms, ploidy = 2)
rows <- as.data.frame(report)
rows <- rows[order(rows$sample), , drop = FALSE]
print(rows[c("sample", "ploh_fraction", "aneuploidy_count",
"hrd_score", "evidence_tier")])
oracle_ok <- TRUE
for (i in seq_len(nrow(rows))) {
sub <- segs[segs$sample == rows$sample[i], , drop = FALSE]
oracle_ok <- oracle_ok &&
identical(rows$ploh_fraction[i], .hrd_ref_ploh(sub, chr_size)) &&
identical(rows$aneuploidy_count[i],
.hrd_ref_aneuploidy_count(sub, arms, 2, rm_black_arms = FALSE))
}
cat("dual-oracle byte-exactness:", oracle_ok, "\n")
stopifnot(isTRUE(oracle_ok))
stopifnot(identical(rows$hrd_score[rows$sample == "H1"], 1L) &&
all(rows$hrd_score[rows$sample != "H1"] == 0L))
ev <- file.path(PKG, "bench", "scenarios", "out", "hrd_report.csv")
if (file.exists(ev)) {
ref <- read.csv(ev, check.names = FALSE)
cat("live run vs committed evidence", ev, ":",
identical(as.integer(ref$hrd_score), rows$hrd_score) &&
identical(as.character(ref$sample), rows$sample), "\n")
}
}
#> sample ploh_fraction aneuploidy_count hrd_score evidence_tier
#> 1 H1 5.00 5 1 standard
#> 3 H2 0.45 0 0 standard
#> 2 H3 0.00 5 0 standard
#> 4 H4 0.00 0 0 standard
#> dual-oracle byte-exactness: TRUE
#> live run vs committed evidence /home/runner/work/msuiter/msuiter/bench/scenarios/out/hrd_report.csv : TRUEA missing arm table degrades to attestation, never fabrication: the
coverage_note ledger names the missing input and
aneuploidy_count stays at zero rather than being invented.
The committed table for this leg:
bench/scenarios/out/hrd_report.csv.
Treatment forensics: association plus clonal timing
No mainstream signature tool offers a treatment-history reporting module. msuiter’s face fits the governed therapy classes (temozolomide, platinum, thiopurine) against the officially recognised signature identities served from the sealed bundled COSMIC v3.6 reference, adjudicates each cell by a bootstrap contrast clearing an effect-size floor, and splits the exposed mean into clonal/subclonal timing shares. The snippet runs the complete face on the committed synthetic cohort; the degradation faces (no treatment records, no clonality) return honest empty tables plus an attestation ledger, never invented rows.
if (is.na(PKG)) {
cat("rendered outside the package root: the cohort generator is not reachable here\n")
} else {
source(file.path(PKG, "tests", "testthat", "helper-treatment.R"))
fx <- .ms_trt_cohort()
sigs <- sort(unique(unlist(.MS_TRT_MAP, use.names = FALSE)))
clon <- expand.grid(sample = fx$samples, signature = sigs)
clon <- clon[c("sample", "signature")]
clon$clonal_fraction <- 0.5
rep1 <- ms_treatment_report(fx$catalog, fx$treatments, clonality = clon,
n_boot = 200, seed = 7)
print(rep1)
tab <- as.data.frame(rep1)
cat("association rows:", nrow(tab), " adjudicated associated:",
sum(tab$is_associated), "\n")
prov <- attr(rep1, "provenance")
cat("provenance: refdb", prov$refdb_version, "build", prov$build,
"method", prov$method, "n_boot", prov$n_boot, "seed", prov$seed, "\n")
truth_ok <- all(tab$is_associated == vapply(seq_len(nrow(tab)), function(i)
fx$planted_truth(tab$treatment[i], tab$signature[i]), logical(1)))
cat("planted-fidelity adjudication:", truth_ok, "\n")
stopifnot(isTRUE(truth_ok))
ev <- file.path(PKG, "bench", "scenarios", "out", "treatment_association.csv")
if (file.exists(ev)) {
ref <- read.csv(ev, check.names = FALSE)
cat("live run vs committed evidence", ev, ":",
nrow(ref) == nrow(tab) &&
identical(as.logical(ref$is_associated), as.logical(tab$is_associated)),
"\n")
}
}
#> <msuiter treatment report>
#> refdb=v3.6 build=GRCh37 method=nnls n_boot=200 seed=7
#> association=full
#> association (12 rows):
#> treatment signature n_exposed n_unexposed prevalence_exposed
#> temozolomide SBS11 4 5 1.00
#> temozolomide SBS31 4 5 1.00
#> temozolomide SBS32 4 5 1.00
#> temozolomide SBS35 4 5 0.50
#> platinum SBS11 4 5 0.50
#> platinum SBS31 4 5 1.00
#> platinum SBS32 4 5 1.00
#> platinum SBS35 4 5 1.00
#> thiopurine SBS11 4 5 0.50
#> thiopurine SBS31 4 5 0.75
#> thiopurine SBS32 4 5 1.00
#> thiopurine SBS35 4 5 0.50
#> prevalence_unexposed mean_exposed mean_unexposed contrast contrast_lower
#> 1 559.035832 18.93254 540.103291 488.832976
#> 1 51.085840 16.21274 34.873097 5.404249
#> 1 47.251189 23.24596 24.005228 -9.867679
#> 1 5.691984 22.60875 -16.916770 -26.006420
#> 1 13.048698 18.93254 -5.883844 -21.063274
#> 1 581.466258 16.21274 565.253515 509.688787
#> 1 28.842356 23.24596 5.596395 -18.146152
#> 1 557.755725 22.60875 535.146970 489.410549
#> 1 20.870731 18.93254 1.938189 -12.035828
#> 1 29.504018 16.21274 13.291275 -13.678608
#> 1 615.612893 23.24596 592.366932 560.465510
#> 1 14.267304 22.60875 -8.341450 -21.037117
#> contrast_upper is_associated
#> 582.719269 TRUE
#> 74.480071 FALSE
#> 65.777400 FALSE
#> -0.290450 FALSE
#> 15.439007 FALSE
#> 619.726828 TRUE
#> 26.798330 FALSE
#> 584.637147 TRUE
#> 27.043682 FALSE
#> 31.981235 FALSE
#> 613.227812 TRUE
#> 8.672979 FALSE
#> timing (12 rows):
#> treatment signature n_exposed exposure_pre exposure_post pre_lower
#> temozolomide SBS11 4 279.517916 279.517916 253.971960
#> temozolomide SBS31 4 25.542920 25.542920 11.596076
#> temozolomide SBS32 4 23.625594 23.625594 7.136747
#> temozolomide SBS35 4 2.845992 2.845992 0.000000
#> platinum SBS11 4 6.524349 6.524349 0.000000
#> platinum SBS31 4 290.733129 290.733129 262.907846
#> platinum SBS32 4 14.421178 14.421178 3.499895
#> platinum SBS35 4 278.877862 278.877862 256.007940
#> thiopurine SBS11 4 10.435365 10.435365 4.135583
#> thiopurine SBS31 4 14.752009 14.752009 4.049051
#> thiopurine SBS32 4 307.806447 307.806447 291.316141
#> thiopurine SBS35 4 7.133652 7.133652 1.719598
#> pre_upper post_lower post_upper
#> 299.13046 253.971960 299.13046
#> 46.16511 11.596076 46.16511
#> 42.53066 7.136747 42.53066
#> 10.17515 0.000000 10.17515
#> 15.87979 0.000000 15.87979
#> 316.87740 262.907846 316.87740
#> 23.15746 3.499895 23.15746
#> 303.90452 256.007940 303.90452
#> 23.94582 4.135583 23.94582
#> 23.97608 4.049051 23.97608
#> 317.81182 291.316141 317.81182
#> 14.82868 1.719598 14.82868
#> association rows: 12 adjudicated associated: 4
#> provenance: refdb v3.6 build GRCh37 method nnls n_boot 200 seed 7
#> planted-fidelity adjudication: TRUE
#> live run vs committed evidence /home/runner/work/msuiter/msuiter/bench/scenarios/out/treatment_association.csv : TRUEThe committed tables for this leg:
bench/scenarios/out/treatment_association.csv,
treatment_timing.csv. The recognition anchors and the
adjudication-floor rationale are documented at the face’s manual page
(?ms_treatment_report) and in the Phase-8 design memo under
vignettes/articles/devlog-*.qmd.