Benchmark Simulation#

Transcriptome-scale hierarchical gamma-Poisson simulator for generating realistic single-cell clinical-trial data with known ground truth, plus the calibration, validation and orchestration tools used to benchmark statistical methods against it.

The generative model, per cell c of participant i at visit t, gene g:

\[\log \mu_{icgt} = \log L_{ic} + \alpha_g + b_{ig} + u_{igt} + \gamma_g \mathrm{Post}_t + \beta_g (T_i \times \mathrm{Post}_t)\]

with \(Y \sim \mathrm{NB}(\mu, \phi_g)\) and \(\mathrm{Var} = \mu + \phi_g \mu^2\). The full transcriptome is simulated; analysis panels are drawn from it as nested subsets of the detectable genes, so normalisation and library-size offsets see the whole transcriptome exactly as a real workflow would.

Simulator#

class sctrial.benchmark.TranscriptomeSimConfig(n_per_arm: int = 6, design: ~typing.Literal['two_arm', 'single_arm'] = 'two_arm', n_genes_transcriptome: int = 20284, arm_ratio: tuple[int, int] | None = None, panel_min_mean_count: float = 0.05, panel_sizes: tuple[int, ...] = (50, 200, 500, 2000), cells_per_pv_mean: float = 5898.0, cells_per_pv_cv: float = 0.911, cells_per_pv_min: int = 147, cells_per_pv_max: int = 27653, cells_scale: float = 1.0, cells_per_pv_fixed: int | None = None, missing_rate: float = 0.0, use_empirical_library: bool = True, use_empirical_cells_per_pv: bool = True, empirical_path: str | None = None, lib_log_mean: float = 7.7333, lib_log_sd: float = 0.8, use_empirical_gene_rates: bool = True, gene_rate_log_mean: float = -4.5338, gene_rate_log_sd: float = 2.6428, use_empirical_dispersion: bool = True, dispersion_median: float = 0.4416, dispersion_mean_slope: float = -0.1753, dispersion_anchor: float = -10.5486, dispersion_residual_sd: float = 1.3224, between_participant_sd: float = 0.5969, prepost_corr: float = 0.3987, time_effect: float = 0.0, effects: dict[str, float] = <factory>, seed: int = 0)[source]#

Bases: object

Configuration for the transcriptome-scale simulator.

Defaults are the TNBC-calibrated values. Every one is an empirical quantity, so a run with defaults is a calibrated run — the previous design’s defaults were arbitrary and silently used whenever calibration was not threaded through.

property participant_sd: float#

the participant component of the total between-participant SD.

Type:: sigma_b

property participant_visit_sd: float#

the participant x visit component.

corr(pre, post) = sigma_b^2 / (sigma_b^2 + sigma_u^2), so sigma_u = sd * sqrt(1 - corr). This is the level whose absence made participant_sd inert in the old simulator.

Type:: sigma_u

sctrial.benchmark.simulate_trial_v2(cfg: TranscriptomeSimConfig) → dict[source]#

Simulate a full transcriptome under the three-level hierarchical NB model.

Returns:

A dictionary with the following keys:

adata: Cell-level raw counts (sparse); obs has participant/visit/arm and the OBSERVED library total.
pseudobulk_counts: Participant x visit summed counts, full transcriptome.
pseudobulk_means: Participant x visit mean counts, full transcriptome.
gene_names: Transcriptome gene names.
panels: Nested panel -> gene names.
truth: Beta_g by gene (the injected effect).
oracle: Per-gene truth on each method class’s own estimand scale (see oracle_estimands()).
latent: VALIDATION ONLY: b_ig, u_igt, alpha_g, phi_g. Never use as an analysis input.
config: The config used.

Return type:

dict

sctrial.benchmark.make_signal(panel_genes: list[str], n_signal: int, architecture: Literal['balanced', 'heterogeneous', 'one_directional'] = 'balanced', magnitude: float = 0.5, rng: Generator | None = None) → dict[str, float][source]#

Effect sizes for the tested panel under one of three signal architectures.

balanced (primary): Half +magnitude, half -magnitude. No net directional shift, so the library size is unmoved and there is no compositional artifact.
heterogeneous (primary): Symmetric mixture of weak/moderate/large effects — the most realistic.
one_directional (stress test): All +magnitude. This is the previous design. Retain it, but label it a COMPOSITION-STRESS scenario: the coordinated shift moves the library-size reference, and roughly two thirds of the dreamlet inflation previously attributed to empirical-Bayes moderation is attributable to it.

sctrial.benchmark.nested_panels(cfg: TranscriptomeSimConfig, sizes: tuple[int, ...] | None = None, rng: Generator | None = None) → dict[int, list[int]][source]#

Nested analysis panels drawn from the DETECTABLE genes.

Nesting is what lets a panel-size effect be separated from a gene-identity effect. With independently drawn panels the two are confounded, and the previous benchmark’s “progressive miscalibration with panel size” could not be distinguished from a change in which genes were tested.

Panels are drawn from eligible_panel_genes(), not from the whole transcriptome. The transcriptome is still simulated in full and still supplies every normalisation denominator and offset; it is only the TESTED set that is restricted, exactly as a real pipeline restricts it.

sctrial.benchmark.oracle_estimands(params: dict) → dict[str, ndarray][source]#

Per-gene POPULATION truth on each METHOD CLASS’s own estimand scale.

Different methods do not estimate the same functional, so scoring them all against the injected beta silently penalises whichever method’s estimand differs most from it. Two published conclusions have already been produced that way (log2-versus-natural-log, then this).

count_link: The log-link coefficient: exactly beta_g. This is the estimand of NEBULA (NB log link) and of log-CPM pseudobulk models (dreamlet, limma-voom, edgeR): with full-transcriptome normalisation the library reference does not move with a panel-restricted effect, so DiD[log CPM_g] = beta_g.
log1p_cpm: The estimand targeted by sctrial and the Wilcoxon change score: the difference-in-differences of log(1 + CPM) at participant level. Because d/dx log(1+x) = 1/(1+x) this equals beta_g only when CPM >> 1, and is attenuated for low-expression genes at realistic depth. The attenuation is a real property of the estimand and is reported, not corrected away.

MARGINALISED, NOT CONDITIONED. The expectation is taken over the participant and participant x visit random effects rather than evaluated at their realised values. Conditioning on the realised draw makes the “truth” a random quantity: under a true null it returns values like +0.15 and -0.86 instead of zero, so a perfectly calibrated method would be scored as biased and a null scenario would have a non-null target. The random effects are exactly the variability the standard error is meant to cover.

The expectation E[log(1 + exp(m + eps))], eps ~ N(0, sigma_b^2 + sigma_u^2), is evaluated by Gauss-Hermite quadrature.

sctrial.benchmark.simulator_v2.load_tnbc_targets(path: str | Path | None = None) → dict[source]#

Empirical TNBC targets measured from the processed h5ad.

These are the quantities a defensible calibration must reproduce. Measured by scripts/calibrate_simulator.py targets from the v5 TNBC object (141,553 cells x 20,284 genes, 12 paired participants, 6 v 6 arms).

sctrial.benchmark.simulator_v2.load_empirical(path: str | Path | None = None) → dict | None[source]#

Empirical TNBC nuisance pools (library sizes, cells per participant-visit).

Returns None if absent, so the simulator falls back to the parametric fits rather than failing — but the fits are known to reproduce the library-size distribution poorly, so the empirical file should normally be present.

sctrial.benchmark.simulator_v2.build_params(cfg: TranscriptomeSimConfig) → dict[source]#

Draw every gene-level and participant-level latent parameter.

Separated from cell generation so that the full simulation (simulate_trial_v2()) and the calibration gates (sctrial.benchmark.gates, which only need summary statistics and must not materialise 141k x 20k count matrices) consume one implementation of the generative model. Two implementations of a generative model is how a calibration silently stops describing the thing it calibrates.

sctrial.benchmark.simulator_v2.gene_baseline_rates(cfg: TranscriptomeSimConfig) → ndarray[source]#

alpha_g: log gene proportions, summing to 1 on the exponential scale.

Reproduced exactly as build_params() draws it – it is the FIRST draw from default_rng(cfg.seed), so a fresh generator reproduces it bit for bit. That lets the analysis panel be chosen before the effects are defined (which the orchestrator needs, since the signal is injected on the tested genes) without simulating anything. tests/test_benchmark.py asserts the two agree.

A gene with zero empirical counts gets alpha = -inf, i.e. it is never expressed. That is faithful rather than a degenerate case, and callers must treat exp(alpha) as the rate rather than assuming alpha is finite.

sctrial.benchmark.simulator_v2.expected_counts_per_cell(cfg: TranscriptomeSimConfig) → ndarray[source]#: Expected counts per cell for each gene: E[L] * exp(alpha_g).

sctrial.benchmark.simulator_v2.eligible_panel_genes(cfg: TranscriptomeSimConfig) → ndarray[source]#: Indices of genes a real analysis would carry into testing.

sctrial.benchmark.simulator_v2.iter_pv_blocks(cfg: TranscriptomeSimConfig, params: dict | None = None)[source]#

Yield one participant-visit block of cell-level counts at a time.

Yielding rather than returning is what makes full-TNBC-scale Monte Carlo calibration tractable: a single replicate is 141k cells x 20,284 genes, which never has to exist at once if the consumer only needs summary statistics.

Yields:

dict with participant, visit, arm, counts (n_cells x G int32),
library (the drawn latent L per cell), pi, ti.

Method contracts#

Different methods do not estimate the same functional, and they require different input representations. Both are declared explicitly rather than inferred at call time.

sctrial.benchmark.prepare_inputs(sim: dict, panel_genes: list[str]) → dict[source]#

Build every method’s contracted input from one simulated dataset.

Returns:: participant_log1p_cpm outcome frame for sctrial / Wilcoxon pseudobulk_counts panel raw counts for the count-based models lib_size full-transcriptome total per pseudobulk sample cell_counts cell-level AnnData restricted to the panel cell_lib_size full-transcriptome library size per cell oracle per-gene truth keyed by estimand name panel_genes the tested panel
Return type:: dict with

sctrial.benchmark.contracts.prepare_inputs_from_adata(adata, panel_genes: list[str], participant_col: str = 'participant', visit_col: str = 'visit', arm_col: str = 'arm', counts_layer: str = 'counts') → dict[source]#

The same contracts, applied to REAL data.

The permutation and subsampling analyses run the identical methods on real cohorts and must therefore hand them the identical representations. They previously built their own pseudobulk and normalised inside the tested panel, so the real-data results were produced under a different normalisation scope from the simulation used to characterise those same methods. Sharing this function is what stops the two drifting apart again.

The CPM denominator is the whole measured transcriptome, not the tested panel, exactly as in prepare_inputs().

sctrial.benchmark.contracts.participant_log1p_cpm(pseudobulk_counts: DataFrame, panel_genes: list[str], gene_cols: list[str] | None = None) → DataFrame[source]#

log(1 + CPM) per participant-visit, normalised on the TRANSCRIPTOME.

The denominator is the sum over gene_cols (the full transcriptome), not over panel_genes. Normalising within the tested panel makes the reference move with the signal: under a coordinated effect the panel total shifts, and every null gene in the panel acquires an offsetting apparent effect. That artifact was previously attributed to empirical-Bayes moderation.

Panel selection happens AFTER normalisation, which is what a real workflow does and what makes panel size separable from normalisation scope.

sctrial.benchmark.METHOD_INPUT: dict[str, str]#

The input representation each method expects.

participant_log1p_cpm — sctrial_did, sctrial_mixed, wilcoxon_paired
pseudobulk_counts — dreamlet, limma_voom
cell_counts — nebula

sctrial.benchmark.METHOD_ESTIMAND: dict[str, str]#

The estimand scale each method is scored against.

log1p_cpm — sctrial_did, sctrial_mixed, wilcoxon_paired
count_link — dreamlet, limma_voom, nebula

Calibration#

Estimators that measure empirical properties of a real dataset (library-size distribution, per-gene dispersion, per-gene baseline rates) used to calibrate the simulator so it reproduces observed nuisance statistics.

class sctrial.benchmark.calibration.SummaryAccumulator(n_genes: int, gene_names: list[str] | None = None, min_mean_count: float = 0.02, min_strata_detected: int = 5, min_cpm_genewise: float = 0.0, cell_umi: list = <factory>, cell_genes: list = <factory>, pv_rows: list = <factory>, n_cells_total: int = 0)[source]#

Bases: object

Accumulates every validation statistic block by block.

A block is one homogeneous stratum’s cell x gene count matrix. Nothing here requires the full matrix to exist, which is what makes full-scale Monte Carlo calibration (141k cells x 20k genes x hundreds of replicates) affordable.

Conditional dispersion is accumulated as pooled within-stratum moments: within stratum k the fitted mean is mu_c = s_c * lambda_k with lambda_k = sum(y) / sum(s), so

sum (y - mu)^2 = sum y^2 - 2 lambda sum(y s) + lambda^2 sum(s^2)

and every term is a per-gene reduction over cells. The NB2 moment estimate is then alpha = (sum (y-mu)^2 - sum mu) / sum mu^2 pooled over strata.

add_block(counts: ndarray, participant: str, visit: str, arm: str = 'NA', stratum: str | None = None) → None[source]#

Add one homogeneous stratum of cells.

counts is n_cells x n_genes. stratum labels the homogeneity unit (participant x visit for simulated data; participant x visit x cell type for real data) and is used only for provenance – the caller is responsible for calling once per stratum.

conditional_alpha() → tuple[ndarray, ndarray][source]#

Pooled within-stratum NB2 moment dispersion and the per-gene mean.

Returns (alpha, mean_count_per_cell); alpha is NaN where the gene is not estimable.

genewise_corr_within_stratum(min_cpm: float = 1.0) → dict[source]#

Gene-wise pre/post correlation computed WITHIN each cell type.

The pooled version differences participants after summing over cell types, so a gene restricted to one cell type inherits that cell type’s ABUNDANCE variation across participants. That alone creates gene-to-gene heterogeneity in participant-level correlation, with no gene-intrinsic biology behind it.

The simulator contains one homogeneous population and no composition variation, so if TNBC’s heterogeneity is compositional then the pooled statistic is not a quantity the simulator could or should reproduce – the same conditional-versus-marginal distinction already applied to dispersion. Measuring within cell type is what makes the two arms of the gate comparable.

pv_frame() → DataFrame[source]#: Participant x visit pseudobulk counts, summed over strata.

statistics() → dict[source]#: Every gate statistic, as a flat dict of scalars plus a few vectors.

variance_components(min_cpm: float = 10.0, within_stratum: bool = False) → dict[source]#

LATENT participant and participant x visit SDs on the log-rate scale.

These are the quantities the simulator is parameterised by, and they are NOT the observable correlation. Conflating the two is the same error as calibrating dispersion on the marginal curve: the observable pre/post correlation of log(1+CPM) is attenuated by pseudobulk sampling noise, so setting the generating parameter equal to it under-disperses the hierarchy (measured: observable 0.348 for a generating 0.466).

Identification, with two visits per participant and one pseudobulk value each, comes from the split halves:

Var(y_A - y_B) = 4 * sigma_e^2          (halves share b and u)
Var(post - pre) = 2 sigma_u^2 + 2 sigma_e^2
Var((post + pre)/2) = sigma_b^2 + sigma_u^2/2 + sigma_e^2/2

Restricted to genes above min_cpm because log(1+x) ~ log(x) only there; below it the transform itself shrinks the variance and the components would be biased low.

sctrial.benchmark.calibration.summarize_blocks(blocks, n_genes: int) → SummaryAccumulator[source]#: Feed an iterable of (counts, participant, visit, arm, stratum) blocks.

sctrial.benchmark.calibration.summarize_adata(adata, participant_col: str = 'participant', visit_col: str = 'visit', arm_col: str | None = 'arm', celltype_col: str | None = 'cell_type', layer: str | None = None) → SummaryAccumulator[source]#

Summarise real data, conditioning on participant x visit x cell type.

celltype_col=None conditions on participant x visit only. That is correct for an already-homogeneous population and WRONG for a real mixed sample: on TNBC it inflates the conditional dispersion 2.8x by loading between-cell-type mean differences onto it.

sctrial.benchmark.calibration.summarize_simulation(cfg) → SummaryAccumulator[source]#

Summarise one simulated replicate through the identical statistic code.

Consumes sctrial.benchmark.simulator_v2.iter_pv_blocks(), so no replicate is ever materialised in full and the simulated arm of every gate is computed by exactly the same code as the real arm.

sctrial.benchmark.calibration.conditional_dispersion(adata, participant_col: str = 'participant', visit_col: str = 'visit', celltype_col: str | None = 'cell_type', layer: str | None = None, alpha_grid: ndarray | None = None, gene_chunk: int = 200, verbose: bool = True) → DispersionFit[source]#

Gamma-Poisson dispersion MLE within homogeneous strata, with EB shrinkage.

The model is Y_cg ~ NB2(mu_cg = s_c * lambda_kg, alpha_g) where k is the stratum (participant x visit x cell type) and s_c the library-size offset. alpha_g is maximised over a log-spaced grid on the Cox-Reid adjusted profile likelihood, refined by parabolic interpolation, then shrunk toward a mean-dependent trend with an empirical-Bayes weight.

Conditioning on cell type is not optional for real mixed samples: without it the estimate absorbs between-cell-type mean differences (TNBC: 0.774 pooled versus 0.275 within cell type).

Measure every simulator target from real data and write the canonical files.

ALL TARGETS ARE MEASURED WITHIN CELL TYPE. The primary simulator represents one homogeneous cell population, so cell-level dispersion, participant and participant-by-visit variance, cells per participant-visit, gene rates and the longitudinal covariance must all be conditioned at that level. Mixing conditional dispersion with pooled cell counts and pooled longitudinal covariance is the incoherence the calibration gates kept detecting.

Scalar targets are the MEDIAN across cell types (a typical cell type), with the inter-cell-type range recorded. Cell-type-pooled values are also written, with a pooled_ prefix, for description only – they are what a whole-sample analysis would see and are NOT calibration targets.

Validation gates#

Monte Carlo envelope tests (Gates A–E) that verify the simulator faithfully reproduces the empirical nuisance statistics it was calibrated against.

sctrial.benchmark.gates.GATE_STATISTICS: dict[str, list[str]]#: Maps each gate label ("A", "B", …) to the list of statistics it checks.

sctrial.benchmark.gates.PINNED_STATISTICS: frozenset[str]#: Statistics that are near-deterministic readbacks of an empirical pool the simulator resamples. Agreement is true by construction; failure indicates an implementation defect, not a fidelity defect.

class sctrial.benchmark.gates.GateResult(gate: str, statistic: str, observed: float, sim_median: float, sim_lo95: float, sim_hi95: float, percentile: float, verdict: str)[source]#

Bases: object

Envelope verdict for one statistic.

sctrial.benchmark.gates.run_gates(cfg, observed: dict, n_mc: int = 200, n_jobs: int = 8, seed0: int = 100000, out_dir: str | Path | None = None, verbose: bool = True, bootstrap: list[dict] | None = None) → DataFrame[source]#

Run every envelope gate.

Parameters:

cfg – A TranscriptomeSimConfig. It is used verbatim except for seed, which is varied across replicates.
observed – Real-data statistics from sctrial.benchmark.calibration.SummaryAccumulator.statistics().
n_mc – Number of simulated replicates. 200 gives a 2.5th-percentile bound with about 5 order statistics below it, which is the practical minimum for a 95% envelope; 500 is preferable when the compute allows.

sctrial.benchmark.gates.composition_ablation(cfg, architectures=('balanced', 'heterogeneous', 'one_directional'), signal_fraction: float = 0.2, magnitude: float = 0.5, panel_size: int = 200, n_rep: int = 5, seed0: int = 900000, verbose: bool = True) → DataFrame[source]#

Isolate normalisation scope from signal direction.

For each signal architecture the same generated cells are analysed three ways: normalised on the full transcriptome, normalised on the tested panel only, and against the noiseless oracle estimand. If the recovery gap tracks normalisation scope rather than signal direction, the compositional explanation is demonstrated rather than merely consistent with the data.

Metrics#

Per-replicate and per-scenario metric functions used to evaluate statistical methods.

sctrial.benchmark.metrics.compute_fpr(pvalues: ndarray, alpha: float = 0.05) → dict[source]#

Per-test false positive rate under the null.

Parameters:

pvalues (array) – P-values from null genes only.
alpha (float) – Significance threshold.

Return type:

dict with “fpr”, “n_tests”, “wilson_ci_lo”, “wilson_ci_hi”

sctrial.benchmark.metrics.compute_fdr_tpr(pvalues: ndarray, is_signal: ndarray, q: float = 0.05) → dict[source]#

FDR and TPR (power) at a given BH threshold.

Parameters:

pvalues (array) – P-values for all genes (null + signal).
is_signal (array of bool) – True for genes with non-zero effect.
q (float) – BH-FDR threshold.

Return type:

dict with “fdr”, “tpr”, “n_discoveries”, “n_true_positives”, “n_signal”

sctrial.benchmark.metrics.compute_bias_rmse(estimated: ndarray, true: ndarray) → dict[source]#

Bias and RMSE of effect-size estimates.

Parameters:

estimated (array) – Estimated β values.
true (array) – True β values.

Return type:

dict with “bias”, “rmse”, “n”

sctrial.benchmark.metrics.compute_ci_coverage(ci_lo: ndarray, ci_hi: ndarray, true: ndarray) → dict | None[source]#

CI coverage rate (95% nominal).

Returns None if no valid intervals are available (e.g., Wilcoxon). This follows locked rule 2: coverage only where intervals are defined.

sctrial.benchmark.metrics.compute_sign_recovery(estimated: ndarray, true: ndarray, threshold: float = 0.05) → dict[source]#: Fraction of correct effect-size signs for |β| > threshold.

sctrial.benchmark.metrics.compute_lambda_gc(pvalues: ndarray) → float[source]#

Genomic inflation factor from null p-values.

Secondary diagnostic (locked rule 3): use alongside FPR and QQ, not as a headline metric.

sctrial.benchmark.metrics.compute_topk_jaccard(ranking_a: Series, ranking_b: Series, k: int = 20) → float[source]#

Jaccard overlap of top-k genes between two rankings.

Rankings are Series indexed by gene name, values are p-values (lower = more significant).

sctrial.benchmark.metrics.compute_failure_rates(results: list[dict]) → dict[source]#

Compute failure rates split by mode (locked rule 4).

Parameters:: results (list of dict) – Each dict has “failure_mode”: None | “convergence” | “numerical” | “timeout”
Returns:: dict with keys
Return type:: “convergence_rate”, “numerical_rate”, “timeout_rate”, “total_failure_rate”, “n”

sctrial.benchmark.metrics.summarize_iteration(results: dict[str, dict], truth: dict[str, float], signal_genes: set[str]) → dict[source]#

Compute all metrics for a single method on a single iteration.

Parameters:

results (dict) – gene_name → {“beta”, “pvalue”, “ci_lo”, “ci_hi”, “converged”, “failure_mode”}
truth (dict) – gene_name → true effect size
signal_genes (set) – Gene names with non-zero true effects

Return type:

dict of all metrics

Endpoints#

Pre-specified statistical endpoints computed from replicate-level results.

sctrial.benchmark.endpoints.ALPHA: float = 0.05#: Locked nominal significance threshold.

sctrial.benchmark.endpoints.Q_FDR: float = 0.05#: Locked Benjamini–Hochberg FDR threshold.

sctrial.benchmark.endpoints.replicate_endpoints(g: DataFrame, alpha: float = 0.05, q: float = 0.05) → dict[source]#

Endpoints for ONE method within ONE replicate.

g must contain every gene ASSIGNED to the panel, including those the method could not evaluate; that is what makes the end-to-end denominator meaningful.

sctrial.benchmark.endpoints.scenario_endpoints(df: DataFrame, alpha: float = 0.05, q: float = 0.05) → DataFrame[source]#

Endpoints per scenario and method, with replicate-level Monte Carlo error.

Each rate is computed within a replicate and then averaged; the reported *_mcse is the standard error over replicates. Averaging genes across replicates instead would understate it, because genes in one replicate are not independent of each other.

Orchestrator#

Monte Carlo orchestration, scenario grid construction, and the top-level entry points for running and sensitivity-testing the benchmark.

sctrial.benchmark.CORE_METHODS: list[str]#: The methods included in every standard benchmark run: sctrial_did, wilcoxon_paired, dreamlet, limma_voom, nebula.

sctrial.benchmark.run_benchmark(designs: list[str] | None = None, methods: list[str] | None = None, n_iterations: int = 200, n_jobs: int = 1, output_dir: str | Path = 'benchmark_results', resume: bool = True, base_config: dict | None = None, manifest: dict | None = None, completion_dir: str | Path | None = None) → DataFrame[source]#: Run the core scenario grid.

sctrial.benchmark.run_sensitivity_benchmark(designs: list[str] | None = None, methods: list[str] | None = None, n_iterations: int = 200, n_jobs: int = 1, output_dir: str | Path = 'benchmark_results/sensitivity', resume: bool = True, base_config: dict | None = None, manifest: dict | None = None, panels: list[int] | None = None, completion_dir: str | Path | None = None) → DataFrame[source]#: Run the panel-size x signal-fraction sensitivity grid.

sctrial.benchmark.orchestrator.build_scenario_grid(design: str = 'two_arm') → list[dict][source]#: Core scenario grid for one design family.

sctrial.benchmark.orchestrator.build_sensitivity_grid(design: str = 'two_arm', panels=None) → list[dict][source]#

Panel size x signal BURDEN sensitivity, declared as integer gene counts.

Panels are NESTED subsets of one simulated transcriptome, so a panel-size effect is separable from a gene-identity effect.

SIGNAL IS AN INTEGER COUNT, and the fractions are chosen to be EXACTLY realisable at every panel size, so the grid is a complete factorial:

fraction     50   200   500  2000
    2%        1     4    10    40
    4%        2     8    20    80
   10%        5    20    50   200
   20%       10    40   100   400

The earlier 1% and 5% are not realisable at 50 genes (0.5 and 2.5 genes). Rounding them would mislabel the lowest-signal condition – one gene out of 50 is 2%, not 1% – and skipping them instead leaves holes in the very panel-size comparison the analysis exists to make. Nothing is scientifically privileged about 1% and 5%; they were sensitivity values. 2% and 4% are more rigorous because the experimental factor is then exactly defined everywhere.

sctrial.benchmark.orchestrator.mc_max_for(scenario: dict) → int[source]#

The replicate cap for one scenario.

Keyed on the SIGNAL GENE COUNT, not on the number of genes tested. A 2,000-gene scenario does not carry 40x more independent information than a 50-gene one: genes within a replicate share participants, libraries, random effects, normalisation and signal architecture. The independent unit is the simulated dataset, which is the same pseudoreplication argument the manuscript itself makes.

sctrial.benchmark.orchestrator.scenario_seed(master_seed: int, scenario_name: str, replicate: int) → int[source]#

A seed addressed by (scenario, replicate index), not by draw order.

With a single sequential RNG stream, adding a batch of replicates under adaptive stopping would renumber every later dataset, so replicate 437 would not be the same data across a resume, a different worker count, or a re-partitioned SLURM job. Hashing the address makes replicate 437 always replicate 437.

hashlib rather than hash(): Python salts str hashing per process, so hash() is not reproducible across runs at all.

Permutation & subsampling#

Real-data hypothesis tests and sensitivity analyses that operate on observed AnnData objects rather than simulated data.

sctrial.benchmark.permutation.run_permutation_test(adata, gene_cols: list[str], design_type: str = 'two_arm', methods: list[str] | None = None, n_permutations: int = 1000, n_jobs: int = 1, participant_col: str = 'participant', arm_col: str = 'arm', visit_col: str = 'visit', output_path: str | Path | None = None, seed: int = 42) → DataFrame[source]#

Run participant-label permutation test on real data.

Parameters:

adata (AnnData) – Real dataset.
gene_cols (list[str]) – Genes to test.
design_type ({"two_arm", "single_arm"})
methods (list[str]) – Methods to run. Default: all core methods.
n_permutations (int) – Number of permutations.
n_jobs (int) – Parallel workers.

Returns:

DataFrame with columns

Return type:

permutation, method, gene, pvalue

sctrial.benchmark.subsample.run_subsampling(adata, gene_cols: list[str], methods: list[str] | None = None, fractions: list[float] | None = None, n_resamples: int = 100, participant_col: str = 'participant', arm_col: str = 'arm', visit_col: str = 'visit', output_path: str | Path | None = None, seed: int = 42) → DataFrame[source]#

Run subsampling reproducibility on a single dataset.

Parameters:

adata (AnnData) – Full dataset.
gene_cols (list[str]) – Genes to test.
methods (list[str]) – Methods to benchmark.
fractions (list[float]) – Participant fractions to subsample. Default: [0.5, 0.7, 0.9].
n_resamples (int) – Number of random subsamples per fraction.

Returns:

DataFrame with

Return type:

fraction, resample, method, spearman_rho, jaccard_top20

Ablation#

Estimator variants used to isolate the contribution of individual design choices (normalisation scope, aggregation level, fixed effects).

sctrial.benchmark.ablation.ABLATION_VARIANTS: dict[str, tuple]#: Registry mapping variant name to (input_type, runner, label).

sctrial.benchmark.ablation.run_ablation(inputs: dict, gene_cols: list[str], variants: list[str] | None = None) → dict[str, dict][source]#

Run ablation variants on a single simulated/real dataset.

Parameters:

inputs (dict) – From sctrial.benchmark.contracts.prepare_inputs(). Supplying the prepared contract rather than a raw sim dict is what guarantees every rung sees the same outcome.
gene_cols (list[str]) – Genes to test.
variants (list[str]) – Which ablation variants to run. Default: all.

Returns:

dict

Return type:

variant_name → {gene → {“beta”, “pvalue”, “ci_lo”, “ci_hi”}}

Infrastructure#

Reproducibility and result-layout utilities shared across the benchmark pipeline.

Result layout

class sctrial.benchmark.paths.ResultLayout(root: Path | str, manifest_sha: str)[source]#

Bases: object

The directory tree for one manifest’s results.

completed_scenarios(grid: str | None = None) → dict[str, dict][source]#

Scenarios with a completion record, keyed by scenario id.

A CSV without a record is NOT complete: it is what a job killed part-way through an adaptive extension leaves behind.

completion_marker(grid: str) → Path[source]#

The marker written by the aggregator for ONE grid, at the end.

Per GRID, not per run. The core and sensitivity grids are aggregated by separate jobs into the same manifest directory, so a single shared marker would be written twice and the survivor would attest to whichever finished last – the same last-writer-wins defect as the combined file itself, reintroduced one level up.

orphan_scenarios(grid: str | None = None) → list[str][source]#: Scenario CSVs with no completion record – truncated or in flight.

publication_marker() → Path[source]#

The WHOLE-BENCHMARK marker, written only by the finalizer.

Grid-level markers say “core finished” or “sensitivity finished”. Neither says “the manuscript benchmark is complete”, and the asymmetric failure is real: core aggregates successfully, sensitivity never finishes, and a figure script finds a perfectly valid core marker and regenerates outputs from half the benchmark.

This project has already met both other members of that class – several jobs writing one combined file, and two grids sharing one scenario directory – so the third is closed here rather than left to discipline.

scenarios_for(grid: str) → Path[source]#

Scenario CSVs for ONE grid.

Separated because the aggregator globs this directory and validates the result against the exact set the grid declares. The core and sensitivity grids write into the same manifest directory, so a shared directory would hand the core aggregator all 36 sensitivity files as UNEXPECTED – and the natural “fix” of filtering to the expected set would destroy the only check that catches a genuinely unexpected scenario.

sctrial.benchmark.paths.preflight_layout(root: Path | str, label: str) → ResultLayout[source]#

A clearly-separated tree for pre-freeze probes.

Uses the same class so a probe exercises the real execution path, but under _preflight/ so it is obvious that nothing here is a manuscript result and a single directory removal disposes of all of it.

sctrial.benchmark.paths.require_layout(root: Path | str, manifest_sha: str) → ResultLayout[source]#

Resolve a layout that must already contain a completed run.

This is the loader entry point. It refuses to guess: no “latest”, no newest CSV, no search across manifests.

Scenario contracts

class sctrial.benchmark.scenario_contract.Violation(field: str, requested: Any, realised: Any, detail: str)[source]#

Bases: object

One requested-versus-realised mismatch.

class sctrial.benchmark.scenario_contract.ContractReport(scenario: 'str', violations: 'list[Violation]', checked: 'list[str]', unchecked: 'list[str]')[source]#: Bases: object

sctrial.benchmark.scenario_contract.check_simulation(scenario: dict, cfg: Any, sim: dict, panel_genes: list, signal_genes: set) → ContractReport[source]#

Validate one realised simulation against the scenario that requested it.

Runs per iteration, immediately after simulate_trial_v2, so a mismatch is caught on the first replicate rather than after a scenario has burned a 72-hour allocation.

sctrial.benchmark.scenario_contract.check_scenario_results(scenario_id: str, scenario: dict, df: DataFrame, manifest_sha: str | None = None) → ContractReport[source]#

Validate a completed scenario’s result table before it is marked complete.

check_simulation guards each replicate as it is produced. This is the second gate: it re-derives the same quantities from the RECORDED columns, so a defect in the recording path – rather than in the simulation path – is also caught. The two are deliberately redundant; the original defect survived precisely because only one layer was checked.

sctrial.benchmark.scenario_contract.evaluability_by_method(df: DataFrame) → dict[str, float][source]#

Fraction of prespecified genes each method actually returned inference for.

Reported, never used to filter. A method that drops genes at low cell yield has a lower end-to-end detection rate, and that is a benchmark result rather than a reason to shrink its denominator.

sctrial.benchmark.scenario_contract.completion_record(scenario_id: str, scenario: dict, df: DataFrame, stop_reason: str, max_replicates: int, manifest_sha: str | None, mcse_target_fpr: float, mcse_target_power: float) → dict[source]#

The per-scenario completion record.

Written ONLY after a scenario finishes, so its absence marks truncation. This matters because adaptive stopping makes replicate count a per-scenario outcome: a scenario killed part-way through an extension still holds more rows than the base batch, so a count threshold cannot distinguish it from one that legitimately stopped early. The record can.

Manifest & reproducibility

sctrial.benchmark.manifest.SOURCE_TREE_PATHS: tuple[str, ...]#: Paths hashed when computing the source-tree digest: src, scripts, pyproject.toml.

sctrial.benchmark.manifest.manifest_hash(manifest: dict) → str[source]#: Stable hash of a manifest, independent of key order.

sctrial.benchmark.manifest.verify_manifest(manifest: dict, artifacts: dict[str, Path] | None = None) → None[source]#: Raise if the recorded hashes no longer describe what is on disk.

sctrial.benchmark.manifest.assert_single_manifest(df: DataFrame, context: str = '') → str[source]#

Refuse a table whose rows came from different runs.

Combining results across manifests is how a corrected run gets silently averaged with the run it was meant to replace. Loaders call this before plotting anything.

sctrial.benchmark.manifest.source_tree_sha256(repo: Path | None = None, paths=('src', 'scripts', 'pyproject.toml')) → str[source]#

Deterministic hash of the source that will actually run.

A commit SHA says what SHOULD be there; this says what IS there. The difference is not hypothetical here: the cluster spent this project with its HEAD pinned at one commit while rsync had overwritten the files with code many commits newer, so the nominal commit described nothing that was executing.

It also needs no git, which matters because git is absent from this cluster’s compute nodes – so a job can verify its own source at run time, which is exactly where verification is worth having.

Sorted relative paths, content-hashed, excluding bytecode and egg-info so the hash is stable across installs.

Seeds

sctrial.benchmark.seeds.SEP: str#: Field separator used when constructing seed strings (ASCII unit separator, \x1f).

sctrial.benchmark.seeds.stable_seed(namespace: str, *parts: object, bits: int = 64) → int[source]#

A deterministic seed derived from a namespace and identifying parts.

Identical across processes, interpreters and PYTHONHASHSEED values, which is the entire point: hash() is not.

Parameters:

namespace – Distinguishes unrelated uses of the same identifiers. Include a version suffix ("..._v1") when a change to the derived quantity should deliberately produce different draws.
parts – Identifying components. Converted with str and joined with an unambiguous separator.
bits – Width of the returned integer. 64 by default; NumPy accepts arbitrary non-negative integers as seeds, and 64 bits is far past any collision concern at this scale.